A multilayer soil approach for seepage process analysis in earthen levees

The analysis of seepage through earthen levees is fundamental to investigate the levee safety condition, since that the seepage phenomenon can often lead to piping causing the collapse of the levee (Saghaee et al., 2012). This work proposes a new approach for identifying the location of the saturation line through the earthen levee by using a multilayer soil approach that assumes the levee body and the levee foundation as two different layers, characterized by using a different hydraulic conductivity parameter. The main purpose is to investigate the difference in the predicted saturation line identified through the proposed method and the one based on the single layer analytical solution of the seepage flow suggested by Marchi (1961). On comparison of the saturation lines provided by two different approaches, the ratio between the river water head and the groundwater thickness or depth increases; allowing the determination of a ratio threshold above which the Marchi equation should no longer be applied. The proposed multilayer (levee and foundation) soil approach is also used to assess the seepage vulnerability of the methodology proposed by Barbetta et al. (2017). The results show that the differences between the seepage vulnerabilities evaluated using the multilayer approach and the single layer method are quite low. Although notably, significant modifications of the seepage probability are observed when the distance between the ground level and the groundwater table changes.


| INTRODUCTION
Floods are among the most widespread and destructive natural hazards worldwide. The severe effects of flood events are often mitigated through structural measures alongside complementary non-structural measures.
Earthen levees represent one of the most common structural measures used to reduce flood risk. The levees structural vulnerability becomes increasingly concerning when combined with a poor perception of the flood hazard by people who live in the floodplain, protected by the levee systems. This misconception of perceivable safety can encourage the growth of human settlements and consequently can lead to an increase in the number of people living and affected by floods in an area at risk from a potential levee failure (Ludy & Kondolf, 2012). In the last decades, the scientific community has faced the flood risk management moving gradually from a structure-oriented approach towards a risk-based approach (Apel et al., 2006;Kellens et al., 2013;Schumann, 2017;Vrijling, 2001;USACE, 1996USACE, , 1999. In the first case, only the efficiency of the mitigation structures was evaluated by considering their capability to remain undamaged during flood event characterized by assigned return periods. On the contrary, the risk-based approach also considers damages to people, human health, public and private goods, cultural goods and ecosystem services, providing the risk assessment as the combination of exposure and vulnerability. This was required by the growing demand from the authorities involved in the control and mitigation of risk to develop comprehensive methods addressed to the verification of effectiveness of flood defenses structures, such as levees, and the quantification of expected damages and residual risk. A proper assessment of the residual risk in the areas protected by the earthen levees is a key element with regard to flood risk mitigation policies, as stated by the "Floods Directive" (Directive 2007/60/EC). This can help decision-makers to assess alternative actions in territorial planning and emergency strategies, to aim of improve public understanding and acceptance of a possible residual risk, involving it in mitigation strategies. For example, in 2013 an organization from six countries (France, Germany, Ireland, Netherlands, United Kingdom and United States) published the International Levee Handbook (ILH) with the key point of accepting the issue of levee failure and collecting a consistent classification and treatment of deterioration and damage mechanisms. The document provides a good practice guidance for the management, design and construction of earthen levees that, even if built according to these guidelines, could be affected by failures. However, flood maps delineation in some countries, like in Italy, typically assumes that levees maintain its original functionality during flood events and, hence, the flood risk does not consider the possibility of failure of the longitudinal defense structures (Mazzoleni, Bacchi, et al., 2014;Nagy & Toth, 2005;Sellmeijer & Koenders, 1991;Saghaee et al., 2012). However, possible damages can compromise the effectiveness in some countries, like in Italy, and induce a "residual risk" to be managed in order to avoid flood risk underestimation (Castellarin et al., 2011;Ludy & Kondolf, 2012;Wilby et al., 2008).
The river levees failure can be due to several processes: scouring and sliding of the foundation, overtopping, seepage/piping through levee's body and foundation (ASCE, 2011;Huang et al., 2014;ICOLD, 2013;Ojha et al., 2001;Serre et al., 2008;Sills et al., 2008). In particular, Foster et al. (2000) showed that piping is one of the most common mechanisms causing levee failure. In fact, when the seepage phenomenon starts the resistance of the soil particles decreases and the internal erosion triggers, causing the removal of the soil particles by the flow and the beginning of the piping process. This mechanism is not clearly visible and although there are sometimes surface level indicators such as wetting, slumping or scouring, on the whole this is very hard to sport and predict.
Various models in the literature have been proposed to identify the location of the saturation line: from more complex models, based on the solution of the Richards equation, to more simple approaches applying the Darcy's law (Bardet & Tobita, 2002;Casagrande, 1940;Cividini & Gioda, 1990;Iterson, 1917;Kozeny, 1931;Mishra & Singh, 2005;Palladino et al., 2020;Schafferank, 1917). However, in many cases these models do not consider the uncertainty related to the hydraulic parameters of the soils, such as hydraulic conductivity or soil porosity. Therefore, their application is limited to single and well-known case studies for which the information regarding permeability, grain distribution, and thickness of soil layer are required (Mazzoleni, Bacchi, et al., 2014;Rice & Polanco, 2012;VanBeek et al., 2011). To improve analysis on catchment scale levee reaches, the development of expeditious procedures for identifying the most vulnerable levees is the desirable solution.
To this end, Vorogushyn et al. (2009) studied the levee vulnerability using fragility curves; incorporating the uncertainty of the soil parameters and defining a simple vulnerability index based on the hydraulic and geometric characteristics of the levee.  randomized the geometry of the levee, by a Monte Carlo approach, to achieve fragility curves estimation of potential failure. Barbetta et al. (2017) proposed a procedure for homogeneous levees (Marchi, 1961), define the levee body vulnerability to seepage, using measured hydraulic and geometric characteristics. This method accounts for the uncertainty in the estimate of soil hydraulic conductivity through a probabilistic approach based on a Monte Carlo sampling method. In particular, it assumes the same level for both the ground and the groundwater table, developing a more sophisticated analysis that assured a higher level of safety. Barbetta et al. (2017) identified the probability of occurrence of a critical condition that would allow the soil particles' erosion through the levee body and potentially lead to structural failure. The critical condition was met when the saturation line intercepts the landward side of the levee.
The analysis was based on two main assumptions, firstly that the groundwater table was at the same level as the ground and secondly, the soil was characterized by the same hydraulic conductivity parameter as the levee body and foundation.
In this context, this work proposes a different approach for identifying the location of the saturation line within the earthen levee body and the foundation by using an equation that considers different hydraulic conductivity parameter value for the levee body than the foundation.
The main purpose of the study is to investigate the difference between the saturation line, exploiting the continuity and flow equations and fixing a different characteristic hydraulic conductivity for each layer (levee and foundation) of those used by Marchi (1961). Marchi derived the solution assuming a steady state condition, that is, stationary water level in the river for a fixed period. The levee body and its foundation are characterized as a uniform layer with uniform properties. Marchi's method is applicable when the thickness of the groundwater depth, H 0 , below the levee is much greater than the water depth in the river, h 1 (see Figures 1 and 2).
This study looks to identify a validity threshold for the ratio h 1 /H 0 , beyond which the Marchi's equation should no longer be used. Moreover, the need to consider the groundwater table below the ground level is further considered. Finally, the impact of using the multilayer F I G U R E 1 (a) Conceptual scheme; (b) sketch of earthen levee and the foundation soil and representation of the variables adopted for the definition of the saturation line with the multilayer soil equation F I G U R E 2 Vulnerability index to seepage for a levee with known geometry (for symbols see text) soil equation on the levee seepage vulnerability assessment is investigated.
Section 2 presents the mathematical basis of the new multilayer soil approach and describes the expeditious procedure for seepage vulnerability assessment . Section 3 presents and discusses the obtained results and Section 4 outlines conclusions from the analysis.

| Multilayer soil approach
In order to calculate the location of the saturation line within an earthen levee foundation and body, the flow equation along with the continuity equation are applied (see Figure 1a): where q is the flow rate, h 1 is the hydraulic head in the river; a is the distance between the groundwater table and the ground level; H 0 is the groundwater thickness ( Figure 1b), n is the soil porosity; x and y represent the reference coordinate system; H is the saturation line height along the horizontal distance x. V LB and V F are the Darcy's flow velocity in the levee body and in the soil foundation, respectively, defined as: K LB and K F are the soil hydraulic conductivities of the levee body and the foundation, respectively (see Figure 1a). Using Equations (3a) and (3b), Equation (1) can be expressed as: Use of Equation (4) in Equation (2) leads to the expression: Considering that ∂H 0 ∂x ¼ 0 and at a permanent state of flow the coefficient is a constant, the saturation line referring to a multilayer soil can be inferred similarly to Marchi's equation for a single layer : where D is the duration of the flood event and erf is the error function, that is, twice the integral of the Gaussian distribution with zero mean and variance equal to 0.5. The Marchi's equation Camici et al., 2017;Marchi, 1961) is given by: This solution for the identification of the saturation line in earthen levees considers the steady state reached after the boundary conditions (i.e., river water head, h 1 ) are kept constant for a period equal to the entire flood duration, D.

| Vulnerability index and seepage probability
The delineation of the saturation line in the levee body and foundation allows an understanding of the levee condition in relation to seepage vulnerability. As previously mentioned, Barbetta et al. (2017) assumed that the critical condition is met when the saturation line intercepts the levee landside. If the saturation line is contained within the levee, the levee failure due to seepage process is assumed avoided. The interception of the landside identifies a necessary but not sufficient condition that would allow the soil particles' erosion through the levee body and potentially could lead to the levee piping and failure.
The levee vulnerability is assessed through a simple vulnerability index that compares the maximum length of the seepage line, x max , representing the distance at which the seepage line intercepts the ground level, with the length at the levee base (L) (see Figure 2).
Specifically, the vulnerability index is defined for a dimensionless geometry by using the quantities x* = x/L and H* = H/H s , where L is the bottom width of the levee and H s is the levee depth, as Palladino et al., 2020): In Equation (8), x Ã max is the dimensionless maximum length of the seepage line computed as x max L ; x 0 Ã ¼ H s Àh 1 L cotα, with H s equal to levee depth and α equal to slope of the levee riverside (see Figure 2 for all details).
The maximum length of seepage, x max , is quantified imposing in Equation (6) H = a, that is, solving the equations to identify the distance at which the saturation line intercepts the ground level.
If I v <0 the seepage line is enclosed within the levee and, as a consequence, the levee is assumed to be safe. If I v = 0 the seepage line crosses the landside levee toe, marking the threshold condition for seepage. When I v >0, the seepage line intercepts the landside slope and the structure is potentially affected by piping. Barbetta et al. (2017) wrote Equation (7) using dimensionless characteristics of the levee as: where δ ¼ H 0 D L 2 n . The description of the seepage flow within the levee is considerably affected by the uncertainty in the soil parameters estimate, such as the soil porosity n and the hydraulic conductivity, K LB . Previous studies showed that the location of the saturation line is much more sensitive to the hydraulic conductivity Camici et al., 2017). For this reason, the seepage probability is evaluated by identifying K LB through a Monte Carlo approach (Pohl, 2000;USACE, 1999). In particular, 1000 values are randomly sampled from a lognormal distribution with mean μ K ¼ 10 À5 ms À1 and standard deviation σ K ¼ 25μ K Vorogushyn et al., 2009).
For a specific value of the δ parameter, I v can be assessed for each of the 1000 values of K LB . The resulting curves ("fragility curves"), one for each δ value, are estimated for a specific ratio h 1 /H s (e.g., Figure 3a). The cumulative probability of I v values allows the assessment of the levee vulnerability to its seepage and is shown in Figure 3b ( Barbetta et al., 2017).

