Formulation, Implementation and Validation of the Horizontal Coupling Method for 1D/2D Shallow Water Flow Models
Chinedu Nwaigwe, Andreas Dedner
TL;DR
This paper introduces the Horizontal Coupling Method (HCM), a novel efficient numerical scheme for coupling 1D and 2D shallow water flow models to improve accuracy in flood simulations without high computational costs.
Contribution
The paper develops a new coupling approach for 1D and 2D shallow water models that accurately captures lateral discharges and is computationally efficient.
Findings
HCM is well-balanced and satisfies no-numerical flooding property.
The method accurately approximates lateral discharges with minimal computational overhead.
Numerical tests confirm high accuracy and efficiency of the proposed coupling scheme.
Abstract
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. %leading to a new coupling term. Finite volume methods are formulated for both the 2D and 1D models including a discrete coupling term in closed form.…
| Channel Grid | Floodplain Grid | No. of time steps | Processor time (in seconds) | |
|---|---|---|---|---|
| Full 2D | 3,669 | 3,110.31 | ||
| HCM | 2,616 | 1,420.4 | ||
| FBM | 2,592 | 1,311.26 |
| Channel Grid | Floodplain Grid | No. of time steps | Processor time (in seconds) | |
|---|---|---|---|---|
| Full 2D | 4,963 | 4,100.47 | ||
| HCM | 3,235 | 1,710.36 | ||
| FBM | 3,178 | 1,555.08 |
| Channel Grid | Floodplain Grid | |
|---|---|---|
| Full 2D | ||
| HCM | ||
| FBM |
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
TopicsComputational Fluid Dynamics and Aerodynamics · Meteorological Phenomena and Simulations · Fluid Dynamics and Turbulent Flows
Formulation, Implementation and Validation of the Horizontal Coupling Method for 1D/2D Shallow Water Flow Models
C. Nwiagwe [email protected] Centre for Scientific Computing and Warwick Mathematics Institute, University of Warwick, United Kingdom
A. S. Dedner [email protected] Centre for Scientific Computing and Warwick Mathematics Institute, University of Warwick, United Kingdom
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 [17] 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 [2], 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 [4]; 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 [18] 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 [14], see also [13], 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 [16]. The superposition approach, proposed in [12], 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 [12] was extended to finite volume methods in [6], 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 [8].
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 [7]. In [8], the exchange terms derived in [12] 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 [12] 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 [12] (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
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
[TABLE]
bounded below by a fixed bottom, and above by the water free-surface elevation, given by
[TABLE]
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, )**
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
[TABLE]
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
[TABLE]
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
[TABLE]
such that
[TABLE]
and the bottom elevation satisfies
[TABLE]
see figure 2(b). Furthermore, we extend the definition of the width functions above the channel top () in the following:
[TABLE]
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
[TABLE]
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
[TABLE]
(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 [10]:
[TABLE]
[TABLE]
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:
[TABLE]
So that the section-averaged velocity, is given as
[TABLE]
First, we note that -independence of the free-surface, means that the sum,
[TABLE]
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:
[TABLE]
Therefore, the FSEE (10) -(16), reduce to the following system:
[TABLE]
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 [12]. Repeating the same process for (22), gives the discharge equation (26) below:
[TABLE]
see [17, 12] for details; and are the coupling terms defined as
[TABLE]
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
[TABLE]
where is any of the integrands appearing in (27)-(28). Therefore, we can conveniently write the coupling terms as follows:
[TABLE]
To proceed, let us define the following quantities:
[TABLE]
Using equations (31) the coupling terms become
[TABLE]
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
[TABLE]
Define the vector , then we write the coupling terms as
[TABLE]
where
[TABLE]
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:
[TABLE]
Hence, the 1D channel models with coupling term, in the presence of friction is
[TABLE]
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 [12]. 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 :
[TABLE]
2.3 Floodplain Flow Model
We describe the flow in the floodplains using the 2D Shallow water equations, namely
[TABLE]
where
[TABLE]
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.
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
[TABLE]
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 [1]:
[TABLE]
where
[TABLE]
The function, is any numerical flux function consistent with the 1D component, of the 2D flux, . Here, we consider the HLL scheme [9], namely
[TABLE]
where and are estimates of the smallest and largest wave speeds in the solution of the associated 1D Riemann problem [20]. There are several choices for [20, 19]. We use the ones given in [3] namely
[TABLE]
where , are the eigenvalues of the Jacobian matrix of the system.
And is a rotation matrix which depends on the normal vector, and is its inverse; they are given by
[TABLE]
see [20].
3.2 Schemes for the Channel Models
Here, we describe the finite volume method to discretize the channel models with coupling terms, equations (40), (41) and (42). To design a method which reuses existing 1D channel solvers, we discretize the purely 1D channel models, (40) and (41) separately from the lateral discharge model, (42). Let be points in the 1D grid, and be a cell centred at in . Where is the number of cells in the 1D grid. Let be a vector of conserved quantities at point, and time, , then the cell average vector, in cell is defined as
[TABLE]
where , and is the time step.
For each 1D cell, , the channel lateral boundaries, and are approximated with straight edges which we call South (S) and North (N) edges (or faces) respectively with unit normals and (see figure 5). This means that the channel normals, and are approximated with the edge normals and respectively, that is
[TABLE]
We start by presenting the scheme of the 1D model without coupling terms, next the scheme for channel lateral discharge model is presented and finally, the discrete coupling term is derived.
3.2.1 Scheme for 1D Model without Coupling Terms
We now focus on the 1D channel model (40)-(41) but without the coupling terms, namely,
[TABLE]
We consider the scheme of [15] as summarised in [14]. The scheme is based on the formulation of the St Venant model as presented in [5] and rewrites the model in the quasi-linear form.
[TABLE]
where , the Jacobian matrix, is given by
[TABLE]
is the top width at the free-surface, is the water depth from the 1D bottom elevation, to the flat free-surface, and is the negative of channel bed slope. Details about this formulation can be found in [14]. The eigenvalues and eigenvectors of are
[TABLE]
respectively.
Define the Roe averages:
[TABLE]
Define for any quantity, . Then, define
[TABLE]
The Roe averaged eigenvalues and eigenvectors are
[TABLE]
The artificial viscosity (entropy fix), to correct the entropy problem associated with the Roe method [14] is given by
[TABLE]
Hence, the numerical scheme of [14, 15] for the 1D channel model without coupling term, (54) is given by
[TABLE]
where
[TABLE]
3.2.2 Approximating Channel Lateral Discharge
Here, the goal to solve the lateral discharge model, (42) along the channel. Consider the 1D channel cell, , with cell average vector and denote by and , the -th 2D floodplain cells adjacent to on its North edge and South edge respectively, see figure 6(a). Let the cell averages in the adjacent 2D cells and be
[TABLE]
respectively. and are the number of the adjacent 2D cells on the North and South edges, respectively of . and are the average water depth, average discharge along -direction and average discharge along -direction respectively, in 2D cell while and are those of cell .
To discretize the lateral discharge model (42) in , we subdivide into two subcells, and and view them as 2D cells within the channel, see figure 6(b). Let and be the 2D cell average vectors in the subcells, and respectively. Then, we define them as
[TABLE]
where and are the 1D cell average water depth and section-averaged velocity in the channel, and are computed at every time step using the scheme presented below.
The motivation to compute is to apply a well-balanced scheme to the model,(42) in the subcells by taking the bottom to be flat across all the edges within the channel, (see figure 6(b)). To this end, we define the following:
[TABLE]
where is the discrete flat free surface elevation in 1D cell, and is the bed elevation of the adjacent 2D cell , see figure 6(b). Therefore, we propose the following hydrostatic reconstruction scheme [1] for the lateral discharge in subcell :
[TABLE]
denotes the 3rd component of numerical flux, . The quantities; and are the lengths of the corresponding edges of and is the outward unit normal of towards (see figure 6(b)).
Similarly, by defining
[TABLE]
we propose the following scheme for the lateral discharge in :
[TABLE]
where is length of edge, between and , and is the bed elevation of 2D cell, . The initial values of are obtained as explained in the following remark.
Remark 3.1
At initial time (), only the lateral discharge, for the full cell, is given. Then we initialize and to be equal to it, namely
[TABLE]
For other times, (), we compute using the scheme described above.
3.2.3 Discrete Coupling Terms
To discretize the coupling term, let us denote by the approximation of a function, at the edge, between 1D cell and 2D cell (see figure 6(b)). Then the approximation, of over the entire South edge, of , is given by averaging over all edges on the South edge, namely
[TABLE]
where is the length of edge, ; is sum of all edges of on the South side. Similarly (f|_{R})_{i}=\sum_{j=1}^{Nn}\bigg{(}f_{i,j}^{N}\frac{|e_{ij}^{N}|}{|e_{i}^{N}|}\bigg{)} for North edge, where is the length of edge between and , and is the sum of all edges on North side of . Hence we can approximate the coupling term as
[TABLE]
where
[TABLE]
are the discrete coupling terms at the edges and respectively (see figure 6(b)). and are respectively, the discrete water depth, first and second components of 2D numerical flux at edge . While and are respectively, the water depth, first and second components of 2D numerical flux at edge, . We now focus on how to compute them.
Given the 1D cell average, in from which we obtain the cell average in the subcell, (using equation (67)). Then we directly approximate and by computing the 2D numerical flux, , at edge, (see figure 6(b)), namely
[TABLE]
Similarly, by using and we approximate and using
[TABLE]
The hydrostatically reconstructed quantities, and are defined in equations (68) and (70) respectively.
To approximate and , we propose to adapt the hydrostatic reconstruction approach [1], namely
[TABLE]
We therefore summarise the discrete coupling term as
[TABLE]
3.3 Summary of the Channel Flow Solver
The complete scheme for the channel flow model with coupling term, (40)-(41) is
[TABLE]
where is the solution of the purely channel model without the coupling terms,(54) which is given in (64). is the discrete coupling term summarised in (79). The channel lateral discharges are computed using the schemes in (69) and (71).
3.4 Summary of the Horizontal Coupling Method
The channel flow model, (40), (41) and (42) is simulated as summarised in section 3.3 while the flood flow model, (43) is solved with the 2D solver, (46). At 2D/1D edge, the 2D numerical flux is computed by using the 2D cell, averages and the averages obtained from the adjacent channel subcell, ( or ) as described in (67). A flow chart for the implementation of the HCM is given in figure 7.
4 Properties of the HCM
We discuss a few properties of the method in this section.
Definition 4.1** (Well Balance of Lake at rest)**
Assuming that the existing numerical schemes for the uncoupled 1D and 2D models are well balanced with respect to lake at rest, then the coupled scheme is said to be well balanced with respect to lake at rest if the coupling term vanishes whenever the lake at rest condition holds.
Theorem 4.1
The coupling term derived in equation (79) leads to a fully well-balanced scheme with respect to lake at rest.
proof 4.1
Assuming that the condition of water at rest holds, then
[TABLE]
where and are the free surface elevation in the adjacent 2D cells, and respectively. Hence,
[TABLE]
Therefore,
[TABLE]
Since all velocities (and discharges) are zero, then
[TABLE]
Hence, by the consistency of the numerical flux, with the physical flux, , we have
[TABLE]
That is
[TABLE]
Similarly,
[TABLE]
So that
[TABLE]
Therefore, using equations (81)-(83), then the discrete coupling term in (79) becomes
[TABLE]
as claimed.
We now introduce the concept of "No-Numerical Flooding".
Definition 4.2** (No Numerical Flooding Property)**
We shall say that a 2D/1D coupling scheme preseves the No Numerical Flooding property, if all its coupling terms vanish whenever there is no flooding or draining.
Theorem 4.2
The scheme (80) preserves the no numerical flooding property.
proof 4.2
If no flooding, then
[TABLE]
Hence,
[TABLE]
Again, since floodplain is dry, we have
[TABLE]
These give the numerical fluxes:
[TABLE]
So all the flux terms in are zero. Finally,
[TABLE]
Therefore, . Which means that no water is gained from or lost to the floodplain as required.
5 Numerical Results
In this section, we present some numerical experiments to investigate the performance of the proposed method. We use full 2D simulation results as the reference solution and compare these results with those of the HCM and of the flux-based method (FBM) of [2]. All the algorithms are implemented in a C++ code and the experiments are run on the Cluster of Workstations (COW) of the Centre for Scientific Computing, University of Warwick, United Kingdom.
5.1 Test Case 1 : Dam-Break Flow into a Flat Floodplain
The first test case is suggested in [14]. The setup consists of a dam break flow in a 19.3 meter long, 0.5 meter constant width flat channel with adjacent flat floodplain, see figure 8. The National Laboratory of Civil Engineering in the IST in Portugal designed and measured this test case [21, 14]. A reservoir is located from the left end of the channel to 6.10 metres (position of dam in figure 8). The initial condition is
[TABLE]
The manning coefficient, (see section 2.2) for both channel and floodplain is s/m1/3 and the boundaries are all closed walls except the right side as indicated in figure 8. The labels are probe points in the flow domain. More about this test case can be found in [14, 21].
Here, a full 2D simulation was run with a grid of cells in the channel and cells in the floodplain, while a simulation with the proposed method was run with a grid of cells for the floodplain and cells in the channel, and the simulation using the FBM was also run with a grid of cells in floodplain and in the channel, see table 1.
Figure 9 displays the free surface elevation for the three simulation methods. We can see that both the proposed method and the FBM capture the behaviour of the full 2D simulation, and from the right end of the channel, one can also see that the proposed method approximates the full 2D result better than the FBM. This is more obvious in figure 10, in which the time evolution of the free surface elevation at the probe points (indicated in figure 8), are plotted for all the simulation methods. Furthermore, figure 11 displays the plots of the time evolution of the - and - velocity components at selected probe points. One can see that the HCM performs better than the FBM at the indicated points. Of particular interest is the -velocity component at the probe point, which is located within the channel. We can see that while the FBM wrongly computed a zero -velocity all the time, the HCM computed the correct none-zero values with very good accuracy. It is thus very clear that the proposed method, HCM significantly outperforms the FBM for this test case.
In terms of efficiency, as shown in table 1, the full 2D simulation took 3,669 time steps and 3,100.31 seconds to complete this simulation while the coupling methods took less number of time steps and more than 50% reduced time to complete the same simulation. With the above observations, we conclude that it is possible to efficiently use coupling methods instead of full 2D simulation, and that the proposed method is capable of reproducing the full 2D solutions with greater accuracy than the FBM.
5.2 Test case 2 : Channel Flow into Elevated 2D Floodplain
This test case involves the same channel as in the previous example but connected to an elevated floodplain located in the region (see figure 12). The channel bed is flat and the floodplain bed is 0.5 meters high. The initial condition is the following.
[TABLE]
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 12, others are closed. Just like the previous test case, here nine probe points, are identified, see figure 12.
All the methods solved this problem with a grid of cells in the floodplain while the channel consists of , and cells for the full 2D, the HCM and the FBM respectively, see table 2. The simulation was run for ten seconds. Figures 13 and 14 show the free surface elevation and velocity magnitude after the last time step for each method. It can be seen that the HCM provides a better approximation of the full 2D results than the FBM. As a further validation of this claim, the time evolution of the free surface elevation is plotted in figures 15, while those of the -velocity and -velocity components are plotted in figure 16 for selected probe points. It can be seen that the horizontal coupling method is more accurate than the FBM at the points for all flow quantities and almost all the time. Again, the HCM really captures the flow structure of the full 2D simulation. This proves the accuracy of the proposed methods over the FBM for this test case.
For efficiency, we see from the processing time in table 2 that the FBM is very efficient but not very accurate while the horizontal coupling method is both very efficient and also has good accuracy.
5.3 Test case 3 : Flooding of an initially dry floodplain
The final test case involves the overflowing of a channel onto an initially dry floodplain. Both the channel and the floodplain are located in the domain, . The channel occupies the region, with flat bottom, , while the floodplain occupies the rest of the domain, , where . The bottom topography of the entire domain is the following
[TABLE]
where
[TABLE]
is the elevation of the channel wall, see figures 17(a) and 17(b) for the plots of and .
The initial condition consists of stationary water of depth, 0.08 meters in the channel and dry floodplain. The boundary conditions are time-dependent water depth at the left boundary of the channel and zero velocity at the right channel boundary, namely
[TABLE]
where
[TABLE]
where and . is a constant initial free-surface elevation inside the channel. The remaining boundaries are closed and the manning coefficients are the same as used in the previous cases. The following probe points are chosen, , , , , , , , , , , , , , and .
Table 3 shows the domain discritization for both the channel and the floodplain for each method being discussed. As before, all methods use the same grid for the floodplain but different grids for the channel. This problem was simulated for seconds. We report, in figures 18 - 21, the results of the simulation after 40 seconds and in figures 22 and 23, we report the results at selected probe points throughout the duration of the simulation.
As can been seen from the pictures, both coupling methods provide very good approximation of full 2D simulation results for both the free surface elevation (figure 18), the velocity components (figures 19 and 20) and the velocity magnitude (figure 21) for this test case. And in terms of accuracy of y-velocity component, the HCM provides better approximations as can be seen in figure 20.
To further understand the results of the simulations, the time evolution of the flow quantities at the probe points, have been examined. Here we report the results at the probe points and which are in the channel and the points, and in the floodplain. Figures 22 and 23 show the results for the selected points in the channel and floodplain respectively. In each figure, the left column displays the water depth, the second (middle) column shows the component of velocity, while the third(right) column shows the component of velocity.
From figure 22, we can see that both coupling methods provide very good approximation of the results of the full 2D simulations, especially for the water depth and component of velocity. However, only the HCM is able to compute the variation in the velocity and it does so with very good accuracy, see in figure 22. This further verifies the ability of of the HCM to compute the lateral discharges within the channel.
From figure 23, we also see that for the points in the floodplain, the coupling methods computed very good approximations of results of the full 2D simulation with the HCM computing more accurate results especially for the velocity. This figure also verify the no-numerical flooding property of the methods. That is, the floodplain initially remained dry until the time when water height rose above the channel banks. This is the reason why, for all points in the floodplain, the water depth and velocity remained at zero for the first several seconds of the simulation. Another thing to note is that due to the time-dependent boundary condition for this problem, water flowed onto the floodplain and after some time the water level in the channel decreased, hence the water in the floodplain drains back into the channel leaving the floodplain dry again. The coupling methods truly capture this phenomenon as one can see in figure 23 where the water depth and velocity return to zero towards the end of the simulation and remain at zero throughout the rest of the simulation. This is true for all the points in the floodplain, even those not reported here.
6 Conclusion
A horizontal coupling method has been proposed, implemented and tested in this paper. It presents a strategy to overcome the difficulty in computing the channel lateral discharges, circumvent the 1D assumption on the channel lateral discharge during flooding and propose a variant of the coupling term of [12] without the use of or imposing any restriction on the channel width variation. Numerical experiments show that the method computes adequate results. Particularly, the channel lateral discharges are properly computed without adopting complicated/iterative procedures. Finally, we note that for all the numerical test cases considered in this paper, the HCM would coincide with the FBM if the lateral discharges were not computed in the HCM. Therefore, the improved solution observed in the HCM over the FBM, for these test cases, is a result of the lateral discharges that are computed in the HCM. We, therefore, conclude that properly computing and restoring the channel lateral discharge, improves the quality of the computed solution and this can be done without introducing much computational overhead.
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. 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.
- 3[3] F. Bouchut. Efficient Numerical Finite Volume Schemes for Shallow Water Models . Elsevier, 2007.
- 4[4] 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.
- 5[5] J. A. Cunge, F. M. Holly, and A. Verwey. Practical aspects of computational river hydraulics . Pitman publishing, 1980.
- 6[6] E. D. Fernandez-Nieto, J. Marin, and J. Monnier. Coupling superposed 1d and 2d shallow-water models: Source terms in finite volume schemes. Computers & Fluids , 39(6):1070–1082, 2010.
- 7[7] R. Ghostine, I. Hoteit, J. Vazquez, A. Terfous, A. Ghenaim, and R. Mose. Comparison between a coupled 1d-2d model and a fully 2d model for supercritical flow simulation in crossroads. Journal of Hydraulic Research , 53(2):274–281, 2015.
- 8[8] N. Goutal, M. Parisot, and F. Zaoui. A 2d reconstruction for the transverse coupling of shallow water models. International Journal for Numerical Methods in Fluids , 75(11):775–799, 2014.