A flexible z-layers approach for the accurate representation of free surface flows in a coastal ocean model (SHYFEM v. 7_5_71)

We propose a discrete multilayer shallow water model based on z-layers, which, thanks to the insertion and removal of surface layers, can deal with an arbitrarily large tidal oscillation independently of the vertical resolution. The algorithm is based on a classical two-step procedure used in numerical simulations with moving boundaries (grid movement followed by a grid topology change, that is, the insertion/removal of surface layers), which avoids the appearance of surface layers with very small or negative thickness. With ad hoc treatment of advection terms at nonconformal edges that may appear owing to insertion/removal operations, mass conservation and the compatibility of the tracer equation with the continuity equation are preserved at a discrete level. This algorithm called z-surface-adaptive, can be reduced, as a particular case when all layers are moving, to the z-star coordinate. With idealized and realistic numerical experiments, we compare the z-surface-adaptive against z-star and we show that it can be used to simulate coastal flows effectively.

1 Introduction

The accuracy of ocean models in reproducing many dynamical processes is highly related to their vertical coordinate system. In the literature, many choices exist covering the spectrum of coordinate systems. There are four main types of vertical coordinates that correspond to different vertical subdivisions of the fluid domain: (1) isopycnal layers with the interfaces that are material surfaces (Lagrangian framework); (2) z-layers with fixed interfaces parallel to geopotentials (Eulerian framework); (3) terrain/surface-following σ or s-layers with interfaces adapted to the ocean surface and bottom boundaries; (4) adaptive coordinate with interfaces that dynamically adapt to better capture different flow features (Lagrangian tendencies, stratification, and shear). The last two types move “arbitrarily” with respect to the flow, either to adapt to the free surface or any other features, and belong to the Arbitrary Lagrangian Eulerian (ALE) framework.

z-layers were used in early ocean and coastal models and are nowadays implemented and used in some ocean models (HAMburg Shelf Ocean Model, HAMSOM, Backhaus, 1985; Tidal, Residual, Intertidal Mudflat-3D, TRIM-3D, Cheng et al., 1993; UNTRIM-3D, Casulli and Walters, 2000; Shallow water HYdrodynamic Finite Element Model, SHYFEM, Umgiesser, 2022). They are attractive when simulating strongly stratified flows (Hordoir et al., 2015) and low-frequency motions (Leclair and Madec, 2011). This occurs because the isopycnals are well aligned to the z-interfaces or they slowly depart from them. At the same time, the truncation error of the internal pressure gradient term remains very small.

A vertical discretization based on fixed interfaces is expected to have issues with the complex and moving boundaries represented by the free surface and by the ocean bottom. In this paper, we focus on z-layers performances relative to the treatment of the free surface boundary. To simplify the boundary condition at the free surface, z-layers were typically coded, allowing the surface layer to vary in thickness (Griffies et al., 2001). However, in such models, the surface layer cannot vanish, which implies that the free surface variation must be smaller than the surface layer thickness. For coastal applications, this is a serious drawback, especially for the vertical resolution in shallow areas with high tidal elevations. In order to overcome this problem, other z-type coordinates have been introduced over the years: they are based on z-layers that move to accommodate the tidal oscillation, but the bottom is not a coordinate surface (they are surface-following but not terrain-following). These coordinates are clearly of ALE-type but in the ocean modeling literature they are classified as z because the deviation from the geopotentials is very small. They combine small diapycnal mixing, especially for internal tides computations, and small truncation error on the pressure gradient term. The z-star of Adcroft and Campin (2004), the quasi-z of Mellor et al. (2002), and the hybrid $z/\mathit{\sigma }$ of Burchard and Petersen (1997) all belong to such z-surface-following systems. An alternative to deal with the moving surface is to keep the vertical grid perfectly aligned to geopotentials, thus working in a truly Eulerian framework, but allowing the surface layer(s) to be removed or inserted. We refer to this system as z-surface-adaptive. Insertion/removal of the surface layer has been discussed in Casulli and Cheng (1992) and it is used for example in Burchard and Baumert (1998). However “both the accuracy and stability are suspect; it is most likely difficult to make the transition of a vanishing layer smooth enough to not generate numerical problems; conservation issues are a major concern and the likelihood of vanishing layers becomes more frequent with increasing vertical resolution” (Adcroft and Campin, 2004).

In this paper we propose an algorithm for the z-surface adaptive coordinate that goes beyond such limitations. We employ a classical grid adaptation strategy for situations in which the adaptation is driven by a moving boundary (Guardone et al., 2011). It combines a first ALE grid movement step (surface interface displacement stretched by the free surface displacement) and a second topology modification step (layer insertion, layer removal). All these operations are easily performed on the one-dimensional vertical grid. If the water depth is positive, the thickness of the surface layers remains positive, avoiding stability issues related to the appearance of small or negative layers. We show that the mass is conserved. Also, the discrete preservation of a constant tracer can be easily accomplished, which guarantees complete consistency at a discrete level of the tracer equation with the continuity equation as shown since the work of Lin and Rood (1996); Gross et al. (2002).

This solution generalizes z-layers in the sense that the same algorithm can be easily reduced to z-star and can be added to a flexible vertical coordinate system. In fact, the grid adaptation has one free parameter that controls the number of moving layers. Tuning such a parameter, so that all the layers along the water column are moving, we show the link of the proposed approach with the z-star.

The algorithm is implemented in the SHYFEM finite-element ocean model of the CNR-ISMAR (Umgiesser et al., 2004, https://github.com/SHYFEM-model/shyfem, last access: 5 November 2023), which implements the multilayer shallow water equations with z and σ layers. SHYFEM uses a semi-implicit finite element discretization on unstructured B-type grids derived from the work of Casulli and Cheng (1992) and Williams and Zienkiewicz (1981).

The paper is organized as follows: in Sect. 2 we introduce the vertical discretization and the multilayer shallow water model. Three different vertical discretizations are considered: the standard multilayer shallow water model based on σ-layers, then the z-star and the z-layers. In Sect. 3 we provide the semi-implicit finite element discretization of the multilayer equations. In Sect. 4 we describe the z-surface-adaptive algorithm, in Sect. 5 we detail the issue of a spatially variable number of surface layers caused by the insertion/removal operations. In Sect. 6 we provide numerical tests and in Sect. 7 we conclude with a discussion.

2 Multilayer shallow water model

We start considering the multilayer (or layer integrated) shallow water model for stratified flows studied in Audusse et al. (2011). The space variable is $(\mathit{x},z)\in {\mathbb{R}}^{\mathrm{3}}$ with $\mathit{x}=(x,y)\in {\mathbb{R}}^{\mathrm{2}}$ that denotes the horizontal space variable. We consider the fluid domain Ω:

$$\mathrm{\Omega }=\mathit{\left\{}\right(\mathit{x},z):\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}},-z_{\mathrm{b}}(\mathit{x})\le z\le \mathit{\zeta }(\mathit{x},t\left)\mathit{\right\}},$$

where Ω_x is the projection of Ω onto the horizontal plane, ζ(x,t) is a function that represents the free surface elevation, and z_b(x) is the bathymetry, which does not depend on time. The water depth is $H(\mathit{x},t)=\mathit{\zeta }(\mathit{x},t)+z_{\mathrm{b}}\left(\mathit{x}\right)$. As depicted in Fig. 1, right panel, the multilayer shallow water model is based on a discretization of the domain Ω with a vertical grid composed of N layers denoted Ω_α with $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$, ordered from the free surface to the bottom. The layers are non-overlapping with $\mathrm{\Omega }={\bigcup }_{\mathit{\alpha }=\mathrm{1}}^N{\mathrm{\Omega }}_{\mathit{\alpha }}$. Each layer Ω_α is delimited laterally by the vertical domain boundary and in the vertical by the time-dependent interfaces ${\mathrm{\Gamma }}_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}\left(t\right)$ defined by the set of points of coordinates (x,z) such that $z=z_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}(\mathit{x},t)$. The free surface Γ^ζ and the bottom interfaces Γ^b are described by the free surface elevation $z_{\mathrm{1}/\mathrm{2}}=\mathit{\zeta }(\mathit{x},t)$ and by the bathymetry function $z_{N+\mathrm{1}/\mathrm{2}}=-z_{\mathrm{b}}\left(\mathit{x}\right)$ respectively. In order to provide the rules for such slicing of the domain, we define a reference domain that is constant in time, with space variables $(\mathit{x},s)\in {\mathbb{R}}^{\mathrm{3}}$ such that:

$${\mathrm{\Omega }}^{\mathrm{0}}=\mathit{\left\{}\right(\mathit{x},s):\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}},-\mathrm{1}\le s\le \mathrm{0}\mathit{\}},$$

and discretized by means of a vertical grid similarly composed of N layers, each denoted ${\mathrm{\Omega }}_{\mathit{\alpha }}^{\mathrm{0}}$. The reference layers are delimited vertically by the fixed-in-time interfaces ${\mathrm{\Gamma }}_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}^{\mathrm{0}}$, which are placed at the vertical coordinate given by the coefficients $s_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}$. Such constants can be ordered:

$$s_{\mathrm{1}/\mathrm{2}}=\mathrm{0}>s_{\mathrm{2}-\mathrm{1}/\mathrm{2}}>\mathrm{\dots }>s_{N+\mathrm{1}/\mathrm{2}}=-\mathrm{1}.$$

Then, the interface position can be obtained by mapping the reference interface ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{0}}$ to the actual or physical interface ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\left(t\right)$. In general we assume that there exists a function, for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(1)}& A:{\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{0}}\to {\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\left(t\right),z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=A(\mathit{x},s_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}},t)\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}}.\end{array}$$

To prescribe this function we use the generalized vertical coordinate transformation, see Mellor et al. (2002):

$$\begin{array}{}\text{(2)}& z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=\mathit{\zeta }(\mathit{x},t)+s_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\left(\mathit{\zeta }(\mathit{x},t)+z_b\left(\mathit{x}\right)\right),\end{array}$$

which assures a surface- and terrain-following grid that is limited by the interfaces ${\mathrm{\Gamma }}_{\mathrm{1}/\mathrm{2}}\left(t\right)={\mathrm{\Gamma }}^{\mathit{\zeta }}\left(t\right)$ and ${\mathrm{\Gamma }}_{N+\mathrm{1}/\mathrm{2}}={\mathrm{\Gamma }}^b$. The reference and the physical domains with their vertical subdivisions are sketched in Fig. 1.

Figure 1 One-dimensional sketch of the reference (a) and physical (b) domains for the multilayer shallow water model.

Using this transformation, the layer thickness can be deduced from the water depth, for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(3)}& {\displaystyle }{h}_{\mathit{\alpha }}(\mathit{x},t)& {\displaystyle }=z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}(\mathit{x},t)-z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}(\mathit{x},t)\text{(4)}& {\displaystyle }& {\displaystyle }=\left(s_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-s_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}\right)H(\mathit{x},t)=l_{\mathit{\alpha }}H(\mathit{x},t),\end{array}$$

where the coefficients $l_{\mathit{\alpha }}=s_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-s_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$ are prescribed after the creation of the reference grid. They are positive and they sum to one ${\sum }_{\mathit{\alpha }=\mathrm{1}}^N{l}_{\mathit{\alpha }}=\mathrm{1}$. The multilayer model is based on a piecewise constant approximation, on the vertical grid, of the horizontal fluid velocity and of a generic tracer. For $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{u}}_{\mathit{\alpha }}(\mathit{x},t)={\displaystyle \frac{\mathrm{1}}{h_{\mathit{\alpha }}}}\underset{z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}}{\overset{z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}{\int }}\mathit{u}(\mathit{x},z,t)\mathrm{d}z,\text{(6)}& {\displaystyle }{T}_{\mathit{\alpha }}(\mathit{x},t)={\displaystyle \frac{\mathrm{1}}{h_{\mathit{\alpha }}}}\underset{z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}}{\overset{z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}{\int }}T(\mathit{x},z,t)\mathrm{d}z.\end{array}$$

