Multilayer shallow water models with locally variable number of layersand semi-implicit time discretization

# Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization

Luca Bonaventura , Enrique D. Fernández-Nieto ,
###### Abstract

We propose an extension of the discretization approaches for multilayer shallow water models, aimed at making them more flexible and efficient for realistic applications to coastal flows. A novel discretization approach is proposed, in which the number of vertical layers and their distribution are allowed to change in different regions of the computational domain. Furthermore, semi-implicit schemes are employed for the time discretization, leading to a significant efficiency improvement for subcritical regimes. We show that, in the typical regimes in which the application of multilayer shallow water models is justified, the resulting discretization does not introduce any major spurious feature and allows again to reduce substantially the computational cost in areas with complex bathymetry. As an example of the potential of the proposed technique, an application to a sediment transport problem is presented, showing a remarkable improvement with respect to standard discretization approaches.

MOX – Modelling and Scientific Computing,

Dipartimento di Matematica, Politecnico di Milano

Via Bonardi 9, 20133 Milano, Italy

luca.bonaventura@polimi.it

IMUS & Departamento de Matemática Aplicada I.

Avda. Reina Mercedes 2, 41012, Sevilla, Spain

edofer@us.es, jgarres@us.es, gnarbona@us.es

Keywords: Semi-implicit method, multilayer approach, depth-averaged model, mass exchange, sediment transport.

AMS Subject Classification: 35F31, 35L04, 65M06, 65N08, 76D33

## 1 Introduction

Multilayer shallow water models have been first proposed in [2] to account for the vertical structure in the simulation of large scale geophysical flows. They have been later extended and applied in [5], [4], [6]. This multilayer model was applied in [3] to study movable beds by adding an Exner equation. A different formulation, to which we will refer in this paper, was proposed in [32], which has several peculiarities with respect to previous multilayer models. The model proposed in [32] is derived from the weak form of the full Navier-Stokes system, by assuming a discontinuous profile of velocity, and the solution is obtained as a particular weak solution of the full Navier-Stokes system. The vertical velocity is computed in a postprocessing step based on the incompressibility condition, but accounting also for the mass transfer terms between the internal layers. In [22], this multilayer approach is applied to dry granular flows, for which an accurate approximation of the vertical flow structure is essential to approximate the velocity-pressure dependent viscosity.

Multilayer shallow water models can be seen as an alternative to more standard approaches for vertical discretizations, such as natural height coordinates, (also known as coordinates in the literature on numerical modelling of atmospheric and oceanic flows), employed e.g. in [10], [13], [15], terrain following coordinates (also known as coordinates in the literature), see e.g. [27], and isopycnal coordinates, see e.g. [8], [14]. Each technique has its own advantages and shortcomings, as highlighted in the discussions and reviews in [1], [10], [11], [28]. Multilayer approaches are appealing, because they share some of the advantages of coordinates, such as the absence of metric terms in the model equations, while not requiring special treatment of the lower boundary. On the other hand, multilayer approaches share one of the main disadvantages of coordinates, since they require, at least in the formulations employed so far, to use the same number of layers independently of the fluid depth. Furthermore, an implicit regularity assumption on the lower boundary is required, in order to avoid that too steeply inclined layers arise, which would contradict the fundamental hydrostatic assumption underlying the model.

In this work, we propose two concurrent strategies to make multilayer models more efficient and fully competitive with their and coordinates counterparts. On one hand, we propose a novel discretization approach, in which the number of vertical layers can vary over the computational domain. We show that, in the typical regimes in which the application of multilayer shallow water models is justified, the resulting discretization does not introduce significant errors and allows to reduce substantially the computational cost in areas with complex bathymetry. Thus making multilayer approach fully competitive with coordinate discretizations for large scale, hydrostatic flows. Furthermore, efficient semi-implicit discretizations are applied for the first time to this kind of models, allowing to achieve the same kind of computational gains that have been obtained for other vertical discretization approaches. In this paper, for simplicity, we have restricted our attention to constant density flows. An extension to variable density problems in the Boussinesq regime will be presented in a forthcoming paper.

