BBL

madsen

madsen(ub_r, omega_r, uc_r, phi_c, phi_wr, z_r, kN, max_iter=10, tol=0.0001)

Grant and Madsen (1979) combined wave-current bottom boundary layer model as implemented in Madsen (1994).

Parameters:

Name	Type	Description	Default
ub_r	float	Wave orbital velocity measured at height z_r (m/s)	required
omega_r	float	Radian wave frequency measured at height z_r (rad/s)	required
uc_r	float	Mean current velocity measured at height z_r (m/s)	required
phi_c	float	Angle of propagation of the mean current (degrees)	required
phi_wr	float	Angle of propagation of the surface waves (degrees)	required
z_r	float	reference height for wave and current measurements (m)	required
kN	float	Sand grain (Nikuradse) roughness height (m)	required
max_iter	int	Maximum number of iterations, by default 10	10
tol	float	Convergence tolerance, by default 1e-4 relative	0.0001

Returns:

Type	Description
dict	Dictionary containing current friction velocity ustar_c, maximum wave friction velocity ustar_wm, combined wave-current friction velocity ustar_wc, and combined wave current friction factor f_wc

References

Grant, W. D., & Madsen, O. S. (1979). Combined wave and current interaction with a rough bottom. Journal of Geophysical Research: Oceans, 84(C4), 1797-1808

Madsen, O. S. (1994). Spectral wave-current bottom boundary layer flows. In Coastal engineering 1994 (pp. 384-398).

Source code in src/pytoast/boundaries/bbl.py

def madsen(ub_r, omega_r, uc_r, phi_c, phi_wr, z_r, kN, max_iter=10, tol=1e-4):
    """
    Grant and Madsen (1979) combined wave-current bottom boundary layer model as implemented in Madsen (1994).

    Parameters
    ----------
    ub_r : float
        Wave orbital velocity measured at height z_r  (m/s)
    omega_r : float
        Radian wave frequency measured at height z_r (rad/s)
    uc_r : float
        Mean current velocity measured at height z_r (m/s)
    phi_c : float
        Angle of propagation of the mean current (degrees)
    phi_wr : float
        Angle of propagation of the surface waves (degrees)
    z_r : float
        reference height for wave and current measurements (m)
    kN : float
        Sand grain (Nikuradse) roughness height (m)
    max_iter : int, optional
        Maximum number of iterations, by default 10
    tol : float, optional
        Convergence tolerance, by default 1e-4 relative

    Returns
    -------
    dict
        Dictionary containing current friction velocity `ustar_c`, maximum wave friction velocity `ustar_wm`,
        combined wave-current friction velocity `ustar_wc`, and combined wave current friction factor `f_wc`

    References
    ----------
    Grant, W. D., & Madsen, O. S. (1979). Combined wave and current interaction with a rough bottom. Journal of
        Geophysical Research: Oceans, 84(C4), 1797-1808

    Madsen, O. S. (1994). Spectral wave-current bottom boundary layer flows. In Coastal engineering 1994 (pp. 384-398).



    """
    phi_wc = np.deg2rad(min_angle(phi_c - phi_wr))
    z0 = kN / 30

    num_iter = 0
    delta_fwc = np.inf
    c_mu = 1.0  # c_mu when mu = tau_c / tau_wr = 0
    while (num_iter < max_iter) and (delta_fwc > tol):
        relative_roughness = (c_mu * ub_r) / (kN * omega_r)

        # Wave-current friction factor
        f_wc = _fwc_m94(relative_roughness, c_mu)

        # All the friction velocities
        ustar_wm = np.sqrt(0.5 * f_wc * ub_r**2)
        ustar_r = np.sqrt(c_mu * ustar_wm**2)

        if relative_roughness > 8:
            delta_wc = 2 * kappa * ustar_r / omega_r
        else:
            delta_wc = kN

        ustar_c = (
            (ustar_r / 2)
            * (np.log(z_r / delta_wc) / np.log(delta_wc / z0))
            * (-1 + np.sqrt(1 + (4 * kappa * np.log(delta_wc / z0) / np.log(z_r / delta_wc) ** 2) * (uc_r / ustar_r)))
        )

        mu = (ustar_c / ustar_wm) ** 2
        c_mu = np.sqrt(1 + 2 * mu * np.abs(np.cos(phi_wc)) + mu**2)
        f_wc_new = _fwc_m94(relative_roughness, c_mu)
        delta_fwc = np.abs(f_wc - f_wc_new) / f_wc
        num_iter += 1

    out = {
        "ustar_c": ustar_c,
        "ustar_wm": ustar_wm,
        "ustar_wc": ustar_r,
        "f_wc": f_wc_new,
    }
    return out

styles

styles(ub, ab, ur, zr, deg, d_median, s=_S, nu=_NU, g=_G, beta=_BETA, Alpha=_ALPHA, Con=_CON, kbr_def=_KBR_DEF, max_iter=50, tol=0.0001)

