Vertical Mixing

The discrete form of the ocean subgrid scale physics has been presented in §5.3 and §6.7. At the surface and bottom boundaries, the turbulent fluxes of momentum, heat and salt have to be defined. At the surface they are prescribed from the surface forcing (see Chap. 7), while at the bottom they are set to zero for heat and salt, unless a geothermal flux forcing is prescribed as a bottom boundary condition ($i.e.$ key_ trabbl defined, see §5.4.3), and specified through a bottom friction parameterisation for momentum (see §10.4).

In this section we briefly discuss the various choices offered to compute the vertical eddy viscosity and diffusivity coefficients, $A_u^{vm}$ , $A_v^{vm}$ and $A^{vT}$ ($A^{vS}$), defined at $uw$-, $vw$- and $w$- points, respectively (see §5.3 and §6.7). These coefficients can be assumed to be either constant, or a function of the local Richardson number, or computed from a turbulent closure model (TKE, GLS or KPP formulation). The computation of these coefficients is initialized in the zdfini.F90 module and performed in the zdfric.F90, zdftke.F90, zdfgls.F90 or zdfkpp.F90 modules. The trends due to the vertical momentum and tracer diffusion, including the surface forcing, are computed and added to the general trend in the dynzdf.F90 and trazdf.F90 modules, respectively. These trends can be computed using either a forward time stepping scheme (namelist parameter ln_zdfexp=true) or a backward time stepping scheme (ln_zdfexp=false) depending on the magnitude of the mixing coefficients, and thus of the formulation used (see §3).

Constant (key_ zdfcst)

!-----------------------------------------------------------------------
&namzdf        !   vertical physics
!-----------------------------------------------------------------------
   rn_avm0     =   1.2e-4  !  vertical eddy viscosity   [m2/s]          (background Kz if not "key_zdfcst")
   rn_avt0     =   1.2e-5  !  vertical eddy diffusivity [m2/s]          (background Kz if not "key_zdfcst")
   nn_avb      =    0      !  profile for background avt & avm (=1) or not (=0)
   nn_havtb    =    0      !  horizontal shape for avtb (=1) or not (=0)
   ln_zdfevd   = .true.    !  enhanced vertical diffusion (evd) (T) or not (F)
   nn_evdm     =    0      !  evd apply on tracer (=0) or on tracer and momentum (=1)
   rn_avevd    =  100.     !  evd mixing coefficient [m2/s]
   ln_zdfnpc   = .false.   !  Non-Penetrative Convective algorithm (T) or not (F)
   nn_npc      =    1            !  frequency of application of npc
   nn_npcp     =  365            !  npc control print frequency
   ln_zdfexp   = .false.   !  time-stepping: split-explicit (T) or implicit (F) time stepping
   nn_zdfexp   =    3            !  number of sub-timestep for ln_zdfexp=T
/

Options are defined through the namzdf namelist variables. When key_ zdfcst is defined, the momentum and tracer vertical eddy coefficients are set to constant values over the whole ocean. This is the crudest way to define the vertical ocean physics. It is recommended that this option is only used in process studies, not in basin scale simulations. Typical values used in this case are:

$$A_u^{vm} = A_v^{vm} = 1.2 10^{-4} m^2.s^{-1}$$

$$A^{vT} = A^{vS} = 1.2 10^{-5} m^2.s^{-1}$$

These values are set through the rn_avm0 and rn_avt0 namelist parameters. In all cases, do not use values smaller that those associated with the molecular viscosity and diffusivity, that is $\sim10^{-6} m^2.s^{-1}$ for momentum, $\sim10^{-7} m^2.s^{-1}$ for temperature and $\sim10^{-9} m^2.s^{-1}$ for salinity.

Richardson Number Dependent (key_ zdfric)