In section 2, the equations defining the multilayer shallow water models of interest will be reviewed. In section 3, the spatial discretization is introduced in a simplified framework, showing how the number of layers can be allowed to vary over the computational domain. In section 4, some semi-implicit time discretizations are introduced for the model with a variable number of layers. Results of a number of numerical experiments are reported in section 5, showing the significant efficiency gains that can be achieved by combination of these two techniques. Some conclusions and perspectives for future work are presented in section 6.

## 2 Multilayer shallow water models

We consider the multilayer shallow water model described pictorially in Figure 1. In this approach, subdivisions of the domain are introduced in the vertical direction. We denote by the height of the layer and by the total height. Note that and that each subdomain is delimited by time dependent interfaces that are assumed to be represented by the one valued functions . For a function and for , we also define, as in [32],

 f−α+12:=(f|Ωα(t))|Γα+12(t)\;\;and\;\;f+α+12:=(f|Ωα+1(t))|Γα+12(t).

Obviously, if the function is continuous,

 fα+12:=f|Γα+12(t)=f+α+12=f−α+12.

Following [22], [32], the equations describing this multilayer approach can be written for as

 ∂thα+∇x⋅(hαuα)=Gα+12−Gα−12,∂t(hαuα)+∇x⋅(hαuα⊗uα)++ghα∇x(b+h)=1ρ0(Kα−12−Kα+12)++12Gα+12(uα+1+uα)−12Gα−12(uα+uα−1). (1)

Here, we consider a fluid with constant density is the velocity in the layer is a function describing the bathymetry (which is assumed to be constant in time, so that ) and the terms model the shear stresses between the layers. Notice that the atmospheric pressure has assumed to be zero. The vertical velocity profile is recovered from the integrated incompressibility condition, obtaining for and ,

 wα(t,x,z)=w+α−12(t,x)−(z−zα−12)∇x⋅uα(t,x), (2)

where

 w+α+12=(uα+1−uα)∇xzα+12+w−α+12,

and

 w−α+12=w+α−12−hα∇x⋅uαwithw+12=u1∇xb−G12.

Since we are focusing in this work mostly on subcritical flows, there is no special reason to choose discharge rather than velocity as a model variable. Therefore, we rewrite the previous system as

 ∂thα+∇x⋅(hαuα)=Gα+12−Gα−12,hα∂tuα+hαuα⋅∇uα+ghα∇x(b+h)==1ρ0(Kα−12−Kα+12)+(Gα+12Δ~uα+12+Gα−12Δ~uα−12), (3)

where From the derivation in [32], it follows that

 Gα+12 = ∂tzα+12+uα+1⋅∇xzα+12−w+α+12=∂tzα+12+uα⋅∇xzα+12−w−α+12, Kα+12 = −μα+12UHZα+12, (4)

where denotes the dynamic viscosity and is an approximation of at . We then define the vertical partition of the domain, setting where, for are positive constants such that

 N∑α=1lα=1.

Note that model (1) consists of equations in the unknowns

 h,{uα}α=1,…,N,    {Gα+12}α=1,…,N−1.

However, the mass transfer terms can be rewritten as

 Gα+12=∂tzα+12+uα+uα+12∂xzα+12−wα+12,wherewα+12=w+α+12+w−α+122.

As a consequence, the system has unknowns, now corresponding to the total height , the velocity in each layer and the averaged vertical velocity at each internal interface . By combining the continuity equations, the system can be rewritten with equations and unknowns. The unknowns of the reduced system are the total height and the velocity in each layer, for Note that can be written, by summing the mass equations from 1 to , as

 Gα+12=G12+α∑β=1(∂thβ+∇x⋅(hβuβ)). (5)

Moreover, for the special case and , the above equation leads to

 ∂th+∇x⋅(hN∑β=1lβuβ)=0.

By introducing this in the mass equation we obtain

 Gα+12= α∑β=1(∇x⋅(hlβuβ)−lβN∑γ=1∇x⋅(lγhuγ)). (6)

