Alongshore transport

Alongshore transport#

First import some necessary packages#

import logging
import pathlib
import sys
import warnings

import bokeh.io
import colorcet as cc
import dotenv
import geopandas as gpd
import holoviews as hv
import hvplot.pandas  # noqa: API import
import numpy as np
import pandas as pd
import panel as pn
import pooch
from bokeh.models import HoverTool, PanTool, WheelZoomTool
from bokeh.resources import INLINE

import coastal_dynamics as cd

# Activate Panel extension to make interactive visualizations
pn.extension()

# # Read questions from cloud storage
questions = cd.read_questions(
    # "./7_alongshore_transport.json"
    "az://coastal-dynamics/questions/7_alongshore_transport_hashed.json",
    storage_options={"account_name": "coclico"},
)

%run initialize/init_7.ipynb

Alongshore sediment transport#

Welcome to the notebook of week 7! The main topic of this notebook is alongshore sediment transport. This notebook is relatively short. We will look into the CERC formula, and how to use it to obtain an (S, $\phi$-curve. Afterwards, we will use the (S, $\phi$-curve to determine coastal evolution around a breakwater.

CERC formula and ($S$,$\phi$)-curve#

Many different formulas exist to calculate bulk longshore sediment transport. One widely used formula is the CERC formula (Equation 8.4 in the book). We will use a version here which uses deep-water values for the wave height and wave angle (Equation 8.10 in the book). This formula, which is applicable for straight, parallel depth contours, is defined as: $$S = \frac{K}{32(s-1)(1-p)} c_b \sin{2 \phi_0} H_0^2 $$ For a complete description of each parameter, see section 8.2.3 of the book.

In line with Example 8.1 in the book, we will assume that the offshore wave height $H_0$ is equal to the deep-water root-mean-square wave height $H_{rms,0}$, and can be represented by a single value of 2 m, with a period T of 7 s. The corresponding value for $K$ is given by $K_{rms} \approx 0.7$. We further use a porosity $p = 0.4$ and relative density $s = 2.65$. We use $\gamma = H_{rms,b}/h_b = 0.8$.

You might notice that, given the offshore wave height, the only variables left in the CERC formula are the wave celerity at breaking $c_b$ and the offshore wave angle $\phi_0$. Starting from our values for $H_0$ and $T$, we can use linear wave theory to compute values for $c_b$ as a function of $\phi_0$. The computation procedure is described in Example 8.1 in the book. The function below takes care of this computation. Each step in the function is explained and it is instructive to see if you can follow the function step-by-step. The result of the function is $c_b$ for a certain deep water wave angle of incidence $\phi_0$. Once $c_b$ is known, we can compute the transport $S$. Note that the function also outputs $\phi_b$.

In the computation, we use an efficient solver for the dispersion relationship (© Gert Klopman, 1994; conversion to Python by M. van der Lugt). You could also use the dispersion relationship that you programmed in Waves or in the notebooks of week 2. If you would like to do so, implement a function by uncommenting the cell below and adding your code.

# def disper(w, h, g):
#     """
#     Input:
#     w = 2*pi/T, were T is wave period
#     h = water depth
#     g = gravity constant

#     Output:
#     k = wave number
#     """

#     k = ...

#     return k

def find_cb_phib(
    phi0, H0, T, g=9.81, gamma=0.8, hb=np.arange(0.1, 5.0, 0.01), print_report=False
):
    """
    Returns breaking wave celerity cb [m/s] and angle of incidence at breaking phib [degrees] for given:
    - phi0 : angle of incidence [degrees]
    - H0   : deep water wave height [m]
    - T    : period [s]

    The parameter hb_guess is used as guessed values for the breaking depth.
    From this array, the best-fitting value is chosen in the end. You can adjust this
    array to make estimates more accurate at the cost of computational efficiency.
    """
    # First convert the angle of incidence to radians
    phi_rad = phi0 / 360 * 2 * np.pi

    # We start with calculating deep water celerity, wavelength, and angular frequency
    c0 = g * T / (2 * np.pi)
    L0 = c0 * T
    w = T / (2 * np.pi)

    # For every value of hb_guess, the wavenumber k is determined using the dispersion relation
    k = disper(w, hb, g=g)  # Feel free to use your own implementation from week 2!

    # Next we calculate the celerity and group celerity for each breaking depth
    c = np.sqrt(g / k * np.tanh(k * hb))
    n = 1 / 2 * (1 + (2 * k * hb) / (np.sinh(2 * k * hb)))
    cg = n * c

    # In order to correctly shoal the waves, we also need the deep water group celerity
    n0 = 1 / 2
    cg0 = n0 * c0

    # And to account for refraction we need the angle of incidence at breaking using Snell's law
    phi = np.arcsin(np.sin(phi_rad) / c0 * c)

    # Shoaling & refraction coefficients
    Ksh = np.sqrt(cg0 / cg)
    Kref = np.sqrt(np.cos(phi_rad) / np.cos(phi))

    # Wave heights Hb at depth hb
    Hb = Ksh * Kref * H0

    # We are looking for an hb where the breaker parameter is 0.8
    # We can determine which value of hb in our array gets closest using the
    # following line of code:
    i = np.argmin(np.abs(Hb / hb - gamma))
    Hb_pred, hb_pred = Hb[i], hb[i]

    # Let's print what we found
    if print_report:
        print(f"predicted breaking depth: {hb_pred:.2f} m")
        print(f"predicted breaking wave height: {Hb_pred:.2f} m")
        print(f"gamma = {Hb_pred / hb_pred:.2f} [-]")

    # And finally return the associated value cb for the celerity at breaking, as well as the angle of incidence at breaking phib
    return c[i], phi[i] / (2 * np.pi) * 360

Let’s check that it works. The cell below shows that for an offshore wave height of 2 m, a period of 7 s, and a deep water angle of incidence of 5 degrees, we get a $c_b$ of 4.92 m/s. You can check that this is in reasonable agreement with Example 8.1 in the book. Also check the other values that the cell below prints.

phi0 = 5  # degrees
H0 = 2  # m
T = 7  # s

cb, phib = find_cb_phib(phi0, H0, T, print_report=True)

print(f"cb: {cb:.2f} m/s")
print(f"phib: {phib:.2f} degrees")

predicted breaking depth: 2.79 m
predicted breaking wave height: 2.23 m
gamma = 0.80 [-]
cb: 4.92 m/s
phib: 2.25 degrees

Now that we have this function, we can calculate a breaking wave celerity for each angle of incidence. We use a range from negative to positive angles here, but this is strictly not necessary since the transport S for negative angles only differs in sign from the result for positive angles .

# Set the range for which we want to calculate cb
phi_array = np.arange(-80, 81, 1)

# Initialize cb array
cb_array = np.zeros(phi_array.shape)
phib_array = np.zeros(phi_array.shape)

# Loop through each phi and compute associated value for cb
for i in range(len(phi_array)):
    cb_array[i], phib_array[i] = find_cb_phib(phi_array[i], H0, T)

With $c_b$ as a function of $\phi_0$, we can calculate our transport $S$ as a function of $\phi_0$, which means we can generate an ($S$,$\phi$)-curve!

Remember that we can use the CERC formulation of Equation 8.10 for this. We use typical values of $K=0.7$, $p=0.4$, and $s=2.65$.

Check the equation below to compute the bulk sediment transport S.

def CERC(cb, phi0, H0, K=0.7, s=2.65, p=0.4):
    """
    cb:   celerity at breaking
    phi0: offshore angle of incidence (degrees)
    H0:   offshore wave height

    K:    coefficient
    s:    relative density
    p:    porosity
    """

    return (
        K / (32 * (s - 1) * (1 - p)) * cb * np.sin(2 * (phi0 / 360 * 2 * np.pi)) * H0**2
    )

First, let’s see what our value for transport is for the previously used angle of incidence of 5 degrees.

S_5 = CERC(cb, 5, H0)

print(f"S for angle of incidence of 5 degrees: {S_5:.4f} m3/s")

S for angle of incidence of 5 degrees: 0.0756 m3/s

This is in reasonable agreement with Example 8.1 from the book! You can of course also change $\phi_0$ in the above computation of $c_b$ and $S$ to check the values for other deep water angles.

Let’s continue with calculating S for a range of values of $\phi$, in order to generate an ($S$, $\phi$)-curve.

S = CERC(cb_array, phi_array, H0)

warnings.filterwarnings("ignore")
logging.getLogger().setLevel(logging.ERROR)

(
    hv.Curve((phi_array, S))
    * hv.HLine(0).opts(color="black")
    * hv.VLine(0).opts(color="black")
).opts(
    xlabel="deep water angle [degrees]",
    ylabel="S [m3/s]",
    title="(S, phi_0)-curve",
    width=800,
    height=400,
    show_grid=True,
)

This plot is comparable to Figure 8.4 in Example 8.1 in the book.

Using this plot and the code above, try to answer the questions below.

q1 = cd.QuestionFactory(questions["Q7-1"]).serve()
q2 = cd.QuestionFactory(questions["Q7-2"]).serve()
q3 = cd.QuestionFactory(questions["Q7-3"]).serve()
q4 = cd.QuestionFactory(questions["Q7-4"]).serve()

pn.Column(q1, q2, q3, q4)

This is the end of the first part of this notebook.

Coastal evolution around a breakwater#

We will now use the ($S$,$\phi$)-curve generated above to predict the coastal evolution updrift of a breakwater using single-line theory. For now, we will consider only a single wave condition, which is thought to be representative:

$H_0=1.2$
$T=7$
$\phi_0=10$

The angle is quite small, but it simplifies our notebook and helps a comparison that we will make below between a numerical and analytical solution.

This condition occurs roughly 10% of the time. Note that that means that for 90% of the time no significant transport occurs. Remember that positive angles result in positive transport (see also the image in Table 8.2).

Let’s determine the average yearly transport for this. We just calculate the transport S, and convert it to m3/year:

H0 = 1.2
T = 7
phi0 = 10

cb, phib = find_cb_phib(phi0, H0, T)
S = CERC(cb, phi0, H0)

# This transports S is given in m3/s. Remember that the wave climate occurs 10% of the year. For 1 year, the contribution is:
S_total = (0.1 * S) * 365.25 * 24 * 3600

print(f"Breaking wave celerity {cb:.2f} m/s")
print(f"Angle at breaking {phib:.2f} degrees")
print(f"Transport per second: {S:.3f} m3/s")
print(f"Total yearly transport: {S_total:.0f} m3/year")

Breaking wave celerity 4.04 m/s
Angle at breaking 3.68 degrees
Transport per second: 0.044 m3/s
Total yearly transport: 138817 m3/year

This total transport is positive, which means the cumulative transport is in the positive direction. But what does this mean for the coastal morphology? According to Equation 8.16, we know the coastline will evolve by: $$\frac{\partial Y}{\partial t} + \frac{1}{d} \frac{\partial S_x}{\partial x}=0$$ with $S_x$ the sediment transport in the $x$-direction. This equation means that the temporal evolution of the coastline position Y is determined by the *alongshore gradient* of the alongshore transport. For any location $x$ along the coastline and a given active profile height $d$, we can determine the coastline evolution through the following steps:

Determine the angle of wave attack relative to the coastline $\phi$
Determine the sediment transport, which is a function of the local angle $\phi$
Determine the alongshore gradient of the alongshore sediment transport
Determine the expected morphological change in a certain short time interval

However, when the coastline orientation changes, the angle of wave attack relative to the coastline changes as well and therewith the transport. We thus need to iterate through these steps. This is called the morphological feedback loop. When repeated, for example through a numerical model, this strategy can help predict morphological changes. This is useful for predicting the shoreline evolution due to for instance human interventions! Let us for example consider a breakwater perpendicular to the shore. We will try to use a numerical model to determine the shoreline evolution. Don’t worry too much about the numerical model, for now, we will focus on the concepts here! Therefore, and because it may take a little while to run, we have run the model already for you, and will simply load the results. You will learn more about the numerical modelling part in the Coastal Modelling unit of the Coastal Engineering track and also do this yourself.

Let’s place the breakwater at $x=0$, and consider the domain $x \leq 0$. The breakwater extends into the water from $y=0$ m to $y=500$ m. An overview of the situation is provided by running the cell below.

breakwater = hv.Curve(((0, 0), (0, 500)), label="Breakwater").opts(color="black")
shoreline = hv.Curve(((-10000, 500), (0, 0)), label="Shoreline").opts(color="#D2B48C")
(breakwater * shoreline).opts(
    width=1200,
    height=400,
    legend_position="top_left",
    xlabel="alongshore position [m]",
    ylabel="cross-shore position [m]",
)

Firstly, some assumptions need to be made. For simplicity, we set $d$ to an assumed depth of closure of $d=7$ m. Secondly, let’s assume that the wave climate is given by the single condition described above (for which we already calculated the total yearly transport). We also need to impose initial and boundary conditions. These are thoroughly described by Equation 8.21 - 8.23. Briefly, for initial conditions, we assume a horizontal coastline (i.e $y=0$ along the coast). The following two boundary conditions are imposed:

$S_x = S$, for $x=-\infty$ and for all $t$
$S_x = 0$, for $x=0$ and for all $t$

We will look at a stretch of coast with a length of order 10 km. As mentioned, we have already run the numerical model for you. Using the cell below, you can simply load the results!

Remember, the angle $\phi$ is defined as positive when it induces positive transport, the same as the image in Table 8.2.

fp_X = pooch.retrieve(
    "https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_X.txt",
    known_hash="908f9b3871098c0446fb3b6c8933bd62cfeb83df33db5469814f6516be596767",
)
fp_T = pooch.retrieve(
    "https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_T.txt",
    known_hash="5e71b53e92b4493c5e8a0b2b29fcd1bcd1620645a561896df309883982155ff9",
)
fp_Y = pooch.retrieve(
    "https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_Y_t.txt",
    known_hash="480eb21ccfad136acc81ced0c915cc4ebcaed36259ec22148f5c0d295ed32f16",
)
fp_S = pooch.retrieve(
    "https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_S_t.txt",
    known_hash="47500101458117c5041cb9a0afd96fd15306fa7d9bb7e48ae25c91ec3ad63970",
)

X = np.loadtxt(fp_X)
T = np.loadtxt(fp_T)
Y_t = np.loadtxt(fp_Y)
S_t = np.loadtxt(fp_S)

Downloading data from 'https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_X.txt' to file '/home/runner/.cache/pooch/9646315e3d0e10c953d5c4ccddff42e2-7_X.txt'.

Downloading data from 'https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_T.txt' to file '/home/runner/.cache/pooch/43d7b53d39a97c47a760f7e91fc1d8be-7_T.txt'.

Downloading data from 'https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_Y_t.txt' to file '/home/runner/.cache/pooch/6853584fb91e61a513b51cd3d118c5e9-7_Y_t.txt'.

Downloading data from 'https://coclico.blob.core.windows.net/coastal-dynamics/7_alongshore_transport/7_S_t.txt' to file '/home/runner/.cache/pooch/7fcf093c61da351fed871d6b6ef6931e-7_S_t.txt'.

Let’s plot the results for a selection of years. You can modify which years to display if you want. We plot the total yearly transport (in the top plot), as well as the coastline (in the bottom plot). Each line represents the coastline after a certain amount of years.

years = [0, 5, 10, 15, 20]
# years = np.arange(0, 21, 1)   # Uncomment this line if you want a visualization for every year!

transport_plot = hv.Curve(((), ()))
coastline_plot = hv.Curve(((), ()))

for year in years:

    id = np.argmin(np.abs(T - year))

    transport_plot *= hv.Curve((X, S_t[id]), label=f"Year = {year}")
    coastline_plot *= hv.Curve((X, Y_t[id]), label=f"Year = {year}")

breakwater = hv.Curve(((0, 0), (0, 500)), label="Breakwater").opts(color="black")
shoreline = hv.Curve(((np.min(X), 500), (0, 0)), label="Shoreline").opts(
    color="#D2B48C"
)

(
    transport_plot.opts(
        width=1200,
        height=400,
        legend_position="bottom_left",
        xlabel="alongshore position [m]",
        ylabel="transport S [m3/year]",
        show_grid=True,
    )
    + (coastline_plot * breakwater * shoreline).opts(
        width=1200,
        height=400,
        legend_position="top_left",
        xlabel="alongshore position [m]",
        ylabel="cross-shore position [m]",
        show_grid=True,
    )
).opts(shared_axes=False).cols(1)

This numerical solution looks quite nice. You can verify the transport magnitude far away from the breakwater and at the breakwater from the boundary conditions and the above-computed value for the total yearly transport. Have a look at how the shoreline and transport changes over the years and compare this also to the more qualitative Figures 8.11 and 8.13 in the book. You can see that the development of the shoreline slows down in time. You can also see the difference between the initial transport at t = 0 and the transport after some time.

In the plot with the morphological evolution through time, we see that the angle the shoreline makes with the breakwater is fairly constant throughout the years. This is the angle that leads to zero alongshore transport for the given conditions at this location! Can you estimate the shoreline angle at the breakwater, explain its value and why it is constant in time?

Note that in the textbook, lectures and exercises, we often only consider the initial transport and the morphological change as a consequence of this initial transport. So that would correspond to the first morphodynamic time step only.

Above we have considered a simplified situation, for which the textbook also presents some analytical solutions, for 1) the accretion length L(t) at the breakwater, and 2) the region of influence of the breakwater X = $5 \sqrt{at}$ (see Figure 8.12 and Equation 8.25). Let’s see how they compare with the numerical solution. In the first cell below we get the accretion length and region of influence for the modelled time. In the second cell, we calculate these same values analytically. Then we plot both to compare.

