Formulation, Implementation and Validation of the Horizontal Coupling Method for 1D/2D Shallow Water Flow Models

# Formulation, Implementation and Validation of the Horizontal Coupling Method for 1D/2D Shallow Water Flow Models

C. Nwiagwe nwaigwe@warwick.ac.uk Centre for Scientific Computing and Warwick Mathematics Institute, University of Warwick, United Kingdom A. S. Dedner a.s.dedner@warwick.ac.uk Centre for Scientific Computing and Warwick Mathematics Institute, University of Warwick, United Kingdom

Formulation, Implementation and Validation of the Horizontal Coupling Method for 1D/2D Shallow Water Flow Models

Received: never; final revision: never; published: never.

One dimensional (1D) simulations of the flow and flooding of open channels are known to be inaccurate as the flow is multi-dimensional in nature, especially at the flooded regions. However, multi-dimensional simulations, even in two dimensions (2D), are computationally expensive, hence the problem of efficiently coupling 2D and 1D simulations for the flow and flooding of open channels has been the subject of much research and is investigated in this paper. We adopt a 1D model with coupling term for the channel flow and the 2D shallow water flow model for the floodplain. The 1D model with coupling term is derived by integrating the 3D Free Surface Euler equations but without imposing any restriction on the channel width variations. Finite volume methods are formulated for both the 2D and 1D models including a discrete coupling term in closed form. Coupling is achieved through the discrete coupling term in the 1D model and the lateral numerical fluxes in the 2D model. Since the lateral discharge in the channel cannot be guaranteed to be zero during flooding, we aim to recover the lateral variation by computing two lateral discharges over each cross section and propose to use an ad-hoc model based on the -discharge equation in the 2D model for this purpose. We then propose the numerical scheme for this ad-hoc model following the hydrostatic reconstruction philosophy. Then, we show that the resulting method, named Horizontal Coupling Method (HCM), is well-balanced; we introduce the no-numerical flooding property and also show that the method satisfies the property. Three numerical test cases are used to verify the performance of the method. The results show that the method performs well in both accuracy and efficiency and also approximates the channel lateral discharges with very good accuracy and little computational overhead.

## 1 Introduction

Flows in open channels, such as rivers, in which the vertical and lateral variations in velocities can be assumed negligible, can be accurately simulated using the 1D Saint Venant Equations. During flooding, the channel overflows and the flow becomes high dimensional, rendering 1D simulation inadequate. These claims have also been observed numerically, see  for example. But then, even a 2D simulation of the entire flow is computationally expensive. This leads to the difficulty of choosing between an expensive but more accurate high dimensional simulations and an inexpensive but less accurate 1D simulation. To tackle this problem, a 1D simulation can be used along the channel while a 2D simulation is used for the floodplains. The problem of how to couple the two simulations then arises. This has led to many research work and also the subject of this paper.

A lot of research has been carried out to propose methods to couple 1D channel model with 2D floodplain flow model. In , a method, which numerically couples the 1D and 2D models by including the lateral numerical fluxes in the 1D numerical scheme, is proposed. They referred to the method as the Flux-Based Method (FBM). The theory of characteristics was employed to couple 1D/2D models in ; matching conditions are defined at the 2D/1D interfaces, then a prediction and correction algorithm was used to ensure that these conditions are satisfied. The 1D river model and the 2D non-inertia model were also coupled in  to simulate the interaction of a sewer system with over-land flow. Here, the water level differences between the flows in the two domains are used to calculate the interacting discharges in the sub-domains.

In , see also , two methods which are based on post-processing of separately computed solutions of the existing 1D and 2D models, are proposed. In these methods, the separately computed solutions are used to calculate the total water volume in a 1D cell and all its adjacent floodplain 2D cells. Then, a common water level is found for all the cells, finally the wetted cross sectional area for the 1D cell and the water height for all the adjacent 2D cells are found. These methods have been applied to Tiber River, Rome in . The superposition approach, proposed in , classically derives the exchange terms in the 1D model from the full 3D Inviscid Euler’s equations; an optimal control process is used to couple the models. The superposition approach of  was extended to finite volume methods in , proposing a discrete exchange term that leads to globally well-balanced scheme. This approach superposes a 2D grid over the 1D channel grid and convergence is achieved using a Schwartz-like iterative algorithm. For practical cases, the iterative algorithm can jeopardise the overall efficiency of the method .