The tracer for us will be the salinity. We assume that the fluid density depends on salinity through the UNESCO equation of state (Millero and Poisson, 1981) at one atmosphere and at a constant potential temperature of 12.4 ^∘C. If the equation of state is of type ρ=ρ(T), the density vertical discretization derives from the tracer one, for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(7)}& {\mathit{\rho }}_{\mathit{\alpha }}(\mathit{x},t)=\mathit{\rho }\left(T_{\mathit{\alpha }}\right(\mathit{x},t\left)\right).\end{array}$$

We introduce the following notation for a generic function f(z):

• To express a function that is discontinuous at the interface, we use the same notation as Fernández-Nieto et al. (2014):

  $$f_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}={\left({\left.f\right|}_{{\mathrm{\Omega }}_{\mathit{\alpha }}}\right)}_{{\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}},f_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}={\left({\left.f\right|}_{{\mathrm{\Omega }}_{\mathit{\alpha }-\mathrm{1}}}\right)}_{{\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}.$$

• If the function is continuous

  $$f_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=f_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}=f_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}={\left.f\right|}_{{\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}.$$

• The difference of the function between the upper and lower interface is

  $$[f{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=f_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-f_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}.$$

Mass conservation reads

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}+\mathrm{\nabla }\cdot \left(\sum _{\mathit{\beta }=\mathrm{1}}^N{h}_{\mathit{\beta }}{\mathit{u}}_{\mathit{\beta }}\right)=\mathrm{0}.\end{array}$$

In this work we consider the multilayer shallow water model for stratified fluid with the Boussinesq assumption. Momentum and tracer equations in the multilayer approach can be written for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial h_{\mathit{\alpha }}{\mathit{u}}_{\mathit{\alpha }}}{\partial t}}& +\mathrm{\nabla }\cdot \left(h_{\mathit{\alpha }}{\mathit{u}}_{\mathit{\alpha }}\otimes {\mathit{u}}_{\mathit{\alpha }}\right)=[\mathit{u}G{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-g{h}_{\mathit{\alpha }}\mathrm{\nabla }\mathit{\zeta }\\ & +[\mathit{K}{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}+{\mathit{B}}_{\mathit{\alpha }},\end{array}\text{(10)}& {\displaystyle \frac{\partial h_{\mathit{\alpha }}{T}_{\mathit{\alpha }}}{\partial t}}+\mathrm{\nabla }\cdot \left(h_{\mathit{\alpha }}{T}_{\mathit{\alpha }}{\mathit{u}}_{\mathit{\alpha }}\right)=[TG{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}+[K_T{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}},\end{array}$$

where $G_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}$ is the mass-transfer function responsible for the vertical mass exchange between the layers, ${\mathit{K}}_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}$ are the vertical viscous fluxes that model the shear stress between the layers, $K_{T,\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}$ are the vertical diffusive fluxes that model the diffusive process between the layers, and B_α models the pressure force related to the buoyancy gradient. The system Eqs. (8), (9) and (10) is implemented in the SHYFEM model, as well as in many other ocean models (Burchard and Petersen, 1997; Klingbeil et al., 2018). If N is the number of vertical layers, the equations are solved for 2N+1 unknown variables, which are the free surface elevation, the layer discharges h_αu_α, and the layer-integrated tracer h_αT_α. The layer thickness is deduced from the water depth through Eq. (4). In the following, we give the details of the SHYFEM implementation of each term on the right-hand side.

From the derivation of Fernández-Nieto et al. (2014), the definition of the mass-transfer function is

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{G}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}& =\left(\mathrm{\nabla }{z}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\cdot {\mathit{u}}_{\mathit{\alpha }}\right)+{\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}\\ & =\left(\mathrm{\nabla }{z}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\cdot {\mathit{u}}_{\mathit{\alpha }-\mathrm{1}}\right)+{\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}\end{array},\end{array}$$

with ${\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ the velocity of the grid interface:

$$\begin{array}{}\text{(12)}& {\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\partial z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}{\partial t}},\end{array}$$

and $w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\pm }$ the vertical fluid velocity at the interface. The vertical velocity is computed from the following relationships:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{w}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}& =-w_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{-}-h_{\mathit{\alpha }}\mathrm{\nabla }\cdot {\mathit{u}}_{\mathit{\alpha }}\mathrm{and}{w}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}=w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}\\ & +\mathrm{\nabla }{z}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\cdot \left({\mathit{u}}_{\mathit{\alpha }}-{\mathit{u}}_{\mathit{\alpha }-\mathrm{1}}\right),\end{array}\end{array}$$

which are evaluated starting from the bottom $\mathit{\alpha }=N,\mathrm{\dots },\mathrm{1}$, where the no-slip condition is imposed $w_{N+\mathrm{1}/\mathrm{2}}^{-}={\mathit{u}}_N\cdot \mathrm{\nabla }{z}_{\mathrm{b}}$. In practice, and as is standard in ocean models, the mass-transfer function is computed directly from the layer-integrated mass equation

$$\begin{array}{}\text{(14)}& G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=G_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}+{\displaystyle \frac{\partial h_{\mathit{\alpha }}}{\partial t}}+\mathrm{\nabla }\cdot \left(h_{\mathit{\alpha }}{\mathit{u}}_{\mathit{\alpha }}\right).\end{array}$$

Summing from N to α as

$$\begin{array}{}\text{(15)}& G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=G_{N+\mathrm{1}/\mathrm{2}}+\sum _{\mathit{\beta }=N}^{\mathit{\alpha }}{\displaystyle \frac{\partial h_{\mathit{\beta }}}{\partial t}}+\sum _{\mathit{\beta }=N}^{\mathit{\alpha }}\mathrm{\nabla }\cdot \left(h_{\mathit{\beta }}{\mathit{u}}_{\mathit{\beta }}\right),\end{array}$$

which implies $G_{\mathrm{1}/\mathrm{2}}=\mathrm{0}$ or no mass loss at the free surface. The vertical velocity at the interfaces $w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\pm }$ no longer appears in the system, but it can be computed from the incompressibility condition Eq. (13) in a post-processing step. With a horizontal velocity and tracer that are discontinuous at the interfaces, the vertical momentum flux in Eq. (9) is computed with a numerical flux. An upwind flux is used in this study, for ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ it reads

$$G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}{\mathit{u}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}{\mathit{u}}_{\mathit{\alpha }}+G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}{\mathit{u}}_{\mathit{\alpha }-\mathrm{1}},$$

with $G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+}=max(\mathrm{0},G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}})$ and $G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}=min(\mathrm{0},G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}})$. For the tracer, a total variation diminishing (TVD) flux is employed (LeVeque, 2002).

The terms ${\mathit{K}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ and $K_{T,\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ are the vertical viscous and diffusive fluxes computed at the interface ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$:

$$\begin{array}{rl}& {\mathit{K}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}={\mathit{\nu }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}{D}_z{\mathit{u}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}},\\ & K_{T,\mathit{\alpha }-\mathrm{1}/\mathrm{2}}={\mathit{\nu }}_{T,\mathit{\alpha }-\mathrm{1}/\mathrm{2}}{D}_z{T}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\end{array},$$

where ${\mathit{\nu }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ is the vertical viscosity and ${\mathit{\nu }}_{T,\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ the vertical diffusivity. D_z(⋅) is an approximation of the vertical derivative evaluated at the interface and resolved with finite differences. The vertical viscosity and diffusivity can be laminar or computed with a turbulence model. The bottom momentum flux is specified by a quadratic friction model. Then, the viscous fluxes read

$${\mathit{K}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=\left\{\begin{array}{ll}\mathrm{0},& \mathit{\alpha }=\mathrm{1}\\ {\mathit{\nu }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\frac{{\mathit{u}}_{\mathit{\alpha }-\mathrm{1}}-{\mathit{u}}_{\mathit{\alpha }}}{(h_{\mathit{\alpha }-\mathrm{1}}+h_{\mathit{\alpha }})/\mathrm{2}},& \mathit{\alpha }=\mathrm{2},\mathrm{\dots },N\\ {\mathit{\tau }}_b=-C_{\mathrm{F}}\left|{\mathit{u}}_N\right|{\mathit{u}}_N,& \mathit{\alpha }=N+\mathrm{1}\end{array}\right.,$$

with C_F the bottom friction coefficient. Similarly, the diffusive fluxes read

$$K_{T,\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=\left\{\begin{array}{ll}\mathrm{0},& \mathit{\alpha }=\mathrm{1}\\ {\mathit{\nu }}_{T,\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\frac{T_{\mathit{\alpha }-\mathrm{1}}-T_{\mathit{\alpha }}}{(h_{\mathit{\alpha }-\mathrm{1}}+h_{\mathit{\alpha }})/\mathrm{2}},& \mathit{\alpha }=\mathrm{2},\mathrm{\dots },N\\ \mathrm{0},& \mathit{\alpha }=N+\mathrm{1}\end{array}\right.,$$

with no tracer fluxes through the free surface and the bottom.

Finally, the term B_α represents the internal pressure gradient force. The layer-integrated pressure gradient term ${\int }_{z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}}^{z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}\mathrm{\nabla }p\left(z\right)\mathrm{d}z$, instead of applying the Leibniz rule (Audusse et al., 2011), has been split into the external pressure gradient, related to the free surface slope, and the internal pressure gradient, related to the buoyancy gradient. The internal pressure gradient term is written in the density Jacobian form of Song (1998):

$${\mathit{B}}_{\mathit{\alpha }}=h_{\mathit{\alpha }}{b}_{\mathrm{1}}\mathrm{\nabla }\mathit{\zeta }+h_{\mathit{\alpha }}\sum _{\mathit{\beta }=\mathrm{1}}^{\mathit{\alpha }}\mathit{J}(b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}},z_{\mathit{\beta }-\mathrm{1}/\mathrm{2}})h_{\mathit{\beta }-\mathrm{1}/\mathrm{2}},$$

where $h_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}$ is the distance between the layer centers, that is, $h_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}=(h_{\mathit{\beta }-\mathrm{1}}+h_{\mathit{\beta }})/\mathrm{2}$ for $\mathit{\beta }=\mathrm{2},\mathrm{\dots },N$ and $h_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}=h_{\mathrm{1}}/\mathrm{2}$ for β=1. The summation over the layers corresponds to vertical integration of the density Jacobian based on the piecewise constant profile of the density with the quadrature points placed at the interfaces. The density Jacobian at the interface is

$$\mathit{J}(b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}},z_{\mathit{\beta }-\mathrm{1}/\mathrm{2}})=\mathrm{\nabla }{b}_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}-D_z\left(b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}\right)\mathrm{\nabla }{z}_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}.$$

If $b_{\mathit{\beta }}=g\frac{{\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathit{\beta }}}{{\mathit{\rho }}_{\mathrm{0}}}$ is the layer buoyancy, the buoyancy at the interface is resolved with an average $b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}\left(b_{\mathit{\beta }-\mathrm{1}}+b_{\mathit{\beta }}\right)$ for $\mathit{\beta }=\mathrm{2},\mathrm{\dots }N$ and $b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}=b_{\mathrm{1}}$ for β=1. The approximation of the vertical derivative evaluated at the interface is resolved with finite differences. It is taken zero for the first interface $D_z\left(b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}\right)=\mathrm{0}$ for β=1 and $D_z\left(b_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}\right)=(b_{\mathit{\beta }-\mathrm{1}}-b_{\mathit{\beta }})/h_{\mathit{\beta }-\mathrm{1}/\mathrm{2}}$ for $\mathit{\beta }=\mathrm{2},\mathrm{\dots },N$. These choices allow us to recover a standard formula that can be found in Shchepetkin and McWilliams (2003) or in Klingbeil et al. (2018).

