CALCULATING EARTH DAM SEEPAGE USING HYDRUS SOFTWARE APPLICATIONS

This paper presents simulations of water seepage within and under the embankment dam of Lake Kowalskie reservoir. The aim of the study was to compare seepage calculation results obtained using analytical and numerical methods. In April 1985, after the first filling of the reservoir to normal storage levels, water leaks was observed at the base of the escarpment, on the air side of the dam. In order to control seepage flow, drainage was performed and additional piezometers installed. To explain the causes of increased pressure in the aquifer under the dam in May 1985 a simplified calculation of filtration was performed. Now, on the basis of archived data from the Department of Hydraulic and Sanitary Engineering using 3D HYDRUS STANDARD software, the conditions of seepage under the dam have been recreated and re-calculated. Piezometric pressure was investigated in three variants of drainage, including drainage before and after modernization.


INTRODUCTION
The use of computer software in the modelling of physical phenomena has become widespread.However, it is important that these tools should be properly calibrated and validated [Simunek et al. 1999, Li et al. 2015], preferably on the basis of values measured in the field [Okamoto et al. 2015] rather than estimates obtained according to some theoretical function.In this context, HYDRUS software, is for example commonly used in the numerical modelling of irrigation systems [Martello et al. 2015, Han et al. 2015], or in the analysis of salt transport dynamics in soil and uptake by roots in different simulation scenarios [Li et al. 2015].
The HYDRUS 3D STANDARD program is a computer package that allows to simulate the flow of water, pollution and heat through saturated and unsaturated porous media, and under steady-state or variable conditions.One of the most common causes of earth dams embankment damages is the phenomenon of internal erosion (suffusion or increase breakdown pressure) that occurs most frequently as a consequence of the water level raising during the transition of the wave [Harris 2015, Vogel et al. 2015, Wosiewicz et al. 2015].
In this study, attempts were made to model the failure of a reservoir earth dam due to seepage at the base of the air-side escarpment.In this paper, the use of archival hydrological documentation and studies has been used.Based on the analytical model, developed in May 1985, a two-dimensional numerical model was built in the HYDRUS 3D Standard program.Numerical simulations were subsequently verified by field measurements.

RESEARCH SITE
The Lake Kowalskie Reservoir was formed after damming the waters of the Główna River (Fig. 1), which is a right tributary of the Warta.The catchment area of the river leading to the dam profile is 189.35 km 2 , while the total catchment area of the river is 235.81 km 2 and the length is 43 km [Kondracki 2002].Lake Kowalskie is characterized by a varied shape and a well-developed coastline, measuring 19.1 km.The length of the reservoir is 7.10 km, its average width is 0.27 km, while its depth varies from 1.5 m to 6.5 m at the dam [Hydroprojekt 2004].It is a multipurpose reservoir, whose most important tasks are to equalize the main river flows, provide water retention for agriculture, and to reduce flood risk.The reservoir measures 162.8 hectares and has a 5.99 million m 3 capacity.