A great difficulty for coupling methods is how to calculate the lateral discharge along the river channel because the 1D model does not have an equation to compute it. The channel lateral discharge is set to zero and used to calculate the 2D numerical fluxes at the 2D/1D interfaces in . In , the exchange terms derived in  were adopted and a strategy to estimate the lateral discharge without superposition or overlapping, was proposed. The approach is an iterative technique which uses the solution of successive Riemann problems to estimate the transverse velocity. Another approach which decomposes the channel 1D discharge into lateral and frontal components, using the angle which the channel axis makes with the -axis, can be found in [2, 14]. However, if this angle is zero, then this approach would be inadequate whenever the channel is full because it would always compute a zero discharge which is unrealistic. Therefore, the problem of computing the channel lateral or transverse discharge remains challenging.

In addition to the above difficulty, another fundamental issue is the 1D assumption on the channel flow, namely that both the free surface elevation and lateral velocity are laterally constant. By physical intuition, during overflow like flooding or draining, water flows out of or into the channel from both of its lateral boundaries. This means that the lateral velocities (or discharges) at both sides are in opposite directions and very likely to differ in magnitude. Therefore, the lateral discharge will rarely be constant across the channel cross sections, even when the free-surface elevation is assumed constant over the cross section. This means that the 1D assumption is inadequate if overflowing. However, most existing coupling methods retain this 1D assumption even during overflowing. We, therefore, propose that different discharges for each lateral boundary of a cross section, need to be computed; we propose to use the 2D -component shallow water equation, as an ad-hoc model, to compute these lateral discharges.

The method we propose here, which we the Horizontal Coupling Method (HCM), follows the lines of  to derive a similar but slightly different coupling terms however, we do not impose or use any restriction on the channel width variations. The essence of this paper is, therefore, to propose a strategy (i) to overcome the difficulty in calculating the lateral discharges, (ii) eliminate the limitations of the 1D assumption on the channel lateral discharge during flooding, (iii) to derive a more general variant of the coupling terms of  (iv) prove the properties of the resulting method and (v) validate the method using some numerical test cases.

The paper is organised as follows. In section 2, we present the background for the problem of the flow and flooding of open channels and derive the channel flow model with coupling term. We also present the 2D shallow water flow models for the floodplain flows. The numerical schemes for the uncoupled 1D and 2D models are presented in section 3; the numerical scheme for the lateral discharge model and the detailed formulation of the discrete coupling term are given in sections 3.2.2 and 3.2.3 respectively. The algorithm for the HCM is summarised in section 3.4; a flow chart for the implementation is also given. The properties of the method are discussed in section 4 where we introduce the no-numerical flooding property and show that the method preserves the property and is also well-balanced. In section 5, we present some numerical test cases to access the performance of the method and discuss the results; we conclude the paper in section 6.

## 2 Mathematical Models

### 2.1 Background Figure 1: Flow over a domain with bottom topography zb(x,y) (dashed line) comprising of a channel and floodplain. The channel length is along the x-axis and the width, along the y-axis; H(x,y,t) is water depth, η(x,y,t), the free surface elevation, t is time variable and (x,y,z)∈R3. (a) Flow cross section at a fixed point x. The bottom topography comprises of the channel with bank elevations, zbl(x) and zbr(x), and the floodplain which occupies the remaining regions.

In this section, we present the model equations for the problem under consideration. Let us begin by considering the flow of water over the fixed horizontal 2D domain, which consists of a channel and floodplains (see figure 1), such that the flow at time, occupies the 3D domain, defined by

 Ωt={(→X,z)∈R3:→X=(x,y)∈ΩH,zb(→X)≤z≤η(→X,t)}, (1)

