An analytical process-based approach to predicting breach width in levees constructed from dilatant soils

An accurate prediction of the breach widening rate after the onset of a levee failure is essential for flood risk assessments. Current state-of-the-art analytical breach growth relations are empirical in nature. The large variety in loading conditions, levee design, and levee construction material, combined with the limited amount of accurate measurements of breach growth, make the development of accurate empirical breach growth relations challenging. In this paper, a process-based breach widening relation is presented for levees constructed of dilatant soils. The process-based relation is derived from the weir flow equation and a process-based erosion equation. The breach widening relation can account for the effects of variations in soil parameters. For those cases for which soil parameters are unknown, a calibrated catch-all-coefficient is provided. The relation is benchmarked against the state-of-the-art empirical breach growth relation used in the Netherlands and validated against data on historical levee failures and experimental data.


Introduction
The EU Floods Directive 2007/60/ EC (European Parliament 2007) was developed with the aim to reduce and manage flood risk with respect to public health, environment, cultural heritage, and economic activity. Flood risk maps were completed for these areas by 2013, and prevention, protection, and preparedness plans were developed by 2015 (European Parliament 2007). Nevertheless, research is ongoing on how to improve the accuracy of flood risk assessments. To accurately model flood spreading, an accurate prediction of the flow rate of water entering a flood area is essential. Any error in predicting the rate of breach formation directly affects the accuracy of breach discharge predictions. The discharge and water depth are both important parameters for determining the extent of damages and threat to life (Jonkman 2007). An error in discharge predictions therefore affects consequence estimates and flood risk assessments. When the flood area serves as a retention area, the volume of water leaving the main water system through a breach also reduces the 1 3 risk of flooding downstream. An improved method for predicting the flood volume through a breach then becomes essential for a more reliable assessment of different flood scenarios.
The breach initiation stage depends on the mechanism of failure. For overtopping failures, this stage is characterized by failure of the levee protection, and subsequent erosion of the embankment soil. For piping, this may entail the formation of a full pipe, growth of the pipe, and subsequent collapse of the pipe roof (Mohamed 2002). During breach growth, the full cross section of the embankment is exposed to erosion and the breach widens until it reaches its final width. This breach growth stage is generic. Since peak discharges often occur during the breach growth stage, this stage is especially important for flood risk analyses. Simple data-based empirical breach relations have been developed in support of flood risk analyses (Verheij and Van der Knaap 2003;Wahl 2004). These relations are not based on a process-based simulation of the breach development but are the outcome of a parametrized analysis of the final breach size and the breach development time (Wahl 2004). The large variety in loading conditions, construction materials, and levee design, combined with limited data available on failed levees, limit the accuracy of empirical breach growth relations (Morris 2011). Semi-physically based, analytical or parametric models improve the predictive capability of models by including a process-based description of the breach without complicating the procedure (Morris et al. 2009). To expand the range of applicability, process-based breach models have been developed which model the breach flow and erosion processes during breach initiation and breach widening. Although the flow rate predictions in semi-physically based models or process-based breach models are based on the process-based description of the flow, the accuracy of these models is often still governed by the accuracy of empirical relations underlying the erosion rate predictions.
Here a process-based, analytical breach width equation is presented. The relations are based on the process-based description of erosion of dilatant soil and a process-based description of the flow. The process-based erosion relations have already been validated against erosion experiments performed on sand (Van Damme 2018; Bisschop et al. 2016). Validation of the process-based erosion relation for other soil types is recommended before applying it to other soil types. The predictions following from the new semi-analytical model solution have been compared against data and benchmarked against the parametrized model of Verheij and Van der Knaap (2003) for sand dikes, which partly forms the basis in flood risk analyses in the Netherlands. Model parameters have been calibrated to extend the applicability of the model to cases for which little information is available on the specific type of sandy material.

Methodology
This section outlines how the developed process-based analytical model accounts for the erosion process and breach flow, and how the developed model is validated.

