Tendency Terms

Tendencies for updating state variables are computed using functors, which are class objects that can be called like functions. These functors are used within Kokkos parallel loops to compute the full tendency arrays. The following tendency terms are currently implemented:

Name	Description
ThicknessFluxDivOnCell	divergence of layer thickness flux, defined at cell center
PotentialVortHAdvOnEdge	horizontal advection of potential vorticity, defined on edges
KEGradOnEdge	gradient of kinetic energy, defined on edges
SSHGradOnEdge	gradient of sea-surface height, multiplied by gravitational acceleration, defined on edges
VelocityDiffusionOnEdge	Laplacian horizontal mixing, defined on edges
VelocityHyperDiffOnEdge	biharmonic horizontal mixing, defined on edges
TracerHorzAdvOnCell	horizontal advection of thickness-weighted tracers
TracerHighOrderHorzAdvOnCell	second order horizontal advection of thickness-weighted tracers
TracerDiffOnCell	horizontal diffusion of thickness-weighted tracers
TracerHyperDiffOnCell	biharmonic horizontal mixing of thickness-weighted tracers
WindForcingOnEdge	forcing by wind stress, defined on edges
BottomDragOnEdge	bottom drag, defined on edges

Among the internal data stored by each functor is a bool which can enable or disable the contribution of that particular term to the tendency. These flags need to be provided in the configuration file:

Tendencies:
   ThicknessFluxTendencyEnable: true

Some functors have member variables that can be set by the user by the configuration file. The following are all the user-configurable parameters for the currently available tendency terms:

Term	Parameter	Description
ThicknessFluxDivOnCell	ThicknessFluxTendencyEnable	enable/disable term
PotentialVortHAdvOnEdge	PVTendencyEnable	enable/disable term
KEGradOnEdge	KETendencyEnable	enable/disable term
SSHGradONEdge	SSHTendencyEnable	enable/disable term
VelocityDiffusionOnEdge	VelDiffTendencyEnable	enable/disable term
	ViscDel2	horizontal viscosity
VelocityHyperDiffOnEdge	VelHyperDiffTendencyEnable	enable/disable term
	ViscDel4	horizontal biharmonic mixing coefficient for normal velocity
	DivFactor	scale factor for the divergence term
TracerHorzAdvOnCell	TracerHorzAdvTendencyEnable	enable/disable term
	HorzTracerFluxOrder	1 for standard linear advection
TracerHighOrderHorzAdvOnCell	TracerHorzAdvTendencyEnable	enable/disable term
	HorzTracerFluxOrder	2 for second order advection algorithm
TracerDiffOnCell	TracerDiffTendencyEnable	enable/disable term
	EddyDiff2	horizontal diffusion coefficient
TracerHyperDiffOnCell	TracerHyperDiffTendencyEnable	enable/disable term
	EddyDiff4	biharmonic horizontal mixing coeffienct for tracers
WindForcingOnEdge	WindForcingTendencyEnable	enable/disable term
BottomDragOnEdge	BottomDragTendencyEnable	enable/disable term
	BottomDragCoeff	bottom drag coefficient

Second Order Horizontal Advection Algorithm

The horizontal advection is done independently within each ocean layer so a two dimensional scheme is required to advect mixing ratios. The second order horizontal advection scheme is described in the paper

William C. Skamarock and Almut Gassmann, “Conservative Transport Schemes for Spherical Geodesic Grids: High-Order Flux Operators for ODE-Based Time Integration” Monthly Weather Review, Vol 139: Issue 9, pp 2962–2975, 2011. doi.org/10.1175/MWR-D-10-05056.1

and only a summary of a few equations from that paper are reproduced here.

The conservative form of the scalar transport equation within a fluid is given in equation (1) of Skamarock and Gassmann,

\[ \frac{\partial(\rho\psi)}{\partial t} = -\nabla\cdot\mathbf{V}\rho\psi \]

