Abstract
This study models a dam-break flow over a bed by using a depth-averaged numerical model based on finite-volume method and computes the dam-break flow and bed morphology characteristics. The generalized shallow water equations considering the sediment transport and bed change on dam-break flow are adopted in the numerical model, and the vegetation effects on the flow and morphological changes are considered. The model is verified against three cases from the laboratory and field data documented in the literature. The numerical results are consistent with the measured results, which show that the model could accurately simulate the evolution of the dam-break flows and the morphology evolution of bed within a computational domain with complex plant distribution. The results show that the riparian vegetation in the waterway narrows the channel and reduces the conveyance capacity of river. The flood flow is diverted away from the vegetation community toward two sides and forms a weak flow region behind the vegetation domain. The resistance of plants markedly reduces the flow velocity, which directly alters the fluvial processes and influences the waterway morphology.
1. Introduction
Increase in catastrophic flood events has attracted increasing attention on the prediction of their dynamics and bed change evolution, especially in relation to riverine plant effect. Usually, grasses, shrubs, and mangroves, growing in watercourses and floodplains, are key members of the water ecosystem. They play important roles in water purification, flood control, and maintenance of bank stability, but they also create nonpositive obstruction effect on flood propagation in waterways. Recently, interactions between flow, morphology, and aquatic plant have been studied in the laboratory and field-scale experiments. Manners et al. (2014) believed that rapid expansion of tamarisk led to a narrower channel, which pushed the Yampa River into a new equilibrium having altered fluvial processes [1]. Some researchers employed flume experiments to investigate the potential effects of living vegetation and large wood on river morphology, specifically aiming to explore how different wood input and vegetation scenarios impact channel patterns and dynamics [2]. Asami et al. (2012) investigated the morphological characteristics of refugium where riverine plants survived large floods [3]. Because numerical method has its advantages of wide application range, repeatability, and low cost compared to other methods, some effective and efficient numerical methods were developed for understanding flood wave propagation in urban areas, complex open channels, and overland flow systems [4]. Dam-break flows were usually formed in mixed flow regimes with discontinuities; the numerical schemes that were often used for such studies include the shock-capture schemes, such as Total Variation Diminishing schemes (TVD) and Godunov type schemes [5–8]. During a flooding hazard involving rapid transients, interaction between dam-break waves and topography changes could be significant; however, the interaction was ignored in the early numerical models for simulating dam-break flows over mobile beds.
During dam-break flow conditions, the flow velocity is large, the sediment concentration is high, and the bed varies so rapidly that their effects on the flow cannot be ignored [9]. To correctly predict the consequences of a dam failure in a complex topography, the interaction between flow and bed morphology must be modeled using a coupled solution [9–12]. However, such models are based on either fixed rectangular grids or uniform curvilinear boundary-fitted grids and, therefore, simulate the large-scale flow features on a structured grid system, which may consequently reduce the modeling accuracy and efficiency [13–15]. One of the major objectives of studies was to investigate the turbulence structure in vegetated environments [16, 17]. However, there are restrictions and limitations about the knowledge of the positive or nonpositive effects on the topography change caused by vegetation. The crucial contributions of plants to the flood hydrodynamics and the bed change need to be recognized.
In the present study, a depth-averaged hydrodynamic and sediment transport model, based on finite-volume method, is developed to simulate the flood waves and bed change due to vegetation effect. For improving simulation accuracy, local mesh at selected locations is refined by using triangular mesh. The proposed model is first used to calculate the dam-break flows over fixed bed and to compare the water levels and velocities with measured data. Then, the bed changes of dam-break flood are computed to investigate the temporal and spatial variability in the bed elevation. Finally, the hydrodynamic variation and the evolution of bed elevation through the rigid vegetation domain are investigated and discussed.
2. Mathematical Equations
2.1. Governing Equations
The depth-averaged shallow water equations were obtained by integrating the Navier–Stokes equations over the flow depth, consisting of continuity equation and momentum equations for depth-averaged free surface flows. The conservative and vector forms of the 2D shallow water equations are written in (1) and (2), respectively, as follows:where is the time, and are Cartesian coordinates describing the horizontal plane, , , and are vectors of the conserved flow variables and the convection fluxes in the and directions, respectively, and is source term [18].where is the flow depth; and are the depth-averaged velocity in the and directions, respectively, and and are the bed shear stress term with and components defined by the velocities and , respectively. , , and is Manning’s coefficient. Here, is the water level, is the gravitational acceleration, and is density of the water and sediment mixture in the water column, determined by . Here, is the volumetric concentration of total-load sediment dimensionless, is density of the water and sediment mixture in the bed surface layer, determined by , is the porosity of the bed material, and and are water and sediment densities. and are the sediment entrainment and deposition fluxes, respectively. and are the vegetation drag forces in the and directions.
The sediment deposition and entrainment are estimated as follows:where is an empirical coefficient and is the settling velocity of sediment particle [9] given as follows:where is kinematic viscosity and is the sediment diameter. in (4) is the depth-averaged sediment transport capacity given as follows:where is the bed-load sediment transport capacity given as follows:In (7), is the modified parameter, is the Shields parameter such that , is the threshold Shields parameter, and .
The bed deformation can be determined as follows:where is the bed surface elevation above a reference datum.
2.2. Vegetation Resistance
The vegetation effect on the flow was included to the momentum equations as an internal source of resistant force per unit fluid mass. Therefore, the drag force exerted on vegetation per unit volume can be expressed as follows [16]:where is the drag force coefficient, is the vegetation density defined as number of vegetation elements per square meter, is diameter of a vegetation element, and is the height of vegetation.
Equation (9) is derived from experimental and field data and it has been extensively used in wave propagation and it could be also employed to describe dam-break wave propagation.
3. Numerical Method
3.1. Finite-Volume Method
A standard finite-volume method with Roe Riemann solver for wet and dry computation has been developed. The discretization of the governing equations was based on the finite-volume method, for which the unstructured triangular mesh is shown in Figure 1. The conserved variables were defined at the cell centers and represented the average value over each cell. Hence, the integral form of (1) over the th control volume can be written as follows:where subscript indicates the element number and subscript indicates the side of the element. is the advective fluxes term of the water and is the diffusive fluxes term of the water. is the control volume of the th element. Using Green’s theorem, (10) becomeswhere is the average value of the conserved variables over the th cell and is stored at the center of each cell, with . is the boundary of and is the area of the th cell. is the outward surface normal vector of , with , where is the included angle between the outward normal vector and the direction.
In a 2D triangular grid system, the line integral term can be approximated and assessed as follows:where is the total number of edges of the control cell ( is 3 for a triangular mesh); the subscript denotes the index of the edge of a triangular cell; is the arc length; and presents the outward normal flux vector.
3.2. Evaluation of Numerical Fluxes
Riemann problems at each cell interface were solved by various Riemann approximations for assessing the interface fluxes. In particular, Roe’s approach was used in this paper. The interface flux of Roe’s solver was computed as follows:where and are reconstructed Riemann state variables on the right and left sides, respectively. The flux Jacobian matrix is expressed as follows:where and are the components of the outward surface normal vector in the and directions, respectively, and is Roe’s averaged wave velocity.
The three distinct eigenvalues of (by hyperbolicity of the system) areEquation (14) can be alternatively expressed such that . Here, is a diagonal matrix of the absolute values of the eigenvalues of given as follows:and and are the right and left eigenvector matrices given as follows:The Riemann state variables , , and on the face boundaries were needed to calculate the flux and they are given by Roe’s average method as follows:The subscripts + and − in (18) refer to the right and left sides of the edge considered, respectively. It should be noted that when the dry-wet interfaces exist, Roe’s average was defined as and .
3.3. Evaluation of Sediment Fluxes
The sediment flux across the interface of two neighboring elements in (2) was solved by a flux interpolation technique, instead of Roe’s approach. The sediment flux value across the interface was determined based on the neighboring values at element centers, such thatwhere is the sediment convection flux; and are the distances between the interface center and the elements center sharing the same interface, respectively; and are the outward surface normals in and directions; and is the length of the interface.
3.4. Treatment of Wetting and Drying Fronts
A treatment technique for wet and dry boundaries was adopted to achieve the purpose of zero mass error [19–21]. A criterion, ε, was used to classify the following four types of edges:(1)Wet edge as shown in Figure 2(a): two adjacent cells are wet, in which and .(2)Partially wet edge (with flux) as shown in Figure 2(b): a wet cell (left) links to a dry cell on the right, and the water level of the wet cell is higher than that of the dry cell, in which , and .(3)Partially wet edge (no flux): a wet cell (left) links to a dry cell on the right, and the water level of the wet cell is lower than that of the dry cell as shown in Figure 2(c), in which , , and . To overcome a nonphysical flux problem produced in the numerical model, the water level and bed level for the dry cell were temporarily replaced by a value which equaled to the water level in the wet cell.(4)Dry edge as shown in Figure 2(d): two adjacent cells are dry, in which and .According to these four types of edges, cells were correspondingly divided into three types as follows:(1)Wet cell: all the edges of this cell consisted of wet edge or partially wet edge (with flux) and all the nodes of the cell were flooded.(2)Dry cell: all the edges of this cell consisted of dry edge or partially wet edge (no flux).(3)Partially wet cell: all other cells were not satisfying the criteria of wet and dry cell as defined above.
(a)
(b)
(c)
(d)
4. Numerical Study
4.1. Dam-Break Flow with Natural Geology
The numerical model has been validated against theoretical solutions (not the focus of the paper). In the following, a series of applications are illustrated to appreciate the model applicability to real cases for which experimental data are available in the literature. First study considered the failed Malpasset dam that was once located in a narrow gorge of the Reyran river valley in Southern France with a width 500 m. The dam site was located 12 km upstream of Frejus estuary [11, 21]. The arch dam had a height of 66.5 m and a crest length of 223 m length. The maximum reservoir capacity of the Malpasset dam was 55 × 106 m3. In the immediate downstream of the dam, the Reyran river valley was very narrow and had two consecutive sharp bends. Then the valley widened as it went downstream and eventually reached a flat plain. The dam failed in 1959 following exceptionally heavy rain. After the dam failure, a field survey was performed to obtain the maximum water level along the Reyran river valley. In addition, a physical model with a scale of 1 : 400 was built to study the maximum water level and the flood wave arrival time at nine points along the river valley in 1964. Because of its complex topography and availability of measured data, the Malpasset dam-break case was selected as a test example for the present dam-break model.
In this computation, 11800 grid cells were composed of the refined mesh sizes at the dam site, main river channel, and downstream valley. The initial water level in the reservoir was set at 100 m above sea level. The rest of the computational domain was considered as dry bed. Manning’s coefficient was 0.019 over the entire computational domain. The time interval was 0.01 s. The values of maximum water levels and flood wave arrival times available from physical model were compared with the present numerical results as shown in Figures 3 and 4. The computation results are in good agreement with the measured data scaled up to prototype scale. Figure 5 shows the computed water depths with the irregularities of the computational domain corresponding to 800 s and 1800 s after the dam failure. The wave fronts reached the location of measuring station 12 and near estuary at 800 s and 1800 s after the dam failure, respectively. This simulation demonstrates that predictions for discontinuous flow with wetting and drying over uneven bed are reasonable based on unstructured triangular mesh, and, thus, the model proposed in this study can be utilized for field-scale applications.
4.2. Idealized Dam-Break Flows over Mobile Beds
For the proposed dam-break model, an idealized test on dam-break flow over a mobile bed was numerically investigated to examine the capability of capturing the wet-dry boundary and bed change. The experiment performed by Capart and Young [22] was adopted in this study. The flume was 1.2 m long, 0.2 m wide, and 0.7 m high. It was initially covered by a 5 to 6 cm thick layer of light artificial pearls (not natural sand), its diameter was 6.1 mm, specific gravity was 1.048, and settling velocity was 0.076 m/s. The bed of the experimental channel was horizontal. In the experiment [22], an idealized dam was located in the middle of the flume separating the upstream static flow region of 10 cm depth from the dry downstream part. At s, the dam was removed rapidly to form the dam-break wave at the downstream. Manning’s coefficient was 0.025 and the time step was 0.001 s; the total simulation period was 0.6 s. The computational domain was divided into 10354 triangular grids. The bed porosity was 0.28, the threshold Shields parameter was set at 0.15, and and were equal to 3 and 6 in this study [12]. Figure 6 shows the comparison of calculated and measured water levels and bed surface elevations at different time series. The results demonstrate that the present model produced predictions in good agreement with Carpart’s experimental data [22]. The erosion magnitudes and the wave front locations were well predicted by the present model. A hydraulic jump was formed near the initial dam site and an obvious scour arose near the hydraulic jump. These discrepancies regarding the location of the hydraulic jump are typical of the numerical studies, as suggested by [9]. In general, the hydrodynamic and sediment model was capable of accurately reproducing physical properties of dam-break waves with wet-dry moving boundaries.
4.3. Dam-Break Flows over Mobile Beds with Suddenly Widening Channel
An experiment of dam-break flow over mobile beds was adopted to test the capability of the present model by calculating the interactions on flow, sediment transport, and bed change [23]. The test involved a 6 m long flume with nonsymmetrical sudden enlargement from 0.25 m to 0.5 m width, located 1.0 m downstream of the gate (Figure 7). The initial conditions consisted of a 0.1 m high horizontal layer of fully saturated sand over the whole flume. The sediment used was uniform coarse sand with the median diameter of 1.82 mm, density of 2680 kg/m3, and porosity of 0.47 and the Manning’s coefficient was 0.022. An initial water depth was 0.25 m at the upstream of the gate. As shown in Figure 8, a nonuniform space step was used for refining the sudden-expanded channel with a total of 5778 triangular meshes. There was a free outflow boundary condition at the downstream outlet. The coordinates of measured nodes (P1–P6) and three cross-sections were as shown in Figure 7. Figure 9 compares the water level-time series at six points from the model prediction and the laboratory data. The predicted water level hydrographs agree quite well with the observed data under the condition of mobile bed. Figure 10 compares the predicted and observed bed level profiles at three measured cross-sections at 3 ( = 4.2 m), 5 ( = 4.3 m), and 9 ( = 4.5 m). At 3, the simulated magnitude of scour depth was identical with the measured data, except for the location of scour hole. A similar trend was observed between the experimental and the numerical results of erosion and of deposition. At 5, the lateral bed profile was predicted correctly and the maximum bed scouring depth was underestimated slightly compared to the experimental measurement. The bed variations included general erosion on the straight bank and a deposition on the sudden enlargement side at 9.
In order to investigate the interaction of flood, bed morphology, and floodplain vegetation for dam-break flow during a dam-break event, the sudden enlargement case in the previous section was recalculated by adding vegetation in a circular zone (the center of the circle had = 4.2 m and = 0.175 m, and the diameter of vegetation domain was 0.15 m) as shown in Figure 11. The plant parameters were as follows: the diameter of trees was 15 mm, the vegetation height was 0.5 m, and the vegetation density was 1500 numbers per square meter. The drag force coefficient was 1.5, and the coefficients and were equal to 10 and 0.5 in this study based on [11]. Figure 12 compares the calculated water levels and absolute velocities between vegetated and nonvegetated scenarios at different stations. A significant difference in the form of the water level and flood velocity was observed: the water level corresponding to vegetation effect was markedly higher compared to the nonvegetation effect at the center of vegetation. The absolute velocity of the presence of vegetation decreased quickly when dam-break wave reached P1 station. The presence of vegetation reduced the conveyance capacity of channel and its floodplain with a consequent raise of local flood levels (P2 station and P4 station); the absolute velocities were accelerated at P2 and P4 before 5.4 s and 4.5 s and decelerated thereafter. At P3 station, the differences of the water level variations were not apparent, but absolute velocity varied markedly. As expected, the existence of vegetation caused the redistribution of the velocity field, and the velocities in the vegetated domain were reduced. Figure 13 compares the calculated final bed levels between the vegetated and nonvegetated scenarios for two cross-sectional profiles ( = 4.2 m and = 4.35 m). In the absence of vegetation, the deposition hump moved to the expansion zone because of the nature of sudden dam-break flow. The erosion occurred near the two sides of vegetation community, and the sediment was mostly deposited inside the vegetated domain and behind the vegetation group. Figure 14 shows the calculated bed variation and flow pattern near the sudden expansion zone at different time series. It also shows the difference in currents and topographical features for scenarios with and without vegetated resistance. The time series showed a very high temporal variability in the velocity field in the sides and around the vegetation groups. Flow was diverted away from the vegetation zone toward two sides and subsequently enhanced sediment deposition around vegetation group. Sediment erosion was increased in nonvegetated areas to satisfy sediment conservation.
5. Conclusions
In this study, an accurate and efficient depth-averaged hydro-morphodynamic model was developed based on finite-volume method using the unstructured triangular grid. In solving the associated equations, a framework of the fully coupled procedure was deployed with flow, sediment, and bed change computed simultaneously in the entire computational domain. The intercell fluxes were evaluated based on Roe’s approximation of Riemann solver with second-order accuracy as obtained by employing a MUSCL reconstruction technique, which provides an accurate description of flow near the moving waterline for dry and wet boundaries. The hydrostatic pressure was put into the source term of momentum equations such that it successfully eliminated the numerical imbalance between the source and the flux terms [19]. The vegetation field was represented by a momentum sink term as a dynamic drag force. The model was validated with several types of flood propagation problems, including dam-break flows over fixed and mobile beds under the conditions with and without aquatic plants. The plants formed a group providing a buffer against the disturbance by reducing floodplain stripping for flood and the flow was restricted because of the narrow pass. The vegetation community also altered the flow direction and then changed the deposition hump by causing variation in flooding velocity.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (51579030), the Program for Liaoning Excellent Talents in University (LJQ2013077), the Liaoning Natural Science Foundation (2014020148), and the Open Fund of the State Key Laboratory of Hydraulics and Mountain River Engineering (SKHL1517).