!-----------------------------------------------------------------------
&namzdf_ric    !   richardson number dependent vertical diffusion       ("key_zdfric" )
!-----------------------------------------------------------------------
   rn_avmri    = 100.e-4   !  maximum value of the vertical viscosity
   rn_alp      =   5.      !  coefficient of the parameterization
   nn_ric      =   2       !  coefficient of the parameterization
   rn_ekmfc    =   0.7     !  Factor in the Ekman depth Equation
   rn_mldmin   =   1.0     !  minimum allowable mixed-layer depth estimate (m)
   rn_mldmax   =1000.0     !  maximum allowable mixed-layer depth estimate (m)
   rn_wtmix    =  10.0     !  vertical eddy viscosity coeff [m2/s] in the mixed-layer
   rn_wvmix    =  10.0     !  vertical eddy diffusion coeff [m2/s] in the mixed-layer
   ln_mldw     = .true.    !  Flag to use or not the mized layer depth param.
/

When key_ zdfric is defined, a local Richardson number dependent formulation for the vertical momentum and tracer eddy coefficients is set through the namzdf_ric namelist variables.The vertical mixing coefficients are diagnosed from the large scale variables computed by the model. In situ measurements have been used to link vertical turbulent activity to large scale ocean structures. The hypothesis of a mixing mainly maintained by the growth of Kelvin-Helmholtz like instabilities leads to a dependency between the vertical eddy coefficients and the local Richardson number ($i.e.$ the ratio of stratification to vertical shear). Following Pacanowski and Philander [1981], the following formulation has been implemented:

$$\left\{ \begin{aligned}A^{vT} &= \frac {A_{ric}^{vT}}{\left( 1+... ...vT} }{\left( 1+ a \;Ri \right) } + A_b^{vm} \end{aligned} \right.$$

where $Ri = N^2 / \left(\partial_z \textbf{U}_h \right)^2$ is the local Richardson number, $N$ is the local Brunt-Vaisälä frequency (see §5.8.2), $A_b^{vT}$ and $A_b^{vm}$ are the constant background values set as in the constant case (see §10.1.1), and $A_{ric}^{vT} = 10^{-4} m^2.s^{-1}$ is the maximum value that can be reached by the coefficient when $Ri\leq 0$, $a=5$ and $n=2$. The last three values can be modified by setting the rn_avmri, rn_alp and nn_ric namelist parameters, respectively.

A simple mixing-layer model to transfer and dissipate the atmospheric forcings (wind-stress and buoyancy fluxes) can be activated setting the ln_mldw =.true. in the namelist.

In this case, the local depth of turbulent wind-mixing or "Ekman depth" $h_{e}(x,y,t)$ is evaluated and the vertical eddy coefficients prescribed within this layer.

This depth is assumed proportional to the "depth of frictional influence" that is limited by rotation:

$$h_{e} = Ek \frac {u^{*}} {f_{0}}$$ (10.2)

where, $Ek$ is an empirical parameter, $u^{*}$ is the friction velocity and $f_{0}$ is the Coriolis parameter.

In this similarity height relationship, the turbulent friction velocity:

$$u^{*} = \sqrt \frac {\vert\tau\vert} {\rho_o}$$ (10.3)

is computed from the wind stress vector $\vert\tau\vert$ and the reference density $\rho_o$. The final $h_{e}$ is further constrained by the adjustable bounds rn_mldmin and rn_mldmax. Once $h_{e}$ is computed, the vertical eddy coefficients within $h_{e}$ are set to the empirical values rn_wtmix and rn_wvmix [Lermusiaux, 2001].

TKE Turbulent Closure Scheme (key_ zdftke)

!-----------------------------------------------------------------------
&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  ("key_zdftke")
!-----------------------------------------------------------------------
   rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
   rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
   rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
   rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
   rn_bshear   =   1.e-20  ! background shear (>0) currently a numerical threshold (do not change it)
   nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                           !                 = 1 bounded by the local vertical scale factor
                           !                 = 2 first vertical derivative of mixing length bounded by 1
                           !                 = 3 as =2 with distinct disspipative an mixing length scale
   nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
   rn_lc       =   0.15    !  coef. associated to Langmuir cells
   nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
                           !        = 0 no penetration
                           !        = 1 add a tke source below the ML
                           !        = 2 add a tke source just at the base of the ML
                           !        = 3 as = 1 applied on HF part of the stress    ("key_oasis3")
   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                           !        = 0  constant 10 m length scale
                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
/

The vertical eddy viscosity and diffusivity coefficients are computed from a TKE turbulent closure model based on a prognostic equation for $\bar {e}$, the turbulent kinetic energy, and a closure assumption for the turbulent length scales. This turbulent closure model has been developed by Bougeault and Lacarrere [1989] in the atmospheric case, adapted by Gaspar et al. [1990] for the oceanic case, and embedded in OPA, the ancestor of NEMO, by Blanke and Delecluse [1993] for equatorial Atlantic simulations. Since then, significant modifications have been introduced by Madec et al. [1998] in both the implementation and the formulation of the mixing length scale. The time evolution of $\bar {e}$ is the result of the production of $\bar {e}$ through vertical shear, its destruction through stratification, its vertical diffusion, and its dissipation of Kolmogorov [1942] type:

$$\frac{\partial \bar{e}}{\partial t} = \frac{K_m}{{e_3}^2 }\;\left... ...\bar{e}}{\partial k}} \right] - c_\epsilon \;\frac{\bar {e}^{3/2}}{l_\epsilon }$$ (10.4)