| Performance measures
The comparison between the saturation line assessed through the Marchi equation (Equation (7)) and the one identified with the proposed multilayer model (Equation (6)) is based on the root mean square error (RMSE) and the Nash-Sutcliffe efficiency coefficient (NSE): where H Marchi and H multi refer to the Marchi solution and the multilayer model solution, respectively; N is the number of values where the saturation line is quantified by using Equation (6) or Equation (7). H mean multi is the mean of H multi values.
Notably, this is not a comparison between observed and computed values, but a comparison between the results of two different approaches for the saturation line estimate.
The comparison is based on the difference between the maximum length of seepage, x max , assessed by the two models, quantified both as absolute value and as a percentage:

| RESULTS AND DISCUSSION
The new approach based on the multilayer soil equation is applied to identify the location of the saturation line through the earthen levee body and the foundation. First, the difference between the saturation line identified through the proposed method and the one based on the analytical solution by Marchi  is estimated and a validity threshold of the ratio h 1 /H 0 is identified. This threshold represents the value beyond which the Marchi equation is characterized by nonnegligible differences with the solution provided by the proposed multilayer approach. This first analysis is carried out assuming a = 2 m, which corresponds to a position of the groundwater table below the ground level.
Second, the effect of adopting the multilayer equation for levee vulnerability assessment is investigated; specifically, two different values of a (a = 0 m and a = 2 m) are considered to enable assessment of the impact of the water table location below ground level.
The analyses are carried out for a levee characterized by: L = 9.3 m, H s = 2.3 m, α = 0.405. The flood duration, D, is first set equal to 12 h for the saturation line analysis and then increased up to 72 h for the levee vulnerability assessment. As previously mentioned, this study is carried out assuming stationary conditions, that is, the water head h 1 in the river is constant for the flood period equal to D, and we investigate the reached steady state. For the first analysis, the river water head, h 1 , is assumed constant and equal to 1.8 m while the groundwater thickness, H 0 , is varied within a range of 1.8-35 m. Specifically, h 1 /H 0 values are changed with a constant interval equal to 0.1 and the relative values of H 0 are defined in the above-mentioned range.
Further analyses are carried out considering a different levee geometry, characterized by L = 40 m, α = 0.405, H s = 10 m, investigating the effect on the h 1 value. In this analysis h 1 values are varied up to 10 m and H 0 ranged between 10 and 200 m.