Assuming also system (3)-(2) is finally re-written as

 ∂tη + ∇x⋅(hN∑β=1lβuβ)=0, ∂tuα + uα⋅∇uα+g∇xη= =Kα−12−Kα+12ρ0hα+Gα+12Δ~uα+12+Gα−12Δ~uα−12hα,

for Here, we have set as customary in the literature . The transport equation for a passive scalar can be coupled to the previous continuity and momentum equation, in such a way as to guarantee compatibility with the continuity equation in the sense of [26]. If denotes the average density of the passive scalar in , it verifies the following tracer equation:

 ∂t(ραhα)+∇x⋅(ραhαuα)=ρα+1/2Gα+12−ρα−1/2Gα−12,

where

 ρα+1/2=ρα+ρα+12+12sgn% (Gα+12)(ρα+1−ρα).

In principle, any appropriate turbulence and friction model can be considered to define the turbulent fluxes , . Here, we have employed a parabolic turbulent viscosity profile and friction coefficients derived from a logarithmic wall law:

 ν=μρ0=κu∗(z−b)(1−z−bh),

where is the von Karman constant, is the friction velocity and denotes the shear stress. In order to approximate this turbulence model we set for :

 Kα+12=μα+12uα+1−uα(hα+hα+1)/2,withμα+12=ρ0κu∗α+12⎛⎝α∑β=1lβh⎞⎠(N∑γ=α+1lγ).

Trivially, . For and , standard quadratic models for bottom and wind stress are considered. We then set

 K1/2=−Cf∥u1∥u1,KN+1/2=−Cw∥uw−uN∥(uw−uN),

where denotes the wind velocity and the friction coefficient between at the free surface. The friction coefficient is defined, according to the derivation in [19], as:

 Cf= κ2(1−Δzrh)(ln(ΔzrΔz0))2, (8)

where is the roughness length and is the length scale for the bottom layer. Under the assumption that it can be seen that

 utu∗≈1κln(z−bΔz0)

where is the tangential velocity. In practice, we identify with , the horizontal velocity of the layer closest to the bottom, in the multilayer model. The definition of given by equation (8) is deduced by using previous relation of the ratio between and (see [19]). Then, we set

 u∗α+12=u1κln(α∑β=1lβh/Δz0),

in the definition of .

## 3 Spatial discretization with variable number of layers

The multilayer shallow water model (2) can be discretized in principle with any spatial discretization approach. For simplicity, we present the proposed discretization approach in the framework of simple finite volume/finite difference discretization on a staggered Cartesian mesh with C-grid staggering. A discussion of the advantages of this approach for large scale geophysical models can be found in [21]. The C-grid staggering also has the side benefit of providing a more compact structure for the system of equations that is obtained when a semi-implicit method is applied for time discretization. In order to further simplify the presentation, we only introduce the discretization for an vertical slice, even though any of the methods presented in the following can be easily generalized to the full three dimensional case. Generalization to structured and unstructured meshes can be obtained e.g. by the approaches proposed in [15] and [12], [17], [18], respectively, but higher order methods such as those of [34], [35] could also be applied. It is to be remarked that the choice of a staggered mesh is by no means necessary and that the approach proposed below to handle a variable number of layers can be easily extended to colocated meshes as well.

The solution domain will then coincide with an interval that is assumed to be subdivided into control volumes The step in the mesh is defined by and , where as usual. The discrete free surface variables are defined at integer locations corresponding to the centers of the control volumes, while the discrete velocities are defined at the half-integer locations As suggested in [26], the value of is taken to be that of the control volume located upwind of the volume edge.

The vertical number of layers employed, in the approach proposed in [32], is a discretization parameter whose choice depends on the desired accuracy in the approximation of the vertical structure of the flow. In order to make this type of model more flexible and more efficient, we propose to allow for a number of vertical layers that is not constant throughout the domain. The transition between regions with different numbers of layers is assumed to take place at the center of a control volume so that one may have different for and as a consequence, the discrete layer thickness coefficients are also defined at the half-integer locations The number of layers considered at the cell center for the purpose of the discretization of the tracer equation are defined as and the discrete layer thickness coefficients at integer locations are taken to be equal to those at the neighbouring half-integer location with larger number of layers. We will also assume that, whenever for some one has, without loss of generality, then for any there exist

 1≤α−i−12(β)≤α+i−12(β)≤Ni−12   such that   lβ,i+12=α+i−12(β)∑α=α−i−12(β)lα,i−12. (9)