$$\begin{split}K_m &= C_k l_k \sqrt {\bar{e}\; } K_\rho &= A^{vm} / P_{rt} \end{split}$$ (10.5)

where $N$ is the local Brunt-Vaisälä frequency (see §5.8.2), $l_{\epsilon }$ and $l_{\kappa }$ are the dissipation and mixing length scales, $P_{rt}$ is the Prandtl number, $K_m$ and $K_\rho$ are the vertical eddy viscosity and diffusivity coefficients. The constants $C_k = 0.1$ and $C_\epsilon = \sqrt {2} /2$ $\approx 0.7$ are designed to deal with vertical mixing at any depth [Gaspar et al., 1990]. They are set through namelist parameters nn_ediff and nn_ediss. $P_{rt}$ can be set to unity or, following Blanke and Delecluse [1993], be a function of the local Richardson number, $R_i$:

$$P_{rt} = \begin{cases} 1 & \text{if $ R_i \leq 0.2$} 5 ... ...xt{if $ 0.2 \leq R_i \leq 2$} 10 & \text{if $ 2 \leq R_i$} \end{cases}$$

Options are defined through the namzdfy_tke namelist variables. The choice of $P_{rt}$ is controlled by the nn_pdl namelist variable.

At the sea surface, the value of $\bar {e}$ is prescribed from the wind stress field as $\bar{e}_o = e_{bb} \vert\tau\vert / \rho_o$, with $e_{bb}$ the rn_ebb namelist parameter. The default value of $e_{bb}$ is 3.75. [Gaspar et al., 1990]), however a much larger value can be used when taking into account the surface wave breaking (see below Eq. (10.10)). The bottom value of TKE is assumed to be equal to the value of the level just above. The time integration of the $\bar {e}$ equation may formally lead to negative values because the numerical scheme does not ensure its positivity. To overcome this problem, a cut-off in the minimum value of $\bar {e}$ is used (rn_emin namelist parameter). Following Gaspar et al. [1990], the cut-off value is set to $\sqrt{2}/2 10^{-6} m^2.s^{-2}$. This allows the subsequent formulations to match that of Gargett [1984] for the diffusion in the thermocline and deep ocean : $K_\rho = 10^{-3} / N$. In addition, a cut-off is applied on $K_m$ and $K_\rho$ to avoid numerical instabilities associated with too weak vertical diffusion. They must be specified at least larger than the molecular values, and are set through rn_avm0 and rn_avt0 (namzdf namelist, see §10.1.1).

Turbulent length scale

For computational efficiency, the original formulation of the turbulent length scales proposed by Gaspar et al. [1990] has been simplified. Four formulations are proposed, the choice of which is controlled by the nn_mxl namelist parameter. The first two are based on the following first order approximation [Blanke and Delecluse, 1993]:

$$l_k = l_\epsilon = \sqrt {2 \bar{e}\; } / N$$ (10.6)

