Temporal discretization: Asynchronous time stepping
FESOM2 uses asynchronous time stepping which means that velocity and scalar fields are shifted by a half time step. It is assumed that velocity is known at integer time steps and scalars at half-integer time steps. The time index is denoted by \(n\). The advantage of this simple arrangement is that the discrete velocity is centered in time for the time step of scalar quantities, and pressure gradient is centered for the velocity time step. All other terms in equations require some interpolations to warrant second-order accuracy in time, and routinely the second-order Adams-Bashforth estimate is used in FESOM2. Other variants of time stepping are explored, and will be added in future.
Thickness and tracers
We advance thickness and scalar quantities as
and
Here \(\tau\) is the time step and \(D_T\) stands for the terms related to diffusion. The time index on \(w\) is not written, for \(w\) is diagnosed by \({\bf u}\) and \(h\).
Note
Note that \(\mathbf{u}\) and \(h\) enter similarly in the equations for thickness and tracer. This warrants consistency: if \(T=\rm{const}\), the tracer equation reduces to the thickness equation.
Since the horizontal velocity is centered in time, these equations will be of the second order for the advective terms if \(h^*\) is centered (\(h^*=h^n\)) or if any other high-order estimate is used. The treatment of \(h^*\) depends on options. At present FESOM allows only the simplest choices when \(h\) vary only slightly with respect to unperturbed thickness and \(\partial_th\) is prescribed by the evolution of \(\eta\). We take \(h^*=h^{n-1/2}\) in this case because this choice can easily be made consistent with elevation, see further.
Although this formally reduces the time order to the first, the elevation is usually computed with the accuracy shifted to the first-order in large-scale ocean models, including this one.
Note
Other options, including those allowing \(h\) to follow isopycnal dynamic will be gradually added to FESOM. They will differ in the way how \(h^*\) is computed.
Elevation
We introduce
where \(H\) is the unperturbed ocean thickness. \(\overline h\) would be identical to the elevation \(\eta\) in the continuous world, but not in the discrete formulation here. The difference between these two quantities is that the elevation is defined at integer time levels. More importantly, it has to be computed so that the fast external mode is filtered. FESOM2 uses implicit method for that. The point is to make \(\overline h\) and \(\eta\) fully consistent.
In order to filter the external mode \(\eta\) is advanced implicitly in time. For \(h^*=h^{n-1/2}\) we write for the elevation
Here \(\alpha\) is the implicitness parameter (\(0.5\le\alpha\le1\)) in the continuity equation. Note that the velocities at different time steps are taken with their respective thicknesses in the same way as they appear in the thickness equation (14). The approach below is inspired by [CAHM04]. The equation for thicknesses can be vertically integrated giving, under the condition that the surface value of \(w^t\) vanishes,
Expressing the rhs in the formula for \(\eta\) (16) through the difference in surface displacements \(\overline{h}\) from the last formula we see that \(\eta\) and \(\overline{h}\) can be made consistent if we require
To eliminate the possibility for \(\eta\) and \(\overline{h}\) to diverge, we compute \(\eta^n\) from the last formula, then estimate \(\eta^{n+1}\) by solving dynamical equations (equation (16) requires \(\mathbf{u}^{n+1}\), so it is solved simultaneously with the momentum equation), and use it only to compute \({\bf u}^{n+1}\). On the new time step a ‘copy’ of \(\eta^{n+1}\) will be created from the respective fields \(\overline{h}\). We commonly select \(\alpha=1/2\), in this case \(\eta^n\) is just the interpolation between the two adjacent values of \(\overline{h}\). Note that (18) will be valid for any estimate \(h^*\) as far as it is used consistently in the product with the horizontal velocity.
Momentum equation
The momentum equation has to be solved together with the elevation equation (16), which is done with a predictor-corrector method. The method is largely an adaptation of pressure correction algorithm from computational fluid dynamics.
Assuming the forms (12) or (13) we write (using \(\partial_z\) for brevity instead of the top-bottom differences)
Here \(\theta\) is the implicitness parameter for the elevation, \({\bf R}^{n+1/2}_u\) includes all the terms except for vertical viscosity and the contribution from the elevation which are treated implicitly. To compute \({\bf R}^{n+1/2}_u\), we use the second-order Adams-Bashforth (AB) method for the terms related to the momentum advection and Coriolis acceleration. The AB estimate of quantity \(q\) is
Here \(\epsilon\) is a small parameter (\(\le0.1\)) needed to ensure stability in the case of advection operators. The contribution of pressure \(P\) does not need the AB interpolation (because it is centered). The horizontal viscosity is estimated on the level \(n\) because this term is commonly selected from numerical, not physical reasons.
We write the predictor equation
(19)\[{\bf u}^{*}-{\bf u}^{n}-\tau\partial_z\nu_v\partial_z({\bf u}^{*}-{\bf u}^n)=\tau({\bf R}^{n+1/2}_u+\partial_z\nu_v\partial_z{\bf u}^{n}-g\nabla\eta^n).\]The operator on the lhs connects three vertical levels, leading to three-diagonal linear problem for \(\Delta {\bf u}=\mathbf{u}_k^*-\mathbf{u}_k^n\) for each vertical column. Solving it we find the predicted velocity update \(\Delta {\bf u}\). (The vertical viscosity contribution on the rhs is added during the assembly of the operator on the lhs.)
The corrector step is written as
(20)\[{\bf u}^{n+1}-{\bf u}^{*}=-g\tau\theta\nabla(\eta^{n+1}-\eta^n).\]Expressing the new velocity from the last equation and substituting it into the equation for the elevation (16), we find
(21)\[\begin{split}\frac{1}{\tau}(\eta^{n+1}-\eta^n)-\alpha\theta g\tau\nabla\cdot(\overline{h}^{n+1/2}+H)\nabla(\eta^{n+1}-\eta^n)dz= \nonumber \\ -\alpha(\nabla\cdot\sum_k{h}_k^{n+1/2}({\bf u}^n+\Delta{\bf u})_k+W^{n+1/2})-(1-\alpha)(\nabla\cdot\sum{h}_k^{n-1/2}{\bf u}_k^n+W^{n-1/2}).\end{split}\]Here, the operator part depends on \(h^{n+1/2}\), which is the current value of thickness. The last term on the rhs is taken from the thickness computations on the previous time step.
The overall solution strategy is as follows.
Compute \(\eta^n\) from (18). Once it is known, compute \(\Delta {\bf u}\) from (19).
Update the matrix of the operator on the lhs of (21). Solve (21) for \(\eta^{n+1}-\eta^n\) using an iterative solver and estimate the new horizontal velocity from (20).
Compute \(\overline{h}^{n+3/2}\) from (17).
Determine layer thicknesses and \(w\) according to the options described below.
Advance the tracers. The implementation of implicit vertical diffusion will be detailed below.
Momentum equation in form (10)
Here an additional difficulty is the presence of \(h\) in the time derivative and on the rhs. The rule is that \(\mathbf{u}\) should appear with the same \(h^*\) as in the thickness or tracer equation. We used thus far the choice \(h^*\mathbf{u}^n=h^{n-1/2}\mathbf{u}^n\) in these equations, which implies that the time derivative will be
\(h^{n-1/2}\) will be used on the rhs with pressure gradient term, and the predictor equation will be written for \(h^{n-1/2}\Delta\mathbf{u}\). In this case \(h^{n-1/2}\) can be factored out of the lhs, which will make predictor solution similar. The corrector step will be modified to
It will lead to the replacement of \(h^{n+1/2}\) in the lhs of (21) by \(h^{n-1/2}\). We stress once again that the expressions in this section are for the particular choice of \(h^*\).
Current options for the vertical coordinate
The following options for the vertical coordinate are available at present:
Linear free surface: If we keep the layer thicknesses fixed, the time derivative in (14) drops out, and it becomes the standard equation to compute \(w\), starting from the bottom and continuing to the top,
\[w^t-w^b+\nabla\cdot({h\bf u})=0.\]
If this option is applied also to the first layer, the freshwater flux cannot be taken into account in the thickness equation. Its contribution to the salinity equation is added through a virtual salinity flux. In this option, \(w\) at the (fixed) ocean surface differs from zero, and so do the tracer fluxes. They do not necessarily integrate to zero over the ocean surface which is why tracer conservation is violated.
Full (nonlinear) free surface: We adjust the thickness of the upper layer, while the thicknesses of all other layers are kept fixed, \(\partial_th_k=0\) for \(k>1\). The thickness equations are used to compute \(w\) on levels \(k=2:K_v\) starting from the bottom. The change in the thickness of the first layer \(h^{n+3/2}_1-h^{n+1/2}_1\) is given by (17) written for the respective time interval. In this case there is no transport through the upper moving surface (the transport velocity \(w_1\) is identically zero). This option requires minimum adjustment with respect to the linear free surface. However, the matrix of the operator in (21) needs to be re-assembled on each time step.
We can distribute the total change in height \(\partial_t\overline h\) between several or all eligible layers. Due to our implementation, at each scalar horizontal location they can only be the layers that do not touch the bottom topography. If all eligible layers are involved we estimate
\[\partial_t h_k=(h_k^0/\tilde H)\partial_t\overline h,\]where \(h_k^0\) are the unperturbed layer thicknesses and \(\tilde H\) is their sum for all eligible layers. The thickness of the layers adjacent to the topography is kept fixed. The equation on thickness, written for each layer, is used to compute transport velocities \(w\) starting from zero bottom value. This variant gives the so-called \(z^*\)-coordinate.
Additional options will be gradually added. Layer thicknesses can vary in many ways provided that their tendencies sum to \(\partial_t\overline h\) over the layers. In particular, requiring that transport velocities \(w\) are zero, isopycnal layers can be introduced. The levels can move with high-pass vertical velocities, leading to the so called \(\tilde z\) coordinate, see [LM11], [PJR+15] or follow density gradients as in [HBB10]. The unperturbed layer thicknesses need not follow the geopotential surfaces and can be terrain following for example.
The ALE vertical coordinate is only a framework where many options are in principle possible. Additional measures may be required in each particular case, such as computations of pressure gradients with reduced errors. Updated transport algorithms may be needed (in the spirit of [LemarieKS+12] to minimize spurious numerical mixing in terrain-following layers. These generalizations are among the topics of ongoing work.
Implicit vertical diffusion
We return to the tracer equation (15). The vertical diffusion in this equation may present a CFL limitation and is treated implicitly.
Because of varying layer thicknesses, the implementation of implicit vertical diffusion needs slight adjustment compared to the common case of fixed layers. We write, considering time levels \(n-1/2\) and \(n+1/2\),
and split it into
and
Here \(R_T\) contains all advection terms and the terms due to the diffusion tensor except for the diagonal term with \(K_{33}\). The preliminary computation of \(T^*\) is necessary to guarantee that a uniform tracer distribution stays uniform (some significant digits will be lost otherwise).