The tracer Eq. (10) admits a trivial solution, which we also want to inherit at the discrete level, the so-called tracer constancy condition: for a constant tracer, Eq. (10) reduces to the layerwise mass Eq. (14). The importance of preserving this property at a discrete level has been discussed extensively in Gross et al. (2002).

2.1 z-star

The multilayer model presented so far is based on a vertical subdivision of the fluid domain through the surface and terrain-following transformation Eq. (2), which leads to the coefficients l_α given in Eq. (4). Other vertical subdivisions can be used leading to different coefficients that, however, must both verify the positivity constraint and sum to one. In the following, we specify a slicing of the domain with both these properties based on a vertical coordinate transformation called z-star (Adcroft and Campin, 2004). The reference domain, with vertical coordinate Z, is

$${\mathrm{\Omega }}^{\mathrm{0}}=\mathit{\left\{}\right(\mathit{x},Z):\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}},-z_{\mathrm{b}}(\mathit{x})\le Z\le \mathrm{0}\mathit{\}}.$$

This domain is discretized by means of a vertical grid composed of N layers, with interfaces ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{0}}$, which are aligned to the geopotential. These interfaces can be described by constant functions:

$$Z_{\mathrm{1}/\mathrm{2}}=\mathrm{0}>Z_{\mathrm{2}-\mathrm{1}/\mathrm{2}}>\mathrm{\dots }>Z_{N+\mathrm{1}/\mathrm{2}}=-max{z}_{\mathrm{b}}\left(\mathit{x}\right).$$

As shown in Fig. 2, there is a substantial difference with respect to the vertical subdivision of the terrain-following grid. The grid interfaces could intersect the bathymetry and should be defined only in the fluid domain. We define the projection of the interface ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{0}}$ onto the horizontal plane as

$$\begin{array}{}\text{(16)}& {\mathrm{\Omega }}_{\mathit{x},\mathit{\alpha }}=\mathit{\{}\mathit{x}:\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x}}\mathrm{and}-z_{\mathrm{b}}(\mathit{x})\le Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\mathit{\}}.\end{array}$$

If a layer is bounded laterally by the bathymetry interface we can denote this lateral land boundary of the layer as

$${\mathrm{\Gamma }}_{\mathit{\alpha }}^b=\mathit{\left\{}\right(\mathit{x},Z):Z=-z_{\mathrm{b}}(\mathit{x})\mathrm{and}{Z}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}\le Z\le Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\mathit{\}}.$$

Each layer ${\mathrm{\Omega }}_{\mathit{\alpha }}^{\mathrm{0}}$ is delimited on the upper and bottom sides by ${\mathrm{\Gamma }}_{\mathit{\alpha }\mp \mathrm{1}/\mathrm{2}}^{\mathrm{0}}$ and laterally by the vertical domain boundary; it could also be delimited by ${\mathrm{\Gamma }}_{\mathit{\alpha }}^b$ (see Fig. 2, right panel). To map the reference interface ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{0}}$ to the physical interface ${\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$, again, we can use a generalized coordinate transformation, for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N$:

$$\begin{array}{}\text{(17)}& z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=\mathit{\zeta }(\mathit{x},t)+S_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\left(\mathit{x}\right)\left(\mathit{\zeta }(\mathit{x},t)+z_{\mathrm{b}}\left(\mathit{x}\right)\right),\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x},\mathit{\alpha }},\end{array}$$

with $S_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}$ a stretching function defined as

$$S_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\left(\mathit{x}\right)={\displaystyle \frac{Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}{z_{\mathrm{b}}\left(\mathit{x}\right)}}.$$

As in the previous Section, the layer thickness can be deduced from the total water depth. After some calculations we get

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{h}_{\mathit{\alpha }}(\mathit{x},t)& =z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}(\mathit{x},t)-max\left(z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}(\mathit{x},t),-z_{\mathrm{b}}\left(\mathit{x}\right)\right)\\ & =\left(S_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\left(\mathit{x}\right)-max\left(S_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}\left(\mathit{x}\right),-\mathrm{1}\right)\right)H(\mathit{x},t)\\ & =l_{\mathit{\alpha }}\left(\mathit{x}\right)H(\mathit{x},t),\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x},\mathit{\alpha }}\end{array}.\end{array}$$

If we define $\mathrm{\Delta }{Z}_{\mathit{\alpha }}\left(\mathit{x}\right)=Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}-max\left(Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}},-z_{\mathrm{b}}\left(\mathit{x}\right)\right)$ we can rewrite the coefficients, for $\mathit{\alpha }=\mathrm{1},N$:

$$l_{\mathit{\alpha }}\left(\mathit{x}\right)={\displaystyle \frac{\mathrm{\Delta }{Z}_{\mathit{\alpha }}\left(\mathit{x}\right)}{z_{\mathrm{b}}\left(\mathit{x}\right),}},\mathit{x}\in {\mathrm{\Omega }}_{\mathit{x},\mathit{\alpha }},$$

which is prescribed once the reference grid is created. The coefficients satisfy both the positivity constraint and locally they sum to one.

Figure 2 One-dimensional sketch of the reference (a) and physical (b) domains for the multilayer shallow water model with z-star layers.

An important property of the z-star transformation is that the horizontal domain Ω_x,α where the layer thickness h_α is defined, does not depend on time, as one can verify after computing the transformation Eq. (17) for $Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=-z_{\mathrm{b}}\left(\mathit{x}\right)$. This is particularly helpful because the number of layers does not depend on time, and the coefficients too. Other z-layers formulations based on similar mappings, such as the quasi-z layers (Mellor et al., 2002) or the hybrid $z/\mathit{\sigma }$ layers (Burchard and Petersen, 1997), do not share this property. For these coordinates a special treatment of the bottom is necessary: either an ad hoc modification of the bottom geometry or more interestingly these coordinates could be coupled with the porosity approach recently proposed by Debreu et al. (2020), where all the layers are present in the computation. For z-star the bottom momentum and tracer fluxes must be properly modified, replacing the maximum number of layers N, with the local number of layers $N_b\left(\mathit{x}\right)=\mathit{\{}\mathit{\alpha }:Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}<-z_{\mathrm{b}}(\mathit{x})\le Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\mathit{\}}$.

2.2 z-layers

For z-layers the actual interfaces do not depend on time and space:

$$z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}.$$

This method is implemented in the ocean models by allowing the top layer to vary in thickness without vanishing (Griffies et al., 2001). For the above transformation with fixed interfaces, the mass-transfer function (Eq. 14) coincides with the vertical velocity:

$$G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=-w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{-}=-w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{+},\mathit{\alpha }=\mathrm{2},N+\mathrm{1}.$$

Replacing the mass-transfer function with the vertical velocity in the multilayer model, we obtain the Eulerian model of Rambaud (2011).

3 Semi-implicit staggered finite element discretization

The discretization for both the z-star and the z-layers shallow water model can proceed in an equivalent fashion. We consider a discretization of the horizontal domain Ω_x∈ℝ² composed of non-overlapping triangular elements. We denote the horizontal grid by 𝒯 with K∈𝒯 the generic triangle, and $\left|K\right|$ its area. The local reference element length is h_K and it is computed as the minimum length of the triangle sides. With i∈𝒯 we denote the nodes of the grid. When no confusion is generated, we will locally number the nodes of the generic triangle as ($j=\mathrm{1},\mathrm{2},\mathrm{3}$ or j∈K). Given a node i in an element K, ${\mathit{n}}_i^K$ denotes the inward vector normal to the edge of K opposite to i, scaled by the length of the edge, see Fig. 3, left panel. For every node of the triangulation, 𝒟_i denotes the subset of triangles containing i. The dual cell C_i is obtained by joining the barycenters of the triangles in 𝒟_i with the midpoints of the edges meeting in i, as illustrated in Fig. 3, middle panel. Its area is

$$\left|C_i\right|=\sum _{K\in {\mathcal{D}}_i}{\displaystyle \frac{\left|K\right|}{\mathrm{3}}},$$

delimited by the boundary ∂C_i. The edge of ∂C_i separating C_i∩K and C_j∩K has an exterior normal called ${\mathit{n}}_{ij}^K$, as illustrated in Fig. 3, left panel. As before it is scaled by the edge length. Moreover, owing to the definition of the dual cell, we have

$$\begin{array}{}\text{(19)}& \sum _{j\in K,j\ne i}{\mathit{n}}_{ij}^K=-{\displaystyle \frac{{\mathit{n}}_i^K}{\mathrm{2}}}.\end{array}$$

After the horizontal discretization, the domain is subdivided into prismatic boxes $K\times [z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}},z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}]$. At the bottom, z-layers models apply a mask to non-existing land boxes that make the bathymetry stepped, as sketched in Fig. 3, right panel. The bottom layer for each element will be denoted as N_K. For a staggered discretization it is also helpful to define a nodal bottom layer $N_i={max}_{K\in {\mathcal{D}}_i}{N}_K$. The projections of the interfaces onto the horizontal plane are still denoted as Ω_x,α and defined with Eq. (16), this time evaluated with the stepwise approximation of the bathymetry. Then, a layer dual cell C_α,i can be defined by considering 𝒟_α,i the subset of elements sharing node i and in Ω_x,α. Its area is

$$\left|C_{\mathit{\alpha },i}\right|=\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}{\displaystyle \frac{\left|K\right|}{\mathrm{3}}}.$$

Figure 3 Grid and notation. (a) Triangle K with nodes and scaled normals. (b) Set 𝒟_i with dual cell area C_i and dual cell boundary ∂C_i. The degrees of freedom are also shown: discharge ▪, tracer, and free surface ◦. (c) Stepped bathymetry with masked boxes in brown, after the horizontal discretization.

On a B-staggered grid the free surface elevation, the discharges, and the tracers are described with basis functions of different order and support (Williams and Zienkiewicz, 1981). The discharge field and the tracer field belong to a finite dimensional space with basis composed by the piecewise constant functions. For the discharges, the space has a basis {ψ_K}_K∈𝒯 composed of the characteristic functions on the triangle, whereas for the tracers we choose {ϕ_i}_i∈𝒯 composed of the characteristic functions on the dual cell. The discharge fields q_α=h_αu_α and the tracers T_α are approximated through (we use an abuse of notation employing the same symbol of the continuous variable):

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathit{q}}_{\mathit{\alpha }}(\mathit{x},t)=\sum _{K\in \mathcal{T}}{\mathit{\psi }}_K\left(\mathit{x}\right){\mathit{q}}_{\mathit{\alpha },K}\left(t\right),\text{(21)}& {\displaystyle }{T}_{\mathit{\alpha }}(\mathit{x},t)=\sum _{i\in \mathcal{T}}{\mathit{\varphi }}_i\left(\mathit{x}\right)T_{\mathit{\alpha },i}\left(t\right),\end{array}$$

with q_α,K(t), defined for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N_K$, being the elemental discharge values and with T_α,i(t), defined for $\mathit{\alpha }=\mathrm{1},\mathrm{\dots },N_i$, the nodal tracer values. The free surface belongs to a space of finite dimension, with a basis {φ_i}_i∈𝒯, which denotes the standard continuous piecewise linear Lagrangian basis. The discrete free surface is given by

$$\begin{array}{}\text{(22)}& \mathit{\zeta }(\mathit{x},t)=\sum _{i\in \mathcal{T}}{\mathit{\phi }}_i\left(\mathit{x}\right){\mathit{\zeta }}_i\left(t\right),\end{array}$$