bounded below by a fixed bottom, and above by the water free-surface elevation, given by

 η(→X,t)=zb(→X)+H(→X,t) (2)

and is the depth of fluid at time, . A cross section of the flow domain is shown in figure 2(a). The length of the channel lies along the -axis (frontal direction) and the width, along the -axis (lateral direction), while and are the left and right bank elevation of the channel, see figure 2(a). An important quantity is the maximum channel wall elevation or simply, channel wall elevation.

###### Definition 2.1 (Channel wall elevation, zwb(x))

The channel wall elevation at cross section x, is the minimum elevation of the channel banks above which flooding is said to have occurred. We denote it by , that is

 zwb(x)=min(zbl(x),zbr(x)), (3)

see figures 2(a) and 2(b).

That is, is the channel top. Figure 2(b) shows the channel cross section, depicting its geometry, including the 1D laterally constant channel bottom topography given by

 Zb(x)=minyzb(x,y).

It also depicts the wall elevation , the top width and the -coordinates, and of the left and right boundaries at the top, where gives the channel width at an elevation above a reference elevation , and and are the coordinates of the left and right lateral boundaries, respectively, at elevation . So that

 B(x,z)=yr(x,z)−yl(x,z)∀z, (4)

such that

 yr(x,z)=yl(x,z),B(x,z)=0 for all z

and the bottom elevation satisfies

 zb(→X)|y=yl(x,z),yr(x,z)=z∀z∈[Zb(x),zwb(x)], (6)

see figure 2(b). Furthermore, we extend the definition of the width functions above the channel top () in the following:

 yl,r(x,z)=ywl,r(x),B(x,z)=B(x,zwb(x))∀z≥zwb(x), (7)

see figure 2(b).

The flow cross section, in 2(a) has been partitioned into (i) the channel cross section, and (ii) the floodplains, and , see figure 2(b). The flow in the floodplains is simulated with the standard 2D shallow water models (see section 2.3), therefore we focus on deriving the model equations for the flow in the channel.

With the channel geometry completely defined, we now consider the initial flow condition. In general, the free surface elevation is 2D, see figures 1 and 2(a). However, we assume that it is always 1D (laterally constant and given by ) within the channel, see figure 2(b)). Hence, we say the channel is full whenever

 ¯η(x,t)>zwb(x). (8)

Note that in general, the channel flow lateral boundaries are at the coordinates, , not , see figures 2(b) and 3(a). These two sets of coordinates are only equal if the channel is full, see (7) and also figure 2(b). They are not equal if the channel is not full, see figure 3(a). If the channel is not full (), then the water height and velocities are zero at the top lateral boundaries, that is

 (u(→X,z,t),v(→X,z,t),w(→X,z,t),H(→X,t))∣∣∣y=ywl,r(x)=0 whenever ¯η(x,t)≤zwb(x), (9)

(see figure 3(a)) where are velocity components along the directions respectively.

### 2.2 Derivation of the channel flow model with coupling term

Under the assumption of compressible and inviscid fluid, the flow of water in the channel is governed by the following 3D Free-Surface Euler Equations :

 ∂xu(→X,z,t)+∂yv(→X,z,t)+∂zw(→X,z,t)=0. (10) ∂tu(→X,z,t)+u(→X,z,t)∂xu(→X,z,t)+v(→X,z,t)∂yu(→X,z,t)+w(→X,z,t)∂zu(→X,z,t)=−1ρ∂xP(→X,z,t). (11) ∂tv(→X,z,t)+u(→X,z,t)∂xv(→X,z,t)+v(→X,z,t)∂yv(→X,z,t)+w(→X,z,t)∂zv(→X,z,t)=−1ρ∂yP(→X,z,t). (12) ∂tw(→X,z,t)+u(→X,z,t)∂xw(→X,z,t)+v(→X,z,t)∂yw(→X,z,t)+w(→X,z,t)∂zw(→X,z,t)=−1ρ∂xP(→X,z,t)−g. (13) (u(→X,z,t)∂xzb(→X)+v(→X,z,t)∂yzb(→X)−w(→X,z,t))∣∣∣z=zb(→X)=0. (14) (∂tη(→X,t)+u(→X,z,t)∂xη(→X,t)+v(→X,z,t)∂yη(→X,t)−w(→X,z,t))∣∣∣z=η(→X,t)=0. (15)
 P(→X,z,t)=Patm on z=η(→X,t), (16)

