Split-explicit external mode solver in the finite volume sea ice–ocean model FESOM2

A novel split-explicit (SE) external mode solver for the Finite volumE Sea ice–Ocean Model (FESOM2) and its sub-versions is presented. It is compared with the semi-implicit (SI) solver currently used in FESOM2. The split-explicit solver for the external mode utilizes a dissipative asynchronous (forward–backward) time-stepping scheme that does not require additional filtering during the barotropic sub-cycling. Its implementation with arbitrary Lagrangian–Eulerian vertical coordinates like Z star (z^*) and Z tilde ($\stackrel{\mathrm{̃}}{z}$) is explored. The comparisons are performed through multiple test cases involving idealized and realistic global simulations. The SE solver demonstrates lower phase errors and dissipation, while maintaining a simulated mean ocean state very similar to the SI solver. The SE solver is also shown to possess better run-time performance and parallel scalability across all workloads tested.

1 Introduction

All version 2 iterations of the Finite volumE Sea ice–Ocean Model (FESOM2; Danilov et al., 2017), the same as its predecessor FESOM1.4, rely on an implicit algorithm for solution of the external mode. Its computational algorithm maintains elementary options of the arbitrary Lagrangian–Eulerian (ALE) vertical coordinate, such as z^* or non-linear free surface, but needs modifications to incorporate more general options beginning from $\stackrel{\mathrm{̃}}{z}$, where information on horizontal divergence in scalar cells is used when making decisions about layer thicknesses on a new time level in the internal (baroclinic) mode. This work aims to present the modified time-stepping algorithm and its extension through a split-explicit option for solution of the external mode. Many modern ocean circulation models rely on the split-explicit method to solve for their external mode. The primary motivation behind such a choice is expectation of better parallel scalability in massively parallel applications. Indeed, as is well known, the need for global communications to calculate certain global dot products in most iterative solvers is a factor that potentially slows down the overall performance (see e.g. Huang et al., 2016; Koldunov et al., 2019). Although there are solutions minimizing the number of global communications per iteration (see e.g. Cools and Vanroose, 2017), as well as solutions where global communications are avoided (e.g. Huang et al., 2016), split-explicit methods are an obvious alternative. This method is followed by the Geophysical Fluid Dynamics Laboratory (GFDL) Global Ocean and Sea Ice Model OM4 whose ocean component uses version 6 of the Modular Ocean Model (MOM; Adcroft et al., 2019), the Nucleus for European Modelling of the Ocean (NEMO; Madec et al., 2019), the Regional Oceanic Modelling System (ROMS; Shchepetkin and McWilliams, 2005), and the Model for Prediction Across Scales – Ocean (MPAS-O; Ringler et al., 2013) to mention some widely used cases. A careful analysis in Shchepetkin and McWilliams (2005) discusses many details of the numerical implementation for a split-explicit external mode algorithm and proposes the AB3–AM4 (Adams–Bashforth and Adams–Moulton) method, which is at present followed by several models, i.e. ROMS (Shchepetkin and McWilliams, 2005), the Coastal and Regional Ocean COmmunity model (CROCO; Jullien et al., 2022), and FESOM-C (Androsov et al., 2019). However, a recent analysis in Demange et al. (2019) suggests a simpler choice of dissipative forward–backward time stepping. The built-in dissipation in this case allows one to avoid filtering the external mode solution (see Shchepetkin and McWilliams, 2005). The simplest commonly used filter requires that an external mode be stepped across two baroclinic time steps to ensure temporal centring. This doubles the computational cost of the external mode solution. In the forward–backward dissipative method by Demange et al. (2019) an external mode is stepped precisely across one baroclinic time step and not beyond. This method was ultimately found to be the best choice around which to build a split-explicit scheme for our purposes. Demange et al. (2019) also showed how dissipation can be added to the AB3–AM4 method of Shchepetkin and McWilliams (2005). We thus also explore dissipative AB3–AM4 for dissipation and phase errors. The rest of the paper is thus structured as follows. We begin with providing a breakdown for individual steps of split-explicit schemes as adopted in FESOM2 (Sect. 2). It is then followed by a comparison of the individual temporal discretizations that are characteristic of these time-stepping schemes (Sect. 3). We then perform various experiments comparing the new solver against the existing one using different test cases (Sects. 4 and 5). Finally, we summarize the results and argue for the new proposed scheme and its solver being a great choice for FESOM2 moving forward (Sect. 6).

2 Split-explicit asynchronous time stepping

This section provides a detailed description of the proposed asynchronous time-stepping scheme for FESOM2 that incorporates a split-explicit barotropic solver whose schematic can be found in Appendix A. FESOM in its standard version relies on a semi-implicit barotropic solver that already uses an asynchronous time stepping (Danilov et al., 2017). The asynchronous time stepping is a variant of forward–backward time stepping, which is formulated by considering scalar and horizontal velocities as being displaced by half a time step $\mathit{\tau }/\mathrm{2}$. The asynchronous time stepping is taken as the simplest option. Other time-stepping options for the baroclinic part, such as a third-order Runge–Kutta method, are under consideration for future versions.

2.1 Momentum equation

The standard set of equations under the Boussinesq and standard approximations is solved. The equations are taken in a layer-integrated form, and the placement of the variables on the mesh is explained in Danilov et al. (2017). The layer-integrated momentum equation in the flux form is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\partial }_{\mathrm{t}}{\mathit{U}}_k& +{\mathrm{\nabla }}_{\mathrm{h}}\cdot (\mathit{U}\mathit{u}{)}_k+(w^{\mathrm{t}}{\mathit{u}}^{\mathrm{t}}-w^b{\mathit{u}}^b{)}_k+f{\mathit{e}}_z\times {\mathit{U}}_k\\ & +h_k({\mathrm{\nabla }}_{\mathrm{h}}p+g\mathit{\rho }{\mathrm{\nabla }}_{\mathrm{h}}Z{)}_k/{\mathit{\rho }}_{\mathrm{0}}=(V_{\mathrm{h}}\mathit{U}+({\mathit{\nu }}_{\mathrm{v}}{\partial }_z\mathit{u}{)}^{\mathrm{t}}\\ & -({\mathit{\nu }}_{\mathrm{v}}{\partial }_z\mathit{u}{)}^b{)}_k,\end{array}\end{array}$$

with U_k=u_kh_k the horizontal transports, u the horizontal velocity, h_k the layer thickness, V_h the horizontal viscosity operator, ν_v the vertical viscosity coefficient, f the Coriolis parameter, e_z a unit vertical vector, and ${\mathrm{\nabla }}_{\mathrm{h}}=({\partial }_x,{\partial }_y)$ with respect to a constant model layer. Here, k is the layer index, starting from 1 ar the surface layer and increasing downward to the available number of levels with maximum value N_l. We ignore the momentum source due to the added water W at the surface. The term with the pressure gradient gρ∇_hZ_k accounts for the fact that layers deviate from geopotential surfaces. The quantity Z_k appearing in this term is the z coordinate of the midplane of the layer with the thickness h_k. The equation for elevation is written as

$$\begin{array}{}\text{(2)}& \partial \mathit{\eta }+{\mathrm{\nabla }}_{\mathrm{h}}\cdot \stackrel{\mathrm{‾}}{\mathit{U}}=W,\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathit{U}}={\sum }_k{\mathit{U}}_k$, and equations for layer thicknesses h_k and tracers will be presented in the following sections. Equations further in this section are for a particular layer k, and the index k will be suppressed. In the implementation described here, the discrete scalar state variables (elevation η, temperature T, salinity S, and layer thicknesses h) are defined at full time steps denoted by the upper index n, whereas the 3D velocities $\mathbf{v}=(\mathit{u},w)$ and horizontal transports U are defined at half-integer time steps ($n+\mathrm{1}/\mathrm{2},\mathrm{\dots }$). Since thicknesses and horizontal velocities are not synchronous, layer transports U are chosen as prognostic variables. Using them, we avoid the question of $h^{n+\mathrm{1}/\mathrm{2}}$ up to the moment of barotropic correction. Note that the flux form of the momentum advection is used in Eq. (1). Adjustments needed for other forms are straightforward and will not be discussed here.

First, we estimate the transport ${\mathit{U}}^{n+\mathrm{1}/\mathrm{2},*}$ assuming that $h^n,T^n,S^n$, ηⁿ, ${\mathit{u}}^{n-\mathrm{1}/\mathrm{2}}$, ${\mathit{U}}^{n-\mathrm{1}/\mathrm{2}}$, and $w^{n-\mathrm{1}/\mathrm{2}}$ are known.

$$\begin{array}{}\text{(3)}& {\mathit{U}}^{n+\mathrm{1}/\mathrm{2},*}-{\mathit{U}}^{n-\mathrm{1}/\mathrm{2}}=\mathit{\tau }\left[({\mathit{R}}_U^A+{\mathit{R}}_U^C{)}^n+({\mathit{R}}_U^{\mathrm{P}}{)}^n+({\mathit{R}}_U^{\mathrm{hV}}{)}^{n-\mathrm{1}/\mathrm{2}}\right]\end{array}$$

The ${\mathit{R}}_U^i$ terms with $i=A,C,P,\mathrm{hV}$ indicate advective, Coriolis, pressure gradient, and horizontal viscosity components estimated at time step n for $i=A,C,P$, and $n-\mathrm{1}/\mathrm{2}$ for the horizontal viscosity. The momentum advection term is