This allows a more straightforward implementation of the numerical approximation of horizontal advection in the velocity and in the tracer equation, which are the only ones involving a horizontal stencil. Finally, again for simplicity of the implementation and without great loss of generality, it is assumed that if one has as well as

A sample configuration of this kind is depicted in figure 2. Notice that a dependence of the number of layers on time could also be introduced, in order to adapt the global maximum number of layers to the flow conditions, but this has not been done in the present implementation.

## 4 Semi-implicit time discretizations

The previous definitions yield a space discretization that can be easily coupled to any time discretization that yields a stable fully discrete space time scheme. For example, a time discretization by a third order Runge Kutta scheme has been employed as a reference in the numerical tests presented in section 5. However, we will focus here on semi-implicit time discretization approaches aimed at reducing the computational cost in subcritical regime simulations.

With this goal, it is immediate to notice that the formal structure of system (2) is entirely analogous to that of the three dimensional hydrostatic system considered in [15], [16], so that we can build semi-implicit time discretizations along the same lines, i.e. by treating implicitly the velocity in the continuity equation and the free surface gradient in the momentum equation. In the following, we present first a more conventional time discretization based on the off-centered trapezoidal rule (or -method, see e.g. [31]) and then a more advanced Implicit-Explicit Additive Runge Kutta second order method (IMEX-ARK2).

### 4.1 A θ-method time discretization

Following [15], we first consider a semi-implicit discretization based on the -method, which can be defined for a generic ODE system as

 yn+1=yn+Δt[θf(yn+1,tn+1)+(1−θ)f(yn,tn)],

where denotes the time step and is a implicitness parameter. If the method is unconditionally stable and the numerical diffusion introduced by the method increases when increasing . For the second order Crank-Nicolson method is obtained. In practical applications, is usually chosen just slightly larger than in order to allow for some damping of the fastest linear modes and nonlinear effects. We then proceed to describe the time discretization of system (2) based on the -method.

For control volume the continuity equation in (2) is then discretized as

 Δxiηn+1i+θΔt⎛⎜ ⎜⎝Ni+12∑β=1lβ,i+12hni+12un+1β,i+12−Ni−12∑β=1lβ,i−12hni−12un+1β,i−12⎞⎟ ⎟⎠ (10) =Δxiηni−(1−θ)Δt⎛⎜ ⎜⎝Ni+12∑β=1lβ,i+12hni+12unβ,i+12−Ni−12∑β=1lβ,i−12hni−12unβ,i−12⎞⎟ ⎟⎠.

It can be noticed that the dependency on has been frozen at time level in order to avoid solving a nonlinear system at each timestep. As shown in [15], [34], this does not degrade the accuracy of the method, even in the case of a full second order discretization is employed. For nodes the momentum equations for in (2) are then discretized as

 un+1α,i+12+gθΔtΔxi+12(ηn+1i+1−ηn+1i) (11) −Δtθlα,i+12hni+12⎛⎜ ⎜⎝νnα+12,i+12un+1α+1,i+12−un+1α,i+12lα+12,i+12hni+12−νnα−12,i+12un+1α,i+12−un+1α−1,i+12lα−12,i+12hni+12⎞⎟ ⎟⎠ =unα,i+12+ΔtAu,nα,i+12−g(1−θ)ΔtΔxi+12(ηni+1−ηni) +Δt(1−θ)lα,i+12hni+12⎛⎜⎝νnα+12,i+12unα+1,i+12−unα,i+12lα+12,i+12hni+12−νnα−12,i+12unα,i+12−unα−1,i+12lα−12,i+12hni+12⎞⎟⎠ +ΔtΔxi+12lα,i+12hni+12(Δ˜unα+12,i+12Gnα+12,i+12+Δ˜unα−12,i+12Gnα−12,i+12),

