Navier-Stokes and tracer conservation equations

Oceananigans.jl solves the incompressible Navier-Stokes equations and an arbitrary number of tracer conservation equations. Physics associated with individual terms in the momentum and tracer conservation equations –- the background rotation rate of the equation's reference frame, gravitational effects associated with buoyant tracers under the Boussinesq approximation^[1], generalized stresses and tracer fluxes associated with viscous and diffusive physics, and arbitrary "forcing functions" –- are determined by the whims of the user.

Coordinate system and notation

Oceananigans.jl is formulated in a Cartesian coordinate system $\bm{x} = (x, y, z)$ with unit vectors $\bm{\hat x}$, $\bm{\hat y}$, and $\bm{\hat z}$, where $\bm{\hat x}$ points east, $\bm{\hat y}$ points north, and $\bm{\hat z}$ points 'upward', opposite the direction of gravitational acceleration. We denote time with $t$, partial derivatives with respect to time $t$ or a coordinate $x$ with $\partial_t$ or $\partial_x$, and denote the gradient operator $\bm{\nabla} \equiv \partial_x \bm{\hat x} + \partial_y \bm{\hat y} + \partial_z \bm{\hat z}$. We use $u$, $v$, and $w$ to denote the east, north, and vertical velocity components, such that $\bm{u} = u \bm{\hat x} + v \bm{\hat y} + w \bm{\hat z}$.

The Boussinesq approximation

In Oceananigans.jl the fluid density $\rho$ is, in general, decomposed into three components:

where $\rho_0$ is a constant 'reference' density, $\rho_*(z)$ is a background density profile typically associated with the hydrostatic compression of seawater in the deep ocean, and $\rho'(\bm{x}, t)$ is the dynamic component of density corresponding to inhomogeneous distributions of a buoyant tracer such as temperature or salinity. The fluid buoyancy, associated with the buoyant acceleration of fluid, is defined in terms of $\rho'$ as

where $g$ is gravitational acceleration.

The Boussinesq approximation is valid when $\rho_* + \rho' \ll \rho_0$, which implies the fluid is approximately incompressible^[2] In this case, the mass conservation equation reduces to the continuity equation

The momentum conservation equation

Under the Boussinesq approximation and in terms of the buoyancy $b$, the equations expressing conservation of momentum in a rotating fluid with an inhomogeneous density distribution are

where $\bm{\tau}$ is the kinematic stress tensor, $\bm{F_u}$ denotes an internal forcing of the velocity field $\bm{u}$, $\phi$ is the potential associated with kinematic and constant hydrostatic contributions to pressure, and $\bm{f}$ is Coriolis parameter, or the background vorticity associated with the specified rate of rotation of the frame of reference.

The tracer conservation equation

The conservation law for tracers in Oceananigans.jl is

where $\bm{q}_c$ is the diffusive flux of $c$ and $F_c$ is an arbitrary source term. Oceananigans.jl permits arbitrary tracers and thus an arbitrary number of tracer equations to be solved simultaneously with the momentum equations.

Buoyancy model and equations of state

The buoyancy model determines the relationship between tracers and the buoyancy $b$ in the momentum equation.

Buoyancy tracer

The simplest buoyancy model uses buoyancy $b$ itself as a tracer: $b$ obeys the tracer conservation equation and is used directly in the momentum equations in the momentum equation.

Seawater buoyancy

For seawater buoyancy is, in general, modeled as a function of conservative temperature $\theta$, absolute salinity $S$, and depth below the ocean surface $d$ via

where $g$ is gravitational acceleration, $\rho_0$ is the reference density. The function $\rho'(\theta, S, d)$ in the seawater buoyancy relationship that links conservative temperature, salinity, and depth to the density perturbation is called the equation of state. Both $\theta$ and $S$ obey the tracer conservation equation.

Linear equation of state

Buoyancy is determined from a linear equation of state via

where $g$ is gravitational acceleration, $\alpha_\theta$ is the thermal expansion coefficient, and $\beta_S$ is the haline contraction coefficient.

Nonlinear equation of state

Buoyancy is determined by the simplified equations of state introduced by Roquet et al (2015).

Coriolis forces

The Coriolis model controls the manifestation of the term $\bm{f} \times \bm{u}$ in the momentum equation.

The "$f$-plane" approximation

Under an $f$-plane approximation^[3] the reference frame in which the momentum and tracer equations are are solved rotates at a constant rate around a vertical axis, such that

where $f$ is constant and determined by the user.

The $\beta$-plane approximation

Under the $\beta$-plane approximation, the rotation axis is vertical as for the $f$-plane approximation, but $f$ is expanded in a Taylor series around a central latitude such that

where $f_0$ is the planetary vorticity at some central latitude, and $\beta$ is the planetary vorticity gradient. The $\beta$-plane model is not periodic in $y$ and thus can be used only in domains that are bounded in the $y$-direction.

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 defned

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

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

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

and diffusive flux divergence

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

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

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

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

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

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

where the eddy diffusivity $\kappa_e$ is

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 \emph{predictors} $\nu_e^\dagger$ and $\kappa_e^\dagger$, such that

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

while the eddy diffusivity predictor for tracer $c$ is

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

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

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:

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

Boundary conditions

In Oceananigans.jl the user may impose \textit{no-penetration}, \textit{flux}, \textit{gradient} (Neumann), and \textit{value} (Dirichlet) boundary conditions in bounded, non-periodic directions. Note that the only boundary condition available for a velocity field normal to the bounded direction is \textit{no-penetration}.

Flux boundary conditions

A flux boundary condition prescribes flux of a quantity normal to the boundary. For a tracer $c$ this corresponds to prescribing

where $\partial \Omega_b$ is an external boundary.

Gradient (Neumann) boundary condition

A gradient boundary condition prescribes the gradient of a field normal to the boundary. For a tracer $c$ this prescribes

Value (Dirichlet) boundary condition

A value boundary condition prescribes the value of a field on a boundary; for a tracer this prescribes

No penetration boundary condition

A no penetration boundary condition prescribes the velocity component normal to a boundary to be 0, so that