$${\mathit{R}}_U^A=-\mathrm{\nabla }\left(\mathit{u}\mathit{U}\right)-w\mathit{u}{|}_{\mathrm{b}}^{\mathrm{t}},$$

where ${|}_{\mathrm{b}}^{\mathrm{t}}$ implies that the difference between the top and bottom interfaces of layer k is taken. The Coriolis term is

$${\mathit{R}}_U^C=-f{\mathbf{e}}_z\times \mathit{U}.$$

Fields entering these advection and Coriolis terms are known at $n-\mathrm{1}/\mathrm{2}$, and the second- or third-order Adams–Bashforth method is used to get an estimate of ${\mathit{R}}_U^A$ and ${\mathit{R}}_U^C$ at n. For any quantity f, $f^{\mathrm{AB}}=(\mathrm{3}/\mathrm{2}+\mathit{\beta })f^n-(\mathrm{1}/\mathrm{2}+\mathrm{2}\mathit{\beta })f^{n-\mathrm{1}}+\mathit{\beta }{f}^{n-\mathrm{2}}$. For classical third-order interpolation (AB3), β is $\mathrm{5}/\mathrm{12}$, while β=0 gives the second-order result (AB2). The pressure gradient term is,

$${\mathit{R}}_U^{\mathrm{P}}=-h{\mathrm{\nabla }}_{z}p/{\mathit{\rho }}_{\mathrm{0}}=-h{\mathrm{\nabla }}_{\mathrm{h}}p/{\mathit{\rho }}_{\mathrm{0}}-hg{\mathrm{\nabla }}_{\mathrm{h}}Z,$$

where ∇_z means differencing at constant z. Since pressure and thicknesses are known at the time level n, no interpolation is needed. This is one of the advantages of asynchronous time stepping. Depending on how much layer thicknesses are perturbed, other algorithms than that written above can be applied to minimize pressure gradient errors. As a default, the approach by Shchepetkin and McWilliams (2003) is used in FESOM. The pressure contains contributions from η, density perturbations in the fluid column, and contributions from atmospheric and ice loading. The contribution from horizontal viscosity is of either the harmonic or biharmonic type. For simplicity, we write it here as

$${\mathit{R}}_U^{\mathrm{hV}}={\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(h{A}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}\mathit{u}\right).$$

The implicit contribution from vertical viscosity is added as

$$\begin{array}{}\text{(4)}& {\mathit{U}}^{n+\mathrm{1}/\mathrm{2},**}-{\mathit{U}}^{n+\mathrm{1}/\mathrm{2},*}=\mathit{\tau }{\mathit{R}}_U^{vV}=\mathit{\tau }\left(A_{\mathrm{v}}{\partial }_z{\mathit{u}}^{n+\mathrm{1}/\mathrm{2},**}\right){|}_{\mathrm{b}}^{\mathrm{t}}.\end{array}$$

The latter equation is rewritten for increments $\mathrm{\Delta }\mathit{u}={\mathit{u}}^{n+\mathrm{1}/\mathrm{2},**}-{\mathit{u}}^{n+\mathrm{1}/\mathrm{2},*}=({\mathit{U}}^{n+\mathrm{1}/\mathrm{2},**}-{\mathit{U}}^{n+\mathrm{1}/\mathrm{2},*})/h^{*}$ and solved for Δu. At this stage of the scheme, a reliable estimate for h^* at $n+\mathrm{1}/\mathrm{2}$ is not available, but, since A_v is a parameterization and this step is first order in time, we use $h^{*}=h^n$. When Δu is obtained, we update $\mathit{\tau }{\mathit{R}}_U^{vV}=\mathrm{\Delta }\mathit{u}{h}^n$ and ${\mathit{U}}^{n+\mathrm{1}/\mathrm{2},**}={\mathit{U}}^{n+\mathrm{1}/\mathrm{2},*}+\mathrm{\Delta }\mathit{u}{h}^n$. In preparation for the barotropic time step, the vertically integrated forcing from baroclinic dynamics is computed as

$$\begin{array}{}\text{(5)}& {\stackrel{\mathrm{‾}}{\mathit{R}}}^n=\sum _k\left[\right({\mathit{R}}_U^A{)}^n+({\stackrel{\mathrm{̃}}{\mathit{R}}}_U^{\mathrm{P}}{)}^n+({\mathit{R}}_U^{\mathrm{hV}}{)}^{n-\mathrm{1}/\mathrm{2}}+({\mathit{R}}_U^{vV}{)}^{n+\mathrm{1}/\mathrm{2}}{]}_k.\end{array}$$

Here, ${\stackrel{\mathrm{̃}}{R}}_U^{\mathrm{P}}$ represents the pressure gradient force excluding the contribution of η, as it will be accounted for explicitly in the barotropic equation. The Coriolis term is also omitted for the same reason. The vertically summed contribution from vertical viscosity is in reality the difference in surface stress and bottom stress. The bottom stress in FESOM is commonly computed as $C_{\mathrm{d}}\left|{\mathit{u}}_{\mathrm{b}}^{n-\mathrm{1}/\mathrm{2}}\right|{\mathit{u}}_{\mathrm{b}}^{n-\mathrm{1}/\mathrm{2}}$, where u_b is the bottom velocity.

2.2 Barotropic time stepping

Next is the barotropic step where η and $\stackrel{\mathrm{‾}}{\mathit{U}}={\sum }_k{\mathit{U}}_k$ are estimated by solving

$$\begin{array}{}\text{(6)}& {\partial }_{\mathrm{t}}\stackrel{\mathrm{‾}}{\mathit{U}}+f{\mathbf{e}}_z\times \stackrel{\mathrm{‾}}{\mathit{U}}+gH\mathrm{\nabla }\mathit{\eta }=\stackrel{\mathrm{‾}}{\mathit{R}},{\partial }_{\mathrm{t}}\mathit{\eta }+{\mathrm{\nabla }}_{\mathrm{h}}\stackrel{\mathrm{‾}}{\mathit{U}}+W=\mathrm{0}.\end{array}$$

Here, $H=H_{\mathrm{0}}+\mathit{\eta }$, W is the freshwater flux (positive out of ocean), and $\stackrel{\mathrm{‾}}{\mathit{R}}$ the is the forcing from the 3D part defined above. These equations are solved from time step n to n+1 as detailed below. Note that the baroclinic forcing term is taken at time level n. Centring it at $n+\mathrm{1}/\mathrm{2}$ would have involved many additional computations and is not implemented at present. This set of equations presents a minimum model. It is sufficient for basins with simple geometry. In realistic applications, it has been found that an additional viscous regularization term is needed to suppress oscillations in narrow straits with irregular coastlines. In such cases we add

$$\begin{array}{}\text{(7)}& {\stackrel{\mathrm{‾}}{\mathit{R}}}^{\mathrm{hV}}={\mathrm{\nabla }}_{\mathrm{h}}\cdot (H{\stackrel{\mathrm{‾}}{A}}_{\mathrm{h}}{\mathrm{\nabla }}_{\mathrm{h}}\stackrel{\mathrm{‾}}{\mathit{U}}/H)\end{array}$$

to the right-hand side of momentum Eq. (6) and subtract the initial value of this term from $\stackrel{\mathrm{‾}}{\mathit{R}}$ at each baroclinic time step. Here, ${\stackrel{\mathrm{‾}}{A}}_{\mathrm{h}}$ is the viscosity coefficient tuned experimentally to ensure stability in narrow, shallow regions. We express it as a combination of some background viscosity and a flow-dependent part that is proportional to the differences in barotropic velocity across cell edges. The subtraction mentioned serves to minimize the inconsistency created by adding the new term. Note that in coastal applications, one generally keeps bottom drag acting on the barotropic flow as well as barotropic momentum advection (Klingbeil et al., 2018). We treat them as slow processes here, but modifications might be needed for possible future applications.

As a default time stepping for the barotropic part, the forward–backward dissipative time stepping by Demange et al. (2019) is used. It is abbreviated as SE (for split-explicit, as it will be the default choice) in subsequent sections.

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+(m+\mathrm{1})/M}-{\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}=(\mathit{\tau }/M)[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}f{e}_z\\ & \times \left({\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+(m+\mathrm{1})/M}+{\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}\right)\\ & -g{H}^{n+m/M}{\mathrm{\nabla }}_{\mathrm{h}}{\mathit{\eta }}^{n+m/M}+{\stackrel{\mathrm{‾}}{\mathit{R}}}^n]\\ & {\mathit{\eta }}^{n+(m+\mathrm{1})/M}-{\mathit{\eta }}^{n+m/M}=(\mathit{\tau }/M)[-{\mathrm{\nabla }}_{\mathrm{h}}\\ & \cdot \left((\mathrm{1}+\mathit{\theta }){\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+(m+\mathrm{1})/M}-\mathit{\theta }{\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}\right)-W]\end{array}\end{array}$$

