1.3. Wave-averaged Equations


Activate wave-current interactions


Activate current effect on waves (2-way interaction)


Analytical (constant) wave parameters (Hs,Tp,Dir)


Activate wave forcing from offline model/data


Activate CROCO’s monochromatic (WKB) model


Activate coupling with spectral wave model (WW3)


Activate bottom friction for WKB model and WAVE_STREAMING


Activate bottom streaming (needs WAVE_FRICTION)


Activate Stokes drift

Preselected options:


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


Activate WKB wave model


Activate wave rollers


Activate bottom friction


Activate additional diffusion to wavenumber field


Active current effect on waves


Activate space filter on ubar,vbar,zeta for CEW


Activate time filter on ubar,vbar,zeta for CEW


Activate wave ramp


Read boundary data from croco.in


Offshore wave forcing at western boundary


Offshore wave forcing at eastern boundary

Preselected options:

# ifdef MRL_CEW
#  undef  WKB_KZ_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


Activate surface breaking acceleration


Activate Church & Thorton (1993) breaking acceleration (default)


Activate Thorton & Guza (1983, 1986)

Preselected options:

# define WAVE_BREAK_CT93
# undef  WAVE_BREAK_TG86

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} \]


Marchesiello, P.; Benshila, R.; Almar, R.; Uchiyama, Y.; McWilliams, J.C., and Shchepetkin, A., 2015. On tridimensional rip current modeling. Ocean Model., 96, 36-48.

McWilliams, J.C., Restrepo, J.M., and Lane, M.R., 2004. An asymptotic theory for the interaction of waves and currents in coastal waters. J. Fluid Mech., 511, 135–178.

Thornton, E.B. & R.T. Guza, 1983: Transformation of wave height distribution, J. Geophys. Res. 88, 5925-5938.

Uchiyama, Y., McWilliams, J., Shchepetkin, A., 2010. Wave-current interaction in an oceanic circulation model with a vortex-force formalism: application to the surf zone. Ocean Modell. 34, 16–35.

Weir, B., Uchiyama, Y., Lane, E. M., Restrepo, J. M., & McWilliams, J. C. (2011). A vortex force analysis of the interaction of rip currents and surface gravity waves. Journal of Geophysical Research: Oceans, 116(C5).