5.1. Vertical mixing parametrizations

CROCO contains a variety of methods for setting the vertical viscous and diffusive coefficients. The choices range from simply choosing fixed values to the KPP and the generic lengthscale (GLS) turbulence closure schemes. See Large (1998) for a review of surface ocean mixing schemes. Many schemes have a background molecular value which is used when the turbulent processes are assumed to be small (such as in the interior).

Related CPP options:

ANA_VMIX

Analytical definition

BVF_MIXING

Brunt-Vaisaleafrequency based

LMD_MIXING

K-profile parametrisation

GLS_MIXING

Generic lengthscale parametrisation

Preselected options:

NONE : default is no mixing scheme

5.1.1. Analytical definition

Related CPP options:

ANA_VMIX

Analytical definition

Preselected options:

NONE

A profile for mixing coeefficient \(K_{m,s}(z)\) can be set in ana_vmix routine for variables Akv (viscosity) and Akt (diffusivity), which is called at each time step. In this case, background coeeficients read in croco.in can be used.

5.1.2. BVF mixing

Related CPP options:

BVF_MIXING

Brunt-Vaisala frequency based

Preselected options:

NONE

It computes diffusivity using a Brunt-Vaisala frequency based vertical mixing scheme. Viscosity is set to its background. In static unstable regime, diffusivity is enhanced.

  • If \(N^2(z) < 0\) :

    \[K_{m,s}(z) = 0.1\;{\rm m}^2\;{\rm s}^{-1}\]

  • If \(N^2(z) > 0\) :

    \[K_{m,s}(z) = 10^{-7} / \sqrt{ N^2(z) }, \qquad K_{m,s}^{\min} \le K_{m,s}(z) \le K_{m,s}^{\max}\]

Default bounds are quite restrictive :

\[K_{m,s}^{\min} = 3 \times 10^{-5}\;{\rm m}^2\;{\rm s}^{-1}, \qquad K_{m,s}^{\max} = 4 \times 10^{-4}\;{\rm m}^2\;{\rm s}^{-1}\]

5.1.3. K-profile parametrization

Large, W., J. McWilliams, and S. Doney, Oceanic vertical mixing: A review and a model with nonlocal boundary layer parameterization, Rev. Geophys., 32, 363-403, 1994.

Related CPP options:

KPP-related options :

LMD_MIXING

K-profile parametrisation

LMD_SKPP

Activate surface boundary layer KPP mixing

LMD_SKPP2005

Activate surface boundary layer KPP mixing (2005 version)

LMD_BKPP

Activate bottom boundary layer KPP mixing

LMD_BKPP2005

Activate bottom boundary layer KPP mixing (2005 version)

LMD_RIMIX

Activate shear instability interior mixing

LMD_CONVEC

Activate convection interior mixing

LMD_DDMIX

Activate double diffusion interior mixing

LMD_NONLOCAL

Activate nonlocal transport for SKPP

LMD_LANGMUIR

Activate Langmuir turbulence mixing

Preselected options:

# ifdef LMD_MIXING
#  define LMD_SKPP
#  define LMD_BKPP
#  define LMD_RIMIX
#  define LMD_CONVEC
#  undef  LMD_DDMIX
#  define LMD_NONLOCAL
#  undef  LMD_LANGMUIR
# endif

#if defined LMD_SKPP # define LMD_SKPP2005 #endif #ifdef LMD_BKPP # undef LMD_BKPP2005 #endif

Surface boundary layer

  • LMD_SKPP (Large et al, 1994)

    • Step 1 : Compute boundary layer depth \(h_{bl} (z_r \rightarrow z_N)\)

      \[{\rm Ri}_b(z) = \frac{g (z_r - z) \left( \rho(z) - \rho_r \right) / \rho_0}{|\mathbf{u}(z) - (\mathbf{u}_h)_r|^2 + V_t^2(z)}, \qquad {\rm Ri}_b(-h_{bl}) = {\rm Ri}_{cr}\]

    • Step 2 : In the stable case math::(B_f > 0) : h_{bl} = min( h_{bl}, h_{ek}, h_{mo} )

      \[h_{ek} = 0.7 u_\star / f, \qquad h_{mo} = u_{\star}^3 / (\kappa B_f).\]

    • step 3 : Compute turbulent viscosity and diffusivity

      \[K_{m,s}(z) = w_{m,s} \; h_{bl} \; G(z/h_{bl}), \qquad w_{m,s} = \kappa \; u_\star \; \psi_{m,s}( z B_f / u_\star^3 )\]