Here, M is the total number of barotropic substeps per baroclinic step τ, and θ controls dissipation. The value of θ=0.14 is mentioned by Demange et al. (2019) as being sufficient¹. This method is first-order accurate for θ≠0. We also used another version that is based on the AB3–AM4 (Adams–Bashforth – Adams–Moulton) approach of Shchepetkin and McWilliams (2005), with dissipative corrections as proposed in Demange et al. (2019) (abbreviated as SESM below). The specific versions of AB3 and AM4 used are

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{f}^{\mathrm{AB}\mathrm{3}}& =(\mathrm{3}/\mathrm{2}+\mathit{\beta })f^m-(\mathrm{1}/\mathrm{2}+\mathrm{2}\mathit{\beta })f^{m-\mathrm{1}}+\mathit{\beta }{f}^{m-\mathrm{2}}\\ f^{\mathrm{AM}\mathrm{4}}& =\mathit{\delta }{f}^{m+\mathrm{1}}+(\mathrm{1}-\mathit{\delta }-\mathit{\gamma }-\mathit{\zeta })f^m+\mathit{\gamma }{f}^{m-\mathrm{1}}+\mathit{\zeta }{f}^{m-\mathrm{2}},\end{array}\end{array}$$

with appropriate values of $\mathit{\beta },\mathit{\delta },\mathit{\gamma },\mathit{\zeta }$ discussed later. The time stepping takes the form²

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{\eta }}^{n+(m+\mathrm{1})/M}-{\mathit{\eta }}^{n+m/M}=(\mathit{\tau }/M)[-{\mathrm{\nabla }}_{\mathrm{h}}\cdot {\stackrel{\mathrm{‾}}{\mathit{U}}}^{\mathrm{AB}\mathrm{3}}-W]\\ & {\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+(m+\mathrm{1})/M}-{\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}=(\mathit{\tau }/M)[-f{e}_z\times {\stackrel{\mathrm{‾}}{\mathit{U}}}^{\mathrm{AB}\mathrm{3}}\\ & -g{H}^{\mathrm{AM}\mathrm{4}}{\mathrm{\nabla }}_{\mathrm{h}}{\mathit{\eta }}^{\mathrm{AM}\mathrm{4}}+{\stackrel{\mathrm{‾}}{\mathit{R}}}^n].\end{array}\end{array}$$

The modified AB3–AM4 scheme should also be moved to the first-order accuracy as concerns amplitudes but remains higher-order with respect to phase errors.

2.3 Reconciliation of barotropic and baroclinic mode

Note that the use of dissipative time stepping in Eq. (8) or Eq. (10) allows one to abandon the filtering of η and $\stackrel{\mathrm{‾}}{\mathit{U}}$ during the barotropic step that would be needed if non-dissipative forward–backward (θ=0) or the original AB3–AM4 schemes were applied instead (see Shchepetkin and McWilliams, 2005). The most elementary form of filtering involves integration to n+2 with subsequent averaging to n+1, which would double the computational expense for the barotropic solver. To be in agreement with the traditional notation (Shchepetkin and McWilliams, 2005), we write $\langle \mathit{\eta }{\rangle }^{n+\mathrm{1}}={\mathit{\eta }}^{n+m/M}$ and $\langle \stackrel{\mathrm{‾}}{\mathit{U}}{\rangle }^{n+\mathrm{1}}={\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}$for m=M (there would be a difference if filtering were needed). By summing the elevation equations over M substeps, one gets for the forward–backward dissipative case (Eq. 8)

$$\begin{array}{}\text{(11)}& \langle \mathit{\eta }{\rangle }^{n+\mathrm{1}}-\langle \mathit{\eta }{\rangle }^n=-\mathit{\tau }{\mathrm{\nabla }}_{\mathrm{h}}\cdot \langle \langle \stackrel{\mathrm{‾}}{\mathit{U}}\rangle {\rangle }^{n+\mathrm{1}/\mathrm{2}},\end{array}$$

where

$$\begin{array}{}\text{(12)}& \langle \langle \stackrel{\mathrm{‾}}{\mathit{U}}\rangle {\rangle }^{n+\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{M}}\sum _{m=\mathrm{1}}^M{\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}+{\displaystyle \frac{\mathit{\theta }}{M}}\left(\langle {\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+\mathrm{1}}\rangle -\langle {\stackrel{\mathrm{‾}}{\mathit{U}}}^n\rangle \right).\end{array}$$

While η is consistently initialized with 〈η〉ⁿ for m=0, there is no good answer for $\stackrel{\mathrm{‾}}{\mathit{U}}$. One can use the last available $\langle \stackrel{\mathrm{‾}}{\mathit{U}}{\rangle }^n$, but because 3D and barotropic velocities are integrated using different methods, this may lead to divergences with time unless some synchronization with 3D velocities is foreseen. We return to this topic below. On time level n+1, the total thickness becomes $H^{n+\mathrm{1}}=H^{\mathrm{0}}+\langle \mathit{\eta }{\rangle }^{n+\mathrm{1}}$. The horizontal transport is finalized by making the vertically integrated transport equal to the value obtained from the barotropic solution.

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}& ={\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2},**}\\ & -{\displaystyle \frac{h_k^{n+\mathrm{1}/\mathrm{2}}}{{\sum }_k{h}_k^{n+\mathrm{1}/\mathrm{2}}}}\left(\sum _k{\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2},**}-\langle \langle \stackrel{\mathrm{‾}}{\mathit{U}}\rangle {\rangle }^{n+\mathrm{1}/\mathrm{2}}\right)\end{array}\end{array}$$

2.4 Finalization of baroclinic mode

The estimate of the thickness at $n+\mathrm{1}/\mathrm{2}$ depends on the option of the ALE vertical coordinate and will be detailed below. In treating the scalar part, we are relying on the V-ALE (vertical ALE) approach in the terminology of Griffies et al. (2020). It is assumed that there is some external procedure to predict $h_k^{n+\mathrm{1}}=h_k^{\mathrm{target}}$ constrained by the condition ${\sum }_k{h}_k^{n+\mathrm{1}}=H^{\mathrm{0}}+\langle \mathit{\eta }{\rangle }^{n+\mathrm{1}}$. In the simplest case, this is the z^* vertical coordinate with $h_k^{n+\mathrm{1}}=h_k^{\mathrm{0}}(H/H^{\mathrm{0}})$. Here, as well as in other cases when the decision regarding h^target does not depend on layer horizontal divergences, $h^{n+\mathrm{1}/\mathrm{2}}$ in Eq. (13) is half the sum of the n and n+1 values. In more complicated cases, such as $\stackrel{\mathrm{̃}}{z}$ (Leclair and Madec, 2011; Petersen et al., 2015; Megann et al., 2022), the horizontal divergence in layers $\mathrm{\nabla }\cdot {\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}$ is needed to predict $h_k^{n+\mathrm{1}}$, and a reliable estimate of $h^{n+\mathrm{1}/\mathrm{2}}$ is not immediately available. Requiring that $h_k^{n+\mathrm{1}}$ is smooth and positive and also satisfies the barotropic constraint ${\sum }_k{h}_k^{n+\mathrm{1}}=H^{\mathrm{0}}+\langle \mathit{\eta }{\rangle }^{n+\mathrm{1}}$ could be a non-trivial task and may require a special procedure (see Hallberg and Adcroft, 2009 and Megann et al., 2022) that simultaneously adjusts ${\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}$ and hⁿ⁺¹. The description of the current implementation of $\stackrel{\mathrm{̃}}{z}$ in FESOM is presented in Appendix C. The potential presence of such complications is the reason why the decision regarding hⁿ⁺¹ is delayed to the end, and the discretization of momentum equation is performed in terms of U. Once the new thickness is determined, the thickness equation,

$$\begin{array}{}\text{(14)}& h_k^{n+\mathrm{1}}=h_k^n-\mathit{\tau }[{\mathrm{\nabla }}_{\mathrm{h}}\cdot \mathit{U}+w{|}_{\mathrm{b}}^{\mathrm{t}}{]}_k-\mathit{\tau }W{\mathit{\delta }}_{k\mathrm{1}},\end{array}$$

is used to estimate the diasurface velocity w. Tracers are then advanced, first taking into account advection and horizontal (isoneutral) diffusion before being trimmed by implicit vertical diffusion.

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& (h^{n+\mathrm{1}}{T}^{*}{)}_k=(hT{)}_k^n-\mathit{\tau }\left[\mathrm{\nabla }\right(UT)+(wT){|}_{\mathrm{b}}^{\mathrm{t}}{]}_k\\ & -\mathit{\tau }W{T}_W{\mathit{\delta }}_{k\mathrm{1}}+\mathit{\tau }\left(\mathrm{\nabla }\right(hK)\mathrm{\nabla }T{)}_k^n\\ & (h^{n+\mathrm{1}}{T}^{n+\mathrm{1}}{)}_k=(h^{n+\mathrm{1}}{T}^{*}{)}_k+\mathit{\tau }(K_v{\partial }_z{T}^{n+\mathrm{1}}{)}_k{|}_{\mathrm{b}}^{\mathrm{t}}\end{array}\end{array}$$

Here, T_W it the value of scalar T in a freshwater flux. In this procedure, if Tⁿ=const, the second equation will return this constant in T^*. The two equations above could have been combined into a single one. We treat them separately to avoid the loss of some significant digits (and to avoid ensuing errors in constancy preservation). Before solving the last equation in Eq. (15), it is rewritten for the increment $\mathrm{\Delta }T=T^{n+\mathrm{1}}-T^{*}$.