where \(\rho\) is the fluid density, \(\psi\) is the mixing ratio, and \(\mathbf{V}\) is the fluid velocity. This is simplified by assuming that \(\rho\) is constant and cancels out. The advection scheme for a finite volume method involves approximating this equation along an element edge at a certain time step and multiplying by the time step length to get an approximate scalar transport from element to element during the time step. The left hand side of this equation is a time derivative and for Omega this is handled by standard Runge—Kutta methods and not discussed here but can be found in the Skamarock and Gassmann paper. The right hand side of this equation is the spatial derivative that needs to be evaluated by the Omega horizontal advection scheme. Omega ocean uses unstructured Voronoi meshes and a finite volume scheme with \(\psi\) defined on cell centers. A first order approximation to the spatial derivative of the left hand termacross a single edge of an element would be the difference of \(\psi\) at element centroids divided by the distance from cell centroid to cell centroid. This is expressed in equation (5) of the paper. For edge \(i\) between elements \(i+1/2\) and \(i-1/2\),

\[ \frac{\partial(u\psi_i)}{\partial x} = \frac{1}{\Delta x}[F_{i+1/2}(u\psi) - F_{i-1/2}(u\psi)] + O(\Delta x^2). \]

Where \(u\) is \(\mathbf{V}\cdot\mathbf{n}\) where \(\mathbf{n}\) is the unit normal along the edge being evaluated and \(F\) is the evaluation of \(\psi\) for the elements on each side of the edge being evaluated. This is used when the TracerHorzAdvOnCell user option is active and HorzTracerFluxOrder is 1. To make the comparison to higher order methods explicit, this first order method is equivalent to defining \(\psi\) as a linear function, \(\psi = c_0 + c_x x + c_y y\), between the two elements sharing the edge \(i\) and taking the directed derivative along the edge.

For higher order flux calculations, higher order derivatives of \(\psi\) are needed and for the user option of TracerHorzAdvOnCell along with HorzTracerFluxOrder set to 2, a quadratic approximation of \(\psi\) is used,

\[ \psi = c_0 + c_x x + c_y y + c_{xx} x^2 + c_{xy} xy + c_{yy} y^2, \]

defined for each edge in the mesh. The six coefficients in this quadratic equation are calculated from a least squares fit of \(\psi\) in a neighborhood of elements around an edge. The least squares fit is defined in equation (12) of the paper,

\[ \mathbf{f = (P^T P)^{-1}P^T s = Bs} \]

where the vector \(\mathbf{f} = [c_0,c_x,c_y,c_{xx},c_{xy},c_{yy}]\) and \(\mathbf{P}\) is an \(m\times 6\) matrix with each row, \((1,x_i,y_i,x_i^2,x_iy_i,y_i^2)\), based on the cell center coordinates where there is one row for each of the \(m\) elements in the neighborhood. This results in a \(6\times m\) matrix \(\mathbf{B}\). The values of \(\mathbf{s}\) are the mixing ratios, \(\mathbf{s}=[\psi_0, \psi_1,...,\psi_m]\), of the elements in the neighborhood at the time step being evaluated. Note that \(\mathbf{B}\) depends on the mesh geometry only, so only needs to be computed once for each edge during initialization. Since only the second order terms are used to determine the second order derivatives needed for the higher order flux corrections to the lower order linear flux given above and the derivatives are taken along a line connecting cell centers, the amount of data that needs to be stored is minimized and \(\mathbf{B}\) reduces to a vector.

The second order derivatives of \(\psi\) are then combined as shown in equation (11) of Skamarock and Gassmann to get a second order horizontal transport algorithm.

The computation of these second order fluxes is complicated by having to do these computations on a sphere. Figure 2 from the Skamarock and Gassmann paper describs the slight modifications needed to compute the entries of the \(\mathbf{P}\) matrix for a spherical geometry where the distances need to be measured as arc lengths along the surface.

In practice, the convergence of this higher order algorithm on a sphere is slightly below second order. For the advection of a \(\psi\) in the shape of a cosine bell being advected around a sphere and compared with the analytic solution, the spatial convergence for grid refinement is about order 1.7 as shown in tracer-higher-order-convergence:

Tracer higer order convergence example of a cosine bell advected on a sphere showing an order 1.71 convergence rate