where ζ_i(t) are the nodal free surface values. Note that the discrete discharges and discrete tracers are discontinuous respectively across the boundaries of the triangles and of the dual cells, whereas the discrete free surface is globally continuous. On a B-grid the layer thickness is naturally computed at the grid nodes h_α,i, where the free surface is available. The element values h_α,K are a conservative average of the nodal values. The element velocities are obtained from ${\mathit{u}}_{\mathit{\alpha },K}=\frac{{\mathit{q}}_{\mathit{\alpha },K}}{h_{\mathit{\alpha },K}}$.

We obtain the weak formulation multiplying mass and momentum Eqs. (8) and (9) by the test functions that belong to the same space of the solution and integrating it on the horizontal domain. The finite element discretization reduces to compute the integrals accounting for the different terms. For the mass flux term, which is integrated by parts we need to compute with a proper quadrature rule the following integral (only x-component shown):

$$a_{iK}^x=\underset{K}{\int }{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}\mathrm{d}\mathit{x}.$$

The boundary term has been neglected as it cancels out except at the lateral domain boundary. Similarly, for the terms that will be treated explicitly in the momentum equation namely the horizontal and vertical advection and the internal pressure gradient, we have

$$f_{\mathit{\alpha },K}^x=-\underset{\partial K}{\int }\widehat{{\mathit{q}}_{\mathit{\alpha }}{u}_{\mathit{\alpha }}}\cdot \mathit{n}\mathrm{d}s+\underset{K}{\int }\left(B_{\mathit{\alpha }}^x+[uG{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\right)\mathrm{d}\mathit{x}.$$

The horizontal advection term is resolved with a first-order upwind flux $\widehat{{\mathit{q}}_{\mathit{\alpha }}{u}_{\mathit{\alpha }}}$ (Umgiesser et al., 2004). To write the scheme in matrix form, exploiting the compactness of the staggered discretization, we introduce “vertical” vectors and matrices, that pile up all the layers for a single element K, and we denote them with bold capital letters. For example, the layer discharges and the layer thickness are regrouped in the following vectors:

$$\begin{array}{rl}{\mathbf{U}}_K& =\left(\begin{array}{c}{q}_{\mathrm{1},K}^x\\ \mathrm{\dots }\\ q_{\mathit{\alpha },K}^x\\ \mathrm{\dots }\\ q_{N_K,K}^x\end{array}\right),{\mathbf{V}}_K=\left(\begin{array}{c}{q}_{\mathrm{1},K}^y\\ \mathrm{\dots }\\ q_{\mathit{\alpha },K}^y\\ \mathrm{\dots }\\ q_{N_K,K}^y\end{array}\right),\\ & {\mathbf{H}}_K=\left(\begin{array}{c}{h}_{\mathrm{1},K}\\ \mathrm{\dots }\\ h_{\mathit{\alpha },K}\\ \mathrm{\dots }\\ h_{N_K,K}\end{array}\right),\end{array}$$

and analogously the explicit terms:

$${\mathit{F}}_K^x=\left(\begin{array}{c}{f}_{\mathrm{1},K}^x\\ \mathrm{\dots }\\ f_{\mathit{\alpha },K}^x\\ \mathrm{\dots }\\ f_{N_K,K}^x\end{array}\right),{\mathit{F}}_K^y=\left(\begin{array}{c}{f}_{\mathrm{1},K}^y\\ \mathrm{\dots }\\ f_{\mathit{\alpha },K}^y\\ \mathrm{\dots }\\ f_{N_K,K}^y\end{array}\right).$$

The vertical viscous term is recast in matrix form via a tridiagonal matrix ${\mathbf{A}}_K^d\in {\mathbb{R}}^{N_K\times N_K}$. The bottom momentum flux has to be integrated into this matrix. Note that all these vectors and matrices are restricted to nonmasked boxes.

Following Casulli and Cattani (1994), we build a semi-implicit time discretization by treating semi-implicitly the mass flux and the free surface gradient in the momentum equation. The vertical viscous term can also cause a restrictive time-step and is handled here implicitly without major computation issues but allowing the CFL condition to be relaxed. We define the variation of a quantity in a time step as $\mathrm{\Delta }u=u^{n+\mathrm{1}}-u^n$, then:

$$u^{n+\mathit{\theta }}=\mathit{\theta }{u}^{n+\mathrm{1}}+(\mathrm{1}-\mathit{\theta })u^n=\mathit{\theta }\mathrm{\Delta }u+u^n.$$

After applying the previous definition to the semi-discrete equations, the semi-implicit momentum equations on an unstructured B-grid read

$$\begin{array}{}\text{(23)}& {\displaystyle }\mathrm{\Delta }{\mathbf{U}}_K=\mathrm{\Delta }{\mathbf{U}}_K^{*}-\mathrm{\Delta }tg{\mathbf{A}}_K^{-\mathrm{1}}{\mathbf{H}}_K^n\sum _{j\in K}{a}_{jK}^x\mathit{\theta }\mathrm{\Delta }{\mathit{\zeta }}_j,\text{(24)}& {\displaystyle }\mathrm{\Delta }{\mathbf{V}}_K=\mathrm{\Delta }{\mathbf{V}}_K^{*}-\mathrm{\Delta }tg{\mathbf{A}}_K^{-\mathrm{1}}{\mathbf{H}}_K^n\sum _{j\in K}{a}_{jK}^y\mathit{\theta }\mathrm{\Delta }{\mathit{\zeta }}_j,\end{array}$$

with ${\mathbf{A}}_K=\left(\mathit{I}\left|K\right|-\mathrm{\Delta }t{\mathbf{A}}_K^d\right)$ a tridiagonal, positive definite, and diagonally dominant matrix. The nonlinear dependence of the external pressure gradient term from H_K has been resolved by using the old value. Also, the viscous matrix has been computed with frozen values at tⁿ. In ${\mathit{F}}_K^n$ all the quantities are computed at tⁿ, including the mass-transfer function. These choices avoid solving a nonlinear system at each time step. The variation $\mathrm{\Delta }(\cdot {)}^{*}=(\cdot {)}^{*}-(\cdot {)}^n$ is the solution of the following Euler step with an explicit external pressure gradient:

$$\begin{array}{}\text{(25)}& {\displaystyle }\mathrm{\Delta }{\mathbf{U}}_K^{*}=\mathrm{\Delta }t{\mathbf{A}}_K^{-\mathrm{1}}({\mathit{F}}_K^{x,n}+{\mathbf{A}}_K^d{\mathbf{U}}_K^n-g{\mathbf{H}}_K^n\sum _{j\in K}{a}_{jK}^x{\mathit{\zeta }}_j^n),\text{(26)}& {\displaystyle }\mathrm{\Delta }{\mathbf{V}}_K^{*}=\mathrm{\Delta }t{\mathbf{A}}_K^{-\mathrm{1}}({\mathit{F}}_K^{y,n}+{\mathbf{A}}_K^d{\mathbf{V}}_K^n-g{\mathbf{H}}_K^n\sum _{j\in K}{a}_{jK}^y{\mathit{\zeta }}_j^n).\end{array}$$

If the expressions for ΔU_K and ΔV_K, Eqs. (23) and (24), are introduced into the discrete mass equation, we obtain a linear system with only the free surface coefficients as unknowns:

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\sum _{K\in {\mathcal{D}}_i}\sum _{j\in K}& \left(m_{ij}^K+g{\mathit{\theta }}^{\mathrm{2}}\mathrm{\Delta }{t}^{\mathrm{2}}\left(a_{iK}^x{\mathbf{1}}^{\mathrm{T}}{\mathbf{A}}_K^{-\mathrm{1}}{\mathbf{H}}_K^n{a}_{jK}^x\right.\right.\\ & \left.\left.+a_{iK}^y{\mathbf{1}}^{\mathrm{T}}{\mathbf{A}}_K^{-\mathrm{1}}{\mathbf{H}}_K^n{a}_{jK}^y\right)\right)\mathrm{\Delta }{\mathit{\zeta }}_j\\ & =\mathrm{\Delta }t\sum _{K\in {\mathcal{D}}_i}\left(a_{iK}^x{\mathbf{1}}^{\mathrm{T}}\left(\mathit{\theta }\mathrm{\Delta }{\mathbf{U}}_K^{*}+{\mathbf{U}}_K^n\right)\right.\\ & \left.+a_{iK}^y{\mathbf{1}}^{\mathrm{T}}\left(\mathit{\theta }\mathrm{\Delta }{\mathbf{V}}_K^{*}+{\mathbf{V}}_K^n\right)\right)\end{array},\end{array}$$

where $m_{ij}^K={\int }_K{\mathit{\phi }}_i{\mathit{\phi }}_j\mathrm{d}\mathit{x}$ is the Galerkin mass matrix based on the piecewise linear Lagrange basis functions. The Galerkin mass matrix, in SHYFEM, is lumped. The vector $\mathbf{1}\in {\mathbb{R}}^{N_K}$ has all components being one.

The hydrodynamic time step flow chart is thus the following: we first perform the Euler step Eqs. (25) and (26). Then we resolve the mass Eq. (27) and we complete the momentum update with the semi-implicit step Eqs. (23) and (24). Finally, we compute the layer thickness at the grid nodes. For z-star, we use the expression Eq. (18) at the grid nodes. For the z-layers, the layer thickness does not change except for the first layer.

3.1 Mass-transfer function

After the hydrodynamic update of the previous paragraph, the discrete mass-transfer function is computed. We employ the same continuous piecewise linear approximation used for the free surface. The nodal values are computed from a finite-element mass-lumped discretization of the layerwise mass Eq. (14). As for the depth-integrated mass equation, the discharge is evaluated semi-implicitly. Starting from the bottom with $G_{N_i+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}=\mathrm{0}$, for $\mathit{\alpha }=N_i,\mathrm{\dots },\mathrm{1}$:

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\left|C_{\mathit{\alpha },i}\right|G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}& =\left|C_{\mathit{\alpha }+\mathrm{1},i}\right|G_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}+\left|C_{\mathit{\alpha },i}\right|{\displaystyle \frac{\mathrm{\Delta }{h}_{\mathit{\alpha },i}}{\mathrm{\Delta }t}}\\ & -\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}(a_{iK}^x{q}_{\mathit{\alpha },K}^{x,n+\mathit{\theta }}+a_{iK}^y{q}_{\mathit{\alpha },K}^{y,n+\mathit{\theta }})\end{array}.\end{array}$$