While ${\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}$, trimmed as given by Eq. (13), ensures by virtue of the first equation in Eq. (15) that ${\sum }_k{h}_k^{n+\mathrm{1}}=h^{\mathrm{0}}+\langle \mathit{\eta }{\rangle }^{n+\mathrm{1}}$ (as required for perfect volume conservation), its vertical sum ${\sum }_k{\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}=\langle \langle \stackrel{\mathrm{‾}}{\mathit{U}}\rangle {\rangle }^{n+\mathrm{1}/\mathrm{2}}$ deviates from the barotropic transport at time level $n+\mathrm{1}/\mathrm{2}$ (i.e for $m=M/\mathrm{2}$). We tried to compensate for this difference by saving ${\stackrel{\mathrm{‾}}{\mathit{U}}}^{n+m/M}$ for $m=M/\mathrm{2}$ in the barotropic step and re-trimming the 3D transports once tracers are advanced, but this has been found to be redundant in practice^{3 4}. The number of barotropic substeps M depends on the quality of meshes with varying resolutions. It can always be estimated based on mesh cell size and local depth. The presence of particularly small cells in deep water may limit the scheme globally even if the rest of the mesh is regular. Such limitations are absent in the present version of FESOM that is based on an implicit barotropic solver. Appendix B summarizes the changes needed to extend it (Danilov et al., 2017) to more general ALE options.

3 Temporal discretization for barotropic solver

This section provides a detailed numerical analysis of the new external mode solver when using SE or SESM time stepping vs. the current semi-implicit solver (Danilov et al., 2017) of FESOM2, which uses a first-order implicit time stepping in global simulations. Although parts of these implementations are already known, we repeat them here for clarity and comparison. A simple prototype system relevant for this analysis is

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\partial }_{\mathrm{t}}\stackrel{\mathrm{̃}}{u}& =-c_{\mathrm{p}}{\partial }_x\stackrel{\mathrm{̃}}{\mathit{\eta }}\\ {\partial }_{\mathrm{t}}\stackrel{\mathrm{̃}}{\mathit{\eta }}& =-c_{\mathrm{p}}{\partial }_x\stackrel{\mathrm{̃}}{u},\end{array}\end{array}$$

where $c_{\mathrm{p}}=(g{H}_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{2}}$ is the phase velocity, $\stackrel{\mathrm{̃}}{u}$ the dimensionless vertically averaged velocity ($\stackrel{\mathrm{‾}}{u}/c_{\mathrm{p}}$), and $\stackrel{\mathrm{̃}}{\mathit{\eta }}$ the dimensionless surface elevation ($\mathit{\eta }/H$). It will be assumed that $\stackrel{\mathrm{̃}}{\mathit{\eta }},\stackrel{\mathrm{̃}}{u}\sim e^{ikx}$, where k is the wave number.

3.1 Characteristic matrix forms

3.1.1 Semi-implicit (SI) method of Danilov et al. (2017)

We begin with the semi-implicit method of Danilov et al. (2017) used in the current model, FESOM2, which was adapted from FESOM 1.4 (Wang et al., 2014) and can be written as

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}& ={\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n-ic\left[\mathit{\alpha }{\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}+(\mathrm{1}-\mathit{\alpha }){\stackrel{\mathrm{̃}}{u}}^n\right]\\ {\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}& ={\stackrel{\mathrm{̃}}{u}}^n-ic\left[\mathit{\theta }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}+(\mathrm{1}-\mathit{\theta }){\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n\right].\end{array}\end{array}$$

Here, $\mathrm{1}/\mathrm{2}\le \mathit{\theta },\mathit{\alpha }\le \mathrm{1}$ are the control parameters, and c=c_pkτ is the Courant number. The characteristic matrix form for this scheme is

$$\begin{array}{}\text{(18)}& \left\{\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}\end{array}\right\}={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }\mathit{\theta }{c}^{\mathrm{2}}+\mathrm{1}}}\left[\begin{array}{cc}\mathit{\alpha }\mathit{\theta }{c}^{\mathrm{2}}-\mathit{\alpha }{c}^{\mathrm{2}}+\mathrm{1}& -ic\\ -ic& \mathit{\alpha }\mathit{\theta }{c}^{\mathrm{2}}-\mathit{\theta }{c}^{\mathrm{2}}+\mathrm{1}\end{array}\right]\left\{\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n\\ {\stackrel{\mathrm{̃}}{u}}^n\end{array}\right\}.\end{array}$$

3.1.2 Split-explicit (SESM) method of Shchepetkin and McWilliams (2005)

The explicit method of Shchepetkin and McWilliams (2005) based on an advanced forward–backward method combining the AB3 and AM4 steps can be expressed as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}& ={\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n-ic\left[\left(\mathrm{3}/\mathrm{2}+\mathit{\beta }\right){\stackrel{\mathrm{̃}}{u}}^n-\left(\mathrm{1}/\mathrm{2}+\mathrm{2}\mathit{\beta }\right){\stackrel{\mathrm{̃}}{u}}^{n-\mathrm{1}}+\mathit{\beta }{\stackrel{\mathrm{̃}}{u}}^{n-\mathrm{2}}\right]\\ {\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}& ={\stackrel{\mathrm{̃}}{u}}^n-ic\left[\mathit{\delta }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}+(\mathrm{1}-\mathit{\delta }-\mathit{\gamma }-\mathit{\zeta }){\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n+\mathit{\gamma }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n-\mathrm{1}}+\mathit{\zeta }{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n-\mathrm{2}}\right].\end{array}\end{array}$$

Here, too, $\mathit{\beta },\mathit{\delta },\mathit{\gamma },\mathit{\zeta }$ are the control parameters. Its characteristic matrix form is then

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& \left\{\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n-\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{u}}^n\\ {\stackrel{\mathrm{̃}}{u}}^{n-\mathrm{1}}\end{array}\right\}\\ & =\left[\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& -ic(\mathrm{3}/\mathrm{2}+\mathit{\beta })& ic(\mathrm{1}/\mathrm{2}+\mathrm{2}\mathit{\beta })& -ic\mathit{\beta }\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -ic(\mathrm{1}-\mathit{\gamma }-\mathit{\zeta })& -ic\mathit{\gamma }& -ic\mathit{\zeta }& \mathrm{1}-c^{\mathrm{2}}\mathit{\delta }(\mathrm{3}/\mathrm{2}+\mathit{\beta })& c^{\mathrm{2}}\mathit{\delta }(\mathrm{1}/\mathrm{2}+\mathrm{2}\mathit{\beta })& -c^{\mathrm{2}}\mathit{\delta }\mathit{\beta }\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\end{array}\right]\\ & \times \left\{\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n-\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n-\mathrm{2}}\\ {\stackrel{\mathrm{̃}}{u}}^n\\ {\stackrel{\mathrm{̃}}{u}}^{n-\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{u}}^{n-\mathrm{2}}\end{array}\right\}.\end{array}\end{array}$$

3.1.3 Split-explicit (SE) method of Demange et al. (2019)

