A two-layer approach for Coupling 1D/2D Shallow Water Flow Models
Andreas Dedner, Chinedu Nwaigwe
TL;DR
This paper introduces a vertical coupling method for 1D/2D shallow water models that improves simulation accuracy and flexibility by modeling the channel as two layers, enabling better flow structure recovery and adaptability.
Contribution
The paper proposes a novel vertical coupling approach that treats the channel as two layers, allowing dynamic switching and improved flow structure preservation over existing methods.
Findings
The VCM recovers 2D channel flow structure without efficiency loss.
The method is well-balanced and preserves conservation properties.
It adapts to different flow situations effectively.
Abstract
In this paper, we propose a novel approach for coupling 2D/1D shallow water flow models. Efficiently coupling these models is vital for simulating the flow and flooding of open channels. Currently, existing methods couple the models either at the channel lateral boundaries (lateral methods) or at the location, along the channel flow direction, where the two sub-domains intersect (frontal methods). We classify these methods as horizontal methods since the coupling points are on the horizontal plane. The limitations of these methods include their inability to recover the 2D channel flow structure during flooding without losing efficiency, and their inability to switch between lateral and frontal methods, based on the problem. Here, we propose a new paradigm which is to think of the channel flow as a two layer flow. This leads to the vertical coupling method (VCM) which is able to recover…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsMeteorological Phenomena and Simulations · Flood Risk Assessment and Management · Hydrology and Watershed Management Studies
A two-layer approach for Coupling 1D/2D Shallow Water Flow Models
Andreas Dedner [email protected] Centre for Scientific Computing and Warwick Mathematics Institute, University of Warwick, United Kingdom
Chinedu Nwaigwe corresponding author, [email protected] Centre for Scientific Computing and Warwick Mathematics Institute, University of Warwick, United Kingdom
Abstract
In this paper, we propose a novel approach for coupling 2D/1D shallow water flows based on a two layer model in the channel, a 1D lower layer model and a 2D upper laayer model. The upper layer is only used in regions where flooding occurs otherwise the model reduces to a standard 1D channel model. To switch between the one layer and the two layer models the user prescribes an elevation above which the channel is considered to be full, i.e., a flooding event may be taking place. In the case of flooding the 2D upper layer model make it strightforward to couple the channel flow solver to a 2D shallow water solver used in the floodplain. We show that the resulting method (i) is well-balanced (ii) preserves a no-numerical flooding property (iii) preserves conservation properties of the underlying 1D and 2D finite volume schemes used for the flow in the channel and the floodplain. Numerical tests show that the method performs well compared to two horizontal coupling methods found in the literature. The results show that the method recovers the 2D flow structure in the channel in flooding regions, retains 1D flow structure in non-flooding regions while maintains good efficiency.
1 Introduction
One dimensional Saint Venant models are often used to simulate open channel flows but they become inadequate once the channel overflows. However, multi-dimensional (even 2D shallow water) simulations are computationally expensive. This has led to the development of methods to couple 1D channel simulations with 2D floodplain simulations. River/floodplain coupling simulations started as quasi 2D models in which floodplains are represented by storage cells, then coupled with an existing 1D river model [9, 4]. These approaches would not allow to simulate the fluid dynamics in the floodplains [10].
In [5], the coupling of a 1D channel model with a full 2D floodplain model was achieved by including the 2D numerical fluxes into the finite volume scheme for the 1D model, while [8] utilized the theory of characteristics to couple 1D/2D models through suitable matching conditions defined at the 2D/1D interface. In [20], the 1D river model and a 2D non-inertia model were also coupled where the water level differences between the flows in the two domains are used to calculate the interacting discharges in the sub-domains.
Methods based on post-processing the separately computed solutions are presented in [15], see also [14]. At 2D/1D interface, each model computes its own solution from which the total water volume in a 1D cell and all its adjacent 2D cells is computed. Then, the water height for 2D cells and the wetted Area for the associated 1D cell are found. In [17], the methods have been applied to Tiber River, Rome. In [13], the 1D model including the coupling terms were classically derived from the full 3D inviscid Euler’s equations and an optimal control process applied to couple the models. These models have been numerically treated with the finite volume method [10] where the discrete exchange term, which lead to globally well-balanced scheme, is proposed. This approach superposes a 2D grid over the 1D channel grid and convergence is achieved using a Schwartz-like iterative algorithm.
A major difficulty in coupling 2D/1D shallow water models is the computation of the channel flow lateral discharge. In [11], this channel lateral discharge was set to zero, while [12] adopted an iterative technique that uses the solution of successive Riemann problems to estimate the transverse velocity. This difficulty in computing the channel lateral discharge remains challenging to compute accurate fluxed between the channel and the floodplain.
A further issue is the assumption that the flow will remain one dimensional in the channel even in the region, where a flooding event is occurring. In [19] a method to compute different lateral discharges at each channel boundary was purposed. However, the free-surface and -velocity component were still assumed to be laterally constant across the channel and although the approach does improve the accuracy of the method, it still does not recover the complete 2D flow structure during flooding.
Efficient methods that recover the 2D flow structure during flooding but revert back to 1D simulation if no flooding is occurring could solve both the issues mentioned above. Frontal coupling methods in which the floodplain extends into the channel in parts of the domain, recover 2D flow structure but they loose efficiency because they compute the 2D solutions at all times. Moreover, most of the existing methods need to know the location of possible flooding a-priori and can not take into account that flooding locations may vary with time. So methods that can adapt to flooding regions a-posteriori are desirable.
The goal of this paper is to propose a method, the vertical coupling method (VCM), that (i) recovers the 2D flow structure during flooding while reverting back to 1D simulation in the channel regions where no flooding occurs. (ii) automatically detects flooding regions. (iii) can be easily added to well established 1D channel and 2D floodplain flow solvers.
The VCM is based on partitioning the flow in an overflowing channel into two vertical layers where the flows in the lower and upper layers are simulated using 1D and 2D models, respectively. We then derive a coupled models for the two layers which represents the channel flow. The upper layer model can be easily coupled to a 2D floodplain model. The method is parameterized by prescribing a height function along the channel which determines at which water height the channel is in danger of overflowing. Only when this water height is reached will the upper layer model be activated, below that layer the channel flow is simulated using a standard 1D model. Different choices for will lead to different method and special choices will lead to some standard lateral or frontal coupling methods found in the literature. In this sense VCM is a superset of some existing methods.
The rest of the paper is organised as follows: The two layer channel model is presented in section 2. We start by presenting the notation used to derive the VCM in section 2.1, derive the lower and upper layer flow models in sections 2.2 and 2.3 respectively, the complete two layer channel flow model is summarised in section 2.4. For the floodplain flow we use a standard shallow water model which is summarized in section 2.5. A numerical algorithm for the coupled channel flow models is presented in detail in section 3 and the properties of the method are considered in section 4 where we prove that the method is well-balanced, preserves no-numerical flooding and is mass conservative. Numerical experiments are presented in section 5 to evaluate the performance of the method compared to simpler approaches found in the literature. Finally, we give a summary in section 6.
2 Mathematical Models for the Fluid Dynamics
This section presents the model equations used for flow in the channel and fllodplain and derive, in detail, the two layer models for the channel flow.
2.1 Background
Let be a fixed 2D horizontal domain with fixed bottom elevation, , . Let denote the depth of water at point, at time, , so that
[TABLE]
is the water free-surface elevation at point, , at time . Then, at time, the flow occupies the 3D domain, defined by
[TABLE]
A cross section of the flow domain, at a fixed x is shown in figure 1(a).
2.1.1 Channel Geometry:
Figure 1(b) shows the flow cross section indicating the floodplain part and the channel cross section. To simplify the derivation of the models we assume its length lies along the -axis (frontal direction) and the width along the -axis (lateral direction). We define by
[TABLE]
The functions and are the -coordinates of the left and right lateral wall boundaries respectively at the elevation, . The function, gives the channel cross sectional lateral width, i.e.,
[TABLE]
and for convenience we set
[TABLE]
see figure 1(b).
2.1.2 The VCM Background
Given the channel geometry above, the formulation of the VCM starts by choosing an elevation, above which the channel is considered full, see figure 1(b). This elevation is completely decided by the user; the only constraint is for it not to be less than . So, becomes the channel top and is referred to as the maximum wall elevation defined below.
Definition 2.1** (Maximum channel wall elevation, )**
The maximum channel wall elevation at cross section x, is the minimum elevation of the channel banks above which flooding is said to have occurred [18].
Once has been chosen, other quantities are derived. The lateral width and the -coordinates of the lateral wall boundaries at the top () becomes and respectively. For the VCM we then partition the flow into the channel flow within the region , and the flow in the floodplain the remaining region, namely and , see figure 1(b). As noted above, the floodplain flows will be simulated using a standard 2D shallow water models (given in section 2.5), hence we concentrate on deriving the models for the flow in the channel region.
Remark 2.1
The idea of the VCM method is to use a 1D channel model for regions in the domain where the channel is not full while using a two layer model otherwise, where only the lower layer is evolved using a 1D model while the upper layer uses a 2D model. As said above the choice of when to switch from a one layer to a two layer model depends on the choice of and thus is up to the user.
The closer the chosen is to , the smaller the channel flow region , hence the larger the floodplain flow region (see figure 1(b)). In particular, choosing in some part of the channel, will result in an approach where the 2D floodplain model is used in parts of the channel (a frontal type coupling approach), while taking will result in a one layer 1D model being used in that region leading to a horizontal type coupling [18, 19]. Thus, different choices of lead to different coupling approaches.
Having identified the channel and floodplain flow regions, we now concentrate on the channel region. First, we note that the following condition holds for fixed :
[TABLE]
We also extend the definition of the width functions, to the region above the top (that is where ) as follows:
[TABLE]
see figure 1(b). This results in a channel with straight vertical walls above the height where the channel is assumed to be full.
For the scheme later on we need to assume that the channel geometry is known such that we can reconstruct the water height given the wetted cross section:
Definition 2.2** (Water height for given wetted cross section)**
Let be a given function describing the wetted area in the channel defined by the bottom topography . We define the height function with
[TABLE]
Remark 2.2
With the above definition is independent of and represent a 1D flow with the wetter area given by . Due to the way we extended for it is easy to see that the above definition provides a unique height function even for large value of , i.e., even in the case where the given wetted cross section leads to a overfull channel.
2.1.3 Important Quantities for the Channel Flow
We now proceed to define other important quantities for the channel flow.
Definition 2.3** (Channel Depth)**
The channel depth, is the laterally varying height between the channel bed and the chosen elevation, , that is
[TABLE]
see figure 1(b).
We also introduce the critical area defined as follows.
Definition 2.4** (Critical Area)**
The critical area, is the wetted cross sectional area of an exactly filled cross section. That is, the wetted cross sectional area when the water level is exactly at the chosen elevation, . It is defined by
[TABLE]
Remark 2.3
If the channel is exactly full (), then the water depth is exactly equal to the channel depth, , and consequently
[TABLE]
Definition 2.5** (Channel Flow Lateral Boundaries)**
Let and denote the -coordinates of the left and right channel flow lateral boundaries, i.e.,
[TABLE]
Note that if there is water everywhere in the channel, then .
Remark 2.4** (Channel Assumption)**
We assume that the channel never goes dry anywhere between the lateral wall boundaries, i.e., there are no islands in the channel:
[TABLE]
Definition 2.6** (Average Free Surface)**
With the flow lateral boundaries known, we can define the laterally averaged free surface elevation, namely
[TABLE]
Next, we define precisely what we mean by a full channel.
Definition 2.7** (Full Channel)**
We say a channel cross section at is full if
[TABLE]
(recall (1) for the definition of the water height and definition 2.3 for the channel depth ).
This means that if the channel is full the following statements hold:
- •
given the total wetted area in the cross section and the critical area as given in definition the following inequality holds (recall definition 2.4):
[TABLE]
- •
also since by definition
[TABLE]
2.1.4 Vertically Partitioned Channel Flow
The main idea of the VCM is to partition the channel flow into two layers, based on the chosen channel top elevation, . Since we want to approximate the flow with a 1D model if the channel is not full, we make the following assumption to allow consistency with a 1D formulations.
Definition 2.8** (1D Consistency Assumption)**
If the channel is not full, then the lateral variation in free surface elevation, is negligible and the free surface can be taken to be its lateral average, .
Remark 2.5
This simply means that we assume the free surface to be laterally flat. This is exactly the 1D assumption and allows to apply 1D modelling whenever the channel is not full (a similar assumption is made for example for the 1D models in [15, 5, 9]).
Next, we define the following elevation:
Definition 2.9** (Time Dependent Interface)**
The time dependent interface, is the elevation defined as,
[TABLE]
Proposition 2.1
.
Proof 2.1**.**
Case 1 : If the channel is full (see figure 2(a)), then by definition 2.7 we have . Hence . Furthermore, by definition 2.5, the lateral walls are .
Case 2 : If the channel is not full (see figure 2(b)), then . And definition 2.8 requires free surface to be constant and equal to the average, , hence the free surface at the lateral walls is . So, the lateral boundaries are .
We can now identify the channel flow domain, at time, as:
[TABLE]
Definition 2.10** **(Flow Partitions)
Given the total water depth, , we partition the flow into the two layers
[TABLE]
see figure 1(b). The water depths in the lower and upper layers are given by
[TABLE]
Note that defines the top of the lower layer and so by our 1D consistency assumptions is independent of .
2.1.5 Fundamental Equations:
Under the assumption of hydrostatic pressure, the flow in the domain, is governed by the incompressible free-surface Euler equations:
[TABLE]
where is the fluid velocity vector at point at time . Furthermore, the following boundary conditions hold
[TABLE]
[TABLE]
2.2 The Lower Layer Flow Model
The lower layer flow is assumed to be always one dimensional and thus the lateral variations in the free-surface elevation do not have any impact on the lower layer flow. Therefore, using the averaged free-surface, , instead of . We thus replace the horizontal moment equation (22) and the kinematic boundary condition (25) by
[TABLE]
We define the area of wetted cross section, , the section averaged volumetric discharge, and section averaged velocity, , for the lower layer, as follows:
[TABLE]
Before we proceed, let us state the following important relations which are easily proven (see [18] for details):
Lemma 2.1
Let and be defined in (9) and (14) respectively, and
[TABLE]
then
[TABLE]
Theorem 2.1
Let , , and be defined in (9), (14), (29) and (31) respectively, then
[TABLE]
2.2.1 Mass Conservation Equation for Lower Layer:
First we derive the model equations for the lower layer. Integrating the equation (21) over the lower layer cross section, and applying the Leibniz rule using the kinematic boundary condition (24), we obtain
[TABLE]
where
[TABLE]
Remark 2.6
Note that is the averaged mass exchange term between the two layers. If the channel is not full, i.e., then using the boundary condition (27)
2.2.2 Momentum Equation for Lower Layer:
To derive the momentum equation for the lower layer, we integrate equation (26) over the lower layer cross section and using standard arguments arrive at
[TABLE]
where the last integral on the right is the momentum exchange term between the two layers.
2.3 Upper Layer Flow Model
The upper layer flow is allowed to remain fully two-dimensional, so the Free-Surface Euler Equations, (21) - (23) and (25) are applicable. However, the kinematic boundary condition on the bottom does not apply here because the bottom of the upper layer is which is not a physical boundary that fluid particles cannot cross. Define the following quantities:
[TABLE]
So that the velocities are
[TABLE]
where is the upper layer discharge vector and is the velocity vector in the upper layer.
In the following, we derive the equations for the quantities. Integrating equation (21) vertically over the upper layer, we have
[TABLE]
Also, integrating equation (22) vertically over , applying the kinematic boundary condition (equation (25)) and simplifying, we have
[TABLE]
Similarly, integrating equation (23) vertically over , applying the kinematic boundary condition (equation (25)) and simplifying, we have
[TABLE]
The mass exchange between the layers is defined in equation (35).
2.4 Summary of Coupled Channel Flow Models
The two layer channel model derived so far consisting of the 1D lower layer model
[TABLE]
and the 2D upper layer model,
[TABLE]
where the fluxes are given by
, , \vec{u}_{\eta_{1}}=(u(\vec{\boldsymbol{X}},t),v(\vec{\boldsymbol{X}},t))^{T}\bigg{|}_{z=\eta_{1}}, , and .
Note that the models (42) and (43) are not closed since the exchange term and interface velocity are not known. Following a similar idea as used in [3] we will solve the system using a two step approach where in the first step we solve the equations without the exchange term
[TABLE]
and
[TABLE]
and then use this intermediate solution to approximate the mass and momentum exchange between the layers. But even then there is a difficulty in solving the 1D lower layer model since the free-surface elevation term, , appearing in (44) is different from the actual top level, of the lower layer flow, which is represented by . This makes directly deriving a well-balance scheme for the lower layer model challenging and more importantly we cannot directly reuse an existing 1D channel solver for the solution. Since this is one of our major goals, we base our numerical scheme on a reformulation of the above coupled channel model:
Integrating (43) over the cross section we arrive at
[TABLE]
where
[TABLE]
assuming that the average horizontal velocities in both layers are very similar, i.e., we can add the equations for and together arriving at a model for the 1D full channel
[TABLE]
So with this further approximation we end up with the following systems to model the flow within the channel:
[TABLE]
[TABLE]
where and .
Remark 2.7
The system (49) is a standard 1D approximation to the channel flow with providing the exchange terms between the channel and the floodplain along the horizontal boundary. Thus a standard model can be used to approximate (49). The advantage of our approach is that the second set of equations (50) provides additional 2D information in the channel that can be used to accurately compute the exchange terms .
2.5 Floodplain Flow Model
We describe the flow in the floodplains using the 2D Shallow water equations, namely
[TABLE]
where
[TABLE]
3 Numerical Algorithm for the Coupled Channel Flow Models
In this section, we describe in detail the numerical schemes and algorithms to solve the two-layer models, (49) and (50) for the channel flow. The 1D model (49) is solved on a 1D (quasi 2D) channel mesh (see figure 3(a)), while the 2D upper layer model (50) is solved on a 2D channel mesh (see figure 3(b)). The results are effectively combined to derive the update values for the channel flow.
3.1 Background
Let be a 1D (quasi 2D) channel mesh with cells , see figure 3(a). To consider the upper layer channel model, (43), we further discretize each 1D channel cell into 2D cells, , see figure 3(b). This forms the 2D channel mesh , so that
[TABLE]
For the cell, , we define the following discrete channel quantities:
[TABLE]
as the channel wall elevation, channel top width, critical area, 1D bed elevation, and channel depth, respectively, evaluated at the center of each 1D cell, (see figure 4(a)). We also define the discrete quantities
[TABLE]
denoting approximations to the average total wetted cross section and section-averaged discharge in .
For the associated cells, in , we define two similar sets of discrete quantities
[TABLE]
being the bed elevation and channel depth at the cell center, of , see figure 4(b). Note the identity, . Furthermore, let
[TABLE]
where and is the size(area) of . Finally we also define full 2D cell averages given by
[TABLE]
which are cell averages for the full 2D data (sum of lower and upper layer) in 2D cells, .
Given a lateral distribution of the free-surface elevation in the 2D cells for we assume that we can compute the wetted cross sectional area, in by
[TABLE]
This requires information about the channel geometry, for instance, if the channel has a rectangular cross-section, with a laterally constant bottom elevation, for all , then is given by
[TABLE]
Next, we formulate a discrete 1D consistency assumption similar to definition 2.8:
Definition 3.1** **(Discrete Consistency requirement)
The solution, at satisfies a discrete consistency requirement, if the following condition holds. If such that , then
[TABLE]
We will later show that if the initial data satisfies this discrete consistency requirement, then the requirement is satisfied for all time by our scheme.
3.2 Solution of the two-layer Channel Models
We now solve the two layer channel models, (49), (50) to evolve the full 2D data, from the initial to the final time. To this end we first initialize the value for the full 1D model and the values for the 2D upper layer model (in the region where the channel is full). To evolve these values to the next time level we first update the 2D upper layer values ignoring the exchange term between the two layers but taking the horizontal fluxes with the floodplain into account. Next we update the values of the full 1D model, using the channel-floodplain lateral fluxes computed by the 2D upper layer solver to completely solve (49). Finally, we correct the solution in the upper layer by estimating the exchange term between the two layers. Note that the first step only leads to computational cost in regions where the channel is full (i.e. ). Otherwise this step can be skipped and only the 1D step has to be carried out. Exchange terms with the floodplain are also correctly taken into account in the first step.
3.2.1 Initializing the two-layer model from given full 2D data
Given full 2D data, at the initial time and in each 2D channel cell , we initialize the 1D model, (49) and the 2D upper layer model, (50) as follows:
[TABLE]
where and is the average lateral width of .
3.3 Step 1: Solution of Upper Layer Model without Exchange Terms
We now propose a method to evolve the solution for the model, (45). Let us emphasis that any appropriate method for a 2D shallow water model can be used here since (45) can be interpreted as a standard 2D shallow water equation with topography, referred to as the apparent topography in [6, 7]. In our simulation we use an apparent topography hydrostatic reconstruction scheme [1, 6, 7]:
Let be a 2D cell in a 2D channel mesh and be its neighbour. Let be the edge between and , and be the unit vector normal outwards to . Furthermore, let and be the area of and length of respectively and let be the set of all edges of . Let and be the lower layer water depth, channel depth, bottom elevation and upper layer cell average vector in 2D channel cell, ; while and are those in 2D cell, . Note that the neighbour cell, could actually be in the floodplain; the only requirement is that the current cell, in which we compute the upper layer solution, must be in the channel.
Define the apparent topography:
[TABLE]
Then reconstruct the apparent topography in the hydrostatic fashion:
[TABLE]
Next, define
[TABLE]
Then, the apparent topography hydrostatic reconstruction scheme [6, 7] for the model reads:
[TABLE]
where is the friction term defined in (52), and
[TABLE]
3.4 Step 2: Complete Solution of the Full 1D Model
We now solve the full 1D model, (49). This involves two sub-steps, namely
3.4.1 Step 2.1: Black-Box Stage
Here, we use any available 1D solver to solve the full 1D model without the lateral flux term, , namely
[TABLE]
with the data, . In this paper, we adopt the scheme from [15, 16], see also [18, 19]. This gives the approximation, .
3.4.2 Step 2.1: Add the floodplain coupling term
The discrete coupling term, is the vector consisting of the first two components of the 2D lateral fluxes already computed by the upper layer solver in section 3.3. Therefore, the final update value for the full 1D model, (49) is then given by
[TABLE]
where is the time step.
3.5 Step 3: Final Update of the Upper Layer Model
With the intermediate solution, known, we now completely solve the upper layer model including the exchange term, (50). The approximate solution of (50) is the approximate solution of the system
[TABLE]
with the initial data, , where
[TABLE]
Then using forward Euler time discretization, the approximate solution of (71) is
[TABLE]
To complete the description of the scheme, we need to define the exchange terms and interface velocities and . But there are no equations for these terms. However they can be determined by requiring that the following conditions be satisfied: (i) The operation from the intermediate solutions to the new update solutions must locally (and globally) conserve both mass and momentum. (ii) The final update values must satisfy the discrete consistency requirement (definition 3.1). (iii) the new update values must satisfy the non-negativity of water heights. These conditions allow us to first obtain the heights , then the exchange terms, , and finally the interface velocities are calculated so that we can compute .
3.5.1 Step 3.1: Approximating the Lower Layer Flow
By (33) we have
[TABLE]
given computed above. Note that we can not directly apply since we do not yet know . However, by mass/momentum conservation, the following must hold
[TABLE]
Hence the intermediate lower layer wetted area is given by
[TABLE]
Using this lower layer update, and the intermediate solutions, and we can now compute the upper layer heights in the next step.
3.5.2 Step 3.2: Upper Layer Heights
Using the definition of from (73) we have to distinguish two cases and :
Case 1: If
(Lower Layer not full at )
Recall, , so we have . Since , we must have
[TABLE]
Case 2: If (Lower Layer full at )
Denoting , and replacing by , then we can write the first equation in (74) in the following form:
[TABLE]
and consider two further cases.
Case 2a: If (Lower Layer full at intermediate state)
Then by an amount , so we add the constant excess height, to the intermediate solution uniformly over all 2D cells, namely
[TABLE]
Case 2b: If (Lower Layer not full at intermediate state)
Then by the amount, . Hence we remove from using the following algorithm.
- •
initialize for all .
- •
.
- •
while( )
- i
for all
- (a)
, 2. (b)
Reduce upper layer height by the gap height :
[TABLE] 3. (c)
Remove area of reduced height from total gap area:
, 2. ii
.
Equations (76), (78) and (79) compute in for all cells.
3.5.3 Step 3.3: Upper Layer Discharge and Lower Layer Discharge
Having computed , we compute the exchange terms, using the first equation in (72), hence
[TABLE]
and compute the interface velocity, following [3] namely
[TABLE]
Thus the interface velocity is equal to the upper layer velocity if mass is flowing from the upper to the lower layer and otherwise equal to the lower layer velocity.
Using the above definitions of the exchange terms and interface velocity, we define the upper layer discharge by using the second equation in (72).
This completes the upper layer solution .
Implicitly we can also compute the lower layer solution by using (73) and local momentum conservation
[TABLE]
which leads to the lower layer discharge given by
[TABLE]
3.6 Obtaining the full data ,
and for
This final step is only required for postprocessing and for studying the properties of the scheme in the following section.
The full water depth, is given by
[TABLE]
where
[TABLE]
The full 2D -discharge is computed as
[TABLE]
Finally, the -discharge is computed as follows :
[TABLE]
where
[TABLE]
This completes the numerical algorithm for the channel flow, for the floodplain flow we use a hydrostatic reconstruction scheme as described for the upper layer model, see [18] for details. At the channel/floodplain interface, the lateral fluxes are effortlessly computed since in the 2D channel cell we have complete 2D value at least if the channel is not full. In this important case, we have an approximation of the vertical discharge which is not the case for most existing methods where only a 1D model is used in the channel which does not provide values for the lateral discharge. If the channel is not full but the surrounding floodplain is not dry, the channel 1D flow data is used to compute the coupling fluxes.
4 Properties of the Scheme
We now consider some important properties of the scheme proposed above.
Definition 4.1** **(Consistency of Distribution Operation)
We require a distribution operation such as the one in section 3.2.1 to be consistent in the sense that whenever the channel is not full ( for all ), then the following conditions must hold: and each equal zero .
This ensure that we do not solve a two-layer problem when the channel is not full.
Definition 4.2** **(No-Numerical Flooding Properties)
In order to satisfy the no-numerical flooding property, the full 2D data must satisfy the following condition for all :
- i
Either the lower layer is full or the upper layer is empty, i.e.,
[TABLE] 2. ii
If , then (a) (b) (which is constant laterally) and (c) ; 3. iii
*If such that , then . *
Condition means that the lower layer is either full or upper layer is dry . In other words, there should not be gap between the two layers. A consequence of this property is also that no overflowing of the channel can occure unless it is full. Condition () states that if the upper layer is dry, the upper layer velocities and discharges must vanish, and the full layer flow velocity, must be laterally uniform. Condition means that if the channel is not full, then the free-surface must be flat, i.e., not vary in the lateral direction. Thus the discrete consistency requirement (definition 3.1) is satisfied at all time steps.
In the following we will show that if the solution at time satisfies the above no numerical flooding properties then this is also true at time . In addition, we also prove that the scheme is well-balanced and conserves mass under suitable conditions on the intermediate solutions.
Theorem 4.1** **(Consistency of Distribution Operation)
The distribution operation proposed in section 3.2.1 is consistent with the problem in the sense of definition 4.1.
Proof 4.1**.**
We need to prove that the ascertion in definition 4.1 is true. Let the channel not be full, that is, . Then, we have and (see (63)).
Theorem 4.2** **(No-Numerical Flooding Property)
The vertical coupling scheme as derived in (84), (3.6) and (87), preserves the no-numerical flooding property in the sense of definition 4.2.
Proof 4.2**.**
- i
We prove (89) on case-by-case bases (see section 3.5.2).
Case 1: , then
[TABLE]
Cases 2a and 2b : , then
[TABLE]
Therefore, we have or in either case. Hence,
[TABLE]
as claimed. 2. ii
*Let , then *
(a)
[TABLE]
(b) Since , then by (3.6), we have
[TABLE]
Hence,
[TABLE]
(c)
[TABLE] 3. iii
Assume that for we have , then this corresponds to case 1 in section 3.5.2 because cases 2a and 2b satisfy
[TABLE]
Since, corresponds to case 1, then it satisfies
[TABLE]
So,
[TABLE]
That is is independent of , hence also equal to . Therefore,
[TABLE]
Hence the free surface is laterally flat. Therefore, the discrete consistency requirement (definition 3.1) holds at , and the scheme satisfies the no-numerical flooding property.
Theorem 4.3** **(Well-balanced property)
If the numerical schemes used to compute the intermediate solutions, and are well-balanced, then the vertical coupling method, (84), (3.6) and (87), is well-balanced. This is true for any kind of well-balance, not only for lake at rest.
Proof 4.3**.**
Let the intermediate solutions be well-balanced, i.e.,
[TABLE]
*We need to show that .
First, we show that . Recall that by (74)*
[TABLE]
*Hence,
Secondly, we show that using (74) and that :
[TABLE]
Hence,
**Thirdly, we use the above results to show that .
Case 1: . We have already shown that so and . Hence, in case 1.
Cases 2a and 2b: , that is using . But**
[TABLE]
Hence by (78) and (79), we have
[TABLE]
**Fourtly, we show that and .
Since , then we have using (80), (72) that and .**
Also,
[TABLE]
Finally, we use the above results to show that . By (84)-(87), we have
[TABLE]
Hence, we have shown that , so the method is well-balanced.
Next, we prove that the proposed schemes for the intermediate solutions are mass conservative and also well-balanced for lake at rest. Since the hydrostatic reconstruction method is mass conservative, [2, 1] and the model, (45) does not introduce any source term to the height equation, then the scheme (67) is mass conservative. We state and prove a theorem below to show that this scheme, like the standard hydrostatic reconstruction [1] method, preserves well-balance of lake at rest.
Theorem 4.4
The upper layer scheme, (67) is well-balanced with respect to lake at rest.
Proof 4.4**.**
*The proof follows the same lines given in [2, 1], for details see [18]. *
Theorem 4.5
If the underlying 1D solver for (48) is well-balanced then the intermediate lower layer scheme, (75) is also well-balanced.
Proof 4.5**.**
Let scheme for (48) be well-balanced, then . Since the solver, (67), is also well-balanced, so . Therefore, the lower layer scheme, (75) becomes
[TABLE]
So, the scheme is well-balanced.
Theorem 4.6** **(Conservation)
If the intermediate solutions are mass conservative, then the vertical coupling solution is also mass conservative.
Proof 4.6**.**
Let the intermediate solutions be mass conservative, then
[TABLE]
By (74), we have
[TABLE]
Hence,
[TABLE]
Theorem 4.7
If the underlying solver for (48) is mass conservative, then the intermediate lower layer scheme, (75) is also mass-conservative.
Proof 4.7**.**
Let scheme for (48) be mass conservative, then . Since the solver, (67), is also mass conservative, so . Therefore, (75) gives
[TABLE]
So, the scheme is mass conservative.
Theorem 4.8
The vertical coupling method (VCM), described in (84)- (88) is well-balanced with respect to lake at rest and is mass conservative.
Proof 4.8**.**
Since the solver [15] is well-balanced for lake at rest and mass conservative, then by theorems 4.5 and 4.7 the same holds for the intermediate lower layer scheme. Since the intermediate upper layer scheme is also well-balanced and conservative the results follows using theorems 4.3, 4.6.
This completes the theoretical aspect of the VCM. In the next section, we present some numerical experiments to evaluate its performance compared to other coupling methods.
5 Numerical Results
We consider two simple test cases to verify the performance of the proposed method. We use the full 2D simulation results as the reference solution, and compare the results of the VCM with those of the Horizontal Coupling Method (HCM) [19] and the Flux-Based Method (FBM) [5]. Recall that in addition to the usual data required for the simulation, e..g, domain size and intitial conditions, the VCM method also requires a choice of the function used to determin when a channel is considered full and flodding might occure.
5.1 Test Case 1 : Dam-Break Flow into a Flat Floodplain
We consider a dam break flow in a 19.3 meter long, 0.5 meter constant width flat channel with adjacent flat floodplain [21, 15, 18], see figure 5(a). The labels are chosen probe points in the flow domain. To compare our simulations with the literature we add bottom friction terms to the channel and the floodplain models using a manning coefficient of s/m1/3. The boundaries are all closed walls except the right side as indicated in the figure. The wall elevation within the channel is meters. For there is no channel wall. The initial flow condition is given by
[TABLE]
To apply the VCM to this problem, we consider the following :
[TABLE]
5.1.1 Result of test 1:
The four simulation methods were all run with a grid of cells in the floodplain. For the channel region, the full 2D simulation used 2D cells, the VCM used 193 1D cells and an upper layer 2D grid of cells (that is, 8 2D upper layer cells per one 1D cell), while the HCM and FBM used 193 1D cells each. Each simulation was run for ten seconds with a CFL number of 0.95.
The free surface elevation at the last time step is shown in figure 5 for the full 2D, VCM and HCM. The figure shows that the coupling methods, VCM and HCM, approximates the full flow field with good accuracy, however, one can also see that the VCM computes better approximation than the HCM. Especially the 2D flow structure in the right part of the channel where the upper layer is active, is correctly captured by the VCM, unlike the HCM. The accuracy of the vertical coupling method is further illustrated in figure 6 which displays the time evolution of the free surface elevation, at the selected probe points. It can be seen that the VCM captures the full 2D results better than both the HCM and FBM, at all the probe points. Again, from figure 5, one can see that the VCM recovers 2D flow structure within the channel unlike the HCM and FBM. For the this test case, the vertical coupling method results in about 48% gain in computational time over the full 2D simulation. The fastest method (FMB) leads to a gain of about 55%. Hence in this test case, the VCM truly has shown good accuracy improving on the simple flux based coupling method while retaining most of the gain in efficiency.
5.2 Test case 2 : Channel Flow into Elevated 2D Floodplain
The second test case was suggested in [18], see also [19]. The set up consists of a channel of length 19.3 metre, width 0.5 metres, zero bottom topography and connected to a 0.5 metre high floodplain which is located in the region, , see figure 5(b). The manning coefficient for both channel and floodplain is taken as 0.009 the boundaries are only open at the sides indicated "exit" in figure 5(b), others are closed. Nine probe points are considered. From the bottom topographies for the channel and floodplain given above, the channel wall elevation is equal to meters within the channel except in a breach region for where the wall height drops to meters. The initial condition is the following.
[TABLE]
To apply the vertical coupling method to this problem, we consider the following smoother version of the wall elevation:
[TABLE]
5.2.1 Result of test 2:
For all four methods we simulated this problem with grid cells in the floodplain, while the grid resolutions in the channel are exactly the same as in the first test case. The simulation was run for ten seconds and CFL of .
The free surface elevation and velocity magnitude are shown in figures 7 and 8 respectively. We can see that that the VCM computes a better approximation of the full 2D solution than the HCM which in turn, is more accurate than the FBM. Moreover, the VCM reproduces a non laterally constant free surface elevation and velocity within the channel, which can not be achieved by either HCM or FBM. To further understand the results, the free surface elevation, the -component and -component velocity are plotted for selected probe points in figures 9. It can be seen that the vertical coupling method is more accurate than the other methods for all three flow quantities at all probe points and almost throughout the duration of simulation. Again, the vertical coupling method really captures the flow structure of the full 2D simulation.
Finally, figures 5 and 8 show that the VCM, unlike the other methods, recovers the 2D flow structure within the channel at the flooding regions. Again, the VCM continues to compute 1D solutions at non flooding regions; this demonstrates the self-adaptive nature of the method.
In this example VCM was about 75% slower then the computationally much simpler FBM method and about 35% more efficient then the full 2D simulation.
For this problem the VCM clearly improved the accuracy of the simulation both within the channel and in the floodplain compared to simpler coupling method but at a increased computational cost. But note that in this examples the percentage of the domain where a 1D assumption for the flow is valid is quite small so that it is not surprising that the gain in computational efficiency between full 2D and VCM is not so large. For a simulation of a very large network of rivers where the 2D region might be very small compared to the entire computational domain, the difference in efficiency of one coupling method over another will be far less significant, so that accuracy becomes the deciding factor. In this regard, the VCM is the best method of the three considered here.
6 Conclusion
In this paper we investigated a new method for coupling 2D and 1D shallow water models. We focused on the need to efficiently recover 2D channel flow structure during a flooding event to achieve high accuracy while maintaining the efficiency of the 1D channel model as much as possible. To this end we presented a vertical coupling approach which adds a second 2D layer inside the channel in regions where overflow can be expected to occur. This makes coupling the channel flow to the floodplain straightforward since 2D information is always available also within the channel when required. Any standard 1D channel flow and 2D floodplain solver can be used as building blocks for the new VCM method and only a slight modification of standard 2D solvers is required for the evolution of the second layer. We proved that the resulting method retains many properties of the 1D and 2D solvers used, e.g., mass conservation and well balancing. In addition we also studied a no-numerical flooding property. Our numerical results show that the VCM, in most cases, outperforms the other methods studied in this paper and does accurately recover 2D flow structures also within the channel when flooding occurs.
More detailed investigations with more complex channel geometries will be the focus of our further research. As mention in this paper, VCM actually consists of a very large family of methods. It depends on the choice of which is used to determine when a channel is considered full and the 2D layer model is used. For example, taking everywhere leads to a horizontal flux coupling method similar to the FBM method. On the other hand taking in regions in danger of flooding and very large away from these regions using a smooth transition, results in a frontal type coupling method where a standard 2D solver is used in the region where and on the boundary of this region the VCM will lead to a blending type approach between the 2D and the 1D regions. The width of the blending region will depend on how changes from [math] to . Further tests are required to understand the influence of different choices of and are the focus of ongoing work.
Acknowledgement
We are grateful to the Petroleum Technology Development Fund (PTDF), Nigeria for funding this study and to the Centre for Scientific Computing, University of Warwick for providing the computing resources.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] E. Audusse, F. Bouchut, M. Bristeau, R. Klein, and B. Perthame. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM Journal Scientific Computing , 25:2050–2065, 2004.
- 2[2] E. Audusse and M. Bristeau. A well-balanced positivity preserving second-order scheme for shallow water flows on unstructured meshes. Journal of Computational Physics , 206(1):311–333, 2005.
- 3[3] E. Audusse, M. Bristeau, B. Perthame, and J. Sainte-Marie. A multilayer saint-venant system with mass exchanges for shallow water flows. derivation and numerical validation. ESAIM: Mathematical Modelling and Numerical Analysis , 45(01):169–200, 2011.
- 4[4] E. Bladé, M. Gómez, and J. Dolz. Quasi-two dimensional modelling of flood routing in rivers and flood plains by means of storage cells. In Modelling of Flood Propagation Over Initially Dry Areas , pages 156–170. ASCE, 1994.
- 5[5] E. Bladé, M. Gómez-Valentín, J. Dolz, J. Aragón-Hernández, G. Corestein, and M. Sánchez-Juny. Integration of 1d and 2d finite volume schemes for computations of water flow in natural channels. Advances in Water Resources , 42:17–29, 2012.
- 6[6] F. Bouchut. Nonlinear stability of finite Volume Methods for hyperbolic conservation laws: And Well-Balanced schemes for sources . Springer Science & Business Media, 2004.
- 7[7] F. Bouchut. Efficient Numerical Finite Volume Schemes for Shallow Water Models . Elsevier, 2007.
- 8[8] Y. Chen, Z. Wang, Z. Liu, and D. Zhu. 1d-2d coupled numerical model for shallow-water flows. Journal of Hydraulic Engineering. , 138:122–132, 2012.