# minimal accretion to include it in region of influence

threshold = 0.0001  # times accretion length at

acc_num = Y_t[:, -1]
roi_num = np.zeros(T.shape)

for it in range(len(T)):
    X_influence = X[Y_t[it, :] > threshold * acc_num[it]]
    roi_num[it] = np.max(X_influence) - np.min(X_influence)

In Equation 8.25, $\phi$’ is the angle of incidence at the closure depth of $h = 7$ m. Here we will use $\phi$’$≈ \phi_0$, where $\phi_0$ is a small angle in line with the assumptions behind the analytical solution.

We also need (yearly) transport in the undisturbed zone (where the coastline still has its original orientation). We already computed this above as $S_{total}$ as follows: $S = CERC(c_b, 10, 1.2)$ This transport S is given in m$^3$/s. Remember that the wave climate occurs 10% of the year. For 1 year, the contribution is: $S_{total} = (0.1 * S) * 365.25 * 24 * 3600$

phi_prime = phi0

# For the depth of closure, we have previously assumed d=7
d = 7  #  m

# We can calculate our accretion length analytically
acc_ana = 2 * np.sqrt((phi_prime / 360 * 2 * np.pi) * S_total * T / (np.pi * d))

# And from that we can calculate our region of influence analytically
roi_ana = 2.5 * np.sqrt(np.pi) * acc_ana / (phi_prime / 360 * 2 * np.pi)