Finally, the SE method by Demange et al. (2019) can be expressed as

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}& ={\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n-ic\left[(\mathrm{1}+\mathit{\theta }){\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}-\mathit{\theta }{\stackrel{\mathrm{̃}}{u}}^n\right]\\ {\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}& ={\stackrel{\mathrm{̃}}{u}}^n-ic{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n,\end{array}\end{array}$$

with θ being the control parameter. Its characteristic matrix form is

$$\begin{array}{}\text{(22)}& \begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^{n+\mathrm{1}}\\ {\stackrel{\mathrm{̃}}{u}}^{n+\mathrm{1}}\end{array}=\left[\begin{array}{cc}\mathrm{1}-c^{\mathrm{2}}(\mathrm{1}+\mathit{\theta })& -ic\\ -ic& \mathrm{1}\end{array}\right]\begin{array}{c}{\stackrel{\mathrm{̃}}{\mathit{\eta }}}^n\\ {\stackrel{\mathrm{̃}}{u}}^n\end{array}.\end{array}$$

3.2 Dissipation and phase analysis

Depending on the control parameters, the schemes above may lead to different dissipation and phase errors. Let the characteristic matrices of Eqs. (18), (20), and (22) be denoted by M_c. For I as an identity matrix of same rank as M_c and λ an eigenvalue of M_c, the characteristic polynomials for each scheme, given by det(M_c−λI), are

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& P(\mathit{\lambda }{)}_{i=\text{SE}}=-c^{\mathrm{2}}{\mathit{\theta }}_i+\mathrm{1}+(c^{\mathrm{2}}{\mathit{\theta }}_i+c^{\mathrm{2}}-\mathrm{2})\mathit{\lambda }+{\mathit{\lambda }}^{\mathrm{2}}\\ & P(\mathit{\lambda }{)}_{i=\text{SI}}={\displaystyle \frac{{\mathit{\alpha }}_i{c}^{\mathrm{2}}{\mathit{\theta }}_i-{\mathit{\alpha }}_i{c}^{\mathrm{2}}-c^{\mathrm{2}}{\mathit{\theta }}_i+c^{\mathrm{2}}+\mathrm{1}}{{\mathit{\alpha }}_i{c}^{\mathrm{2}}{\mathit{\theta }}_i+\mathrm{1}}}\\ & +\left({\displaystyle \frac{-\mathrm{2}{\mathit{\alpha }}_i{c}^{\mathrm{2}}{\mathit{\theta }}_i+{\mathit{\alpha }}_i{c}^{\mathrm{2}}+c^{\mathrm{2}}{\mathit{\theta }}_i-\mathrm{2}}{{\mathit{\alpha }}_i{c}^{\mathrm{2}}{\mathit{\theta }}_i+\mathrm{1}}}\right)\mathit{\lambda }+{\mathit{\lambda }}^{\mathrm{2}}\\ & P(\mathit{\lambda }{)}_{i=\text{SESM}}={\mathit{\beta }}_i{\mathit{\zeta }}_i{c}^{\mathrm{2}}+c^{\mathrm{2}}\left({\mathit{\beta }}_i({\mathit{\gamma }}_i-\mathrm{2}{\mathit{\zeta }}_i)-{\displaystyle \frac{{\mathit{\zeta }}_i}{\mathrm{2}}}\right)\mathit{\lambda }\\ & +-c^{\mathrm{2}}\left(({\mathit{\delta }}_i+\mathrm{3}{\mathit{\gamma }}_i-\mathrm{1}){\mathit{\beta }}_i+{\displaystyle \frac{{\mathit{\gamma }}_i}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{3}{\mathit{\zeta }}_i}{\mathrm{2}}}\right){\mathit{\lambda }}^{\mathrm{2}}\\ & +{\displaystyle \frac{c^{\mathrm{2}}}{\mathrm{2}}}\left({\mathit{\beta }}_i(\mathrm{6}{\mathit{\delta }}_i+\mathrm{6}{\mathit{\gamma }}_i+\mathrm{4}{\mathit{\zeta }}_i-\mathrm{4})+{\mathit{\delta }}_i+\mathrm{4}{\mathit{\gamma }}_i+{\mathit{\zeta }}_i-\mathrm{1}\right){\mathit{\lambda }}^{\mathrm{3}}+\\ & \left(\mathrm{1}+{\displaystyle \frac{c^{\mathrm{2}}}{\mathrm{2}}}\left({\mathit{\beta }}_i\right(-\mathrm{6}{\mathit{\delta }}_i-\mathrm{2}{\mathit{\gamma }}_i-\mathrm{2}{\mathit{\zeta }}_i+\mathrm{2})-\mathrm{4}{\mathit{\delta }}_i-\mathrm{3}{\mathit{\gamma }}_i-\mathrm{3}{\mathit{\zeta }}_i+\mathrm{3})\right){\mathit{\lambda }}^{\mathrm{4}}\\ & +\left(-\mathrm{2}+{\displaystyle \frac{c^{\mathrm{2}}{\mathit{\delta }}_i}{\mathrm{2}}}(\mathrm{2}{\mathit{\beta }}_i+\mathrm{3})\right){\mathit{\lambda }}^{\mathrm{5}}+{\mathit{\lambda }}^{\mathrm{6}}.\end{array}\end{array}$$

Here, the index i serves to distinguish between the SE, SI, and SESM schemes and their control parameters. The equations were obtained using the Maple symbolic solver. If λ is an eigenvalue, then P(λ)_i=0. Given that the physical eigenvalue should closely resemble the continuous solution λ=e^ic5, we can expand it for small c as

$$\begin{array}{}\text{(24)}& \mathit{\lambda }=\mathrm{1}+mc+n{c}^{\mathrm{2}}+q{c}^{\mathrm{3}}+O\left(c^{\mathrm{4}}\right).\end{array}$$

If the schemes are to be at least second-order dissipative with respect to c (see Demange et al., 2019), λ must also obey the relationship

$$\begin{array}{}\text{(25)}& \left|\mathit{\lambda }\right|=\mathrm{1}-\mathit{\chi }{c}^{\mathrm{2}}+O\left(c^{\mathrm{4}}\right),\end{array}$$

where χ is a parameter characterizing dissipation. Similarly, the phase ${\tan }^{-\mathrm{1}}\left(\mathrm{\Im }\right({\mathit{\lambda }}_{\mathrm{p}})/\mathfrak{R}({\mathit{\lambda }}_{\mathrm{p}}\left)\right)$ must also closely resemble the ideal phase c. The two conditions (24) and (25) then tie all control parameters together. From the requirement to remain formally second-order dissipative, one gets the conditions

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& \text{SE}\left(\mathit{\theta }\right)\text{, }\mathrm{2}\mathit{\chi }=\mathit{\theta }\\ & \text{SI}(\mathit{\theta },\mathit{\alpha })\text{, }\mathrm{2}\mathit{\chi }=\mathit{\theta }+\mathit{\alpha }-\mathrm{1}\\ & \text{SESM}(\mathit{\delta },\mathit{\gamma },\mathit{\zeta })\text{, }\mathrm{2}\mathit{\chi }=\mathit{\delta }-\mathit{\gamma }-\mathrm{2}\mathit{\zeta }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

with the requirement that m=i. Note that here $i=\sqrt{-\mathrm{1}}$ and that each equation in Eq. (26) admits its own set of parameters, i.e. θ_SI≠θ_SE. We also reject the possibility of $m=-i$ as it immediately gives the wrong phase.

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& \text{SE}\left(\mathit{\theta }\right)\text{, }n=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}+\mathit{\theta }),q=-{\displaystyle \frac{i}{\mathrm{8}}}(\mathrm{1}+\mathit{\theta }{)}^{\mathrm{2}}\\ & \text{SI}(\mathit{\theta },\mathit{\alpha })\text{, }n=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathit{\alpha }+\mathit{\theta }),q=-{\displaystyle \frac{i}{\mathrm{8}}}({\mathit{\alpha }}^{\mathrm{2}}+{\mathit{\theta }}^{\mathrm{2}}+\mathrm{6}\mathit{\alpha }\mathit{\theta })\\ & \text{SESM}(\mathit{\delta },\mathit{\gamma },\mathit{\zeta },\mathit{\beta })\text{, }n=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}+\mathit{\delta }-\mathit{\gamma }-\mathrm{2}\mathit{\zeta }\right),q\\ & =-{\displaystyle \frac{i}{\mathrm{8}}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[\mathit{\zeta }\right(\mathrm{16}\mathit{\gamma }-\mathrm{16}\mathit{\delta }+\mathrm{24})+\mathit{\gamma }(-\mathrm{8}\mathit{\delta }+\mathrm{4})-\mathrm{7}]\right.\\ & \left.+\mathrm{4}{\mathit{\zeta }}^{\mathrm{2}}+{\mathit{\gamma }}^{\mathrm{2}}+{\mathit{\delta }}^{\mathrm{2}}+\mathrm{3}\mathit{\delta }+\mathrm{4}\mathit{\beta }\right)\end{array}\end{array}$$

Note that here too as in Eq. (26) the parameters will be different for each scheme, i.e. θ_SI≠θ_SE. At this point, only the SE scheme is fully defined. The other schemes still have free parameters in need of optimization – α for SI and $\mathit{\beta },\mathit{\gamma },\mathit{\zeta }$ for the SESM scheme. As in Shchepetkin and McWilliams (2005), β can be set to 0.281105 for the largest stability limit. Given that dissipation is now the same (up to the second order), one can seek to optimize for phase errors. If third-order phase accuracy is desirable, then for the SI scheme, it is only possible if $\mathit{\alpha }=\mathit{\chi }+\mathrm{1}/\mathrm{2}\pm (\mathrm{1}/\mathrm{6})\sqrt{\mathrm{6}+\mathrm{72}{\mathit{\chi }}^{\mathrm{2}}}$. For the SESM scheme, this gives $\mathit{\gamma }=-{\mathit{\chi }}^{\mathrm{2}}-\mathrm{3}\mathit{\zeta }+\mathrm{1}/\mathrm{3}-\mathit{\beta }$. This still leaves ζ open for optimization. It can be obtained through further optimizing for either phase accuracy or stability limit. If optimizing for stability limit, the limit can be pushed much higher, as in Demange et al. (2019) if one relaxes the third-order phase accuracy constraint. The results of both optimizations are as follows,

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& \text{With}O\left(c^{\mathrm{3}}\right)\text{phase accuracy, }\mathit{\zeta }\approx -\mathrm{0.123}c+\mathrm{0.223}\\ & -\mathrm{0.169}\mathit{\beta }-\mathrm{0.169}{\mathit{\chi }}^{\mathrm{2}}\text{, }\mathit{\gamma }=\mathrm{1}/\mathrm{3}-\mathit{\beta }-\mathrm{3}\mathit{\zeta }-{\mathit{\chi }}^{\mathrm{2}}.\\ & \text{Without}O\left(c^{\mathrm{3}}\right)\text{phase accuracy, }\mathit{\zeta }\approx \mathrm{0.010}-\mathrm{0.135}\mathit{\chi }\text{, }\\ & \mathit{\gamma }=\mathrm{0.083}-\mathrm{0.514}\mathit{\chi }.\end{array}\end{array}$$