Choice of the critical Richardson number \({\rm Ri}_{cr}\) : \(\qquad {\rm Ri}_{cr} \in [ 0.15, 0.45 ]\)

  • LMD_SKPP2005 (Shchepetkin et al, 2005)

    • Criteria for \(h_{bl}\) : integral layer where production of turbulence by shear balances dissipation by the stratification

      \[{\rm Cr}(z) = \int_{z}^{\zeta} {\cal K}(z') \left\{ \left| \partial_{z'} \mathbf{u}_h \right|^2 - \frac{N^2}{{\rm Ri}_{cr}} - C_{Ek} \; f^2 \right\} dz' + \frac{V_t^2(z)}{(\zeta - z)}, \; {\rm Cr}(-h_{bl}) = 0\]

    • Consistent with the original KPP

      \[{\rm Cr}(-h_{bl}) = 0 \Rightarrow \frac{ (\zeta - z) \int_{z}^{\zeta} {\cal K}(z') N^2(z') dz'}{ (\zeta - z) \int_{z}^{\zeta} {\cal K}(z') \left\{ \left| \partial_z \mathbf{u}_h \right|^2 - C_{Ek} \; f^2 \right\} dz' + V_t^2(z)} = {\rm Ri}_{cr}\]
Advantages :
-> consistent with Ekman problem
-> tends to give deeper boundary layers : \((\zeta - z) \int_{z}^{\zeta} \left| \partial_{z'} \mathbf{u}_h \right|^2 dz' \ge |\mathbf{u}_h(z) - \mathbf{u}_h(\zeta)|^2\).

  • cpp key LMD_LANGMUIR (McWilliams & Sullivan, 2000)

Following the work of McWilliams and Sullivan (2000), we introduce in KPP an enhancement factor E to the turbulent velocity scale as a function of the turbulent Langmuir number \(La_t=\sqrt{u_\star / u_{Stokes}}\), but this function is taken as in Van Roekel et al. (2012) which gives good results in Li et al. (2016) – still assuming that Stokes drift is aligned with the surface wind stress:

\[\begin{split}w_{m,s} = \frac{\kappa \; u_\star}{\phi_{m,s}} \; E , \qquad E = \sqrt{1+0.104 La_t^{-2} + 0.034 La_t^{-4}} \\\end{split}\]

Interior scheme

\[K_{m,s}(z) = K_{m,s}^{\rm sh}(z) + K_{m,s}^{\rm iw}(z) + K_{m,s}^{\rm dd}(z)\]

  • cpp key LMD_RIMIX, RI_(H-V)SMOOTH (Large et al., 1994)

    \[{\rm Ri}_g = N^2 / \left[ ( \partial_z u )^2 + ( \partial_z v )^2 \right]\]
    \[\begin{split} K_{m,s}^{\rm sh}(z) = \left\{ \begin{array}{lrl} K_{0,c} & {\rm Ri}_g < 0 & {\rm { \leftarrow [LMD\_CONVEC]}} \\ K_0 \left[ 1 - ( \frac{{\rm Ri}_g}{{\rm Ri}_0} )^3 \right] & 0 < {\rm Ri}_g < {\rm Ri}_0 \\ 0 & {\rm Ri}_0 < {\rm Ri}_g \\ \end{array} \right.\end{split}\]
    \[K_0 = 5 \times 10^{-3} \; {\rm m}^2\;{\rm s}^{-1}, {\rm Ri}_0 = 0.7\]

  • cpp key LMD_NUW_GARGETT (Gargett & Holloway)

    \[K_{m}^{\rm iw}(z) = \frac{10^{-6}}{ \sqrt{\max(N^2(z),10^{-7})} } , \qquad K_{s}^{\rm iw}(z) = \frac{10^{-7}}{ \sqrt{\max(N^2(z),10^{-7})} }\]

  • cpp key LMD_DDMIX (cf Large et al., 1994, eqns (31))

Bottom boundary layer

  • cpp key LMD_BOTEK : Bottom Ekman layer

    \begin{eqnarray} h_{\rm Ek} &=& \min \left\{ \frac{0.3 u_{\star,b}}{| f |}, h \right\} \nonumber \\ \sigma_{k+\frac{1}{2}} &=& (z_{k+\frac{1}{2}}-h) / h_{\rm Ek} \nonumber \\ K^{\rm Ek}_{k+\frac{1}{2}} &=& \max \left\{ 4 \; \kappa \; u_{\star,b} \; h_{\rm Ek} \; \sigma ( 1 - \sigma ), K_{\min} \right\} \nonumber \\ {\rm AKv}_{k+\frac{1}{2}} &=& {\rm AKv}_{k+\frac{1}{2}} + K^{\rm Ek}_{k+\frac{1}{2}} \nonumber \\ {\rm AKt}_{k+\frac{1}{2}} &=& {\rm AKt}_{k+\frac{1}{2}} + K^{\rm Ek}_{k+\frac{1}{2}} \nonumber \end{eqnarray}

  • cpp key LMD_BKPP (Bottom KPP 1994)

Same rationale than surface KPP but this time we search for the critical value \({\rm Ri}_{\rm cr}\;(\approx 0.3)\) starting from the bottom

\[h_{\rm bbl} = \min \left( h_{\rm bbl}, \frac{0.7 u_{\star,b}}{| f |} \right) K_{m,s}(z) = \kappa \; u_{\star,b} \; h_{bbl} \; G( \sigma ), \qquad \sigma = \frac{(z-h)}{h_{bbl}}\]

5.1.4. Generic length scale

GLS-related options :

GLS_MIXING

Activate Generic Length Scale scheme, default is k-epsilon (see below)

GLS_KOMEGA

Activate K-OMEGA (OMEGA=frequency of TKE dissipation) originating from Kolmogorov (1942)

GLS_KEPSILON

Activate K-EPSILON (EPSILON=TKE dissipation) as in Jones and Launder (1972)

GLS_GEN

Activate generic model of Umlauf and Burchard (2003)

CANUTO_A

Option for CANUTO A stability function (default, see below)

GibLau_78

Option for Gibson & Launder, 1978 stability function

MelYam_82

Option for Mellor & Yamada, 1982 stability function

KanCla_94

Option for Kantha & Clayson, 1994 stability function

Luyten_96

Option for Luyten, 1996 stability function

CANUTO_B

Option for CANUTO B stability function

Cheng_02

Option for Cheng, 2002 stability function

Preselected options for GLS:

#ifdef GLS_MIXING
# if   defined GLS_KOMEGA
# elif defined GLS_KEPSILON
# elif defined GLS_GEN
# else
#  define GLS_KEPSILON
# endif
# if   defined CANUTO_A
# elif defined GibLau_78
# elif defined MelYam_82
# elif defined KanCla_94
# elif defined Luyten_96
# elif defined CANUTO_B
# elif defined Cheng_02
# else
#  define CANUTO_A
# endif
#endif

The objective of this section is to describe the current implementation of a Generic Length Scale (GLS) turbulence scheme in CROCO that computes the turbulent viscosity \(K_m\) and diffusivity \(K_s\). First of all, as usually done in most implementations, the assumption of an horizontally homogeneous flow is made. Following Umlauf & Burchard (2003), the equations satisfied by the two prognostic variables \(k\) (the kinetic energy) and \(\psi\) (the generic length scale) are

\begin{eqnarray} \partial_t k &=& \partial_z ( K_k \partial_z k ) + P + B - \varepsilon, \qquad \qquad \qquad \qquad \qquad K_k = K_m / {\rm Sc}_{k} \nonumber \\ \partial_t \psi &=& \partial_z ( K_\psi \partial_z \psi ) + \psi k^{-1} \left( \beta_1 P + \beta_3^{\pm} B - \beta_2 \varepsilon \right), \qquad \qquad K_\psi = K_m / {\rm Sc}_{\psi} \nonumber \end{eqnarray}

where the \(\beta_j\) (j=1,3) are constants to be defined, \(P\) represents the TKE production by vertical shear \(P = K_m \left[ (\partial_z u)^2 + (\partial_z v)^2 \right]\) and \(B\) the TKE destruction by stratification \(B = - K_s N^2\) (with \(N^2\) the local Brunt-Vaisala frequency). The dissipation rate \(\varepsilon\) is related to the generic length scale \(\psi\) following

\[\varepsilon = (c_\mu^0)^{3+p/n} k^{3/2+m/n} \psi^{-1/n}, \qquad \psi = (c_\mu^0)^p k^m l^n, \qquad l = (c_\mu^0)^{3} k^{3/2} \varepsilon^{-1}\]

with \(l\) a mixing length and \(c_\mu^0\) a constant (whose value is between \(0.526\) and \(0.555\)) to be defined. Depending on the parameter values for the triplet \((m,n,p)\) the GLS scheme will either correspond to a \(k-\varepsilon\), a \(k-\omega\) or the so-called generic (Umlauf & Burchard, 2003) turbulence scheme (to simplify the code and because this scheme do not generally outperform other schemes, the possibility to use the so-called \(k\)-\(k\)\(l\) scheme is not implemented in Croco). Since the equations for \(e\) and \(\psi\) bear lots of similarities, to avoid excessive code duplication, a unique equation is solved for a quantity \({\cal T}_i\) encompassing \(k\) (when \(i = i_{\rm tke}\)) and \(\psi\) (when \(i = i_{\rm gls}\), \(i_{\rm gls} = i_{\rm tke} + 1\)) such that

\[\partial_t {\cal T}_i = \partial_z ( K_{{\cal T}_i} \partial_z {\cal T}_i ) + ( c_i^1 P + c_i^{3,\pm} B - c_i^2 \varepsilon), \qquad K_{{\cal T}_i} = K_m / {\rm Sc}_{{\cal T}_i}\]

where

\[{\rm Sc}_{{\cal T}_{i_{\rm tke}}} = {\rm Sc}_{k}, \qquad {\rm Sc}_{{\cal T}_{i_{\rm gls}}} = {\rm Sc}_{\psi}\]

and

\begin{eqnarray} c_i^1 &=& (i_{\rm gls} - i) + (i - i_{\rm tke}) \beta_1 e^{-1} \psi \nonumber \\ c_i^2 &=& (i_{\rm gls} - i) + (i - i_{\rm tke}) \beta_2 e^{-1} \psi \nonumber \\ c_i^{3,\pm} &=& (i_{\rm gls} - i) + (i - i_{\rm tke}) \beta_3^{\pm} e^{-1} \psi \nonumber \end{eqnarray}

In practice this explains why in the code the two prognostic quantities \(k\) and \(\psi\) are stored in a single array \({\rm trb(i,j,k,ntime,ngls)}\) avec \({\rm ngls} = 2\), \(i_{\rm tke} = 1\) and \(i_{\rm gls} = 2\). Once the quantities \(k\) and \(\psi\) (hence \(\varepsilon\)) are known, the turbulent viscosity/diffusivity are given by

\[K_m = c_\mu \left( \frac{k^2}{\varepsilon} \right) = \frac{c_\mu}{(c_\mu^0)^{3}} (l \sqrt{k}), \qquad K_s = c_\mu^{'} \left(\frac{k^2}{\varepsilon}\right) = \frac{c_\mu^{'}}{(c_\mu^0)^{3}} (l \sqrt{k}).\]

where \(c_\mu\) and \(c_\mu^{'}\) are determined through so-called stability functions (see below).

Choice of parameter values and stability functions

A particular GLS occurence is defined by the following parameters :

  • The exponents \((m,n,p)\) in the definition of \(\varepsilon\)

  • The Schmidt numbers \({\rm Sc}_{k}\) and \({\rm Sc}_{\psi}\)

  • The coefficients \(\beta_j\) (j=1,3)

  • The constant \(c_\mu^0\)

  • The stability functions which are generally function of

\[\alpha_M = \left(\frac{k}{\varepsilon}\right)^2 \left[ (\partial_z u)^2 + (\partial_z v)^2 \right], \qquad \alpha_N = \left(\frac{k}{\varepsilon}\right)^2 N^2\]

Where \((m,n,p)\), \({\rm Sc}_{k}\), \({\rm Sc}_{\psi}\), \(\beta_j\) (j=1,3) are tied to a particular choice of GLS scheme (see table below) while \(c_\mu^0\), \(c_\mu\) and \(c_\mu^{'}\) are tied to a particular choice of stability function. The formulation of numerous stability functions can be reconciled when written using the generic form

\begin{eqnarray} c_\mu &=& \frac{n_0 + n_1 \alpha_N + n_2 \alpha_M}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2} \nonumber \\ c_\mu^{'} &=& \frac{n_0^{'} + n_1^{'} \alpha_N + n_2^{'} \alpha_M}{d_0 + d_1 \alpha_N + d_2 \alpha_M + d_3 \alpha_N \alpha_M + d_4 \alpha_N^2 + d_5 \alpha_M^2} \nonumber \end{eqnarray}

where a given choice of stability function will define the parameter values for \(n_i\), \(d_j\), and \(n_{k}^{'}\). In Croco, 7 options are available, these are referred to CANUTO-A, CANUTO-B, Gibson & Launder (1978), Mellor & Yamada (1982), Kantha & Clayson (1994), Luyten (1996), Cheng (2002).

Table: parameter values corresponding to each particular GLS model

GLS model

\(m\)

\(n\)

\(p\)

\(\beta_1\)

\(\beta_2\)

\(\beta_3^{-}\)

\(\beta_3^{+}\)

\({\rm Sc}_{e}\)

\({\rm Sc}_{\psi}\)

\(k-\omega\)

0.5

-1

-1

0.555

0.833

-0.6

1

0.5

0.5

\(k-\varepsilon\)

1.5

-1

3

1.44

1.92

-0.4

1

1

0.7692

Gen

1

-0.67

0

1

1.22

0.05

1

1.25

0.9345

The quantities \(\alpha_N\) and \(\alpha_M\) in the formulation of \(c_\mu\) and \(c_\mu^{'}\) must satisfy some constraints to guarantee the regularity of numerical solutions. In CROCO, the following steps are done:

1. Apply the Galperin (1988) limitation i.e. \(l \le l_{\rm lim} = \beta_{\rm galp} \sqrt{2 k / N^2}\) on \(\psi\) with \(\beta_{\rm galp} = 0.53\). The first step is to use this mixing length \(l_{\rm lim}\) to compute \(\psi_{\rm min} = (c_\mu^0)^p k^m (l_{\rm lim})^n\) and to correct \(\psi\) to satisfy the constraint

\[\psi = \max \left( \psi, \psi_{\rm min} \right)\]

here the max function is used since the exponent \(n\) is negative whatever the GLS scheme.

2. Compute the dissipation rate \(\varepsilon = (c_\mu^0)^{3+p/n} k^{3/2+m/n} \psi^{-1/n}\) and correct it

\[\varepsilon = \max \left( \varepsilon, \varepsilon_{\rm min} \right), \qquad \varepsilon_{\rm min} = 10^{-12} \;{\rm m}^2\;{\rm s}^{-3}\]

3. Compute \(\alpha_N\) and \(\alpha_M\), and apply “stability and realisability” constraints following Umlauf & Burchard (2003) (their Sec. 4). A first constraint applies on \(\alpha_N\) to ensure that \(- \partial_{\alpha_N}(c_\mu^{'}/\alpha_N) > 0\) to prevent the occurence of oscillations in \(c_\mu^{'}\). This translates into the following limiter

\[\alpha_N^{\min} = \frac{ -( d_1 + n_0^{'} ) + \sqrt{ ( d_1 + n_0^{'} )^2 - 4 d_0 ( d_4 + n_1^{'} ) } }{ 2 ( d_4 + n_1^{'} ) }, \qquad \alpha_N = \min \left( \max( 0.73 \alpha_N^{\min} ) , 10^{10} \right)\]

where the coefficient 0.73 is used to ensure the so-called realisability and has been empirically computed thanks to Table 3 in Umlauf & Burchard (2003) in order to satisfy their constraint (48). Then an upper limit is applied on \(\alpha_M\) to ensure that \(\partial_{\alpha_M}(c_\mu \sqrt{\alpha_M}) \ge 0\) which is also a prerequisite for stability reasons

\[\alpha_M^{\max} = \frac{d_0 n_0 + ( d_0 n_1 + d_1 n_0 ) \alpha_N + (d_1 n_1 + d_4 n_0) \alpha_N^2 + d_4 n_1 \alpha_N^3}{ d_2 n_0 + (d_2 n_1 + d_3 n_0) \alpha_N + (d_3 n_1) \alpha_N^2}, \qquad \alpha_M = \min \left( \alpha_M, \alpha_M^{\max} \right)\]

Once those quantities are computed, stability functions are evaluated as well as the turbulent viscosity/diffusivity.

Surface and bottom boundary conditions

In current version of Croco, both \(k\) and \(\psi\) are formulated with Neumann boundary conditions at the top and at the bottom. However the nature of those boundary conditions also requires the determination of bottom and surface values for \(k\) and \(\psi\).

  • For turbulent kinetic energy, the “diagnostic” surface and bottom values are given by

\[k_{\rm sfc} = ( u_\star^{s} / c_{\mu}^0 )^2, \qquad k_{\rm bot} = ( u_\star^{b} / c_{\mu}^0 )^2\]

and simple homogeneous Neumann boundary conditions are applied

\[\left. K_k \partial_z k \right|_{\rm sfc} = 0, \qquad \left. K_k \partial_z k \right|_{\rm bot} = 0\]

In practice, due to the placement of \(k\) and \(\psi\) on the computational grid, the Neumann boundary condition is not applied strictly at the surface (resp. at the bottom) but at \(z = z_N\) (resp. \(z = z_1\)) whereas the surface (resp. bottom) is located at \(z=z_{N+1/2}\) (resp. \(z=z_{1/2}\)) with \(N\) the number of vertical levels (i.e. the number of cells in the vertical).

  • For the generic length scale, a roughness is defined as

\[z_{0,s} = \max\left\{ 10^{-2}\;{\rm m}, \frac{C_{\rm ch}}{g} (u_\star^{s})^2 \right\}, \qquad C_{\rm ch} = 1400\]

at the surface and

\[z_{0,b} = \max\left\{ 10^{-4}\;{\rm m}, {\rm Zob} \right\}\]

at the bottom with \({\rm Zob}\) a user defined roughness length (usually \({\rm Zob} = 10^{-2}\;{\rm m}\)). Again, the boundary conditions are applied at the center of the shallowest and deepest grid cells and not at their interfaces which means that the relevant length scales are

\[L_{\rm sfc} = \kappa \left( \frac{\varDelta z_N}{2} + z_{0,s} \right), \qquad L_{\rm bot} = \kappa \left( \frac{\varDelta z_1}{2} + z_{0,b}\right)\]

with \(\kappa\) the von Karman constant. Moreover TKE values are interpolated at \(z = z_N\) and \(z = z_1\)

\[\widetilde{k}_{\rm sfc} = \frac{1}{2}\left( k_{\rm sfc} + k_{\rm N-1/2} \right), \qquad \widetilde{k}_{\rm bot} = \frac{1}{2}\left( k_{\rm bot} + k_{\rm 3/2} \right)\]

where \(k_{\rm sfc}\) and \(k_{\rm bot}\) are the diagnostic values given above. The “diagnostic” surface and bottom values for \(\psi\) are thus given by

\[\psi_{\rm sfc} = (c_\mu^0)^p (L_{\rm sfc})^n (\widetilde{k}_{\rm sfc})^m, \qquad \psi_{\rm bot} = (c_\mu^0)^p (L_{\rm bot})^n (\widetilde{k}_{\rm bot})^m\]

Then the surface and bottom flux are defined as

\[{\cal F}_{\psi}^{\rm sfc} = \left. K_{\psi} \partial_z \psi \right|_{\rm sfc} = - n (c_\mu^0)^{p+1} \frac{\kappa}{{\rm Sc}_{\psi}} (\widetilde{k}_{\rm sfc})^{m+1/2} (L_{\rm sfc})^n\]
\[{\cal F}_{\psi}^{\rm bot} = \left. K_{\psi} \partial_z \psi \right|_{\rm bot} = - n (c_\mu^0)^{p+1} \frac{\kappa}{{\rm Sc}_{\psi}} (\widetilde{k}_{\rm bot})^{m+1/2} (L_{\rm bot})^n\]

which correspond to the Neumann boundary conditions applied in the code.