Determining breach growth rates
To account for the continuous change in breach shape due to erosion, the shape of the breach is simplified. A breach widens due to a combination of erosion and mass instability. Breach side slopes are undercut by the flow to the point at which overhanging soil sections become unstable and collapse into the breach (Mohamed 2002). Morris (2011) noted from evaluating the impact of block failures on the breach flow that the impact of block failures 1 3 on the flow rate through the breach is small as the material is rapidly removed from the breach. This was supported by the analysis of video material of large-scale breach experiments performed as a part of the IMPACT project (Morris 2005(Morris , 2011Morris and Hassan 2005). Therefore, a rectangular breach shape which deepens and/or widens (Temple et al. 2005;van Damme et al. 2012) is considered a practical yet reasonable assumption to describe the breach formation stage. Here a rectangular breach shape is assumed which has developed over the full height and therefore solely widens.
The widening rate of a breach is often determined from erosion relations. Many erosion relations relate the erosion rate to the excess shear stress, defined as the operating shear stress minus some critical shear stress, via an erodibility coefficient. How this erodibility coefficient depends on soil characteristics has been difficult to quantify. Hanson and Cook (2004), Hanson and Hunt (2007) and Hunt et al. (2005) noted that material texture, compaction moisture content, and compaction energy all significantly affect the erodibility of soil and thereby significantly influence levee performance (Morris et al. 2009). Of these, material texture and compaction energy have been difficult to account for. Van Rhee (2010) and Bisschop et al. (2016) noted that erosion of dilatant material due to high velocities can be described as the shear failure of thin layers of soil. A soil fails due to shear over a depth d over which the shear stresses exceed the shear resistance ( Fig. 1). Van Rhee (2010) and Bisschop et al. (2016) highlighted the significant effects of dilation on the shear resistance when the ratio between the erosion velocity and the hydraulic conductivity exceeds a factor of 3. During dilation, shear stresses induce an increase in pore volume. The porosity thereby increases from the initial porosity to the critical porosity, defined by that porosity for which shear is possible without a change in pore volume. By solving the mass and momentum balance equations that describe the process of dilation, Van Damme (2018) showed that the rate of displacement of a soil surface due to erosion under high flow velocities corresponds with a maximum in shear resistance induced by dilation. This theory, which inherently accounts for the influence of material texture, compaction energy, and compaction moisture content, underlies the developed breach growth relation. Below this theory is outlined briefly. "Appendix" provides more detail on the full derivation of the process-based method.
For a rectangular-shaped breach profile, the breach sides consist of vertical faces. At a vertical breach face, shear stresses are induced by gravity and the flow of water. The shear stresses induce dilatancy of the soil, which leads to an increase in pore volume. This increase is accommodated by an inflow of water through the surface of the vertical breach face. The pressure gradient required for the inflow results in a reduction in pore water pressures and an increase in shear resistance. During the exchange of water and particles at the soil surface, mass and momentum are exchanged and the breach face displaces. The mass balance equation that describes the displacement rate c (m/s) of a breach face in the normal direction is given by where n 0 is the porosity of the embankment soil prior to dilation, n loose denotes the critical porosity of the embankment soil, v w is the normal velocity component of the water into the wall, and v p is the normal velocity component of particles into the wall. With mass, also momentum is exchanged. Equation 2 describes the momentum balance in the direction of the breach flow parallel to the wall.
Here u p (m/s) and u w (m/s) equal the shear velocity of the breach flow at the vertical breach wall, p (kg/m 3 ) and w (kg/m 3 ), respectively, denote the density of the particles and pore water, g (m/s 2 ) denotes the gravitational constant, x (N/m 2 ) denotes the shear stress component in x-direction, acting on the vertical plane, and K s (m/s) denotes the hydraulic conductivity of the embankment material. The normal velocity components of the water v w (m/s) and the particles v p (m/s) are related to the displacement rate c via The displacement rate c consists of the displacement rate c s induced by erosion by surface shear stresses and the displacement rate c g due to gravity-induced internal shear stresses. A body of soil with a vertical slope is inherently unstable. For a shear failure to occur, soil needs to dilate. The dilation requires water to flow towards the shear plane. The reduction in pore pressures required to induce the water flow temporarily stabilizes the soil (Van der Schrieck 2006; Van Rhee 2010). Due to this gravity-induced dilation, the vertical breach face displaces horizontally with a rate c g . At this displacement rate, the Coulomb friction matches the gravity-induced shear stresses (Van Rhee 2010). The Coulomb friction follows from where c denotes the shear strength, C denotes the cohesion, ′ denotes the effective stress, and is the angle of internal friction. Setting the Coulomb friction equal to the gravityinduced shear stresses now gives Shear stresses applied to the breach wall are transmitted to the deeper soil layers. Integrating Eq. 2 over the depth d, over which the shear stresses x exceed the shear resistance c , gives after substituting which relates the failure depth d to the rate of displacement of the boundary due to erosion c.
The absolute value signs in Eq. 6 make it challenging to develop a closed-form relation for the erosion rates c s and c g as a function of the wall shear stress . For that reason, an empirical approximation has been developed on the basis of the general behaviour Equation 6 shows that the depth d varies with c. For only one displacement rate c, the depth-averaged effective stress ′ , and thus shear resistance, is maximum (Eq. 7 and Fig. 2). This displacement rate corresponds with the most stable situation of the soil and the and the most efficient transfer of momentum.
Predictions for the displacement rate c, for which the shear depth average effective stress is maximum, match the measurements performed by Bisschop et al. (2016) (Van Damme 2018. The effects of dilation on the displacement rate are thereby in line with the observations of Hanson and Cook (2004), Hanson and Hunt (2007), Hunt et al. (2005) that material texture, moisture content, and compaction energy influence soil erodibility.
Here the accuracy of the breach growth rate prediction has been illustrated against a breach widening experiment performed by Hunt et al. (2005). Hunt et al. (2005) performed a large-scale breach widening test on a homogeneous levee constructed of silty sand. The levee was 1.3 m high and had a crest width of 1.8 m. The landside and waterside slope gradients were 1:3 (V:H). Prior to testing, a 0.3-m-wide notch was cut over the full height of the test levee and backfilled with sand. The water level in the test basin was increased until the sand in the breach notch overtopped. With the subsequent rapid removal of the sand, the breach started to widen. During the experiment, the water level upstream of the levee was kept at constant elevation such that the flow velocities 1 3 through the breach remained constant. The average wall shear stress exerted on to the breach side slopes was 20 N/m 2 (Hunt et al. 2005). These levees contained a significant degree of fines. The d 10 was estimated from the distribution of fines at 2 μ m. The density of the particles was 2670 kg/m 3 . The initial porosity was obtained from the dry specific weight and the grain density and was n 0 = 0.35 . The critical porosity was estimated at n loose = 0.46 . The Kozeny-Carman equation (Eq. 9) now gives a hydraulic conductivity of K s = 3.31 × 10 −8 m/s. For = 0 and = 20 N/m 2 , the lateral erosion rate c for which the average shear resistance is maximum is c = 1.02 × 10 −4 m/s (Fig. 3, which corresponds with a breach widening rate of 2.04 × 10 −4 m/s). This prediction was thereby in exact agreement with the measured breach widening rate. Due to the absolute value signs present in Eq. 6, it is difficult to develop a closedform relation for the breach widening rate. For that reason, a closed-form empirical fit to the process-based relation was created. Figure 4 shows the relation between the bed shear stress 0 (N/m 2 ) and the displacement rate c (m/s) alongside the empirical fit given by where m = 0.0003253 m 2 s∕kg and c 1 = 0.00625 m∕s are, respectively, the displacement factor and displacement coefficients for this specific fit. As indicated by Eq. 8, the relation between the displacement rate and the shear stress approximates a square root relation. This corresponds with findings of Bisschop et al. (2016). Considering the accuracy of breach growth measurements, a sufficiently close match is found.
The square root relation between the shear stress and displacement rate, given by Eq. 8, has been applied in the analytical breach growth model. Besides the shear stress, Eqs. 6 to 7 indicate that the displacement rate is a function of the hydraulic conductivity K s , initial porosity n 0 , and critical porosity n loose . The impact of these parameters on the magnitude of the factor m has therefore been determined for the empirical fit. The hydraulic conductivity has been related to the initial porosity n 0 via the Kozeny-Carman equation given by where (m 2 /s 2 ) denotes the kinematic viscosity and d 10 (m) denotes the particle diameter exceeded by 90% of the particles. The particle diameters, initial porosity, and critical porosity are independent parameters. In order to include the effects of the hydraulic conductivity in the displacement relation, it was determined how c 0.5 0 differs with n 0 , n loose , and  Figure 5 shows the initial porosity n 0 set out against c c fit . The dependence of the displacement rate on n 0 is given by The effect of n 0 thereby goes beyond the influence on the hydraulic conductivity. Also the effects of n 0 on the exchange of momentum and soil weight have been accounted for. The procedure was repeated to determine the effects of n loose . This effect is depicted in Fig. 6. This relation is given by Finally, the influence of the particle diameter d 10 on the final solution was determined. This effect is depicted in Figure 7 and is given by

3
A fit to the process-based displacement relation now follows from Eq. 8 whereby m and c 1 are given by the product of Eqs. 10, 11 and 12 multiplied by, respectively, 0.0003253 and 0.00625. The correspondence between the process-based relation and the fitted relation has been evaluated by drawing values for the parameters n 0 , n loose , d 10 , and from uniform distributions. The ranges for these parameters are given in Table 1. Figure 8 depicts the correspondence between the fitted relation and the process-based relation. The coefficient of determination corresponding to the data depicted in Fig. 8 is 0.97. This match was considered to be sufficiently accurate for inclusion into an analytical breach growth relation.

Breach model
This section deals with the description of the flow through the breach, the means of estimating the shear stresses from the flow, and the means of using these stresses to determine the rate of breach growth. Mohamed (2002) noted that the weir equation outperforms the shallow water equations in terms of accuracy of the discharge prediction because it is possible to take the non-hydrostatic effects due to flow curvature into account by calibrating the weir coefficient. Morris (2011) found that use of a variable, instead of a constant weir coefficient, leads to slight improvements in accuracy, but concluded that the primary source of error lies in representing the erosion process. Here the weir equation has been (12) c c fit f (n 0 )f (n loose ) = f (d 10 )502.1(e −414.6d 10 − e −428.7d 10 ) 1 3 combined with the empirical fit to the process-based displacement relation with the aim to arrive at an accurate analytical semi-process-based breach model. This has been done for a situation of a free breach flow and a drowned breach flow.

Free breach flow
In the analytical model, water is assumed to flow through the breach and enter an enclosed polder with surface area A polder . Before the breach flow drowns, the spatial averaged rate at which the water level in the polder rises is given by here b (m) denotes the breach width, Q (m 3 /s) denotes the discharge through the breach, h out (m) denotes the outside water level, h breach (m) denotes the invert level of the breach, and h polder (m) denotes the water level in the polder. For the analysis, the outside water level is assumed to be relative to the breach invert level. Equation 13 shows that the breach discharge is linearly dependent on the breach width. The rate of breach growth is a function of the wall shear stress, which is assumed to be quadratically dependent on the flow velocity through the breach. The flow velocity through the breach has been calculated on the basis that the breach invert level equals the polder level. The average flow velocity through the breach (u) is now given by The wall shear stress relation now follows from Mannings' equation. Substituting the velocity in Mannings' equation by Eq. 14 and substituting the hydraulic radius in Mannings' equation by the depth of water in the breach give where n (s/m 1 3 ) denotes the Manning roughness coefficient. The lateral breach growth rate is given by 2 times the displacement rate of the breach faces. This finally leads to the following expression for the rate of breach growth where m c = √ 0.7m and c 1 follows from the empirical fit to the process-based erosion relation (Eq. 8). For a free breach flow, the breach growth rate is constant for constant hydrodynamic loading conditions. Integrating Eq. 16 with respect to time now gives the breach width as a function of time b(t) before drowning (Eq. 17).
Substituting the solution for Eq. 17 in Eq. 13 now gives the relation for the rate of water level rise, which is used to determine the onset of drowning.
where b 0 refers to the initial breach width. The following section deals with the determination of the breach growth rate during drowning.

Drowned breach flow
Once the water level in the polder h polder (Eq. 18) equals 2∕3h out , the flow becomes affected by the downstream water level. The breach flow is no longer super-critical, and the breach no longer widens linearly with time. During the drowning stage, it has been assumed that the analytical displacement relation is still valid. Due to the steep side walls of the breach, gravity-induced dilation would still affect the process of breach widening even when the flow has stopped. This process continues until a large slump failure of the side walls leads to a stable soil configuration. Studying what happens when the breach flow has seized is, however, of little interest as it has no effect on the flow through the breach. Drowning of the breach flow causes a decrease in flow velocities and breach discharge. The breach is assumed to be present over the full height of the embankment. The corresponding reduced rate at which the water level in the polder rises due to the flow through the breach is now approximated by A rise in the water level in the polder h polder gives a drop in the driving head of water h . Equation 19 can now be written in terms of a change in water level difference by defining h polder = h out − h (m). h out is thereby assumed not to vary during breach formation. This leads to the following differential equation for the head difference over the breach.
If the water level h out varies with time, then a time discretized outside water level could be applied. The maximum value for h is given by 1 3 rd of the outside water level as this is the level at which the breach flow starts to drown. From that point onwards, the water level difference only grows smaller. In order to develop an analytical solution for the rate at which the breach grows and the water level rises, the water depth in the breach has been given by h out − h(t) ≈ 0.83h out . This assumption gives a maximum error of 17%. Inherent to the applied weir equation is the assumption that the pressure distribution in the breach is hydrostatic. Due to flow contraction, the mean pressures in the flow could become lower than hydrostatic leading to an error in prediction of the discharge. This error could increase to as much as 40% (Nortier and De Koning 1991). The 17% error is of the same order as the error due to flow contraction and therefore considered allowable to solve the equations.
As an in-between step to further linearizing the equation, a new variable u was defined as u 2 = h(t) leading to here b(t)0.83h out denotes the cross-sectional area of the flow through the breach. As shown in Eqs. 20 and 21, the rate in which the water level rises is a function of the breach growth rate. The breach growth rate follows from twice the displacement rate given by Eq. 8 after substituting the Manning equation for the wall shear stress. The hydraulic radius in the Manning equation has been replaced by the estimated depth in the breach given by 0.83h out which leads to where u = The analytical solution for u and b follows from the analytical solution of Eq. 24 which is given by The coefficients a and b thereby have the dimensional units m − 1 12 . The first condition for a and b follows from Eq. 22 after substituting that at t = 0 h = 1∕3h out The second condition is given by Eq. 23 and is

3
With u(0) = √ 1 3 h out , this gives The breach width predictions that follow from Eqs. 17 for the non-drowned flow and 27 for a drowned breach flow have been validated. Details on this validation are given in the next section.

Validation
Although the sophistication and predictive capability of breach models have improved over the past years, the performance of breach models is often misinterpreted. Morris et al. (2008) mentioned that breach models are usually validated and calibrated against a single breach event because of the lack of high-quality datasets available for calibration purposes. Consequently, the ability of the model to reproduce the breach event is inherently quite good. An effort has been initiated as a part of the Dutch research project SAFELevee to collect and organize all available data with respect to levee performance into an online accessible database in support of the development of breach models. The database contains data on real levee failures and experimental data. To evaluate the accuracy of the approach, breach width predictions from the analytical model have been validated against data from this database. The polder area, outside water levels, and displacement parameters have been entered into Eq. 18 to determine the onset of drowning. When the flow starts to drown the breach width prediction shifted from the use of Eqs. 17 to 27. The breach width prediction from Eq. 17 at the onset of drowning thereby forms the initial breach width in Eq. 27. Model coefficients were derived to extend the usability to cases for which little information is known on the construction of the levee. Tables 2 and 5 together contain information on 66 time steps during various breach events which have been used for validation of the analytical breach growth equation. These cases address breach formation due to drowned flows and non-drowned flows. The data have been selected on the basis that the levees were constructed of dilatant sandy material such that the process-based displacement relation could be applied. The new model has furthermore been benchmarked against the current standard for sand dikes given by the HIS-OM model of Verheij and Van der Knaap (2003) which states that the breach width b(t) is given by

3
here H (m) denotes the initial difference in water level between the polder and the outside water level, and t denotes the time. The critical velocity u c has been set to 0.2 m/s (Verheij and Van der Knaap 2003). Table 2 contains experimental data of large and small-scale breach experiments that have been performed. Table 3 contains information on those cases Here b w (m) denotes the breach width, T (h) is the time at which the breach width has been reached, T 0 (h) is the start of breach development, A polder (m 2 ) is the known size of the polder which is flooded, Z w (m) is the water level driving the breach flow A.D. (Above Datum), Z p (m) is the polder level A.D., b 0 (m) is the initial breach width, and Z b (m) is the bed level AD 1 3 from Table 2 of which the soil parameters could be estimated. The displacement coefficients m and c 1 here follow from the product of Eqs. 10, 11, and 12. The other cases from Table 2 have been used to calibrate an overall displacement factor m = 0.0002253 and displacement coefficient c 1 = 0.008 (Eq. 8) to extend the use of the model to those cases for which insufficient information is available to determine these parameters. The displacement coefficient c 1 thereby equals the average of the displacement coefficients used for those cases for which sufficient information was available. Table 4 contains the displacement coefficients that have been used in the model validation alongside with the sources of the data. Labels thereby indicate the link with the other tables. Due to the higher level of accuracy, data from Table 2 have been used for the probabilistic analysis of the model. Table 5 contains data on the back-analysis of real breach events. Often insufficient  Tables 2 and 3 information was available on the hydrodynamic loading or soil parameters for these cases. The data in Table 5 are therefore expected to have a lower level of accuracy. The datasets contain information on breaches that initiated due to piping, internal erosion, or overtopping. For those cases in Table 2 for which the outside water level varied during breach formation, the outside water level has been discretized. The corresponding breach width has thereby been set as initial breach width at the start of the new water level. Here b w (m) denotes the breach width, T (h) is the time at which the breach width has been reached, T 0 (h) is the start of breach development, A polder (m 2 ) is the size of the polder which is flooded, Z w (m) is the water level A.D. driving the breach flow, Z p (m) is the polder level A.D., b 0 (m) is the initial breach width, and Z b (m) is the bed level A.D

3
For those cases from Table 2 for which the water level remained relatively constant, several moments in time were taken to evaluate whether the predicted breach growth matches the measured breach growth. For the data in Table 5, the initial water level was used to estimate the breach growth as no more information was available. More information on the breaches is given in the following sections. For those cases for which no polder size is given, either no drowning of the breach flow occurred, or the downstream water levels have been given for later stages of breach development in the tables. Data from Tables 2 and 5 have been used to verify the applicability of the process-based displacement equation to predicting breach growth rates. The displacement parameters used for validation are given in Table 4. For those cases for which sufficient soil parameters were known (Table 3), the ratio of the predicted and measured breach width is depicted in Figure 9. The figure was obtained from sorting the ratios from low to high and normalizing the range of data points by the total number of data points. The R 2 coefficient for the match between the breach width predictions following from the model and the data from Table 2 is 0.62, indicating that the model captures 62% of the variability of the dependent variables. For the predictions with the HIS-OM method (Eq. 32), the R 2 value is −60 indicating that the prediction is worse than when only the averaged value of the data would have been used. When both the historical data and the experiments are used, then both the new method and the HIS-OM method give negative R 2 values. For breach widths smaller than 160 m, the new method thereby does give a slightly better prediction than the HIS-OM method R 2 = −18 versus R 2 = −29 . In Fig. 10 the predictions of the new model are set out against the values from Tables 2 and 5. Applying a factor of 1.5 to the breach width predictions results in a 95% certainty that the predicted breach width exceeds the real breach width as is indicated in Figs. 10 and 11

Discussion and conclusions
The process-based analytical breach growth relation presented herein is based on a process-based descriptions of both the flow and the erosion process. Figure 10 shows that by combining the two process-based descriptions the new method provides a significant improvement in accuracy over the HIS-OM method when comparing the predictions against the output of experiments (Eq. 32). This is also supported by the R 2 coefficient of, respectively, 0.62 in comparison to the current method whose R 2 coefficient falls well outside the range of [0, 1]. The best fits are obtained when comparing against experimental data. When the historical failures are included in the R 2 analysis, both methods perform poorly. The source of this poor fit is likely to be related to the model assumption that the water level remains constant during the breach formation. The variations in driving head  Fig. 11 Data versus predicted values after applying a factor of 1.5 to the predictions by the new model outlined in this paper . Here "det" refers to the use of the displacement parameters which are based on soil parameters, and "calib" refers to the use of the overall calibrated displacement parameters 1 3 of water are often unknown for the historical cases but have a significant impact on the prediction. It is recommended to apply the model in a time discretized manner for cases whereby the outside water level is likely to change significantly with time during breach formation, as was done during the validation against the experimental data. Unfortunately, for these historic data, this information is rarely available. For validation of breach models, it is therefore essential that the data are of sufficient quality as indicated by the difference in results between the experimental data and the historical data. The process-based erosion model underlying the breach growth relation accounts for the impact of the permeability and degree of compaction in the breach growth prediction. At the moment, the applicability of the model is limited to predominantly sandy levees as the erosion relation has not yet been validated for other materials. However, the impact of other materials could theoretically be accounted for via the soil parameters in the erosion model. For the expansion of the applicability of the breach model, further validation of the presented erosion relations is recommended. For those cases for which no information is known on the soil condition in the levee, an displacement factor m = 0.0002253 and displacement coefficient c 1 = 0.008 have been shown to provide reliable initial estimates for rate of displacement of sandy breach faces due to erosion (Fig. 9). Upon extension of the model to different types of soil, it is recommended to recalibrate these parameters for the various soil types. In general, the new method gives more accurate predictions than the HIS-OM model against experimental data, or for breach widths smaller than 160 m. By applying a factor of 1.5 to the predicted breach width values, only a 5% chance remains that the model under-predicts the breach width while maintaining its improved accuracy (Figs. 9 and 11). Further research is recommended to the impact of a more accurate breach growth model on flood risk predictions.

3
where ′ (N/m 2 ) denotes the effective stress, is the angle of internal friction, and C (N/ m 2 ) is the cohesion of the material subject to erosion. Shear deformation of dilatant soil is accommodated by an increase in pore volume. The pressure gradient required to accommodate the associated inflow of water gives a decrease in pore pressures in the soil, an increase in effective stresses ′ (N/m 2 ), and an increase in shear resistance. Dissipation of the low pore pressures by the inflow of water induces continuous shallow shear failures of soil, resembling the process of erosion. The constitutive equations that describe this process have been derived for the case of a flow past a vertical wall. The x-direction thereby coincides with the main flow direction past the x, z-plane. The z-coordinate has been taken parallel to the soil-water interface pointing upwards, and the positive y-direction is pointing into the soil (Fig. 12). The positive velocity vectors u (m/s), v (m/s) and w (m/s), respectively, coincide with the positive x-, y-and z-coordinate directions.

Mass balance equation
For the case of incompressible particles and incompressible water, the mass balance equation for water in a porous medium is given by where n denotes the porosity and subscript w refers to water. Similarly, the mass balance equation for solid particles in a porous medium is given by where subscript p refers to particles. Adding Eqs. 34 and 35 gives The last three terms in Eq. 36 give a change in porosity due to volumetric strain and equal 1 1−n n t (Verruijt 2006). The terms n(u w − u p ) , n(v w − w p ), and n(w w − w p ), respectively, denote the specific discharges q x (m/s), q y (m/s), and q z (m/s). Substituting the expressions for the specific discharges in Eq. 36 gives Flow-and gravity-induced shear stresses acting on the x, z-plane are transferred into the y-direction. The thickness of shear layers is assumed small compared to the characteristic length scales of the surface on which the shear stresses act. This allows for the problem of pick-up of sediment to be regarded as the problem of a continuous displacement of a soil plane in the y-coordinate direction. To better describe this displacement, a coordinate system was defined by y = − ct (Jensen and Finlayson 1980;Verruijt 2006), where c (m/s) is the rate of displacement of the x, z-surface in the y-coordinate direction. No significant spatial variations in shear stress are assumed to occur in the x-and z-coordinate directions. The depth-integrated mass balance equation for water (Eq. 34) now becomes Over the depth d the increase in porosity is assumed constant (Verruijt 2006;Van der Schrieck 2006). This corresponds with the conclusions of Casagrande (1936) who observed that samples subjected to large deformations reach the same porosity. It is also a requirement for integrating the momentum balance equations as will be discussed in Section A.2. At depth d shear stresses equal the shear resistance. Over the short time duration required to remove a soil layer with depth d, denoted by t = d∕c , the degree of deformation of the soil skeleton below the failure depth d is assumed negligible compared to the deformation over d. Under this assumption v w | =d = v p | =0 = 0 . Dilation is thereby assumed to be solely accommodated by the inflow of water through the soil surface given by the x, z-plane. Integrating Eq. 38 over from = 0 to = d gives Rewriting Eq. 35 in terms of a moving coordinate system and integrating it from = 0 to = d similarly give the depth-integrated mass balance equation for the particles whereby it has been assumed that v p | d = 0 . The depth-integrated mass balance equation of the total porous medium similarly follows from rewriting and integrating Eq. 37, giving here n| =d has been assumed to equal the initial porosity n 0 , whereas n| =0 is assumed to equal the critical porosity of the fully dilated material n loose (Van der Schrieck 2006).

Momentum balance equations
The general momentum balance equations for a flow in soil comprise the momentum balance equations of the particles and the pore water. Below the momentum balance equations for a wall moving at rate c have been given under the assumption of an infinitely long slope and The contribution of the surface shear stress to the lateral erosion rate follows from the horizontal momentum balance equation whereby x | =0 > 0 . Due to the presence of a horizontal shear stress component acting on a vertical soil surface, a failure depth d (m) needs to be defined. The lateral erosion rate is related to the rate of inflow of water into soil. A high rate of inflow requires a high pressure gradient. The corresponding low pore water pressures give a high shear resistance and a low failure depth. Consequently, c s is inversely related to the failure depth d. The average excess shear stress available for accelerating the sheared soil layers depends on the average resistance against shearing over d. As the failure depth and erosion rate are inversely related, a failure depth d exists for which the average dilatancy-induced shear resistance is maximum, and the momentum available for accelerating the particles is minimum. This situation corresponds with a maximum in effective stresses, averaged over the failure depth. The linear decrease in velocity profiles v w (m/s) and v p (m/s), associated with the spatially constant rate of dilation, is given by Substituting Eq. 56 in Darcy's law and twice integrating with respect to and dividing by d give the average shear induced drop in pore water pressure over d (m). For dP d | =d = 0 and P| =0 = 0, the integral of the drop in pore water pressure, denoted by P f (N/m 2 ), is given by.
Effective stresses are given by the soil pressure minus the water pressure. The failure depth-averaged contribution of pore water pressures to the increase in effective stress,  , and hence shear resistance, follows from substituting = d and dividing by d. The failure depth-averaged effective stress is now given by The depth d for which the averaged effective stress is maximum follows from including the effects of c z (Eq. 55). The shear resistance provided by the gravity-induced dilation equals = ( s − w )cos( )d , for = 0 (Eq. 55). The constant inflow thereby also induces an additional exchange of momentum between the flow and the soil. Accounting for the effects of c z in Eq. 52 gives for = 0 where the contribution of c z (m/s) on the shear resistance is accounted for in the first term of the denominator. Substituting Eq. 59 in Eq. 58 gives a relationship for the average effective stress as a function of the rate of displacement of a vertical soil surface due to erosion c (m/s). Figure 2 shows this relationship for a shear stress of 0 = 100 N/m 2 and K s = 1E − 4 m/s, for a non-cohesive dilatant soil. As illustrated in Fig. 2 a maximum average effective stress is found for a specific value of the boundary displacement rate c (m/s). This maximum corresponds with a minimum in momentum available for the acceleration of soil particles at the surface and resembles the most stable situation.