Vertical discretization: Layer thicknesses and layer equations
FESOM2 uses Arbitrary Lagrangian Eulerian (ALE) vertical coordinate. This implies that level surfaces are allowed to move. ALE vertical coordinate on its own is only the framework enabling moving level surfaces. The way how they are moving depends on one’s particular goal and may require additional algorithmic steps. Two limiting cases obviously include the case of fixed \(z\)-levels and the case when levels are isopycnal surfaces. At present only vertical coordinates that slightly deviate from \(z\)-surfaces are supported in FESOM, but many other options will follow.
The implementation of ALE vertical coordinate in FESOM2 basically follows [RPH+13]. An alternative approach, used by MOM6, see, e.g., [AH06] [AAB+19], is in the exploratory phase.
The essential step toward the ALE vertical coordinate lies in confining equations of section Main equations to model layers.
Introduce layer thicknesses \(h_k=h_k(x,y,t)\), where \(k=1:K\) is the layer index and \(K\) the total number of layers. They are functions of the horizontal coordinates and time in a general case. Each layer consists of prisms defined by the surface mesh but partly masked by bottom topography.
Layers communicate via the transport velocities \(w_{kv}\) through the top and bottom boundaries of the prisms. The transport velocities are the differences between the physical velocities in the direction normal to the layer interfaces and the velocities due to the motion of the interfaces. These velocities are defined at the interfaces (the yellow points in Fig. 3). For layer \(k\) the top interface has index \(k\) and the bottom one is \(k+1\). Note that \(w_{kv}\) coincides with the vertical velocity only if the level surfaces are flat.
All other quantities - horizontal velocities \({\bf u}\), temperature \(T\), salinity \(S\) and pressure \(p\) are defined at mid-layers. Their depths will be denoted as \(Z_k\), and the notation \(z_k\) is kept for the depths of mesh levels (the layer interfaces). They are functions of horizontal coordinates and time in a general case.
The equations of motion, continuity and tracer balance are integrated vertically over the layers. We will use \(T\) as a representative of an arbitrary tracer.
The continuity equation becomes the equation on layer thicknesses
and the tracer equation becomes
Here, \(W\) is the water flux leaving the ocean at the surface, it contributes to the first layer only (hence the delta-function); \(T_W\) is the property transported with the surface water flux and the indices \(t\) and \(b\) imply the top and the bottom of the layer.
The right hand side of (8) contains the 3 by 3 diffusivity tensor \({\bf K}\). We still use \(\nabla\) in (8) for the 3D divergence (the outer \(\nabla\)) for brevity, but assume the discrete form \(\nabla_h(...)+((...)^t-(...)^b)/h_k\), where \((...)\) are the placeholders for the horizontal and vertical components of 3D vector it acts on. A correct discretization of the diffusivity term is cumbersome and will be explained below.
Vertical sum of (7) over layers with account that \(w^t=0\) at the free surface and \(w_b=0\) at the bottom gives the ‘layered’ version of the elevation equation
(9)\[\partial_t\eta+\nabla_h\cdot\sum_kh_k{\bf u}_k+W=0.\]The layer-integrated momentum equation in the flux form is
(10)\[\begin{split}\partial_t(h{\bf u})+\nabla_h\cdot(h{\bf u u})+w^t{\bf u}^t-w^b{\bf u}^b+ f{\bf k}\times{\bf u}h +h(\nabla_h p+g\rho\nabla Z)/\rho_0= \nonumber \\ D_{uh}{\bf u}+(\nu_v\partial_z{\bf u})^t-(\nu_v\partial_z{\bf u})^b,\end{split}\]with \(D_{uh}{\bf u}\) the horizontal viscosity operator for the flux form (to be specified later), \(\nu_v\) the vertical viscosity coefficient, \(f\) the Coriolis parameter and \({\bf k}\) a unit vertical vector. We ignore the momentum source due to the added water \(W\) at the surface. Note that it could be more natural to formulate the solution procedure in terms of the horizontal layer transport velocities \({\bf U}=h{\bf u}\) in this case, but the present implementation in FESOM deals with separate \(h\) and \(\mathbf{u}\).
The pressure field is expressed as
(11)\[p=g\rho_0\eta+P, \quad P_{1}=p_a+g\rho_1h_1/2, \quad P_k=P_{k-1}+g(\rho_{k-1}h_{k-1}+ \rho_kh_k)/2.\]with \(p_a\) the atmospheric pressure, \(\rho\) the deviation of density from its reference value \(\rho_0\), and \(P\) is the variable hydrostatic pressure due to \(\rho\). The pressure gradient in continuous equations (1) has to be computed at constant \(z\). The model levels deviate from surfaces \(z=\rm{const}\). The term \(g\rho\nabla Z\), appearing together with the horizontal pressure gradient in (10) compensates for the deviation. The quantity \(Z\) appearing in this term is the \(z\)-coordinate of the midplane of the layer with the thickness \(h\).
Note
Although \(\nabla p+g\rho\nabla Z\) gives a formally correct estimate of pressure gradient at constant \(z\), the errors of discretization of the two terms in this expression become an issue if level surfaces deviate from \(z\)-surfaces. They are known as pressure gradient errors and require special care. FESOM will propose a selection of known algorithms, including the finite-volume algorithms of pressure gradient force that follows [EKM17] but is adapted to the triangular prisms of FESOM mesh.
Instead of using the flux form of momentum equation (10) representing momentum balance in the layer one can work without layer integration. Of particular interest is the vector-invariant form written as
(12)\[\begin{split}\partial_t{\bf u}+\frac{\omega+f}{h}{\bf k}\times{\bf u}h+((w\partial_z{\bf u})^t+(w\partial_z{\bf u})^b)/2 +\nabla (p/\rho_0+{\bf u}^2/2)+g\rho\nabla Z/\rho_0= \nonumber \\ D_u{\bf u}+((\nu_v\partial_z{\bf u})^t-(\nu_v\partial_z{\bf u})^b)/h.\end{split}\]Here, the identity
\[{\bf u}\cdot\nabla{\bf u}=\omega{\bf k}\times{\bf u}+\nabla({\bf u}^2/2),\quad \omega={\bf k}\cdot(\nabla\times{\bf u})\]was used.
The second term on the lhs of (12) includes division and multiplication with the layer thickness, and in doing so, it introduces the layer potential vorticity (PV), \(q=(\omega+f)/h\) and its transport \({\bf u}h\). The layer thickness formally drops out from the equation (12) which is still continuous in the horizontal direction. However, in the discrete case, the location of vorticity points (vertices) and velocity points is different. By keeping separate \(h\) the equation will then operate on the same horizontal transports as the thickness equations. This is the prerequisite for developing discretizations that conserve potential vorticity.
One more form is possible where the vector-invariant representation is not used
(13)\[\begin{split}\partial_t({\bf u})+\nabla\cdot({\bf u u})+(w^t{\bf u}^t-w^b{\bf u}^b)/h+ f{\bf k}\times{\bf u} +(\nabla p+g\rho\nabla Z)/\rho_0= \nonumber \\ D_{u}{\bf u}+(A_v\partial_z{\bf u})^t-(A_v\partial_z{\bf u})^b/h.\end{split}\]
The default version in FESOM2 is (10). Although the versions are derived from the same continuous equations, they are not equivalent in the discrete case.