1. The Aronnax Model

1.1. The Physics

Aronnax can run in two different modes:

• n + 1/2 layer mode, in which the bottom layer is quiescent and infinitely deep
• n layer mode, in which the bottom topography is specified and all layers have finite thickness

1.1.1. n + 1/2 layer mode

The n + 1/2 layer mode is very cheap to run: Each time-step takes time linear in the number of grid points.

1.1.2. n layer mode

The n layer mode is more expensive, because the integration method is party implicit (to wit, solves a system of equations to ensure non-divergence of the barotropic flow). This tends to cost an extra factor of one grid side length in the run time, but admits a more realistic simulation, in that the effect of the ocean floor is included.

1.2. Governing Equations

The model solves the hydrostatic Boussinesq equations within a finite number of discrete isopycnal layers. At the layer interfaces there is a discrete jump in velocity and density, but not in pressure. The general forms of the equations are shown in Section 1.2.1 and Section 1.2.2, while a derivation of these equations for the 1.5 layer version of the model is shown in Section 1.2.3.

1.2.1. Continuity equation

(1.1)\[\frac{\partial h_{n}}{\partial t} = \mathbf{\nabla} \cdot \left(h_{n} \mathbf{v_{n}} \right).\]

in which \(h_{n}\) is the thickness of layer \(n\), and \(\mathbf{v_{n}}\) represents the vertically averaged horizontal velocity in layer \(n\).

1.2.2. Momentum equations