which is valid in a stable stratified region with constant values of the Brunt- Vaisälä frequency. The resulting length scale is bounded by the distance to the surface or to the bottom (nn_mxl = 0) or by the local vertical scale factor (nn_mxl = 1). Blanke and Delecluse [1993] notice that this simplification has two major drawbacks: it makes no sense for locally unstable stratification and the computation no longer uses all the information contained in the vertical density profile. To overcome these drawbacks, Madec et al. [1998] introduces the nn_mxl = 2 or 3 cases, which add an extra assumption concerning the vertical gradient of the computed length scale. So, the length scales are first evaluated as in (10.6) and then bounded such that:

$$\frac{1}{e_3 }\left\vert {\frac{\partial l}{\partial k}} \right\vert \leq 1 l = l_k = l_\epsilon$$ (10.7)

(10.7) means that the vertical variations of the length scale cannot be larger than the variations of depth. It provides a better approximation of the Gaspar et al. [1990] formulation while being much less time consuming. In particular, it allows the length scale to be limited not only by the distance to the surface or to the ocean bottom but also by the distance to a strongly stratified portion of the water column such as the thermocline (Fig. 10.1). In order to impose the (10.7) constraint, we introduce two additional length scales: $l_{up}$ and $l_{dwn}$, the upward and downward length scales, and evaluate the dissipation and mixing length scales as (and note that here we use numerical indexing):

Figure 10.1: Illustration of the mixing length computation.

$$\begin{aligned}l_{up }^{(k)} &= \min \left( l^{(k)} , l_... ...} \right) \quad &\text{ from $k=jpk$ to $1$ } \end{aligned}$$

where $l^{(k)}$ is computed using (10.6), $i.e.$ $l^{(k)} = \sqrt {2 {\bar e}^{(k)} / {N^2}^{(k)} }$.

In the nn_mxl = 2 case, the dissipation and mixing length scales take the same value: $l_k= l_\epsilon = \min \left( l_{up} \;,\; l_{dwn} \right)$, while in the nn_mxl = 3 case, the dissipation and mixing turbulent length scales are give as in Gaspar et al. [1990]:

$$\begin{aligned}& l_k = \sqrt{ l_{up} l_{dwn} } & l_\epsilon = \min \left( l_{up} \;,\; l_{dwn} \right) \end{aligned}$$

At the ocean surface, a non zero length scale is set through the rn_mxl0 namelist parameter. Usually the surface scale is given by $l_o = \kappa z_o$ where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness parameter of the surface. Assuming $z_o=0.1$ m [Craig and Banner, 1994] leads to a 0.04 m, the default value of rn_mxl0. In the ocean interior a minimum length scale is set to recover the molecular viscosity when $\bar {e}$ reach its minimum value ( $1.10^{-6}= C_k l_{min} \sqrt{\bar{e}_{min}}$ ).

Surface wave breaking parameterization

Following Mellor and Blumberg [2004], the TKE turbulence closure model has been modified to include the effect of surface wave breaking energetics. This results in a reduction of summertime surface temperature when the mixed layer is relatively shallow. The Mellor and Blumberg [2004] modifications acts on surface length scale and TKE values and air-sea drag coefficient. The latter concerns the bulk formulea and is not discussed here.

Following Craig and Banner [1994], the boundary condition on surface TKE value is :

$$\bar{e}_o = \frac{1}{2} \left( 15.8 \alpha_{CB} \right)^{2/3} \frac{\vert\tau\vert}{\rho_o}$$ (10.10)

where $\alpha_{CB}$ is the Craig and Banner [1994] constant of proportionality which depends on the ”wave age”, ranging from 57 for mature waves to 146 for younger waves [Mellor and Blumberg, 2004]. The boundary condition on the turbulent length scale follows the Charnock's relation:

$$l_o = \kappa \beta \frac{\vert\tau\vert}{g \rho_o}$$ (10.11)