| Saturation lines analysis (a = 2 m)
Six different configurations are considered for the identification of the saturation lines. Table 1 summarizes these for each configuration used (i.e., Marchi or multilayer model) and provides the values of the hydraulic conductivity. Notably, K can assume different values, equal to 10 À7 and 10 À5 ms À1 , for the levee body and the foundation in the multilayer approach, whereas for the Marchi solution a single parameter K is considered.
In order to apply Equations (6) and (7) for the levee characterized by H s = 2.3 m, the distance between the groundwater table and the ground level is assumed equal to 2 m (a = 2 m), the soil porosity n = 0.2, the flood duration D = 12 h, the river water head h 1 = 1.8 m and the groundwater depth H 0 variable in the range 1.8-35 m.
The comparison of the saturation lines estimated by the six configurations is shown in Figures 4-7 for different values of the ratio h 1 /H 0 .
When h 1 is much lower than H 0 , the saturation lines provided by Equations (6) and (7) are very similar (see Figure 4) with the dashed black line with circles (C1) almost overlying the green line (C4). Similarly, when K = 10 À7 ms À1 is considered, the black dashed line (C2) and the cyan line (C6) almost overlap. When the multilayer soil approach is applied with different K value for the levee body and the foundation (C3 and C5), the saturation line is found to approach the levee landside more than when compared with the simulation C6 (lower and constant K = 10 À7 ms À1 ). When C3 and C5 results are compared with the simulation C4 (higher and constant K = 10 À5 ms À1 ), the green line (C4) is found closer to the levee landside.