Note that the semi-implicit discretization ensures vertical mass-conservation. Summing up Eq. (28) for all the layers and using Eq. (27) with a lumped Galerkin mass-matrix to cancel the right-hand side, we get the impermeability condition at the free surface $G_{\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}=\mathrm{0}$. With standard z-layers, the contribution related to the grid velocity is zero $\mathrm{\Delta }{h}_{\mathit{\alpha },i}=\mathrm{\Delta }t[{\mathit{\sigma }}_i{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=\mathrm{0}$, except for the first layer.

3.2 Tracers

The semi-implicit update is completed with the time-stepping of the tracer. Vertical diffusion is treated implicitly and the remaining advection terms are explicit. The spatial discretization of the explicit terms implies the computation of the following integrals, which account for the horizontal and vertical advection terms:

$$f_{\mathit{\alpha },i}=-\underset{\partial C_{\mathit{\alpha },i}}{\int }\widehat{T_{\mathit{\alpha }}{\mathit{q}}_{\mathit{\alpha }}}\cdot \mathit{n}\mathrm{d}s+\underset{C_{\mathit{\alpha },i}}{\int }[TG{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\mathrm{d}\mathit{x},$$

where $\widehat{T_{\mathit{\alpha }}{\mathit{q}}_{\mathit{\alpha }}}$ is an appropriate numerical tracer flux across the dual cell boundary. At the lateral boundary $\partial {\mathrm{\Omega }}_{\mathit{x},\mathit{\alpha }}$, the tracer flux is zero for land boundaries, whereas it is determined by the boundary conditions at the domain boundary. In the discussion that follows we consider only nodes that do not lie on the domain boundary. On a triangular grid the two terms read

$$\begin{array}{}\text{(29)}& {\displaystyle }\begin{array}{rl}\underset{\partial C_{\mathit{\alpha },i}}{\int }& \widehat{T_{\mathit{\alpha }}{\mathit{q}}_{\mathit{\alpha }}}\cdot \mathit{n}\mathrm{d}s=\sum _{K\in {\mathcal{D}}_{\mathit{\alpha },i}}\sum _{j\in K,j\ne i}\widehat{T_{\mathit{\alpha }}{\mathit{q}}_{\mathit{\alpha }}}\cdot {\mathit{n}}_{ij}^K\\ & =\sum _{K\in {\mathcal{D}}_{\mathit{\alpha },i}}\sum _{j\in K,j\ne i}{\widehat{H}}_{\mathit{\alpha }}(T_{\mathit{\alpha },i},T_{\mathit{\alpha },j})\end{array},\text{(30)}& {\displaystyle }\begin{array}{rl}\underset{C_{\mathit{\alpha },i}}{\int }& [TG{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\mathrm{d}\mathit{x}=|C_{\mathit{\alpha },i}|T_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}{G}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}\\ & -\left|C_{\mathit{\alpha }+\mathrm{1},i}\right|T_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}{G}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}\end{array},\end{array}$$

with ${\widehat{H}}_{\mathit{\alpha }}(T_{\mathit{\alpha },i},T_{\mathit{\alpha },j})$ being the numerical flux in the horizontal direction and $T_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}{G}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}$ the numerical flux in the vertical direction. The SHYFEM model implements second-order consistent TVD fluxes in both directions.

Using the notation with bold capital letters denoting “vertical” vectors, the tracer values and the explicit term at the nodes are regrouped in the following:

$${\mathit{T}}_i=\left(\begin{array}{c}{T}_{\mathrm{1},i}\\ \mathrm{\dots }\\ T_{\mathit{\alpha },i}\\ \mathrm{\dots }\\ T_{N_i,i}\end{array}\right),{\mathit{F}}_i=\left(\begin{array}{c}{f}_{\mathrm{1},i}\\ \mathrm{\dots }\\ f_{\mathit{\alpha },i}\\ \mathrm{\dots }\\ f_{N_i,i},\end{array}\right).$$

Vertical diffusion can also be assembled in matrix form through the discrete matrix ${\mathbf{A}}_i^d\in {\mathbb{R}}^{N_i\times N_i}$. Then, the discretization of the layerwise tracer Eq. (10) read

$$\begin{array}{}\text{(31)}& {\mathbf{A}}_i{\mathit{T}}_i^{n+\mathrm{1}}=\mathrm{Diag}\mathit{\left\{}\right|C_{\mathit{\alpha },i}\left|h_{\mathit{\alpha },i}^n\mathit{\right\}}{\mathit{T}}_i^n+\mathrm{\Delta }t{\mathit{F}}_i^n,\end{array}$$

with ${\mathbf{A}}_i=\left(\mathrm{Diag}\mathit{\left\{}\right|C_{\mathit{\alpha },i}\left|h_{\mathit{\alpha },i}^{n+\mathrm{1}}\mathit{\right\}}-\mathrm{\Delta }t{\mathbf{A}}_i^d\right)$ the vertical tracer matrix. Although the advection terms are explicit, it should be noted that the horizontal numerical fluxes are computed with the discharges evaluated at ${\mathit{q}}_{\mathit{\alpha }}^{n+\mathit{\theta }}$ whereas the vertical numerical flux uses the last available mass-transfer function $G_{\mathit{\alpha }\pm \mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}$ from Eq. (28). This choice is important to maintain the consistency of the discrete tracer equation with the layerwise mass equation. In fact inserting a constant tracer in Eq. (31), yields exactly the discrete layerwise mass Eq. (28). The proof is left in the Appendix.

To conclude, we summarize the whole time step flow chart: after the hydrodynamic update described in Sect. 3, we compute the mass-transfer function Eq. (28) and, last, we update the tracers with Eq. (31).

4 z-surface-adaptive layers

In this section, we enhance the z-layers shallow water model by introducing a new algorithm that allows for the dynamic insertion and removal of surface boxes or, with an abuse of language, of surface layers. To differentiate it from the standard z-layers, we refer to this enhanced version as z-surface-adaptive layers. The key idea is to interpret the area swept by the layer interface in the time step $\mathrm{\Delta }t\in \left[t^n,t^{n+\mathrm{1}}\right)$ as the sum of two contributions: one owing to the mesh movement driven by the free surface oscillation (grid movement) and one owing to the collapse/expansion of the layer (topology change). These topology changes, in fact, can be seen as a continuous deformation of the layer interfaces performed within the time step. With this in mind, the final position of the interfaces at the grid nodes ${\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}={\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}({\mathit{x}}_i,t^{n+\mathrm{1}})$ is

$${\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}=z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}+\mathrm{\Delta }{\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i},$$

where $z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}=z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}({\mathit{x}}_i,t^{n+\mathrm{1}})$ is the interface position after the grid movement and $\mathrm{\Delta }{\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}$ is the contribution of the interface collapse/expansion, basically a correction term. Similarly, the grid velocity in the time step can be decomposed as

$${\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}={\displaystyle \frac{{\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^n}{\mathrm{\Delta }t}}={\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{mov}}+{\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{top}},$$

with

$${\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{mov}}={\displaystyle \frac{z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^n}{\mathrm{\Delta }t}},{\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{top}}={\displaystyle \frac{\mathrm{\Delta }{\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}}{\mathrm{\Delta }t}}.$$

In the solution of the multilayer shallow water equations we employ a splitting procedure, where the two aforementioned contributions are treated in two steps. In the first step (grid movement) we solve the multilayer model on a vertical grid where the surface layers adjust locally in order to maintain a positive thickness. In the subsequent step, we locally remove surface fluid boxes with a minimal thickness or split fluid boxes that are excessively thick. The evolution of the vertical grid and of the tracer in one time step is shown in Fig. 4. The top row shows the case of a layer removal and the bottom row the case of a layer insertion. As a remark, we stress that the above interpretation of the interface displacement reveals many beneficial aspects with respect to the direct insertion and removal of a layer. Without the grid movement step, it would be more complicated to time step the tracers on a grid with positive layer thickness, with all the related stability issues. In fact in the tracer update Eq. (24) the layer thickness at tⁿ⁺¹ is needed. One may think to compute the tracer after the insertion/removal operations have been performed (thus having positive layer thickness both at tⁿ and tⁿ⁺¹), but then the configuration on which the discrete tracer equation is solved is ambiguous and it seems hard to ensure the consistency with the continuity or to verify the tracer constancy property.

Figure 4 Grid and tracer evolution during one time step. The process is interpreted as four stages, which bring from the pair $(T^n,{\mathit{\zeta }}_{\mathrm{h}}^n)$ to $({\stackrel{\mathrm{̃}}{T}}^{n+\mathrm{1}},{\mathit{\zeta }}^{n+\mathrm{1}})$. The vector $T=\mathit{\{}{T}_{\mathrm{1}},T_{\mathrm{2}}\mathit{\}}$ collects the layer values of the tracer. The dashed line indicates the removed interface. Top: case of surface layer removal. Bottom: case of surface layer insertion.

In the following we provide the technical details to realize such adaptation to the free surface with the z-layers. First, we notice that, since the beginning of the simulation, the index of the surface layer may change spatially at the element boundaries. Given the initial free surface elevation ζ⁰(x), we define a set of active indices and the surface layer index, by element, as

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{active},K}& =\mathit{\{}\mathit{\alpha }\in {\mathit{\alpha }}_K:Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}+{\mathit{ϵ}}_{\mathrm{top}}<\underset{\mathit{x}\in K}{min}{\mathit{\zeta }}^{\mathrm{0}}(\mathit{x}\left)\mathit{\right\}},\\ & {\mathit{\alpha }}_{\mathrm{top},K}=min{\mathit{\alpha }}_{\mathrm{active},K},\end{array}\end{array}$$

with ${\mathit{\alpha }}_K=\mathit{\{}\mathrm{1},\mathrm{\dots },N_K\mathit{\}}$. Owing to the staggering of the grid, it is convenient to define also at each node:

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{active},i}& =\mathit{\{}\mathit{\alpha }\in {\mathit{\alpha }}_i:Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}+{\mathit{ϵ}}_{\mathrm{top}}<{\mathit{\zeta }}_i^{\mathrm{0}}\mathit{\}},\\ & {\mathit{\alpha }}_{\mathrm{top},i}=min{\mathit{\alpha }}_{\mathrm{active},i},\end{array}\end{array}$$

with ${\mathit{\alpha }}_i=\mathit{\{}\mathrm{1},\mathrm{\dots },N_i\mathit{\}}$. The parameter ϵ_top is a small positive constant that fixes the minimum allowable depth for a surface layer to exist. Below this threshold, the layer is removed. We have fixed it as ϵ_top=0.2ΔZ_α. It turns out that this parameter is quite important as it avoids the presence of very small layers. Such layers can lead to a restrictive time step owing to the explicit discretization of vertical advection terms. In Fig. 5 we illustrate the spatial variation of the top layer index for a one-dimensional example.

Figure 5 This one-dimensional example shows the grid for the z-surface-adaptive layers. Elemental surface layer indices are shown at the bottom, nodal surface layer indices are shown at the top.

4.1 Vertical grid movement

We evolve the discrete multilayer shallow water equations using the semi-implicit finite element method detailed in Sect. 3. The vertical vectors and matrices are restricted to the layers with an active index. Moreover, to account for the movement of the surface layers, the layer thickness is updated as follows:

• We identify the indices associated with the layers that, locally, undergo a deformation. They are defined as the layers of the reference grid whose top interface is above the free surface or by the set of indices:

  $$\begin{array}{}\text{(34)}& {\mathit{\alpha }}_{\mathrm{mov},i}=\mathit{\{}\mathit{\alpha }\in {\mathit{\alpha }}_i:Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}+{\mathit{ϵ}}_{\mathrm{mov}}>{\mathit{\zeta }}_i^{n+\mathrm{1}}\mathit{\}}.\end{array}$$

  ϵ_mov is a small and positive constant that we have added. Below this threshold, the vertical grid movement is deployed. As seen for ϵ_top, it avoids the presence of very small layers that can be dangerous from a numerical point of view. The bottom-most layer is denoted by $N_{\mathrm{mov},i}=max{\mathit{\alpha }}_{\mathrm{mov},i}$. The depth of the moving layers is

  $$z_{\mathrm{mov},i}=max\left(Z_{N_{\mathrm{mov},i}+\mathrm{1}/\mathrm{2}},-z_{b,i}\right).$$

• We compute the new layer thickness after a local grid deformation that absorbs the free surface movement. To move the interfaces of the layers contained in the set, we use the generalized coordinates transformation Eq. (1), which takes the form:

  $$\begin{array}{}\text{(35)}& z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}={\mathit{\zeta }}_i^{n+\mathrm{1}}+S_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}\left({\mathit{\zeta }}_i^{n+\mathrm{1}}+z_{\mathrm{mov},i}\right),\end{array}$$

  with $S_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}$ a stretching function. Then, the nodal layer thickness reads:

  $$\begin{array}{}\text{(36)}& h_{\mathit{\alpha },i}^{n+\mathrm{1}}=l_{\mathit{\alpha },i}\left({\mathit{\zeta }}_i^{n+\mathrm{1}}+z_{\mathrm{mov},i}\right),\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i},\mathrm{\dots },N_{\mathrm{mov},i}.\end{array}$$

  For the proportionality coefficients, we have tried different definitions allowing a smooth movement on the interfaces between the time steps, without experiencing any major impact on the results. For simplicity we have thus implemented a z-star definition $l_{\mathit{\alpha },i}=\frac{\mathrm{\Delta }{Z}_{\mathit{\alpha }}}{z_{\mathrm{mov},i}}$ (see Sect. 2).

