Preprint
Article
Numerical Simulation of Sediment Transport in Unsteady Open Channel Flow
Altmetrics
Downloads
161
Views
58
Comments
0
A peer-reviewed article of this preprint also exists.
This version is not peer-reviewed
Submitted:
06 June 2023
Posted:
07 June 2023
You are already at the latest version
Alerts
Abstract
This paper presented a two-dimensional, well-balanced hydrodynamic and sediment transport model based on the solutions of variable density shallow water equations (VDSWEs) for sediment-laden flows, and the Exner equation for bed changes. Those equations are solved in a coupled way by the first-order Godunov-type finite volume method. The Harten-Lax-van Leer-Contact (HLLC) Riemann solver is extended to find the local Riemann fluxes in order to maintain the exact balance between the momentum term and the bed slope term. The advantage of a well-balanced scheme over an unbalanced scheme is demonstrated by the synthetic standing contact-discontinuity test case. Then, the model is employed to simulate two laboratory experiments. At last, a field case, the 1996 Lake Ha! Ha! flood event (Canada), is simulated. Results of cross sectional geometries and profiles of longitudinal thalweg agree well with measurements. The accuracy and simplicity of the numerical model, together with the robust implementation, make the model a good candidate for practical engineering applications.
Keywords:
Subject: Engineering - Civil Engineering
1. Introduction
For dam-break flows over mobile beds, spatial variations of densities of sediment-laden flow is considerable and may impose significant impacts on flow fields (Cao et al. 2004). Since flow densities are not constants but spatial functions, dam-break flows over mobile beds are often modeled by the variable density shallow water equations (VDSWEs) (Cao et al. 2004; Wu & Wang 2007; Rosatti et al. 2008; Leighton et al. 2010; Begnudelli & Rosatti 2011; Wu et al. 2011; Li et al. 2013; Guan et al. 2014). As the name implies, the VDSWEs is an extension of the shallow water equations (SWEs) and is also based on the shallow water assumption (Vreugdenhil 1994). Besides, the following simplifications are employed in the derivation of the VDSWEs (Fraccarollo & Capart 2002; Cao et al. 2004): (1) the flow is dominated by suspended sediment load; (2) there are no obvious lags between sediment particles and water.
To account for its hyperbolic nature, the VDSWEs is often solved by the Godunov-type finite volume method (LeVeque 2002; Toro 2009). A thorough review about the application of the method to the numerical solution of the SWEs is given by Toro and Garcia-Navarro (2007). Being different from the SWEs, density of sediment-laden flow in the VDSWEs is an independent variable. To simplify the solution of the VDSWEs, Brufau et al. (2000) and Cao et al. (2004) reformulated the VDSWEs by moving the density-related terms to the right-hand side of the equations. After the reformulation, the left-hand side of the equations has the same form as the SWEs, while the right-hand side has two additional terms than the SWEs. Cao et al. (2004) pointed out that one additional term represents the effect of spatial variation of flow densities and the other one described the momentum exchange between flow and bed. Since the reformulated VDSWEs has the same left-hand side as the SWEs, it is a common practice to apply the previously well-developed Godunov-type SWEs solvers to solve the VDSWEs and to discretize the additional terms by a central difference scheme (Cao et al. 2004; Wu et al. 2011; Li et al. 2013; Guan et al. 2014).
However, one of the two additional terms on the right-hand side involves the gradient of density/concentration. Mathematically, the term represents a product of the Heaviside function by the Dirac function. Across discontinuous fronts of a density field, its mathematical meaning is not well-defined (Pares 2006; Castro et al. 2007). Numerically, this term originates from the advection term of the VDSWEs and it should to be discretized by an upwind scheme. Cozzolino et al. (2013) argued that it is difficult to maintain the conservative property (C-property) (Bermudez & Vazquez 1994) of numerical models based on the reformulated VDSWEs. Cozzolino et al. (2013) also suggested that the better way to create a well-balanced model is to solve the original formulation of the VDSWEs.
The objective of this paper is to present a novel well-balanced numerical model for the simulation of dam-break flows over mobile bed. The model is based on the original formulation of the VDSWEs. In the model, the Godunov-type finite volume method is applied to solve the system of equations simultaneously. As far as we know, the model is the first two-dimensional model that solves the VDSWEs using a well-balanced numerical scheme.
The rest of the paper is organized as follows: in Section 2, the governing equations, i.e. the two-dimensional variable density shallow water equations, and its two formulations are presented. In Section 3, the well-balanced numerical scheme is developed for the VDSWEs. In Section 4, the proposed model is tested against three test cases and the results are discussed and compared with observations. Finally, the main conclusions are drawn in Section 5.
2. Governing Equations
In this study, the dam-break flows over mobile beds are described by the variable density shallow water equations (VDSWEs). In the VDSWEs, the mixture of water and sediment is treated as a continuum, and both phases are moving together at the same velocity. The system of governing equations consists of volume, mass and momentum conservation equations for flow and sediment. For a sediment-laden flow, the bulk volume conservation equation is written as:
(1)where 
is time; 
and 
are the spatial coordinates; 
is the flow depth; 
and 
are the depth-averaged flow velocities in 
and 
directions, respectively; 
is the sediment exchange rate between flow and bed. Since the density of a sediment-laden flow is a spatial function, the mass conservation equation of is expressed as:
is time; and are the spatial coordinates; is the flow depth; and are the depth-averaged flow velocities in and directions, respectively; is the sediment exchange rate between flow and bed. Since the density of a sediment-laden flow is a spatial function, the mass conservation equation of is expressed as:
(2)where 
is the density of the sediment-laden flow; 
is the density of the saturated mobile bed; 
is the porosity of the mobile bed; 
is the density of the water; 
is the density of the bed material. The concentration of sediment-laden flow can be calculated by 
. The depth-averaged momentum conservation equations in 
and 
directions are given as:
is the density of the sediment-laden flow; is the density of the saturated mobile bed; is the porosity of the mobile bed; is the density of the water; is the density of the bed material. The concentration of sediment-laden flow can be calculated by . The depth-averaged momentum conservation equations in and directions are given as:
(3)
(4)where 
is the gravitational acceleration; 
and 
are the 
and 
components of bed slope; 
is the elevation of the immobile bed; 
is the depth of the mobile bed; 
is the drag coefficient; 
is Manning’s roughness coefficient; 
is the magnitude of the flow velocity; 
is the velocity vector.
is the gravitational acceleration; and are the and components of bed slope; is the elevation of the immobile bed; is the depth of the mobile bed; is the drag coefficient; is Manning’s roughness coefficient; is the magnitude of the flow velocity; is the velocity vector.The system of governing equations can also be written in a vectorial form as:
(5)where 
and 
are the vectors of primitive variables and conservative variables:
and are the vectors of primitive variables and conservative variables:
(6), 
and 
are the vectors of advective fluxes:
and are the vectors of advective fluxes:
(7), and 
, 
and 
are the bed slope, bed friction, and bed material source terms:
, and are the bed slope, bed friction, and bed material source terms:
(8)Equation (5) represents a system of time dependent, nonlinear hyperbolic partial differential equations. This system may result in sharp and discontinuous solutions even if starting from continuous initial conditions.
2.1. Nonequilibrium Sediment Transport Model
The evolution of mobile bed is described by the following conservation equation of bed material:
(9)The mass exchange between sediment-laden flow and mobile bed is evaluated by the nonequilibrium sediment transport equation (Armanini & Di Silvio 1988; Wu & Wang 2007). It reads:
(10)where 
is the nonequilibrium adaption length of sediment transport; 
is the actual flux of sediment transport; and 
is the sediment transport capacity or the equilibrium sediment transport rate. On account of the recommendations of El Kadi Abdeerezzak and Paquier (2010) and Van Emelen et al. (2014), the sediment transport capacity is calculated by the modified Meyer-Peter-Müller (MPM) formula:
is the nonequilibrium adaption length of sediment transport; is the actual flux of sediment transport; and is the sediment transport capacity or the equilibrium sediment transport rate. On account of the recommendations of El Kadi Abdeerezzak and Paquier (2010) and Van Emelen et al. (2014), the sediment transport capacity is calculated by the modified Meyer-Peter-Müller (MPM) formula:
(11)where 
is the Shields number; 
is the friction velocity; 
is the critical Shields number.
is the Shields number; is the friction velocity; is the critical Shields number.The nonequilibrium adaption length denotes the spatial lag between actual sediment transport rate and its saturation/equilibrium rate. It is calculated by Wu (2004):
(12)where 
is the adaption length of bed-load; 
is the adaption coefficient of suspended-load; 
is the settling velocity of a single sediment particle. Since there have been no consensus on 
and 
, both parameters are treated as calibration parameters in the proposed model. The settling velocity of sediment particle is calculated by Cao (1999):
is the adaption length of bed-load; is the adaption coefficient of suspended-load; is the settling velocity of a single sediment particle. Since there have been no consensus on and , both parameters are treated as calibration parameters in the proposed model. The settling velocity of sediment particle is calculated by Cao (1999):
(13)where 
is the kinematic viscosity of water; 
is the specific gravity of submerged sediment particle; 
is the median diameter of sediment particle.
is the kinematic viscosity of water; is the specific gravity of submerged sediment particle; is the median diameter of sediment particle.2.2. Numerical Methods
2.2.1. Godunov-type Finite Volume Method
The integral form of the VDSWEs over a computational cell can be written as:
(14)where 
represents a computational cell; 
is the boundary of the cell; 
is the matrix of fluxes; 
is the outward unit normal vector; 
and 
are the components of the normal vector in 
and 
directions, respectively; and 
is the fluxes normal to the cell boundary designated by 
. In this study, the fully coupled system of governing equations is solved by the first-order Godunov-type finite volume method over a regular Cartesian mesh. For a rectangular cell, the fluxes across the cell boundaries are calculated by
represents a computational cell; is the boundary of the cell; is the matrix of fluxes; is the outward unit normal vector; and are the components of the normal vector in and directions, respectively; and is the fluxes normal to the cell boundary designated by . In this study, the fully coupled system of governing equations is solved by the first-order Godunov-type finite volume method over a regular Cartesian mesh. For a rectangular cell, the fluxes across the cell boundaries are calculated by
(15)where subscript 
is the edge index; subscript 
and 
denote the left and right edges of a cell boundary, respectively; 
is the Riemann fluxes at the edge 
; 
and 
are the vectors of conservative variables reconstructed at the left and right edges, respectively; 
is the unit normal vector for the kth edge; 
is the length of the kth edge.
is the edge index; subscript and denote the left and right edges of a cell boundary, respectively; is the Riemann fluxes at the edge ; and are the vectors of conservative variables reconstructed at the left and right edges, respectively; is the unit normal vector for the kth edge; is the length of the kth edge.A first-order, two-step fractional scheme is employed in the numerical model to update the solution at each time step:
Step one:
(16)Step two:
(17)where superscript 
denotes the solution at time step 
; superscript 
denotes the time step 
; superscript 
denotes the intermediate solution between time step 
and 
.
denotes the solution at time step ; superscript denotes the time step ; superscript denotes the intermediate solution between time step and .Time steps are limited by two stability conditions. The first is the Courant-Friedrichs-Lewy (CFL) condition that requires the maximum CFL number should not be greater than 0.5 (LeVeque 2002) at each time step. The second is the bed change condition that limits the change of mobile bed at each time step up to 10% of the local bed depth (Wu et al. 2011).
2.2.2. HLLC Approximated Riemann Solver
The HLLC (Harten-Lax-van Leer-Contact) approximate Riemann solver (Toro 2009) is extended to calculate numerical fluxes across cell boundaries. The wave structure of a typical HLLC solution is demonstrated in Cozzolino et al. (2013). The solution is separated by three waves: the slowest wave (
), the middle wave (
), and the fastest wave (
). There are four constant states in the solution: the left state (
), the right state (
), the left star state (
), and the right star state (
). The corresponding HLLC numerical flux 
is defined as:
), the middle wave (), and the fastest wave (). There are four constant states in the solution: the left state (), the right state (), the left star state (), and the right star state (). The corresponding HLLC numerical flux is defined as:
(18)where 
represents the supercritical flow from the left to the right; 
represents the supercritical flow from the right to the left; 
and 
are the fluxes at the star regions. By applying the Rankine-Hugoniot conditions across the waves, 
and 
, the numerical fluxes at the star regions can be calculated by:
represents the supercritical flow from the left to the right; represents the supercritical flow from the right to the left; and are the fluxes at the star regions. By applying the Rankine-Hugoniot conditions across the waves, and , the numerical fluxes at the star regions can be calculated by:
(19)The middle wave speed is estimated by Cozzolino et al. (2013):
(20)where 
and 
are the hydrostatic pressures at the left and right edges, respectively; 
and 
are the normal velocities at the left and right edges, respectively. The left and right states at the star regions can be determined by:
and are the hydrostatic pressures at the left and right edges, respectively; and are the normal velocities at the left and right edges, respectively. The left and right states at the star regions can be determined by:
(21)
(22)The left and right wave speeds are estimated by:
(23)where 
and 
are the Roe averaged flow depth and velocity.
and are the Roe averaged flow depth and velocity.2.2.3. Well-balanced Scheme for the VDSWEs
Without any loss of generality, only one-dimensional equations are considered in this section. For a stationary flow (
), the momentum equation of SWEs becomes:
), the momentum equation of SWEs becomes:
(24)If a numerical scheme can maintain this stationary state, the scheme is called a well-balanced scheme (Greenberg & Leroux 1996), and also it is said to satisfy the conservation property (or C-property) (Bermudez & Vazquez 1994).
As to the VDSWEs, the momentum equation changes into:
(25)It is obvious that the definition of a well-balanced scheme or the conservative property can be extended to the VDSWEs, if a numerical scheme can maintain Equation (25). Rearranging Equation (25) into a similar form as Equation (24):
(26)Since the density gradient term on the right hand side comes from the advection term of Equation (25), the term bears an upwind characteristics. It is difficult to maintain the conservative property of a numerical scheme if this term is discretized by a central difference method.
To create a well-balanced scheme, the hydrostatic reconstruction approach (Audusse et al. 2004) for the SWEs is incorporated into the proposed model. According to the approach, the bed elevation at the cell interface (
) between the cell (
) and the cell (
) is evaluated by
) between the cell () and the cell () is evaluated by
(27)where 
is the bed elevation; subscript 
and 
are the cell indices in 
and 
directions; subscript 
represents the cell boundary. Then, flow depths at the left and right sides of a cell boundary are reconstructed as:
is the bed elevation; subscript and are the cell indices in and directions; subscript represents the cell boundary. Then, flow depths at the left and right sides of a cell boundary are reconstructed as:
(28)The reconstructed conservative variables at the cell boundary equal to:
(29)Accordingly, the bed slope term at the cell center is discretized as:
(30)where subscript E, W, N, and S represent the east, west, north, and south sides of a cell.
To prove the well-balanced property of the scheme, the numerical discretization to a stationary flow is presented. For a flow in hydrostatic equilibrium state, the wave speeds are:
(31)For the cell i, the left-hand side (LHS) of Equation (25) is discretized as:
(32)while the right-hand side (RHS) of Equation (25) equals:
(33)Since the left hand side is exactly equal to the right hand side, the proposed numerical scheme is a well-balanced scheme.
3. Results
The numerical method was tested by using two cases: one is a synthetic case of standing contact discontinuity, and the other is dam break flow with sediment concentration gradient.
3.1. Standing Contact-Discontinuity
The aim of this synthetic case is to demonstrate the well-balanced property of the model. The case is proposed by Cozzolino et al. (2013): assuming a stationary fluid in a flume (L = 500 m) with a flat bed and both sides of the flume are solid walls, the initial conditions is given as:
(34)with:
(35)For the given initial condition, the hydrostatic pressure is constant through the flume. The discontinuity at x0 = 250 m is a standing contact-discontinuity and the corresponding solution is of fluid at stationary state.
The numerical simulation is carried out over a Cartesian mesh with Δx =1 m. For the simulation, the time step is Δt = 0.02 s and the simulation time is tmax =10 s. The simulation results, including surface profile, density profile, and velocity profile, are shown in Figure 1. The profiles of the results are identical to the stationary solution of this test. Therefore, the well-balance property of the numerical model is confirmed by this test case.
Besides, the results of an unbalanced scheme (Cao et al. 2004) are also plotted in Figure 1. It is apparent that there are two oscillations propagate toward the upstream and downstream ends of the flume and the unbalanced scheme cannot keep the stationary state solution. These comparisons demonstrate the advantage of a well-balanced scheme over an unbalanced scheme: a well-balanced scheme can maintain the exact equilibrium between the hydrostatic pressure term and the bed slope term and there are no oscillations in the solution of a well-balanced scheme.
3.2. Case 1: One-dimensional Dam-break Flow over Mobile Bed
In this test case, the model is applied to a one-dimensional laboratory experiment (Spinewine & Zech 2007). The objective of this experiment is to investigate the erosional behavior of dam-break flow over mobile bed. The setup of the experiment is shown in Figure 2. The flume is 6 m long, 0.25 m wide, and 0.7 m high. The flume is equipped with a fast downward-moving gate, which is used to simulate idealized dam-break events. The experiment uses a flat, loose granular sediment bed with following parameters: the density of sediment particle, ps = 2683 kg/m3; the median diameter of sand, d50 = 1.82 mm; the settling velocity, 0 = 0.16 m/s; and the porosity, = 0.47. The Manning’s roughness coefficient was estimated to be n = 0.0165 s/m1/3.
Initially, the reservoir at the upstream is filled with water to a depth, h = 0.35 m. The downstream of the flume is dry bed. The mobile bed layer consists of fully saturated sands with an initial thickness, b = 0.1 m. As to the boundary conditions, the upstream boundary condition is wall boundary, and the transmissive boundary condition is applied to the downstream boundary. At time t = 0 s, the gate was suddenly removed to create a dam-break flow. During the experiment, the dam-break flow is recorded by several high speed digital cameras. The total experimental time, tmax = 1.5 s.
The numerical simulation is carried out with the following parameters: Δx = 0.06 m and Δt = 0.001 s. In Figure 3, the calculation results are compared with the measured data at T= 0.25 s, 0.50 s, 0.75s, 1.00 s, 1.25 s, and 1.50 s. An overall comparison shows reasonably well agreements between the calculated results and the measurements. Not only the propagation of the shock-front is accurately predicted, but also the temporal evolution of free surface agrees well with the interface tracking results of the digital cameras. Although the bed changes are not significant in the experiment, the magnitudes of the simulated bed changes are consistent with the measurements.
To quantitatively evaluate the accuracy of the developed model, the RMSEs and NRMSEs of the simulated results are calculated and shown in Table 1. As shown in Table 1, the RMSEs of the simulated surface profiles are around 0.01 m, while the RMSEs of the simulated bed profiles are less than 0.05 m. However, the NRMSEs of the simulated bed profiles are much greater than the simulated surface profiles. This reveals that the simulated surface profiles are more accurate than the simulated bed profiles.
The simulated concentration profiles at different time instants are plotted in Figure 4. The existence of high concentration at the shock front (about 20%) justifies the significance of the variable density effect. This is also a challenge to the numerical modeling of sediment transport because the concentration profile is discontinuous at the shock wave front. Figure 4 demonstrated that the proposed model is able to capture this concentration discontinuity without introducing spurious oscillations.
3.3. Case 2: Two-dimensional dam-break flow over mobile bed
This laboratory experiment was conducted at the Hydraulics Unit of the Mechanical and Civil Engineering Laboratory, Universite catholique de Louvain, Belgium (Soares-Frazao et al. 2012). The experiment aims to provide a benchmark test case to validate numerical models for the simulation of dam-break flow over mobile bed (Wu et al. 2011; Canelas et al. 2013; Li et al. 2013). The experiment is a dam-break flow from an upstream reservoir flowing into a flume over a mobile bed made of uniform coarse particles. A schematic view of the flume is shown in Figure 5. The length of the flume is 36 m; the width of the flume is 3.6 m. The breached dam is located in a narrow reach (1 m long and 1 m wide) between two impervious blocks. The upstream of the flume was a reservoir, whereas the downstream is a flooded channel. The origin of the coordinates is situated at the center of the dam. The Manning’s coefficient for the sandy bed is given as .
The bottom of the flume was covered with a layer of saturated sand, which expanded from 1 m upstream of the dam to 9 m downstream. The thickness of mobile sand layer is 8.5 cm. The related properties of the sand layer are: the mean diameter of the particles, d50 = 1.61 mm; the porosity of bed material, ε = 0.42; the density, ps = 2630 kg/m3; and the critical Shields parameter, = 0.047 . The Manning’s roughness coefficient for this sandy channel is measured as n = 0.0165 s/m1/3.
The dam-break flow was triggered by rapidly lifting the gate separating the reservoir and the channel. The initial water level in the reservoir is 0.47 m, while the downstream channel is initially dry (h=0 m). The wall boundary condition is applied to the upstream end of the reservoir, and similarly, the transmissive boundary condition is imposed at the downstream outlet of the flume. The experiment lasted for 20 s. During the experiment, eight sonic water level gauges were used to record water levels. The bed profiles were measured after the experiment along three longitudinal lines: y = 0.2 m, y = 0.7 m, and y = 1.45 m.
For the numerical simulation, the flume is discretized using a Cartesian mesh (Δx = Δy = 0.1 m), with a total of 10,602 rectangular cells. The time step is set to be Δt= 0.01 s. The comparisons of simulated and measured water level hydrographs at eight gauges are shown in Figure 6. The quantitative evaluation indices, i.e., the RMSEs and the NRMSEs, are shown in Table 2. The comparisons showed that, at all eight gauges, both the arrival times of the dam-break waves and the peak water levels are accurately predicted. At the early stage of experiment, the water levels at the gauges, G2 and G3, are overpredicted, while the water levels are underpredicted at the gauges, G5 and G8. At the late stage, the water levels are underpredicted at G1, G4, G6, G7, and G8. For the symmetric gauge pairs (G1/G4, G2/G3, G5/G8, and G6/G7), the simulated hydrographs are nearly symmetrical except for the results at G5 and G8 at the late stage. The RMSEs of the predicted water level hydrographs are between 0.015 m and 0.02 m, while the NRMSEs are between 12% and 22%. The NRMSE results of the upstream gauges (G1, G2, G3, and G4) are better than the downstream gauges (G5, G6, G7, and G8). A possible explanation is that, as the dam-break flow propagating from upstream to downstream, three dimensional flow effects becomes more dominant.
The comparisons between the calculated and measured bed profiles along the longitudinal lines at the end of the experiment are shown in Figure 7. The RMSEs and the NRMSEs are shown in Table 3. The magnitudes of RMSEs are between 0.009 m and 0.017 m, while the NRMSEs’ are between 14% and 33%. The simulated mean bed profiles matched the measurements at all three measured longitudinal lines. Although the calculated and the measured bed profiles have the same trends, the calculated bed profiles are smoother than the measured bed profiles. This implies that the developed model cannot capture the local perturbations generated by 3D flow effects (Canelas et al. 2013). These perturbations are caused by the flow in the vertical direction, and cannot be captured by a 2D depth-averaged model. On the other hand, the energy dissipations caused by these perturbations can be calculated by turbulence modeling method (Yu & Duan 2012). It thus follows that further researches are necessary for incorporating turbulence model into the proposed model.
3.4. Case 3: 1996 Lake Ha! Ha! catastrophic flood event
In this case, the performance of the proposed model is tested against a field event: the 1996 Lake Ha! Ha! catastrophic flood event, in the Saguenay region of Quebec, Canada (Figure 8). A detailed description of this flood event and extensively documented data are provided by Capart et al. (2007). Brooks and Lawrence (1999) gave a detailed geomorphic description about the event. Ferreira et al. (2009) and El Kadi Abderrezzak and Paquier (2009) carried out a one-dimensional numerical simulation of this event.
From July 18 to 21, 1996, an extreme precipitation event affected the Saguenay region of Quebec, Canada. At the Ha! Ha! Lake, an earthfill dyke was being overtopped by up to 0.26 m of water, and a new outlet channel formed. Overall, of water was estimated to have drained from the lake. The failure of the dyke resulted in a peak discharge of 8 times the 100-year flood. The Ha! Ha! River was severely damaged by the resulting flood flow (Brooks & Lawrence 1999).
The numerical simulation started with the digital elevation model (DEM) of the Ha! Ha! River, which was surveyed in May 1994 (Capart et al. 2007). The spatial data is based on the Modified Transverse Mercator (MTM) projection, zone 7 coordinates (NAD83). The spatial ranges of the DEM data are: 275,000 m x 282,000 m on the east-west direction, and 5318000 mN y 5354000 mN. In addition to the DEM data, the geometry data of evenly spaced cross sections are also provided. These 363 cross sections spaced at 100 m interval, and are oriented approximately normal to the central line of the Ha! Ha! River. Each cross section is identified by its streamwise distance measured from the failed dyke. In spite of abundant data, it has to be pointed out that the average error of the preflood riverbed elevations is estimated to be about 2 m (El Kadi Abderrezzak & Paquier 2009).
Based on the field observation and photo interpretation, superficial materials exposed along the channel and floodplains are classified as noncohesive sediment (sand and gravel), cohesive sediment (glaciomarine or glaciodiamicton), or bedrock (Capart et al. 2007). Due to the lack of sufficient field data, the influence of cohesion on sediment transport was not considered in the numerical model. According to El Kadi Abderrezzak and Paquier (2009), the median diameter of bed material for the whole Ha! Ha! River is equal to 0.5 mm, the density of bed material is ps = 2650 kg/m and the porosity equals to 0.4. Besides, the influence of bedrock on the channel profile was studied in Capart et al. (2007), the bedrock constraint was also imposed in this simulation for which the outcrops of bedrocks or coarse glacial deposits are assumed to act as rigid, non-erodible beds.
Initially, the Ha! Ha! River is assumed to be dry everywhere. The discharge hydrograph shown in Figure 9 is used as the upstream inflow boundary condition, and the inflow is assumed to be clear water (El Kadi Abderrezzak & Paquier 2009). The downstream boundary at the Ha! Ha! Bay is set to satisfy the transmissive boundary condition. The overtopping occurred at 14:00 on July 19, 1996. The lake is estimated to be emptied in 18 hours. During the event, the water level of the lake dropped from 381 m to 370 m. The total simulation time is set to be 24 hours. For the numerical simulation, the river is discretized by a Cartesian mesh with Δx = Δy = 20m. The time step size is Δt = 0.5 s.
Comparisons of the selected cross sections along the river are shown in Figure 10. The cross sectional changes are reasonably well predicted by the model. The RMSEs and NRMSEs of bed elevation changes are calculated for all cross sections, and are shown in Figure 11. The average of RMSEs is 4.4 m, and the average of NRMSEs is 49.53%. As shown in Figure 12, among all the data, 50% of the RMSEs are less than 3.3 m, and 90% of the RMSEs are less than 7.9 m. As to the NRMSEs, about 50% cross sections have NRMSE values less than 40%, and 90% of the NRMSEs are less than 70%.
In the longitudinal direction, the measured and calculated thalwegs are plotted in Figure 13. The RMSE of thalweg change is 5.3 m, and the NRMSE is 20.18%. It is obvious that the model performed better in the longitudinal direction than in the cross sectional direction. Also, the simulated thalweg profile at the middle reach is overestimated, especially at the 23 km from the channel mouth. According to El Kadi Abderrezzak and Paquier (2009), this part of the channel is controlled by outcrops of bedrocks. During the 1996 flood event, a new reach formed at the right floodplain and the bed was eroded up to 20 m there. The proposed model failed to predict this geomorphic change.
For this real-world flood event, although the RMSEs of bed elevation (in meter scale) are two orders of magnitude greater than the laboratory cases (in centimeter scale), the NRMSEs of bed elevation (49.53% in horizontal direction and 20.18% in longitudinal direction) are at the same order of magnitude as the laboratory cases (24.68% for 1D dam-break flow case, 20.96% for 2D dam-break flow case). Since there is no consensus on which sediment transport equation is most suitable for a specific river reach, an appropriate calibration procedure by using different bed load transport equation and varying Manning’s roughness coefficient is needed to reach better matches of modeling results with observations. The focus of this paper is the new well-balanced numerical scheme for solving the VDSWEs, this calibration procedure was not performed for each testing case. Nevertheless, the simulating results of bed elevation changes satisfactorily match both the laboratory and field observations. Therefore, the developed well-balanced numerical scheme for solving VDSWEs is proved to be a robust and accurate method for simulating dam-break flows over mobile beds.
4. Conclusions
A two-dimensional, well-balanced, nonequilibrium sediment transport model for the simulation of dam-break flows over mobile beds is developed and tested. Based on the VDSWEs, the model simulates the development of dam-break flows, sediment transport, and bed elevation changes. In comparison with previously developed 2D models, the governing equations of the proposed model are based on the original formulation of VDSWEs, thus it is simple to preserve the conservative property of sediment-laden flow. In the model, the system of governing equations is solved by a fully coupled first-order Godunov-type finite volume method. The HLLC Riemann solver is extended to the two-dimensional VDSWEs to calculate the Riemann fluxes across cells’ boundaries. The model also features a well-balanced scheme that maintains the exact balance between the momentum term and the bed slope term.
The performance of the model is verified by four test cases. A synthetic case with analytic solution is used to verify the well-balance property of the model. Two laboratory experimental studies of 1D and 2D dam-break flows over mobile beds show the accuracy of the model for reproducing not only bed profiles but also the flooding processes. The last case is the 1996 Lake Ha! Ha! flood event. The model predicts reasonable bed cross sections and thalwegs for this field case. Also, the robustness of the model is demonstrated by this field case. The accuracy and simplicity of the proposed model, together with the robust implementation of well-balanced numerical scheme, make this model suitable for practical hydraulic engineering applications.
Acknowledgments
The research is partially funded by NSF Award EAR-0846523 to the University of Arizona. The funding support is essential for authors to complete this research.
References
- Armanini, A. & Di Silvio, G. 1988 A one-dimensional model for the transport of a sediment mixture in non-equilibrium conditions. Journal of Hydraulic Research 26 (3), 275-292.
- Audusse, E. , Bouchut, F., Bristeau, M. O. Klein, B. & Perthame, B, T. 2004 A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. Siam Journal on Scientific Computing 25 (6) 2050-2065.
- Begnudelli, L. & Rosatti, G. 2011 Hyperconcentrated 1D shallow flows on fixed bed with geometrical source term due to a bottom step. Journal of Scientific Computing 48 (1), 319-332.
- Bermudez, A. & Vazquez M. E. 1994 Upwind methods for hyperbolic conservation-laws with source terms. Computers & Fluids 23 (8), 1049-1071.
- Brooks, G. R. & Lawrence, D. E. 1999 The drainage of the Lake Ha!Ha! reservoir and downstream geomorphic impacts along Ha!Ha! River, Saguenay area, Quebec, Canada. Geomorphology 28 (1) 141-168.
- Brufau, P. , Garcia-Navarro, P., Ghilardi, P., Natale, L. & Savi, F. 2000 1D mathematical modelling of debris flow. Journal of Hydraulic Research 38 (6), 435-446.
- Canelas, R. , Murillo, J. & Ferreira, R. 2013 Two-dimensional depth-averaged modelling of dam-break flows over mobile beds. Journal of Hydraulic Research 51 (4) 392-407.
- Cao, Z. 1999 Equilibrium near-bed concentration of suspended sediment. Journal of Hydraulic Engineering 125 (12), 1270-1278.
- Cao, Z. , Pender, G., Wallis, S. & Carling, P. 2004 Computational dam-break hydraulics over erodible sediment bed. Journal of Hydraulic Engineering 130 (7), 689-703.
- Capart, H. , Spinewine, B., Young, D. L., Zech, Y., Brooks, G. R., Leclerc, M. & Secretan, Y. 2007 The 1996 Lake Ha! Ha! breakout flood, Quebec: Test data for geomorphic flood routing methods. Journal of Hydraulic Research 45 (S1) 97-109.
- Castro, M. J. , Pardo Milanes, A. & Pares, C. 2007 Well-balanced numerical schemes based on a generalized hydrostatic reconstruction technique. Mathematical Model and Methods in Applied Sciences 17 (12), 2055-2113.
- Cozzolino, L. , Cimorelli, L., Covelli, C., Morte, R. D. & Pianese, D. 2013 Novel numerical approach for 1D variable density shallow flows over uneven rigid and erodible beds. Journal of Hydraulic Engineering 140 (3), 254-268.
- El Kadi Abderrezzak, K. & Paquier, A. 2009 One-dimensional numerical modeling of sediment transport and bed deformation in open channels. Water Resources Research 45 (5) W05404.
- El Kadi Abderrezzak, K. & Paquier, A. 2010 Applicability of sediment transport capacity formulas to dam-break flows over movable beds. Journal of Hydraulic Engineering 137 (2), 209-221.
- Ferreira R., M. , Franca, M. J., Leal, J. G. & Cardoso A. H. 2009 Mathematical modelling of shallow flows: Closure models drawn from grain-scale mechanics of sediment transport and flow hydrodynamics. Canadian Journal of Civil Engineering 36 (10) 1605-1621.
- Fraccarollo, L. & Capart, H. 2002 Riemann wave description of erosional dam-break flows. Journal of Fluid Mechanics, 461, 183-228.
- Greenberg, J. M. & Leroux, A. Y. 1996 A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. Siam Journal on Numerical Analysis 33 (1), 1-16.
- Guan, M. , Wright, N. G. & Sleigh, P. A. 2014 2D Process-Based Morphodynamic Model for Flooding by Noncohesive Dyke Breach. Journal of Hydraulic Engineering 140 (7), 04014022. 0140. [Google Scholar]
- Leighton, F. Z. , Borthwick, A. G. L. & Taylor, P. H. 2010 1-D numerical modelling of shallow flows with variable horizontal density. International Journal for Numerical Methods in Fluids 62 (11), 1209-1231.
- LeVeque, R. J. 2002 Finite volume methods for hyperbolic problems. Cambridge University Press, Cambridge, New York.
- Li, W. , de Vriend, H. J., Wang, Z. & van Maren, D. S. 2013 Morphological modeling using a fully coupled, total variation diminishing upwind-biased centered scheme. Water Resources Research 49 (6), 3547-3565.
- Pares, C. 2006 Numerical methods for nonconservative hyperbolic systems: A theoretical framework. Siam Journal on Numerical Analysis 44 (1), 300-321.
- Rosatti, G. , Murillo, J. & Fraccarollo, L. 2008 Generalized Roe schemes for 1 D two-phase, free-surface flows over a mobile bed. Journal of Computational Physics 227 (24), 10058-10077.
- Soares-Frazao, S. , Canelas, R., Cao, Z., Cea, L., Chaudhry, H. M., Die Moran, A., El Kadi, K., Ferreira, R., Cadorniga, I. F., Gonzalez-Ramirez, N., Greco, M., Huang, W., Imran, J., Le Coz, J., Marsooli, R., Paquier, A., Pender, G., Pontillo, M., Puertas, J., Spinewine, B., Swartenbroekx, C., Tsubaki, R., Villaret, C., Wu, W., Yue, Z. & Zech, Y. 2012 Dam-break flows over mobile beds: experiments and benchmark tests for numerical models. Journal of Hydraulic Research 50 (4) 364-375.
- Spinewine, B. & Zech, Y. 2007 Small-scale laboratory dam-break waves on movable beds. Journal of Hydraulic Research 45 (S1) 73-86.
- Toro, E. F. 2009 Riemann solvers and numerical methods for fluid dynamics: a practical introduction, 3rd edn. Springer, Dordrecht, New York.
- Toro, E. F. 2009 Riemann solvers and numerical methods for fluid dynamics: a practical introduction, 3rd edn. Springer, Dordrecht, New York.
- Van Emelen, S. , Zech, Y. & Soares-Frazao, S. 2014 Impact of sediment transport formulations on breaching modelling. Journal of Hydraulic Research 53 (1), 60-72.
- Vreugdenhil, C. B. 1994 Numerical methods for shallow-water flow. Kluwer Academic Publishers, Dordrecht, Boston.
- Wu, W. 2004 Depth-averaged two-dimensional numerical modeling of unsteady flow and nonuniform sediment transport in open channels. Journal of Hydraulic Engineering 130 (10), 1013-1024.
- Wu, W. , Marsooli, R. & He, Z. 2011 Depth-averaged two-dimensional model of unsteady flow and sediment transport due to noncohesive embankment break/breaching. Journal of Hydraulic Engineering 138 (6), 503-516.
- Wu, W. & Wang S. S. 2007 One-dimensional Modeling of dam-break flow over movable beds. Journal of Hydraulic Engineering 133 (1), 48-58.
- Yu, C. & Duan, J. 2012 Two-dimensional depth-averaged finite volume model for unsteady turbulent flow. Journal of Hydraulic Research 50 (6), 599-611.
Figure 1.
Solution of test case 1: (a) surface profile; (b) density profile; and (c) velocity profile.
Figure 1.
Solution of test case 1: (a) surface profile; (b) density profile; and (c) velocity profile.
Figure 2.
Experimental setup of test case 2.
Figure 3.
Measured and calculated results of test case 2: (a) t = 0.25 s; (b) t = 0.50 s; (c) t = 0.75 s; (d) t = 1.00 s; I t = 1.25 s; and (f) t = 1.50 s.
Figure 3.
Measured and calculated results of test case 2: (a) t = 0.25 s; (b) t = 0.50 s; (c) t = 0.75 s; (d) t = 1.00 s; I t = 1.25 s; and (f) t = 1.50 s.
Figure 4.
Concentration profiles of 1D dam-break flow.
Figure 5.
Experimental setup of 2D dam-break flow.
Figure 6.
Measured and calculated water levels of 2D dam-break flow at: (a) Gauge 1; (b) Gauge 2; (c) Gauge 3; (d) Gauge 4; (e) Gauge 5; (f) Gauge 6; (g) Gauge 7; and (h) Gauge 8.
Figure 6.
Measured and calculated water levels of 2D dam-break flow at: (a) Gauge 1; (b) Gauge 2; (c) Gauge 3; (d) Gauge 4; (e) Gauge 5; (f) Gauge 6; (g) Gauge 7; and (h) Gauge 8.
Figure 7.
Measured and calculated bed profiles of 2D dam-break flow: (a) y = 0.20 m; (b) y = 0.70 m; (c) y = 1.45 m.
Figure 7.
Measured and calculated bed profiles of 2D dam-break flow: (a) y = 0.20 m; (b) y = 0.70 m; (c) y = 1.45 m.
Figure 8.
DEM map of the Ha! Ha! River, Canada.
Figure 9.
Inlet discharge hydrograph of 1996 Lake Ha! Ha! flood event.
Figure 10.
Measured and calculated bed cross sections of Lake Ha! Ha! flood event: (a) cross section #20; (b) cross section #70; (c) cross section #120; (d) cross section #170; (e) cross section #220; (f) cross section #270; and (g) cross section #320.
Figure 10.
Measured and calculated bed cross sections of Lake Ha! Ha! flood event: (a) cross section #20; (b) cross section #70; (c) cross section #120; (d) cross section #170; (e) cross section #220; (f) cross section #270; and (g) cross section #320.
Figure 11.
Analysis of calculated bed cross sections of Lake Ha! Ha! flood event: (a) RMSEs; (b) NRMSEs.
Figure 11.
Analysis of calculated bed cross sections of Lake Ha! Ha! flood event: (a) RMSEs; (b) NRMSEs.
Figure 12.
Cumulative percentage distributions of Lake Ha! Ha! flood event: (a) RMSE; and (b) NRMSE.
Figure 12.
Cumulative percentage distributions of Lake Ha! Ha! flood event: (a) RMSE; and (b) NRMSE.
Figure 13.
Measured and calculated thalwegs of Lake Ha! Ha! flood event: (a) 0 – 12 km; (b) 12 – 24 km; (c) 24 – 36 km.
Figure 13.
Measured and calculated thalwegs of Lake Ha! Ha! flood event: (a) 0 – 12 km; (b) 12 – 24 km; (c) 24 – 36 km.
Table 1.
RMSE and NRMSE of 1D dam-break flow.
| Time (s) | Water Level | Bed Elevation | ||
|---|---|---|---|---|
| RMSE (cm) | NRMSE (%) | RMSE (cm) | NRMSE (%) | |
| 0.25 | 1.13 | 3.35 | 0.45 | 23.20 |
| 0.50 | 0.57 | 1.80 | 0.38 | 25.09 |
| 0.75 | 0.78 | 2.79 | 0.39 | 26.47 |
| 1.00 | 0.84 | 3.37 | 0.33 | 22.24 |
| 1.25 | 0.82 | 3.58 | 0.35 | 23.75 |
| 1.50 | 0.87 | 4.83 | 0.32 | 27.35 |
| Average | 0.84 | 3.29 | 0.37 | 24.68 |
Table 2.
RMSE and NRMSE of water levels of 2D dam-break flow.
| Gauge # | RMSE (cm) | NRMSE (%) |
|---|---|---|
| 1 | 1.71 | 15.11 |
| 2 | 2.07 | 12.66 |
| 3 | 2.17 | 12.72 |
| 4 | 1.59 | 13.81 |
| 5 | 1.39 | 18.90 |
| 6 | 1.49 | 20.73 |
| 7 | 1.58 | 21.91 |
| 8 | 1.98 | 21.50 |
| Average | 1.75 | 17.17 |
Table 3.
RMSE and NRMSE of bed elevations of 2D dam-break flow.
| Bed Profile (m) | RMSE (cm) | NRMSE (%) |
|---|---|---|
| 0.20 | 1.24 | 13.03 |
| 0.70 | 1.67 | 32.65 |
| 1.45 | 0.91 | 17.20 |
| Average | 1.27 | 20.96 |
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content. |
© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
MDPI Initiatives
Important Links
© 2024 MDPI (Basel, Switzerland) unless otherwise stated