Configuration ID Equation
Hydraulic conductivity (ms À1 ) Moreover, it is worth noting that the safety assumption assumed by Marchi, that is, the same hydraulic conductivity for levee body and foundation, is shown in Figure 4 for the condition represented by the magenta line and the green line (configurations C5 and C4).
When the value of h 1 /H 0 increases, the saturation lines provided by the two approaches for K = 10 À5 ms À1 (C1 and C4) become more distant from each other and the difference increases as the height value of h 1 /H 0 increases (see Figures 5 and 6), up to a maximum value equal to 1 (see Figure 7 for h 1 /H 0 = 1). Specifically, the dashed black line with circles (C1), representing Marchi's saturation line, identifies a maximum seepage length which is lower than the one indicated by the green line (C4) (Figure 7).
As is expected, when the hydraulic conductivity is low, that is, K = 10 À7 ms À1 (C2 and C6), both saturation lines are almost vertical (Figures 4-6) and when the hydraulic conductivity increases, the saturation line is located much closer to the levee landside (see the green line referring to C4 and the black dashed line referring to C2) indicating a higher vulnerability condition. The comparison between the saturation lines achieved using the two approaches also points out that the difference between Marchi and multilayer solutions increases when the hydraulic conductivity is increased, that is, equal to 10 À5 ms À1 .
The differences between the saturation lines identified through the two methods are quantified by the RMSE, the NSE coefficient and the difference between the maximum length of seepage. The performance measures are computed by comparing the configurations C1 and C4 (K = 10 À5 ms À1 ) considering h 1 = 1.8 m and a = 2 m, and H 0 varying in the range 1.8-35 m. Figure 8a shows that the RMSE increases with the value of the ratio h 1 /H 0 . Figure 8b demonstrates how the NSE value decreases while the ratio value is approaching to 1 (i.e., h 1 = H 0 ). Both RMSE and NSE measures' behavior indicates that the differences are significant for high values of the ratio, while they become negligible when the river water head is much lower than the groundwater thickness, as assumed by Marchi. The results depicted in Figure 9 indicate that the difference between the two maximum seepage lengths increases while the ratio increases (in reference to Δx max ).
The selected performance measures are also estimated for other combinations of K, a, and h 1 values (K LB F I G U R E 6 As for Figure 4, but for h 1 =H 0 ¼ 0:5 F I G U R E 4 Comparison between the saturation lines estimated through Marchi equation and the multilayer soil method: h 1 =H 0 ¼ 0:1 (a = 2 m, L = 9.30 m, H s = 2.30, D = 12 h). The colored lines are all referred to solutions obtained with the multilayer model (Equation (6)) F I G U R E 5 As for Figure 4, but for h 1 =H 0 ¼ 0:3 F I G U R E 7 As for Figure 4, but for h 1 =H 0 ¼ 1 and K F varying within the range of 10 À9 to 10 À4 ms À1 , h 1 variable from 1.8 to 10 m and a ranging from 2 to 5 m), providing similar results here omitted for sake of brevity. This last analysis is carried out for the levee characterized by H s = 10 m.