where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant. Mellor and Blumberg [2004] suggest $\beta = 2.10^{5}$ the value chosen by Stacey [1999] citing observation evidence, and $\alpha_{CB} = 100$ the Craig and Banner's value. As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} \vert\tau\vert / \rho_o$, with $e_{bb}$ the rn_ebb namelist parameter, setting rn_ebb = 67.83 corresponds to $\alpha_{CB} = 100$. Further setting ln_mxl0 to true applies (10.11) as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value. Note that a minimal threshold of rn_emin0 $=10^{-4} m^2.s^{-2}$ (namelist parameters) is applied on surface $\bar {e}$ value.

Langmuir cells

Langmuir circulations (LC) can be described as ordered large-scale vertical motions in the surface layer of the oceans. Although LC have nothing to do with convection, the circulation pattern is rather similar to so-called convective rolls in the atmospheric boundary layer. The detailed physics behind LC is described in, for example, Craik and Leibovich [1976]. The prevailing explanation is that LC arise from a nonlinear interaction between the Stokes drift and wind drift currents.

Here we introduced in the TKE turbulent closure the simple parameterization of Langmuir circulations proposed by [Axell, 2002] for a $k-\epsilon$ turbulent closure. The parameterization, tuned against large-eddy simulation, includes the whole effect of LC in an extra source terms of TKE, $P_{LC}$. The presence of $P_{LC}$ in (10.4), the TKE equation, is controlled by setting ln_lc to true in the namtke namelist.

By making an analogy with the characteristic convective velocity scale ($e.g.$, D'Alessio et al. [1998]), $P_{LC}$ is assumed to be :

$$P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}}$$ (10.12)

where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth. With no information about the wave field, $w_{LC}$ is assumed to be proportional to the Stokes drift $u_s = 0.377 \vert\tau\vert^{1/2}$, where $\vert\tau\vert$ is the surface wind stress module ^10.1. For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at a finite depth $H_{LC}$ (which is often close to the mixed layer depth), and simply varies as a sine function in between (a first-order profile for the Langmuir cell structures). The resulting expression for $w_{LC}$ is :

$$w_{LC} = \begin{cases}c_{LC} u_s \sin(- \pi z / H_{LC} ) & \text{if $-z \leq H_{LC}$} 0 & \text{otherwise} \end{cases}$$ (10.13)

where $c_{LC} = 0.15$ has been chosen by [Axell, 2002] as a good compromise to fit LES data. The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second. The value of $c_{LC}$ is set through the rn_lc namelist parameter, having in mind that it should stay between 0.15 and 0.54 [Axell, 2002].

The $H_{LC}$ is estimated in a similar way as the turbulent length scale of TKE equations: $H_{LC}$ is depth to which a water parcel with kinetic energy due to Stoke drift can reach on its own by converting its kinetic energy to potential energy, according to

$$- \int_{-H_{LC}}^0 { N^2\;z \;dz} = \frac{1}{2} u_s^2$$ (10.14)

Mixing just below the mixed layer

Vertical mixing parameterizations commonly used in ocean general circulation models tend to produce mixed-layer depths that are too shallow during summer months and windy conditions. This bias is particularly acute over the Southern Ocean. To overcome this systematic bias, an ad hoc parameterization is introduced into the TKE scheme Rodgers et al. [2014]. The parameterization is an empirical one, $i.e.$ not derived from theoretical considerations, but rather is meant to account for observed processes that affect the density structure of the ocean’s planetary boundary layer that are not explicitly captured by default in the TKE scheme ($i.e.$ near-inertial oscillations and ocean swells and waves).

When using this parameterization ($i.e.$ when nn_etau = 1), the TKE input to the ocean ($S$) imposed by the winds in the form of near-inertial oscillations, swell and waves is parameterized by (10.10) the standard TKE surface boundary condition, plus a depth depend one given by:

$$S = (1-f_i) \; f_r \; e_s \; e^{-z / h_\tau}$$ (10.15)

