Articles | Volume 16, issue 22
Development and technical paper
28 Nov 2023
Development and technical paper |  | 28 Nov 2023

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

Luca Arpaia, Christian Ferrarin, Marco Bajo, and Georg Umgiesser

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, Backhaus1985; Tidal, Residual, Intertidal Mudflat-3D, TRIM-3D, Cheng et al.1993; UNTRIM-3D, Casulli and Walters2000; Shallow water HYdrodynamic Finite Element Model, SHYFEM, Umgiesser2022). They are attractive when simulating strongly stratified flows (Hordoir et al.2015) and low-frequency motions (Leclair and Madec2011). 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/σ 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 Campin2004).

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,, 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 (x,z)R3 with x=(x,y)R2 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 zb(x) is the bathymetry, which does not depend on time. The water depth is H(x,t)=ζ(x,t)+zb(x). 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 α=1,,N, ordered from the free surface to the bottom. The layers are non-overlapping with Ω=α=1NΩα. Each layer Ωα is delimited laterally by the vertical domain boundary and in the vertical by the time-dependent interfaces Γα±1/2(t) defined by the set of points of coordinates (x,z) such that z=zα±1/2(x,t). The free surface Γζ and the bottom interfaces Γb are described by the free surface elevation z1/2=ζ(x,t) and by the bathymetry function zN+1/2=-zb(x) 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 (x,s)R3 such that:


and discretized by means of a vertical grid similarly composed of N layers, each denoted Ωα0. The reference layers are delimited vertically by the fixed-in-time interfaces Γα±1/20, which are placed at the vertical coordinate given by the coefficients sα±1/2. Such constants can be ordered:


Then, the interface position can be obtained by mapping the reference interface Γα-1/20 to the actual or physical interface Γα-1/2(t). In general we assume that there exists a function, for α=1,,N:

(1) A : Γ α - 1 / 2 0 Γ α - 1 / 2 ( t ) , z α - 1 / 2 = A ( x , s α - 1 / 2 , t ) x Ω x .

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

(2) z α - 1 / 2 = ζ ( x , t ) + s α - 1 / 2 ζ ( x , t ) + z b ( x ) ,

which assures a surface- and terrain-following grid that is limited by the interfaces Γ1/2(t)=Γζ(t) and ΓN+1/2=Γb. The reference and the physical domains with their vertical subdivisions are sketched in Fig. 1.

Figure 1One-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 α=1,,N:


where the coefficients lα=sα-1/2-sα+1/2 are prescribed after the creation of the reference grid. They are positive and they sum to one α=1Nlα=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 α=1,,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 Poisson1981) at one atmosphere and at a constant potential temperature of 12.4C. If the equation of state is of type ρ=ρ(T), the density vertical discretization derives from the tracer one, for α=1,,N:

(7) ρ α ( x , t ) = ρ ( T α ( x , t ) ) .

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):

  • If the function is continuous

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


Mass conservation reads

(8) ζ t + β = 1 N h β u β = 0 .

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 α=1,,N:


where Gα±1/2 is the mass-transfer function responsible for the vertical mass exchange between the layers, Kα±1/2 are the vertical viscous fluxes that model the shear stress between the layers, KT,α±1/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 Petersen1997; 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

(11) G α - 1 / 2 = z α - 1 / 2 u α + σ α - 1 / 2 - w α - 1 / 2 + = z α - 1 / 2 u α - 1 + σ α - 1 / 2 - w α - 1 / 2 - ,

with σα-1/2 the velocity of the grid interface:

(12) σ α - 1 / 2 = z α - 1 / 2 t ,

and wα-1/2± the vertical fluid velocity at the interface. The vertical velocity is computed from the following relationships:

(13) w α - 1 / 2 + = - w α + 1 / 2 - - h α u α and w α - 1 / 2 - = w α - 1 / 2 + + z α - 1 / 2 u α - u α - 1 ,

which are evaluated starting from the bottom α=N,,1, where the no-slip condition is imposed wN+1/2-=uNzb. In practice, and as is standard in ocean models, the mass-transfer function is computed directly from the layer-integrated mass equation

(14) G α - 1 / 2 = G α + 1 / 2 + h α t + h α u α .

Summing from N to α as

(15) G α - 1 / 2 = G N + 1 / 2 + β = N α h β t + β = N α h β u β ,

which implies G1/2=0 or no mass loss at the free surface. The vertical velocity at the interfaces wα-1/2± 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 Γα-1/2 it reads