| Applicability of Marchi equation: Identification of the threshold ratio value
On the basis of the results in Figure 9, a threshold value of the ratio h 1 /H 0 above which the saturation lines identified by Marchi equation are significantly different from those provided by the multilayer soil approach can be identified. Specifically, for fixed values of K, h 1 and a, where we consider the difference as acceptable when Δx max is lower than 20% and we assume the corresponding ratio as the applicability threshold ( Figure 10). The identified ratio value is equal to 0.57, for K = 10 À5 ms À1 , a = 2 m and two values of river water head, h 1 = 1.8 m and h 1 = 6 m. Moreover, the same analysis carried out for different values of h 1 , a and K provides always the same critical value equal to 0.57 (Table 2). Therefore, when h 1 /H 0 <0.57, the two F I G U R E 8 Comparison of the saturation lines estimated with the two method: (a) Root mean square error; (b) Nash-Sutcliffe efficiency coefficient F I G U R E 9 Difference between the x max , the maximum length of the saturation line, computed through the two models: (a) in meters; (b) as a percentage F I G U R E 1 0 Identification of the threshold value of the ratio h 1 =H 0 saturation lines are similar and both methods are adequate; if the ratio is higher than 0.57, the results are no longer negligible in their difference and the multilayer soil approach should be used.

