## Abstract

In this research, SEEP2D and SEEP/W numerical models are used to simulate seepage through earth dams with internal cores. In order to evaluate the two models' performance, they were compared in cases with no, vertical, and wedge-shaped cores. SEEP/W was then used to study further cases due to its accuracy in drawing the phreatic line within the core zone. The effect of the core's characteristics on the amount of discharge, and the phreatic line's levels at the core's upstream and downstream faces were investigated. Four core types – vertical, wedge-shaped, upstream inclined, and downstream inclined – were considered. Different hydraulic conductivities, upper widths, and core slopes were also evaluated. The wedge-shaped core is the most effective of those investigated in reducing seepage discharge and the phreatic line's level at the core's downstream face, the vertical core came second. Design equations are provided for all the core shapes considered in the study.

## HIGHLIGHTS

Compares the performance of SEEP2D and SEEP/W numerical models in simulating earth dams with internal cores.

Evaluates the most effective earth dam core shape.

Covers downstream inclined cores, which have not been considered previously in the literature.

Investigates the effect of core characteristics on seepage discharge and the phreatic line.

Suggests design equations for the core types studied.

### Graphical Abstract

## ABBREVIATIONS

- b
core upper width (L).

- B
dam crest width (L).

- H
upstream water head acting on the earth dam (L).

- h
total head (L).

- h
_{d} datum (L).

- h
_{el} elevation head (L).

- h
_{p} pressure head (L).

- i
hydraulic gradient (dimensionless).

- k
hydraulic conductivity (LT

^{−1}).- k
_{s} saturated hydraulic conductivity (LT

^{−1}).- k
_{r} relative hydraulic conductivity (LT

^{−1}).- K
_{shell} dam shell's hydraulic conductivity (LT

^{−1}).- K
_{core} dam core's hydraulic conductivity (LT

^{−1}).- q
seepage discharge rate per unit length of the dam (L

^{2}T^{−1}).- R
^{2} coefficient of determination (dimensionless).

- S
core slope/inclination (horizontal:vertical) (dimensionless).

- V
Darcy's velocity (LT

^{−1}).- y
_{1} phreatic level at the core's upstream face (taking the dam's base as the datum) (L).

- y
_{2} phreatic level at the core's downstream face (taking the dam's base as the datum) (L).

## INTRODUCTION