where and are the fluid density, velocity vector and pressure at point at time, and is the atmospheric pressure, which is usually conveniently taken to be zero.

The flow quantities of interest in the 1D channel model are the wetted cross sectional area, and the section averaged discharge, given by the following averages:

 Q(x,t) =∫yr(x,¯η(x,t))yl(x,¯η(x,t))∫¯η(x,t)zb(→X)u(→X,z,t)dzdy, (17) A(x,t) =∫yr(x,¯η(x,t))yl(x,¯η(x,t))∫¯η(x,t)zb(→X)dzdy=∫yr(x,¯η(x,t))yl(x,¯η(x,t))H(→X,t)dy. (18)

So that the section-averaged velocity, is given as

 u––(x,t)=Q(x,t)A(x,t)=1A(x,t)∫yr(x,¯η(x,t))yl(x,¯η(x,t))∫¯η(x,t)zb(→X)u(→X,z,t)dzdy. (19)

First, we note that -independence of the free-surface, means that the sum,

 H(→X,t)+zb(→X)=¯η(x,t),∀yl(x,¯η(x,t))≤y≤yr(x,¯η(x,t)), (20)

is constant in , even though each of and depends on , see figures 2(b) and 3(a).

The shallow water assumption that water depth is small compared to horizontal length leads to neglect vertical acceleration, hence the z-momentum equation, (13) and the dynamic boundary condition, (16), lead to the following hydrostatic pressure:

 P(→X,z,t)=ρg(¯η(x,t)−z)⟹∂xP(→X,z,t)=ρg∂x¯η(x,t).

Therefore, the FSEE (10) -(16), reduce to the following system:

 ∂xu(→X,z,t)+∂yv(→X,z,t)+∂zw(→X,z,t)=0. (21) ∂tu(→X,z,t)+∂x(u2(→X,z,t))+∂y(u(→X,z,t)v(→X,z,t))+∂z(u(→X,z,t)w(→X,z,t))=−g∂x¯η(x,t). (22) (u(→X,z,t)∂xzb(→X)+v(→X,z,t)∂yzb(→X)−w(→X,z,t))∣∣∣z=zb(→X)=0. (23) (∂t¯η(x,t)+u(→X,z,t)∂x¯η(x,t)−w(→X,z,t))∣∣∣z=¯η(x,t)=0. (24)

Integrating (21) vertically (over ), and laterally (over ), applying the Leibnitz rule and using the kinematic boundary conditions, (23), (24), we have the following mass equation (25) below see . Repeating the same process for (22), gives the discharge equation (26) below:

 ∂tA(x,t)+∂xQ(x,t) =ΦA(x,t), (25) ∂tQ(x,t)+∂x(Q2(x,t)A(x,t)) =−gA(x,t)∂x¯η(x,t)+ΦQ(x,t), (26)