with Gα-1/2+=max(0,Gα-1/2) and Gα-1/2-=min(0,Gα-1/2). For the tracer, a total variation diminishing (TVD) flux is employed (LeVeque2002).

The terms Kα-1/2 and KT,α-1/2 are the vertical viscous and diffusive fluxes computed at the interface Γα-1/2:


where να-1/2 is the vertical viscosity and νT,α-1/2 the vertical diffusivity. Dz(⋅) 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 CF 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 layer-integrated pressure gradient term zα+1/2zα-1/2p(z)dz, 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β-1/2 is the distance between the layer centers, that is, hβ-1/2=(hβ-1+hβ)/2 for β=2,,N and hβ-1/2=h1/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β=gρ0-ρβρ0 is the layer buoyancy, the buoyancy at the interface is resolved with an average bβ-1/2=12bβ-1+bβ for β=2,N and bβ-1/2=b1 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 Dz(bβ-1/2)=0 for β=1 and Dz(bβ-1/2)=(bβ-1-bβ)/hβ-1/2 for β=2,,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).


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 Campin2004). The reference domain, with vertical coordinate Z, is


This domain is discretized by means of a vertical grid composed of N layers, with interfaces Γα-1/20, 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 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 Γα-1/20 onto the horizontal plane as

(16) Ω x , α = { x : x Ω x and - z b ( x ) Z α - 1 / 2 } .

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


Each layer Ωα0 is delimited on the upper and bottom sides by Γα1/20 and laterally by the vertical domain boundary; it could also be delimited by Γαb (see Fig. 2, right panel). To map the reference interface Γα-1/20 to the physical interface Γα-1/2, again, we can use a generalized coordinate transformation, for α=1,,N:

(17) z α - 1 / 2 = ζ ( x , t ) + S α - 1 / 2 ( x ) ζ ( x , t ) + z b ( x ) , x Ω x , α ,

with Sα-1/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

(18) h α ( x , t ) = z α - 1 / 2 ( x , t ) - max z α + 1 / 2 ( x , t ) , - z b ( x ) = S α - 1 / 2 ( x ) - max S α + 1 / 2 ( x ) , - 1 H ( x , t ) = l α ( x ) H ( x , t ) , x Ω x , α .

If we define ΔZα(x)=Zα-1/2-maxZα+1/2,-zb(x) we can rewrite the coefficients, for α=1,N:


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

Figure 2One-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α-1/2=-zb(x). 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/σ layers (Burchard and Petersen1997), 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 Nb(x)={α:Zα+1/2<-zb(x)Zα-1/2}.


For z-layers 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 mass-transfer function (Eq. 14) coincides with the vertical velocity:


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∈ℝ2 composed of non-overlapping triangular elements. We denote the horizontal grid by 𝒯 with K∈𝒯 the generic triangle, and |K| its area. The local reference element length is hK 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=1,2,3 or jK). Given a node i in an element K, niK 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 Ci 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 Ci. The edge of Ci separating CiK and CjK has an exterior normal called nijK, 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

(19) j K , j i n i j K = - n i K 2 .

After the horizontal discretization, the domain is subdivided into prismatic boxes K×[zα+1/2,zα-1/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 NK. For a staggered discretization it is also helpful to define a nodal bottom layer Ni=maxKDiNK. 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


Figure 3Grid and notation. (a) Triangle K with nodes and scaled normals. (b) Set 𝒟i with dual cell area Ci and dual cell boundary Ci. 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 Zienkiewicz1981). 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 α=1,,NK, being the elemental discharge values and with Tα,i(t), defined for α=1,,Ni, 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