Dams are built to store water for domestic and irrigation use, to regulate flow, reduce the dangers of floods and droughts, generate electrical energy, and so on. Earth-fill dams include fully homogeneous embankment dams and inhomogeneous embankment dams (zoned and diaphragm). A core is used to reduce seepage through the dam and lower the phreatic level, and, accordingly, reduce the height of the wetted portion on the dam's downstream face. The core material should have low permeability – for example, clay, silt, or silty clay – while the dam's external shell is usually pervious – for example, coarse sand. Earth dams have low rigidity, making them relatively likely to collapse because of hydraulic failure (e.g., overtopping and erosion of upstream face by waves' action), seepage failure (e.g., piping and sloughing), or structural failure (e.g., upstream and downstream slope slides), for instance (Garg 2005). A statistical study of dam collapses showed that 80% of such incidents were caused by water seepage through the dam body or foundations, due to the differential head between the dam's upstream and downstream sides (Garg 2005; Modi 2011).

Numerous studies have considered the problem of earth dams with cores, many being case studies of specific dams and sites. Ismaeel & Noori (2011),**s**tudied Duhok zoned earth dam using SEEP2D. The model's accuracy was tested using experimental and field data, and the results showed that it was acceptable. Hasani *et al.* (2013) studied Ilam earth dam, which has a central clay core, using SEEP/W, and evaluated the effect of the type and size of its mesh on the total flow rate. Mohammadi *et al.* (2013) studied Birjand Hesar Sangi earth dam, which has a core, using SEEP/W to optimize its geometry. Fattah *et al.* (2014) used SEEP/W to study Al-Adhaim zoned earth dam, evaluating the effect of the shell's permeability, and the low permeability core's location and thickness, on seepage. Abdulkareem *et al.* (2014) used SEEP/W to assess the safety of Duhok and Bawashaswar dams. Khassaf & Madhloom (2017) studied the seepage parameters of Khassa Chai earth dam, which has a central core, using SLIDE V.5.0 for different core permeabilities and thicknesses. Li *et al.* (2019) performed a case study on the Shizihe earth-rock dam using ABAQUS.

Other researchers have investigated the dam core effects experimentally. Rezk (1995) studied the effect of changing the relative permeability between the dam's core and body on seepage through an earth dam with a vertical internal core. Kanchana & Prasanna (2015) used different combinations of materials to form zoned earthen dams with central vertical cores and evaluated the behavior of the phreatic line. Salem *et al.* (2019) studied the effect of core permeability, width, base thickness, and penetration, for an earth dam with and without an internal core, both experimentally and, using SEEP/W, numerically. The numerical analysis was verified using the experimental model. Rezk & Senoon (2011) studied an earth dam with internal core on an impervious base using an analytical method, investigating the effect of the core's relative permeability on seepage and the drop of the phreatic surface.

Numerical models have also been used widely in studies considering seepage through earth dams with cores. Shakir (2011) studied a zoned earth dam with two core types; a vertical low permeability core and an inclined upstream core, using a finite element model to study seepage quantities and the location of the free surface, in relation to different core permeability, thickness and location scenarios. Fakhari & Ghanbari (2013) investigated an embankment dam considering two types of dam cores, a vertical core and an oblique core towards the upstream side, using SEEP/W, studying the effect of water height in the reservoir, dam crest width and the core's central angle. Li & Jianjun (2010) used Seep/W to analyze seepage through a cored dam under both steady and transient conditions. Zahedi & Aghazani (2018) used GeoStudio software to study earth-fill dams, considering vertical, inclined and diaphragm core types, with different slopes and geometries, and investigate the seepage behavior through the dam's body and the effect of the core's anisotropy. Jamel (2018) studied seepage volumes through a cored earth dam using SEEP/W to investigate different upstream and downstream slopes of the dam and core. Kheiri *et al.* (2020) used SEEP/W to study an embankment dam with a core, cutoff wall and horizontal drain, and evaluate the effect of the cutoff wall's position and depth on seepage under the dam. Salmasi *et al.* (2020) used SEEP/W to simulate the seepage and hydraulic gradients of a zoned earth dam with a clay core, evaluating and comparing vertical and upstream-inclined cores. Sazzad & Alam (2021) implemented SEEP/W to study the seepage characteristics of different types of earth dam – that is, homogeneous, zoned with transition filter, and diaphragm with impervious rock core.

In this study, SEEP2D and SEEP/W finite element numerical models were used to simulate seepage through earth dams with internal cores. They were used first to study the cases of no, vertical, and wedge-shaped cores, and the results were compared to evaluate the models' performance. SEEP/W was then used to complete the study due to its observed accuracy in locating the phreatic line within the core. The effect of the core's characteristics on seepage discharge rate and the phreatic line's levels (taking the dam's base as the datum) at the core's upstream and downstream faces was investigated using different core shapes and geometric characteristic scenarios. Four core geometries were studied – vertical, wedge-shaped (vertical core with sloping sides and a narrower upper width that increases gradually towards the dam's base), upstream inclined (located on the upstream side and thus inclined downstream), and downstream inclined (on the downstream side and inclined upstream). In each case, the effect of changing the core's upper width and slope was also evaluated. The literature review shows that the performances of SEEP2D and SEEP/W were not compared in previous studies, and neither was the downstream inclined core. Design equations are also provided.

## DIMENSIONAL ANALYSIS

## DAM DIMENSIONS

The dam design studied is formed on an impervious foundation to minimize seepage beneath it. The dam's total height is 40 m and its crest width (B) is 8 m. The upstream head (H) is 37 m, with a 3 m freeboard, and the downstream side is dry. The upstream side slope is 3:1, and the downstream slope 2.5:1. Figure 2 shows the dam cross section and dimensions. The crest width and minimum freeboard were selected following U.S. Bureau of Reclamation recommendations (Garg 2005; USBR 2015), while the upstream and downstream slopes were selected on the basis of Terzaghi's recommended values for earth dam side slopes (Garg 2005).

## NUMERICAL MODELING METHODS

### SEEP2D

_{p}). In order to model flow in the unsaturated zone (negative h

_{p}), the hydraulic conductivity can be expressed as the product of the saturated hydraulic conductivity

*k*and the relative hydraulic conductivity

_{s}*k*(which is less than unity and depends on the soil type) – Equation (6):

_{r}Two methods are available to determine *k _{r}*, the linear frontal function method and the Van Genuchten model method (GMS 2008). Trial runs were performed using the two methods and the results compared. The difference was very slight, so the study was done using the linear frontal function method. Also, the amount of unsaturated flow is much smaller than the saturated flow and, hence, has minimal effects on the total seepage discharge results.

