Bottom Boundary Layer model
===========================

**Related CPP options:**

=============== =================================================================================
BBL             Activate bottom boundary layer parametrization
ANA_WWAVE       Set analytical (constant) wave forcing (hs,Tp,Dir).
ANA_BSEDIM      Set analytical bed parameters (if SEDIMENT is undefined)
Z0_BL           Compute bedload roughness for ripple predictor and sediment purposes
Z0_RIP          Determine bedform roughness ripple height and ripple length for sandy bed
Z0_BIO          Determine (biogenic) bedform roughness ripple height and ripple length for silty beds
=============== =================================================================================

*Preselected options:*
::

 #ifdef BBL
 # ifdef OW_COUPLING
 # elif defined WAVE_OFFLINE
 # elif defined WKB_WWAVE
 # else
 #  define ANA_WWAVE
 # endif
 # ifdef SEDIMENT
 #  undef  ANA_BSEDIM
 # else
 #  define ANA_BSEDIM
 # endif
 # ifdef SEDIMENT
 #  define Z0_BL
 # else
 #  undef  Z0_BL
 # endif
 # ifdef Z0_BL
 #  define Z0_RIP
 # endif
 # undef  Z0_BIO
 #endif


**DESCRIPTION**

Reynolds stresses, production and dissipation of turbulent kinetic energy, and 
gradients in velocity and suspended-sediment concentrations vary over short 
vertical distances, especially near the bed, and can be difficult to resolve with 
the vertical grid spacing used in regional-scale applications. CROCO provides 
algorithms to parameterize some of these subgrid-scale processes in the water column 
and in the bottom boundary layer (BBL). Treatment of the BBL is important for the 
circulation model solution because it determines the stress exerted on the flow by 
the bottom, which enters the Reynolds-averaged Navier-Stokes equations as a boundary 
conditions for momentum in the x and y directions:

.. math::

   K_m \frac{\partial u}{\partial s} = \tau_{bx} \\
   K_m \frac{\partial v}{\partial s} = \tau_{by}

Determination of the BBL is even more important for the sediment-transport formulations 
because bottom stress determines the transport rate for bedload and the resuspension 
rate for suspended sediment.

CROCO implements either of two methods for representing BBL processes: (1) simple 
drag-coefficient expressions or (2) more complex formulations that represent the 
interactions of wave and currents over a moveable bed. The drag-coefficient methods 
implement formulae for linear bottom friction, quadratic bottom friction, or a 
logarithmic profile. The other, more complex  wave-current BBL model is described 
by :cite:t:`blaas_sediment-transport_2007` with an example of its use on the Southern 
California continental shelf. The method uses efficient wave-current BBL computations 
developed by :cite:t:`soulsby1995bed` in combination with sediment and bedform roughness 
estimates of :cite:t:`grant_movable_1982`, :cite:t:`nielsen_suspended_1986` 
and :cite:t:`li_sedtrans96_2001`. 

**Linear/quadratic drag**

The linear and/or quadratic drag-coefficient methods depend only on velocity components 
u and v in the bottom grid cell and constant, spatially-uniform coefficients :math:`\gamma_1` 
and :math:`\gamma_2` specified as input:

.. math::

    \tau_{bx} = (\gamma_1 + \gamma_2 \sqrt{u^2+v^2}) ~u  \\
    \tau_{by} = (\gamma_1 + \gamma_2 \sqrt{u^2+v^2}) ~v

where :math:`\gamma_1` is the linear drag coefficient and :math:`\gamma_2` is the 
quadratic drag coefficient. The user can choose between linear or quadratic drag by 
setting one of these coefficients to zero. The bottom stresses computed from these 
formulae depend on the elevation of u and v (computed at the vertical mid-elevation of the 
bottom computational cell). Therefore, in this s-coordinate model, the same drag coefficient 
will be imposed throughout the domain even though the vertical location of the velocity 
is different.

**Logarithmic drag** (with roughness length :math:`z_0`)


To prevent this problem, the quadratic drag :math:`\gamma_2` can be computed assuming 
that flow in the BBL has the classic vertical logarithmic profile defined by a shear 
velocity :math:`u_*` and bottom roughness length :math:`z_0` (m) as:

.. math::

   \left | u \right | = \frac{u_*}{\kappa} \ln \left ( \frac{z}{z_0} \right ) 