where $z$ is the depth, $e_s$ is TKE surface boundary condition, $f_r$ is the fraction of the surface TKE that penetrate in the ocean, $h_\tau$ is a vertical mixing length scale that controls exponential shape of the penetration, and $f_i$ is the ice concentration (no penetration if $f_i=1$, that is if the ocean is entirely covered by sea-ice). The value of $f_r$, usually a few percents, is specified through rn_efr namelist parameter. The vertical mixing length scale, $h_\tau$, can be set as a 10 m uniform value (nn_etau = 0) or a latitude dependent value (varying from 0.5 m at the Equator to a maximum value of 30 m at high latitudes (nn_etau = 1).

Note that two other option existe, nn_etau = 2, or 3. They correspond to applying (10.15) only at the base of the mixed layer, or to using the high frequency part of the stress to evaluate the fraction of TKE that penetrate the ocean. Those two options are obsolescent features introduced for test purposes. They will be removed in the next release.

TKE discretization considerations (key_ zdftke)

Figure 10.2: Illustration of the TKE time integration and its links to the momentum and tracer time integration.

The production of turbulence by vertical shear (the first term of the right hand side of (10.4)) should balance the loss of kinetic energy associated with the vertical momentum diffusion (first line in (2.34)). To do so a special care have to be taken for both the time and space discretization of the TKE equation [Marsaleix et al., 2008, Burchard, 2002].

Let us first address the time stepping issue. Fig. 10.2 shows how the two-level Leap-Frog time stepping of the momentum and tracer equations interplays with the one-level forward time stepping of TKE equation. With this framework, the total loss of kinetic energy (in 1D for the demonstration) due to the vertical momentum diffusion is obtained by multiplying this quantity by $u^t$ and summing the result vertically:

$$\begin{split}\int_{-H}^{\eta} u^t \partial_z &\left( {K_m}^... ...^t \partial_z{u^t} \partial_z u^{t+\rdt} dz } \end{split}$$ (10.16)

Here, the vertical diffusion of momentum is discretized backward in time with a coefficient, $K_m$, known at time $t$ (Fig. 10.2), as it is required when using the TKE scheme (see §3.3). The first term of the right hand side of (10.16) represents the kinetic energy transfer at the surface (atmospheric forcing) and at the bottom (friction effect). The second term is always negative. It is the dissipation rate of kinetic energy, and thus minus the shear production rate of $\bar {e}$. (10.16) implies that, to be energetically consistent, the production rate of $\bar {e}$ used to compute $(\bar{e})^t$ (and thus ${K_m}^t$) should be expressed as ${K_m}^{t-\rdt} (\partial_z u)^{t-\rdt} (\partial_z u)^t$ (and not by the more straightforward $K_m \left( \partial_z u \right)^2$ expression taken at time $t$ or $t-\rdt$).

A similar consideration applies on the destruction rate of $\bar {e}$ due to stratification (second term of the right hand side of (10.4)). This term must balance the input of potential energy resulting from vertical mixing. The rate of change of potential energy (in 1D for the demonstration) due vertical mixing is obtained by multiplying vertical density diffusion tendency by $g z$ and and summing the result vertically:

$$\begin{split}\int_{-H}^{\eta} g z \partial_z &\left( {K_\rh... ...\rho^{t+\rdt} {K_\rho}^t (N^2)^{t+\rdt} dz } \end{split}$$ (10.17)

where we use $N^2 = -g \partial_k \rho / (e_3 \rho)$. The first term of the right hand side of (10.17) is always zero because there is no diffusive flux through the ocean surface and bottom). The second term is minus the destruction rate of $\bar {e}$ due to stratification. Therefore (10.16) implies that, to be energetically consistent, the product ${K_\rho}^{t-\rdt} (N^2)^t$ should be used in (10.4), the TKE equation.

Let us now address the space discretization issue. The vertical eddy coefficients are defined at $w$-point whereas the horizontal velocity components are in the centre of the side faces of a $t$-box in staggered C-grid (Fig.4.1). A space averaging is thus required to obtain the shear TKE production term. By redoing the (10.16) in the 3D case, it can be shown that the product of eddy coefficient by the shear at $t$ and $t-\rdt$ must be performed prior to the averaging. Furthermore, the possible time variation of $e_3$ (key_ vvl case) have to be taken into account.

The above energetic considerations leads to the following final discrete form for the TKE equation:

$$\begin{split}\frac { (\bar{e})^t - (\bar{e})^{t-\rdt} } {\rdt... ...}}{l_\epsilon}\right)^{t-\rdt} (\bar {e})^{t+\rdt} \end{split}$$ (10.18)