PHENOMENON OF SEEPAGE DURING THE FIRST FILLING
Lake Kowalskie is formed within the valley of the Główna River earthen dam.It is homogeneous, in the form of embankments made of fine and medium-sized sand (Fig. 2).The width of the dam at the crown is 8.50 m, the slope 1:3.The embankment on its water side is strengthened upstream with a screen made of concrete slabs of 0.15 m thickness, placed on a layer of 0.10 m lean concrete.There is no vertical core inside the dam body.The elevation of the dam crest were: 89.00 m above sea level, maximum water level (MWL) 87.30 m, normal water level (NWL) 87.00 m, minimum water level 83.00 m and the level of the bottom of the watercourse below the dam are 80.50 m above sea level.The air-side escarpment is cut with a bench to a width of 3.00 m at 85.50 m above sea level, under which a drainage layer is installed.Moreover, the drainage layer is made of stoneware perforated pipes, placed in a gravel-stone layer.Control wells of φ = 1.00 m are located along the drainage and the drainage outlets are located 81.50 m a.s.l.
The base of the dam front is built of loamy sand originating from Baltic glaciations.These loamy sands are characterized by multiple and irregular inter-beddings and sandy veins.Organic soils occurring in the immediate vicinity of the Główna River have been designed for total removal.
Inter-beddings of sand and gravel at the base of the dam, not separated by a vertical screen from the reservoir, were the cause of the unfavourable effects of seepage at the first filling.In April 1985, after damming up to the elevation of 86.50 m, seepage was found on the air side of the slope, at the front of the dam, on the left bank over a length of about 5 m.The outflow from the earth dam took place therefore at the ground surface between the embankment base and into the band-shaped ditch localized parallel to the dam.Also quite intense seepage occurred in the ditch effluent.Consequently, at the bottom of the ditch and on both of its slopes, an outflow of water containing soil particles was observed.
In the summer of 1985 symptoms of soil particle transportation in the form of very small craters was also seen in the surface below the dam embankment [Lewandowski and Rembeza 1987].The geotechnical investigations completed in May 1986 during the installation of the additional piezometers near the observed water outflow on the surface slope fully confirmed the assumption that the observed seepages resulted from water filtration through the ground under the dam.These studies showed that a continuous aquifer with a free surface occurred above low-permeable loams.The aquifer layer contained saturated medium and fine native sands and a sand-gravel mix.
In order to collect and discharge seepage water, in June 1986 a ditch was made along the air-side slope of the dam on the left bank.To reach the aquifer layer built of gravel sand a relatively deep trench was made, which cutting the low permeability layer.The researchers, Lewandowski and Rembeza [1987], also attempted to conduct seepage calculations to explain the occurrence of excessive pressure in the aquifer layer.However, due to the lack of sufficient data describing the geological composition of the base, results of calculations were approximated by the above-mentioned authors.
The seepage area was divided by a vertical line into two parts, each of a different thickness.Assuming that the dividing line is a line of equal potential, it was assumed that the thickness of the aquifer increased from 0.30-0.50m at the drainage ditch to 4.0 m below the dam bench on the air-side escarpment which further extends under the entire body of the dam.It was also assumed that the low-permeability clay covering the aquifer under the whole body of the dam is impermeable.Calculations were made for the two levels of water in the reservoir.Despite a significant simplification of the geological structure of the substrate, the results were consistent with field measurement [Lewandowski andRembeza 1987, Lewandowski et al. 1990].

EXPERIMENTAL PROCEDURES
Analytical solution Lewandowski and Rembeza [1987] calculated seepage under the dam to clarify the occurrence of excessive pressure in the aquifer (indicated as O).However, as the authors mentioned, in the absence of sufficient data describing the geological conditions, the calculations are approximate and the locations of layers are shown only schematically.The seepage area was divided identically as mentioned above, assuming that the dividing line is a line of equal potential (I and II, Fig. 2).The thickness of the aquifer layer in Part I (D 1 ), at the drainage ditch, was initially assumed to be 0.3 m.In Part II -under the air-side escarpment bench, and under the whole body of the dam -the thickness of the aquifer layer (D) was assumed to be 4.0 m.The capacity of water q 1 flowing in Part I of the thickness D 1 was calculated from the relationship: where: -width of the ditch, m, D -thickness under the dam, m, D 1 -thickness under the ditch, m.Symbols related to the equations presented above are shown in Fig. 2. Capacity q 2 in the second part of the area was determined from the relationship: (2) Lewandowski and Rembeza [1987] assumed that under steady-state flow conditions, capacities in both parts of the filtration area equal (q 1 = q 2 ) and the elevations of the piezometric pressure line h can be determined from the relationship: where: M 1 and M 2 -the denominators of fractions in equations ( 1) and ( 2).
The calculations of seepage under the dam were carried out at two levels of filling of the reservoir: the maximum level of water: 87.30 m a.s.l.(H = 5.0 m) and partial filling: 86.24 m a.s.l.(H = 3.94 m).For the calculations it was assumed that the path of filtration L 1 is variable.The other values were assumed to be: L = 48.0m, D = 4.0 m, D 1 = 0.3 m and B = 0.8 m.The calculated values for the elevation of piezometric head h in the control piezometer at the border of zones I and II are summarized in Table 1.The value of Δh represents the difference between the elevation of the piezometric head and the terrain elevation at the border of zones I and II.The calculation results show the strong influence of the distance L 1 on the location of the piezometric pressure.The change of L 1 from 1.5 m to 3.0 m caused an increase in the pressure level of 0.60 ÷ 0.75 m.This was caused by the small thicknesses of the aquifer level D 1 = 0.30 m.The assessment of the impact of D 1 thickness on the elevation of the piezometric head line h was performed using formulas (1) and ( 2), after taking into account the fixed value L 1 = 5.0 m.The calculation results are presented in Table 2. Based on the results presented in Table 2 it can be seen that the change in thickness D 1 of 0.3 m to 0.9 m causes a reduction of the pressure line of 0.94 ÷ 1.20 m, depending on the height of the water level H.
The results of calculations were verified by measurements of ground water level in the open piezometer located at the assumed separation line of the two areas of filtration.Depending on the height H, the measured pressure levels were equal to 1.85 m at the maximum water level (H = 5.0 m) and 1.45 m at the reduced water level in the reservoir (H = 3.94 m).Lewandowski and Rembeza [1987] emphasized that full compliance of measured and calculated elevation of the piezometric head was obtained with the following assumptions: a) thickness D 1 = 0.3 m and length of the distance of filtration L 1 = 2.30 m, b) thickness D 1 = 0.7 m and length L 1 = 5.0 m.