where  :math:`\left | u \right | =\sqrt{u^2+v^2}`, friction velocity  :math:`u_*=\sqrt{\tau_b}`, 
z is the elevation above the bottom (vertical mid-elevation point of the bottom cell), 
:math:`\kappa=0.41` is von Kármán’s constant. :math:`z_0` is an empirical parameter. 
It can be constant (default) or spatially varying. Kinematic stresses are calculated as`

.. math::

    \tau_{bx} = \frac{\kappa^2}{ \ln^2 \left ( z/z_0 \right )} \sqrt{u^2+v^2} u  \\
    \tau_{by} = \frac{\kappa^2}{ \ln^2 \left ( z/z_0 \right )} \sqrt{u^2+v^2} v

The advantage of this approach is that the velocity and the vertical elevation of that 
velocity are used in the equation. Because the vertical elevation of the velocity in the 
bottom computational cell will vary spatially and temporally, the inclusion of the 
elevation provides a more consistent formulation.

**Combined wave-current drag (BBL)**

To provide a more physically relevant value of :math:`z_0`, especially when considering 
waves and mobile sediments, a more complex formulation is available (BBL). 

The short (order 10-s) oscillatory shear of wave-induced motions in a thin (a few cm) 
wave-boundary layer produces turbulence and generates large instantaneous shear stresses. 
The turbulence enhances momentum transfer, effectively increasing the bottom-flow coupling 
and the frictional drag exerted on the wave-averaged flow. The large instantaneous shear 
stresses often dominate sediment resuspension and enhance bedload transport. Sediment 
transport can remold the bed into ripples and other bedforms, which present roughness 
elements to the flow. Bedload transport can also induce drag on the flow, because momentum 
is transferred to particles as they are removed from the bed and accelerated by the flow. 
Resuspended sediments can cause sediment-induced stratification and, at high concentrations, 
change the effective viscosity of the fluid.

The BBL parameterization implemented in CROCO requires inputs of velocities u and v at 
reference elevation z, representative wave-orbital velocity amplitude :math:`u_b`, wave period 
T, and wave propagation direction :math:`\theta` (degrees, clockwise from north). The wave 
parameters may be the output of a wave model such as WKB or WW3 or simpler calculations based 
on specified surface wave parameters. Additionally the BBL models require bottom sediment 
characteristics (median grain diameter 
:math:`D_{50}`, mean sediment density :math:`\rho_s`, and representative settling 
velocity :math:`w_s`); these are constant (ANA_BSEDIM) or based on the composition of 
the uppermost active layer of the bed sediment during the previous time step if the sediment 
model is used. 

The wave-averaged, combined wave–current bottom stress is expressed as function of :math:`\tau_w` 
and :math:`\tau_c` (i.e., the stress due to waves in the absence of currents and due to currents 
in the absence of waves, respectively) according to :cite:t:`soulsby1995bed`:

.. math::

   \bar{\tau}_{wc}=\tau_c \left ( 1+1.2 \left ( \frac{\tau_w}{\tau_w + \tau_c} \right )^{3.2}  \right ) 

The maximum wave–current shear stress within a wave cycle is obtained by adding :math:`\bar{\tau}_{wc}` 
and :math:`\tau_w` (with :math:`\phi` the angle between current and waves):

.. math::

   \tau_{wc}=  \left ( (\bar{\tau}_{wc}+\tau_w \cos{\phi})^2 + (\tau_w \sin{\phi})^2 \right )^{1/2}

The stresses :math:`\tau_c` and :math:`\tau_w` are determined using:

.. math::

   \tau_{c} = \frac{\kappa^2}{ \ln^2 \left ( z/z_0 \right )}  { \left | u \right | }^2 \\
   \tau_{w} = 0.5 \rho f_w u_b^2

:math:`u_b`, the bottom orbital velocity, is determined from the significant wave height :math:`H_s` 
and peak frequency :math:`\omega_p` using the Airy wave theory:

.. math::

   u_b=\omega_p \frac{H_s}{2 \sinh{kh}}

with h the local depth and k the local wave number from the dispersion relation. The wave-friction 
factor :math:`f_w` is, according to :cite:t:`soulsby1995bed`:

.. math::

   f_w=1.39(u_b/\omega_p z_0)^{-0.52}

The wave–current interaction in the BBL is taken into account only if :math:`u_b>1` cm/s; otherwise, 
current-only conditions apply. 

**Shear stress for sediment resuspension and roughness length due to bed form**

To determine the shear stress relevant for sediment resuspension and the roughness length due to 
bed forms, we follow the concept of :cite:t:`li_sedtrans96_2001` briefly summarized here. First, 
the maximum wave–current skin friction :math:`\tau_{s}` is computed from the equations above, 
using the Nikuradse roughness :math:`z_0= D_{50}/12`. 

A bed-load layer develops as soon as the maximum wave–current skin friction :math:`\tau_{s}` 
exceeds the critical stress :math:`\tau_{cr}`. This layer affects the stress effective for ripple 
formation and sediment resuspension. Subsequently, for sandy locations, ripple height and length 
are computed, leading to a  space- and time-dependent ripple roughness length :math:`z_0= z_{rip}`, 
which is used to compute the drag on the flow (instead of a constant value when BBL is not 
activated). This drag provides boundary conditions to the momentum and turbulence equations 
(KPP or GLS).