the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A flexible zlayers approach for the accurate representation of free surface flows in a coastal ocean model (SHYFEM v. 7_5_71)
Christian Ferrarin
Marco Bajo
Georg Umgiesser
We propose a discrete multilayer shallow water model based on zlayers, 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 twostep 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 zsurfaceadaptive, can be reduced, as a particular case when all layers are moving, to the zstar coordinate. With idealized and realistic numerical experiments, we compare the zsurfaceadaptive against zstar and we show that it can be used to simulate coastal flows effectively.
 Article
(6963 KB)  Fulltext XML
 BibTeX
 EndNote
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) zlayers with fixed interfaces parallel to geopotentials (Eulerian framework); (3) terrain/surfacefollowing σ or slayers 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.
zlayers 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 Mudflat3D, TRIM3D, Cheng et al., 1993; UNTRIM3D, 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 lowfrequency motions (Leclair and Madec, 2011). This occurs because the isopycnals are well aligned to the zinterfaces 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 zlayers performances relative to the treatment of the free surface boundary. To simplify the boundary condition at the free surface, zlayers 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 ztype coordinates have been introduced over the years: they are based on zlayers that move to accommodate the tidal oscillation, but the bottom is not a coordinate surface (they are surfacefollowing but not terrainfollowing). These coordinates are clearly of ALEtype 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 zstar of Adcroft and Campin (2004), the quasiz of Mellor et al. (2002), and the hybrid $z/\mathit{\sigma}$ of Burchard and Petersen (1997) all belong to such zsurfacefollowing 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 zsurfaceadaptive. 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 zsurface 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 onedimensional 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 zlayers in the sense that the same algorithm can be easily reduced to zstar 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 zstar.
The algorithm is implemented in the SHYFEM finiteelement ocean model of the CNRISMAR (Umgiesser et al., 2004, https://github.com/SHYFEMmodel/shyfem, last access: 5 November 2023), which implements the multilayer shallow water equations with z and σ layers. SHYFEM uses a semiimplicit finite element discretization on unstructured Btype 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 zstar and the zlayers. In Sect. 3 we provide the semiimplicit finite element discretization of the multilayer equations. In Sect. 4 we describe the zsurfaceadaptive 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.
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 Ω:
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 nonoverlapping 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 timedependent 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:
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 fixedintime 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:
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$:
To prescribe this function we use the generalized vertical coordinate transformation, see Mellor et al. (2002):
which assures a surface and terrainfollowing 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.
Using this transformation, the layer thickness can be deduced from the water depth, for $\mathit{\alpha}=\mathrm{1},\mathrm{\dots},N$:
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$:
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$:
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ándezNieto et al. (2014):
$${f}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}}^{+}={\left({\left(\right)}_{f}{\mathrm{\Omega}}_{\mathit{\alpha}}\right)}_{}{\mathrm{\Gamma}}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}},\phantom{\rule{1em}{0ex}}{f}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}}^{}={\left({\left(\right)}_{f}{\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(\right)}_{f}{\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
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$:
where ${G}_{\mathit{\alpha}\pm \mathrm{1}/\mathrm{2}}$ is the masstransfer 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 layerintegrated 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 righthand side.
From the derivation of FernándezNieto et al. (2014), the definition of the masstransfer function is
with ${\mathit{\sigma}}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}}$ the velocity of the grid interface:
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:
which are evaluated starting from the bottom $\mathit{\alpha}=N,\mathrm{\dots},\mathrm{1}$, where the noslip 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 masstransfer function is computed directly from the layerintegrated mass equation
Summing from N to α as
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 postprocessing 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
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}}$:
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
with C_{F} the bottom friction coefficient. Similarly, the diffusive fluxes read
with no tracer fluxes through the free surface and the bottom.
Finally, the term B_{α} represents the internal pressure gradient force. The layerintegrated pressure gradient term ${\int}_{{z}_{\mathit{\alpha}+\mathrm{1}/\mathrm{2}}}^{{z}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}}}\mathrm{\nabla}p\left(z\right)\phantom{\rule{0.125em}{0ex}}\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):
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
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 socalled 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 zstar
The multilayer model presented so far is based on a vertical subdivision of the fluid domain through the surface and terrainfollowing 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 zstar (Adcroft and Campin, 2004). The reference domain, with vertical coordinate Z, is
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:
As shown in Fig. 2, there is a substantial difference with respect to the vertical subdivision of the terrainfollowing 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
If a layer is bounded laterally by the bathymetry interface we can denote this lateral land boundary of the layer as
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$:
with ${S}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}}$ a stretching function defined as
As in the previous Section, the layer thickness can be deduced from the total water depth. After some calculations we get
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$:
which is prescribed once the reference grid is created. The coefficients satisfy both the positivity constraint and locally they sum to one.
An important property of the zstar 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 zlayers formulations based on similar mappings, such as the quasiz 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 zstar 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 zlayers
For zlayers the actual interfaces do not depend on time and space:
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 masstransfer function (Eq. 14) coincides with the vertical velocity:
Replacing the masstransfer function with the vertical velocity in the multilayer model, we obtain the Eulerian model of Rambaud (2011).
The discretization for both the zstar and the zlayers shallow water model can proceed in an equivalent fashion. We consider a discretization of the horizontal domain Ω_{x}∈ℝ^{2} composed of nonoverlapping triangular elements. We denote the horizontal grid by 𝒯 with K∈𝒯 the generic triangle, and $\leftK\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
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 deﬁnition of the dual cell, we have
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, zlayers models apply a mask to nonexisting 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
On a Bstaggered 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):
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
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 Bgrid 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 xcomponent shown):
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
The horizontal advection term is resolved with a firstorder 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:
and analogously the explicit terms:
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 semiimplicit time discretization by treating semiimplicitly the mass flux and the free surface gradient in the momentum equation. The vertical viscous term can also cause a restrictive timestep 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:
After applying the previous definition to the semidiscrete equations, the semiimplicit momentum equations on an unstructured Bgrid read
with ${\mathbf{A}}_{K}=\left(\mathit{I}\leftK\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^{n}. In ${\mathit{F}}_{K}^{n}$ all the quantities are computed at t^{n}, including the masstransfer 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:
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:
where ${m}_{ij}^{K}={\int}_{K}{\mathit{\phi}}_{i}{\mathit{\phi}}_{j}\phantom{\rule{0.125em}{0ex}}\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 semiimplicit step Eqs. (23) and (24). Finally, we compute the layer thickness at the grid nodes. For zstar, we use the expression Eq. (18) at the grid nodes. For the zlayers, the layer thickness does not change except for the first layer.
3.1 Masstransfer function
After the hydrodynamic update of the previous paragraph, the discrete masstransfer function is computed. We employ the same continuous piecewise linear approximation used for the free surface. The nodal values are computed from a finiteelement masslumped discretization of the layerwise mass Eq. (14). As for the depthintegrated mass equation, the discharge is evaluated semiimplicitly. 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}$:
Note that the semiimplicit discretization ensures vertical massconservation. Summing up Eq. (28) for all the layers and using Eq. (27) with a lumped Galerkin massmatrix to cancel the righthand side, we get the impermeability condition at the free surface ${G}_{\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}=\mathrm{0}$. With standard zlayers, 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 semiimplicit update is completed with the timestepping 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:
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
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 secondorder 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:
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
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 masstransfer 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 masstransfer function Eq. (28) and, last, we update the tracers with Eq. (31).
In this section, we enhance the zlayers 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 zlayers, we refer to this enhanced version as zsurfaceadaptive 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{\u0303}}{z}}_{\mathit{\alpha}\mathrm{1}/\mathrm{2},i}^{n+\mathrm{1}}={\stackrel{\mathrm{\u0303}}{z}}_{\mathit{\alpha}\mathrm{1}/\mathrm{2}}({\mathit{x}}_{i},{t}^{n+\mathrm{1}})$ is
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{\u0303}}{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
with
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^{n+1} 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^{n} and t^{n+1}), 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.
In the following we provide the technical details to realize such adaptation to the free surface with the zlayers. 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 ζ^{0}(x), we define a set of active indices and the surface layer index, by element, as
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:
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 onedimensional example.
4.1 Vertical grid movement
We evolve the discrete multilayer shallow water equations using the semiimplicit 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{\u03f5}}_{\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 bottommost 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}}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}{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),\phantom{\rule{1em}{0ex}}\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 zstar 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 reevaluating Eqs. (32) and (33) with the new free surface elevation ζ^{n+1}. 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{\u0303}}{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:
with an upwind flux:
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
with the new nodal layer thickness:
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 depthintegrated 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 depthintegrated 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 zstar
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 zsurfaceadaptive 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 zlayers having all the interfaces aligned to the geopotentials.
Interestingly, we can obtain other grids by increasing r_{mov}. We define
with ${\mathit{\zeta}}_{\mathrm{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

zstar 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.

zstar + 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 zsurface adaptive layers against zstar.
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^{n} and t^{n+1}. 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)).
Consider the one dimensional example in Fig. 7, where two contiguous elements with a different toplayer 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 toplayer 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:
with ${\mathit{\alpha}}_{\mathrm{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}}_{\mathrm{min},K},\mathrm{\dots},{\mathit{\alpha}}_{\mathrm{top},K}$:
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 massflux term (only the xcomponent shown) reads, for $\mathit{\alpha}={\mathit{\alpha}}_{\mathrm{top},i},\mathrm{\dots},{\mathit{\alpha}}_{\mathrm{top},K}$:
with
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 massflux term in the layerwise mass equation. It allows the masstransfer 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 massconserving. Summing the masstransfer 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}$:
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.
The tests have been run with implicitness parameter θ=0.5. We check discrete massconservation at t^{n+1} 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}:
To quantify the tracer constancy error, we use the L^{1}−norm:
with T_{0} the initial tracer value.
6.1 Impulsive wave
As the first test, we check the accuracy of the zsurfaceadaptive layers with an increasing vertical resolution. We use a closed basin $[\mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{m},\mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{m}]\times [\mathrm{5}\phantom{\rule{0.125em}{0ex}}\mathrm{m},\mathrm{5}\phantom{\rule{0.125em}{0ex}}\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:
with A=0.5 m, τ=0.5 m^{2} 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_{1}=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_{1}=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 depthintegrated velocities of the shallow water equations (onelayer case). Of course, this is not a property of the zlayers (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 onelayer discrete solution is considered here as a reference solution against which we compare our implementation of the zlayers. 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 24layer grid, up to six layers are progressively removed (and then reinserted). 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 singleprecision). 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 threelayer case. Tracer constancy, as expected, is also preserved up to machine precision (see Fig. 10).
6.2 1d 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 semiimplicit update and are treated in a massconserving manner as described in Umgiesser (2022). The test that we propose, presented in Oey (2005), is a benchmark for wettingdrying 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}}/L\phantom{\rule{0.125em}{0ex}}x$ and H_{0}=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}\mathrm{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 onelayer models (Warner et al., 2013). Here we use the onelayer solution (1L) as a reference and we test the fivelayer solution with surfaceadaptation and the fivelayer solution with zstar. In the 5L zsurfaceadaptive 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 zstar instead, the number of layers remains constant over time.
In Fig. 11 we check the alongchannel solution profiles. Despite the different vertical resolution in the wetdry and dry regions for the 5L zsurfaceadaptive and 5L zstar simulations, quite good agreement is observed for the free surface. Larger differences are found for the barotropic velocity, where both the fivelayer simulations appear noisier at the wetdry interface. In Fig. 12 we check volume conservation for this case, which involves an uneven bathymetry and wettingdrying. Although in correspondence with wetdry nodes the relative volume error is much larger, we can verify that the zsurface adaptive has the same level of relative error of zstar, which we accept to be within the roundoff errors. The same argument applies to the error for the tracer constancy (see Fig. 12).
6.3 Po delta idealized test
We test the different zlayers 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 timestepping 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 zsurface adaptive with fixed interfaces and zstar 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 semidiurnal 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^{3} s^{−1}, 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 24layer grid has one layer with ΔZ_{1}=1 m. This choice avoids the drying of the first layer. The 27layer 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 zlayers (24L z), one with 27 zsurfaceadaptive layers (27L zsurfadapt), and one with 27 zstar layers (27L zstar).
Given the fine vertical resolution at the surface and the tidal amplitude of 0.4 m, the 27L zsurfadapt 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.
To diagnose the river plume we look at the minimum surface salinity during the simulation. From Fig. 15, it is clear that both 27layer 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 smallscale internal structures that are present under the surface to be captured. The zsurface 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 zlayers and zstar 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 zstar masstransfer 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.
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 zstar), 1 998 969 s (24L zsurfadapt) 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.
In this work, we have studied the performances of multilayer shallow water models based on zlayers for the simulation of free surface coastal flows. We have investigated a wellknown issue of zlayers 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 zsurfaceadaptive. 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 twostep 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 massconservation. As a particular case, the algorithm can be reduced to the popular zstar.
Without the limitation on the surface resolution, we have been able to compare the zlayers with insertion/removal (surfaceadaptive) against zstar for typical coastal applications of semienclosed shallow seas with a tidal signal imposed at the openings and wettingdrying at intertidal flats. The comparison has been carried out with idealized and realistic numerical experiments. We show that zsurfaceadaptive layers can be used to simulate wetting and drying without a significant loss of accuracy with respect to zstar. We found that zlayers and zstar exhibit differences when simulating large, lowfrequency 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 zlayers 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 zsurface adaptive layers, we experienced some stability issues in the computation of the tracers. This occurred for nonconformal boxes undergoing wettingdrying and is under current investigation. We are trying a simpler treatment of the nonconformal surface boxes as in Bonaventura et al. (2018).
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
with
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):
In the last step we used the fact that for piecewise linear basis functions we have $\frac{{\mathit{n}}_{i}^{K}}{\mathrm{2}}=\leftK\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:
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 massfluxes. According to Eq. (44), the horizontal tracer fluxes in nonconformal boxes should be computed with
For a constant tracer, it can be rewritten for $\mathit{\alpha}={\mathit{\alpha}}_{\mathrm{top},i},\mathrm{\dots}{\mathit{\alpha}}_{\mathrm{top},K}$:
where we have also used the definition (Eq. 43). Thus,
This gives exactly the contribution from nonconformal boxes to the masstransfer (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
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{\u0303}}{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{\u0303}}{h}}_{\mathit{\alpha},i}^{n+\mathrm{1}}={\stackrel{\mathrm{\u0303}}{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.
The SHYFEM hydrodynamic model is open source (GNU General Public License as published by the Free Software Foundation) and freely available through GitHub at https://github.com/SHYFEMmodel (last access: 5 November 2023). The current developments have been implemented in a branch of the SHYFEM code that can be accessed from Zenodo (https://doi.org/10.5281/zenodo.8356398, Arpaia, 2023). Configuration files and data used to run each test case are also available at the same Zenodo repository.
LA: Conceptualization, Methodology, Software, Validation, Writing, Formal analysis. CF: Conceptualization, Methodology, Funding acquisition, Writing, Resources, Validation. MB: Methodology, Writing. GU: Conceptualization, Methodology, Writing, Software.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
All the developments presented have been implemented in the Finite Element Model for Coastal Seas SHYFEM (https://github.com/SHYFEMmodel/shyfem, last access: 8 November 2023) developed at the CNRISMAR. The authors acknowledge Debora Bellafiore for fruitful discussions about the implementation of the present work. The corresponding author expresses his gratitude to Luca Bonaventura and Giacomo Capodaglio, the two reviewers, for their valuable comments and feedback that contributed to improve the precision and clarity of the manuscript during the revision process.
This research has been supported by the project AdriaClim (Climate change information, monitoring and management tools for adaptation strategies in Adriatic coastal areas; project ID 10252001) funded by the European Union under the VA Interreg ItalyCroatia CBC programme.
This paper was edited by Vassilios Vervatis and reviewed by Giacomo Capodaglio and L. Bonaventura.
Adcroft, A. and Campin, J.M.: Rescaled height coordinates for accurate representation of freesurface flows in ocean circulation models, Ocean Model., 7, 269–284, 2004. a, b, c
Arpaia, L.: SHYFEM version with surfaceadaptive zcoordinates (7_5_71zlay1.2), Zenodo [code], https://doi.org/10.5281/zenodo.8356398, 2023. a
Audusse, E., Bristeau, M.O., Pelanti, M., and SainteMarie, J.: Approximation of the hydrostatic NavierStokes system for density stratified flows by a multilayer model: Kinetic interpretation and numerical solution, J. Comput. Phys., 230, 3453–3478, 2011. a, b
Backhaus, J. O.: A threedimensional model for the simulation of shelf sea dynamics., Dt. Hydrogr. Z., 38, 165–187, 1985. a
Bellafiore, D., Ferrarin, C., Maicu, F., Manfè, G., Lorenzetti, G., and et al., G. U.: Saltwater intrusion in a Mediterranean delta under a changing climate, J. Geophys. Res.Oceans, 126, 6945–6975, 2021. a
Bonaventura, L., FernandezNieto, E. D., GarresDiaz, J., and NarbonaReina, G.: Multilayer shallow water models with locally variable number of layers and semiimplicit time discretization, J. Comput. Phys., 364, 209–234, 2018. a, b
Buchard, H., Bolding, K., and Villareal, M. R.: GOTM, a General Ocean Turbulence Model. Theory, implementation and test cases, GOTM Report, 2001. a
Burchard, H. and Baumert, H.: The formation of estuarine turbidity maxima due to density effects in the salt wedge. A hydrodynamic process study, J. Phys. Oceanogr., 28, 309–321, 1998. a
Burchard, H. and Petersen, O.: Hybridization between sigma and zcoordinates for improving the internal pressure gradient calculation in marine models with steep bottom slopes, Int. J. Numer Meth. Fl., 25, 1003–1023, 1997. a, b, c
Casulli, V. and Cattani, E.: Stability, accuracy and efficiency of a semiimplicit method for threedimensional shallow water flow, Comput. Math. Appl., 27, 99–112, 1994. a
Casulli, V. and Cheng, R.: Semiimplicit finite difference methods for threedimensional shallow water flow, Int. J. Numer. Meth. Fluids, 15, 629–648, 1992. a, b
Casulli, V. and Walters, R. A.: An unstructured grid, threedimensional model based on the shallow water equations, Int. J. Numer. Meth. Fluids, 32, 331–348, 2000. a
Cheng, R., Casulli, V., and Gartner, J. W.: Tidal, Residual, Intertidal Mudflat (TRIM) model and its applications to San Francisco Bay, California, Estuar., Coast. Shelf S., 36, 235–280, 1993. a
Debreu, L., Kevlahan, N.R., and Marchesiello, P.: Brinkman volume penalization for bathymetry in threedimensional ocean models, Ocean Model., 145, 101530, https://doi.org/10.1016/j.ocemod.2019.101530, 2020. a
FernándezNieto, E., Koné, E., and Rebollo, T. C.: A Multilayer Method for the Hydrostatic NavierStokes Equations: A Particular Weak Solution, J. Sci. Comput., 60, 408–437, 2014. a, b
Fofonova, V., Kärnä, T., Klingbeil, K., Androsov, A., Kuznetsov, I., Sidorenko, D., Danilov, S., Burchard, H., and Wiltshire, K. H.: Plume spreading test case for coastal ocean models, Geosci. Model Dev., 14, 6945–6975, https://doi.org/10.5194/gmd1469452021, 2021. a
Griffies, S., Pacanowski, R., Schmidt, M., and Balaji, V.: Tracer conservation with an explicit freesurface method for zcoordinate ocean models, Mon. Weather Rev., 129, 1081–1098, 2001. a, b
Gross, E., Bonaventura, L., and Rosatti, G.: Consistency with continuity in conservative advection schemes for freesurface models, J. Comput. Phys., 38, 307–327, 2002. a, b
Guardone, A., Isola, D., and Quaranta, G.: Arbitrary Lagrangian Eulerian formulation for twodimensional flows using dynamic meshes with edge swapping, J. Comput. Phys., 230, 7706–7722, 2011. a
Hordoir, R., Axell, L., Loptien, U., Dietze, H., and Kuznetsov, I.: Influence of sea level rise on the dynamics of salt inflows in the Baltic Sea, J. Geophys. Res.Oceans, 120, 6653–6668, 2015. a
Kleptsova, O., Stelling, G., and Pietrzak, D.: An accurate momentum advection scheme for a zlevel coordinate models, Ocean Dynam., 60, 1447–1461, 2010. a
Klingbeil, K., Lemarié, F., Debreu, L., and Burchard, H.: The numerics of hydrostatic structuredgrid coastal ocean models: state of the art and future perspectives, Ocean Model., 125, 80–105, 2018. a, b
Leclair, M. and Madec, G.: zCoordinate, an Arbitrary Lagrangian–Eulerian coordinate separating high and low frequency motions, Ocean Model., 37, 139–152, 2011. a
LeVeque, R. J.: Finite Volume Methods for Hyperbolic Problems, Cambridge University Press, https://doi.org/10.1017/CBO9780511791253, 2002. a
Lin, S. J. and Rood, R. B.: Multidimensional flux form semiLagrangian transport schemes, Mon. Weather Rev., 124, 2046–2070, 1996. a
Mellor, G., Hakkinen, S., Ezer, T., and Patchen, R.: A generalization of a sigma coordinate ocean model and an intercomparison of model vertical grids, in: Ocean Forecasting: Conceptual Basis and Applications, edited by: Pinardi, N. and Woods, J., Springer, New York, 55–72, ISBN 9783642087547, 2002. a, b, c
Millero, F. J. and Poisson, A.: International oneatmosphere equation of state of seawater, DeepSea Res., 28, 625–629, 1981. a
Oey, L.Y.: A wetting and drying scheme for POM, Ocean Model., 2, 133–150, 2005. a
Rambaud, A.: Modélisation, analyse mathématique et simulations numériques de quelques problèmes aux dérivées partielles multiéchelles, PhD thesis, Université Claude Bernard – Lyon I, https://theses.hal.science/tel00656013 (last access: 6 April 2023), 2011. a
Shchepetkin, A. and McWilliams, J.: A method for computing horizontal pressuregradient force in an oceanic model with a nonaligned vertical coordinate, J. Geophys. Res., 108, https://doi.org/10.1029/2001JC001047, 2003. a
Song, Y. T.: A general pressure gradient formulation for ocean models: scheme design and diagnostic analysis, Mon. Weather Rev., 126, 3213–3230, 1998. a
Umgiesser, G.: SHYFEM Finite Element Model for Coastal Seas – User Manual, Tech. Rep., Oceanography, ISMARCNR Arsenale Tesa 104, Castello 2737/F 30122 Venezia, Italy, 2022. a, b
Umgiesser, G., Canu, D. M., Cucco, A., and Solidoro, C.: A finite element model for the Venice Lagoon. Development, set up, calibration and validation, J. Marine Syst., 51, 123–145, 2004. a, b
Verri, G., Barletta, I., Pinardi, N., Federico, I., Alessandri, J., and Coppini, G.: Shelf slope, estuarine dynamics and river plumes in a z^{*} vertical coordinate, unstructured grid model, Ocean Model., 184, 102235, https://doi.org/10.1016/j.ocemod.2023.102235 2023. a
Warner, J., Defne, Z., Haas, K., and Arango, H.: A wetting and drying scheme for ROMS, Comput. Geosci., 58, 54–61, 2013. a
Williams, R. T. and Zienkiewicz, O. C.: Improved finite element forms for the shallowwater wave equations, Int. J. Numer. Meth. Fluids, 1, 81–97, 1981. a, b
 Abstract
 Introduction
 Multilayer shallow water model
 Semiimplicit staggered finite element discretization
 zsurfaceadaptive layers
 Advection with a spatially variable number of layers
 Numerical tests
 Conclusions
 Appendix A: Tracer constancy
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Multilayer shallow water model
 Semiimplicit staggered finite element discretization
 zsurfaceadaptive layers
 Advection with a spatially variable number of layers
 Numerical tests
 Conclusions
 Appendix A: Tracer constancy
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References