see [17, 12] for details; and are the coupling terms defined as

 ΦA(x,t)=∂xyr(x,¯η(x,t))[∫¯η(x,t)zb(→X)u(→X,z,t)dz]y=yr(x,¯η(x,t))−[∫¯η(x,t)zb(→X)v(→X,z,t)dz]y=yr(x,¯η(x,t))−∂xyl(x,¯η(x,t))[∫¯η(x,t)zb(→X)u(→X,z,t)dz]y=yl(x,¯η(x,t))+[∫¯η(x,t)zb(→X)v(→X,z,t)dz]y=yl(x,¯η(x,t)). (27) ΦQ(x,t)=∂xyr(x,¯η(x,t))[∫¯η(x,t)zb(→X)u2(→X,z,t)dz]y=yr(x,¯η(x,t))−[∫¯η(x,t)zb(→X)u(→X,z,t)v(→X,z,t)dz]y=yr(x,¯η(x,t))−∂xyl(x,¯η(x,t))[∫¯η(x,t)zb(→X)u2(→X,z,t)dz]y=yl(x,¯η(x,t))+[∫¯η(x,t)zb(→X)u(→X,z,t)v(→X,z,t)dz]y=yl(x,¯η(x,t)). (28)

Note that if the channel is not full, then we have and because non-full channel means by (6). Hence, both limits in all the integrals in (27)-(28), so the coupling terms vanish. In this case, the model, (25)-(26) reduces to the standard 1D Saint-Venant Models.

It is straight forward to show that

 ∂xyl,r(x,¯η(x,t))[∫¯η(x,t)zb(→X)θ(→X,z,t)dz]y=yl,r(x,¯η(x,t))=∂xywl,r(x)[∫¯η(x,t)zb(→X)θ(→X,z,t)dz]y=ywl,r(x)

where is any of the integrands appearing in (27)-(28). Therefore, we can conveniently write the coupling terms as follows:

 (29) ΦQ(x,t)=∂xywr(x)[∫¯η(x,t)zb(→X)u2(→X,z,t)dz]y=ywr(x)−[∫¯η(x,t)zb(→X)u(→X,z,t)v(→X,z,t)dz]y=ywr(x)−∂xywl(x)[∫¯η(x,t)zb(→X)u2(→X,z,t)dz]y=ywl(x)+[∫¯η(x,t)zb(→X)u(→X,z,t)v(→X,z,t)dz]y=ywl(x). (30)

To proceed, let us define the following quantities:

 qx(→X,t)=∫η(→X,t)zb(→X)u(→X,z,t)dz,qy(→X,t)=∫η(→X,t)zb(→X)v(→X,z,t)dz.So that q2xH≈∫η(→X,t)zb(→X)u2(→X,z,t)dz,qxqyH≈∫η(→X,t)zb(→X)u(→X,z,t)v(→X,z,t)dz(see \@@cite[cite]{[\@@bibref{Number}{chineduthesis, % toroshock}{}{}]}). (31)

Using equations (31) the coupling terms become

 ΦA(x,t)=qx|y=ywr(x)∂xywr(x)−qy|y=ywr(x)−qx|y=ywl(x)∂xywl(x)+qy|y=ywl(x). (32) ΦQ(x,t)=q2xH∣∣∣y=ywr(x)∂xywr(x)−qxqyH∣∣∣y=ywr(x)−q2xH∣∣∣y=ywl(x)∂xywl(x)+qxqyH∣∣∣y=ywl(x). (33)

#### Notational Simplification

We now express the coupling terms as functions of the fluxes at the channel lateral boundaries which are easier to compute.

Let and be the outward unit normal vectors to the lateral boundaries at and respectively (see figure 3(b)). Since these normal vectors are perpendicular to the tangent lines to their respective lateral boundaries, we have

 nylnxl∂xywl(x)=−1,nyrnxr∂xywr(x)=−1⟹∂xywl(x)=−nxlnyl,∂xywr(x)=−nxrnyr,nyl,nyr≠0.

Define the vector , then we write the coupling terms as

 ΦA(x,t)=ΦAL(x,t)+ΦAR(x,t) and ΦQ(x,t)=ΦQL(x,t)+ΦQR(x,t). (34)

