4.3. Advection Schemes#

4.3.1. Lateral Momentum Advection#

Related CPP options:


Activate 3rd-order upstream biased advection scheme


Activate 5th-order upstream biased advection scheme


Activate 2nd-order centred advection scheme
(should be used with explicit momentum mixing)


Activate 4th-order centred advection scheme
(should be used with explicit momentum mixing)


Activate 6th-order centred advection scheme
(should be used with explicit momentum mixing)


Activate WENO 5th-order advection scheme


Activate Total Variation Diminushing scheme

Preselected options:

# define UV_HADV_UP3
# undef  UV_HADV_UP5
# undef  UV_HADV_C2
# undef  UV_HADV_C4
# undef  UV_HADV_C6
# undef  UV_HADV_WENO5
# undef  UV_HADV_TVD

These options are set in set_global_definitions.h as the default UV_HADV_UP3 is the only one recommended for standard users.

4.3.2. Lateral Tracer advection#

Related CPP options:


3rd-order upstream biased advection scheme


Split and rotated 3rd-order upstream biased advection scheme


5th-order upstream biased advection scheme


Split and rotated 5th-order upstream biased advection scheme
with reduced dispersion/diffustion


4th-order centred advection scheme


Activate 6th-order centred advection scheme


5th-order WENOZ quasi-monotonic advection scheme for all tracers


5th-order WENOZ quasi-monotone advection scheme for
passive tracers (including biology and sediment tracers)

Preselected options:

# undef  TS_HADV_UP3
# define TS_HADV_RSUP3
# undef  TS_HADV_UP5
# undef  TS_HADV_RSUP5
# undef  TS_HADV_C4
# undef  TS_HADV_C6
# undef  TS_HADV_WENO5
# if defined PASSIVE_TRACER || defined BIOLOGY || defined SEDIMENT
#  define BIO_HADV_WENO5
# endif

TS_HADV_RSUP3 is recommended for realistic applications with variable bottom topography as it strongly reduces diapycnal mixing. It splits the UP3 scheme into 4th-order centered advection and rotated bilaplacian diffusion with grid-dependent diffusivity. It calls for CPP options in set_global_definitions.h for the explicit treatment of bilaplacian diffusion (see below). TS_HADV_RSUP3 is expensive in terms of computational cost and requires more than 30 sigma levels to perform properly. Therefore, for small domains dominated by open boundary fluxes, TS_HADV_UP5 may present a cheaper alternative and good compromise. TS_HADV_RSUP5 is still experimental but allows a decrease in numerical diffusivity compared to TS_HADV_RSUP3 by using 6th order rather than 4th-order centered advection (it resembles in spirit a split-rotated UP5 scheme but the use of bilaplacian rather than trilaplacian diffusion keeps it 3rd order). TS_HADV_C4 has no implicit diffusion and is thus accompanied by rotated Smagorinsky diffusion defined in set_global_definitions.h; it is not recommended for usual applications. For RSUP family, by default the diffusive part is oriented along geopotential.

4.3.3. Vertical Momentum advection#

Related CPP options:


4th-order compact advection scheme


2nd-order centered advection scheme


5th-order WENOZ quasi-monotone advection scheme


Total Variation Diminushing (TVD) scheme

Preselected options:

#elif defined UV_VADV_WENO5
#elif defined UV_VADV_C2
#elif defined UV_VADV_TVD
# undef  UV_VADV_WENO5
# undef  UV_VADV_C2
# undef  UV_VADV_TVD

4.3.4. Vertical Tracer advection#

Related CPP options:


4th-order compact advection scheme


4th-order centered advection scheme with harmonic averaging


2nd-order centered advection scheme


5th-order WENOZ quasi-monotone advection scheme

Preselected options:

#elif defined TS_VADV_AKIMA
#elif defined TS_VADV_WENO5
#elif defined TS_VADV_C2
# define TS_VADV_AKIMA
# undef  TS_VADV_WENO5
# undef  TS_VADV_C2

4.3.5. Adaptively implicit vertical advection#

Related CPP options:


Activate adapative, Courant number dependent implicit advection scheme


Adaptive treatment at both predictor and corrector steps

Preselected options:

# undef  UV_VADV_C2
# undef  UV_VADV_WENO5
# undef  UV_VADV_TVD
# undef   TS_VADV_AKIMA
# undef   TS_VADV_WENO5
# undef   TS_VADV_C2