| Saturation lines analysis (a = 0 m)
A further analysis is undertaken to provide a preliminary investigation into the impact of the distance between the water table level and the ground level (a). For this, a comparison between the saturation lines is carried out also considering the water table at the ground level (a = 0 m). The results are similar to the ones provided by the previous analysis developed with a = 2 m in terms of saturation lines comparison. By way of example, the black dashed line with circles (C1) and the green line (C4) are nearly overlapping when h 1 /H 0 is equal to 0.1 (Figure 11a), whereas the two lines are different when the ratio is equal to 1 (Figure 11b). Similar results become apparent when comparing the other lines as for the study carried out for a = 2 m.
It is worth noting that the same analysis carried out for different geometric characteristics of the levee (L in

| Vulnerability index and seepage probability
A further analysis for comparing the results of the two models is based on the assessment of vulnerability index to seepage (Equation (8)). This study is necessary to understand the impact of using the two different methods in the quantification of the seepage occurrence probability. The comparison is based on the vulnerability index evaluation and on the fragility curves delineation . We consider 1000 values of K randomly sampled in the range of 10 À14 to 10 À3 ms À1 , following the lognormal distribution with mean μ K ¼ 10 À5 ms À1 and standard deviation σ K ¼ 25μ K . Using δ = 8.7 Â 10 4 ms À1 (L = 9.30 m; n = 0.2; D = 12 h; H 0 = 35 m), the fragility curves are derived by using both the methods. The comparison is carried out for two different locations of the groundwater table, a = 0 and a = 2 m. The results shown in Figure 12 indicate that the Marchi equation provides lower seepage probabilities. Specifically, when a = 2 m, where the continuous blue line in Figure 12a provides lower vulnerability index values and hence a lower seepage occurrence probability (3%Þ (see Figure 11b) compared with the outcomes of the multilayer equation 4% ð Þ (continuous red line in Figure 11). However, when a = 0 m a significant difference in the seepage vulnerability assessment can be observed. In this case, the seepage occurrence probability is quite high (≈30%) (see Table 3) and the difference between the two approaches increases up to 5%. It is important to point out that when the event duration D increases, the failure probability increases as well (Table 3). Moreover, the increase in the seepage probability provided by the proposed approach is not negligible for the longest duration, 72 h, and when h 1 /H 0 = 1.0, that is, when the difference between the two saturation lines is more evident.

| CONCLUSION
The "Floods Directive" (Directive 2007/60/EC) requires to develop methods to mitigate the flood risk, assuming the possible failure of the flood defenses structures, such as levees. In this context, this work proposes a new method to assess the saturation line location within earthen levee bodies and the foundation where a levee is sited, considering a multilayer soil approach. The approach allows the assignment of a different value of the hydraulic conductivity parameter to the levee body and the levee foundation.
The results of the proposed method are compared with the ones derived from the application of Marchi equation. The analysis considers a single soil,  homogenous layer, with a uniform hydraulic conductivity; which is assumed to be applicable when the groundwater thickness is much higher than the water head in the river. The analysis indicates that: (1) the difference between the two different methods saturation lines increases with the ratio h 1 /H 0 up to its maximum value (equal to 1). Specifically, the saturation line estimated by the proposed approach is always found to be closer to the levee land side, indicating that the multilayer soil method depicts the levee as more vulnerable; (2) the difference between the two saturation lines increases when the hydraulic conductivity is higher. Based on these comparisons, the study identifies for the ratio h 1 /H 0 a threshold value equal to 0.57, above which the assumption of the Marchi model cannot be considered valid for the identification of the saturation line.
Finally, the impact of using the proposed approach in the estimation of levee vulnerability to seepage is investigated. The Marchi equation always provides lower seepage probabilities and the difference between methods is higher when the groundwater table is at ground level. Specifically, the seepage probability is found very sensitive to the a values, equal to 2 and 0 m, with a mean increase of probability equal to about 30%.
It is suggested that future investigations could corroborate the results shown in this work, both considering different geometrical characteristics and thoroughly testing the saturation lines estimation through the use of reliable field data. Field data could be used to benchmark and test the capability of the proposed multilayer approach and further help identify the saturation line location in earthen levees.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.