(22) ζ ( x , t ) = i T φ i ( x ) ζ i ( t ) ,

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 uα,K=qα,Khα,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):


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 first-order upwind flux qαuα^ (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 AKdRNK×NK. 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 Δu=un+1-un, then:


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


with AK=I|K|-ΔtAKd a tridiagonal, positive definite, and diagonally dominant matrix. The nonlinear dependence of the external pressure gradient term from HK has been resolved by using the old value. Also, the viscous matrix has been computed with frozen values at tn. In FKn all the quantities are computed at tn, including the mass-transfer function. These choices avoid solving a nonlinear system at each time step. The variation Δ()*=()*-()n is the solution of the following Euler step with an explicit external pressure gradient:


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

(27) K D i j K m i j K + g θ 2 Δ t 2 a i K x 1 T A K - 1 H K n a j K x + a i K y 1 T A K - 1 H K n a j K y Δ ζ j = Δ t K D i a i K x 1 T θ Δ U K * + U K n + a i K y 1 T θ Δ V K * + V K n ,

where mijK=Kφiφjdx is the Galerkin mass matrix based on the piecewise linear Lagrange basis functions. The Galerkin mass matrix, in SHYFEM, is lumped. The vector 1RNK 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 GNi+1/2,in+1=0, for α=Ni,,1:

(28) | C α , i | G α - 1 / 2 , i n + 1 = | C α + 1 , i | G α + 1 / 2 , i n + 1 + | C α , i | Δ h α , i Δ t - K D α i ( a i K x q α , K x , n + θ + a i K y q α , K y , n + θ ) .

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 G1/2,in+1=0. With standard z-layers, the contribution related to the grid velocity is zero Δhα,i=Δt[σi]α+1/2α-1/2=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:


where Tαqα^ is an appropriate numerical tracer flux across the dual cell boundary. At the lateral boundary Ωx,α, 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 H^α(Tα,i,Tα,j) being the numerical flux in the horizontal direction and Tα+1/2,iGα+1/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:


Vertical diffusion can also be assembled in matrix form through the discrete matrix AidRNi×Ni. Then, the discretization of the layerwise tracer Eq. (10) read

(31) A i T i n + 1 = Diag { | C α , i | h α , i n } T i n + Δ t F i n ,

with Ai=Diag{|Cα,i|hα,in+1}-ΔtAid 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 qαn+θ whereas the vertical numerical flux uses the last available mass-transfer function Gα±1/2n+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).

4z-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 Δttn,tn+1 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 z̃α-1/2,in+1=z̃α-1/2(xi,tn+1) is


where zα-1/2,in+1=zα-1/2(xi,tn+1) is the interface position after the grid movement and Δz̃α-1/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




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 tn+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 tn and tn+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.

Figure 4Grid and tracer evolution during one time step. The process is interpreted as four stages, which bring from the pair (Tn,ζhn) to (T̃n+1,ζn+1). The vector T={T1,T2} 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 ζ0(x), we define a set of active indices and the surface layer index, by element, as

(32) α active , K = { α α K : Z α + 1 / 2 + ϵ top < min x K ζ 0 ( x ) } , α top , K = min α active , K ,

with αK={1,,NK}. Owing to the staggering of the grid, it is convenient to define also at each node:

(33) α active , i = { α α i : Z α + 1 / 2 + ϵ top < ζ i 0 } , α top , i = min α active , i ,

with αi={1,,Ni}. 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 5This 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:

    (34) α mov , i = { α α i : Z α - 1 / 2 + ϵ mov > ζ i n + 1 } .

    ϵ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 Nmov,i=maxαmov,i. The depth of the moving layers is

  • 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:

    (35) z α + 1 / 2 , i n + 1 = ζ i n + 1 + S α + 1 / 2 , i ζ i n + 1 + z mov , i ,

    with Sα+1/2,i a stretching function. Then, the nodal layer thickness reads:

    (36) h α , i n + 1 = l α , i ζ i n + 1 + z mov , i , α = α top , i , , N mov , i .

    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α,i=ΔZαzmov,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 ζn+1. We get the new top layer indices αtop,Kn+1 and αtop,in+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 T̃αn+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αn+1. The remapped value is the solution of the following advection equation integrated over the layer thickness:

(37) h ̃ α T ̃ α t = [ σ top T ̃ ] α + 1 / 2 α - 1 / 2 ,

with an upwind flux:

(38) σ α - 1 / 2 top T α - 1 / 2 n + 1 = σ α - 1 / 2 top + T α n + 1 + σ α - 1 / 2 top - T α - 1 n + 1 .

We consider the discrete case. After integration on the dual cell and with a simple forward Euler time stepping (with initial condition Tαn+1) we have

(39) h ̃ α , i n + 1 T ̃ α , i n + 1 = h α , i n + 1 T α , i n + 1 + Δ t σ α - 1 / 2 , i top T α - 1 / 2 , i n + 1 - σ α + 1 / 2 , i top T α + 1 / 2 , i n + 1 ,

with the new nodal layer thickness:


In the case of an element removal (αtop,in+1>αtop,in), we identify the layer that should disappear and we proceed with a collapse of the lower interface to the upper one. For α=αtop,in,,αtop,in+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 (αtop,in+1<αtop,in), we identify the layer that should appear and we expand the interface. Then the remap for α=αtop,in+1,,αtop,in 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=rmovΔZα in the relationship Eq. (34). In order to obtain the z-surface-adaptive grid we have chosen rmovrtop; in practice, we have set rmov=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 rmov. We define

(40) R α = ζ max - Z α - 1 / 2 Δ Z α ,

with ζmaxmaxx,tζ(x,t) an estimate of the maximum free surface height during the simulation. We get

  • z-star if rmovRN 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 rmov=RM 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 rmov we compare the newly developed z-surface adaptive layers against z-star.

Figure 6The 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 tn and tn+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)).

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


Figure 8Treatment of nonconformal box for the one-dimensional case. (a) Splitting with fictitious layers. (b) The mass-transfer function G1+1/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 αtop,i+1/2>αtop,i-1/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:


with αmin,K=minjKαtop,j. Then the boxes called hereinafter for simplicity “nonconformal” can be identified by the pair of indices αtop,K,K*. 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α,K* for the fictitious layers. Then we distribute the discharge of the top layer among the fictitious layers, for α=αmin,K,,αtop,K:

(41) q α , K * = l α , K * q α top , K , K ,

with lα,K*=hα,K*hαtop,KK. 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 (αtop,K,K) with node iK, as illustrated in one dimension in Fig. 7. After the splitting Eq. (41), the mass-flux term (only the x-component shown) reads, for α=αtop,i,,αtop,K:

(42) K φ i x q α * d x = a i K c α , i * q α top , K , K ,


(43) c α , i * = β = α top , i α min , K l β , K * if α = α top , i and α min , K < α top , i l α , K * otherwise (hanging layer) ,

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α-1/2,in+1 for α=αtop,i,α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 α=αtop,i,,αtop,K:

(44) H ^ α = β = α top , i α min , K l β , K * H ^ α top , K ( T β * , i , T β * , j ) if α = α top , i and α min , K < α top , i l α , K * H ^ α top , K ( T α * , i , T α * , j ) otherwise (hanging layer) .

Again, we have separated the cases of a node with or without hanging layers. Note that the subscript α*=max(α,α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 tn+1 by computing the following relative volume error for the dual cell area, which results from the sum of Eq. (28) from Ni to αtop,i:


To quantify the tracer constancy error, we use the L1norm:


with T0 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 [-5m,5m]×[-5m,5m] with constant depth zb=1m. The basin is initially at rest and the free surface is perturbed by the following Gaussian hump:


with A=0.5m, τ=0.5m2 and r=x2+y2. 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 hK=0.25m. We compare different vertical resolutions with variable layer thickness. The coarsest grid has three layers: a first surface layer with a thickness of ΔZ1=0.2m, the second and the third layers have thicknesses of ΔZ2=ΔZ3=0.4m. 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 ΔZ1=0.025m.

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 9Impulsive 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α+1/2 are traced with dashed lines. In the regions where the free surface crosses the interface Zα+1/2 it means that the layer α locally has been removed from the computation.


Figure 10Impulsive 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=25km. The bathymetry is represented by the following function zb(x)=H0/Lx and H0=10m. The horizontal element size is uniform and equal to hK=250m. A periodic water level is imposed at the seaward boundary as ζ(t)=A(1-sin2πtT) with amplitude A=10m, period T=1d and the time t ranging from 0 to 0.5d. 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 111 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α+1/2. In the regions where the free surface crosses the interface Zα+1/2 it means that the layer α locally has been removed from the computation.


Figure 121 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 CF=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.4m and period 12 h. The salinity at the sea boundary is constant and fixed to 38PSU. A weak freshwater flow with a discharge of 500m3 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 38PSU. 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 hK=2km at the sea boundary, to hK=100m in the inner shelf close to the lagoons and river branches, and to hK=50m 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=-1m) is equal for the two grids and it is composed of 23 levels with variable thicknesses from ΔZ=0.5m near the surface up to ΔZN=4m for the last layer. The resolution of the upper part of the water column differs: the 24-layer grid has one layer with ΔZ1=1m. This choice avoids the drying of the first layer. The 27-layer grid, in the upper part, has four layers with constant thickness, ΔZ1=ΔZ2=ΔZ3=ΔZ4=0.25m. 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.4m, 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 13Po river. (a) Horizontal grid. (b) Zoom of the horizontal grid with tidal stations and the transect in the Pila branch.

