FESOM-C v.2: coastal dynamics on hybrid unstructured meshes

Androsov, Alexey; Fofonova, Vera; Kuznetsov, Ivan; Danilov, Sergey; Rakowsky, Natalja; Harig, Sven; Brix, Holger; Wiltshire, Karen Helen

doi:https://doi.org/10.5194/gmd-12-1009-2019

Articles | Volume 12, issue 3

https://doi.org/10.5194/gmd-12-1009-2019

Articles | Volume 12, issue 3

Model description paper

21 Mar 2019

Model description paper |

| 21 Mar 2019

FESOM-C v.2: coastal dynamics on hybrid unstructured meshes

Alexey Androsov, Vera Fofonova, Ivan Kuznetsov, Sergey Danilov, Natalja Rakowsky, Sven Harig, Holger Brix, and Karen Helen Wiltshire

Abstract

We describe FESOM-C, the coastal branch of the Finite-volumE Sea ice – Ocean Model (FESOM2), which shares with FESOM2 many numerical aspects, in particular its finite-volume cell-vertex discretization. Its dynamical core differs in the implementation of time stepping, the use of a terrain-following vertical coordinate, and the formulation for hybrid meshes composed of triangles and quads. The first two distinctions were critical for coding FESOM-C as an independent branch. The hybrid mesh capability improves numerical efficiency, since quadrilateral cells have fewer edges than triangular cells. They do not suffer from spurious inertial modes of the triangular cell-vertex discretization and need less dissipation. The hybrid mesh capability allows one to use quasi-quadrilateral unstructured meshes, with triangular cells included only to join quadrilateral patches of different resolution or instead of strongly deformed quadrilateral cells. The description of the model numerical part is complemented by test cases illustrating the model performance.

Download & links

Article (PDF, 9639 KB)

Download & links

How to cite

How to cite.

Dates

Received: 23 Apr 2018 – Discussion started: 25 Jul 2018 – Revised: 08 Dec 2018 – Accepted: 15 Feb 2019 – Published: 21 Mar 2019

1 Introduction

Many practical problems in oceanography require regional focus on coastal dynamics. Although global ocean circulation models formulated on unstructured meshes may in principle provide local refinement, such models are as a rule based on assumptions that are not necessarily valid in coastal areas. The limitations on dynamics coming from the need to resolve thin layers, maintain stability for sea surface elevations comparable to water layer thickness, or simulate the processes of wetting and drying make the numerical approaches traditionally used in coastal models different from those used in large-scale models. For this reason, combining coastal and large-scale functionality in a single unstructured-mesh model, although possible, would still imply a combination of different algorithms and physical parameterizations. Furthermore, on unstructured meshes, the numerical stability of open boundaries, needed in regional configurations, sometimes requires masking certain terms in motion equations close to open boundaries. This would be an unnecessary complication for a large-scale unstructured-mesh model that is as a rule global.

The main goal of the development described in this paper was to design a tool, dubbed FESOM-C, that is close to FESOM2 (Danilov et al., 2017) in its basic principles, but can be used as a coastal model. Its routines handling the mesh infrastructure are derived from FESOM2. However, the time stepping, vertical discretization, and particular algorithms, detailed below, are different. FESOM-C relies on a terrain-following vertical coordinate (vs. the arbitrary Lagrangian–Eulerian (ALE) vertical coordinate of FESOM2), but takes a step further with respect to the mesh structure. It is designed to work on hybrid meshes composed of triangles and quads. Some decisions, such as the lack of the ALE at the present stage, are only motivated by the desire to keep the code as simple as possible through the initial phase of its development and maintenance. The code is based on the cell-vertex finite-volume discretization, the same as FESOM2 (Danilov et al., 2017) and FVCOM (Chen et al., 2003). It places scalar quantities at mesh vertices and the horizontal velocities at cell centroids.