where the last two terms in (10.18) (vertical diffusion and Kolmogorov dissipation) are time stepped using a backward scheme (see§3.3). Note that the Kolmogorov term has been linearized in time in order to render the implicit computation possible. The restart of the TKE scheme requires the storage of $\bar {e}$, $K_m$, $K_\rho$ and $l_\epsilon$ as they all appear in the right hand side of (10.18). For the latter, it is in fact the ratio $\sqrt{\bar{e}}/l_\epsilon$ which is stored.

GLS Generic Length Scale (key_ zdfgls)

!-----------------------------------------------------------------------
&namzdf_gls                !   GLS vertical diffusion                   ("key_zdfgls")
!-----------------------------------------------------------------------
   rn_emin       = 1.e-7   !  minimum value of e   [m2/s2]
   rn_epsmin     = 1.e-12  !  minimum value of eps [m2/s3]
   ln_length_lim = .true.  !  limit on the dissipation rate under stable stratification (Galperin et al., 1988)
   rn_clim_galp  = 0.267   !  galperin limit
   ln_sigpsi     = .true.  !  Activate or not Burchard 2001 mods on psi schmidt number in the wb case
   rn_crban      = 100.    !  Craig and Banner 1994 constant for wb tke flux
   rn_charn      = 70000.  !  Charnock constant for wb induced roughness length
   rn_hsro       =  0.02   !  Minimum surface roughness
   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met=2)
   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2)
   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum)
   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum)
   nn_stab_func  =     2   !  stability function (0=Galp, 1= KC94, 2=CanutoA, 3=CanutoB)
   nn_clos       =     1   !  predefined closure type (0=MY82, 1=k-eps, 2=k-w, 3=Gen)
/

The Generic Length Scale (GLS) scheme is a turbulent closure scheme based on two prognostic equations: one for the turbulent kinetic energy $\bar {e}$, and another for the generic length scale, $\psi$ [Umlauf and Burchard, 2003, Umlauf and Burchard, 2005]. This later variable is defined as : $\psi = {C_{0\mu}}^{p} {\bar{e}}^{m} l^{n}$, where the triplet $(p, m, n)$ value given in Tab.10.1 allows to recover a number of well-known turbulent closures ($k$-$kl$ [Mellor and Yamada, 1982], $k$-$\epsilon$ [Rodi, 1987], $k$-$\omega$ [Wilcox, 1988] among others [Kantha and Carniel, 2005, Umlauf and Burchard, 2003]). The GLS scheme is given by the following set of equations:

$$\frac{\partial \bar{e}}{\partial t} = \frac{K_m}{\sigma_e e_3 }\;... ... \left[ \frac{K_m}{e_3} \frac{\partial \bar{e}}{\partial k} \right] - \epsilon$$ (10.19)

$$\begin{split}\frac{\partial \psi}{\partial t} =& \frac{\psi}{... ...e_3 } \;\frac{\partial \psi}{\partial k}} \right]\; \end{split}$$ (10.20)