Here, β,δ retain their earlier description. To demonstrate the benefit in terms of phase accuracy for these split-explicit schemes, we analyse their net amplitude and phase errors per baroclinic time step, assuming that it consists of M=30 barotropic steps in Fig. 1, together with the errors in the SI scheme, vs. the baroclinic Courant–Friedrichs–Lewy (CFL) number c=c_pτk, where τ is the baroclinic time step. In practical cases, where the largest k is defined by the mesh resolution Δx (i.e. $k\to \mathit{\pi }/\mathrm{\Delta }x$), this baroclinic CFL number can become rather large. Because of the sub-cycling, the barotropic CFL is M times smaller and stays within the stability bounds of the explicit schemes. Figure 1 shows how all tested schemes in reality are able to maintain low dissipation even for high baroclinic CFLs. The choice is then made based on phase accuracy, which is very different between the implicit and explicit schemes. It is seen that both the SE and SESM schemes have orders-of-magnitude-lower phase errors compared to the SI scheme. For high c, the SI scheme has to be used with parameters α and θ to ensure strong damping of wavenumbers with large dispersive errors. Also, between the SESM and SE schemes, SESM seems to be the most accurate, even for same dissipation. To conclude the tests on phase accuracy, we report that irrespective of dissipation level, the split-explicit schemes will always by design provide orders-of-magnitude-better phase accuracy compared to the semi-implicit implementation, especially in the range of high CFL numbers, i.e. smaller wavelengths.

Figure 1 Comparison of amplitude and phase error for different schemes using the same dissipation χ=0.07 and maintaining at least second-order accuracy, i.e. $\mathit{\zeta }\approx \mathrm{0.010}-\mathrm{0.135}\mathit{\chi },\mathit{\gamma }=\mathrm{0.083}-\mathrm{0.514}\mathit{\chi }$, for the SESM scheme as per Eq. (28). Additionally, the SI scheme is also plotted using its recommended configuration $\mathit{\theta },\mathit{\alpha }=\mathrm{1}$, which successfully damps the high phase error solutions corresponding to large CFLs.

4 Numerical experiments

This section compares measurements from the new external mode solver to the existing one in FESOM2. The tests are done for both an idealized case and a realistic global setup. The idealized case is expected to highlight threshold performance of the new solver compared to the global case, where its impact will also be governed by mesh non-uniformity, the presence of external forcing, complicated boundaries, and bottom topography. The global case will, however, crucially assess the practicality of this new solver.

4.1 Idealized channel

In Sect. 3 (see Fig. 1), the primary characteristics of the new schemes were explored. In this idealized case, the solvers are tested for correct representation of the mean dynamics. We use the zonally re-entrant channel described in Soufflet et al. (2016). It is 2000 km long (north–south), 500 km wide (east–west), and 4 km deep. We test 10 km meshes of different types (triangular, quadrilateral) and unequally spaced vertical levels (40, 60). The baroclinic time step is τ=720 s, and surface gravity wave speed $c_{\mathrm{p}}=(g{H}_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{2}}=\mathrm{200}$ m s⁻¹ so that a mesh cell is crossed by waves in less than 50 s. We take M=30 for the split-explicit solver, which means $\mathit{\tau }/M=\mathrm{24}$ s. The initial density stratification due to temperature corresponds to a zonal jet. The zonally mean stratification and velocity are relaxed to their initial distributions. Some initial temperature perturbation leads to an onset of baroclinic instability, which is maintained through the relaxation of the zonal mean profiles. Simulations are run for 20 years, and the last 13 years are used to compute means. Figure 2 shows that mean depth profiles for this case are not affected by implementation of the new solver regardless of mesh configuration, i.e. a different mesh structure and number of vertical layers. This can be attributed to the predominantly baroclinic character of the flow, so the lower dissipation in the new solver is not necessarily seen. This is further supported by the observation that when the bottom drag is reduced, i.e. when the flow is made more barotropic, we see the separation in mean eddy kinetic energy between the two solvers grow, as shown in Fig. A2.

Figure 2 Comparison of area-averaged mean depth profiles for eddy kinetic energy (m² s⁻²), root-mean-square vertical velocity (m s⁻¹), and buoyancy flux (m² s⁻³), with the number of barotropic sub-cycles at M=30. Here, Δ and [] mean triangular and quadrilateral meshes, respectively. The meshes used have a fixed horizontal resolution of 10 km but varying vertical resolutions (40 or 60 layers).

4.2 Surface gravity wave

In this test, we further show how the SE solver for the barotropic mode being proposed is less dissipative and how the differences between it and the current SI solver of FESOM2 become more obvious when the barotropicity of the flow is dominant. For that, we use a simple surface gravity wave (SGW) setup where we simulate a channel (of the same geometry as in Sect. 4.1) with an initial elevation distribution that is meridionally Gaussian, i.e. $\ln (\mathit{\eta }/A)=-(y-y_{\mathrm{mid}}{)}^{\mathrm{2}}/{\mathit{\sigma }}^{\mathrm{2}}$, where A=3 m is the amplitude and σ=200 km is the half-bell width. The temperature is set at T=20 °C, the velocities are initialized to 0, and the simulation is run for 3 d with a baroclinic time step of τ=5 min. As seen in Fig. 3, which shows the final elevation for both cases after 3 d as well as time series of their available potential energy (APE), the APE for the SI solver asymptotes to 0 within the first 10 h, with the SI solver heavily damping all the SGWs, whereas the SE solver still maintains strong SGWs along with most of its APE. In contrast to the results for the baroclinic test case of Sect. 4.1, here we see a clear difference between the two solvers because of the flow being purely barotropic. As such, we see the gains in terms of the better energetic representation that the new SE solver promises.

Figure 3 Comparison of elevation (m) distribution snapshots for both solvers (after 3 d) and their total available potential energy density (m⁵ s⁻²) time series. The number of barotropic sub-cycles are M=30 and θ=0.14 for the SE solver, while for the SI solver, $\mathit{\alpha }=\mathit{\theta }=\mathrm{1}$ as per the recommended configuration. Here, triangular mesh with sides of 10 km and a baroclinic time step of τ=5 min is used.

4.3 Realistic global ocean

For this case, we now test a more complicated case of a global ocean–sea ice simulation similar to the one used by Scholz et al. (2022). We use the standard coarse mesh of FESOM2 with a minimum resolution of 25 km north of 25° N and a coarse resolution of around 1.5° in the interior of the ocean, with further moderate refinements in the equatorial belt and around Antarctica. The mesh configuration consists of 47 vertical levels with a minimum layer thickness of 10 m near the surface, up to 250 m near the abyssal depth. The baroclinic time step is τ=2700 s, and we take M=50 for the split-explicit solver, which means that $\mathit{\tau }/M=\mathrm{54}$ s. The simulations were forced with the JRA-55do v1.4.0 reanalysis data covering the period from 1958 to 2019. To show the differences in the simulations carried out with the SI and the SE barotropic solvers, we only show mean elevation, surface temperature, and kinetic energy over the last 20 years (1999–2019) of the simulation period. Due to the high similarity of SE and SESM results (as seen earlier in Fig. 2 for the idealized case), only SE results are shown in Fig. 4. The differences in sea surface elevation are found to be rather small. The pattern of difference in the sea surface temperature is most likely associated with the transient variability that is different in the two setups. The eddy kinetic energy increases everywhere outside the equatorial belt. This increase could be associated with the reduction in overall dissipation due to use of the SE barotropic solver and the observation that the barotropic kinetic energy contributes most to the overall kinetic energy budget at middle and high latitudes, as shown in Aiki et al. (2011).

Figure 4 Comparison of biases in sea surface heights (m), temperatures (°C), and eddy kinetic energies (m² s⁻²) using M=50 barotropic sub-cycles. The SE solver uses dissipation parameter θ=0.14, and depth-dependent fields are at 100 m.

In summary, no significant difference in terms of time-averaged measurements from the new SE external mode solver was observed. For both the idealized and the global test cases, the new SE external mode solver maintained mean dynamics close to those reported by the current SI solver.

5 Runtime performance and parallel scalability

This section compares the parallel scalability of the new external mode solver to the existing one of FESOM2. As in Sect. 4, we again utilize the two test cases – idealized and global – described earlier. Additionally, the two cases are also executed on different computer clusters providing for even better estimation of their general performance. Again, because of high similarities between SE and SESM parallel scalability in comparison to SI, the plots of SESM are omitted. In reality, SESM was found to be slightly less scalable than SE. The simulations are run for many model steps, and the mean total time per task the model spends for the barotropic solver is measured. For the idealized case, the simulations were performed using the Ollie HPC at the Alfred Wegener Institute equipped with Intel Xeon E5-2697 v4 (Broadwell) CPUs (308 nodes with 36 cores). To ensure a sufficient workload, we used a fine 2 km triangular mesh with 60 vertical layers on the same channel setup as in Soufflet et al. (2016). The mesh contains approximately 2.5×10⁵ vertices, so the setup is expected to scale almost linearly to about 10³ cores, according to our previous experience (Koldunov et al., 2019). The baroclinic time step has been reduced to τ=144 s, and M was left without changes. For this case, the simulations were run for 600 steps, i.e. 1 simulated day. As seen from Fig. 5 (left panel), the new external mode solver (SE) scales significantly better and is faster than the current SI solver of FESOM across all workloads. In reality, the relative speed of the SE versus SI solver will depend on M and on the efficiency of the preconditioner in SI and may change. Since the barotropic solver takes only a part of the total time step (10 %–20 %), the improved scalability of the SE solver contributes noticeably to the reduction in total computing time only after approximately 400 surface vertices per core, as shown in the left panel of Fig. 5. Its impact becomes significant only when parallelized beyond this limit. For the global case, the measurements were performed using the Albedo HPC at the Alfred Wegener Institute with 2xAMD Epyc-7702 CPUs (240 nodes with 128 cores). It uses the same setup and mesh from the global case in Sect. 4. The mesh contains approximately 1.27×10⁵ vertices. The baroclinic time step and M have been left unchanged (τ=2700 s and M=50, respectively). Here, simulations were run for 11 680 steps, i.e. 1 simulated model year. Similar to the findings from the idealized channel case, Fig. 5 (right panel) shows the new SE solver scaling similarly faster and further for the global case also. We again observe a perceivable difference across all workloads. Similar to the idealized case, these performance improvements only become significant for highly parallelized workflows, i.e. for fewer than approximately 400 vertices per core. In summary, the performance of the SE solver shows visible improvement in parallelization and computing time over the SI solver across all tested workloads. The behaviour remained the same over different test cases (idealized and global) and different computer resources (Ollie HPC, Albedo HPC). For less-parallel workloads, the benefits are marginal, but they become significant for highly parallelized workflows, i.e. in cases with fewer than approximately 400 vertices/core.