This is shown in Fig. 4, first and second columns. The new layer thickness is used in the update of the tracers, Eq. (31). We stress the fact that the vertical configuration is taken constant, i.e., the number of layers at each element remains constant during the time stepping of the discharges and of the tracers.

4.2 Removal/insertion of surface layers

Then we perform the insertion/removal operation based on:

• An update of the active layers and of the top layer index by re-evaluating Eqs. (32) and (33) with the new free surface elevation ζⁿ⁺¹. We get the new top layer indices ${\mathit{\alpha }}_{\mathrm{top},K}^{n+\mathrm{1}}$ and ${\mathit{\alpha }}_{\mathrm{top},i}^{n+\mathrm{1}}$.

• Once we have identified the index that should be inserted/removed in the active set, we proceed with the collapse/expansion of the surface boxes. A conservative remapping step is necessary to pass the unknowns from the old vertical grid to the new one.

We use the tilde ${\stackrel{\mathrm{̃}}{T}}_{\mathit{\alpha }}^{n+\mathrm{1}}$ to distinguish a generic layer variable (the tracer in this case) remapped onto the new grid from the solution time stepped on the old grid $T_{\mathit{\alpha }}^{n+\mathrm{1}}$. The remapped value is the solution of the following advection equation integrated over the layer thickness:

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{h}}_{\mathit{\alpha }}{\stackrel{\mathrm{̃}}{T}}_{\mathit{\alpha }}}{\partial t}}=[{\mathit{\sigma }}^{\mathrm{top}}\stackrel{\mathrm{̃}}{T}{]}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}^{\mathit{\alpha }-\mathrm{1}/\mathrm{2}},\end{array}$$

with an upwind flux:

$$\begin{array}{}\text{(38)}& {\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{top}}{T}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}={\left({\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{top}}\right)}^{+}{T}_{\mathit{\alpha }}^{n+\mathrm{1}}+{\left({\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}^{\mathrm{top}}\right)}^{-}{T}_{\mathit{\alpha }-\mathrm{1}}^{n+\mathrm{1}}.\end{array}$$

We consider the discrete case. After integration on the dual cell and with a simple forward Euler time stepping (with initial condition $T_{\mathit{\alpha }}^{n+\mathrm{1}}$) we have

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{h}}_{\mathit{\alpha },i}^{n+\mathrm{1}}& {\stackrel{\mathrm{̃}}{T}}_{\mathit{\alpha },i}^{n+\mathrm{1}}=h_{\mathit{\alpha },i}^{n+\mathrm{1}}{T}_{\mathit{\alpha },i}^{n+\mathrm{1}}\\ & +\mathrm{\Delta }t\left({\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{top}}{T}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-{\mathit{\sigma }}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{\mathrm{top}}{T}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}\right)\end{array},\end{array}$$

with the new nodal layer thickness:

$${\stackrel{\mathrm{̃}}{h}}_{\mathit{\alpha },i}^{n+\mathrm{1}}={\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-{\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}.$$

In the case of an element removal (${\mathit{\alpha }}_{\mathrm{top},i}^{n+\mathrm{1}}>{\mathit{\alpha }}_{\mathrm{top},i}^n$), we identify the layer that should disappear and we proceed with a collapse of the lower interface to the upper one. For $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i}^n,\mathrm{\dots },{\mathit{\alpha }}_{\mathrm{top},i}^{n+\mathrm{1}}$, the discrete remap Eq. (39) with Eq. (38) reduces trivially to transfer the depth-integrated tracer that belongs to the removed layers to the upper active layer. In the case of an element insertion (${\mathit{\alpha }}_{\mathrm{top},i}^{n+\mathrm{1}}<{\mathit{\alpha }}_{\mathrm{top},i}^n$), we identify the layer that should appear and we expand the interface. Then the remap for $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i}^{n+\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_{\mathrm{top},i}^n$ reduces to distribute the depth-integrated variable across the existing and inserted layers with a weighted average. This is shown in Fig. 4, third and fourth columns. All the unknowns must be remapped. The discharges that are defined on the elements Eq. (37) should be integrated on the element. This completes the time step.

4.3 Connection to z-star

We have a small parameter ϵ_mov to fix. It is convenient to express this constant as a percentage of the reference layer thickness, ϵ_mov=r_movΔZ_α in the relationship Eq. (34). In order to obtain the z-surface-adaptive grid we have chosen r_mov≤r_top; in practice, we have set r_mov=0.15. The grid deformation is localized to the free surface. As long as the surface fluid boxes are deformed, they are recognized as too small and immediately removed in the grid topology step. This implies working, at the next time step, with z-layers having all the interfaces aligned to the geopotentials.

Interestingly, we can obtain other grids by increasing r_mov. We define

$$\begin{array}{}\text{(40)}& R_{\mathit{\alpha }}={\displaystyle \frac{{\mathit{\zeta }}_{\max }-Z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }{Z}_{\mathit{\alpha }}}},\end{array}$$

with ${\mathit{\zeta }}_{\max }\ge \underset{\mathit{x},t}{max}\mathit{\zeta }(\mathit{x},t)$ an estimate of the maximum free surface height during the simulation. We get

• z-star if r_mov≥R_N and no insertion/removal. The whole water column is subjected to the grid movement, whereas the number of layers does not change.

• z-star + z if r_mov=R_M and no insertion/removal. The upper part of the water column, at minimum M layers, is subjected to the grid movement, whereas the lower part is fixed.

Figure 6 shows a sketch of the different possibilities. Tuning properly r_mov we compare the newly developed z-surface adaptive layers against z-star.

Figure 6 The different vertical grids outlined in Sect. 4.3.

5 Advection with a spatially variable number of layers

We have used an approach where the grid topology does not change during the time step of the conserved variables, i.e., the numerical scheme of Sect. 3 works on the deforming grid of Sect. 4.1, with a temporally constant number of layers between tⁿ and tⁿ⁺¹. However, in the previous time step, a layer insertion/removal may occur (to remove very thin surface layers or to split a thicker layer) on a certain element and not on its neighbors. This results in a vertical discretization with a spatially variable number of layers (see Fig. 7), which slightly complicates the treatment of advection terms (see on this topic Bonaventura et al. (2018)).

Figure 7 Nonconformal box for the one-dimensional case. The nonconformal box is in gray. Discharges, layer thickness, and tracers are shown.

Figure 8 Treatment of nonconformal box for the one-dimensional case. (a) Splitting with fictitious layers. (b) The mass-transfer function $G_{\mathrm{1}+\mathrm{1}/\mathrm{2},i}$ at hanging point is represented by a red arrow.

Consider the one dimensional example in Fig. 7, where two contiguous elements with a different top-layer index ${\mathit{\alpha }}_{\mathrm{top},i+\mathrm{1}/\mathrm{2}}>{\mathit{\alpha }}_{\mathrm{top},i-\mathrm{1}/\mathrm{2}}$ exist. In correspondence with node i a change of the element top layer index takes place. Borrowing the vocabulary from the literature on nonconformal meshes, we have a vertical edge with a hanging point. We call it a hanging layer, a layer for which at least one interface ends with a hanging point. The boxes that have vertical edges across which the element top-layer index varies deserve special treatment. In our case, with only insertion/removal of surface layers, we can easily flag such boxes by checking, for each element, that the nodal top layer index is different from the elemental one. The elements of the grid with a nonconformal surface box are indicated by an asterisk:

$$\mathrm{if}{\mathit{\alpha }}_{\min ,K}<{\mathit{\alpha }}_{\mathrm{top},K}\text{then}K=K^{*},$$

with ${\mathit{\alpha }}_{\min ,K}={min}_{j\in K}{\mathit{\alpha }}_{\mathrm{top},j}$. Then the boxes called hereinafter for simplicity “nonconformal” can be identified by the pair of indices $\left({\mathit{\alpha }}_{\mathrm{top},K},K^{*}\right)$. As both mass and tracer fluxes need communication with the neighbors' boxes, they have to be treated differently. Moreover, for the tracer discrete update, we have to take care of preserving the constancy property.

In the case of a nonconformal box we proceed as follows. We split the box vertically in fictitious layers through planar interfaces passing through the hanging points of nonconformal edges and some fraction of the conformal edge length (see Fig. 8, left panel). From this geometrical configuration we compute the element layer thickness $h_{\mathit{\alpha },K}^{*}$ for the fictitious layers. Then we distribute the discharge of the top layer among the fictitious layers, for $\mathit{\alpha }={\mathit{\alpha }}_{\min ,K},\mathrm{\dots },{\mathit{\alpha }}_{\mathrm{top},K}$:

$$\begin{array}{}\text{(41)}& {\mathit{q}}_{\mathit{\alpha },K}^{*}=l_{\mathit{\alpha },K}^{*}{\mathit{q}}_{{\mathit{\alpha }}_{\mathrm{top},K},K},\end{array}$$

with $l_{\mathit{\alpha },K}^{*}=\frac{h_{\mathit{\alpha },K}^{*}}{h_{{\mathit{\alpha }}_{\mathrm{top},K}K}}$. These values are used to complete both mass and tracer fluxes for the missing layers of nonconformal boxes. We consider the case of a nonconformal box $({\mathit{\alpha }}_{\mathrm{top},K},K)$ with node i∈K, as illustrated in one dimension in Fig. 7. After the splitting Eq. (41), the mass-flux term (only the x-component shown) reads, for $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i},\mathrm{\dots },{\mathit{\alpha }}_{\mathrm{top},K}$:

$$\begin{array}{}\text{(42)}& \underset{K}{\int }{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}{q}_{\mathit{\alpha }}^{*}\mathrm{d}\mathit{x}=a_{iK}{c}_{\mathit{\alpha },i}^{*}{q}_{{\mathit{\alpha }}_{\mathrm{top},K},K},\end{array}$$

with

$$\begin{array}{}\text{(43)}& c_{\mathit{\alpha },i}^{*}=\left\{\begin{array}{ll}\sum _{\mathit{\beta }={\mathit{\alpha }}_{\mathrm{top},i}}^{{\mathit{\alpha }}_{\min ,K}}{l}_{\mathit{\beta },K}^{*}& \mathrm{if}\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i}\mathrm{and}{\mathit{\alpha }}_{\min ,K}<{\mathit{\alpha }}_{\mathrm{top},i}\\ l_{\mathit{\alpha },K}^{*}& \text{otherwise (hanging layer)}\end{array}\right.,\end{array}$$

where the two cases account for the contribution of element K to nodes with or without hanging layers, node i and i+1 respectively in Fig. 8. Such a contribution from the nonconformal box is added to the mass-flux term in the layerwise mass equation. It allows the mass-transfer function to be computed at the hanging points $G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}$ for $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i},\mathrm{\dots }{\mathit{\alpha }}_{\mathrm{top},K}$, as shown in Fig. 8, right panel. One can check that this treatment is mass-conserving. Summing the mass-transfer function for all the layers, even in the presence of nonconformal boxes, still yields the discrete mass Eq. (27).

The horizontal advection scheme (Eq. 29) on the nonconformal box can be applied straightforwardly to the fictitious layers. Then, the numerical flux in nonconformal boxes reads for $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i},\mathrm{\dots },{\mathit{\alpha }}_{\mathrm{top},K}$:

$$\begin{array}{}\text{(44)}& {\widehat{H}}_{\mathit{\alpha }}=\left\{\begin{array}{l}\sum _{\mathit{\beta }={\mathit{\alpha }}_{\mathrm{top},i}}^{{\mathit{\alpha }}_{\min ,K}}{l}_{\mathit{\beta },K}^{*}{\widehat{H}}_{{\mathit{\alpha }}_{\mathrm{top},K}}(T_{{\mathit{\beta }}^{*},i},T_{{\mathit{\beta }}^{*},j})\\ \mathrm{if}\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i}\mathrm{and}{\mathit{\alpha }}_{\min ,K}<{\mathit{\alpha }}_{\mathrm{top},i}\\ l_{\mathit{\alpha },K}^{*}{\widehat{H}}_{{\mathit{\alpha }}_{\mathrm{top},K}}(T_{{\mathit{\alpha }}^{*},i},T_{{\mathit{\alpha }}^{*},j})\\ \text{otherwise (hanging layer)}\end{array}\right..\end{array}$$

Again, we have separated the cases of a node with or without hanging layers. Note that the subscript ${\mathit{\alpha }}^{*}=max(\mathit{\alpha },{\mathit{\alpha }}_{\mathrm{top},j})$ avoids selecting tracer values in removed layers. In the Appendix we show that, when a constant tracer is imposed, the horizontal tracer flux reduces to the mass flux even in the case of a nonconformal box.

6 Numerical tests

The tests have been run with implicitness parameter θ=0.5. We check discrete mass-conservation at tⁿ⁺¹ by computing the following relative volume error for the dual cell area, which results from the sum of Eq. (28) from N_i to α_top,i:

$$e_i^{n+\mathrm{1}}=\mathrm{\Delta }t\left|{\displaystyle }\sum _{\mathit{\alpha }=N_i}^{{\mathit{\alpha }}_{\mathrm{top},i}}\left|C_{\mathit{\alpha },i}\right|G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}\right|,e^{n+\mathrm{1}}=\underset{i\in \mathcal{T}}{max}\left({\displaystyle \frac{e_i^{n+\mathrm{1}}}{{\displaystyle }\sum _{\mathit{\alpha }=N_i}^{{\mathit{\alpha }}_{\mathrm{top},i}}\left|C_{\mathit{\alpha },i}\right|h_{\mathit{\alpha },i}}}\right).$$

To quantify the tracer constancy error, we use the L¹−norm:

$$e^{n+\mathrm{1}}={\displaystyle \frac{{\sum }_{\mathit{\alpha },i}\left|C_{\mathit{\alpha },i}\right|h_{\mathit{\alpha },i}^{n+\mathrm{1}}|T_{\mathit{\alpha },i}^{n+\mathrm{1}}-T_{\mathrm{0}}|}{{\sum }_{\mathit{\alpha },i}\left|C_{\mathit{\alpha },i}\right|h_{\mathit{\alpha },i}^{n+\mathrm{1}}{T}_{\mathrm{0}}}},$$

with T₀ the initial tracer value.

6.1 Impulsive wave

As the first test, we check the accuracy of the z-surface-adaptive layers with an increasing vertical resolution. We use a closed basin $[-\mathrm{5}\mathrm{m},\mathrm{5}\mathrm{m}]\times [-\mathrm{5}\mathrm{m},\mathrm{5}\mathrm{m}]$ with constant depth z_b=1 m. The basin is initially at rest and the free surface is perturbed by the following Gaussian hump:

$$\mathit{\zeta }(\mathit{x},t=\mathrm{0})=A\exp (-r^{\mathrm{2}}/\mathit{\tau }),$$

with A=0.5 m, τ=0.5 m² and $r=\sqrt{x^{\mathrm{2}}+y^{\mathrm{2}}}$. A constant passive tracer is prescribed in the background and such a constant state should be preserved along the simulation. The mesh has a uniform horizontal element size of h_K=0.25 m. We compare different vertical resolutions with variable layer thickness. The coarsest grid has three layers: a first surface layer with a thickness of ΔZ₁=0.2 m, the second and the third layers have thicknesses of $\mathrm{\Delta }{Z}_{\mathrm{2}}=\mathrm{\Delta }{Z}_{\mathrm{3}}=\mathrm{0.4}$ m. The other vertical grids are obtained by halving each of these layers. The finest grid has 24 layers with a minimum layer thickness at the surface of ΔZ₁=0.025 m.

Without bottom and surface forcing, if the initial velocities are constant over the layers, they must remain barotropic and equal to the depth-integrated velocities of the shallow water equations (one-layer case). Of course, this is not a property of the z-layers (but the scheme should converge to a barotropic solution refining the resolution). It is however desirable that the results of 2D and 3D models are similar for the typical resolution of an ocean simulation (Kleptsova et al., 2010). The one-layer discrete solution is considered here as a reference solution against which we compare our implementation of the z-layers. The coarse grid with three layers is also used for comparison as the free surface is contained in the first layer and no insertion/removal is necessary. For the 24-layer grid, up to six layers are progressively removed (and then re-inserted). In Fig. 9, all resolutions show a good agreement for both the water level and the barotropic velocity. We can check some conservation properties of the scheme. As usual for such an adaptation strategy, mass is conserved up to machine precision (SHYFEM is coded in single-precision). This is what we check in Fig. 10. Except for a small additional noise associated with the insertion/removal operations, no significant source of mass error is present with respect to the three-layer case. Tracer constancy, as expected, is also preserved up to machine precision (see Fig. 10).

Figure 9 Impulsive wave. Comparison of the free surface elevation and barotropic velocity at different time instants. Vertical grids with different resolutions are compared. For each grid the reference interfaces $Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$ are traced with dashed lines. In the regions where the free surface crosses the interface $Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$ it means that the layer α locally has been removed from the computation.

Figure 10 Impulsive Wave. (a) Relative mass conservation error for the dual cell. (b) Relative tracer constancy error at the final time t=3s.

6.2 1-d tidal flow in a sloping channel

Coastal applications include extensive intertidal flats. As with many ocean models, SHYFEM handles wetting and drying processes in a simplified manner, applying special treatments in dry cells. An extrapolation algorithm for the free surface is used to track the shoreline and identify dry and wet elements. Then, the dry elements are taken out from the semi-implicit update and are treated in a mass-conserving manner as described in Umgiesser (2022). The test that we propose, presented in Oey (2005), is a benchmark for wetting-drying algorithms used in ocean models. The domain consists of a 1 d sloping channel that ranges from x=0 at the landward end to x=L at the seaward boundary, with L=25 km. The bathymetry is represented by the following function $z_{\mathrm{b}}\left(x\right)=H_{\mathrm{0}}/Lx$ and H₀=10 m. The horizontal element size is uniform and equal to h_K=250 m. A periodic water level is imposed at the seaward boundary as $\mathit{\zeta }\left(t\right)=A(\mathrm{1}-\sin \left(\frac{\mathrm{2}\mathit{\pi }t}{T}\right))$ with amplitude A=10 m, period T=1 d and the time t ranging from 0 to 0.5 d. At the beginning of the simulation, the channel is dry. Typically, this test is run with one-layer models (Warner et al., 2013). Here we use the one-layer solution (1L) as a reference and we test the five-layer solution with surface-adaptation and the five-layer solution with z-star. In the 5L z-surface-adaptive simulation, only one layer is present at the beginning of the simulation and then, as the free surface is tilted by the boundary signal, more levels are inserted and then removed during the drying phase. With z-star instead, the number of layers remains constant over time.

In Fig. 11 we check the along-channel solution profiles. Despite the different vertical resolution in the wet-dry and dry regions for the 5L z-surface-adaptive and 5L z-star simulations, quite good agreement is observed for the free surface. Larger differences are found for the barotropic velocity, where both the five-layer simulations appear noisier at the wet-dry interface. In Fig. 12 we check volume conservation for this case, which involves an uneven bathymetry and wetting-drying. Although in correspondence with wet-dry nodes the relative volume error is much larger, we can verify that the z-surface adaptive has the same level of relative error of z-star, which we accept to be within the round-off errors. The same argument applies to the error for the tracer constancy (see Fig. 12).

Figure 11 1 d tidal channel flow. Comparison between the one-layer and five-layer runs. (a, c, e) Free surface elevation. (b, d, f) Barotropic velocity. Dashed gray lines represent the reference interfaces $Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$. In the regions where the free surface crosses the interface $Z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$ it means that the layer α locally has been removed from the computation.

Figure 12 1 d tidal channel flow. (a) Relative mass conservation error for the dual cell. (b) Relative tracer constancy error at the final time.

6.3 Po delta idealized test

We test the different z-layers in a realistic coastal environment forced by the tidal oscillation: the Po delta. We study both the river plume and the penetration of the salt water into the river branches. The reproduction of such phenomena for numerical models is a very delicate issue. Specifically, spurious mixing related to the horizontal and vertical numerical fluxes, the vertical grid and the time-stepping can destroy stratification and frontal characteristics, potentially modifying the plume dynamics (Fofonova et al., 2021). In this discussion we solely focus on the impact of the vertical discretization: the resolution at the surface and the comparison between the z-surface adaptive with fixed interfaces and z-star with moving interfaces.

The vertical eddy viscosity and the vertical tracer eddy diffusivity are computed with the turbulence module General Ocean Turbulence Model (GOTM) (Buchard et al., 2001). The bottom friction is fixed to C_F=0.002. Because of their fundamental role in the plume dynamics, two more terms have been added to the multilayer shallow water model of Sect. 2: the Coriolis force, which is time stepped with an implicitness parameter of 0.5 and a horizontal diffusion term for the salinity equation, treated explicitly. The horizontal viscosity is taken as the Smagorinsky eddy viscosity. The sea boundary is forced with a semi-diurnal tidal signal with amplitude 0.4 m and period 12 h. The salinity at the sea boundary is constant and fixed to 38 PSU. A weak freshwater flow with a discharge of 500 m³ s⁻¹, which is characteristic of the summer season, is enforced at the Pontelagoscuro river boundary. The initial solution corresponds to water at rest and salinity equal to the boundary value of 38 PSU. The simulation lasts 1 month, after which the salinity shows periodic behavior modulated by the tidal cycle.

The computational domain encompasses the entire river network of the delta, stretching from Pontelagoscuro to the sea, including all delta lagoons, as well as a portion of the adjacent shelf sea (Bellafiore et al., 2021). Horizontal resolution ranges from h_K=2 km at the sea boundary, to h_K=100 m in the inner shelf close to the lagoons and river branches, and to h_K=50 m in the inner delta system. The horizontal grid, composed of 38 884 nodes and 69 364 elements, is in Fig. 13. We consider two vertical resolutions, one with N=24 layers and one with N=27 layers. The deeper part (from the bottom to $Z=-\mathrm{1}$ m) is equal for the two grids and it is composed of 23 levels with variable thicknesses from ΔZ=0.5 m near the surface up to ΔZ_N=4 m for the last layer. The resolution of the upper part of the water column differs: the 24-layer grid has one layer with ΔZ₁=1 m. This choice avoids the drying of the first layer. The 27-layer grid, in the upper part, has four layers with constant thickness, $\mathrm{\Delta }{Z}_{\mathrm{1}}=\mathrm{\Delta }{Z}_{\mathrm{2}}=\mathrm{\Delta }{Z}_{\mathrm{3}}=\mathrm{\Delta }{Z}_{\mathrm{4}}=\mathrm{0.25}$ m. Three simulations have been performed: one with 24 standard z-layers (24L z), one with 27 z-surface-adaptive layers (27L z-surf-adapt), and one with 27 z-star layers (27L z-star).

