1.3. Wave-averaged Equations

MRL_WCI	Activate wave-current interactions
MRL_CEW	Activate current effect on waves (2-way interaction)
ANA_WWAVE	Analytical (constant) wave parameters (Hs,Tp,Dir)
WAVE_OFFLINE	Activate wave forcing from offline model/data
WKB_WWAVE	Activate CROCO’s monochromatic (WKB) model
OW_COUPLING	Activate coupling with spectral wave model (WW3)
WAVE_FRICTION	Activate bottom friction for WKB model and WAVE_STREAMING
WAVE_STREAMING	Activate bottom streaming (needs WAVE_FRICTION)
STOKES_DRIFT	Activate Stokes drift

Preselected options:

# define STOKES_DRIFT

A vortex-force formalism for the interaction of surface gravity waves and currents is implemented in CROCO [Marchesiello et al., 2015, Uchiyama et al., 2010]. Eulerian wave-averaged current equations for mass, momentum, and tracers are included based on an asymptotic theory by McWilliams et al. [2004] plus non-conservative wave effects due to wave breaking, associated surface roller waves, bottom streaming, and wave-enhanced vertical mixing and bottom drag especially for coastal and nearshore applications. The wave information is provided by either a spectrum-peak WKB wave-refraction model that includes the effect of currents on waves, or, alternatively, a spectrum-resolving wave model (e.g., WAVEWATCH3) can be used. In nearshore applications, the currents’ cross-shore and vertical structure is shaped by the wave effects of near-surface breaker acceleration, vertical component of vortex force, and wave-enhanced pressure force and bottom drag.

1.3.1. Equations in Cartesian coordinates

In the Eulerian wave-averaged current equations, terms for the wave effect on currents (WEC) are added to the primitive equations. Three new variables are defined:

\[ \begin{align}\begin{aligned}\xi^c & = \xi + \hat{\xi}\\\phi^c & = \phi + \hat{\phi}\\\vec{\textbf v}_L & = \vec{\textbf v} + \vec{\textbf v}_S\end{aligned}\end{align} \]

where \(\xi^c\) is a composite sea level, \(\phi^c\) absorbs the Bernoulli head \(\hat{\phi}\), \(\vec{\textbf v_L}\) is the wave-averaged Lagrangian velocity, sum of Eulerian velocity and Stokes drift \(\vec{\textbf v_S}\). The 3D Stokes velocity is non-divergent and defined for a monochromatic wave field (amplitude A, wavenumber vector \(\vec{\textbf k}=(k_x,k_y)\), and frequency \(\sigma\)) by:

\[ \begin{align}\begin{aligned}u_S & = \frac{A^2 \sigma}{2 \sinh^2 \left ( k D \right ) } \cosh \left ( 2 k (z+h) \right ) k_x\\v_S & = \frac{A^2 \sigma}{2 \sinh^2 \left ( k D \right ) } \cosh \left ( 2 k (z+h) \right ) k_y\\w_S & = - \int_{-h}^{z} \left ( \frac{\partial u_S}{\partial x} + \frac{\partial v_S}{\partial y}\right ) ~dz'\end{aligned}\end{align} \]

Where \(D=h+\xi^c\). The quasi-static sea level and Bernouilli head are:

\[ \begin{align}\begin{aligned}\hat{\xi} & = - \frac{A^2 k}{2 \sinh \left ( 2 k D \right ) }\\\hat{\phi} & = \frac{A^2 \sigma}{4 k \sinh^2 \left ( k D \right ) } \int_{-h}^{z} \frac{\partial^2 \vec{\textbf k}.\vec{\textbf v} } {\partial z'^2} \sinh \left ( 2 k (z-z') \right ) ~dz'\end{aligned}\end{align} \]

The primitive equations become (after re-organizing advection and vortex force terms):

\[ \begin{align}\begin{aligned}\frac{\partial u}{\partial t} + \vec{\bf \nabla} . \left ( \vec{\textbf v}_L u \right ) - f v_L & = - \frac{\partial \phi^c}{\partial x} + \left ( u_S \frac{\partial u}{\partial x} + v_S \frac{\partial v}{\partial x} \right ) + \mathcal{F}_u + \mathcal{D}_u + \mathcal{F^W}_u\\\frac{\partial v}{\partial t} + \vec{\bf \nabla} . \left ( \vec{\textbf v}_L v \right ) + f u_L & = - \frac{\partial \phi^c}{\partial y} + \left ( u_S \frac{\partial u}{\partial y} + v_S \frac{\partial v}{\partial y} \right ) + \mathcal{F}_v + \mathcal{D}_v + \mathcal{F^W}_v\\\frac{\partial \phi^c}{\partial z} + \frac{\rho g}{\rho_0} & = \vec{\textbf v}_S. \frac{\partial \vec{\textbf v}}{\partial z}\\ \frac{\partial C}{\partial t} + \vec{\bf \nabla} . \left ( \vec{\textbf v}_L C \right ) & = \mathcal{F}_C + \mathcal{D}_C + \mathcal{F^W}_C\\ \vec{\bf \nabla} . \vec{\textbf v}_L & = 0\\ \rho & = \rho(T,S,P)\end{aligned}\end{align} \]

The variables used are :

\(\mathcal{D}_u, \mathcal{D}_v, \mathcal{D}_C\) : diffusive terms (including wave-enhaced bottom drag and mixing)

\(\mathcal{F}_u, \mathcal{F}_v, \mathcal{F}_C\) : forcing terms

\(\mathcal{F^W}_u, \mathcal{F^W}_v, \mathcal{F^W}_C\) : wave forcing terms (bottom streaming, breaking acceleration)

\(f(x,y)\) : Traditional Coriolis parameter \(2 \Omega sin \phi\)

\(g\) : acceleration of gravity

\(\phi(x,y,z,t)\) : dynamic pressure \(\phi=P/\rho_0\), with P the total pressure

\(\rho_0+\rho(x,y,z,t)\) : total in situ density

\(u,v,w\) : the (x,y,z) components of vector velocity \(\vec{\textbf v}\)

1.3.2. Embedded wave model

WKB_WWAVE	Activate WKB wave model
WAVE_ROLLER	Activate wave rollers
WAVE_FRICTION	Activate bottom friction
WKB_ADD_DIFF	Activate additional diffusion to wave number field
MRL_CEW	Active current effect on waves
WKB_KZ_FILTER	Activate space filter on ubar, vbar, zeta for CEW
WKB_TIME_FILTER	Activate time filter on ubar, vbar, zeta for CEW
WAVE_RAMP	Activate wave ramp
ANA_BRY_WKB	Read boundary data from croco.in
WKB_OBC_WEST	Offshore wave forcing at the western boundary
WKB_OBC_EAST	Offshore wave forcing at the eastern boundary

Preselected options:

# ifdef MRL_CEW
#  undef  WKB_KZ_FILTER
#  undef  WKB_TIME_FILTER
# endif
# define WKB_ADD_DIFF
# if defined SHOREFACE || defined SANDBAR || (defined RIP && !defined BISCA)
#  define ANA_BRY_WKB
# endif

A WKB wave model for monochromatic waves is embedded in CROCO following Uchiyama et al. [2010]. It is based on the conservation of wave action \(\mathcal{A} = E/\sigma\) and wavenumber \(\textbf k\) – wave crest conservation – and is particularly suitable for nearshore beach applications, allowing refraction from bathymetry and currents (but no diffraction or reflection), with parametrizations for wave breaking and bottom drag:

\[ \begin{align}\begin{aligned}\frac{\partial \mathcal{A}}{\partial t} + \vec{\bf \nabla} . \mathcal{A} \vec{\bf c}_g = -\frac{\epsilon^w}{\sigma}\\ \frac{\partial \vec{\bf k}}{\partial t} + \vec{\bf c}_g . {\bf \nabla} \vec{\bf k} = - \vec{\bf k} . {\bf \nabla} \vec{\bf V} - \frac{k\sigma}{\sinh 2kD} {\bf \nabla} D\end{aligned}\end{align} \]

\(\vec{\bf V}\) is the depth-averaged velocity vector and \(\sigma\) is the intrinsic frequency defined by the linear dispersion relation \(\sigma^2 = gk \tanh kD\). Current effects on waves are noticeable in the groupe velocity \(c_g\) which gets two components: the doppler shift due to currents on waves and the groupe velocity of the primary carrier waves :

\[\vec{\bf c}_g = \vec{\bf V} + \frac{\sigma}{2k^2} \left ( 1+\frac{2kD}{\sinh 2kD} \right ) \vec{\textbf k}\]

The currents may need filtering before entering the wave model equations because the current field should evolve slowly with respect to waves in the asymptotic regime described by McWilliams et al. [2004]. By default, this filtering is turned off (WKB_KZ_FILTER, WKB_TIME_FILTER).

\(\epsilon^w\) is the depth-integrated rate of wave energy dissipation due to depth-induced breaking \(\epsilon^b\) (including white capping) and bottom friction \(\epsilon^{wd}\), both of which must be parameterized (in WKB, WW3 or CROCO if defined WAVE_OFFLINE):

\[\epsilon^w = \epsilon^b + \epsilon^{wd}\]

1.3.3. Breaking acceleration and bottom streaming

A formulation for \(\epsilon^{b}\) is needed in both the wave model (dissipation term) and the circulation model (acceleration term). In the wave-averaged momentum equations of the circulation model, the breaking acceleration enters as a body force through \(\mathcal{F^W}\):

\[\vec{\bf F^b} = \frac{\epsilon^b}{\rho \sigma} \vec{\bf k} ~f_b(z)\]

where \(f_b(z)\) is a normalized vertical distribution function representing vertical penetration of momentum associated with breaking waves from the surface. The penetration depth is controlled by a vertical length-scale taken as \(H_{rms}\).

The wave model can also include a roller model with dissipation \(\epsilon^r\). In this case:

\[\vec{\bf F^b} = \frac{(1-\alpha_r)\epsilon^b+\epsilon^r}{\rho \sigma} \vec{\bf k} ~f_b(z)\]

The idea is that some fraction \(\alpha_r\) of wave energy is converted into rollers that propagate toward the shoreline before dissipating, while the remaining fraction \(1-\alpha_r\) causes local dissipation (hence current acceleration). It can be useful for correcting \(\epsilon^b\) with some flexibility to depict different breaking wave and beach forms (e.g., spilling or plunging breakers, barred or plane beaches), although the parameter \(B_b\) can also be used for that. See Uchiyama et al. [2010] for the roller equation and \(\epsilon^r\) formulation.

Wave-enhanced bottom dissipation enters in the momentum equations through a combined wave-current drag formulation (see parametrizations) and bottom streaming. The latter is due to dissipation of wave energy in the wave boundary layer that causes the instantaneous, oscillatory wave bottom orbital velocities to be slightly in phase from quadrature; this causes a wave stress (bottom streaming) in the wave bottom boundary layer along the direction of wave propagation [Longuet-Higgins, 1953]. The effect of bottom streaming in momentum balance is accounted for by using the wave dissipation due to bottom friction with an upward decaying vertical distribution:

\[\vec{\bf F^{st}} = \frac{\epsilon^{wd}}{\rho \sigma} \vec{\bf k} ~f_{st}(z)\]

where \(f_{st}(z)\) is a vertical distribution function.

1.3.4. Formulation of wave energy dissipation

WAVE_SFC_BREAK	Activate surface breaking acceleration
WAVE_BREAK_CT93	Activate Church and Thornton [1993] breaking acceleration (default)
WAVE_BREAK_TG86	Activate Thornton and Guza [1983], Thornton and Guza [1986]

Preselected options:

# define WAVE_BREAK_CT93
# undef  WAVE_BREAK_TG86
# undef  WAVE_SFC_BREAK

While a few formulations for \(\epsilon^b\) are implemented in CROCO, the one by Church and Thornton [1993] is generally successful for nearshore beach applications:

\[\epsilon^b = \frac{3}{16} \sqrt{\pi} \rho g B^3_b \frac{H^3_{rms}}{D} \left \{ 1 + \tanh \left [ 8 \left ( \frac{H_{rms}}{\gamma_b D} -1 \right ) \right ] \right \} \left \{ 1 - \left [ 1 + \left ( \frac{H_{rms}}{\gamma_b D} \right )^2 \right ]^{-2.5} \right \}\]

where \(B_b\) and \(\gamma_b\) are empirical parameters related to wave breaking. \(\gamma_b\) represents the wave height-to-depth ratio for which all waves are assumed to be breaking and \(B_b\) is the fraction of foam on the face, accounting for the type of breaker. \(H_{rms}\) is the RMS wave height. For the DUCK94 experiment, Uchiyama et al. [2010] suggest \(\gamma_b=0.4\) and \(B_b = 0.8\), while for Biscarrosse Beach, Marchesiello et al. [2015] use \(\gamma_b=0.3\) and \(B_b = 1.3\) from calibration with video cameras.

For \(\epsilon^{wd}\), the dissipation caused by bottom viscous drag on the primary waves, we use a parameterization for the realistic regime of a turbulent wave boundary layer, consistent with the WKB spectrum-peak wave modeling:

\[\epsilon^{wd} = \frac{1}{2 \sqrt\pi} \rho f_w u^3_{orb}\]

where \(u_{orb}\) is the wave orbital velocity magnitude and \(f_w\) is a wave friction factor, function of roughness length \(z_0\):

\[ \begin{align}\begin{aligned}u_{orb} & = \frac{\sigma H_{rms}}{2 \sinh kD}\\f_w & = 1.39 \left ( \frac{\sigma z_0}{u_{orb}} \right ) ^{0.52}\end{aligned}\end{align} \]