# Vertical Mixing Parameterizations

Vertical Mixing Parameterizations

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 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 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 . 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 ( ) are expressed as the product of the length scale , a turbulent velocity scale and a non-dimensional shape function.

 $nu Subscript x Baseline equals h Subscript s b l Baseline w Subscript x Baseline left-parenthesis sigma right-parenthesis upper G Subscript x Baseline left-parenthesis sigma right-parenthesis$ (1)

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

#### Surface Boundary layer depth

The boundary layer depth is calculated as the minimum of the Ekman depth, estimated as,

 $h Subscript e Baseline equals 0.7 u Subscript asterisk Baseline slash f$ (2)

(where is the friction velocity ), the Monin-Obukhov depth:

 $upper L equals u Subscript asterisk Superscript 3 Baseline slash left-parenthesis kappa upper B Subscript f Baseline right-parenthesis$ (3)

(where is von Karman's contant and is the surface buoyancy flux), and the shallowest depth at which a critical bulk Richardson number is reached. The critical bulk Richardson number () is typically in the range 0.25--0.5. The bulk Richardson number () is calculated as:

 $upper R i Subscript b Baseline left-parenthesis z right-parenthesis equals StartFraction left-parenthesis upper B Subscript r Baseline minus upper B left-parenthesis d right-parenthesis right-parenthesis d Over vertical-bar ModifyingAbove upper V With right-arrow Subscript r Baseline minus ModifyingAbove upper V With right-arrow left-parenthesis d right-parenthesis vertical-bar squared plus upper V Subscript t Baseline Superscript 2 Baseline left-parenthesis d right-parenthesis EndFraction$ (4)

where is distance down from the surface, is the buoyancy, is the buoyancy at a near surface reference depth, is the mean horizontal velocity, the velocity at the near surface reference depth and is an estimate of the turbulent velocity contribution to velocity shear.

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

 $upper V Subscript t Superscript 2 Baseline left-parenthesis d right-parenthesis equals StartFraction upper C Subscript v Baseline left-parenthesis minus beta Subscript upper T Baseline right-parenthesis Superscript 1 slash 2 Baseline Over upper R i Subscript c Baseline kappa EndFraction left-parenthesis c Subscript s Baseline epsilon right-parenthesis Superscript negative 1 slash 2 Baseline d upper N w Subscript s Baseline$ (5)

where is the ratio of interior to at the entrainment depth, is ratio of entrainment flux to surface buoyancy flux, and are constants, and 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 (where is - momentum or - any scalar) throughout the boundary layer, surface layer similarity theory is utilized. Following an argument by Troen and Mahrt (1986), Large et al. estimate the velocity scale as

 $w e Subscript x Baseline equals StartFraction kappa u Subscript asterisk Baseline Over phi Subscript x Baseline left-parenthesis zeta right-parenthesis EndFraction$ (6)

where is the surface layer stability parameter defined as . 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 (where is set at 0.10) . Beyond this depth is set equal to its value at .

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

 $phi Subscript x Baseline equals 1 plus 5 zeta$ (7)

In near-neutral unstable conditions common Businger-Dyer forms are used which match with the formulation for stable conditions at . Near neutral conditions are defined as

for momentum and,

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

 $StartLayout 1st Row 1st Column phi Subscript m 2nd Column equals left-parenthesis 1 minus 16 zeta right-parenthesis Superscript 1 slash 4 Baseline 2nd Row 1st Column phi Subscript s 2nd Column equals left-parenthesis 1 minus 16 zeta right-parenthesis Superscript 1 slash 2 EndLayout$ (8)

In more unstable conditions is chosen to match the Businger-Dyer forms and with the free convective limit. Here the flux profiles are

 $StartLayout 1st Row 1st Column phi Subscript m 2nd Column equals left-parenthesis 1.26 minus 8.38 zeta right-parenthesis Superscript 1 slash 3 Baseline 2nd Row 1st Column phi Subscript s 2nd Column equals left-parenthesis negative 28.86 minus 98.96 zeta right-parenthesis Superscript 1 slash 3 EndLayout$ (9)

#### The shape function

The non-dimensional shape function 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.

 $upper G left-parenthesis sigma right-parenthesis equals a Subscript o Baseline plus a 1 sigma plus a 2 sigma squared plus a 3 sigma cubed$ (10)

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

 $StartLayout 1st Row 1st Column a Subscript o 2nd Column equals 0 2nd Row 1st Column a 1 2nd Column equals 1 3rd Row 1st Column a 2 2nd Column equals negative 2 plus 3 StartFraction nu Subscript x Baseline left-parenthesis h Subscript s b l Baseline right-parenthesis Over h w Subscript x Baseline left-parenthesis 1 right-parenthesis EndFraction plus StartFraction normal partial-differential Subscript x Baseline nu Subscript x Baseline left-parenthesis h right-parenthesis Over w Subscript x Baseline left-parenthesis 1 right-parenthesis EndFraction plus StartFraction nu Subscript x Baseline left-parenthesis h right-parenthesis normal partial-differential Subscript sigma Baseline w Subscript x Baseline left-parenthesis 1 right-parenthesis Over h w Subscript x Superscript 2 Baseline left-parenthesis 1 right-parenthesis EndFraction 4th Row 1st Column a 3 2nd Column equals 1 minus 2 StartFraction nu Subscript x Baseline left-parenthesis h Subscript s b l Baseline right-parenthesis Over h w Subscript x Baseline left-parenthesis 1 right-parenthesis EndFraction minus StartFraction normal partial-differential Subscript x Baseline nu Subscript x Baseline left-parenthesis h right-parenthesis Over w Subscript x Baseline left-parenthesis 1 right-parenthesis EndFraction minus StartFraction nu Subscript x Baseline left-parenthesis h right-parenthesis normal partial-differential Subscript sigma Baseline w Subscript x Baseline left-parenthesis 1 right-parenthesis Over h w Subscript x Superscript 2 Baseline left-parenthesis 1 right-parenthesis EndFraction EndLayout$ (11)

where is the viscosity calculated by the interior parameterization at the boundary layer depth.

The second term of the LMD scheme's surface boundary layer formulation is the non-local transport term 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

 $minus StartFraction normal partial-differential Over normal partial-differential z EndFraction upper K gamma period$ (12)

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

 $gamma Subscript upper T Baseline equals upper C Subscript s Baseline StartFraction ModifyingAbove w upper T 0 With bar plus ModifyingAbove w upper T Subscript upper R Baseline With bar Over w Subscript upper T Baseline left-parenthesis sigma right-parenthesis h EndFraction$ (13)

for temperature and

 $gamma Subscript upper S Baseline equals upper C Subscript s Baseline StartFraction ModifyingAbove w upper S 0 With bar Over w Subscript upper S Baseline left-parenthesis sigma right-parenthesis h EndFraction$ (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 .

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

 $nu Subscript x Baseline left-parenthesis d right-parenthesis equals nu Subscript x Superscript s Baseline plus nu Subscript x Superscript d Baseline plus nu Subscript x Superscript w$ (15)

#### Shear generated mixing

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

 $nu Subscript x Superscript s Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column nu 0 2nd Column upper R i Subscript g Baseline less-than 0 comma 2nd Row 1st Column nu 0 left-bracket 1 minus left-parenthesis upper R i Subscript g Baseline slash upper R i 0 right-parenthesis squared right-bracket cubed 2nd Column 0 less-than upper R i Subscript g Baseline less-than upper R i 0 comma 3rd Row 1st Column 0 2nd Column upper R i Subscript g Baseline greater-than upper R i 0 period EndLayout$ (16)

where is , .

#### 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 ()

 $StartLayout 1st Row 1st Column nu Subscript s Superscript d Baseline left-parenthesis upper R Subscript rho Baseline right-parenthesis 2nd Column equals StartLayout Enlarged left-brace 1st Row 1st Column 1 times 10 Superscript negative 4 Baseline left-bracket 1 minus left-parenthesis StartFraction left-parenthesis upper R Subscript rho Baseline minus 1 Over upper R Subscript rho Superscript 0 Baseline minus 1 EndFraction right-parenthesis squared right-parenthesis cubed 2nd Column for 1.0 less-than upper R Subscript rho Baseline less-than upper R Subscript rho Superscript 0 Baseline equals 1.9 comma 2nd Row 1st Column 0 2nd Column otherwise period EndLayout 2nd Row 1st Column Blank 3rd Row 1st Column nu Subscript theta Superscript d Baseline left-parenthesis upper R Subscript rho Baseline right-parenthesis 2nd Column equals 0.7 nu Subscript s Superscript d EndLayout$ (17)

For diffusive convection () LMD suggest several formulations from the literature and choose the one with the most significant impact on mixing (Fedorov 1988).

 $nu Subscript theta Superscript d Baseline equals left-parenthesis 1.5 times 10 Superscript negative 6 Baseline right-parenthesis left-parenthesis 0.909 exp left-parenthesis 4.6 exp left-bracket minus 0.54 left-parenthesis upper R Subscript rho Superscript negative 1 Baseline minus 1 right-parenthesis right-bracket right-parenthesis$ (18)

for temperature. For other scalars,

 $nu Subscript s Superscript d Baseline equals StartLayout Enlarged left-brace 1st Row 1st Column nu Subscript theta Superscript d Baseline left-parenthesis 1.85 minus 0.85 upper R Subscript rho Superscript negative 1 Baseline right-parenthesis upper R Subscript rho Baseline 2nd Column for 0.5 <= upper R Subscript rho Baseline less-than 1.0 comma 2nd Row 1st Column nu Subscript theta Superscript d Baseline 0.15 upper R Subscript rho Baseline 2nd Column otherwise period EndLayout$!--\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 Ledwell et al. (1993), while 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

 $StartLayout 1st Row 1st Column nu Subscript m Superscript w 2nd Column equals 1.0 times 10 Superscript negative 4 Baseline m squared s Superscript negative 1 Baseline 2nd Row 1st Column nu Subscript s Superscript w 2nd Column equals 1.0 times 10 Superscript negative 5 Baseline m squared s Superscript negative 1 EndLayout$ (20)

One of the more popular closure schemes is that of Mellor and Yamada (1982). They actually present a hierarchy of closures of increasing complexity. ROMS provides only the "Level 2.5" closure with the Galperin et al. (1988) modifications as described in Allen et al. (1995). This closure scheme adds two prognostic equations, one for the turbulent kinetic energy () and one for the turbulent kinetic energy times a length scale ().

The turbulent kinetic energy equation is:

 $StartFraction upper D Over upper D t EndFraction left-parenthesis StartFraction q squared Over 2 EndFraction right-parenthesis minus StartFraction normal partial-differential Over normal partial-differential z EndFraction left-bracket upper K Subscript q Baseline StartFraction normal partial-differential Over normal partial-differential z EndFraction left-parenthesis StartFraction q squared Over 2 EndFraction right-parenthesis right-bracket equals upper P Subscript s Baseline plus upper P Subscript b Baseline minus xi Subscript d$ (21)

where is the shear production, is the buoyant production and is the dissipation of turbulent kinetic energy. These terms are given by

 $StartLayout 1st Row 1st Column upper P Subscript s 2nd Column equals upper K Subscript m Baseline left-bracket left-parenthesis StartFraction normal partial-differential u Over normal partial-differential z EndFraction right-parenthesis squared plus left-parenthesis StartFraction normal partial-differential v Over normal partial-differential z EndFraction right-parenthesis squared right-bracket comma 2nd Row 1st Column upper P Subscript b 2nd Column equals minus upper K Subscript s Baseline upper N squared comma 3rd Row 1st Column xi Subscript d 2nd Column equals StartFraction q cubed Over upper B 1 l EndFraction EndLayout$ (22)

where is a constant. One can also add a traditional horizontal Laplacian or biharmonic diffusion () to the turbulent kinetic energy equation. The form of this equation in the model coordinates becomes

 $StartFraction 2 upper H Subscript z Baseline upper K Subscript m Baseline Over m n EndFraction left-bracket left-parenthesis StartFraction normal partial-differential u Over normal partial-differential z EndFraction right-parenthesis squared plus left-parenthesis StartFraction normal partial-differential v Over normal partial-differential z EndFraction right-parenthesis squared right-bracket plus StartFraction 2 upper H Subscript z Baseline upper K Subscript s Baseline Over m n EndFraction upper N squared minus StartFraction 2 upper H Subscript z Baseline q cubed Over m n upper B 1 l EndFraction plus StartFraction upper H Subscript z Baseline Over m n EndFraction script upper D Subscript q Baseline period$ (23)

The vertical boundary conditions are:

top (:

and bottom ():

There is also an equation for the turbulent length scale :

 $StartFraction upper D Over upper D t EndFraction left-parenthesis l q squared right-parenthesis minus StartFraction normal partial-differential Over normal partial-differential z EndFraction left-bracket upper K Subscript l Baseline StartFraction normal partial-differential l q squared Over normal partial-differential z EndFraction right-bracket equals l upper E 1 left-parenthesis upper P Subscript s Baseline plus upper P Subscript b Baseline right-parenthesis minus StartFraction q cubed Over upper B 1 EndFraction upper W overTilde$ (24)

where is the wall proximity function:

 $StartLayout 1st Row 1st Column upper W overTilde 2nd Column equals 1 plus upper E 2 left-parenthesis StartFraction l Over k upper L EndFraction right-parenthesis squared 2nd Row 1st Column upper L Superscript negative 1 2nd Column equals StartFraction 1 Over zeta minus z EndFraction plus StartFraction 1 Over upper H plus z EndFraction EndLayout$ (25)

The form of this equation in the model coordinates becomes

 $StartFraction upper H Subscript z Baseline Over m n EndFraction l upper E 1 left-parenthesis upper P Subscript s Baseline plus upper P Subscript b Baseline right-parenthesis minus StartFraction upper H Subscript z Baseline q cubed Over m n upper B 1 EndFraction upper W overTilde plus StartFraction upper H Subscript z Baseline Over m n EndFraction script upper D Subscript q l Baseline period$ (26)

where is the horizontal diffusion of the quantity . 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 and , the vertical viscosity and diffusivity coefficients are:

 $StartLayout 1st Row 1st Column upper K Subscript m 2nd Column equals q l upper S Subscript m Baseline plus upper K Subscript m Sub Subscript normal b normal a normal c normal k normal g normal r normal o normal u normal n normal d Baseline 2nd Row 1st Column upper K Subscript s 2nd Column equals q l upper S Subscript h Baseline plus upper K Subscript s Sub Subscript normal b normal a normal c normal k normal g normal r normal o normal u normal n normal d Baseline 3rd Row 1st Column upper K Subscript q 2nd Column equals q l upper S Subscript q Baseline plus upper K Subscript q Sub Subscript normal b normal a normal c normal k normal g normal r normal o normal u normal n normal d EndLayout$ (27)

and the stability coefficients , and are found by solving

 $upper S Subscript s Baseline left-bracket 1 minus left-parenthesis 3 upper A 2 upper B 2 plus 18 upper A 1 upper A 2 right-parenthesis upper G Subscript h Baseline right-bracket equals upper A 2 left-bracket 1 minus 6 upper A 1 upper B 1 Superscript negative 1 Baseline right-bracket$ (28)
 $upper S Subscript m Baseline left-bracket 1 minus 9 upper A 1 upper A 2 upper G Subscript h Baseline right-bracket minus upper S Subscript s Baseline left-bracket upper G Subscript h Baseline left-parenthesis 18 upper A 1 squared plus 9 upper A 1 upper A 2 right-parenthesis upper G Subscript h Baseline right-bracket equals upper A 1 left-bracket 1 minus 3 upper C 1 minus 6 upper A 1 upper B 1 Superscript negative 1 Baseline right-bracket$ (29)
 $upper G Subscript h Baseline equals min left-parenthesis minus StartFraction l squared upper N squared Over q squared EndFraction comma 0.028 right-parenthesis period$ (30)
 $upper S Subscript q Baseline equals 0.41 upper S Subscript m$ (31)

The constants are set to . The quantities and are both constrained to be no smaller than while is set to be no larger than .

## Generic Length Scale

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 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 . Their dissipation is given by

 $epsilon equals left-parenthesis c Subscript mu Superscript 0 Baseline right-parenthesis Superscript 3 plus p slash n Baseline k Superscript 3 slash 2 plus m slash n Baseline psi Superscript negative 1 slash n$ (32)

where is a generic parameter that is used to establish the turbulence length scale. The equation for is:

Coefficients and are chosen to be consistent with observations of decaying homogeneous, isotropic turbulence. The parameter has differing values for stable () and unstable () stratification. Also

 $StartLayout 1st Row 1st Column psi 2nd Column equals left-parenthesis c Subscript mu Superscript 0 Baseline right-parenthesis Superscript p Baseline k Superscript m Baseline l Superscript n Baseline 2nd Row 1st Column l 2nd Column equals left-parenthesis c Subscript mu Superscript 0 Baseline right-parenthesis cubed k Superscript 3 slash 2 Baseline epsilon negative 1 EndLayout$ (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.

 $k minus k l$ $k minus epsilon$ $k minus omega$ gen $psi equals k Superscript 1 Baseline l Superscript 1$ $psi equals left-parenthesis c Subscript mu Superscript 0 Baseline right-parenthesis cubed k Superscript 3 slash 2 Baseline l Superscript negative 1$ $psi equals left-parenthesis c Subscript mu Superscript 0 Baseline right-parenthesis Superscript negative 1 Baseline k Superscript 1 slash 2 Baseline l Superscript negative 1$ $psi equals left-parenthesis c Subscript mu Superscript 0 Baseline right-parenthesis squared k Superscript 1 Baseline l Superscript negative 2 slash 3$ $p$ 0.0 3.0 -1.0 2.0 $m$ 1.0 1.5 0.5 1.0 $n$ 1.0 -1.0 -1.0 -0.67 $sigma Subscript k Baseline equals StartFraction upper K Subscript upper M Baseline Over upper K Subscript k Baseline EndFraction$ 1.96 1.0 2.0 0.8 $sigma Subscript psi Baseline equals StartFraction upper K Subscript upper M Baseline Over upper K Subscript psi Baseline EndFraction$ 1.96 1.3 2.0 1.07 $c 1$ 0.9 1.44 0.555 1.0 $c 2$ 0.52 1.92 0.833 1.22 $c 3 Superscript minus$ 2.5 -0.4 -0.6 0.1 $c 3 Superscript plus$ 1.0 1.0 1.0 1.0 $k Subscript m i n$ 5.0e-6 7.6e-6 7.6e-6 1.0e-8 $psi Subscript m i n$ 5.0e-6 1.0e-12 1.0e-12 1.0e-8 $c Subscript mu Superscript 0$ 0.5544 0.5477 0.5477 0.5544