MAPPING OF NATURAL CONDITIONS IN HYDRUS 3D
Using the archival data, the authors of this article decided to simulate the seepage phenomena occurring under the dam at Lake Kowalskie Reservoir, using the program HYDRUS 3D STANDARD.It was assumed that the simulation would be carried out for three cases at a water level in the reservoir H = 5.0 m: a) simulation of seepage under the dam according to assumptions made by Rembeza and Lewandowski [1987] -identification of the model parameters (before modernization); b) calculation of seepage without simplification of the geological structure of the dam substrate (before modernization) with a functioning drainage system; c) calculation of seepage through the dam and ground substrate after drainage system enlargement (after modernization).
Usually simulations of seepage flow on numerical models use proprietary programs [Sroka et al. 2004, Chalfen andCzamara 2007] or commercial software (GEO-SLOPE Office [Kajewski et al. 2006]).The HYDRUS 3D STANDARD software is a computer package that allows the modelling of flow of water, mass and heat in three dimensions in variably saturated media.
The calculation of seepage flow was performed for a two-dimensional model with a mixed flow regime, using the finite element method [Sroka et al. 2004].In the first step the seepage area was defined, and then the discretization mesh was built.Triangular elements were used [Małecki 2006, Hydrus… 2011].The boundary conditions in the conceptual model were presented in three cases: • before modernization (a and b).
The second kind boundary condition q n = 0 (no flow) was assumed on the embankment of the dam, on the base of the aquifer layer, and on the left and right border of the filtration area.Along the air-side embankment from the level of the bench to the end of the ditch a seepage face condition was used.In this case, the length of the seepage face is not known a priori and its length is determined iteratively by the software.In the HYDRUS software it is assumed that the pressure head is constant and equal to zero along the seepage face.The boundary condition of type I (marked in red) defined the constant water level in the reservoir.The initial condition, defined for the entire seepage area, was established using the hydrostatic pressure at the level of -1 m.Flow under unsaturated conditions was calculated in the presented model using the Richards equation: where: θ -volumetric water content, m 3 • m -3 , h -pressure head, m, K -unsaturated hydraulic conductivity, m • d -1 , K ij A -components of a dimensionless anisotropy tensor K A (which reduces to the unit matrix when the medium is isotropic), S -general sink term, m 3 • m -3 • d -1 , accounting for root water uptake t -time, d, x i -spatial coordinate, m.
The unsaturated hydraulic properties, θ (h) and K (h), in (4) are in general highly nonlinear functions of the pressure head [Hydrus… 2011].The HYDRUS software also implements the soil-hydraulic functions of van Genuchten [Simunek 1999], who used [Małecki et al. 2006] the statistical pore-size distribution model of Mualem [1976] to obtain a predictive equation for the unsaturated hydraulic conductivity function in terms of soil water retention parameters.The expressions of van Genuchten [1980] are given by: where: Six independent parameters: θ r , θ s , α, n, K s , and l are presented in Table 3.They were defined using indirect methods using RETC software [Van Genuchten et al. 1991] (a computer program used to analyze soil water retention and hydraulic conductivity functions of unsaturated soils).The parameters were obtained based on a percentage of the size fractions and bulk density.These hydraulic properties are key parameters in any quantitative description of water flow into and through the unsaturated soil zone.Generally it is important to estimate appropriate hydraulic conductivity using for example the field test method appropriate for different soils, as reported by Nieć and Spychała [2014].Using an inappropriate method for estimation of hydraulic conductivity can give a value that differs by as much as one order of magnitude.
Layers were simulated based on hydro-geological drillings and, in addition, tight (impermeable) reinforcement of the upstream concrete membrane.Calculations were performed with a varying time step in the range of 1 • 10 -5 to 5 days.The model was run for 100 days.

