Spatial discretization of equations (mostly in horizontal directions)
Finite-volume basics
To obtain the finite-volume discretization, the governing equations are integrated over the selected control volumes. The flux divergence terms are then, by virtue of Gauss theorem, transformed to the net fluxes leaving the control volumes. Basic discrete fields such as horizontal velocities and scalar tracers are considered to be mean over the respective finite volumes. We already described layered equations, i.e., the basic part of vertical discretization, where fluxes through top and bottom layer surfaces were appearing. In this section we discuss only the horizontal part. The discrete fields are understood as
and similarly for the temperature and other scalars,
Here \(A_c\) and \(A_{kv}\) are the horizontal areas of cells and scalar prisms. The scalar areas vary with depth, hence the index \(k\) in \(A_{kv}\) in the formula above (the index \(k\) will be suppressed in some cases). For layer \(k\), \(A_{kv}\) is the area of the prism \(kv\) including its top face. The area of bottom face is \(A_{(k+1)v}\) and may differ from that of the top one if the bottom is encountered. To be consistent in spherical geometry, we us
where \(\overline C(v)\) is the set of wet prisms containing \(v\) in layer \(k\).
Since the horizontal velocity is at centroids, its cell-mean value \({\bf u}_c\) can be identified with the value of the field \({\bf u}\) at the centroid of cell \(c\) with the second order of spatial accuracy. For scalar quantities a similar rule is valid only on uniform meshes, but even in this case it is violated in the vicinity of boundaries or topography. This has some implications for the accuracy of transport operators.
Horizontal operators
Scalar gradient
Scalar gradient takes vertex values of a scalar field and returns the gradient at the cell center:
where \({\bf n}_e\) is the outer normal to cell \(c\). Clearly \(l_e{\bf n}_e=-{\bf k}\times{\bf l}_e\) if \(c\) is the first (left) cell of \(c(e)\). This procedure introduces \({\bf G}_{cv}=(G^x_{cv},G^y_{cv})\) with the \(x\)- and \(y\) component matrices \(G^x_{cv}\) and \(G^y_{cv}\). They have three non-zero entries for each cell (triangle) which are stored in the order of triangle vertices, first for \(x\) and then for \(y\) component in the array gradient_sca(1:6,1:myDim\_elem2D) in the code. In contrast to FESOM1.4, where similar arrays are stored for each tetrahedron (and for 4 vertices and 3 directions), here only surface cells are involved.
Vector gradient
Vector gradient takes the values of velocity components at neighbor triangles of a given triangle and returns their derivatives at this triangle. They are computed through the least squares fit based on the velocities on neighboring cells sharing edges with cell \(c\). Their set is \(N(c)\). The derivatives \((\alpha_x, \alpha_y)\) of the velocity component \(u\) are found by minimizing
Here \({\bf r}_{cn}=(x_{cn}, y_{cn})\) is the vector connecting the center of \(c\) to that of its neighbor \(n\). The solution of the minimization problem can be represented as two matrices \(g_{cn}^x\) and \(g_{cn}^y\), acting on velocity differences \(u_n-u_c\) and returning the derivatives. Computations for \(v\)-component result in the same matrices. The explicit expressions for matrix entries are:
Here \(d=X^2Y^2-(XY)^2\), \(X^2=\sum_{n\in N(c)} x_{cn}^2\),
\(Y^2=\sum_{n(c)} y_{cn}^2\) and \(XY=\sum_{n\in N(c)} x_{cn}y_{cn}\). The matrix entries stored in the array gradient_vec(1:6,1:myDim_elem2D) in the order the neighbor elements are listed, first for \(x\) and then for \(y\) component. They are computed only once.
On the cells touching the lateral walls or bottom topography we use ghost cells (mirror reflections with respect to boundary edge). Their velocities are computed either as \({\bf u}_{n}=-{\bf u}_{c}\) or \({\bf u}_{n}={\bf u}_{c}-2({\bf u}_{c}\cdot{\bf n}_{nc}){\bf n}_{nc}\) for the no-slip or free-slip cases respectively. Here \(n\) is the index of the ghost cell, and \({\bf n}_{nc}\) is the vector of unit normal to the edge between cells \(c\) and \(n\). Note that filing ghost cells takes additional time, but allows using matrices \(g_{cn}^x\) and \(g_{cn}^y\) related to the surface cells only. Otherwise separate matrices will be needed for each layer. Note also that ghost cells are insufficient to implement the free-slip condition. In addition, the tangent component of viscous stress should be eliminated directly.
We stress that matrices \(g_{cn}^x\) and \(g_{cn}^y\) return derivatives of velocity components, and not the components of the tensor of velocity derivatives. The latter includes additional metric terms that need to be taken into account separately.
Flux divergence
Flux divergence takes fluxes defined at boundaries of scalar control volume (and in the simplest case just on cells) and returns their divergence on scalar control volumes:
where \({\bf n}_{ec}\) is the outer normal to control volume \(v\). Clearly, if \(v\) is the first vertex in the set \(v(e)\), \({\bf n}_{ec}d_{ec}=-{\bf k}\times{\bf d}_{ec}\) if \(c\) is the first in the set \(c(e)\) (signs are changed accordingly in other cases). While these rules may sound difficult to memorize, in practice computations are done in a cycle over edges, in which case signs are obvious.
In contrast to the scalar gradient operator, the operator of divergence depends on the layer (because of bottom topography), which is one of the reasons why it is not stored in advance. Besides, except for the simplest cases such as \(\mathbf{F}=\mathbf{u}\), the fluxes \({\bf F}\) involve estimates of the scalar quantity being transported. Computing these estimates requires a cycle over edges in any case, so there would be no economy even if the matrices of the divergence operator were introduced.
Velocity curl
Velocity curl takes velocities at cells and returns the relative vorticity at vertices using the circulation theorem:
where \({\bf t}_{ec}\) is the unit vector along \({\bf d}_{ec}\) oriented so as to make an anticlockwise turn around vertex \(v\). If \(v\) is the first in the set \(v(e)\) and \(c\) is the first in the set \(c(e)\), \({\bf t}_{ec}d_{ec}={\bf d}_{ec}\). This operator also depends on the layer and is not stored.
Mimetic properties
It can be verified that the operators introduced above are mimetic, which means that they satisfy the properties of their continuous analogs. For example, the scalar gradient and divergence are negative adjoint of each other in the sense of scalar products that provide energy norm, and the curl operator applied to the scalar gradient operator gives identically zero. The latter property implies that the curl of discrete pressure gradient is identically zero, which is a prerequisite of a PV conserving discretization.
Momentum advection
FESOM2.0 has three options for momentum advection. Two of them are different implementations of the flux form and the third one relies on the vector invariant form. In spherical geometry the flux form acquires an additional term \(M{\bf k}\times{\bf u}\) on the lhs, where \(M=u\tan\lambda/r_E\) is the metric frequency, with \(\lambda\) the latitude and \(r_E\) the Earth radius. All the options are based on the understanding that the cell-vertex discretization has an excessive number of velocity degrees of freedom on triangular meshes. The implementation of momentum advection must contain certain averaging in order to suppress the appearance of grid-scale noise.
Flux form with velocities averaged to vertices
The momentum flux is computed based on vertex velocities. We compute vertex velocities by averaging from cell to vertex locations
and use them to compute the divergence of horizontal momentum flux:
Here \({\bf n}_e\) is the external normal and \(l_e{\bf n}_e=-{\bf k}\times{\bf l}_e\) if \(c\) is the first one in the set \(c(e)\). Since the horizontal velocity appears as the product with the thickness, the expressions here can be rewritten in terms of transports \({\bf U}={\bf u}h^*\).
The fluxes through the top and bottom faces are computed with \(w_c=\sum_{v\in V(c)}w_v/3\) using either the second or fourth order centered, or high-order upwind algorithms.
Flux form relying on scalar control volumes
Instead of using vector control volumes, we assemble the flux divergence on the scalar control volumes and then average the result from the vertices to the cells. Here the same idea of averaging as in the previous case is applied to the momentum advection term instead of velocities. For the horizontal part,
with the same rule for the normals as in the computations of the divergence operator. The contributions from the top and bottom faces of scalar control volume are obtained by summing the contributions from the cells:
for the top surface, and similarly for the bottom one. The estimate of \({\bf u}^t\) can be either centered or upwind as above.
This option of momentum advection is special in the sense that the continuity is treated here in the same way as for the scalar quantities.
Vector-invariant form
The relative vorticity in the cell-vertex discretization is defined on vertices, and so should be the Coriolis parameter. We use the following representation
The representation with the thicknesses,
is reserved for future. The gradient of kinetic energy should be computed in the same way as the pressure gradient, which necessitates computations of \({\bf u}^2\) at vertices. This is done as
The vertical part follows (12),
for the top surface and similarly for the bottom. Note that the contributions from the curl of horizontal velocity, the gradient of kinetic energy and the vertical part involve the same stencil of horizontal velocities.
The three options above behave similarly in simple tests on triangular meshes, but their effect on flow-topography interactions or eddy dynamics still needs to be studied.
Horizontal viscosity operators
Because the cell placement of velocities, there are more velocity degrees of freedom than needed for vertex scalars. This leads to spurious grid-scale oscillations, and the task of horizontal viscosity operator is to eliminate them.
The derivatives of horizontal velocity can be estimated at cell locations and then averaged to edges of triangles enabling computation of the viscous stress tensor \(\sigma_{ij}=\nu_hs_{ij}\), \(s_{ij}=(\partial_iu_j+\partial_ju_i)/2\), where the indices \(i,j\) imply the horizontal directions, \(s_{ij}\) is the strain rate tensor and \(\nu_h\) is the harmonic horizontal viscosity coefficient. Its divergence will give the viscous force. This would be the standard way of introducing viscosity. It turns out that for cell velocities such a viscous force is insensitive to the difference in the nearest velocities (see [DK19]). To eliminate grid-scale fluctuations FESOM is bound to use other discretizations.
FESOM relies on a simplified stresses \(\sigma_{ij}=\nu_h\partial_iu_j\). As discussed by [Gri04], their divergence still ensures energy dissipation, but is nonzero for solid-body rotations if \(\nu_h\) is variable. In spite of this drawback, using the simplified form is much more convenient for numerical reasons: since the divergence of stresses reduces to fluxes over vertical faces of triangular prisms, only contraction of stresses with normal vector appears, i.e., \(\nu_hn_i\partial_iu_j\), which is the derivative in the direction of \({\bf n}\).
‘Canonical’ operators
‘Canonical’ harmonic viscosity operator is written as
The operator \(D_{uh}\) appearing in (10) will only differ by the absence of \(h_c\) on the lhs, and will not be written separately.
Let \(n\) be the cell sharing edge \(e\) with cell \(c\). We formally write for the edge normal vector \({\bf n}_e={\bf r}_{cn}/|{\bf r}_{cn}|+({\bf n}-{\bf r}_{cn}/|{\bf r}_{cn}|)\), where \({\bf r}_{cn}={\bf d}_{en}-{\bf d}_{ec}\) is the vector connecting the centroids of cells \(c\) and \(n\). Then
The first term on the rhs of (24) depends on the velocity difference across the edge, and this ensures that the viscous operator (23)-(24) will act to reduce the difference. The viscosity \(\nu_h\) has to be estimated at edges.
The operations are implemented in two cycles. The first one estimates velocity derivatives at cells. The second one is over edges, and edge values of velocity derivatives are obtained by averaging between the cells sharing the edge.
The ‘canonical’ biharmonic viscous operator is obtained by computing first the field
and then applying the operator \(D_{u}\) (23)-(24) to this field under the agreement that \(\nu_h=(\nu_{bh})^{1/2}\), where \(\nu_{bh}\) is positive biharmonic viscosity coefficient.
FESOM relies on an alternative version of biharmonic operator, where the viscosity coefficient is taken at \(c\) locations and the \(\mathbf{B}\) field becomes
In this case \(D_{u}\) with \(\nu_h=1\) is applied to \(\mathbf{B}\). The advantage of this form is that the combination \((\nu_{bh})_c/A_c\) in (25) has the dimension of harmonic viscosity \(\nu_h\) and thus can be replaced by \((\nu_h)_c\), avoiding the need to specify \((\nu_{bh})_c\). Even if viscosity is flow-dependent, one needs a single routine to compute \(\nu_h\).
Selection of viscosity coefficient
[FKM08] review common recipes for the viscosity coefficient and provide necessary references. Below we list the options available in FESOM. They can be used with both harmonic and biharmonic operators as explained above.
Simple viscosity A frequent practice in ocean modeling community is that \(\nu_h\) is selected as \(\nu_h=Vl\), where \(V\) is a velocity scale (commonly about 1 cm/s) and \(l\) is the cell size. In FESOM this rule is replaced by \((\nu_h)_c=VA_c^{1/2}\), and the edge values are obtained by averaging of values at neighboring cells. Note that FESOM choice automatically provides scaling.
Smagorinsky viscosity In this case \(\nu_h\) is taken as
\[(\nu_h)_c=C_{Smag}(A_c/\pi^2)(s_{ij}^2)^{1/2},\]where \(C_{Smag}\) is a dimensionless factor about one, and \(A_c/\pi^2\) is assumed to be an estimate of maximum wavenumber squared. In reality \(C_{Smag}\) is a tunable parameter available through FESOM namelists. The velocity derivatives are needed at \(c\) locations to compute \(s_{ij}^2\), see the section on velocity gradients.
Leith viscosity and its modification In this case
\[(\nu_h)_c=C_{Leith}(A_c/\pi^2)^{3/2}(|\nabla \omega|+C_{div}|\nabla\nabla\cdot\mathbf{u}|).\]Here \(C_{Leith}\) and \(C_{div}\) are dimensionless factors of order one. The term with \(C_{div}\) was added by [FKM08]. The entire construction is called the modified Leith viscosity. Here vorticity and divergence should be computed first at vertices. Then their scalar gradients are estimated, giving values at cells.
Note that the Smagorinsky and (modified) Leith viscosities rely on additional computations which take time. Viscosities are computed at \(c\) points and averaged to edges whenever needed.
Numerical viscosities (viscosity filters)
The canonical viscosity, especially its biharmonic version, proves to be costly. Its expensive part involving the computation of velocity derivatives at cells is only needed on deformed meshes, but even there it contributes little to penalizing differences between the nearest velocities. This leads to the idea of simplified numerical operators based only on the nearest neighbors. We refer to them as viscosity filters. They resemble canonic operators, but deviate from them on irregular meshes.
Quasi-harmonic viscosity filter
An analog of harmonic viscosity operator in FESOM is introduced as
Here \(l_{nc}\) is the length of the edge between cells \(n\) and \(c\), \(h_{nc}\) the layer thickness interpolated to the edge and \({\bf r}_{nc}\) is a vector connecting the cell centroids (we use \(nc\) to identify edges here). If mesh cells (triangles) are equilateral, \(({\bf u}_n-{\bf u}_c)/|{\bf r}_{nc}|\) is the velocity gradient in the direction of the normal to the edge between \(n\) and \(c\). The expression above is then the sum of viscous fluxes leaving cell \(c\), which gives a harmonic viscosity operator \(D_u\mathbf{u}=\nabla\nu\nabla{\bf u}\). This interpretation fails on general meshes, but the formula above for viscous operator provides the most efficient penalty of velocity difference between the nearest neighbors.
The appearance of \(A_ch_c\) on the left hand side is critical and is needed for conservation on general meshes (see below). Since \(l_{nc}/|{\bf r}_{nc}\) and the mean viscosity are related to the edge between \(n\) and \(c\), they can be incorporated in generalized viscosity \(\nu_{nc}\) associated to edges (hence \(\nu_{nc}=\nu_{cn}\)), which gives
Since \(\sum_cA_ch_c(D_u\mathbf{u})_c=0\) (the difference between \(n\) and \(c\) velocities appear with opposite signs in equations for \(n\) and \(c\)), the viscosity operator does not violate momentum conservation. Furthermore, taking \(\sum_c{\bf u}_cA_ch_c(V)_c\), we see that the contribution from \(n\) and \(c\) comes twice, one time as \(h_{nc}\nu_{nc}{\bf u}_c({\bf u}_n-{\bf u}_c)\) and the other time as \(\nu_{nc}h_{nc}{\bf u}_n({\bf u}_c-{\bf u}_n)\), which sums to \(-\nu_{nc}({\bf u}_n-{\bf u}_c)^2\). This proves that area mean kinetic energy dissipation is non-positive.
Quasi-biharmonic viscosity filter
It is only a bit more complicated with biharmonic implementation. First, we define \(b_{nc}=((\nu_{bh})_{nc}^b)^{1/2}\). We write
where
The momentum is conserved for the same reason as above. To see why the form above leads to kinetic energy dissipation we introduce matrix \(B\) with the dimension of the number of cells, such that its entries are \(B_{nc}=b_{nc}h_{nc}^{1/2}\) for those \(n\) and \(c\) that have common edges. It is a symmetric matrix. Then, we put the sum of row entries with the opposite sign at its diagonal. Obviously, in terms of this matrix \({\bf L}_c=B_{cn}{\bf u}_n\) (we assume summation over the repeating indices in matrix-vector multiplications). The energy dissipation is \(\sum_c{\bf u}_cA_ch_c(D_u^b\mathbf{u})_c\) can be written as \(-{\bf u}_mB_{mc}B_{cn}{\bf u}_n=-{\bf L}^T{\bf L}\le 0\). Note that area appears here only at the last step of the procedure, which makes the two steps different.
Similarly to the canonic biharmonic viscosity, the split of biharmonic viscosity between the two harmonic steps can be avoided. In this case we write
and
so that \({\bf L}_c=B_{cn}{\bf u}_n\) where \(B_{cn}=1\) if \(c\) and \(n\) are neighbors, and \(B_{cc}=-\sum_{n\ne c}B_{nc}\). We introduce a diagonal matrix \(N\) such that \(N_{cc}=(\nu_{bh})_ch_c\). The biharmonic operator can be written then as \(A_ch_c(D_u^b\mathbf{u})_c=-B_{cn}N_{nm}B_{mj}{\bf u}_j\), and making a dot product with \({\bf u}_c\) we will have a nonpositive energy dissipation. This alternative form is default in FESOM.
In all viscosity filter cases the main assembly is in the cycle(s) over edges.
Viscosities in viscosity filters
Although viscosity filters can be using the standard viscosities specified above, there are less expensive and more convenient options. For the quasi-harmonic viscosity filter (26) the edge harmonic viscosity \((\nu_h)_{nc}\) can be estimated in the same cycle where velocity differences are assembled as
where \(C_h\) is a small factor determined experimentally (about 1/20). Since it depends on velocity differences, it is an analog of the Smagorinsky viscosity. For the quasi-biharmonic viscosity filter (27),:eq:eq_viscQbh2, once \(\mathbf{L}_c\) is computed, the cell biharmonic viscosity can be estimated as
Since \(\mathbf{L}_c\) contains double differences (on uniform meshes), it is an analog of the Leith viscosity. The advantage of these simplified expressions is that they use the already computed terms, leading to a more economical implementation.
Transport of scalar quantities. Horizontal part
Horizontal and vertical fluxes are taken into account together, without operator splitting commonly used on structured meshes. However, the mesh is unstructured in the horizontal direction, the computation of the horizontal and vertical fluxes is different. They are therefore presented separately.
FESOM will provide the following horizontal advection schemes:
A third-fourth order scheme based on gradient estimate (GE34),
A third-fourth order scheme based on quadratic polynomial reconstruction (QR34)
A compact scheme (C34)
GE34 is available at present, and the two other will be made available in 2020 (their description will be added).
All of them can be run with flux corrected transport (FCT) limiter. The FCT operates on full three-dimensional fluxes.
GE34
Consider edge \(e\). Let \(V(e)=(v_1,v_2)\) and \(C(e)=(c_1, c_2)\). The advective flux of scalar quantity \(T\) through the face of scalar volume associated to this edge is
The quantity \(Q_e\) is the volume flux associated with edge \(e\) which leaves the control volume \(v_1\) and enters the control volume \(v_2\). To compute the mid-edge tracer estimate \(T_e\), for each edge \(e\) the indices of the cells up or down this edge in the edge direction \({\bf l}_e\) are stored. Two estimates
and
are computed. In these expressions, \((\nabla T)^c{\bf l}_e=T_{v2}-T_{v1}\) (\(c\) for centered), while \(u\) and \(d\) imply the gradients on up- and down-edge cells. They are computed and stored in a separate cycle by applying the scalar gradient operator.
The estimate
provides the standard third-order upwind method, and the estimate
provides the fourth-order centered method. The combination
takes the fourth-order part with the weight \(\gamma\) and the third order part, with \(1-\gamma\). \(\gamma=0.75\) will reducing the upwind dissipation by a factor of 4 compared to the third order method. The question whether this reduction is needed depends on applications.
The high order of the scheme above is only achieved on uniform meshes. Since \(T_e\) is computed through linear reconstruction, the second order is warranted on general meshes.
The scheme requires preliminary computation of scalar gradients on cells before the main cycle over edges. An extended halo exchange is needed to make these gradients available during flux assembly.
Edges touching the topography lack either \(u\) or \(d\) cells. In this case the simplest choice is either to use the central estimate or the estimate based on the mean vertex gradient \(A_v(\nabla T)_v=\sum_{c\in \overline C(v)}A_c(\nabla T)_c/3\). The associated logistics is expensive and increases the CPU cost of the scheme.
Vertical advection of scalars
The approach advocated by the ROMS community is to avoid dissipation in vertical advection, i. e. use high-order centered schemes. Although vertical grids are not uniform, this does not affect the convergence rate very much if level spacing is smoothly varying. For discussion, see Treguier et al. (1998). What is affected, of course, is the magnitude of residual errors. Many high-order algorithms essentially rely on mesh uniformity (some error cancellation). Nevertheless, they are applied on nonuniform meshes. We describe the selection of schemes available now or soon in FESOM. One-dimensional approach is used.
Centered second order
Let \(T\) be the scalar to be advected. A bar will be used to denote the interface value: \(\overline{T}_{kv}\) is the value reconstructed to the level \(k\) separating layer \(k\) above from layer \(k+1\) below. The ‘vertical’ flux \(F\) through the level \(k\) for the centered scheme is
Vertical GE34
We write
where \(G^c=(T_{(k-1)v}-T_{kv})/(Z_{(k-1)v}-Z_{kv})\) is the central gradient estimate and \(G^u=(T_{kv}-T_{(k+1)v})/(Z_{kv}-Z_{(k+1)v})\) for positive and \(G^u=(T_{(k-2)v}-T_{(k-1)v})/(Z_{(k-2)v}-Z_{(k-1)v})\) for negative \(w_{kv}\). Note that our estimates of gradients are based on values that are mean over control volume. So the estimates themselves are not very accurate. It is the combination (of central and upwind) values that is accurate.
Using
will give a fourth-order scheme on a uniform mesh if \(\gamma=1\). A blended third-fourth order scheme follows for \(0\le\gamma<1\).
Compact scheme (also the Parabolic Spline Method
We need scalar values at interfaces. An elegant way to find them is to use splines, requiring continuity of reconstruction and first derivatives at level locations. The result is
The boundary conditions are those of natural spline, i. e.,
This method requires three-diagonal solve, which takes the same time as two vertical loops. The name compact reflects the fact that the equation above involves stencil of minimum size. It becomes the PSM method if used with semi-Lagrangian time stepping, as in PPM.
The result is more accurate than PPM (see further). It is of the fourth order as PPM on uniform grid, but has a smaller residual term. Those who learned piecewise linear finite elements may see some analogies in the reconstruction procedure. ROMS uses this method for vertical advection of both tracers and momentum.
Piecewise Parabolic Method
To be written
FCT
The FCT limiter in FESOM2 uses the first-order upwind method as the low-order monotonic method and a combination of methods above as the high-order one. The low-order solution and the antidiffusive fluxes (the difference between the high-order and low-order fluxes) are assembled in one pass (in a cycle over edges for the horizontal part and over vertices for the vertical part). We experimented with separate pre-limiting of horizontal and vertical antidiffusive fluxes and found that commonly this leads to an increased dissipation, for the horizontal admissible bounds are in many cases too tight. For this reason, the computation of admissible bounds and limiting is three-dimensional. As a result, it will not necessarily fully eliminate non-monotonic behavior in the horizontal direction. The basic difference from the FCT algorithm used in FESOM1.4 is the construction of low-order solution. In FESOM1.4 the low-order solution is obtained by adding an artificial diffusion to the high-order right hand side. Using the FCT roughly doubles the cost of transport algorithm, but makes the code more stabe in practice.
Vertical velocity splitting
As demonstrated in [LemarieDM+15], the strongest practical Courant number limitation is imposed by vertical advection in isolated patches adjacent to the coast. The code numerical efficiency can be improved if measures are taken to stabilize it with respect to sporadic events with large vertical velocities. Unstructured meshes may even be more vulnerable to such events because mesh irregularity can easily provoke a noisy pattern in \(w\) just on its own. FESOM offers the approach proposed by [Shc15] according to which the vertical transport velocity is split into two contributions \(w=w_{ex}+w_{im}\) where the first one is determined by the maximum admissible Courant number, and the second one takes the rest. The advection with \(w_{ex}\) is done explicitly using schemes mentioned above. The advection with \(w_{im}\) is implicit. It uses the first-order upwind (backward Euler in time). This method leads to an operator that is diagonally dominant. The implicit advective terms are added to the implicit vertical mixing terms and the resulting three-diagonal system of equations is solved with the standard sweep algorithm. Because of this, only very small additional costs incur if this algorithm is used. Although the first order upwind scheme is dissipative, it is applied only in critical cases to excessively large velocities.
Operator splitting
FESOM2 does not use operator splitting at present and takes the horizontal and vertical fluxes in a single step. However, from the viewpoint of increasing admissible time steps it is worthwhile to provide the implementation of advection in which tracers are updated separately for horizontal and vertical contributions. As is well known, the sequence horizontal-vertical should alternate with vertical-horizontal in this case. This work is planned, and this section will be updated in due course.
GM and isoneutral operators
The eddy-induced transport
FESOM2 follows the algorithm proposed by [FNV10] to implement the Gent-McWilliams (GM) parameterization [GM90],:cite:Gent1995. FESOM1.4 operates with skewsion (see [Gri04] for mathematical detail). While working with skewsion is convenient in FESOM1.4 due to its variational formulation, it is less straightforward in FESOM2. Besides, the algorithm by [FNV10] provides an explicit expression for the eddy bolus velocity streamfunction.
The bolus velocity \({\bf v}^*=({\bf u}^*,w^*)\) is expressed in terms of eddy-induced streamfunction \(\boldsymbol{\Psi}\),
where \(\boldsymbol{\gamma}\) is a two-dimensional vector. In agreement with [FNV10], it is computed by solving
with boundary conditions \(\boldsymbol{\gamma}=0\) at the surface and ocean bottom. In this expression, \(c\) is the speed of the first baroclinic mode, \(\sigma\) the isoneutral density, \(\kappa\) the thickness diffusivity, \(N\) the Brunt–Väisälä frequency, and the index \(z\) means that the gradient is computed for fixed \(z\) (it differs from the gradient along layers, \(\nabla_z\sigma=\nabla\sigma-\partial_z\sigma\nabla Z\)). In terms of the vector \(\boldsymbol{\gamma}\) the components of eddy-induced velocity are computed as
It is easy to see that solving (29) plays a role of tapering, for the solution is a smooth function satisfying boundary conditions. The residual velocity \({\bf u}_r={\bf u}+{\bf u}^*\), \(w_r=w+w^*\) which is the sum of the eddy-induced velocity and the mean velocity \(({\bf u},w)\) is consistent with \(\overline h\) because the vertically integrated divergence of \({\bf u}^*\) is zero. The inclusion of eddy-induced velocity implies that the thickness and tracer equations are now written for the residual velocity \({\bf u}_r\).
Although the natural placement for \(\boldsymbol{\gamma}\) is at the cell centroids, it is moved to the mesh vertices in order to reduce the amount of computations. The vertical location is at full levels (layer interfaces). The horizontal bolus velocities are then computed at cell centroids as
The vertical bolus velocity \(w^*\) is then found together with \(w\) at the end of the ALE step and the full residual velocity is used to advect tracers.
We compute the speed \(c\) in the WKB approximation as
Among other factors, the magnitude of the thickness diffusivity \(\kappa\) depends on the resolution \(r\) and the local Rossby radius \(L_R=c/f\):
where \(f_{\kappa}\) is a cut-off function that tends to 0 if \(r/L_R<1\) and to 1 otherwise. The resolution is defined as a square root of the area of the scalar control volume. On general meshes it may exhibit substantial local variations, so smoothing over the neighbor vertices is done.
Isoneutral diffusion
Assuming that the slope of isopycnals is small, the diffusivity tensor can be written as
Here \(K_i\) and \(K_d\) are the isoneutral and diapycnal diffusivities, and \({\bf s}\) is the isoneutral slope vector. Its derivatives are computed along layers,
If layer interfaces deviate substantially from geopotential surfaces, for example, if layers follow the bottom topography, the slope vector can be substantially larger than typically found on \(z\)-coordinate meshes. Mixed derivatives in \(\nabla_3 h {\bf K}\nabla_3\) operator in this case can limit the time step [LemarieDSM12]. To maintain stability, the term \(h\partial_z(s^2K_i+K_d )\partial_z\) is treated implicitly, as suggested by [LemarieDSM12]. Appendix \(app:isoneutral\) shows the details of the numerical discretization of isoneutral diffusion.
Equation of state
FESOM still works with potential temperature. The conservative temperature and TEOS10 will be made available soon. The estimates of density by the equation of state are made columnwise. To facilitate these estimates, for each column the arrays are computed of quantities appearing with different powers in \(z\). Then they are combined to estimate the in-situ density and pressure as well as to compute the Brunt–Väisälä frequency in the same routine.