Figure 5 Scaling results for the idealized test case with a 2 km uniform mesh on the Ollie HPC cluster and the global test case with a 60–25 km unstructured mesh on the Albedo HPC cluster. The black line indicates linear scaling and the coloured lines give the mean computing time over the parallel partitions for the solver part of the code. Here, the wall-clock time measured corresponds to the model runtime per baroclinic step.

A preliminary implementation of the new split-explicit external mode solver within the sea ice model FESOM2 as proposed in this paper, including the test cases, can be found in the public repository at https://doi.org/10.5281/zenodo.10040944 (Banerjee et al., 2023).

6 Conclusions

The new split-explicit external mode solver proposed in this paper is more phase accurate, faster, and more scalable than the SI solver used in FESOM (Danilov et al., 2017). The dissipative asynchronous time-stepping scheme (SE) of Demange et al. (2019) is able to deliver phase accuracy orders of magnitude higher than the first-order SI scheme used before. It also provides comparable phase accuracy and dissipation to the dissipatively modified AB3–AM4 scheme (SESM) of Shchepetkin and McWilliams (2005). No filtering of fast dynamics is required due to the dissipative character of the SE solver. It is easier to implement compared to SESM, and leads to very similar results as SESM in practice.

The new SE solver is one part of the adjusted time stepping of FESOM that facilitates the use of the arbitrary Lagrangian–Eulerian vertical coordinate. As a demonstration, we extended z^* in FESOM2 to the $\stackrel{\mathrm{̃}}{z}$ vertical coordinate in the development version of FESOM2. The implementation of $\stackrel{\mathrm{̃}}{z}$ is outlined in Appendix C but still needs to be tested in realistic simulations. Across different test cases using different mesh geometries and computing resources, the new solver is shown to represent mean dynamics similarly to the existing solver with no significant difference. In the case of runtime performance and parallel scalability, it is shown to improve across all workloads. The improvements are shown to be especially significant for highly parallelized workloads. Kang et al. (2021) presents a semi-implicit solver for MPAS, showing that in contrast to the present work, it is more computationally efficient than their split-explicit solver. While a detailed answer to the question why the opposite conclusion was reached needs a separate study, here we can only mention that by following Demange et al. (2019), we perform fewer external time steps per baroclinic time step than in MPAS (Ringler et al., 2013).

We note that on unstructured meshes, a semi-implicit method can be more forgiving than a split-explicit one toward the size of mesh elements. A small element in deep water will hardly affect the solution of the semi-implicit solver but may require an increased number of barotropic substeps in a split-explicit method. This is why the semi-implicit option will be maintained in FESOM alongside the novel split-explicit option. It will however, be modified to allow more general ALE options as described in Appendix B. To conclude, this work suggests that the new split-explicit external mode solver is a good alternative to the existing solver of FESOM2.

Appendix A: Supplementary figures

Figure A1 Schematic diagram of the control flow for the split-explicit asynchronous time stepping proposed in this paper for FESOM2.

Figure A2 Comparison of area-averaged depth profiles for eddy kinetics (m² s⁻²) over a period of 10 years simulated with a baroclinic time step of τ=10 min on a 10 km triangular mesh with different linear bottom drags.

Appendix B: Adaptation of the semi-implicit scheme in FESOM2

The main difference from the split-explicit method is that the elevation η has to be defined at the same time levels as the horizontal velocity. The elevation is therefore detached from the thicknesses, which creates some conceptual difficulty. We consider quantities at $n-\mathrm{1}/\mathrm{2}$ and n to be known (to start from velocity).

• The predictor step is as follows

  $$\begin{array}{rl}{\mathit{U}}^{*}& ={\mathit{U}}^{n-\mathrm{1}/\mathrm{2}}+\mathit{\tau }({\mathit{R}}_U^A+{\mathit{R}}_U^C+{\stackrel{\mathrm{̃}}{\mathit{R}}}_U^{\mathrm{PGF}}{)}^n+\mathit{\tau }({\mathit{R}}_U^{\mathrm{hV}}{)}^{\mathrm{v}}\\ & -\mathit{\tau }g{h}^n{\mathrm{\nabla }}_{\mathrm{h}}{\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}\mathit{\tau }.\end{array}$$

  Here, tilde implies that the contribution from the elevation to the pressure gradient force (PGF) is omitted. It is taken into account explicitly (the last term). However, since ${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}$ is unknown, we take the value from the current time level $n-\mathrm{1}/\mathrm{2}$. Our intention is to get the Semi-Implicit form $\mathit{\theta }{\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}+(\mathrm{1}-\mathit{\theta }){\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}$, $\mathrm{1}/\mathrm{2}\le \mathit{\theta }\le \mathrm{1}$ in the end. The momentum advection and Coriolis terms are AB2 or AB3 interpolated to n. Implicit vertical viscosity is taken into account by solving

  $${\mathit{U}}^{**}={\mathit{U}}^{*}+\mathit{\tau }\left(A_{\mathrm{v}}{\partial }_z{\mathit{u}}^{**}\right){|}_{\mathrm{b}}^{\mathrm{t}}.$$

  It is solved similarly to the split-explicit asynchronous case.

• The corrector step is as follows

  $${\mathit{U}}^{n+\mathrm{1}/\mathrm{2}}={\mathit{U}}^{**}-\mathit{\tau }\mathit{\theta }g{h}^n{\mathrm{\nabla }}_{\mathrm{h}}({\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}-{\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}).$$

  This step is only written but is evaluated after ${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}$ is available.

• We write the elevation step as

  $${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}-{\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}=-\mathit{\tau }{\mathrm{\nabla }}_{\mathrm{h}}\cdot \sum _k(\mathit{\alpha }{\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}+(\mathrm{1}-\mathit{\alpha }\left){\mathit{U}}_k^{n-\mathrm{1}/\mathrm{2}}\right).$$

  Here, $\mathrm{1}/\mathrm{2}\le \mathit{\alpha }\le \mathrm{1}$, which is needed for stability. This equation has to be solved together with the corrector equation. We express ${\mathit{U}}^{n+\mathrm{1}/\mathrm{2}}$ from the corrector equation and insert the corrector step into the elevation equation to get

  $$\mathit{\delta }\mathit{\eta }=g\mathit{\theta }\mathit{\alpha }{\mathit{\tau }}^{\mathrm{2}}{\mathrm{\nabla }}_{\mathrm{h}}\cdot H^n{\mathrm{\nabla }}_{\mathrm{h}}\mathit{\delta }\mathit{\eta }-\mathit{\tau }{\mathrm{\nabla }}_{\mathrm{h}}\cdot \sum _k(\mathit{\alpha }{\mathit{U}}_k^{**}+(\mathrm{1}-\mathit{\alpha }\left){\mathit{U}}_k^{n-\mathrm{1}}\right).$$

  This equation is solved for $\mathit{\delta }\mathit{\eta }={\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}-{\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}$, giving ${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}={\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}+\mathit{\delta }\mathit{\eta }$.

• The corrector step is used to compute ${\mathit{U}}^{n+\mathrm{1}/\mathrm{2}}$.

• We write the thickness equation for the ALE step as

  $$h_k^{n+\mathrm{1}}-h_k^n=-\mathit{\tau }[{\mathrm{\nabla }}_{\mathrm{h}}\cdot {\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}+w{|}_{\mathrm{b}}^{\mathrm{t}}].$$

  These equations are summed vertically to give

  $$H^{n+\mathrm{1}}-H^n=-\mathit{\tau }{\mathrm{\nabla }}_{\mathrm{h}}\cdot \sum _k{\mathit{U}}_k.$$

  The quantity $H^{n+\mathrm{1}}-H^{\mathrm{0}}$ is the elevation at time step n+1. It is used to define hⁿ⁺¹ for the z^* vertical coordinate. The extension to $\stackrel{\mathrm{̃}}{z}$ follows similarly to the SE case. After hⁿ⁺¹ is defined, w is found from the thickness equation.