Let’s plot these solutions in the same plot. Note that the numerical solution in the right plot is influenced by the threshold chosen above as 0.01% of the accretion length at the breakwater at each time. This is the minimum accretion length before we count it as change. You can change it if you’d like and see how it affects the results.

acc_plot = (
    hv.Curve((T, acc_num), label="numerical")
    * hv.Curve((T, acc_ana), label="analytical")
).opts(
    width=600,
    height=400,
    legend_position="top_left",
    xlabel="time [years]",
    ylabel="Accretion length at breakwater L [m]",
    show_grid=True,
)
roi_plot = (
    hv.Curve((T, roi_num), label="numerical")
    * hv.Curve((T, roi_ana), label="analytical")
).opts(
    width=600,
    height=400,
    legend_position="bottom_right",
    xlabel="time [years]",
    ylabel="Region of influence X of breakwater [m]",
    show_grid=True,
)

(acc_plot + roi_plot).opts(shared_axes=False)

Using these plots, try to answer the questions below.

q5 = cd.QuestionFactory(questions["Q7-5"]).serve()
q6 = cd.QuestionFactory(questions["Q7-6"]).serve()
q7 = cd.QuestionFactory(questions["Q7-7"]).serve()
q8 = cd.QuestionFactory(questions["Q7-8"]).serve()

