Overtopping breaching of river levees constructed with cohesive sediments

Experiments were conducted in a bend flume to study the overtopping breaching process and the corresponding overflow rates of river levees constructed with cohesive sediments. The river and land regions were separated by the constructed levee in the bend flume. Results showed that the levee breaching process can be subdivided into a slope erosion stage, a headcut retreat stage and a breach widening stage. Mechanisms such as flow shear erosion, impinging jet erosion, side slope erosion and cantilever collapse were discovered in the breaching process. The erosion characteristics were determined by both flow and soil properties. Finally, a depth-averaged 2-D flow model was used to simulate the levee breaching flow rates, which is well expressed by the broad-crested weir flow formula. The deduced discharge coefficient was smaller than that of common broad-crested rectangular weirs because of the shape and roughness of the breach.


Introduction
River levees, as a kind of embankment structure, are constructed around rivers and parallel to the main flow to constrain flow and protect local residents from flooding disasters (Schmocker and Hager, 2009).Levees constructed with cohesive sediments are the most common type due to their low cost and the convenience with which construction materials can be acquired locally.When lacking a protection layer, this kind of levee can easily be breached by overtopping flow if the water level exceeds the design water level, including the freeboard (Pickert et al., 2011).Consequently, the protected area will be submerged, threatening lives and properties.Under the circumstances, predicting flood propagation processes and repairing the levee breach as soon as possible is crucial for diminishing losses.For this, a profound understanding of the cohesive levee breaching process and overflow rates is necessary.
Much research has been performed on overtopping breaching of the embankments constructed at a normal angle to river flow, and quite a few studies are about embankments constructed with non-cohesive materials.For example, Coleman et al. (2002) found that the breach channel of overtopped embankments under constant water level has a curved shape, and the breach development obeyed the minimum energy dissipation rate rule for streams.Schmocker and Hager (2009) studied the required minimum sediment size, dike width, dike height and unit discharge in their laboratory experiments of non-cohesive dike breach to avoid the side wall effect, scale effects and cohesion.Pontillo et al. (2010) applied a 1-D two-phase model to simulate flow propagation and the breaching process of a trapezoidal-shaped sediment dike.Pickert et al. (2011) found that embankments composed of finer materials exhibited discontinuous erosion affected by cohesion due to porewater pressure.Based on a large amount of experiment data, Schmocker and Hager (2012) studied the effect of sediment size, dike height and inflow discharge on plane dike-breach due to overtopping, and proposed three formulae that express the relationship between the final dike height, final dike volume, maximum breach discharge and the above three factors.Case Q (L s −1 ) w ( %) e (%) ρ (kg m −3 ) ρ d (kg m −3 ) c (KPa) ϕ (  Many studies on the cohesive imbankments have also been conducted.Compared to non-cohesive embankments, overtopping breaching of cohesive embankments is a more complex phenomenon involving impinging jet flows, reverse roller structures and headcut erosion (Zhu, 2006).Many researchers (Ralston, 1987;Powledge et al., 1989;Hanson et al., 1999;Wahl, 2004) have found that headcut retreat is a predominant mode during the overtopping breaching process of cohesive embankments, and many prediction models of headcut retreat rate have been put forward (Hanson et al., 2001;Stein and LaTray, 2002;Zhao et al., 2013).Zhu (2006) and Zhao (2016) respectively developed a model for overtopping breaching of cohesive embankments in which the headcut erosion process was included.Hanson (1999), the IM-PACT project (Morris et al., 2005), Zhang et al. (2009) and Zhu (2011) all studied the cohesive embankment breach process based on laboratory experiments or large-scale experiments.The IMPACT project (Morris et al., 2005) was also concerned with the influence of soil grading, soil water content, compaction and embankment geometry, while Zhang et al. (2009) studied the cause of headcut erosion, double spiral flow at the dam crest and the effect of soil cohesion on the breach process.
However, research on the overtopping breaching of river levees parallel to the main flow is limited.Kakinuma and Shimizu (2014) conducted several groups of large-scale experiments on breaching of non-cohesive river levees and subdivided the breach development into four stages.Yu et al. (2013) studied the influencing factors on the overtopping breaching process of non-cohesive levees and the breach dis-charge properties.Liang et al. (2002) built a numerical model to simulate the breach of the Yellow River dike built with a combination of non-cohesive and cohesive soils.
River levee breach is quite different from that of embankments constructed normal to the flow in morphology, hydraulics and inflow variation characteristics (Kakinuma and Shimizu, 2014).Moreover, measured data of cohesive levee breach have not been reported until now.Hence, four groups of experiments on cohesive levee breach were performed in a bend flume with varied inflow discharge, soil water content and porosity.In these experiments, levees were constructed in the flume with an initial breach.Different stages of the levee breach process and flow characteristics near the breach were analysed.The process of the breaching overflow rates, simulated by a depth averaged 2-D flow model, was also studied in detail.

Experiment set-up
In a bend, the water surface at the concave bank is higher than that at the convex bank; so levees at the concave bank are more likely to be breached by overflow than that at the convex bank or in straight channels.Thus, we conducted the experiments in a bend flume, and Fig. 1a shows the experiment layout.The cross section of the levee is shown in Fig. 1d.The initially dry area, which was separated from the wet part and will be submerged after the levee breaching, is named as the land region, while the remaining wet reach is the river region.The initial breach was located at the place where overtopping breaching is most likely to appear according to the character- istics of bend flow, and its scale is shown in Fig. 1c.The transverse water-surface gradient of the bend flow produces the bend circulation, with flow directing toward the concave bank in the water surface.This provides the initial velocity of levee breach flow.River flow can flow into the land region easily.At the end of the flume, there were two sluice gates, and the one at the river region was used to adjust the water level.
Four automatic water-level gauges were placed at the points numbered from S1 to S4 shown in Fig. 1a.S2 and S4 were in the land region, while the other two were in the river region.A topography meter was placed above the initial breach, measuring the breach profile variation in time by moving back and forth.Velocity changes near the breach were measured by an acoustic Doppler velocimeter (ADV).The location of the velocity monitoring points (M 1 and M 2 ) is shown in Fig. 1b and d.The breach development was recorded by a camera fixed at the flume wall near the breach.
Experiment cases are listed in Table 1.Such main factors as upstream inflow discharge Q, soil porosity e and water content w of levee materials were considered.The levee was constructed using lean clay (according to USCS), with its grain size distribution shown in Fig. 2. As can be seen, the soil consists of 10 % clay particles (d < 0.005 mm), 70 % silt particles (0.005 mm < d < 0.075 mm) and 10 % sand particles (d > 0.075 mm).The liquid and plastic limit of the soil is 28.6 and 18.1 % respectively.Soil density ρ, dry density ρ d , soil cohesion c and initial friction angle ϕ are all affected by e and w.
The levee was constructed by gradually adding soil and compacting it, layer by layer.Before experiments began, the levee surface in the river region was covered with thin films to prevent infiltration.Soil samples were selected from the levee to test w, e, c and ϕ, the values of which are listed in Table 1.Flow entering the flume was adjusted until the discharge reached the scheduled value in Table 1.Then the sluice gate at the river region was adjusted to ensure a very slow rise of the water level.The thin films were removed just before the flow overflowed the levee top.

General description
The process of overtopping breaching can be subdivided into three stages according to the breach erosion characteristics shown in Fig. 3.The initial stage, characterised by flow shear erosion on the levee slope at the land side, is named as the slope erosion stage (shown in Fig. 3a).In this stage, some small scour holes appeared on the land-side slope.Enlargement of the scour holes steepened the slope and then a largescale scarp known as headcut developed (shown in Fig. 3b).This is the beginning of the second stage defined as the headcut retreat stage, shown in Fig. 3b, c, d and e.At the end of this stage, the breach cross section was almost washed out and a deep gully formed (shown in Fig. 3f).Then flow in the gully began to erode side slopes of the breach, which is the third-stage erosion defined as the breach widening stage.

Slope erosion stage
Initial erosion of this stage usually occurred near the levee toe at the land side, due to large flow shear stress there (shown in Fig. 4), which shows the water surface and flow shear stress distribution in a certain moment after overtopping of a levee calculated by Briaud et al. (2008).The large negative shear stress at the toe is caused by flow separation due to sharp corners between the levee and bed surface (Briaud et al., 2008).If the shear stress there surpassed the soil critical shear stress, cohesive soil blocks would be eroded and a scour hole (shown in Fig. 3a) would appear.The small scour hole increased bed roughness and flow turbulence, which in turn accelerated local scour, making the scour hole enlarge.Scour holes may have also appeared at other weak places but they developed more slowly and finally merged into the large scour hole at the bottom.

Headcut retreat stage
According to its shape, the headcut can be categorised into single-step form (Fig. 5a and d) and multiple-step form (Fig. 5b and c).Initially, an incomplete single-step headcut appeared due to the erosion at the bottom (shown in Figs.3b  and 5a).After the headcut retreated to the brink of the levee crest, a multi-step headcut appeared due to the layer structure of the levee (shown in Fig. 3c, d, 5b and c) and at last, all the steps disappeared and merged into a complete singlestep headcut (shown in Figs.3e and 5d).
When the overflow velocity was small, the flow just streamed down along the vertical headcut surface where flow shear erosion mainly occurred.If the flow velocity increased, the overflow departed from the headcut surface, forming an   Briaud et al. (2008) impinging jet.For multi-step headcuts, when the flow velocity was not large enough, multiple small impinging jets may appear due to the steps (shown in Figs.3d and 5c).For singlestep headcuts or multiple-step headcuts with large-velocity overflow, a single impinging jet existed (shown in Figs.3b,  e, 5a, b and d).The single jet and small bottom jet directly impinged the non-erodible foundation, with part reflected toward the headcut as a reverse roller which undermined the headcut face.The other small upper jets, however, impinged the erodible platforms of the headcut steps, exerting both normal and horizontal shear stress on the platforms.The platforms were eroded or collapsed and finally disappeared.The above erosion mode, caused by impinging jets, was defined as jet impinging erosion.
It should be noted that in this stage, flow shear erosion and jet impinging erosion usually appeared simultaneously, with the common style that flow shear erosion happened above and the other happened below; or they distributed alternatively along the river flow direction.Except for the two erosion modes, discrete mass failure may also appear during the erosion process.
The breach height reduction, caused by flow shear erosion, was with a very small rate.This made the overflow velocity quite small (about 0.2 m s −1 ), hardly influenced by inflow discharge.So when soil properties were similar and inflow discharge the same (cases 3 and 4, for example), headcut retreat rate differed slightly (shown in Fig. 6).Furthermore, there is a negative correlation between headcut retreat rate and the soil cohesion (shown in Fig. 6 and Table 1).A reasonable explanation is that for cohesive soil, large cohesion between particles can enhance soil erosion resistance decisively.

Breach widening stage
During this stage, the breach side slopes under water were eroded by flow shear stress and retreated, causing the remaining part above the water to be suspended as a cantilever, shown in Fig. 3f.
The breach widening process under water is shown in Fig. 7. Initially, the slope migration rate of both sides was similar.However, later, the downstream side eroded faster because breaching main flow concentrates near the downstream side slope, as shown in Fig. 3f.In addition, the flow velocity near the downstream side is larger than that near the upstream side.(This can be seen in Fig. 8 showing the calculated flow field near the breach of case 2 when the breach width was about 40 cm using the following 2-D numerical model.)Meanwhile, around the downstream side, the flow velocity near the toe of the land region is larger than that of    pared to case 3. The same as in the headcut retreat stage, decrease of soil cohesion can accelerate the breach widening.
The cantilever above the water surface can sustain itself for a while before collapse, attributed to the inner tensile stress of soil.Assuming that the cantilever is cuboid and sustains flexural deformation, the forces acting on the cantilever in a critical fracture state are shown in Fig. 10.The inner tensile stress on top of the fracture surface approaches the soil tensile strength σ t , and the moment equilibrium equation of the cantilever per unit width is expressed as in Fukuoka (1994): where G is the cantilever weight per unit width expressed as G = γ L c H c ; γ is the unit weight of the cantilever expressed as γ = ρg.L c and H c are the length and height of the cantilever when fracturing, respectively.The critical cantilever length can be deduced from Eq. ( 1) as (2) Assuming that the soil tensile strength is 0.7 times of the soil cohesion (Zhu, 2008), the average critical cantilever length L c can be calculated by Eq. ( 2), and the results are shown in Table 2.It can be seen that the calculated L c rather approximates to the measured value, which proves the reasonability of the above assumption about the cantilever fracture.
4 Flow characteristics near the breach

Flow velocity and water level variation
The velocity monitoring point of case 1 is located at M 1 while that of case 2 is located at M 2 , as shown in Fig. 11.x is horizontal and perpendicular to the levee axis pointing to the land region, direction y is horizontal and parallel to the levee axis pointing to the downstream, while direction z is vertically upward.U x , U y and U z denote velocity in the x , y and z direction, respectively.U is the resultant velocity.
The variation of water level in the land region corresponded well with the breach height H , increasing gradually initially but sharply along with a sharp drop of the breach height.Then it kept almost constant, although the breach went on widening.
The flow velocity variation was also related to the breach height changes.For case 1 (shown in Fig. 12a), before 75 min, U increased slowly with gradual decrease of the breach height.From the time of 75 min, the breach height H dropped sharply, and correspondingly, U raised sharply with 1 min lag.At the time of 78 min, H reduced to the minimum value of zero and U reached the maximum value of 0.7 m s −1 .After that, along with the increasing water level in the land region, U dropped gradually and stabilised itself at the value of about 0.5 m s −1 .
General flow direction and its variation trend can be deduced from the component flow velocities.Initially, U z was almost zero, demonstrating the 2-D characteristics of the breach flow.U x was almost the same as U y initially but surpassed it afterwards, indicating that the flow began to be toward the breach due to decrease of the breach height.With the sharp decrease of the breach height, there was a different variation trend of U y .For case 1 (shown in Fig. 12a), M 1 is located at the upstream side of the breach and it is inevitable that there is component velocity in the downstream direction, so U y increased and stayed at a relatively larger value; while for case 2 (shown in Fig. 12b), M 2 is located in the middle of the breach and the flow velocity there was directly toward the land region, so U y decreased to almost zero.

Overflow rates
It is rather difficult to acquire the flow rates over the levee breach directly by measurements; so a depth-averaged 2-D flow model was established to compute them.In the calculations, the breach geometry was modified at each time step according to the above measured breach developing process.

Numerical methods
The governing equations, written as vector form, are given as follows: where U, F, G and S are expressed as where x and y are Cartesian coordinates shown in Fig. 13a; h is water depth, u and v are flow velocity in the x and y direction, respectively.S bx , S by , S f x and S fy are respectively the bed slope and friction slope in the x and y directions.S f x and S fy are expressed by the Manning equation as follows: where n is the Manning roughness coefficient.Equation ( 3) was discretised by the finite volume method in space.The calculating grids are shown in Fig. 13b.At any calculating cell (i, j ), the equation was discretised as where x and y are grid size in the x and y direction respectively.The numerical fluxes F and G were calculated by the by weighted essentially non-oscillatory (WENO) method together with Roe method based on Riemman solvers that can capture flow discontinuity (Dou et al., 2014).The thirdorder Runge-Kutta method was used to discretise Eq. ( 5).
As the scheme is explicit, time steps should satisfy Courant-Friedrichs-Lewy (CFL) criteria (Toro, 1999): where N CFL is Courant number; t is computational time step; N CFL was valued much smaller than unity because of large bed elevation variation.Here, it was proved that the value of 0.03 for N CFL can achieve a reasonable result.Rectangular grids with size set as 0.5 cm × 0.5 cm were deployed.To avoid a zigzag boundary caused by rectangular grids, the calculating area was extended to a large rectangle ABCD shown in Fig. 13a.The elevations of the area outside the flume were all assigned a larger value than the maximum water level in the flume.A simplified balancing-point method (Zhou, 1988) was used to simulate solid boundary conditions.This is on the assumption that the boundary node is in the middle of the two nodes and the boundary is vertical to the x direction at the point.When grid size is sufficiently small, the method is viable, without changing the real boundary obviously.During the propagation process of levee breach flow, temporarily dry nodes without water may exist.To handle the problem, a minimum water depth of 0.001 m was defined and initially all dry nodes were assumed to be 0.001 m.At a given time, if water depth at a node was less than 0.001 m, then the node is regarded dry and the velocity was set as zero.
The flow rates at the inflow boundary CE were set as that listed in Table 1, and the measured water level process at the point S4 was set as outflow boundary condition.Before overflow calculation, with the given inflow rate and corresponding water level at the point S3, a steady-state flow condition in the river region was calculated, which was set as the initial condition of overflow.The Manning roughness coefficient (0.02) was acquired by comparing the calculated and measured water level.
The model was used by Dou et al. (2014) in simulating the non-cohesive levee breaching process in the same flume.The calculated water level variation process matched well with the measured data, which verified the reliability of the model.
The flow parameters for calculating the overflow rates were defined at the initial land-side brink of the levee top.According to the computed velocity and water depth at each grid between the breach, the average velocity and water depth at the brink can be acquired, by which the overflow rates were calculated.

Calculated results and analysis
The calculated overflow rates Q b are shown in Fig. 14.Initially, the overflow rate was very small and increased gradually due to the large and slowly reducing breach height.Then it increased sharply, following a rapid decrease of breach height.Just after the breach height decreased to zero, the overflow rate reached its maximum value and then decreased.Finally, the overflow rate kept almost stable at a certain value.It can be also seen that the maximum overflow rate for case 4, with a larger inflow discharge (28.53 L s −1 ) than the other three (14.64L s −1 ), is also larger.
To verify the modelled results, we make use of the measured data of case 2, of which the velocity monitoring point    is located at the middle of the breach and the flow direction there was directly toward the land region.Then the measured velocity here can represent the average breaching flow velocity.In Table 3, we list the breach geometry parameters (the breach width B and the breach height H ), the hydraulic parameters (the water head above the breach crest h and the measured velocity U ), the estimated breach discharge Q bm (calculated by Q bm = Bh U ) and the numerically calculated breach discharge Q b .It can be seen that the numerically cal-culated breach discharge matches well with the estimated values overall.The breach overflow rates can also be simulated by broadcrested weir flow formula, which can be expressed as in Hager and Schwalt (1994): where B is the levee breach width and h 0 is the approaching energy head.h 0 = h + v 2 a / 2g, with v a the approaching velocity and h the water head above the breach crest.C d is the dimensionless discharge coefficient.When the effect of approaching velocity can be omitted, Eq. ( 7) can be expressed as For simplicity, Q b , B and h are non-dimensionalised by the gravity acceleration g and the initial breach height H 0 to Q b * , B * and h * respectively (Coleman et al., 2004). . ; B * = B / H ; h * = h / H . Then Eq. ( 8) can be replaced by * can be expressed as From Eq. ( 9) and Eq. ( 10), C d is deduced as 0.24, smaller than the discharge coefficient of common rectangular weirs (Hager and Schwalt, 1994) and smooth breach (Zhao et al., 2015).This may result from both the initial trapezoid shape of the breach cross section and the large surface roughness which appeared later on the land side (Pařílková et al., 2012).

Discussion
The main factor resisting flow erosion of non-cohesive sediments is the effective gravity of sediment particles, while for compacted cohesive soils, the main factor is the cohesive force between soil particles, which is usually much larger than the effective gravity of non-cohesive sediment particles.This makes the affecting factors for overtopping breaching of non-cohesive and cohesive levees quite different.Noncohesive levee breaching is mainly affected by particle diameter and density, while for cohesive levee breach, the main factors are those affecting the cohesive force such as soil composition, compaction degree and soil water content.The breaching processes of the two kinds of levees are also quite different.Compared with non-cohesive levees, the overtopping breaching of cohesive levees is much slower and it consists of two special stages -the slope erosion stage and the headcut retreat stage.This can be attributed to the strong anti-erodibility of cohesive soils.The cohesive levee crest and land-side levee slope are eroded with slower rate and hardly form erosion gully; thus the slope erosion stage exists for a period.Additionally, because the erosion rate of the levee crest is smaller and the land-side levee slope does not collapse easily, the morphology of headcut forms.For non-cohesive levees, the fast-eroded sediment materials from the levee crest deposit at the land-side levee slope toe, and frequent collapse in the land-side slope occurs.This makes the breach a gentler slope that seldom forms abrupt elevation changes like headcut.

Conclusions
Overtopping breaching of levees constructed with cohesive sediments can be subdivided into three stages according to breach erosion characteristics.The initial stage was defined as the slope erosion stage, characterised by flow shear erosion on the land-side slope and small scour holes at the bank toe.The next stage was headcut retreat stage, in which single-step and multi-step headcut can be discovered.Both flow shear erosion and jet impinging erosion existed in this stage.Headcut migration rate was hardly influenced by inflow discharge but will decrease obviously with the increase of soil cohesion.The final stage was the breach widening stage, including erosion of side slopes under water and the cantilever collapse above water surface.The migration rate of side slopes was affected by both soil cohesion and nearby flow characteristics.The maximum cantilever length, how-ever, was influenced by soil density, average cantilever height and soil tensile strength.Both the water level and flow velocity variation near the breach corresponded well with the breach height changes.The magnitude of the component flow velocities near the breach can indicate general direction of the breaching flow.
The calculated overflow rates varied the same as the flow velocity near the breach and increased along with the inflow discharge.By substituting the calculated overflow rates and measured breach size into the broad-crested weir flow formula, the discharge coefficient was deduced, which was smaller than that of common rectangular broad-crested weirs.

Figure 1 .
Figure 1.Levee model: (a) top view of experiment layout, (b) velocity monitoring point, (c) longitudinal section of initial breach and (d) transverse section of initial breach.

Figure 2 .
Figure 2. Grain size distribution of experiment materials.

Figure 4 .
Figure 4.A typical flow shear stress distribution along the dike calculated byBriaud et al. (2008)

Figure 10 .
Figure 10.Critical fracture state of a cantilever.

Fig. 12
Fig.12shows variation with time of the water level Z at the point S2 and the water surface velocity at monitoring point M 1 (M 2 ).Directions of the velocity are shown in Fig.11.Direction x is horizontal and perpendicular to the levee axis pointing to the land region, direction y is horizontal and parallel to the levee axis pointing to the downstream, while direction z is vertically upward.U x , U y and U z denote velocity in the x , y and z direction, respectively.U is the resultant velocity.The variation of water level in the land region corresponded well with the breach height H , increasing gradually initially but sharply along with a sharp drop of the breach height.Then it kept almost constant, although the breach went on widening.The flow velocity variation was also related to the breach height changes.For case 1 (shown in Fig.12a), before 75 min, U increased slowly with gradual decrease of the breach height.From the time of 75 min, the breach height H dropped sharply, and correspondingly, U raised sharply with 1 min lag.At the time of 78 min, H reduced to the minimum value of zero and U reached the maximum value of 0.7 m s −1 .After that, along with the increasing water level in
10) Depending on the simulated overflow rates and the measured breach hydraulic parameters, the relation between Q b * and B * h 1.5 * is shown in Fig. 15.With a correlation coefficient of 0.93, the fitting relation of Q b * and B * h 1.5

Table 1 .
Experiment cases and parameters.

Table 3 .
Numerically calculated breach discharge and estimated values t (min)