In SEEP2D, the problem's geometry is drawn, the boundary conditions are set, and the material or soil properties – hydraulic conductivity and the unsaturated zone parameters – are established for all elements. Finally, the finite element mesh is generated and solution determination starts. The results are displayed as equipotential head contours, pore pressure contours, velocity vectors, and flow lines (GMS 2008).

In this study, the total head boundary condition was used at the upstream dam face with a value of 37 m, while the downstream face was taken as the exit with zero head. Mesh cells of different sizes were then investigated until one was found that yielded accurate results within a suitable calculation time. Good results were obtained using a mesh cell size ranging from 2 to 1 m, which, accordingly, was used in the study. The soil of the dam's shell and core are taken as isotropic, k_{y} = k_{x}, and the linear frontal function method was used to model the saturated/unsaturated flow.

SEEP2D has previously been verified by Ismaeel & Noori (2011) using experimental and field data, and Refaiy *et al.* (2021) using experimental data, and shown to be reliable and offer acceptable accuracy in modeling earth dams.

### SEEP/W

The hydraulic conductivity in Equation (7) is kept constant for each saturated soil zone, but, when unsaturated soil is modeled, k takes different values that depend on the soil's water content (GEO-SLOPE 2012).

The first step with SEEP/W is to define the units and scale, and draw a system cross-section. The soil regions are then defined by connecting the points that define areas with different material properties and the materials assigned to the regions. The boundary conditions can then be assigned. In earth dam simulations, a zero-pressure boundary condition should be assigned to the downstream dam toe, a potential seepage face boundary condition selected for the downstream face, and the reservoir total head boundary condition applied to the upstream face.

The saturated/unsaturated material model was used for the core and shell materials in this study. The sample functions method with saturated water content 0.35 m^{3}/m^{3} was used for the shell (sand) and 0.5 m^{3}/m^{3} for the core (clay) for the volumetric water content data point function. For the hydraulic conductivity function, the Van Genuchten method was used to model the unsaturated zone. The shell and core soils were taken as isotropic, k_{y} = k_{x.} The total head boundary condition was used on the upstream face (H = 37 m), while the potential seepage face boundary condition was used in relation to the downstream face, and the zero pressure boundary condition at the downstream toe.

The results from SEEP/W are plotted as velocity vectors that show the flow direction, the location of the zero-pressure contours and total head contours.

SEEP/W has been verified by Kachare & Jagtap (2017) using field data, and Salem *et al.* (2019) using experimental data, and was found to be accurate and give good results in modeling earth dams.

### Numerical modeling scenarios

Initially, a coreless dam was studied, then the four different core shape and position scenarios (noted above – Figure 2). For each scenario, four different core upper width ratios (b/B) were investigated (b/B = 0.25, 0.5, 0.75, 1), with five core hydraulic conductivity ratios in each case (K_{core}/K_{shell} = 1, 0.1, 0.01, 0.001, 0.0001). For wedge-shaped and inclined cores, three slopes were evaluated (S (H:V) = 0.5:1, 1:1, 1.5:1).

Initially, both SEEP2D and SEEP/W were used. Their results were compared to evaluate their performance with respect to scenarios involving no core, vertical core, and wedge-shaped core (S = 1:1) 41 runs were done with each model for the comparison, then SEEP/W was used to complete the study (160 runs). The scenario schemes are shown in Figure 3.

## RESULTS AND DISCUSSION

### Comparison of SEEP2D and SEEP/W

The modeling results were compared for three cases no core, vertical core, and wedge-shaped core (S = 1:1) to check whether they agreed and decide which was best to use in the study, using R^{2}. R^{2} is a statistical measure of the closeness of the data fit to a regression line and ranges from 0 to 1, with the latter representing a perfect fit relationship (Glantz & Slinker 1990). The model that provides a more reasonable shape for the phreatic line was also evaluated, by visual inspection.

The results from the comparison show that the seepage discharge ratio (q/K_{shell}H) values from SEEP/W match but are slightly higher than those from SEEP2D. SEEP/W calculates the discharge more accurately because its inputs include more soil property details. R^{2} in this case – matching between the SEEP/W and SEEP2D models (y = x) – ranged from 0.9937 to 0.9997, indicating extremely good agreement between them in calculating (q/K_{shell}H), (y_{1}/H) and (y_{2}/H) – see Figure 4.

Figures 5 and 6 show samples of the comparison between the phreatic lines obtained from SEEP2D and SEEP/W for various cases. The phreatic line given by the two models almost coincided for K_{core}/K_{shell} ≥0.01 but, after that, the phreatic line from SEEP2D became distorted and inaccurate within the core zone, because SEEP2D reaches its limit at K_{core}/K_{shell} = 0.01. The phreatic line obtained from SEEP/W was more reliable for K_{core}/K_{shell} <0.01, as its deflection at the soil interface matches the shape and flow line angle of deflection reported by the US Army Corps of Engineers (1993) and Arora (2004). Accordingly, SEEP/W was selected to complete the study, and evaluate more core shapes and characteristics.