pn.Column(q5, q6, q7, q8)

We have been working with very simple wave conditions and boundary conditions. We could improve and extend this model further by using more wave conditions and different boundary conditions. The questions below concern this extension.

q9 = cd.QuestionFactory(questions["Q7-9"]).serve()
q10 = cd.QuestionFactory(questions["Q7-10"]).serve()

pn.Column(q9, q10)

We could extend this model with all sorts of additional functionality. Perhaps a first step would be to include the down-drift zone into the model. You will learn all about this in the Coastal Modelling unit!

You have reached the end of this section of the notebook.

Transport for wave climate as a function of coastline orientation#

In the numerical coastline model, a full wave climate is often taken into account. Further, the transport is computed by changing the coastline orientation relative to the waves. First, the transport for the full climate is computed for zero coastline orientation. Then, the coastline orientation changes a bit and the transport for the full climate is recomputed. This is equivalent to constructing an ($S$,$\phi$)-curve.

As a final step, let’s use an entire wave climate and plot the transport as a function of the shoreline orientation.

Note: previously, we plotted the ($S$,$\phi$)-curve for a single wave condition only. We will now include multiple wave conditions (the same as Intermezzo 8.2) and will determine the transport as a function of the coastline orientation. This will result in a different curve, as you will see. For the wave climate, see Table 8.2. Each wave height and offshore angle of incidence is associated with a period assumed to be 7 seconds.

