Turbulence closures

The turbulence closure selected by the user determines the form of stress divergence $\bm{\nabla} \bm{\cdot} \bm{\tau}$ and diffusive flux divergence $\bm{\nabla} \bm{\cdot} \bm{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 \Sigma_{ij} \, ,\]

where $\nu$ is a constant viscosity and $\Sigma_{ij} \equiv \tfrac{1}{2} \left ( u_{i, j} + u_{j, i} \right )$ is the strain-rate tensor. The divergence of $\bm{\tau}$ is then

\[\bm{\nabla} \bm{\cdot} \bm{\tau} = -\nu \bm{\nabla}^2 \bm{u} \, .\]

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

\[\bm{\nabla} \bm{\cdot} \bm{q}_c = - \kappa \bm{\nabla}^2 c \, .\]

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

Constant anisotropic diffusivity

In Oceananigans.jl, a constant anisotropic diffusivity implies a constant tensor diffusivity $\nu_{j k}$ and stress $\bm{\tau}_{ij} = \nu_{j k} u_{i, k}$ with non-zero components $\nu_{11} = \nu_{22} = \nu_h$ and $\nu_{33} = \nu_v$. With this form the kinematic stress divergence becomes

\[\bm{\nabla} \bm{\cdot} \bm{\tau} = - \left [ \nu_h \left ( \partial_x^2 + \partial_y^2 \right ) + \nu_v \partial_z^2 \right ] \bm{u} \, ,\]

and diffusive flux divergence

\[\bm{\nabla} \bm{\cdot} \bm{q}_c = - \left [ \kappa_{h} \left ( \partial_x^2 + \partial_y^2 \right ) + \kappa_{v} \partial_z^2 \right ] c \, .\]

in terms of the horizontal viscosities and diffusivities $\nu_h$ and $\kappa_{h}$ and the vertical viscosity and diffusivities $\nu_v$ and $\kappa_{v}$. Each tracer may have a unique diffusivity components $\kappa_h$ and $\kappa_v$.

Constant anisotropic biharmonic diffusivity

In Oceananigans.jl, a constant anisotropic biharmonic diffusivity implies a constant tensor diffusivity $\nu_{j k}$ and stress $\bm{\tau}_{ij} = \nu_{j k} \partial_k^3 u_i$ with non-zero components $\nu_{11} = \nu_{22} = \nu_h$ and $\nu_{33} = \nu_v$. With this form the kinematic stress divergence becomes

\[\bm{\nabla} \bm{\cdot} \bm{\tau} = - \left [ \nu_h \left ( \partial_x^2 + \partial_y^2 \right )^2 + \nu_v \partial_z^4 \right ] \bm{u} \, ,\]

and diffusive flux divergence

\[\bm{\nabla} \bm{\cdot} \bm{q}_c = - \left [ \kappa_{h} \left ( \partial_x^2 + \partial_y^2 \right )^2 + \kappa_{v} \partial_z^4 \right ] c \, .\]

in terms of the horizontal biharmonic viscosities and diffusivities $\nu_h$ and $\kappa_{h}$ and the vertical biharmonic viscosity and diffusivities $\nu_v$ and $\kappa_{v}$. Each tracer may have a unique diffusivity components $\kappa_h$ and $\kappa_v$.

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} = \nu_e \Sigma_{ij} \, ,\]

where $\Sigma_{ij} = \tfrac{1}{2} \left ( u_{i, j} + u_{j, i} \right )$ is the resolved strain rate. The eddy viscosity is given by

\[ \nu_e = \left ( C \Delta_f \right )^2 \sqrt{ \Sigma^2 } \, \Upsilon(Ri) + \nu \, , \tag{eq:smagorinsky-viscosity}\]

where $\Delta_f$ is the "filter width" associated with the finite volume grid spacing, $C$ is a user-specified model constant, $\Sigma^2 \equiv \Sigma_{ij} \Sigma{ij}$, and $\nu$ is a constant isotropic background viscosity. The factor $\Upsilon(Ri)$ reduces $\nu_e$ in regions of strong stratification where the resolved gradient Richardson number $Ri \equiv N^2 / \Sigma^2$ is large via

\[ \Upsilon(Ri) = \sqrt{1 - \min \left ( 1, C_b N^2 / \Sigma^2 \right )} \, ,\]

where $N^2 = \max \left (0, \partial_z b \right )$ is the squared buoyancy frequency for stable stratification with $\partial_z b > 0$ and $C_b$ is a user-specified constant. Roughly speaking, the filter width for the Smagorinsky-Lilly closure is taken as

\[\Delta_f(\bm{x}) = \left ( \Delta x \Delta y \Delta z \right)^{1/3} \, ,\]

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

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

\[\bm{q}_c = \kappa_e \bm{\nabla} c \, ,\]

where the eddy diffusivity $\kappa_e$ is

\[\kappa_e = \frac{\nu_e - \nu}{Pr} + \kappa \, ,\]

where $Pr$ is a turbulent Prandtl number and $\kappa$ is a constant isotropic background diffusivity. Both $Pr$ and $\kappa$ may be set independently for each tracer.

Anisotropic minimum dissipation (AMD) turbulence closure

Oceananigans.jl uses the anisotropic minimum dissipation (AMD) model proposed by Verstappen18 and described and tested by Vreugdenhil18. 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^\dagger$ and $\kappa_e^\dagger$, such that

\[ \nu_e = \max \left ( 0, \nu_e^\dagger \right ) + \nu \quad \text{and} \quad \kappa_e = \max \left ( 0, \kappa_e^\dagger \right ) + \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

\[ \tag{eq:nu-dagger} \nu_e^\dagger = -(C \Delta_f)^2 \frac {\left( \hat{\partial}_k \hat{u}_i \right) \left( \hat{\partial}_k \hat{u}_j \right) \hat{\Sigma}_{ij} + C_b \hat{\delta}_{i3} \left( \hat{\partial}_k \hat{u_i} \right) \hat{\partial}_k b} {\left( \hat{\partial}_l \hat{u}_m \right) \left( \hat{\partial}_l \hat{u}_m \right)}\]

while the eddy diffusivity predictor for tracer $c$ is

\[ \tag{eq:kappa-dagger} \kappa_e^\dagger = -(C \Delta_f)^2 \frac {\left( \hat{\partial}_k \hat{u}_i \right) \left( \hat{\partial}_k c \right) \hat{\partial}_i c} {\left( \hat{\partial}_l c \right) \left( \hat{\partial}_l c \right)} \, .\]

In the definitions of the eddy viscosity and eddy diffusivity predictor, $C$ and $C_b$ are user-specified model constants, $\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 \Delta_i \partial_i, \qquad \hat{u}_i(x, t) \equiv \frac{u_i(x, t)}{\Delta_i}, \quad \text{and} \quad \hat{\delta}_{i3} \equiv \frac{\delta_{i3}}{\Delta_3} \, .\]

A velocity gradient, for example, is therefore $\hat{\partial}_i \hat{u}_j(x, t) = \frac{\Delta_i}{\Delta_j} \partial_i u_j(x, t)$, while the normalized strain tensor is

\[ \hat{\Sigma}_{ij} = \frac{1}{2} \left[ \hat{\partial}_i \hat{u}_j(x, t) + \hat{\partial}_j \hat{u}_i(x, t) \right] \, .\]

The filter width $\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}{\Delta_f^2} = \frac{1}{3} \left( \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2} \right) \, .\]

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