where denotes a discretization of the mass transfer term and denotes some spatial discretization of the velocity advection term. In the present implementation, an upstream based second order scheme is employed for this term. Notice that, to define this advection term, velocity values from different layers may have to be employed, if some of the neighbouring volumes has a number of layers different from that at For example, assuming again without loss of generality and and using the notation in (9), values

 u∗β,i−12=1lβ,i+12α+i−12(β)∑α=α−i−12(β)lα,i−12unα,i−12

will be used to compute the approximation of the velocity gradient at . Clearly, this may result in a local loss of accuracy, but the numerical results reported show that this has limited impact on the overall accuracy of the proposed method.

The discretization of the mass transfer term is defined as

 Gnα+12,i+12 = α∑β=1lβ,i+12⎛⎝(huβ−N∑γ=1lγhuγ)ni+1−(huβ−N∑γ=1lγhuγ)ni⎞⎠,

where is the upwind value depending on the averaged velocity . For the tracer equation, this term appears at the center of the control volume, so that we set instead

 Gnα+12,i = 1Δxiα∑β=1(lβ,i+12hni+12unβ,i+12−lβ,i−12hni−12unβ,i−12 (12) − lβ,iNi∑γ=1(lγ,i+12hni+12unγ,i+12−lγ,i−12hni−12unγ,i−12)).

The above formulas are to be modified appropriately for cells in which by summing all the contributions on the cell boundary with more layers that correspond to a given term on the cell boundary with less layers, according to the definitions in the previous section.

###### Remark 4.1

The time discretization of the mass transfer terms could be easily turned into an implicit one, by taking instead

 un+1α,i+12+gθΔtΔxi+12(ηn+1i+1−ηn+1i) (13) −Δtθlα,i+12hni+12(γnα+12,i+12(un+1α+1,i+12−un+1α,i+12)−δnα−12,i+12(un+1α,i+12−un+1α−1,i+12)) =unα,i+12+ΔtAu,nα,i+12−g(1−θ)ΔtΔxi+12(ηni+1−ηni) +Δt(1−θ)lα,i+12hni+12(γnα+12,i+12(unα+1,i+12−unα,i+12)−δnα−12,i+12(unα,i+12−unα−1,i+12)),

where now

 γnα+12,i+12=νnα+12,i+12lα+12,i+12hni+12+Gnα+12,i+122Δxi+12   δnα−12,i+12=νnα−12,i+12lα−12,i+12hni+12−Gnα−12,i+122Δxi+12.

This approach might be helpful to relax stability restrictions if large values of arise. In the implementation employed to obtain the numerical results of section 5, however, only the discretization (11) was applied so far.

At the bottom () and at the free surface () layers, the viscous terms are modified by the friction and drag terms at and respectively. We have then

 un+11,i+12+gθΔtΔxi+12(ηn+1i+1−ηn+1i) (14) −Δtθl1,i+12hni+12⎛⎜ ⎜⎝νn1+12,i+12un+12,i+12−un+11,i+12l1+12,i+12hni+12−Cnf,i+12∣∣∣un1,i+12∣∣∣un+11,i+12⎞⎟ ⎟⎠ =un1,i+12+ΔtAu,n1,i+12+ΔtΔxi+12l1,i+12hni+12Δ˜un32,i+12Gn32,i+12−g(1−θ)ΔtΔxi+12(ηni+1−ηni) +Δt(1−θ)l1,i+12hni+12⎛⎜⎝νn1+12,i+12un2,i+12−un1,i+12l1+12,i+12hni+12−Cnf,i+12∣∣∣un1,i+12∣∣∣un1,i+12⎞⎟⎠,
 un+1Ni+12,i+12+gθΔtΔxi+12(ηn+1i+1−ηn+1i) (15) +ΔtθlNi+12,i+12hni+12⎛⎜ ⎜ ⎜ ⎜⎝νnNi+12−12,i+12un+1Ni+12,i+12−un+1Ni+12−1,i+12lNi+12−12,i+12hni+12+˜Cwn,i+12un+1Ni+12⎞⎟ ⎟ ⎟ ⎟⎠ =unNi+12,i+12+ΔtAu,nNi+12,i+12+Δt˜Cwn,i+12lNi+12,i+12