First, we define a dataset containing the entire wave climate.

Hs_array = np.array([0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 2, 2, 2, 2, 4, 4, 4, 4])

phi_array = np.array(
    [-45, -15, 15, 45, -45, -15, 15, 45, -45, -15, 15, 45, -45, -15, 15, 45]
)

days_array = np.array([10, 10, 25, 8, 8, 15, 35, 16, 4, 10, 21, 11, 1, 1, 4, 2])

df = pd.DataFrame(
    data=zip(Hs_array, phi_array, days_array), columns=["Hs", "phi", "days"]
)

And we define a function that calculates the transport S for that wave climate, given a coastline orientation.

def get_S_coastline(coastline_orientation, wave_climate):
    """
    Returns yearly transport for angle phi [degrees]

    Transport is already scaled for the relative occurrence of the conditions.
    """
    total_transport = 0

    for index, row in df.iterrows():

        Hs, angle, days = row

        # Our formulation of the CERC formula (the choice of K) was based on Hrms, so we determine Hrms from Hs
        Hrms = Hs / np.sqrt(2)
        cb, phib = find_cb_phib(angle - coastline_orientation, Hrms, 7)
        S = CERC(cb, angle - coastline_orientation, Hrms)

        total_transport += days / 365.25 * S
    return total_transport * 365.25 * 24 * 3600

The cell below loops through different values of coastal orientation, and calculates the total yearly transport for each of these orientations.

orientations = np.linspace(-45, 45, 200)
transports = np.zeros(orientations.shape)

for i in range(len(orientations)):
    transports[i] = get_S_coastline(orientations[i], df)

(
    hv.Curve((orientations, transports))
    * hv.HLine(0).opts(color="black")
    * hv.VLine(0).opts(color="black")
).opts(
    xlabel="coastline orientation [degrees]",
    ylabel="S [m3/year]",
    title="(S, phi)-curve",
    width=800,
    height=400,
    show_grid=True,
)

The above figure is quite similar to Figure 8.10 in the book. Now you can try to answer the following questions. You can use the cell below for this if you’d like.

coastline_orientation = 15

print(
    f"Total yearly transport: {get_S_coastline(coastline_orientation, df):.0f} [m3/year]"
)

Total yearly transport: -52221 [m3/year]

q11 = cd.QuestionFactory(questions["Q7-11"]).serve()
q12 = cd.QuestionFactory(questions["Q7-12"]).serve()

pn.Column(q11, q12)

You will obtain more insight in how to use such a transport curve is used in a coastline model in the B1 Coastal Engineering module.

You have reached the end of this notebook.