View source for Vertical Mixing Parameterizations

<div class="title">Vertical Mixing Parameterizations</div>

__TOC__

ROMS contains a variety of methods for setting the vertical viscous and diffusive coefficients. The choices range from simply choosing fixed values to the KPP, generic lengthscale (GLS) and Mellor-Yamada turbulence closure schemes.  See [[Bibliography#LargeWG_1994a | 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).

==K-Profile Parameterization==

The vertical mixing parameterization introduced by  [[Bibliography#LargeWG_1994a | Large, McWilliams
and Doney (1994)]] is a versatile first order scheme which has been shown to perform well in open ocean settings. Its design facilitates experimentation with additional or modified representations of specific turbulent processes.

===Surface boundary layer===

The Large, McWilliams and Doney scheme (LMD)
matches separate parameterizations for vertical mixing of the surface boundary layer and the ocean interior. A formulation based on boundary layer similarity theory is applied in the water column above a calculated boundary layer depth <math>h_{sbl}</math>. This parameterization is then matched at the interior with schemes to account for local shear, internal wave and double diffusive mixing effects.

Viscosity and diffusivities at model levels above a calculated surface boundary layer depth (<math>h_{sbl}</math> ) are expressed as the product of the length scale <math>h_{sbl}</math>, a turbulent velocity scale <math>w_x</math> and a non-dimensional shape function.

{| class="eqno"
|<math display="block">\nu_x = h_{sbl} w_x(\sigma)G_x(\sigma)</math><!--\eqno{(1)}-->
|(1)
|}

where <math>\sigma</math> is a non-dimensional coordinate ranging from 0 to 1 indicating depth within the surface boundary layer. The <math>x</math> subscript stands for one of momentum, temperature and salinity.

====Surface Boundary layer depth====

The boundary layer depth <math>h_{sbl}</math> is calculated as the minimum of the Ekman depth, estimated as,
{| class="eqno"
|<math display="block">h_e=0.7u_*/f</math><!--\eqno{(2)}-->
|(2)
|}  

(where <math>u_*</math> is the friction velocity <math>u_*=\sqrt{\tau_x^2+\tau_y^2}/\rho</math> ), the Monin-Obukhov depth:

{| class="eqno"
|<math display="block">L=u_*^3/(\kappa B_f)</math><!--\eqno{(3)}-->
|(3)
|}

(where <math>\kappa = 0.4</math> is von Karman's contant and <math>B_f</math> is the surface buoyancy flux), and the shallowest depth at which a critical bulk Richardson number is reached. The critical bulk Richardson number (<math>Ri_c</math>) is typically in the range 0.25--0.5. The bulk Richardson number (<math>Ri_b</math>) is calculated as:

{| class="eqno"
|<math display="block">
   Ri_b(z)=\frac{(B_r-B(d))d}{|\vec{V}_r-\vec{V}(d)|^2+{V_t}^2(d)}
</math><!--\eqno{(4)}-->
|(4)
|}

where <math>d</math> is distance down from the surface, <math>B</math> is the buoyancy, <math>B_r</math> is the buoyancy at a near surface reference depth, <math>\vec{V}</math> is the mean horizontal velocity, <math>\vec{V}_r</math> the velocity at the near surface reference depth and <math>V_t</math> is an estimate of the turbulent velocity contribution to velocity shear.

The turbulent velocity shear term in this equation is given by LMD as,

{| class="eqno"
|<math display="block">
  V_{t}^{2}(d)=\frac{C_v(-\beta_T)^{1/2}}{Ri_c
  \kappa}(c_s\epsilon)^{-1/2}dNw_s
</math><!--\eqno{(5)}-->
|(5)
|}

where <math>C_v</math> is the ratio of interior <math>N</math> to <math>N</math> at the entrainment depth, <math>\beta_T</math> is ratio of entrainment flux to surface buoyancy flux, <math>c_s</math> and <math>\epsilon</math> are constants, and <math>w_s</math> is the turbulent velocity scale for scalars. LMD derive (5) based on the expected behavior in the pure convective limit.  The empirical rule of convection states that the ratio of the surface buoyancy flux to that at the entrainment depth be a constant.  Thus the entrainment flux at the bottom of the boundary layer under such conditions should be independent of the stratification at that depth.  Without a turbulent shear term in the denominator of the bulk Richardson number calculation, the estimated boundary layer depth is too shallow and the diffusivity at the entrainment depth is too low to obtain the necessary entrainment flux.  Thus by adding a turbulent shear term proportional to the stratification in the denominator, the calculated boundary layer depth will be deeper and will lead to a high enough diffusivity to satisfy the empirical rule of convection.

====Turbulent velocity scale====

To estimate <math>w_x</math> (where <math>x</math> is <math>m</math> - momentum ''or'' <math>s</math> - any scalar) throughout the boundary layer, surface layer similarity theory is utilized.  Following an argument by [[Bibliography#Troen_1986 | Troen and Mahrt (1986)]], [[Bibliography#LargeWG_1994a | Large et al.]] estimate the velocity scale as

{| class="eqno"
|<math display="block">we_x=\frac{\kappa u_*}{\phi_x(\zeta)}</math><!--\eqno{(6)}-->
|(6)
|}

where <math>\zeta</math> is the surface layer stability parameter defined as <math>z/L</math>.  <math>\phi_x</math> is a non-dimensional flux profile which varies based on the stability of the boundary layer forcing.  The stability parameter used in this equation is assumed to vary over the entire depth of the boundary layer in stable and neutral conditions.  In unstable conditions it is assumed only to vary through the surface layer which is defined as <math> \epsilon h_{sbl} </math> (where <math>\epsilon</math> is set at 0.10) .  Beyond this depth <math>\zeta</math> is set equal to its value at <math> \epsilon h_{sbl} </math>.

The flux profiles are expressed as analytical fits to atmospheric surface boundary layer data.  In stable conditions they vary linearly with the stability parameter <math>\zeta</math>  as

{| class="eqno"
|<math display="block">\phi_x=1+5\zeta</math><!--\eqno{(7)}-->
|(7)
|}

In near-neutral unstable conditions common Businger-Dyer forms are used which match with the formulation for stable conditions at <math>\zeta=0</math>.  Near neutral conditions are defined as

<math display="block">-0.2 \leq \zeta < 0</math>

for momentum and,

<math display="block">-1.0 \leq \zeta < 0</math>

for scalars.  The non dimensional flux profiles in this regime are,

{| class="eqno"
|<math display="block">\begin{align}
\phi_m & =(1-16\zeta)^{1/4} \\
\phi_s & =(1-16\zeta)^{1/2} \end{align}
</math><!--\eqno{(8)}-->
|(8)
|}

In more unstable conditions <math>\phi_x</math> is chosen to match the Businger-Dyer forms and with the free convective limit.  Here the flux profiles are

{| class="eqno"
|<math display="block">\begin{align}
\phi_m & =(1.26-8.38\zeta)^{1/3}    \\
\phi_s & =(-28.86-98.96\zeta)^{1/3} \end{align}
</math><!--\eqno{(9)}-->
|(9)
|}

====The shape function====

The non-dimensional shape function <math>G(\sigma)</math> is a third order polynomial with coefficients chosen to match the interior viscosity at the bottom of the boundary layer and Monin-Obukhov similarity theory approaching the surface.  This function is defined as a 3rd order polynomial.

{| class="eqno"
|<math display="block">G(\sigma)=a_o+a_1\sigma+a_2\sigma^2+a_3\sigma^3</math><!--\eqno{(10)}-->
|(10)
|}

with the coefficients specified to match surface boundary conditions and to smoothly blend with the interior,

{| class="eqno"
|<math display="block">\begin{align}
  a_o & =0         \\
  a_1 & =1        \\
  a_2 & =-2+3\frac{\nu_{x}(h_{sbl})}{hw_x(1)}+\frac{\partial_x
  \nu_{x}(h)}{w_{x}(1)}+\frac{\nu_{x}(h)
  \partial_{\sigma}w_x(1)}{hw_{x}^{2}(1)}  \\
  a_3 & =1-2\frac{\nu_{x}(h_{sbl})}{hw_x(1)}-\frac{\partial_x
  \nu_{x}(h)}{w_{x}(1)}-\frac{\nu_{x}(h)
  \partial_{\sigma}w_x(1)}{hw_{x}^{2}(1)} \end{align}
</math><!--\eqno{(11)}-->
|(11)
|}

where <math>\nu_{x}(h)</math> is the viscosity calculated by the interior parameterization at the boundary layer depth.

===Countergradient flux term===

The second term of the LMD scheme's surface boundary layer formulation is the non-local transport term <math>\gamma</math> which can play a significant role in mixing during surface cooling events.  This is a redistribution term included in the tracer equation separate from the diffusion term and is written as

{| class="eqno"
|<math display="block">-\frac{\partial}{\partial z}K\gamma\,.</math><!--\eqno{(12)}-->
|(12)
|}

LMD base their formulation for non-local scalar transport on a parameterization for pure free convection from [[Bibliography#Mailhot_1982 | Mailhôt and Benoit (1982)]]. They extend this parameterization to cover any unstable surface forcing conditions to give

{| class="eqno"
|<math display="block">
  \gamma_{T}=C_s\frac{\overline{wT_0}+
  \overline{wT_R}}{w_T(\sigma)h}
</math><!--\eqno{(13)}-->
|(13)
|}

for temperature and

{| class="eqno"
|<math display="block">\gamma_S=C_s \frac{\overline{wS_0}}{w_S(\sigma)h}</math><!--\eqno{(14)}-->
|(14)
|}   

for salinity (other scalar quantities with surface fluxes can be treated similarly). LMD argue that although there is evidence of non-local transport of momentum as well, the form the term would take is unclear so they simply specify <math>\gamma_m=0</math>.

===The interior scheme===

The interior scheme of Large, McWilliams and Doney estimates the viscosity coefficient by adding the effects of several generating mechanisms:  shear mixing, double-diffusive mixing and internal wave generated mixing.

{| class="eqno"
|<math display="block">\nu_{x}(d)=\nu_{x}^s+\nu_{x}^d+\nu_{x}^w</math><!--\eqno{(15)}-->
|(15)
|}

====Shear generated mixing====

The shear mixing term is calculated using a gradient Richardson number formulation, with viscosity estimated as:

{| class="eqno"
|<math display="block">
   \nu^s_x=\begin{cases}
      \nu_0&   Ri_g<0\,, \\
      \nu_0[1-(Ri_g/Ri_0)^2]^3& 0< Ri_g<Ri_0\,,  \\
      0&   Ri_g>Ri_0\,.\end{cases}
</math><!--\eqno{(16)}-->
|(16)
|}

where <math>\nu_0</math> is <math>5.0 \times 10^{-3}</math>, <math>Ri_0 = 0.7</math>.

====Double diffusive processes====

The second component of the interior mixing parameterization represents double diffusive mixing.  From limited sources of laboratory and field data LMD parameterize the salt fingering case (<math>R_{\rho}>1.0</math>)

{| class="eqno"
|<math display="block">\begin{align}
\nu_{s}^{d}(R_{\rho}) & =
        \begin{cases}
      1\times10^{-4}[1-(\frac{(R_{\rho}-1}{R_{\rho}^0-1})^2)^{3}&
      \text{for} \; 1.0<R_{\rho}<R_{\rho}^0=1.9,\\
           0& \text{otherwise.}      \end{cases}\\      \\
\nu_{\theta}^{d}(R_{\rho}) & =0.7\nu_{s}^{d} \end{align}
</math><!--\eqno{(17)}-->
|(17)
|}

For diffusive convection (<math>0<R_{\rho}<1.0</math>) LMD suggest several formulations from the literature and choose the one with the most significant impact on mixing ([[Bibliography#Fedorov_1988 |Fedorov 1988]]).

{| class="eqno"
|<math display="block">
   \nu_{\theta}^{d}=(1.5 \times 10^{-6})(0.909 \exp(4.6 \exp[-0.54(R_{\rho}^{-1}-1)])
</math><!--\eqno{(18)}-->
|(18)
|}

for temperature.  For other scalars,

{| class="eqno"
|<math display="block">
   \nu_{s}^{d}=
        \begin{cases}
             \nu_{\theta}^{d}(1.85-0.85R_{\rho}^{-1})R_{\rho}& \text{for} 0.5<=R_{\rho}<1.0,\\
             \nu_{\theta}^{d}0.15R_{\rho}&  \text{otherwise.} \end{cases}
</math>!--\eqno{(19)}-->
|(19)
|}

====Internal wave generated mixing====

Internal wave generated mixing serves as the background mixing in the LMD scheme.  It is specified as a constant for both scalars and momentum. Eddy diffusivity is estimated based on the data of [[Bibliography#Ledwell_1993 | Ledwell et al. (1993)]], while [[Bibliography#Peters_1988 | Peters et al. (1988)]] suggest eddy viscosity should be 7 to 10 times larger than diffusivity for gradient Richardson numbers below approximately 0.7. Therefore LMD use

{| class="eqno"
|<math display="block">\begin{align}
   \nu_{m}^w & =1.0 \times 10^{-4} m^2 s^{-1} \\
   \nu_{s}^w & =1.0 \times 10^{-5} m^2 s^{-1} \end{align}
</math><!--\eqno{(20)}-->
|(20)
|}

==Mellor-Yamada 2.5==

One of the more popular closure schemes is that of [[Bibliography#MellorGL_1982a | Mellor and Yamada (1982)]]. They actually present a hierarchy of closures of increasing complexity. ROMS provides only the "Level 2.5" closure with the [[Bibliography#GalperinB_1988a | Galperin et al. (1988)]] modifications as described in [[Bibliography#AllenJS_1995 | Allen et al. (1995)]]. This closure scheme adds two prognostic equations, one for the turbulent kinetic energy (<math>{1 \over 2} q^2</math>) and one for the  turbulent kinetic energy times a length scale (<math>q^2l</math>).

The turbulent kinetic energy equation is: 

{| class="eqno"
|<math display="block">
  {D \over Dt} \left( {q^2 \over 2} \right) -
  {\partial \over \partial z} \left[ K_q {\partial \over \partial z} 
  \left( {q^2 \over 2} \right) \right] = P_s + P_b - \xi_d
</math><!--\eqno{(21)}-->
|(21)
|}

where <math>P_s</math> is the shear production, <math>P_b</math> is the buoyant production and <math>\xi_d</math> is the dissipation of turbulent kinetic energy. These terms are given by

{| class="eqno"
|<math display="block">\begin{align}
   P_s &= K_m \left[ \left( {\partial u \over \partial z }\right)^2 +
   \left( {\partial v \over \partial z} \right)^2 \right],  \\
   P_b &= -K_s N^2, \\
   \xi_d &= {q^3 \over B_1 l} \end{align}
</math><!--\eqno{(22)}-->
|(22)
|}

where <math>B_1</math> is a constant. One can also add a traditional horizontal Laplacian or biharmonic diffusion (<math>{\cal D}_q</math>) to the turbulent kinetic energy equation. The form of this equation in the model coordinates becomes

<math display="block">
  {\partial \over \partial t} \left( {H_z q^2 \over mn} \right) +
  {\partial \over \partial \xi} \left( {H_z u q^2 \over n} \right) +
  {\partial \over \partial \eta} \left( {H_z v q^2 \over m} \right) +
  {\partial \over \partial s} \left( {H_z \Omega q^2 \over mn} \right) -
  {\partial \over \partial s} \left( {K_q \over mnH_z}
  {\partial q^2 \over \partial s} \right) =
</math>

{| class="eqno"
|<math display="block">
  {2H_z K_m \over mn} \left[ \left({\partial u \over \partial z}
  \right)^2 + \left( {\partial v \over \partial z} \right)^2 \right] +
  {2H_z K_s \over mn} N^2 - {2H_z q^3 \over mnB_1 l} +
  {H_z \over mn} {\cal D}_q\,.
</math><!--\eqno{(23)}-->
|(23)
|}

The vertical boundary conditions are:
:top (<math>z = \zeta(x,y,t))</math>:
<math display="block">\begin{align} & {H_z \Omega \over mn} = 0 \\
  & {K_q \over mn H_z} \, \frac{\partial q^2}{\partial s}
    = {B_1^{2/3} \over \rho_o} \left[ \left( \tau_s^\xi \right)^2
    + \left( \tau_s^\eta \right)^2  \right] \\
  & H_z K_m \left( {\partial u \over \partial z},
    {\partial v \over \partial z} \right) = {1 \over \rho_o}
    \left( \tau_s^\xi, \tau_s^\eta \right) \\
  & H_z K_s N^2 = {Q \over \rho_o c_P} \end{align}
</math>

:and bottom (<math>z = -h(x,y)</math>):

<math display="block">\begin{align}
    &{H_z \Omega \over mn} = 0 \\
  & {K_q \over mnH_z} {\partial q^2 \over \partial s}
    = {B_1^{2/3} \over \rho_o} \left[ \left( \tau_b^\xi \right)^2
    + \left( \tau_b^\eta \right)^2  \right] \\
  & H_z K_m \left( {\partial u \over \partial z},
    {\partial v \over \partial z} \right) = {1 \over \rho_o}
    \left( \tau_b^\xi, \tau_b^\eta \right) \\
  & H_z K_s N^2 = 0 \end{align}
</math>

There is also an equation for the turbulent length scale <math>l</math>:

{| class="eqno"
|<math display="block">
  {D \over Dt} \left( {lq^2} \right) -
  {\partial \over \partial z} \left[ K_l
  {\partial lq^2 \over \partial z} 
  \right] = lE_1 ( P_s + P_b ) - {q^3 \over B_1} \tilde{W}
</math><!--\eqno{(24)}-->
|(24)
|}

where <math>\tilde{W}</math> is the wall proximity function:

{| class="eqno"
|<math display="block">\begin{align}
  \tilde{W} &= 1 + E_2 \left( {l \over kL} \right) ^2 \\
  L^{-1} &= {1 \over \zeta -z} + {1 \over H+z} \end{align}
</math><!--\eqno{(25)}-->
|(25)
|}

The form of this equation in the model coordinates becomes

<math display="block">
  {\partial \over \partial t} \left( {H_z q^2l \over mn} \right) +
  {\partial \over \partial \xi} \left( {H_z u q^2l \over n} \right) +
  {\partial \over \partial \eta} \left( {H_z v q^2l \over m} \right) +
  {\partial \over \partial s} \left( {H_z \Omega q^2l \over mn} \right) -
  {\partial \over \partial s} \left( {K_q \over mnH_z}
  {\partial q^2l \over \partial s} \right) =
</math>

{| class="eqno"
|<math display="block">
  {H_z \over mn} lE_1 ( P_s + P_b) - 
  {H_z q^3 \over mnB_1 } \tilde{W} +
  {H_z \over mn} {\cal D}_{ql}\,.
</math><!--\eqno{(26)}-->
|(26)
|}

where <math>{\cal D}_{ql}</math> is the horizontal diffusion of the quantity <math>q^2l</math>. Equations (23) and (26) are timestepped much like the model tracer equations, including an implicit solve for the vertical operations and options for centered second or fourth-order advection. They are timestepped with a predictor-corrector scheme in which the predictor step is only computing the advection.

Given these solutions for <math>q</math> and <math>l</math>, the vertical viscosity and diffusivity coefficients are:

{| class="eqno"
|<math display="block">\begin{align}
  K_m &= qlS_m + K_{m_{\rm background}} \\
  K_s &= qlS_h + K_{s_{\rm background}} \\
  K_q &= qlS_q + K_{q_{\rm background}} \end{align}
</math><!--\eqno{(27)}-->
|(27)
|}

and the stability coefficients <math>S_m</math>, <math>S_h</math> and <math>S_q</math> are found by solving

{| class="eqno"
|<math display="block">
  S_s \left[ 1 - (3A_2 B_2 + 18 A_1 A_2) G_h \right] =
  A_2 \left[ 1 - 6A_1 B_1^{-1} \right]
</math><!--\eqno{(28)}-->
|(28)
|}

{| class="eqno"
|<math display="block">
  S_m \left[ 1 - 9A_1 A_2 G_h \right] - S_s \left[ G_h ( 18 A_1^2 +
  9A_1 A_2 ) G_h \right] =
  A_1 \left[ 1 - 3C_1 - 6A_1 B_1^{-1} \right]
</math><!--\eqno{(29)}-->
|(29)
|}

{| class="eqno"
|<math display="block">
  G_h = \min ( -{l^2N^2 \over q^2}, 0.028 )\,.
</math><!--\eqno{(30)}-->
|(30)
|}

{| class="eqno"
|<math display="block">S_q = 0.41 S_m</math><!--\eqno{(31)}-->
|(31)
|}

The constants are set to <math>(A_1, A_2, B_1, B_2, C_1, E_1, E_2) = (0.92, 0.74, 16.6, 10.1, 0.08, 1.8, 1.33)</math>. The quantities <math>q^2</math> and <math>q^2l</math> are both constrained to be no smaller than <math>10^{-8}</math> while <math>l</math> is set to be no larger than <math>0.53q/N</math>.


==Generic Length Scale==

[[Bibliography#UmlaufL_2003a | Umlauf and Burchard (2003)]] have come up with a generic two-equation turbulence closure scheme which can be tuned to behave like several of the traditional schemes, including that of Mellor and Yamada 2.5 (above). This is known as the Generic Length Scale, or GLS vertical mixing scheme and was introduced to ROMS in [[Bibliography#WarnerJC_2005a | Warner et al. (2005)]]. Its parameters are set in the ROMS input file.

The first of Warner et al.'s equations is the same as equation (21) with <math>k=\frac{1}{2} q^2</math>. Their dissipation is given by

{| class="eqno"
|<math display="block">
  \epsilon = (c^0_\mu ) ^{3+p/n} k^{3/2+m/n} \psi ^{-1/n}
</math><!--\eqno{(32)}-->
|(32)
|}

where <math>\psi</math> is a generic parameter that is used to establish the turbulence length scale. The equation for <math>\psi</math> is: 

{| class="eqno"
|<math display="block">
  {D \psi \over D t} = {\partial \over \partial z} \left( K_\psi  
  {\partial \psi \over \partial z}\right) + {\psi \over k}
  (c_1 P_s + c_3 P_b - c_2 \epsilon F_{wall})
</math><!--\eqno{(33)}-->

Coefficients <math>c_1</math> and <math>c_2</math> are chosen to be consistent with observations of decaying homogeneous, isotropic turbulence. The parameter <math>c_3</math> has differing values for stable (<math>c^+_3</math>) and unstable (<math>c^-_3</math>) stratification. Also

{| class="eqno"
|<math display="block">\begin{align}
   \psi &= (c^0_\mu)^p k^m l^n \\
   l &= (c^0_\mu)^3 k^{3/2} \epsilon{-1} \end{align}
</math><!--\eqno{(34)}-->
|(34)
|}

Depending on the choice of the various parameters, these two equations can
be made to solve a variety of traditional two-equation turbulence closure
models. The list of parameters is shown in the following table and is also
given inside the comments section of the ROMS input file.

{| border="1" cellspacing="0" cellpadding="5" align="center"
|
| <math>k-kl</math> || <math>k-\epsilon</math> || <math>k-\omega</math> || gen 
|-
|
| <math>\psi = k^1 l^1</math> || <math>\psi = (c^0_\mu)^3 k^{3/2} l^{-1}</math>
| <math>\psi = (c^0_\mu)^{-1} k^{1/2} l^{-1}</math> || <math>\psi = (c^0_\mu)^{2} k^{1} l^{-2/3}</math> 
|-
|  <math>p</math> || 0.0 || 3.0 || -1.0 || 2.0
|-
|  <math>m</math> || 1.0 || 1.5 || 0.5 || 1.0
|-
|  <math>n</math> || 1.0 || -1.0 || -1.0 || -0.67
|-
|  <math>\sigma_k = {K_M \over K_k}</math> || 1.96 || 1.0 || 2.0 || 0.8
|-
|  <math>\sigma_\psi = {K_M \over K_\psi}</math> || 1.96 || 1.3 || 2.0 || 1.07
|-
| <math>c_1</math> || 0.9 || 1.44 || 0.555 || 1.0
|-
|  <math>c_2</math> || 0.52 || 1.92 || 0.833 || 1.22
|-
|  <math>c^-_3</math> || 2.5 || -0.4 || -0.6 || 0.1
|-
|  <math>c^+_3</math> || 1.0 || 1.0 || 1.0 || 1.0 
|-
|  <math>k_{min}</math> || 5.0e-6 || 7.6e-6 || 7.6e-6 || 1.0e-8 
|-
|  <math>\psi_{min}</math> || 5.0e-6 || 1.0e-12 || 1.0e-12 || 1.0e-8 
|-
|  <math>c^0_\mu</math> || 0.5544 || 0.5477 || 0.5477 || 0.5544 
|-
|}