Turbulence closures

The turbulence closure selected by the user determines the form of stress divergence $\mathbf{\nabla }\mathbf{\cdot }\mathit{\tau }$ and diffusive flux divergence $\mathbf{\nabla }\mathbf{\cdot }{\mathit{q}}_c$ in the momentum and tracer conservation equations.

Constant isotropic diffusivity

In a constant isotropic diffusivity model, the kinematic stress tensor is defined

$${\tau }_{ij}=-\nu {\mathrm{\Sigma }}_{ij},$$

where $\nu$ is a constant viscosity and ${\mathrm{\Sigma }}_{ij}\equiv \frac{1}{2}(v_{i,j}+v_{j,i})$ is the strain-rate tensor. The divergence of $\mathit{\tau }$ is then

$$\mathbf{\nabla }\mathbf{\cdot }\mathit{\tau }=-\nu {\mathrm{\nabla }}^2\mathit{v}.$$

Similarly, the diffusive tracer flux is ${\mathit{q}}_c=-\kappa \mathbf{\nabla }c$ for tracer diffusivity $\kappa$, and the diffusive tracer flux divergence is

$$\mathbf{\nabla }\mathbf{\cdot }{\mathit{q}}_c=-\kappa {\mathrm{\nabla }}^{2}c.$$

Each tracer may have a unique diffusivity $\kappa$.

Constant anisotropic diffusivity

A constant anisotropic diffusivity implies a constant tensor diffusivity ${\nu }_{jk}$ and stress ${\mathit{\tau }}_{ij}={\nu }_{jk}{u}_{i,k}$ with non-zero components ${\nu }_{11}={\nu }_{22}={\nu }_h$ and ${\nu }_{33}={\nu }_z$. With this form the kinematic stress divergence becomes

$$\mathbf{\nabla }\mathbf{\cdot }\mathit{\tau }=-[{\nu }_h({\partial }_x^2+{\partial }_y^2)+{\nu }_v{\partial }_z^2]\mathit{v},$$

and diffusive flux divergence

$$\mathbf{\nabla }\mathbf{\cdot }{\mathit{q}}_c=-[{\kappa }_h({\partial }_x^2+{\partial }_y^2)+{\kappa }_v{\partial }_z^2]c,$$

in terms of the horizontal viscosities and diffusivities, ${\nu }_h$ and ${\kappa }_h$, and the vertical viscosity and diffusivities, ${\nu }_z$ and ${\kappa }_z$. Each tracer may have a unique diffusivity components ${\kappa }_h$ and ${\kappa }_v$.

Scalar biharmonic diffusivity

A constant biharmonic diffusivity implies a constant tensor diffusivity ${\nu }_{jk}$ and stress${\mathit{\tau }}_{ij}={\nu }_{jk}{\partial }_k^3{u}_i$ with non-zero components ${\nu }_{11}={\nu }_{22}={\nu }_h$ and ${\nu }_{33}={\nu }_z$.

With this form the kinematic stress divergence becomes

$$\mathbf{\nabla }\mathbf{\cdot }\mathit{\tau }=-[{\nu }_h{({\partial }_x^2+{\partial }_y^2)}^2+{\nu }_v{\partial }_z^4]\mathit{v},$$

and diffusive flux divergence

$$\mathbf{\nabla }\mathbf{\cdot }{\mathit{q}}_c=-[{\kappa }_h{({\partial }_x^2+{\partial }_y^2)}^2+{\kappa }_v{\partial }_z^4]c,$$

in terms of the horizontal biharmonic viscosities and diffusivities, ${\nu }_h$ and ${\kappa }_h$, and the vertical biharmonic viscosity and diffusivities, ${\nu }_z$ and ${\kappa }_z$. Each tracer may have a unique diffusivity components ${\kappa }_h$ and ${\kappa }_z$.

Smagorinsky-Lilly turbulence closure

In the turbulence closure proposed by Lilly (1962) and Smagorinsky (1963), the subgrid stress associated with unresolved turbulent motions is modeled diffusively via

$${\tau }_{ij}=-2{\nu }_e{\mathrm{\Sigma }}_{ij},$$

where ${\mathrm{\Sigma }}_{ij}=\frac{1}{2}(v_{i,j}+v_{j,i})$ is the resolved strain rate. The eddy viscosity is given by

$$\begin{array}{}\text{(1)}& {\nu }_e={\left(C{\mathrm{\Delta }}_f\right)}^2\sqrt{{\mathrm{\Sigma }}^2}\varsigma (N^2/{\mathrm{\Sigma }}^2),\end{array}$$

where ${\mathrm{\Delta }}_f$ is the "filter width" associated with the finite volume grid spacing and $C$ is a user-specified model constant, ${\mathrm{\Sigma }}^2\equiv {\mathrm{\Sigma }}_{ij}{\mathrm{\Sigma }}_{ij}$. The factor $\varsigma (N^2/{\mathrm{\Sigma }}^2)$ reduces ${\nu }_e$ in regions of strong stratification via

$$\varsigma (N^2/{\mathrm{\Sigma }}^2)=\sqrt{1-min(1,C_b{N}^2/{\mathrm{\Sigma }}^2)},$$