Styles et al. (2017) combined wave-current bottom boundary layer model. This is a python port of the Matlab source code (https://cirp.usace.army.mil/products/bblm.php). It has been edited to remove unnecessary input parameters and some variable names have been changed for clarity. The custom bisection and secant solvers have also been replaced with scipy.optimize.brentq and scipy.optimize.newton, respectively. The output on the provided test data is within an acceptable tolerance despite these changes (see py-tests/). Finally, all inputs and outputs use SI units (m) rather than CGS units (centimeters, grams, seconds) like the source.

Parameters:

Name	Type	Description	Default
ub	float	Bottom wave orbital velocity amplitude (m/s)	required
ab	float	Bottom wave excursion amplitude (m)	required
ur	float	Mean current speed at height zr (m/s)	required
zr	float	Height above bed where ur is measured (m)	required
deg	float	Angle between wave propagation and current direction (degrees)	required
d_median	float	Median grain diameter (m)	required
s	float	Relative sediment density (default 2.65)	_S
nu	float	Kinematic viscosity (m^2/s; default 1.19e-6 for seawater at 15 C)	_NU
g	float	Gravitational acceleration (m/s^2; default 9.81)	_G
beta	float	Closure constant (default 0.7)	_BETA
Alpha	float	Closure constant (default 0.3)	_ALPHA
Con	float	Closure constant (default 6.4)	_CON
kbr_def	float	Default ripple roughness used when Psi < psicr (m; default 0.03)	_KBR_DEF
max_iter	int	Maximum bisection iterations (default 50)	50
tol	float	Bisection convergence tolerance (default 1e-4)	0.0001

Returns:

Type

Description

dict with keys:

Ro : internal friction Rossby number (Ab / (z0 * ub / ustar_wc)) mu : sqrt(ub_over_ustar_wc * phi); wave/combined stress ratio parameter epsilon : ustar_c / ustar_wc z1_over_z0 : inner wave BBL height / hydraulic roughness z2_over_z0 : outer wave BBL height / hydraulic roughness zr_over_z1 : measurement height / inner BBL height zr_over_z2 : measurement height / outer BBL height kbs : suspended sediment roughness (m) kbr : ripple roughness (m) z0 : hydraulic roughness length (m) ustar_wm : maximum wave friction velocity (m/s) ustar_c : time-averaged current friction velocity (m/s) ustar_wc : combined wave-current friction velocity (m/s)

References

Styles, R., Glenn, S. M., & Brown, M. E. (2017). An optimized combined wave and current bottom boundary layer model for arbitrary bed roughness (No. ERDCCHLTR1711).

Source code in src/pytoast/boundaries/bbl.py

def styles(
    ub,
    ab,
    ur,
    zr,
    deg,
    d_median,
    s=_S,
    nu=_NU,
    g=_G,
    beta=_BETA,
    Alpha=_ALPHA,
    Con=_CON,
    kbr_def=_KBR_DEF,
    max_iter=50,
    tol=1e-4,
):
    """Styles et al. (2017) combined wave-current bottom boundary layer model. This is a python port of the Matlab
    source code (https://cirp.usace.army.mil/products/bblm.php). It has been edited to remove unnecessary input
    parameters and some variable names have been changed for clarity. The custom bisection and secant solvers have also
    been replaced with scipy.optimize.brentq and scipy.optimize.newton, respectively. The output on the provided test
    data is within an acceptable tolerance despite these changes (see py-tests/). Finally, all inputs and outputs use SI
    units (m) rather than CGS units (centimeters, grams, seconds) like the source.

    Parameters
    ----------
    ub : float
        Bottom wave orbital velocity amplitude (m/s)
    ab : float
        Bottom wave excursion amplitude (m)
    ur : float
        Mean current speed at height zr (m/s)
    zr : float
        Height above bed where ur is measured (m)
    deg : float
        Angle between wave propagation and current direction (degrees)
    d_median : float
        Median grain diameter (m)
    s : float
        Relative sediment density (default 2.65)
    nu : float
        Kinematic viscosity (m^2/s; default 1.19e-6 for seawater at 15 C)
    g : float
        Gravitational acceleration (m/s^2; default 9.81)
    beta : float
        Closure constant (default 0.7)
    Alpha : float
        Closure constant (default 0.3)
    Con : float
        Closure constant (default 6.4)
    kbr_def : float
        Default ripple roughness used when Psi < psicr (m; default 0.03)
    max_iter : int
        Maximum bisection iterations (default 50)
    tol : float
        Bisection convergence tolerance (default 1e-4)

    Returns
    -------
    dict with keys:

        * `Ro` : internal friction Rossby number (Ab / (z0 * ub / ustar_wc))
        * `mu` : sqrt(ub_over_ustar_wc * phi); wave/combined stress ratio parameter
        * `epsilon` : ustar_c / ustar_wc
        * `z1_over_z0` : inner wave BBL height / hydraulic roughness
        * `z2_over_z0` : outer wave BBL height / hydraulic roughness
        * `zr_over_z1` : measurement height / inner BBL height
        * `zr_over_z2` : measurement height / outer BBL height
        * `kbs`        : suspended sediment roughness (m)
        * `kbr`        : ripple roughness (m)
        * `z0`         : hydraulic roughness length (m)
        * `ustar_wm`   : maximum wave friction velocity (m/s)
        * `ustar_c`    : time-averaged current friction velocity (m/s)
        * `ustar_wc`   : combined wave-current friction velocity (m/s)

    References
    ----------
    Styles, R., Glenn, S. M., & Brown, M. E. (2017). An optimized combined wave and current bottom boundary layer model
        for arbitrary bed roughness (No. ERDCCHLTR1711).
    """
    # Skin friction Shields parameter (Madsen formula)
    arg_ole = ab / d_median
    fwc_skn = np.exp(5.61 * arg_ole ** (-0.109) - 7.30)
    psi_norm = (s - 1) * g * d_median
    psi = 0.5 * fwc_skn * (1.42 * ub) ** 2 / psi_norm

    # Critical Shields parameter
    star = d_median / (4 * nu) * np.sqrt(g * d_median * (s - 1))
    psicr = _shields_critical(star)

    # Ripple geometry and bottom roughness
    if psi - psicr <= 0:
        kbr = kbr_def
    else:
        chi = 4.0 * nu * ub**2 / (d_median * ((s - 1) * g * d_median) ** 1.5)
        if chi < 2:
            eta = ab * 0.30 * chi ** (-0.39)
        else:
            eta = ab * 0.45 * chi ** (-0.99)
        kbr = Con * eta

    kbs = ab * 0.0655 * (ub**2 / ((s - 1) * g * ab)) ** 1.4
    kb = d_median + kbr + kbs
    z0 = kb / 30.0

    # Derived input parameters
    omega = ub / ab
    theta = np.deg2rad(deg)
    ab_over_z0 = ab / z0
    zr_over_z0 = zr / z0
    ub_over_kappa_ur = ub / (kappa * ur)

    # Nondimensional heights -- using notation for zeta = z / L_cw from the paper rather than z0p from source code
    alpha_loc = Alpha * (1.0 + beta * kb / ab)
    zeta_1 = alpha_loc
    delta = 1.0 / np.sqrt(2.0 * zeta_1)
    mp = delta + 1j * delta

    # Bisection to find ub_over_ustar_wc

    lower_bound = 1e-6
    # Upper bound: pure-wave limit (pwave.m), uses empirical friction factor parameterization to estimate ub / ustar_wm
    t1 = -alpha_loc * kappa * zeta_1
    if ab_over_z0 < 6.25:
        ub_over_ustar_wm = 1.0 / abs(t1 * mp)
    elif ab_over_z0 < 10:
        ub_over_ustar_wm = np.exp(1.488) * ab_over_z0 ** (-0.653) * ab_over_z0 ** (0.185 * np.log(ab_over_z0))
    elif ab_over_z0 < 100:
        ub_over_ustar_wm = np.exp(0.4599) * ab_over_z0**0.1977 * ab_over_z0 ** (0.0085 * np.log(ab_over_z0))
    else:
        ub_over_ustar_wm = np.exp(0.13996) * ab_over_z0**0.3539 * ab_over_z0 ** (-0.0106 * np.log(ab_over_z0))

    upper_bound = _pwave(ab_over_z0, ub_over_ustar_wm, zeta_1, mp)

    # Wrapper function for the bisection with a cache for the extra variables we want
    _cache: list[Any] = [None]

    def _brent_f(x):
        result = _bstress2(x, ab_over_z0, zr_over_z0, ub_over_kappa_ur, theta, alpha_loc, zeta_1, mp)
        _cache[0] = result
        return result[-1]  # fx

    brentq(_brent_f, a=upper_bound, b=lower_bound, maxiter=max_iter, xtol=tol)
    Ro, mu, epsilon, z1_over_z0, z2_over_z0, zr_over_z1, zr_over_z2, _ = _cache[0]

    ustar_wc = Ro * z0 * omega
    ustar_c = epsilon * ustar_wc
    ustar_wm = mu * ustar_wc

    out = {
        "Ro": Ro,
        "mu": mu,
        "epsilon": epsilon,
        "z1_over_z0": z1_over_z0,
        "z2_over_z0": z2_over_z0,
        "zr_over_z1": zr_over_z1,
        "zr_over_z2": zr_over_z2,
        "kbs": kbs,
        "kbr": kbr,
        "z0": z0,
        "ustar_wm": ustar_wm,
        "ustar_c": ustar_c,
        "ustar_wc": ustar_wc,
    }

    return out