### Effect of core shape on the phreatic level and seepage discharge

To evaluate the effect of core shape on the phreatic line's level and seepage discharge, four different core geometries were studied with the different core hydraulic conductivities, as noted above. The upper core width ratio, b/B, and core slope, S, were kept constant (b/B = 1 and S = 1:1) so that only the core shape would affect the results. Figure 7 shows a sample of the model output for the four core shapes (K_{core}/K_{shell} = 0.01).

Figure 8 shows the effect of core shape on the phreatic line level at the core's upstream (y_{1}/H) and downstream (y_{2}/H) faces, and the seepage discharge ratio, q/K_{shell}H, for different K_{core}/K_{shell} ratios, when b/B = 1 and S = 1:1. As can be seen, for K_{core}/K_{shell} from 1 to 0.001, the downstream inclined core gives the lowest phreatic line level at the core's upstream face (y_{1}/H). This may be due to the longer path taken by the water to reach the core's upstream face. The longer path causes more head loss before the water reaches the core and, consequently, a lower value for y_{1}/H.

The wedge-shaped core gave the lowest phreatic line levels at the core's downstream face, y_{2}/H, as well as the lowest seepage discharge ratio (q/K_{shell}H). This can be attributed to the greater width of this core shape, compared to the other cases studied, leading to more flow resistance and head loss.

As the core's hydraulic conductivity ratio approaches K_{core}/K_{shell} = 0.001, the seepage discharge is almost blocked and the head loss is high, which makes the effect of all core shapes almost the same – Figure 8.

### Effect of core slope on the phreatic level and seepage discharge

Figure 9 shows that the core's wedge slope (S) has very little effect on y_{1}/H and q/K_{shell}H, especially when the core's hydraulic conductivity ratio is below K_{core}/K_{shell} = 0.01. The maximum increase in y_{1}/H arising from changing the slope is approximately 4%, while the maximum reduction in q/K_{shell}H is about 16%. The lowest seepage discharge ratio (q/K_{shell}H) and phreatic line level at the core's downstream face, y_{2}/H, occur at S = 1.5:1 for the wedge-shaped core, which is attributed to its greater width.

For the two forms of inclined core, the slope of S = 0.5:1 gives the lowest seepage discharge ratio (q/K_{shell}H) and phreatic line level at the core's downstream face (y_{2}/H). The phreatic line level at the core's downstream face (y_{2}/H) is the parameter most affected by the core's slope, especially for the two types of inclined core, as changing the slope from 1.5:1 to 0.5:1 reduced y_{2}/H by up to 30% for the upstream inclined core and 45% for the downstream inclined one.

### Effect of core upper width and hydraulic conductivity on the phreatic level and seepage discharge

Figures 10–13 show the effect of the core's upper width (b/B) and hydraulic conductivity (K_{core}/K_{shell}) ratios on the seepage discharge (q/K_{shell}H) and phreatic level at the core's upstream (y_{1}/H) and downstream (y_{2}/H) faces. As can be seen, y_{1}/H increases with decreasing K_{core}/K_{shell} ratio in all cases, as reducing the core's hydraulic conductivity inhibits flow, causing the phreatic level to rise at the upstream face.

Figure 10 shows that, for the vertical core, the seepage discharge ratio, q/K_{shell}H, through the dam decreases as the core's hydraulic conductivity ratio, K_{core}/K_{shell}, decreases and as its upper width ratio, b/B, increases. Increasing b/B and reducing K_{core}/K_{shell}, also raises the phreatic level at the core's upstream face (y_{1}/H) while lowering the phreatic level on the downstream face (y_{2}/H). This applies for all K_{core}/K_{shell} and b/B ratios studied, except for K_{core}/K_{shell} = 1, where the core and shell are of the same material, so width has no effect, and K_{core}/K_{shell} = 0.0001, where the core's hydraulic conductivity is very low and seepage discharge, q/K_{shell}H, is almost blocked.

Figure 11 shows that, for wedge-shaped cores, lowering the hydraulic conductivity ratio, K_{core}/K_{shell,} reduces the seepage discharge ratio, q/K_{shell}H, and the phreatic level at the core's downstream face (y_{2}/H), and raises the phreatic level at the upstream face (y_{1}/H). It is also clear that increasing the core's upper width ratio (b/B) has no effect on q/K_{shell}H or y_{2}/H, but causes a very slight increase in y_{1}/H. This may be because changing the wedge-shaped core's upper width (b/B) has minimal effect on the value of its bigger lower width, so its effect can be neglected.