4.3.6. Numerical details on advection schemes#


Fig: variable location on an Arakawa C-grid. Tracer values are cell centered while velocities are defined on interfaces.#

\[\left. \partial_x( u q ) \right|_{x=x_i} = \frac{1}{\varDelta x_i} \left\{ u_{i+\frac{1}{2}} \widetilde{q}_{i+\frac{1}{2}} - u_{i-\frac{1}{2}} \widetilde{q}_{i-\frac{1}{2}} \right\}\] Linear advection schemes#

\begin{eqnarray} \widetilde{q}_{i-\frac{1}{2}}^{\rm C2} &=& \frac{q_i + q_{i-1}}{2} \\ \widetilde{q}_{i-\frac{1}{2}}^{\rm C4} &=& \left(\frac{7}{6}\right) \widetilde{q}_{i-\frac{1}{2}}^{\rm C2} - \left(\frac{1}{12}\right) \left(q_{i+1}+q_{i-2}\right) \\ \widetilde{q}_{i-\frac{1}{2}}^{\rm UP3} &=& \widetilde{q}_{i-\frac{1}{2}}^{\rm C4} + {\rm sign}\left( \frac{1}{12} , u_{i-\frac{1}{2}} \right) \left(q_{i+1} - 3 q_i + 3 q_{i-1} - q_{i-2}\right) \\ \widetilde{q}_{i-\frac{1}{2}}^{\rm C6} &=& \left(\frac{8}{5}\right) \widetilde{q}_{i-\frac{1}{2}}^{\rm C4} - \left(\frac{19}{60}\right) \widetilde{q}_{i-\frac{1}{2}}^{\rm C2} + \left(\frac{1}{60}\right) \left(q_{i+2}+q_{i-3}\right) \\ \widetilde{q}_{i-\frac{1}{2}}^{\rm UP5} &=& \widetilde{q}_{i-\frac{1}{2}}^{\rm C6} - {\rm sign}\left( \frac{1}{60} , u_{i-\frac{1}{2}}\right) \left( q_{i+2} - 5 q_{i+1} + 10 q_i - 10 q_{i-1} + 5 q_{i-2} - q_{i-3} \right) \end{eqnarray}

Fig: amplification errors (left) and phase errors (right) for linear advection of order 2 to 6.# Split upwind schemes#

Because odd-ordered advection schemes can be formulated as the sum of the next higher-order (centered) advection scheme with a dissipation term it is possible to split the purely centered and dissipative parts of UP3 and UP5 schemes. In this case the centered part is treated within the predictor-corrector framework while the flow-dependent dissipative part is treated with a one-step Euler scheme. Such splitting has two advantages:

  1. It allows better stability for SUP3 and SUP5 schemes comapared to UP3 and UP5 schemes.

  2. Isolating the dissipative part allows to rotate it in the neutral direction to reduce spurious diapycnal mixing (RSUP3 scheme). Splines reconstruction and Akima 4th-order schemes#

Similar to a 4th-order compact scheme, the interfacial values for the splines reconstruction scheme are obtained as a solution of a tridiagonal problem

\[{\rm Hz}_{k+1} \widetilde{q}_{k-\frac{1}{2}} + 2 ({\rm Hz}_k + {\rm Hz}_{k+1}) \widetilde{q}_{k+\frac{1}{2}} + {\rm Hz}_k \widetilde{q}_{k+\frac{3}{2}} = 3 ( {\rm Hz}_k \overline{q}_{k+1} + {\rm Hz}_{k+1} \overline{q}_{k} )\]

where \(\overline{q}_{k}\) values should be understood in a finite-volume sense (i.e. as an average over a control volume).


Fig: amplification errors (left) and phase errors (right) for linear advection of order 5 and 6 and for Splines reconstruction.#

The AKIMA scheme corresponds to a 4th-order accurate scheme where an harmonic averaging of the slopes is used instead of the algebraic average used for a standard C4 scheme

