Interactions Study of Hydrodynamic-Morphology-Vegetation for Dam-Break Flows

1 School of Ocean Science and Environment, Dalian Ocean University, Dalian, Liaoning 116023, China 2State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, Chengdu, Sichuan 610065, China 3State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian, Liaoning 116025, China 4Yuanyanggou National Ocean Park, Panjin, Liaoning 124010, China


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][6][7][8].During a flooding hazard involving rapid transients, interaction between dambreak 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.

Mathematical Problems in Engineering
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][10][11][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][14][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 dambreak 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.

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, U, F, and G are vectors of the conserved flow variables and the convection fluxes in the  and  directions, respectively, and S is source term [18].
where ℎ is the flow depth;  and V 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 V, respectively.
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  =  2 * /,   is the threshold Shields parameter, and  = (  /  − 1).
The bed deformation  can be determined as follows: where  is the bed surface elevation above a reference datum.

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,  V is diameter of a vegetation element, and ℎ V 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.

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.E adv  is the advective fluxes term of the water and E dif  is the diffusive fluxes term of the water.  is the control volume of the th element.Using Green's theorem, (10) becomes where U  is the average value of the conserved variables over the th cell and U  is stored at the center of each cell, with is the boundary of   and   is the area of the th cell.n is the outward surface normal vector of   , with n = (  ,   ) = (cos , sin ), 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 n  presents the outward normal flux vector.

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 U  and U  are reconstructed Riemann state variables on the right and left sides, respectively.The flux Jacobian matrix J is expressed as follows: where   and   are the components of the outward surface normal vector in the  and  directions, respectively, and c is Roe's averaged wave velocity.The three distinct eigenvalues of J (by hyperbolicity of the system) are Equation ( 14) can be alternatively expressed such that | J| = |Λ|.Here, |Λ| is a diagonal matrix of the absolute values of the eigenvalues of J given as follows: and  and  are the right and left eigenvector matrices given as follows: ) . ( The Riemann state variables ũ, Ṽ, and c 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 ũ = ( + +  − )/2 and Ṽ = (V + + V − )/2.

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 that where    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.

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][20][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 ℎ  > .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.

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 × 10 6 m 3 .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.

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  = 0 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.

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/m 3 , 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  s 3 ( = 4.2 m),  s 5 ( = 4.3 m), and  s 9 ( = 4.5 m).At  s 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  s 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  s 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.    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   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.

Conclusions
In this study, an accurate and efficient depth-averaged hydromorphodynamic model was developed based on finitevolume 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

Figure 1 :
Figure 1: Control volume as represented by triangular mesh.

Figure 3 :Figure 4 :
Figure 3: Calculated and measured water front arrival times at each measuring point.

Figure 9 :
Figure 9: Comparison between the observed and calculated water levels.

Figure 10 :Figure 11 :
Figure 10: Comparison between the observed and calculated final bed levels for cross-sectional profiles.

Figure 12 :Figure 13 :Figure 14 :
Figure 12: Comparison between the observed and calculated final bed levels for cross-sectional profiles.
), 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 =   (1 −   ) +     .Here,   is the volumetric concentration of total-load sediment dimensionless,  0 is density of the water and sediment mixture in the bed surface layer, determined by  0 =   (1 −   ) +     ,   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.