$$\begin{split}K_m &= C_{\mu} \sqrt {\bar{e}} l K_\rho &= C_{\mu'} \sqrt {\bar{e}} l \end{split}$$ (10.21)

$${\epsilon} = C_{0\mu} \frac{\bar {e}^{3/2}}{l} \;$$ (10.22)

where $N$ is the local Brunt-Vaisälä frequency (see §5.8.2) and $\epsilon$ the dissipation rate. The constants $C_1$, $C_2$, $C_3$, ${\sigma_e}$, ${\sigma_{\psi}}$ and the wall function ($Fw$) depends of the choice of the turbulence model. Four different turbulent models are pre-defined (Tab.10.1). They are made available through the nn_clo namelist parameter.

Table 10.1: Set of predefined GLS parameters, or equivalently predefined turbulence models available with key_ zdfgls and controlled by the nn_clos namelist variable in namzdf_gls .

	$k-kl$	$k-\epsilon$	$k-\omega$	generic
nn_clo	0	1	2	3
$( p , n , m )$	( 0 , 1 , 1 )	( 3 , 1.5 , -1 )	( -1 , 0.5 , -1 )	( 2 , 1 , -0.67 )
$\sigma_k$	2.44	1.	2.	0.8
$\sigma_\psi$	2.44	1.3	2.	1.07
$C_1$	0.9	1.44	0.555	1.
$C_2$	0.5	1.92	0.833	1.22
$C_3$	1.	1.	1.	1.
$F_{wall}$	Yes	-	-	-

In the Mellor-Yamada model, the negativity of $n$ allows to use a wall function to force the convergence of the mixing length towards $K z_b$ ($K$: Kappa and $z_b$: rugosity length) value near physical boundaries (logarithmic boundary layer law). $C_{\mu}$ and $C_{\mu'}$ are calculated from stability function proposed by Galperin et al. [1988], or by Kantha and Clayson [1994] or one of the two functions suggested by Canuto et al. [2001] (nn_stab_func = 0, 1, 2 or 3, resp.). The value of $C_{0\mu}$ depends of the choice of the stability function.

The surface and bottom boundary condition on both $\bar {e}$ and $\psi$ can be calculated thanks to Dirichlet or Neumann condition through nn_tkebc_surf and nn_tkebc_bot, resp. As for TKE closure , the wave effect on the mixing is considered when ln_crban = true [Mellor and Blumberg, 2004, Craig and Banner, 1994]. The rn_crban namelist parameter is $\alpha_{CB}$ in (10.10) and rn_charn provides the value of $\beta$ in (10.11).

The $\psi$ equation is known to fail in stably stratified flows, and for this reason almost all authors apply a clipping of the length scale as an ad hoc remedy. With this clipping, the maximum permissible length scale is determined by $l_{max} = c_{lim} \sqrt{2\bar{e}}/ N$. A value of $c_{lim} = 0.53$ is often used [Galperin et al., 1988]. Umlauf and Burchard [2005] show that the value of the clipping factor is of crucial importance for the entrainment depth predicted in stably stratified situations, and that its value has to be chosen in accordance with the algebraic model for the turbulent fluxes. The clipping is only activated if ln_length_lim=true, and the $c_{lim}$ is set to the rn_clim_galp value.

The time and space discretization of the GLS equations follows the same energetic consideration as for the TKE case described in §10.1.4 [Burchard, 2002]. Examples of performance of the 4 turbulent closure scheme can be found in Warner et al. [2005].

K Profile Parametrisation (KPP) (key_ zdfkpp)

!------------------------------------------------------------------------
&namzdf_kpp    !   K-Profile Parameterization dependent vertical mixing  ("key_zdfkpp", and optionally:
!------------------------------------------------------------------------ "key_kppcustom" or "key_kpplktb")
   ln_kpprimix = .true.    !  shear instability mixing
   rn_difmiw   =  1.0e-04  !  constant internal wave viscosity [m2/s]
   rn_difsiw   =  0.1e-04  !  constant internal wave diffusivity [m2/s]
   rn_riinfty  =  0.8      !  local Richardson Number limit for shear instability
   rn_difri    =  0.0050   !  maximum shear mixing at Rig = 0    [m2/s]
   rn_bvsqcon  = -0.01e-07 !  Brunt-Vaisala squared for maximum convection [1/s2]
   rn_difcon   =  1.       !  maximum mixing in interior convection [m2/s]
   nn_avb      =  0        !  horizontal averaged (=1) or not (=0) on avt and amv
   nn_ave      =  1        !  constant (=0) or profile (=1) background on avt
/

The KKP scheme has been implemented by J. Chanut ... Options are defined through the namzdf_kpp namelist variables.

Note that KPP is an obsolescent feature of the NEMO system. It will be removed in the next release (v3.7 and followings).

Footnotes

... module^10.1

Following Li and Garrett [1993], the surface Stoke drift velocity may be expressed as $u_s = 0.016 \vert U_{10m}\vert$. Assuming an air density of $\rho_a=1.22 Kg/m^3$ and a drag coefficient of $1.5 10^{-3}$ give the expression used of $u_s$ as a function of the module of surface stress