(1.2)\[\frac{D \mathbf{v_{n}}}{D t} + \mathbf{f} \times \mathbf{v_{n}} + g^{'}\mathbf{\nabla}h_{n} = \mathbf{F_{n}},\]

in which \(g^{'}\) is the reduced gravity given by \({g(\rho_{2} - \rho_{1})}/{\rho_{1}}\). The reduced gravity is dynamically equivalent to gravity, but is scaled to take into account the density difference between the two layers.

This can be rewritten in terms of the Bernoulli Potential to give,

(1.3)\[\frac{\partial\mathbf{v_{n}}}{\partial t} - (f+\zeta_{n}) \times v_{n} + \nabla \Pi_{n} + = \kappa \nabla^{2}v_{n} + \frac{\mathbf{F_{n}}}{\rho_{0}}\]

where \(\Pi_{n}\) is the Bernoulli potential, \(\left(\mathbf{v_{n}}\cdot\mathbf{v_{n}}\right)/2 + p/\rho_{0}\), and \(p\) is the hydrostatic pressure. In this form the non-linearity from the material derivative has been moved into the Bernoulli Potential and the vorticity term.

The model can be used in either reduced gravity mode, with a quiescent abyss, or in n-layer mode with bathymetry. In the n-layer case the model can either be run with a rigid lid, or with a free surface. In n-layer simulations the following equation is also solved

(1.4)\[\frac{\partial \eta}{\partial t} + \mathbf{\nabla} \cdot (H \mathbf{V}) = 0,\]

where \(\eta\) is the free surface height, \(H\) is the depth from the free-surface to the bathymetry, and \(V\) is the vertically averaged flow, the barotropic flow. With a rigid lid \(\eta\) represents the pressure field required to keep the vertically integrated horizontal flow divergence free - the result is not carried from one timestep to the next.

1.2.3. Reduced gravity equations

The most idealised version of Aronnax is a 1.5 layer or ‘reduced gravity’ model. Reduced gravity models are often referred to as 1.5 layer models because of the limited dynamics in the lowest layer. The bottom layer is assumed to have no motion, which is mathematically equivalent to an infinite thickness. Since there is no motion in the lowest layer, the model is unable to support vertically homogeneous, or barotropic, motions. This means that the gravest mode supported by the model is the first baroclinic mode. This version of Aronnax also features a rigid lid to remove surface gravity waves from the model solutions. This allows for longer time steps in the numerical solver. Reduced gravity models are often used to model the dynamics of the upper ocean, with the interface between the layers taken as an approximation for the oceanic thermocline. A schematic of a 1.5 layer model is shown in Figure 1.1.

The equations for a model with one active layer above a quiescient abyss will be derived here. Extending the derivation to a model with \(n\) active layers requires substantially more algebra but gives little additional insight.

Figure 1.1 Schematic of a reduced gravity model in a rectangular domain. The fluid within the infinitely deep, quiescent abyss is assumed to be at rest.

If we take the conservation of mass and recast it as conservation of mass in a layer of thickness \(h\) we get

(1.5)\[\frac{D\left(\rho h\right)}{Dt} + \rho h \mathbf{\nabla} \cdot \mathbf{V} = 0,\]

in which \(\mathbf{V}\) represents the vertically averaged horizontal velocity in the layer, and \(\rho h\) is the areal density. By expanding the material derivative, combining terms and assuming that \(\rho\) is constant within an isopycnal layer, (1.5) can be rewritten as

(1.6)\[\frac{\partial h}{\partial t} = \mathbf{\nabla} \cdot \left(h\mathbf{V}\right).\]

With a rigid lid the surface is maintained at a constant height of zero, but the pressure is unconstrained. If we let the pressure at the rigid lid be given by \(P(x,y,t)\), then the pressure in the upper layer at depth \(z\) is

(1.7)\[P_{1}\left(x,y,z,t\right) = P\left(x,y,t\right) - \rho_{1}gz,\]

where \(\rho_{1}\) is the density of the upper layer, \(z\) is the vertical coordinate which becomes more negative with depth. A defining feature of reduced gravity models is the absence of motion in the lowest layer. This means that the horizontal pressure gradients in layer 2 are identically zero, which we can use to solve for the interface displacement. The pressure in layer 2 is given by

(1.8)\[P_{2}(x,y,z,t) = P_{1}(x,y,h,t) - \rho_{2}g(z+h) = P + \rho_{1}gh + \rho_{2}g(z+h),\]

where \(h\) is the thickness of the upper layer. Since a central assumption of the reduced gravity framework is that the horizontal gradients of \(P_{2}\) are zero we can now solve for the horizontal pressure gradient in the upper layer. Taking the gradient of equation (1.8) and setting the left-hand side to zero gives

(1.9)\[0 = \mathbf{\nabla} P + \mathbf{\nabla}\rho_{1}gh + \mathbf{\nabla}\rho_{2}g(-h),\]

which can be rearranged to give

(1.10)\[\mathbf{\nabla}P = g(\rho_{2} - \rho_{1}) \mathbf{\nabla}h,\]

which relates the horizontal pressure gradients in the upper layer to displacements of the interface. The momentum equation for the upper layer is therefore

(1.11)\[\frac{D\mathbf{V}}{Dt} + \mathbf{f} \times \mathbf{V} + g^{'}\mathbf{\nabla}h = \mathbf{F},\]

in which \(g^{'}\) is the reduced gravity given by \({g(\rho_{2} - \rho_{1})}/{\rho_{1}}\). The reduced gravity is dynamically equivalent to gravity, but is scaled to take into account the density difference between the two layers.

1.3. Discretisation

Aronnax is discretised on an Arakawa C grid.

Figure 1.2 A single grid cell from an Arakawa C grid.

1.4. Numerical algorithm

The model solves for two horizontal velocity components and layer thickness in an arbitrary number of layers. The model supports two sets of physics: either a reduced gravity configuration in which the horizontal pressure gradient is set to zero in a quiescent abyss below the lowest active layer; or an n-layer configuration in which bathymetry must be specified.

Aronnax is discretised on an Arakawa C-grid, with the velocity and thickness variables in different locations on the grid cell.

The choice of quiescent abyss or n-layer physics is made by a runtime parameter in the input file. The numerical algorithm for calculating the values at the next time level, \(n+1\), is as follows:

• The Bernoulli Potential is calculated using values from time-level \(n\)
  • The function used depends on whether the model is running in reduced gravity mode or n-layer mode
• The relative vorticity is calculated using values from time-level \(n\)
• The layer thickness tendencies are calculated using the velocities and layer thicknesses from time-level \(n\)
• the velocity tendencies are calculated using values from time-level \(n\)
• the layer thicknesses and velocities are stepped forward in time to \(n+1\) using a third-order Adams-Bashforth algorithm and the stored time derivatives from the previous two timesteps. N.B. for the n-layer version these velocities are not strictly at time \(n+1\), let’s call it time level \(n+*\).
• For the n-layer version:
  • The no-normal flow boundary condition is applied (perhaps unnecessary?)
  • The barotropic velocity required to keep the vertically integrated flow non-divergent in the horizontal is calculated and added to the baroclinic velocities calculated previously. To do this:
    • the barotropic velocities are calculated from the velocities at time-level \(n+*\).
    • the divergence of these velocities is used to solve for the free surface elevation at time-level \(n+1\) that makes the barotropic flow non-divergent
      • This is the step that requires the linear system solve, since we solve the equation implicitly to sidestep the issue of requiring a very short \(\delta t\).
    • the barotropic correction is applied to the velocity fields
    • consistency between the sum of the layer thicknesses and the depth of the ocean is forced by applying a uniform inflation/deflation to the layers. (the model currently prints a warning if the discrepancy is larger than a configurable threshold, which defaults to 1%)
• The no normal flow and tangential (no-slip or free-slip) boundary conditions are applied
• The layer thicnkesses are forced to be larger than a configurable minimum. This is for numerical stability and is probably only necessary for the layer receiving the wind forcing. This is discussed in ticket #26
• the arrays are shuffled to prepare for the next timestep.

N.B. To get the Adams-Bashforth method going, two time steps are initially performed using Runge-Kutta 4th order time stepping.