Our special focus is on using hybrid meshes. In essence, the capability of hybrid meshes is built on the finite-volume method. Indeed, computations of fluxes are commonly implemented as cycles over edges, and the edge-based infrastructure is immune to the polygonal type of mesh cells. However, because of staggering, it is still convenient to keep some computations on cells, which then depend on the cell type. Furthermore, high-order transport algorithms might also be sensitive to the cell geometry. We limit the allowed polygons to triangles and quads. Although there is no principal limitation on the polygon type, triangles and quads are versatile enough in practice for the cell-vertex discretization. Our motivation for using quads is twofold (Danilov and Androsov, 2015). First, quadrilateral meshes have 1.5 times fewer edges than triangular meshes, which speeds up computations because cycles over edges become shorter. The second reason is the intrinsic problem of the triangular cell-vertex discretization – the presence of spurious inertial modes (see, e.g., Le Roux, 2012) and decoupling between the nearest horizontal velocities. Although both can be controlled by lateral viscosity, the control leads to higher viscous dissipation over the triangular portions of the mesh. The hybrid meshes can be designed so that triangular cells are included only to optimally match the resolution or are even absent altogether. For example, FESOM-C can be run on curvilinear meshes combining smooth changes in the shape of quadrilateral cells with smoothly approximated coastlines. One can also think of meshes for which triangular patches are only used to provide transitions between quadrilateral parts of different resolution, implementing an effective nesting approach.

Many unstructured-mesh coastal ocean models were proposed recently (e.g., Casulli and Walters, 2000; Chen et al., 2003; Fringer et al., 2006; Zhang and Baptista, 2008; Zhang et al., 2016). It will take some time for FESOM-C to catch up in terms of functionality. The decision on the development of FESOM-C was largely motivated by the desire to fit in the existing modeling infrastructure (mesh design, analysis tools, input–output organization), and not by any deficiency in existing models. The real workload was substantially reduced through the use or modification of the existing FESOM2 routines.

We formulate the main equations and their discretization in the three following sections. Section 5 presents results of test simulations, followed by a discussion and conclusions.

2 Model formulation

2.1 The governing equations

We solve standard primitive equations in the Boussinesq, hydrostatic, and traditional approximations (see, e.g., Marshall et al., 1997). The solution is sought in the domain $\hat{Q} = Q \times [0, t_{f}]$ , where t_f is the time interval. The boundary ∂Q of domain Q is formed by the free water surface, the bottom boundaries, and lateral boundaries composed of the solid part ∂Q₁ and the open boundary ∂Q₂, $Q = {x, y, z; x, y \in Ω, - h (x, y) \leq z < ζ (x, y, t)}$ , $0 \leq t \leq t_{f}$ . Here ζ is the surface elevation and h the bottom topography. We seek the vector of unknown $q = (u, w, ζ, T, S)$ , where $u = (u, v)$ is the horizontal velocity, w the vertical velocity, T the potential temperature, and S the salinity, satisfying the equations

\begin{array}{l} \frac{\partial u}{\partial t} + \frac{\partial}{\partial x_{i}} (u u_{i}) + \frac{\partial}{\partial z} (u w) + \frac{1}{ρ_{0}} \nabla p + f k \times u \\ (1) & = \frac{\partial}{\partial z} ϑ \frac{\partial u}{\partial z} + \nabla \cdot (K \nabla) u, \\ (2) & \nabla \cdot u + \frac{\partial}{\partial z} w = 0, \\ (3) & \frac{\partial p}{\partial z} = - g ρ, \\ \frac{\partial Θ_{j}}{\partial t} + \frac{\partial}{\partial x_{i}} (u_{i} Θ_{j}) + \frac{\partial}{\partial z} (w Θ_{j}) = \frac{\partial}{\partial z} ϑ_{Θ} \frac{\partial Θ_{j}}{\partial z} \\ (4) & + \nabla (K_{Θ} \nabla) Θ_{j} . \end{array}

Here $i = 1, 2$ , x₁=x, x₂=y, u₁=u, and u₂=v, and summation is implied over the repeating indices i; p is the pressure; and $j = 1, 2$ with Θ₁=T and Θ₂=S represents the potential temperature and salinity, respectively. The seawater density is determined by the equation of state $ρ = ρ (T, S, p)$ , and ρ₀ is the reference density; f is the Coriolis parameter; k is the vertical unit vector; ϑ and K are the coefficients of vertical and horizontal turbulent momentum exchange, respectively; ϑ_Θ and K_Θ are the respective diffusion coefficients; and g is the acceleration due to gravity.