Given the fine vertical resolution at the surface and the tidal amplitude of 0.4 m, the 27L z-surf-adapt simulation should undergo extensive insertion/removal of the surface fluid boxes. In the left picture of Fig. 14 we have reported the time evolution of the number of boxes inserted and removed during two tidal periods. Almost 4000 surface boxes happened to be inserted or removed in a single time step. As is customary we have reported mass conservation and tracer constancy error in Fig. 14. These latest results refer to a shorter simulation that lasted 4 d with a constant salinity obtained by imposing the river salinity equal to the interior one.

Figure 13 Po river. (a) Horizontal grid. (b) Zoom of the horizontal grid with tidal stations and the transect in the Pila branch.

Figure 14 Po river. (a) Time evolution of the total number of layers inserted and removed per time step for the 24L z−surf-adapt simulation. (b) Relative mass conservation error for the dual cell. (c) Relative tracer constancy error after 4 d.

To diagnose the river plume we look at the minimum surface salinity during the simulation. From Fig. 15, it is clear that both 27-layer simulations allow a stronger gravitational circulation with a more extended freshwater plume. Also, the opposite bottom circulation penetrates more upstream, with stronger salinity recorded at the stations G2 and G5, as shown in Fig. 16. To inspect the extension of the saltwater intrusion we have extracted a section of the salinity field in the Pila branch when saltwater reaches the maximum extent, during a flood tide. This is shown in Fig. 17. The higher resolution at the surface also allows some small-scale internal structures that are present under the surface to be captured. The z-surface adaptive simulation exhibits a stronger plume and a more extended salt wedge as well as a sharper surface structure. Similar results can be observed in a recent work (Verri et al., 2023), where standard z-layers and z-star are compared for an analogous river plume experiment. A possible explanation could be related to the fact that, owing to the strong internal motion, the vertical velocity is not in phase with the time derivative of the free surface and it may have an opposite sign with respect to the grid velocity. For this particular case, the z-star mass-transfer function Eq. (11) could be larger than the vertical velocity. In turn, this can be related to a larger multiplicative constant in the truncation error associated with the vertical advection scheme.

Figure 15 Po river. Minimum of the surface salinity (for the 24-layer grid the minimum is computed at the first layer, for the 27-layer grid at the second layer). (a) 24L z. (b) 27L z-surf-adapt. (c) 27L z-star.

Figure 16 Po river. Salinity profile at G2 (a–c) and G5 (d–f). (a, d) 24L z. (b, e) 27L z-surf-adapt. (c, f) 27L z-star.

Figure 17 Salinity section along the Pila branch during the flood tide of day 29 at 16:00. (a) 24L z. (b) 27L z-surf-adapt. (c) 27L z-star.

All the tests have been accomplished with a serial run. We report the CPU time of the serial simulations, which have been run on a modern workstation with a AMD EPYC 7643 Processor 2 073 005 s (24L z-star), 1 998 969 s (24L z-surf-adapt) showing an overhead of around 3.6 % for the insertion/removal operations. Although we have not covered parallel implementation aspects, we mention that the algorithm (grid movement, insertion/removal) mainly operates on the vertical grid, and the parallel execution of these tasks should not encounter any issues. The stencil of the numerical scheme is not enlarged with respect to the standard method. However, some variables have been introduced only for the insertion/removal operations. This is the case of the nodal top layer index, which must be exchanged between the domains.

7 Conclusions

In this work, we have studied the performances of multilayer shallow water models based on z-layers for the simulation of free surface coastal flows. We have investigated a well-known issue of z-layers when incorporating the free surface: the limitation on the resolution of the surface layer thickness. We have proposed a flexible algorithm based on a vertical mesh adaptation to the tidal oscillation called z-surface-adaptive. With a dynamic insertion and removal of surface layers, the grid at the new timestep is aligned to the geopotentials, reducing the pressure gradient error. Thanks to a two-step procedure (vertical grid movement of surface layers followed by the insertion/removal operations), we have been able to evolve the multilayer model on a grid with a temporally constant number of layers in the time step, which allowed a simple implementation. Moreover, this leads to a consistency, at a discrete level, of the tracer equation with the continuity equation as well as to a simple verification of mass-conservation. As a particular case, the algorithm can be reduced to the popular z-star.

Without the limitation on the surface resolution, we have been able to compare the z-layers with insertion/removal (surface-adaptive) against z-star for typical coastal applications of semi-enclosed shallow seas with a tidal signal imposed at the openings and wetting-drying at intertidal flats. The comparison has been carried out with idealized and realistic numerical experiments. We show that z-surface-adaptive layers can be used to simulate wetting and drying without a significant loss of accuracy with respect to z-star. We found that z-layers and z-star exhibit differences when simulating large, low-frequency internal motions combined with a barotropic tide, such as the gravitational circulation in the Po Delta. These differences deserve further investigation. We speculate that for such cases, keeping z-layers may be convenient to reduce truncation errors in the computation of both the internal pressure gradient term and the vertical advection terms.

We conclude by mentioning that the overhead related to insertion/removal operation should be further assessed in realistic applications. With the actual implementation of the z-surface adaptive layers, we experienced some stability issues in the computation of the tracers. This occurred for nonconformal boxes undergoing wetting-drying and is under current investigation. We are trying a simpler treatment of the nonconformal surface boxes as in Bonaventura et al. (2018).

Appendix A: Tracer constancy

We start with the case without nonconformal boxes. We impose a constant tracer vector T_i=1 in the discrete tracer Eq. (31). Each row reduces to

$$\left|C_{\mathit{\alpha },i}\right|h_{\mathit{\alpha },i}^{n+\mathrm{1}}=\left|C_{\mathit{\alpha },i}\right|h_{\mathit{\alpha },i}^n+\mathrm{\Delta }t{f}_{\mathit{\alpha },i}^n,$$

with

$$\begin{array}{rl}{f}_{\mathit{\alpha },i}^n& =-\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}\sum _{j\in K,j\ne i}{\widehat{H}}_{\mathit{\alpha }}\left(\mathrm{1},\mathrm{1}\right)\\ & +\left(\left|C_{\mathit{\alpha },i}\right|G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-\left|C_{\mathit{\alpha }+\mathrm{1},i}\right|G_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}\right)\end{array}.$$

Using, first, the numerical flux consistency ${\widehat{H}}_{\mathit{\alpha }}\left(\mathrm{1},\mathrm{1}\right)={\mathit{q}}_{\mathit{\alpha }}^{n+\mathit{\theta }}\cdot {\mathit{n}}_{ij}^K$ and then the relationship between the element normals and the dual cell ones (Eq. 19):

$$\begin{array}{rl}\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}& \sum _{j\in K,j\ne i}{\widehat{H}}_{\mathit{\alpha }}\left(\mathrm{1},\mathrm{1}\right)=\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}\sum _{j\in K,j\ne i}{\mathit{q}}_{\mathit{\alpha }}^{n+\mathit{\theta }}\cdot {\mathit{n}}_{ij}^K\\ & =-\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}{\mathit{q}}_{\mathit{\alpha }}^{n+\mathit{\theta }}\cdot {\displaystyle \frac{{\mathit{n}}_i^K}{\mathrm{2}}}\\ & =-\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}\left(a_{iK}^x{q}_{\mathit{\alpha },K}^{x,n+\mathit{\theta }}+a_{iK}^y{q}_{\mathit{\alpha },K}^{y,n+\mathit{\theta }}\right)\end{array}.$$

In the last step we used the fact that for piecewise linear basis functions we have $\frac{{\mathit{n}}_i^K}{\mathrm{2}}=\left|K\right|{\left.\mathrm{\nabla }{\mathit{\phi }}_i\right|}_K$. For each element in the subset 𝒟_α,i, the horizontal tracer flux has been reduced to the mass flux. We can write the discrete tracer update:

$$\begin{array}{rl}\left|C_{\mathit{\alpha },i}\right|& {\displaystyle \frac{\mathrm{\Delta }{h}_{\mathit{\alpha },i}}{\mathrm{\Delta }t}}=\sum _{K\in {\mathcal{D}}_{\mathit{\alpha }i}}\left(a_{iK}^x{q}_{\mathit{\alpha },K}^{x,n+\mathit{\theta }}+a_{iK}^y{q}_{\mathit{\alpha },K}^{y,n+\mathit{\theta }}\right)\\ & +\left|C_{\mathit{\alpha },i}\right|G_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-\left|C_{\mathit{\alpha }+\mathrm{1},i}\right|G_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}\end{array},$$

which corresponds to the discrete layerwise mass Eq. (28).

In the case of a nonconformal box, we have to show that the modified horizontal tracer fluxes still reduces to the mass-fluxes. According to Eq. (44), the horizontal tracer fluxes in nonconformal boxes should be computed with

$${\widehat{H}}_{\mathit{\alpha }}=\left\{\begin{array}{l}\sum _{\mathit{\beta }={\mathit{\alpha }}_{\mathrm{top},i}}^{{\mathit{\alpha }}_{\min ,K}}{l}_{\mathit{\beta },K}^{*}{\widehat{H}}_{{\mathit{\alpha }}_{\mathrm{top},K}}(T_{{\mathit{\beta }}^{*},i},T_{{\mathit{\beta }}^{*},j})\\ \mathrm{if}\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i}\mathrm{and}{\mathit{\alpha }}_{\min ,K}<{\mathit{\alpha }}_{\mathrm{top},i}\\ l_{\mathit{\alpha },K}^{*}{\widehat{H}}_{{\mathit{\alpha }}_{\mathrm{top},K}}(T_{{\mathit{\alpha }}^{*},i},T_{{\mathit{\alpha }}^{*},j})\\ \text{otherwise (hanging layer)}\end{array}.\right.$$

For a constant tracer, it can be rewritten for $\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{top},i},\mathrm{\dots }{\mathit{\alpha }}_{\mathrm{top},K}$:

$${\widehat{H}}_{\mathit{\alpha }}=c_{\mathit{\alpha },i}^{*}{\widehat{H}}_{{\mathit{\alpha }}_{\mathrm{top},K}}\left(\mathrm{1},\mathrm{1}\right),$$

where we have also used the definition (Eq. 43). Thus,

$$\sum _{j\in K,j\ne i}{c}_{\mathit{\alpha },i}^{*}{\widehat{H}}_{{\mathit{\alpha }}_{\mathrm{top},K}}\left(\mathrm{1},\mathrm{1}\right)=-c_{\mathit{\alpha },i}^{*}\left(a_{iK}^x{q}_{{\mathit{\alpha }}_{\mathrm{top},K},K}^{x,n+\mathit{\theta }}+a_{iK}^y{q}_{{\mathit{\alpha }}_{\mathrm{top},K},K}^{y,n+\mathit{\theta }}\right).$$

This gives exactly the contribution from nonconformal boxes to the mass-transfer (Eq. 42).

Finally, the tracer remap (Eq. 39) preserves the constancy property. It is enough to verify that with a constant solution it reduces to

$${\stackrel{\mathrm{̃}}{h}}_{\mathit{\alpha },i}^{n+\mathrm{1}}=h_{\mathit{\alpha },i}^{n+\mathrm{1}}+\mathrm{\Delta }t\left({\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{top}}-{\mathit{\sigma }}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{\mathrm{top}}\right),$$

which, thanks to the definition provided in Sect. 4.2 of grid velocity ${\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{\mathrm{top}}=\frac{{\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-z_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}}{\mathrm{\Delta }t}$ and layer thickness ${\stackrel{\mathrm{̃}}{h}}_{\mathit{\alpha },i}^{n+\mathrm{1}}={\stackrel{\mathrm{̃}}{z}}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}-z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}$, is an identity.