Figure 14Po 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 15Po 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 16Po 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 17Salinity 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 Ti=1 in the discrete tracer Eq. (31). Each row reduces to




Using, first, the numerical flux consistency H^α1,1=qαn+θnijK 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 niK2=|K|φiK. 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 mass-fluxes. According to Eq. (44), the horizontal tracer fluxes in nonconformal boxes should be computed with

H^α=β=αtop,iαmin,Klβ,K*H^αtop,K(Tβ*,i,Tβ*,j)ifα=αtop,iandαmin,K<αtop,ilα,K*H^αtop,K(Tα*,i,Tα*,j)otherwise (hanging layer).

For a constant tracer, it can be rewritten for α=αtop,i,αtop,K:


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


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


which, thanks to the definition provided in Sect. 4.2 of grid velocity σα-1/2,itop=z̃α-1/2,in+1-zα-1/2,in+1Δt and layer thickness h̃α,in+1=z̃α-1/2,in+1-zα+1/2,in+1, is an identity.

Code and data availability

The SHYFEM hydrodynamic model is open source (GNU General Public License as published by the Free Software Foundation) and freely available through GitHub at (last access: 5 November 2023). The current developments have been implemented in a branch of the SHYFEM code that can be accessed from Zenodo (, Arpaia2023). Configuration files and data used to run each test case are also available at the same Zenodo repository.

Author contributions

LA: Conceptualization, Methodology, Software, Validation, Writing, Formal analysis. CF: Conceptualization, Methodology, Funding acquisition, Writing, Resources, Validation. MB: Methodology, Writing. GU: Conceptualization, Methodology, Writing, Software.

Competing interests

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 (, last access: 8 November 2023) developed at the CNR-ISMAR. 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.

Financial support

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 V-A Interreg Italy-Croatia CBC programme.

Review statement

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 free-surface flows in ocean circulation models, Ocean Model., 7, 269–284, 2004. a, b, c

Arpaia, L.: SHYFEM version with surface-adaptive z-coordinates (7_5_71-zlay1.2), Zenodo [code],, 2023. a

Audusse, E., Bristeau, M.-O., Pelanti, M., and Sainte-Marie, J.: Approximation of the hydrostatic Navier-Stokes 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 three-dimensional 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., Fernandez-Nieto, E. D., Garres-Diaz, J., and Narbona-Reina, G.: Multilayer shallow water models with locally variable number of layers and semi-implicit 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 z-coordinates 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 semi-implicit method for three-dimensional shallow water flow, Comput. Math. Appl., 27, 99–112, 1994. a

Casulli, V. and Cheng, R.: Semi-implicit finite difference methods for three-dimensional 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 three-dimensional ocean models, Ocean Model., 145, 101530,, 2020. a

Fernández-Nieto, E., Koné, E., and Rebollo, T. C.: A Multilayer Method for the Hydrostatic Navier-Stokes 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,, 2021. a

Griffies, S., Pacanowski, R., Schmidt, M., and Balaji, V.: Tracer conservation with an explicit free-surface method for z-coordinate 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 free-surface models, J. Comput. Phys., 38, 307–327, 2002. a, b

Guardone, A., Isola, D., and Quaranta, G.: Arbitrary Lagrangian Eulerian formulation for two-dimensional 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 z-level coordinate models, Ocean Dynam., 60, 1447–1461, 2010. a

Klingbeil, K., Lemarié, F., Debreu, L., and Burchard, H.: The numerics of hydrostatic structured-grid coastal ocean models: state of the art and future perspectives, Ocean Model., 125, 80–105, 2018.  a, b

Leclair, M. and Madec, G.: z-Coordinate, 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,, 2002. a

Lin, S. J. and Rood, R. B.: Multidimensional flux form semi-Lagrangian 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 978-3-642-08754-7, 2002. a, b, c

Millero, F. J. and Poisson, A.: International one-atmosphere equation of state of seawater, Deep-Sea 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, (last access: 6 April 2023), 2011. a

Shchepetkin, A. and McWilliams, J.: A method for computing horizontal pressure-gradient force in an oceanic model with a nonaligned vertical coordinate, J. Geophys. Res., 108,, 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, ISMAR-CNR 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, 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 shallow-water wave equations, Int. J. Numer. Meth. Fluids, 1, 81–97, 1981. a, b

Short summary
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 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 very thin surface layers.