• The tracers are as follows

  $$\begin{array}{rl}(\mathrm{Th}{)}_k^{n+\mathrm{1}}& -(\mathrm{Th}{)}_k^n=-\mathit{\tau }\left[\mathrm{\nabla }\cdot \left({\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}{T}_k^{n+\mathrm{1}/\mathrm{2}}\right)+w{T}_k^{n+\mathrm{1}/\mathrm{2}}{|}_{\mathrm{b}}^{\mathrm{t}}\right]\\ & +\mathit{\tau }(\mathrm{\nabla }{h}^n\mathbf{K}\cdot {\mathrm{\nabla }}_{\mathrm{3}}{T}^n{)}_k+\mathit{\tau }(K_{\mathrm{v}}{\partial }_z{T}^{n+\mathrm{1}}){|}_{\mathrm{b}}^{\mathrm{t}}.\end{array}$$

• By virtue of the thickness equation above,

  $${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}-{\mathit{\eta }}^{n-\mathrm{1}/\mathrm{2}}=\sum _k\left[\mathit{\alpha }\right(h_k^{n+\mathrm{1}}-h_k^n)+(\mathrm{1}-\mathit{\alpha }\left)\right(h_k^n-h_k^{n-\mathrm{1}}\left)\right].$$

  We ignore the freshwater flux for simplicity, but it can be added. The solution is

  $${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}=\sum _k(\mathit{\alpha }{h}_k^{n+\mathrm{1}}+(\mathrm{1}-\mathit{\alpha }\left)h_k^n\right)-H^{\mathrm{0}}.$$

  If satisfied initially on cold start by formally taking ${\mathit{\eta }}^{-\mathrm{1}/\mathrm{2}}=\mathrm{0}$ and $h_k^{-\mathrm{1}}=h_k^{\mathrm{0}}$, this relationship will persist with time. However, to avoid accumulation of round-off errors, we reset ${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}$ to the right-hand side of the last expression after the computations of $h_k^{n+\mathrm{1}}$. This new ${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}$ will be used only in the next time step. The point here is that the original ${\mathit{\eta }}^{n+\mathrm{1}/\mathrm{2}}$ is computed by an iterative solver, whereby some significant digits are lost. The reset compensates for that. $\mathit{\alpha }=\mathrm{1}/\mathrm{2}$ provides centring in time.

Both split-explicit and semi-implicit asynchronous schemes are relatively straightforward to implement. The semi-implicit method with $\mathit{\theta }=\mathrm{1}/\mathrm{2}$ and $\mathit{\alpha }=\mathrm{1}/\mathrm{2}$ is non-dissipative, and dissipation is added by shifting θ toward 1. As explained above, even though the dissipation can be controlled well by offsetting $\mathit{\theta }=\mathrm{1}/\mathrm{2}$ only slightly, there are dispersive errors. Since the SI method is used with large Courant numbers for surface gravity waves, the contributions from such waves will come with large phase errors and should be damped. To keep the centring of η, we may take $\mathit{\alpha }=\mathrm{1}/\mathrm{2}$ and $\mathit{\theta }>\mathrm{1}/\mathrm{2}$. FESOM in most applications uses θ=1 and α=1, which implies more dissipation.

Appendix C: Implementation of $\stackrel{\mathrm{̃}}{z}$ in FESOM2

In the case of the $\stackrel{\mathrm{̃}}{z}$ vertical coordinate (Leclair and Madec, 2011), the horizontal divergence in a layer is split into fast and slow contributions. The fast one modifies layer thickness, and the slow one leads to diasurface w. Examples of practical implementation are provided by Petersen et al. (2015) and Megann et al. (2022). Our implementation presents a simplified version of both. The desired layer thickness is computed as

$$h_k^{n+\mathrm{1}}=h_k^{\mathrm{target}}=h_k^{*}+h_k^{\mathrm{hf}},$$

where $h_k^{*}$ corresponds to the z^* coordinate, and $h_k^{\mathrm{hf}}$ is the high-frequency component that augments z^* to $\stackrel{\mathrm{̃}}{z}$. They will be defined below. The bottom depth in FESOM is a cell-wise constant, whereas elevation and layer thicknesses are defined at vertices. For this reason, for a given vertex v, we modify thicknesses of $K^{\prime }=K^{\prime }\left(v\right)$ layers that do not touch topography (see Danilov et al., 2017). The total number of layers under vertex v will be denoted K=K(v). We take

$$(h_k^{*}{)}^{n+\mathrm{1}}=h_k^{\mathrm{0}}(\mathrm{1}+{\mathit{\eta }}^{n+\mathrm{1}}/H^{\prime }),H^{\prime }=\sum _{\mathrm{1}}^{K^{\prime }}{h}_k^{\mathrm{0}}.$$

An alternative definition would be to stretch the layers proportionally to their actual thickness, but Megann et al. (2022) warn that some drift in h^* may be present in such a case. Excluding the fixed layers, we split the divergence $D_k=\mathrm{\nabla }\cdot \left({\mathit{U}}_k\right)$ into a quasi-barotropic part that corresponds to $h_k^{*}$ and the remaining quasi-baroclinic part (“quasi” because we are limited to K^′ layers),

$$\begin{array}{}\text{(C1)}& D_k=D_k^{*}+D_k^{\prime },D_k^{*}=h_k^{\mathrm{0}}D/H^{\prime },\end{array}$$

where $D={\sum }_{k=\mathrm{1}}^K{D}_k$ is the vertically integrated divergence (note that all layers contribute to D). We will be interested in $D_k^{\prime }$, which is computed as the difference between D_k and $D_k^{*}$. We use the available thicknesses $h_k^n$ for $h_k^{n+\mathrm{1}/\mathrm{2}}$ in Eq. (13) to determine transports ${\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}$ featuring in D_k. After $h_k^{n+\mathrm{1}}$ is fully specified, we re-trim ${\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}$ using $h^{n+\mathrm{1}/\mathrm{2}}$ defined as a half sum of the thicknesses at full steps. Our treatment of the barotropic part is admittedly less accurate than in Petersen et al. (2015) and Megann et al. (2022), and some barotropic waves will contaminate $h_k^{\mathrm{hf}}$. However, because of the fixed bottom layers, we have already introduced uncertainty from the very beginning. Since ${\partial }_{\mathrm{t}}{h}_k^{*}=-D_k^{*}$, ${\sum }_{k=\mathrm{1}}^K{D}_k^{\prime }=\mathrm{0}$. The high-frequency thickness $h_k^{\mathrm{hf}}$ will be related to $D_k^{\prime }$ and should sum to zero vertically. D^′ is split into low- and high-frequency parts,

$$D_k^{\prime }=D_k^{\mathrm{lf}}+D_k^{\mathrm{hf}}.$$

The low-frequency part is nudged to $D_k^{\prime }$ as

$${\partial }_{\mathrm{t}}{D}_k^{\mathrm{lf}}=(\mathrm{2}\mathit{\pi }/{\mathit{\tau }}_{\mathrm{lf}})(D_k^{\prime }-D_k^{\mathrm{lf}}),$$

where τ_lf is the timescale (about 5 d in Petersen et al., 2015, but larger values can be of interest according to Megann et al., 2022). The fast-frequency part is obtained by subtracting the low-frequency part from $D_k^{\prime }$. The high-frequency contribution to thickness is

$$\begin{array}{}\text{(C2)}& {\partial }_{\mathrm{t}}{h}_k^{\mathrm{hf}}=-D_k^{\mathrm{hf}}-(\mathrm{2}\mathit{\pi }/{\mathit{\tau }}_{\mathrm{hf}})h^{\mathrm{hf}}+{\mathrm{\nabla }}_{\mathrm{h}}\left(K_{\mathrm{hf}}{\mathrm{\nabla }}_{\mathrm{h}}{h}_k^{\mathrm{hf}}\right).\end{array}$$

The second term on the right-hand side damps $h_k^{\mathrm{hf}}$ to zero over the timescale τ_hf (about 30 d). The last term will smooth the thickness, and the diffusivity K_hf is determined experimentally. If K_hf is vertically constant, ${\sum }_{k=\mathrm{1}}^K{h}_k^{\mathrm{hf}}=\mathrm{0}$ if it was initially so. A potential difficulty with Eq. (C2) is that $h_k^{\mathrm{hf}}$ is not bounded. A simple procedure is implemented at present. Equation (C2) is stepped implicitly with respect to the relaxation term, and diffusion is applied in a separate step. If $(h_k^{\mathrm{hf}}{)}^{n+\mathrm{1}}$ is outside the bounds for any k in the column at vertex v, τ_hf is adjusted for the entire column at this time step to ensure that $(h_k^{\mathrm{hf}}{)}^{n+\mathrm{1}}$ will be within the bounds, and computations of $(h_k^{\mathrm{hf}}{)}^{n+\mathrm{1}}$ are repeated. While this procedure is sufficient for the simple channel test case, it remains to be seen whether it will be sufficient in more realistic cases or whether the solutions reported by Megann et al. (2022) will be needed. The field h^hf is always damped stronger at locations close to topography to eliminate possible inconsistencies with h^hf=0 in cells touching bottom topography. After $(h_k^{\mathrm{hf}}{)}^{n+\mathrm{1}}$ is estimated, hⁿ⁺¹ is available; the transports ${\mathit{U}}_k^{n+\mathrm{1}/\mathrm{2}}$ can be re-trimmed, and diasurface velocities can be estimated from the thickness equations.