Laboratory measurements
For the analysed computational schemes and boundary conditions the analysis of seepage under the dam was performed.In Fig. 3 the velocity distributions of the three analysed simulation are shown.Calculations show that for a period longer than 20 days the flow conditions become steady.The highest velocity values were observed around the drain and the base of the drainage trench.
Hydrostatic pressure distributions for t = 20 days reflecting steady flow conditions for different simulation scenarios are presented in Fig. 5. Pressure head indicated the appearance of seepage on the embankment at the base of the bench in the case before modernization (a).
The obtained numerical solution (case a) was comparable to the analytical solution (case 0) [Lewandowski and Rembeza 1987] and with the measurements taken after the deepening of the ditch (after modernization).
In sections 1-2 (Fig. 2) leaks appeared, both in the case with simplified geometry and that with a naturally occurring fault (hydrostatic pressure equal to 0).In the case after modernization a significant reduction of hydrostatic pressure observed -below -0.1 m (Fig. 6).The graphs in Figures 6 and 7 show hydrostatic pressure on the boundary line where x = 0 was assumed for the start of sections 1-2 at point 1, and the end of the x-axis was at the end of the ditch cross section -at point 2 (Fig. 2).Subsequently the flow velocity at the border of sections 1-2 (with the appearance of seepage) was analysed (Fig. 7).In the region of seepage a significant flow rate was observed, measuring about 2.25 m • d -1 and 0.33 m • d -1 , 0.07 m • d -1 , for variants a, b and c, respectively.The modernizations also affected the velocity of water seepage into the ditch.The enlargement of ditch dimensions (variant c) caused the highest seepage water velocity (Fig. 7), while relieving the drain.
The calculations confirmed that due to the performance of drainage before the geometric changes of the substrate, the seepage velocity in the drainage did not change.The outflow to drain velocities are practically equal in cases a and b.Simplification of the geometry of the substrate did not significantly affect the velocities.However, the velocities decreased in the case of modernization related to the deepening and widening of the banding ditch.
Hydrostatic pressure at observation point O (in Fig. 2) is shown in Fig. 8.The red line marks the ground level.Simulations confirmed that in case a the water level was above the ground level and different than after modernizations in variants b and c when the water level declined.

CONCLUSIONS
Numerical analysis enabled the quantitative and qualitative description of seepage under the earth dam in different cases.Velocities and pressure distributions near the seepage region, drainage and banding ditch have all subsequently been analysed in detail.Based on the performed calculations, the following conclusions can be drawn: • There is satisfactory convergence of the results obtained by the analytical method, the numerical solution (using HYDRUS 3D software) and the results of field studies.• The analytical methods require significant simplifications related to at least the geometry of the analysed domain.The calculations are laborious and time-consuming compared to numerical methods.• The presented software tool is very useful in the design of earth dams.Numerical analysis software enable the description of unsteady seepage flow with time varying