where

 ΦAL(x,t)=1nyl(→q.→nl)∣∣∣y=ywl(x),ΦAR(x,t)=−1nyr(→q.→nr)∣∣∣y=ywr(x), (35) ΦQL(x,t)=1nyl(nxl[q2x(→X,t)H(→X,t)+g2H2(→X,t)]+nylqy(→X,t)qy(→X,t)H(→X,t))∣∣∣y=ywl(x)−nxlnylg2H2(→X,t)∣∣∣y=ywl(x), (36) ΦQR(x,t)=−1nyr(nxr[q2x(→X,t)H(→X,t)+g2H2(→X,t)]+nyrqy(→X,t)qy(→X,t)H(→X,t))∣∣∣y=ywr(x)+nxrnyrg2H2(→X,t)∣∣∣y=ywr(x). (37)

Let and denote the first and second components, respectively, of the outgoing flux in the direction of at , and and be those in the direction of , then the coupling terms (35)-(37) can be written in the flowing forms:

 ΦAL(x,t)=1nylf1L(x,t).ΦAR(x,t)=−1nyrf1R(x,t). (38) ΦQL(x,t)=1nylf2L(x,t)−nxlnylg2H2(→X,t)∣∣∣y=ywl(x), ΦQR(x,t)=−1nyrf2R(x,t)+nxrnyrg2H2(→X,t)∣∣∣y=ywr(x). (39)

Hence, the 1D channel models with coupling term, in the presence of friction is

 ∂tA(x,t)+∂xQ(x,t)=1nylf1L(x,t)−1nyrf1R(x,t), (40) ∂tQ(x,t)+∂x(Q2(x,t)A(x,t))=−gA(x,t)∂x¯η(x,t)+1nylf2L(x,t)−nxlnylg2H2(→X,t)∣∣∣y=ywl(x)−1nyrf2R(x,t)+nxrnyrg2H2(→X,t)∣∣∣y=ywr(x)+gA(x,t)Sf. (41)

where is the channel friction slope, is the conveyance, is the wetted perimeter of channel cross-section, and is the Manning coefficient, see [5, 11].

###### Remark 2.1

We obtained the above coupling terms without using or imposing any restriction on the channel width variation as done in . And our coupling term clearly differs from theirs.

#### 2.2.1 Channel Flow Lateral Discharge Model

To compute the lateral discharges in the channel, we use the following -discharge equation in the 2D Shallow water equations :

 ∂tqy(→X,t)+∂xfx(Π)+∂yfy(Π)=−gH(→X,t)∂yzb(→X),Π=(H,qx,qy)T,fx(Π)=qxqyH,fy(Π)=q2yH+12gH2. (42)

### 2.3 Floodplain Flow Model

We describe the flow in the floodplains using the 2D Shallow water equations, namely

 ∂tΠ+∇⋅F(Π)=S(Π,zb)+Sb(Π), (43)

where

 (44)

where is the manning coefficient, is the friction term and is the source term due to bottom topography term.

## 3 Numerical Schemes

In this section, we detail the numerical schemes for the models presented in previous sections. To begin, we partition the channel into a 1D grid, made of cross sections and the floodplains, into a 2D grid , see figure 4(a). We first present the scheme for the 2D flood model, then the schemes for the channel flow model is presented. (a) Grid of the entire domain consisting of the 1D grid Ω1Dh at the middle and the 2D grids Ω2Dh for the floodplains.

### 3.1 Scheme For 2D Model

In this section, we present the scheme for the flood flow model, (43). Let be an element of the 2D mesh in figure 4(a), and be its neighbour cell, see figure 4(b). Let be the edge between and , while is a unit vector normal to edge and outward to . Furthermore, let and be the area of and length of respectively and let be the set of all edges of . And let be the approximate cell averages of the true solution in , namely

 Πnj=1|Tj|∫TjΠ(→X,tn)d→x. (45)

Similarly, let be the cell average vector , while are the cell averages in respectively.

Then, we consider the following 2D hydrostatic reconstruction finite volume scheme :

 Πn+1j=Πnj−Δt|Tj|∑ejk∈Ej|ejk|(T−1→njkϕ(˜T→njkΠnj,˜T→njkΠnk)+T−1→njkShrm(Hnj,~Hn