4.3. Advection Schemes#

4.3.1. Lateral Momentum Advection#

Related CPP options:

UV_HADV_UP3

Activate 3rd-order upstream biased advection scheme

UV_HADV_UP5

Activate 5th-order upstream biased advection scheme

UV_HADV_C2

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

UV_HADV_C4

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

UV_HADV_C6

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

UV_HADV_WENO5

Activate WENO 5th-order advection scheme

UV_HADV_TVD

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:

TS_HADV_UP3

3rd-order upstream biased advection scheme

TS_HADV_RSUP3

Split and rotated 3rd-order upstream biased advection scheme

TS_HADV_UP5

5th-order upstream biased advection scheme

TS_HADV_RSUP5

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

TS_HADV_C4

4th-order centred advection scheme

TS_HADV_C6

Activate 6th-order centred advection scheme

TS_HADV_WENO5

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

BIO_HADV_WENO5

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:

UV_VADV_SPLINES

4th-order compact advection scheme

UV_VADV_C2

2nd-order centered advection scheme

UV_VADV_WENO5

5th-order WENOZ quasi-monotone advection scheme

UV_VADV_TVD

Total Variation Diminushing (TVD) scheme

Preselected options:

#ifdef UV_VADV_SPLINES
#elif defined UV_VADV_WENO5
#elif defined UV_VADV_C2
#elif defined UV_VADV_TVD
#else
# define UV_VADV_SPLINES
# undef  UV_VADV_WENO5
# undef  UV_VADV_C2
# undef  UV_VADV_TVD
#endif

4.3.4. Vertical Tracer advection#

Related CPP options:

TS_VADV_SPLINES

4th-order compact advection scheme

TS_VADV_AKIMA

4th-order centered advection scheme with harmonic averaging

TS_VADV_C2

2nd-order centered advection scheme

TS_VADV_WENO5

5th-order WENOZ quasi-monotone advection scheme

Preselected options:

#ifdef TS_VADV_SPLINES
#elif defined TS_VADV_AKIMA
#elif defined TS_VADV_WENO5
#elif defined TS_VADV_C2
#else
# undef  TS_VADV_SPLINES
# define TS_VADV_AKIMA
# undef  TS_VADV_WENO5
# undef  TS_VADV_C2
#endif

4.3.5. Adaptively implicit vertical advection#

Related CPP options:

VADV_ADAPT_IMP

Activate adapative, Courant number dependent implicit advection scheme

VADV_ADAPT_PRED

Adaptive treatment at both predictor and corrector steps

Preselected options:

#ifdef VADV_ADAPT_IMP
# undef  VADV_ADAPT_PRED
# define UV_VADV_SPLINES
# undef  UV_VADV_C2
# undef  UV_VADV_WENO5
# undef  UV_VADV_TVD
#endif
#ifdef VADV_ADAPT_IMP
# define  TS_VADV_SPLINES
# undef   TS_VADV_AKIMA
# undef   TS_VADV_WENO5
# undef   TS_VADV_C2
#endif

4.3.6. Numerical details on advection schemes#

../_images/grid_adv.png

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

4.3.6.1. 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}
../_images/linearADV.png

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

4.3.6.2. 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).

4.3.6.3. 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).

../_images/linearADV_Co4.png

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

4.3.6.4. 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.

4.3.6.5. Total variation bounded scheme (WENO5)#

../_images/weno5z.png

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).

4.3.6.6. Total variation diminishing scheme#

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