where $N^2=max(0,{\partial }_{z}b)$ is the squared buoyancy frequency for stable stratification with ${\partial }_{z}b>0$ and $C_b$ is a user-specified constant. Lilly (1962) proposed $C_b=1/Pr$, where $Pr$ is a turbulent Prandtl number. The filter width for the Smagorinsky-Lilly closure is

$${\mathrm{\Delta }}_f(\mathit{x})={\left(\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\right)}^{1/3},$$

where $\mathrm{\Delta }x$, $\mathrm{\Delta }y$, and $\mathrm{\Delta }z$ are the grid spacing in the $\hat{\mathit{x}}$, $\hat{\mathit{y}}$, and $\hat{\mathit{z}}$ directions at location $\mathit{x}=(x,y,z)$.

The effect of subgrid turbulence on tracer mixing is also modeled diffusively via

$${\mathit{q}}_c=-{\kappa }_e\mathbf{\nabla }c,$$

where the eddy diffusivity ${\kappa }_e$ is

$${\kappa }_e=\frac{{\nu }_e}{Pr}+\kappa .$$

Both $Pr$ and $\kappa$ may be set independently for each tracer.

Anisotropic minimum dissipation (AMD) turbulence closure

The anisotropic minimum dissipation (AMD) model proposed by Verstappen (2018) and was described and tested by Vreugdenhil and Taylor (2018). The AMD model uses an eddy diffusivity hypothesis similar the Smagorinsky-Lilly model. In the AMD model, the eddy viscosity and diffusivity for each tracer are defined in terms of eddy viscosity and diffusivity predictors ${\nu }_e^{†}$ and ${\kappa }_e^{†}$, such that

$${\nu }_e=max(0,{\nu }_e^{†})+\nu \text{and}{\kappa }_e=max(0,{\kappa }_e^{†})+\kappa ,$$

to ensure that ${\nu }_e\ge 0$ and ${\kappa }_e\ge 0$, where $\nu$ and $\kappa$ are the constant isotropic background viscosity and diffusivities for each tracer. The eddy viscosity predictor is

$$\begin{array}{}\text{(2)}& {\nu }_e^{†}=C{\mathrm{\Delta }}_f^2\frac{({\hat{\partial }}_k{\hat{v}}_i)({\hat{\partial }}_k{\hat{v}}_j){\hat{\mathrm{\Sigma }}}_{ij}+C_b{\hat{\delta }}_{i3}({\hat{\partial }}_k\widehat{v_i})({\hat{\partial }}_{k}b)}{({\hat{\partial }}_l{\hat{v}}_m)({\hat{\partial }}_l{\hat{v}}_m)},\end{array}$$

while the eddy diffusivity predictor for tracer $c$ is

$$\begin{array}{}\text{(3)}& {\kappa }_e^{†}=C{\mathrm{\Delta }}_f^2\frac{({\hat{\partial }}_k{\hat{v}}_i)({\hat{\partial }}_{k}c)({\hat{\partial }}_{i}c)}{({\hat{\partial }}_{l}c)({\hat{\partial }}_{l}c)}.\end{array}$$

In the definitions of the eddy viscosity and eddy diffusivity predictor, $C$ and $C_b$ are user-specified model constants, ${\mathrm{\Delta }}_f$ is a "filter width" associated with the finite volume grid spacing, and the hat decorators on partial derivatives, velocities, and the Kronecker delta ${\hat{\delta }}_{i3}$ are defined such that

$${\hat{\partial }}_i\equiv {\mathrm{\Delta }}_i{\partial }_i,{\hat{v}}_i(x,t)\equiv \frac{v_i(x,t)}{{\mathrm{\Delta }}_i},\text{and}{\hat{\delta }}_{i3}\equiv \frac{{\delta }_{i3}}{{\mathrm{\Delta }}_3}.$$

A velocity gradient, for example, is therefore ${\hat{\partial }}_i{\hat{v}}_j(x,t)=\frac{{\mathrm{\Delta }}_i}{{\mathrm{\Delta }}_j}{\partial }_i{v}_j(x,t)$, while the normalized strain tensor is

$${\hat{\mathrm{\Sigma }}}_{ij}=\frac{1}{2}[{\hat{\partial }}_i{\hat{v}}_j(x,t)+{\hat{\partial }}_j{\hat{v}}_i(x,t)].$$

The filter width ${\mathrm{\Delta }}_f$ in that appears in the viscosity and diffusivity predictors is taken as the square root of the harmonic mean of the squares of the filter widths in each direction:

$$\frac{1}{{\mathrm{\Delta }}_f^2}=\frac{1}{3}(\frac{1}{\mathrm{\Delta }{x}^2}+\frac{1}{\mathrm{\Delta }{y}^2}+\frac{1}{\mathrm{\Delta }{z}^2}).$$

The constant $C_b$ permits the "buoyancy modification" term it multiplies to be omitted from a calculation. By default we use the model constants $C=1/12$ and $C_b=0$.

Convective adjustment vertical diffusivity

This closure aims to model the enhanced mixing that occurs due to convection. At every point and for every time instance, the closure diagnoses the gravitational stability of the fluid and applies the vertical diffusivities (i) background_νz to u, v and background_κz to all tracers if the fluid is gravitationally neutral or stable with $\partial b/\partial z\ge 0$, or (ii) convective_νz and convective_κz if $\partial b/\partial z<0$.

This closure is a plausible model for convection if convective_κz $\gg$ background_κz and convective_νz $\gg$ background_νz.