Figure 12 shows that, for the upstream inclined core, q/K_{shell}H and y_{2}/H fall when the core's hydraulic conductivity ratio (K_{core}/K_{shell}) is reduced, and its upper width ratio (b/B) is increased. Lowering K_{core}/K_{shell} and increasing b/B also raise y_{1}/H. As above and for the same reasons, this applies for all the K_{core}/K_{shell} and b/B ratios studied, except for K_{core}/K_{shell} = 1.

Figure 13 shows that, for the downstream inclined core, reducing the core's hydraulic conductivity ratio (K_{core}/K_{shell}) causes q/K_{shell}H to fall. Increasing the core's upper width ratio (b/B) for the same K_{core}/K_{shell} also causes q/K_{shell}H to fall for all hydraulic conductivity ratios except K_{core}/K_{shell} = 1 and 0.0001. Furthermore, reducing K_{core}/K_{shell}, reduces y_{2}/H, except for K_{core}/K_{shell} = 0.1, which yielded the same result as K_{core}/K_{shell} = 1. Reducing K_{core}/K_{shell} and increasing b/B also raise y_{1}/H. The core's upper width ratio (b/B) has a very slight effect on y_{2}/H in this case and it can be neglected especially for b/B >0.5.

### Design equations for dam cores

_{1}/H, y

_{2}/H, and q/K

_{shell}H, for K

_{core}/K

_{shell}ranging from 0.1 to 0.0001. For the vertical cored earth dam, equations (8), (9), and (10) can be used to calculate y

_{1}/H, y

_{2}/H, and q/K

_{shell}H. R

^{2}for the equations is 0.763 for y

_{1}/H, 0.969 for y

_{2}/H, and 0.97 for q/K

_{shell}H.

## CONCLUSIONS

SEEP2D and SEEP/W were used to simulate seepage through cored earth dams. The performance of the two models was evaluated for cases involving no, vertical, and wedge-shaped cores. The models showed excellent agreement in calculating the seepage discharge ratio (q/K_{shell}H), and the phreatic levels at the core's upstream and downstream faces, y_{1}/H and y_{2}/H, respectively. The shape of the phreatic surface within the core given by SEEP/W was more reliable for K_{core}/K_{shell} <0.01, so SEEP/W was used to complete the rest of the study.

Of the four core forms studied, the wedge-shape is the most effective in reducing q/K_{shell}H and y_{2}/H, because it has a large total width compared to the others. The vertical core comes next, followed by the upstream and downstream inclined cores, which have almost the same effect. The downstream inclined core gives the lowest y_{1}/H values because the seeping water follows the longest path to reach the core's upstream face. When K_{core}/K_{shell} is below 0.001, the core almost blocks the seepage discharge and causes substantial head loss, and all core shapes have almost the same effect.

The core slope S has a very slight effect on y_{1}/H and q/K_{shell}H, especially when the core's hydraulic conductivity ratio, K_{core}/K_{shell}, is below 0.01. For the wedge-shaped core, the lowest values of q/K_{shell}H and y_{2}/H occur at S = 1.5:1, due to its substantial total width especially near the dam's base. For both inclined core types, the slope S = 0.5:1 yields the lowest values of q/K_{shell}H and y_{2}/H. The phreatic level on the core's downstream face, y_{2}/H, is the parameter most affected by core slope, especially for the inclined cores, with maximum percentages of 30 and 45% for the upstream and downstream inclined cores, respectively.

The core's upper width ratio, b/B, is effective for all core types except the wedge-shaped core, which shows no response to changes in the upper width for the same hydraulic conductivity because changing b/B has minimal effect on the value of its bigger lower widths. The core's hydraulic conductivity ratio, K_{core}/K_{shell}, is also very important. The phreatic level at the core's upstream face, y_{1}/H, increases as K_{core}/K_{shell} decreases for all core types studied due to their resistance to the flow. Reducing K_{core}/K_{shell} also reduces the seepage discharge ratio, q/K_{shell}H, and phreatic level on the downstream core face, y_{2}/H, reducing the height of the downstream face's wetted portion and, consequently, increasing the safety against seepage effects.

These results will help designers determine the most suitable core shape, slope, and upper width, according to the dam's purpose and design requirements. The suggested design equations can help to predict the phreatic line levels and the amount of seepage discharge accurately, for earth dams with internal cores.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.