\[\begin{split}\widetilde{q}_{k+\frac{1}{2}} = \frac{q_{k+1}+q_k}{2} - \frac{\overline{\delta q}_{k+1} - \overline{\delta q}_{k}}{6} \qquad \overline{\delta q}_{k} = \left\{ \begin{array}{cc} \displaystyle 2 \frac{\delta q_{k+\frac{1}{2}} \delta q_{k-\frac{1}{2}}}{\delta q_{k+\frac{1}{2}}+\delta q_{k-\frac{1}{2}}}, & \mbox{if}\;\delta q_{k+\frac{1}{2}} \delta q_{k-\frac{1}{2}} > 0 \\ 0, & \mbox{otherwise}\end{array} \right.\end{split}\] Adaptively implicit vertical advection#

Idea: the vertical velocity \(\Omega\) is split between an explicit and implicit contribution depending on the local Courant number

\[\begin{split}\Omega = \mathbf{\Omega^{(e)}} + \mathbf{\Omega^{(i)}}, \qquad \mathbf{\Omega^{(e)}} = \frac{\Omega}{f(\alpha_{\rm adv}^z,\alpha_{\max})}, \quad f(\alpha_{\rm adv}^z,\alpha_{\max}) = \left\{ \begin{array}{ll} 1, & \alpha_{\rm adv}^z \le \alpha_{\max} \\ \alpha / \alpha_{\max}, & \alpha_{\rm adv}^z > \alpha_{\max}\end{array} \right.\end{split}\]
  • \(\Omega^{(e)}\) is integrated with an explicit scheme with CFL \(\alpha_{\max}\).

  • \(\Omega^{(i)}\) is integrated with an implicit upwind Euler scheme.

  • \(f(\alpha_{\rm adv}^z,\alpha_{\max})\) is a function responsible for the splitting of \(\Omega\) between an explicit and an implicit part.

This approach has the advantage to render vertical advection unconditionally stable and to maintain good accuracy in locations with small Courant numbers. The current implementation is based on the SPLINES scheme for the explicit part. Total variation bounded scheme (WENO5)#


Fig: different stencils used to evaluate the interfacial value \(\widetilde{q}_{k+\frac{1}{2}}\) with WENO5 scheme#

Nonlinear weighting between 3 evaluations of interfacial values based on 3 different stencils

\[\widetilde{q}_{k-\frac{1}{2}} = w_0 \widetilde{q}_{k-\frac{1}{2}}^{(0)} + w_1 \widetilde{q}_{k-\frac{1}{2}}^{(1)} + w_2 \widetilde{q}_{k-\frac{1}{2}}^{(2)}\]

where the weights are subject to the following constraints:

  1. Convexity \(\sum_{j=0}^2 w_j = 1\).

  2. ENO property (Essentially non oscillatory).

  3. 5th-order if \(q(x)\) is smooth.

The resulting scheme is not monotonicity-preserving but instead it is Total Variation Bounded (TVB). Total variation diminishing scheme# Upwinding of nonlinear terms#

In CROCO the nonlinear advection terms are formulated as in Lilly (1965) :

\begin{eqnarray} \partial_t ( {\rm Hz} u ) &+& \partial_x\left( ( {\rm Hz}\; u) u \right) + \partial_y\left( ( {\rm Hz} \; v)\; u \right) + \; ... \\ \partial_t ( {\rm Hz} v ) &+& \partial_x\left( ( {\rm Hz}\; u) v \right) + \partial_y\left( ( {\rm Hz} \; v)\; v \right) + \; ... \end{eqnarray}

which are discretised with third order accuracy as

\begin{eqnarray} \left( \widetilde{( {\rm Hz}\; u) u} \right)_{i,j} &=& ( \widetilde{{\rm Hz}\; u})_{i,j}^{\rm C4} \; \widetilde{u}_{i,j}^{\rm UP3} \\ \left( \widetilde{( {\rm Hz}\; v) u} \right)_{i+\frac{1}{2},j+\frac{1}{2}} &=& ( \widetilde{{\rm Hz}\; v})_{i+\frac{1}{2},j+\frac{1}{2}}^{\rm C4} \; \widetilde{u}_{i+\frac{1}{2},j+\frac{1}{2}}^{\rm UP3} \end{eqnarray}

where the direction for upwinding is selected considering

\[u^{\rm upw}_{i,j} = u_{i+\frac{1}{2},j}+u_{i-\frac{1}{2},j}, \qquad v^{\rm upw}_{i+\frac{1}{2},j+\frac{1}{2}} = ({\rm Hz}\; v)_{i,j+\frac{1}{2}} + ({\rm Hz}\; v)_{i+1,j+\frac{1}{2}}\]