Writing

\begin{matrix} (5) & ρ (x, y, z, t) = ρ_{0} + ρ^{'} (x, y, z, t), \end{matrix}

where ρ^′ is the density fluctuation, we obtain, integrating Eq. (3),

p - p_{atm} = \int_{z}^{ζ} ρ g d z = g ρ_{0} (ζ - z) + g \int_{z}^{ζ} ρ^{'} d z,

where p_atm is the atmospheric pressure. The horizontal pressure gradient is then expressed as the sum of barotropic, baroclinic, and atmospheric pressure gradients:

\begin{array}{l} ρ_{0}^{- 1} \nabla p = g \nabla ζ + g ρ_{0}^{- 1} \nabla I + ρ_{0}^{- 1} \nabla p_{atm}, \\ (6) & I = \int_{z}^{ζ} ρ^{'} d z . \end{array}

Note that horizontal derivatives here are taken at fixed z.

2.2 Turbulent closures

The default scheme to compute the vertical viscosity and diffusivity in the system of Eqs. (1)–(4) is based on the Prandtl–Kolmogorov hypothesis of incomplete similarity. According to it, the turbulent kinetic energy b, the coefficient of turbulent mixing ϑ, and dissipation of turbulent energy ε are connected as $ϑ = l \sqrt{b}$ , where l is the scale of turbulence, ϑ_Θ=c_ρϑ, $ε = c_{ε} b^{2} / ϑ$ , and c_ε=0.046 (Cebeci and Smith, 1974). Prandtl's number c_ρ is commonly chosen as 0.1 and sets the relationship between the coefficients of turbulent diffusion and viscosity. The equation describing the balance of turbulent kinetic energy is obtained by parameterizing the energy production and dissipation in the equation for turbulent kinetic energy b as

\begin{matrix} (7) & \frac{\partial b}{\partial t} - ϑ ({|u_{z}|}^{2} + c_{ρ} g ρ_{0}^{- 1} \frac{\partial ρ}{\partial z}) + c_{ε} b^{2} / ϑ = α_{b} \frac{\partial}{\partial z} ϑ \frac{\partial b}{\partial z}, \end{matrix}

with the boundary conditions

\begin{array}{l} b |_{h} = B_{1} | u |^{2}, \\ ϑ b_{z} |_{ζ} = γ_{ζ} u_{* ζ}^{3}, \end{array}

where α_b=0.73, B₁=16.6, and $γ_{ζ} = 0.4 \times 10^{- 3}$ ; $u_{* ζ} = (ρ / ρ_{a})^{1 / 2} u_{*}$ is the dynamical velocity in water near the surface, ρ_a the air density, and u_* the dynamic velocity of water on the interface between air and water.

Dissipative term is written as

b^{2} / ϑ = (2 b^{ν + 1} b^{ν} - (b^{2})^{ν}) / ϑ^{ν},

where ν is the index of iterations. Equation (7) is solved by a three-point Thomas scheme in the vertical direction with the boundary conditions given above. Iterations are carried out until convergence determined by the condition

max |(b^{ν + 1} - b^{ν}) / b^{ν}| < ϖ,

where ϖ is a small value O(10⁻⁶). More details on the solution of this equation are given in Voltzinger et al. (1989).

To determine the turbulence scale l in the presence of surface and bottom boundary layers we use the Montgomery formula (Reid, 1957):

l = \frac{κ}{H} Z_{h} Z_{ζ},

where $H = h + ζ$ is the full water depth, $Z_{h} = z + h + z_{h}$ , $Z_{ζ} = - z + ζ + z_{ζ}$ , κ≃0.4 is the von Kármán constant, z the layer depth, and z_h and z_ζ are the roughness parameters for the bottom and free surface, respectively. To remove turbulent mixing in layers that are distant from interfaces we modify the Montgomery formula by introducing the cutoff function $Z_{0} = 1 - β_{1} H^{- 2} Z_{h} Z_{ζ}$ , $0 \leq β_{1} \leq 4$ (Voltzinger, 1985):

l = \frac{κ}{H} Z_{h} Z_{ζ} Z_{0} .

In addition to the default scheme, one may select a scheme provided by the General Ocean Turbulence Model (GOTM) (Burchard et al., 1999) implemented into the FESOM-C code for computing vertical eddy viscosity and diffusion for momentum and tracer equations. GOTM includes large number of well-tested turbulence models with at least one member of every relevant model family (empirical models, energy models, two-equation models, algebraic stress models, K-profile parameterizations, etc.) and treats every single water column independently. An essential part of GOTM includes one-point second-order schemes (Umlauf and Burchard, 2005; Umlauf et al., 2005, 2007).

2.3 Bottom friction parameterization

The model uses either a constant bottom friction coefficient C_d, or it is computed through the specified bottom roughness height z_h. The first option is preferable if the vertical resolution everywhere in the domain does not resolve the logarithmic layer or when the vertically averaged equations are solved. In the second option the bottom friction coefficient is computed according to Blumberg and Mellor (1987) and has the following form:

C_{d} = {(\ln ((0.5 h_{b} + z_{h}) / z_{h}) / κ)}^{- 2},

where h_b is the thickness of the bottom layer. It is also possible to prescribe C_d or z_h as a function of the horizontal coordinate at the initialization step.

2.4 Boundary conditions

The boundary conditions for the dynamical Eqs. (1)–(2) are those of no-slip on the solid boundary ∂Q₁,

u |_{\partial Q_{1}} = 0 .

As is well known, the formulation of open boundary conditions faces difficulties. They are related to either the lack or incompleteness of information demanded by the theory, for example on velocity components at the open boundary. Furthermore, whatever the external information, it may contradict the solution inside the computational domain, leading to instabilities that are frequently expressed as small-scale vortex structures forming near the open boundary. The procedure reconciling the external information with the solution inside the domain becomes of paramount importance. We use two approaches. The first one is to use a function whereby advection and horizontal diffusion are smoothly tapered to zero in the close vicinity of open boundary ∂Q₂. Such tapering makes the equations quasi-hyperbolic at the open boundary so that the formulation of one condition (for example, for the elevation, $ζ |_{\partial Q_{2}} = ζ_{Γ}$ ) is possible (Androsov et al., 1995).

The other approach is to adapt the external information. It is applied to scalar fields and will be explained further.

Note that despite simplifications, barotropic and baroclinic perturbations may still disagree at the open boundary, leading to instabilities in its vicinity. In this case an additional buffer zone is introduced with locally increased horizontal diffusion and bottom friction.

Dynamic boundary conditions on the top and bottom specify the momentum fluxes entering the ocean. Neglecting the contributions from horizontal viscosities, we write

\begin{array}{l} ϑ \frac{\partial u}{\partial z} |_{ζ} = τ_{ζ} / ρ_{0}, \\ ϑ \frac{\partial u}{\partial z} |_{- h} = τ_{h} / ρ_{0} = C_{d} | u_{h} | u_{h} . \end{array}

The first one sets the surface momentum flux to the wind stress at the surface (τ_ζ), and the second one sets the bottom momentum flux to the frictional flux at the bottom (τ_h), with u_h the bottom velocity.

Now we turn to the boundary conditions for the scalar quantities obeying Eq. (4). This is a three-dimensional parabolic equation and the boundary conditions are determined by its leading (diffusive) terms. We impose the no-flux condition on the solid boundary ∂Q₁ and the bottom $z = - h$ .

The conditions at the open boundary ∂Q₂ are given for outflow and inflow; see, e.g., Barnier et al. (1995) and Marchesiello et al. (2001):

Θ_{t} + a Θ_{x} + b Θ_{y} = - \frac{1}{τ} (Θ - Θ_{Γ}),

where Θ_Γ is the given field value, usually a climatological one or relying on data from a global numerical model or observations. If the phase velocity components and $a = - Θ_{t} Θ_{x} / G$ , $b = - Θ_{t} Θ_{y} / G$ , and $G = [(\partial Θ / \partial x)^{2} + (\partial Θ / \partial y)^{2}]^{- 1}$ (Raymond and Kuo, 1984), Θ propagates out of the domain, and one sets τ=τ₀. If it propagates into the domain, a and b are set to zero and τ=τ_Γ, with τ_Γ≪τ₀. The parameter τ is determined experimentally and commonly is from hours to days. In the FESOM-C such an adaptive boundary condition is routinely applied for temperature and salinity, yet it can also be used for any components of a solution.

At the surface the fluxes are due to the interaction with the atmosphere:

\begin{array}{l} (8) & ϑ_{Θ} \frac{\partial T}{\partial z} |_{ζ} = \hat{Q} (x, y, t) / ρ_{0} c_{p}, \\ (9) & ϑ_{Θ} \frac{\partial S}{\partial z} |_{ζ} = 0, \end{array}

where $\hat{Q}$ is the heat flux excluding shortwave radiation, which has been included as a volume heat source in the temperature equation, with c_p the specific heat of seawater. The impact of the precipitation–evaporation has been included as a volume source in the continuity equation. In the presence of rivers, their discharge is added either as a prescribed inflow at the open boundary in the river mouth or as volume sources of mass, heat, and momentum distributed in the vicinity of the open boundary. In the first case it might create an initial shock in elevation, so the second method is safer.

3 Temporal discretization

As is common in coastal models, we split the fast and slow motions into, respectively, barotropic and baroclinic subsystems (Lazure and Dumas, 2008; Higdon, 2008; Gadd, 1978; Blumberg and Mellor, 1987; Deleersnijder and Roland, 1993). The reason for this splitting is that surface gravity waves (external mode) are fast and impose severe limitations on the time step, whereas the internal dynamics can be computed with a much larger time step. The time step for the external mode τ_2-D is limited by the speed of surface gravity waves, and that for the internal mode, τ_3-D, by the speed of internal waves or advection. The ratio $M_{t} = τ_{3 - D} / τ_{2 - D}$ depends on applications, but is commonly between 10 and 30. In practice, additional limitations are due to vertical advection or wetting and drying processes. We will further use the indices k and n to enumerate the internal and external time steps, respectively.

The numerical algorithm passes through several stages. In the first stage, based on the current temperature and salinity fields (time step k), the pressure is computed from hydrostatic equilibrium Eq. (6) and then used to compute the baroclinic pressure gradient $ρ_{0}^{- 1} \nabla p = g \nabla ζ^{k} + g ρ_{0}^{- 1} \nabla I^{k} + ρ_{0}^{- 1} \nabla p_{atm}$ . We use an asynchronous time stepping, assuming that integration of temperature and salinity is half-step shifted with respect to momentum. The index k on I implies that it is centered between k and k+1 of momentum integration. The elevation in the expression above is taken at time step k, which makes the entire estimate for ∇p only first-order accurate with respect to time.

At the second stage, the predictor values of the three-dimensional horizontal velocity are determined as

\begin{array}{l} {\tilde{u}}^{k + 1} - u^{k} = τ_{3 - D} {(- f k \times u - \nabla \cdot u u + \nabla \cdot (K \nabla u))}^{AB 3} \\ - τ_{3 - D} ρ_{0}^{- 1} \nabla p + τ_{3 - D} \partial_{z} ϑ \partial_{z} {\tilde{u}}^{k + 1} - τ_{3 - D} \partial_{z} {(w u)}^{AB 3} . \end{array}

Here K is the coefficient of horizontal viscosity, and AB3 implies the third-order Adams–Bashforth estimate. The horizontal viscosity operator can be made biharmonic or replaced with filtering as discussed in the next chapter.

To carry out mode splitting, we write the horizontal velocity as the sum of the vertically averaged one $\overline{u}$ and the deviation thereof (pulsation) u^′:

\begin{array}{l} u = \overline{u} + u^{'}, \\ \overline{u} = \frac{1}{H} \int_{- h}^{ζ} u d z, \\ \int_{- h}^{ζ} u^{'} = 0 . \end{array}

By integrating the system in Eqs. (1)–(3) vertically between the bottom and surface, with regard for the kinematic boundary conditions $\partial_{t} ζ + u \nabla ζ = w$ on the surface, $- u \nabla h = w$ at the bottom, and time discretization, we get

\begin{array}{l} (ζ^{n + 1} - ζ^{n}) + τ_{2 - D} \nabla {(H \overline{u})}^{AB 3} = 0, \\ {\overline{u}}^{n + 1} - {\overline{u}}^{n} = τ_{2 - D} {(- f k \times \overline{u} - \nabla \cdot \overline{u} \overline{u} + \nabla \cdot (K \nabla \overline{u}))}^{AB 3} \\ - τ_{2 - D} (g \nabla ζ)^{AM 4} + τ_{2 - D} (τ_{ζ} / ρ_{0} - τ_{h}) - τ_{2 - D} R_{3 - D}^{'} \\ - τ_{2 - D} g ρ_{0}^{- 1} \nabla {\overline{I}}^{k} . \end{array}

Here a specific version of AB3 is used: ${\overline{u}}^{AB 3} = (3 / 2 + β) {\overline{u}}^{k} - (1 / 2 + 2 β) {\overline{u}}^{k - 1} + β {\overline{u}}^{k - 2}$ , with β=0.281105 for stability reasons (Shchepetkin and McWilliams, 2005); AM4 implies the Adams–Multon estimate $ζ^{AM 4} = δ ζ^{n + 1} + (1 - δ - γ - ϵ) ζ^{n} + γ ζ^{n - 1} + ϵ ζ^{n - 2}$ , taken with δ=0.614, γ=0.088, and ϵ=0.013 (Shchepetkin and McWilliams, 2005). In the equations above τ_ζ and τ_h are the surface (wind) and bottom stresses, respectively, and $\nabla \overline{I}$ is the vertically integrated gradient of baroclinic pressure. The term $R_{3 - D}^{'}$ contains momentum advection and horizontal dissipation of the pulsation velocity integrated vertically:

R_{3 - D}^{'} = \frac{1}{H^{k}} {[\int_{- h}^{ζ} \nabla \cdot u^{'} u^{'} - \int_{- h}^{ζ} \nabla \cdot (K \nabla u^{'})]}^{k} .

In this expression H^k is the total fluid depth at time step k, $H^{k} = h + ζ^{k}$ . This term is computed only on the baroclinic time step and kept constant through the integration of the internal mode.

The bottom friction is taken as

τ_{h} = C_{d} {|\overline{u}|}^{n} ({\overline{u}}^{n + 1} - {\overline{u}}^{n}) / H^{n + 1} + C_{d} | {\tilde{u}}_{h} |^{k + 1} {\tilde{u}}_{h}^{k + 1} / H^{n + 1} .

The first part of bottom friction is needed to increase stability, while the second part estimates the correct friction, with ${\tilde{u}}_{h}^{k + 1}$ the horizontal velocity vector in the bottom cell on the predictor time step.

The system of vertically averaged equations is stepped explicitly (except for the bottom friction) through M_t time steps of duration τ_2-D (index n) to “catch up” the k+1 baroclinic time step. The update of elevation is made first, followed by the update of vertically integrated momentum equations.

At the “corrector” step, the 3-D velocities are corrected to the surface elevation at k+1:

u^{k + 1} = \frac{△_{i}^{k}}{△_{i}^{k + 1}} {\tilde{u}}^{k + 1} + ({\overline{u}}^{k + 1} - {\overline{u}}^{P}),

with ${\overline{u}}^{P} = \frac{1}{H^{k + 1}} \sum_{- h}^{ζ} ({\tilde{u}}^{k + 1} △_{i}^{k})$ , where i is the vertical index. Here $△_{i}^{k}$ and $△_{i}^{k + 1}$ are the thicknesses of the ith layer calculated on respective baroclinic time steps. The layer thickness is $△_{i}^{k} = △_{i} H^{k}$ , where △_i is the unperturbed vertical grid spacing. This correction removes the barotropic component of the predicted velocity and combines the result with the computed barotropic velocity. We will suppress the layer index i where it is unambiguous.

The final step in the dynamical part calculates the transformed vertical velocity w^k+1 from 3-D continuity Eq. (2). It is used in the next predictor step. Note that in the predictor step the computations of vertical viscosity are implicit.

New horizontal velocities, the so-called “filtered” ones, are used for advection of a tracer. They are given by the sum of the filtered depth mean and the baroclinic part of the “predicted” velocities (Deleersnijder, 1993),

u_{F}^{k + 1} = \frac{△_{i}^{k}}{△_{i}^{k + 1}} {\tilde{u}}^{k + 1} + ({\overline{u}}^{F} - {\overline{u}}^{P}),

with ${\overline{u}}^{F} = \frac{1}{M_{t} H^{k + 1}} \sum_{n = 1}^{n = M_{t}} ({\overline{u}}^{n} H^{n})$ . The procedure of “filtering” removes possible high-frequency components in the barotropic velocity. It also improves accuracy because it in essence works toward centering the contribution of the elevation gradient. Once the filtered velocity is computed, the vertical velocity is updated to match it.

The equation for temperature is taken in the conservation form

\begin{array}{l} △^{k + 1} T^{k + 1} \\ = △^{k} T^{k} - τ_{3 - D} [\nabla \cdot (u_{F}^{k + 1} △^{k} T^{*}) + w_{Ft} T_{t}^{*} - w_{Fb} T_{b}^{*}] \\ + D + τ_{3 - D} R + τ_{3 - D} C, \end{array}

where D combines the terms related to diffusion, w_Ft,w_Fb, T_t, and T_b are the vertical transport velocity and temperature on top and bottom of the layer, and T^* is computed through the second-order Adams–Bashforth estimate. R is the boundary thermal flux (either from the surface or due to river discharge). The last term in the equation above is

C = T \frac{\partial △_{i}}{\partial t} = - T (\nabla {(△_{i} u_{F})}^{k + 1} + \partial_{z} {(△_{i} w_{F})}^{k + 1}) .

Its two constituents combine to zero because of continuity. Keeping this term makes sense if computations of advection are split into horizontal and vertical substeps. The salinity is treated similarly.

In simulations of coastal dynamics it is often necessary to simulate flooding and drying events. Explicit time stepping methods of solving the external mode are well suited for this (Luyten et al., 1999; Blumberg and Mellor, 1987; Shchepetkin and McWilliams, 2005). The algorithm to account for wetting and drying will be presented in the next section. We only note that computations are performed on each time step of the external mode.

4 Spatial discretization

In the finite-volume method, the governing equations are integrated over control volumes, and the divergence terms, by virtue of the Gauss theorem, are expressed as the sums of respective fluxes through the boundaries of control volumes. For the cell-vertex discretization the scalar control volumes are formed by connecting cell centroids with the centers of edges, which gives the so-called median-dual control volumes around mesh vertices. The vector control volumes are the mesh cells (triangles or quads) themselves, as schematically shown in Fig. 1.

https://www.geosci-model-dev.net/12/1009/2019/gmd-12-1009-2019-f01

Figure 1Schematic of mesh structure. Velocities are located at centroids (red circles) and elevation at vertices (blue circles). A scalar control volume associated with vertex v₁ is formed by connecting neighboring centroids to edge centers. The control volumes for velocity are the triangles and quads themselves. The lines passing through two neighboring centroids (e.g., c₁ and c₂) are broken in a general case at edge centers. Their fragments are described by the left and right vectors directed to centroids (s_l and s_r for edge e). Edge e is defined by its two vertices v₁ and v₂ and is considered to be directed to the second vertex. It is also characterized by two elements c₁ and c₂ to the left and to the right, respectively.

FESOM-C v.2: coastal dynamics on hybrid unstructured meshes

2.1 The governing equations

2.2 Turbulent closures

2.3 Bottom friction parameterization

2.4 Boundary conditions

4.1 Divergence and gradients

4.2 Momentum advection

4.3 Tracer advection

4.4 Viscosity and filtering

4.5 Wetting and drying algorithm

5.1 Sylt–Rømø experiment

5.2 Southeast North Sea circulation

6.1 Triangles vs. quads: numerical performance

6.2 Triangles vs. quads: open boundaries