Evaluation of Earth Dam Leakage Considering the Uncertainty in Soil ‎Hydraulic Parameters

Analysis of earth dams is generally conducted in three stages of stability, deformability and water penetration analysis. Lack of sufficient attention to leakage, as one of the most important issues, leads to erosion and destruction of slope stability. The aim of current paper is to analyze the earth dam leakage with respect to the existing uncertainty in soil hydraulic parameters. In this research, the Monte Carlo (MC) method was used to calculate soil hydraulic parameters. Using these parameters, analysis of Alborz earth dam leakage by means of SEEP/W model based on the finite elements method was investigated. Due to the hydraulic conditions of the core soil, the total head value, pore water pressure, and water flux in core region will change. The results indicate that uncertainty in the hydraulic parameters of Alborz earth dam are significant, thus risk is important in this dam. The application of the proposed methodology in estimation of leakage from Alborz earth dam in Mazandaran province reveals its efficiency and proper accuracy in predicting the amount of leakage flow in earth dams with respect to the possible changes in the hydraulic parameters of the soil. Moreover, it was found that the quantity of seepage increases considerably when the dam is without core, therefore, the core is very necessary to decrease the value of seepage through the earth dam.


Introduction
By development of urbanization and human activities, drought has increased. About 13% of the world's population is currently deprived from drinking water. Water shortages will increase at current pace of water demand. The earth dams are one of the most important and crucial structures for the storage of water and are used in water resources management [1]. Earth dams are of natural materials without mortar. Therefore, the dam body is not integrated due to the bending and tensile forces and all forces are tolerated by the weight of the dam body and with the help of shear strength of the earth materials. The high cost of construction and the severity of damages caused by the failure of earth dams and uncertainties in predicting the behavior of geotechnical structures in terms of the nature of the earth materials, highlight the needs for care and measurement of earth dams. Therefore, the analysis, calculation and design of earth dams should be performed precisely. Analyses of earth dams are generally conducted in three orders of stability, deformability and water penetration analysis [2]. Nowadays, most of the dams worldwide are facing the problem of leakage. This phenomenon threatens their storage function and sometimes causes unpredicted disasters. In unstable areas, a small penetration may grow and will eventually cause the erosion and often the overall destruction of the whole structure [3]. In case of the avoidance of some issues such as the sudden drop of water levels of the reservoir due to releasing and watering and also the uncertainty of soil hydraulic parameters, lead to sudden changes of the boundary conditions of earth dams in terms of pore water pressure and degree of saturation and hence enable instability of the earth dams [3]. The studies conducted on numerous destructed earth dams indicated that around 30 % of these destructions have been caused by the uncontrolled leakage from dam body and foundation [1].
Since leakage in earth dams is one of the most important issues that the lack of attention to it leads to erosion and loss of girder stability, it is necessary to calculate the exact amount of leakage from the body and the dam due to the technical and economic issues to prevent possible financial and property damages and life threats. The existence of the leakage in earth dams is inevitable. However, availability of proper conditions for soil erosion, would lead to soil washing in proper points, followed eventually by dam destruction [4]. According to the performed studies, one of the main and important causes of the destruction of earth dams has been overtopping phenomenon. However, the leakage phenomenon is one of the issues that has gained the attention of the designers of earth dams.
In probabilistic analysis, an explicit or implicit relationship is related to the operation of the input parameters and output response, which is sometimes difficult except for simple cases. Under such circumstances, the concept of uncertainty can be considered. Lohr and Mueller [5] used Monte Carlo (MC) simulation to model hydrograph and to simulate the dam reservoir system. They found that by performing probability analysis, the effects of the available uncertainties can be reduced. Jamel [6] studied the quantity of seepage through a homogenous earth dam without filter resting on impervious foundation using SEEP/W (Geo-Studio). Three different downstream and three different upstream slopes, three different heights and three different top width of earth dam were considered in the evaluations. SEEP/W results were verified with artificial neural networks (ANN) and analytical methods. This study finally resulted in developing an empirical equation in order to determine the quantity of seepage. In a study by Peng and Tian [7] the amount of leakage from the foundation of earth dam was estimated using the neural network technique claiming that the best model in terms of accuracy was selected. Also Fredlund and Gitirana [8] applied the combination of probability methods and Monte Carlo simulation to investigate the leakage in unsaturated soils and to analyse the slope stability. Using the proximal alternating direction method, minimized the number of runs of finite elements model. Montazeri Namin and Boroomand [9] presented a new method to numerically solve the Richard's equation for estimating water flow into porous media like dam bodies.
Ebrahimzadeh et al. [10] studied the risk of earth dam overtopping using system dynamics and Monte Carlo simulation in Hajilarchay dam, north-western of Iran. They considered different uncertainties in return periods of flood and wind and used system dynamics to assess the overtopping risk on earth dams. For evaluating risk in their model, the authors used Monte Carlo simulation method with different iterations. Their findings showed that overtopping risk obtained from the Monte Carlo simulation has less simulation running time and system dynamic is a very effective tool for evaluation and calculation of overtopping risk of earth dams. Rohaninejad and Zarghami [11] combined Monte Carlo simulation and finite difference method to develop a model for predicting the behaviour of earth dams under the postwatering conditions of dam reservoir. Their results showed efficiency of the model and its great convergence to the results made by dam's instrumentations. Mouyeaux et al. [12] studied the probabilistic stability of an existing earth fill dam using a Stochastic Finite Element Method (SFEM) and considering the spatial variability of soil properties based on field data. They implemented a probabilistic field data to analyse the stability of an embankment dam. Random variables and random fields showing the variability of dam materials were integrated into an FE model by performing Monte Carlo simulation. They finally characterized the variability of the sliding safety factor by using probabilistic analysis. Aboelela [13] studied different control methods such as flat slope and toe drainage systems, and a catch drain in the tail water were considered to avoid seepage through earth dams. The hydraulic performance of each control method were evaluated using the analytical solutions. This study also conducted on a physical model. Their findings showed that the control measures have great effect on the characteristics of seepage through earth dams based on pervious foundations. Mohamed [14] developed an empirical relationship to estimate the peak outflow discharge of embankment dam breach due to overtopping flow. The author employed the Monte Carlo method to calculate the area under the flow duration curve and development of the new formula. This formula gave more accurate results and reasonable outcomes in comparison to the actual previous case studies and obtained results from previously developed formulae.
The procedure of dam design and construction in civil engineering with respect to uncertainties in the structure is a very complicated task. The uncertainty parameters in structures include the loads and strengths. By applying this analysis, the parameters are used in analysing as random variables. The reliability of a structure is normally measured with its probability of failure (PF). The application of mathematical relations to compute PF seems to be very simple. However, in practical problems, using these relations becomes very complex. To practically use these relations, the practical methods can be used. One of the best techniques that will help in this regard is using Monte Carlo analysis [11].
In this paper, with regard to the uncertainty of soil hydraulic parameters, the amount of leakage in Alborz earth dam located in Mazandaran province of Iran has been investigated using Monte Carlo method.

Scope of Study
The aerial map and geographic location of the Alborz earth dam is illustrated in Figure 1. This earth dam with a height of 78 meters from the foundation, a crest length of 831 meters and a usable capacity of the reservoir of 150 million cubic meters is located over the Babolrood River at the intersection of Gazu and Chakhani branches, about 40 kilometers south of Babol city in Mazandaran province, Iran. The longitude of the dam is 52.80 E and its latitude is 36.23 N. The crest level and the normal level of water in the reservoir of the dam are 127 and 123 meters above sea level, respectively. The characteristics of the materials of this dam is given in Table 1.

Figure 1. Geographic location of Alborz earth dam area (Mazandaran province, Iran)
After data collection for this dam, the uncertain parameters affecting the leakage were specified by Monte Carlo method and the likely incidental uncertain parameters were generated along with the selection of the proper probability distribution function. Then, using the produced parameters, the finite element method based on the Monte Carlo simulation method was engaged and volume of water leaked into the dam due to variations of dam reservoir levels was analyzed. The leakage analysis was performed using Geo-Studio (SEEP/W) software with respect to the uncertainty of soil hydraulic parameters.

Leakage Phenomenon in Earth Dams
Earth dams are very huge artificial structures whose height often reaches more than 300 meters. These structures are usually generated from the compaction of such materials as clay, sand, soil or rock. These include two main parts, a natural cover and an impermeable core. The stability of these structures is created by the interface friction between the constituent materials [5].
The design and construction of these structures should be in a way that there would be no leakage to the underlying water or around the structure. Hydraulic fracture can occur in the fine-grained core of earth dams and water leakage from the dam leads to terrible consequences [3]. Equations governing leakage flow are as follows: , {H} and <N> are gradient matrix, matrix of hydraulic conditions of the element, pressure vector in each node and the vector of interpolation functions, respectively. Additionally, T is the time of passing discharge per width of each element, q is the amount of passing discharge per width of each element, τ is element thickness, A is a function to compute the sum of elemental areas and λ and L are storage term for analysis in unsteady leakage mode and a function to compute the length of each element, respectively. In axially symmetrical state, equation 1 can be rewritten as ( Equation 2): It can be summarized as (Equation 3):

Uncertainty of Soil Hydraulic Parameters
In designing structures, the applied loads, structure geometry, material properties, construction process and the surrounding environment have uncertainty and all of them are coordinated by probability functions. However, today, in engineering world, structural analysis is more based on certain methods. Thus, in designing, based on code methods, all the effective variables in the problem including the geometry, mechanical properties of materials and certain loading are presumed definite and uncertainties are covered by applying safety factors. These factors guarantee the stability of the structures by the way of regulations and codes. In practice, most of the parameters have much intrinsic uncertainties that presuming them as certain; and applying high safety factors to cover uncertainties causes the plan to be non-economical [4].
Despite the existing uncertainties in the problems of structural analysis and the advantages of using the probability analysis, applying such techniques has not gained the attention of the public. Lack of attention to the application of probability methods by engineers including the current training of statistics and possibilities for engineers has only led to the presentation of basic information. Consequently, instead of dealing with possibilities, they apply safety factors in their analysis and there is a wrong public presumption that probability analysis requires a lot of information and time and is more difficult than certain methods. Eventually, limited studies have been practically published on the performance and the advantages of probability analysis method. The above cases have caused the safety factors rather than probability of minimum failure to have much penetration into the field of civil engineering. However, better and more economical designs of certain methods is not sufficient and merit in today's world full of competition. And, today, the necessity of using statistical methods is increasing more and more [15].

Monte Carlo Method
One of the important and practical methods in probability analysis and reliability of the structures is Monte Carlo simulation method. The abundant abilities of the Monte Carlo method have caused its extensive usage not only in structural reliability but also in other engineering and mathematical fields. The abilities of the Monte Carlo method have caused its increasing application in solving various problems of structural reliability as well. These methods that establish third class of structural reliability are an attempt to obtain the best estimation and probability of structural destruction. The high accuracy of the obtained results, capability of error estimation and control, simplicity of the method, capability of precise solving of any problem including the problems of part or the whole system's reliability, time-coordinated problems or independent ones and not having limitation in confronting the complex forms of limit state function are some of the outstanding points of Monte Carlo methods [3].
Monte Carlo simulation is known as a simple random sampling method or statistical testing method, based on the creation of a random sample for variables having uncertainty. Using Monte Carlo method for probability analysis of structural problems is only practical using computers. It is a powerful mathematical tool for determining the approximate probability of a special incident which is the result of a series of random procedures [3]. The only disadvantage of Monte Carlo method may be the high volume of operations. When using this method, achieving a pre-determined accuracy requires a specific number of simulations and it has an inverse relationship with the value of failure probability. Since, normally, the failure probability in structures is about one ten thousandth or less, to solve a simple problem, a great number of simulation and sampling are usually required. In summary, the Monte Carlo method is generally a simple method which could achieve any desirable levels of accuracy. However, high costs of calculations cause the application of ordinary Monte Carlo methods in structure reliability analysis to be impractical and infeasible and it is only used in research problems to control the result and the accuracy of other methods in cases that exact results are not available.

Generation of Random Samples for Variables in Monte Carlo Simulation Method
In this method, to generate random samples from uncertain variables of the problem, a set of random numbers are taken from computer being between zero and one. Each of the variables has a probability distribution in a manner that the vector of its accumulative distribution function has a range between zero and one. Therefore, Monte Carlo with respect to the probability distribution of each variable produces random samples by establishing a proportion between the random numbers selected by computer and accumulative distribution function of a variable. Initially, m random number like z1, z2, …, zm is taken from computer, where there is n number of uncertain variables and as described above, the value of each variable is calculated (X = x1, x2, …, xm). If the variable X has a normal probability distribution with an average of µ and standard deviation of σ, the probability density function is as follows (Equations 5 to 7): If Z = F(x) then given the Equation 7: And if variable x has a log normal probability distribution with an average of ' and standard deviation of σ', then we have: Using the equations of 9 and 10, the unknown values of  and σ are calculated (Equations 11 to 13).
In this research, soil hydraulic parameters are initially generated in MATLAB for simulating and analyzing the phenomenon of Alborz earth dam leakage using Monte Carlo method and was imported into Geo-studio as primary data for simulation (Table 2). Figure 2 shows a summary of the methodology and steps of calculations.

Hydraulic Conductivity
The ability of soil in transferring water under both saturated and non-saturated conditions is expressed by hydraulic conductivity function. As can be seen in Figure 3, in saturated soil, all the pores are filled by water and as soon as a negative pressure higher than AEV is occurred, air goes inside the biggest pores and these pores filled with air will prevent the flow of water inside soil. Therefore, the ability of soil in transforming water is reduced. With an increase in negative pore water pressure, more pores are filled with air and hence, hydraulic conductivity is more reduced. Figure 3 illustrates that the capability of water flow inside a soil profile depends on the amount of the available water in the soil, specified by water content function. In this article SEEP/W software, which uses the method of Ferdland et al. [15] and that of Van Genuchten [16] to estimate hydraulic conductivity function in unsaturated zones, were employed by importing the generated hydraulic parameters of the soil.

Calculation of the problem variables
Generation of the soil hydraulic parameters in MATLAB Application of the data into Geo-Studio software

Software Verification
To ensure the correctness of software outputs, some outputs must be initially verified with analytical methods provided in valid scientific resources. Therefore, in the first step, some calculation methods for water leakage from the earth dam's body is investigated and compared with the outputs of SEEP/W.

Changing the Length of the Drain and the Distance of the Saturation Line from the Downstream Slope
Chahar [17] proposed the following equation (Equation 16) to determine the distance between downstream slope and free flow line using analytical terms for a homogenous earth dam over an impermeable foundation. The parameters used are shown below. In this section, the given numerical model will be assessed using this equation.

Comparing the Distance between Downstream Slope and Free Line of Flow Obtained from Software and Analytical Method
In both analyses, with an increase in drain length, the distance between downstream slope and free flow line is also increased. The results of the calculated distance at both methods are very close; and in all modes, the data gained from the analytical method are higher than the software (Figure 4).

Results
In this research, Monte Carlo method was used to calculate soil hydraulic parameters. Using these parameters, leakage analysis was conducted using Geo-Studio (SEEP/W) software. The results of simulation including the value of the leaked discharge in various likely hydraulic parameters inside the body of Alborz dam have been illustrated. Pore water pressure in each node obtained from steady state analysis was imported into transitory analysis as boundary conditions. Soil properties like permeability and mass moisture defined in steady state analysis was also imported in transitory analysis. A view of 2D simulation of Alborz earth dam is shown in Figure 5.

Figure 5. A view of the simulated model of Alborz dam in Geo-Studio
Meshing the model was performed to conduct analysis and with respect to this, the given model had 832 nodes and 766 elements. The dimension of each element was approximately 5 meters.

Figure 6. The generated overloads in the body of Alborz earth dam
In this study, parameters of soil hydraulic conductivity were considered as uncertainty factors and variations of the water level behind the dam, pore water pressure, flow speed, and water head inside the dam were studied. According to Figure 6, under the effect of water level behind dam, the body of the dam is pressurized. The generated total head in the initial 200 meters of the earth dam reaches over 45 meters and it gradually reduces towards the behind of the dam. In the core of the dam, it varies from 40 to 20 meters and it reaches zero at the end of the dam. This means that the downstream head of water remained zero during the seepage. In this part of the study, pore water pressure inside the body of the earth dam, including the clay core and covers were studied. As shown in the Figure 7a, the pore water pressure inside the body is 45 kPa in the distance of 233 meters. It was constant in clay core. The generated total head up to the initial 280 meters of the earth dam was 45 kPa and is gradually reduced towards the back of the dam and reached less than 5 kPa in 380 meters, which can be due to the reduction of water level.
According to Figure 7b, pressure head inside the body up to water flow level (i.e., the initial 233 meters) was equal to 200 meters. In other words, the amount of the generated pressure at the lowest point of dam, which is in contact with the water behind the dam, has its highest value and reaches its lowest value when reaches the back of the dam. The pressure head in clay core where the water is not flowing was equal to 200 meters and it became negative at the back of the dam.

Figure 8. Contour of a) the water flow speed and b) output gradient in the body of the dam
To investigate the flow speed and gradient, as can be seen in Figure 8a, water flow speed until a distance of approximately 200 meters was less than 0.02 meters per second. Then, this speed was gradually increased until the clay core of the dam and reached 0.08 meters per second in clay core. Whereas, by increase in dam height, the water flow speed became higher which in the height of about 75 meters, the water flow speed reached its maximumhigher than 0.24 meters per second. The water flow along the dam length after passing the core at the distance of 280 meters from the beginning of dam decreased again and reached 0.02 meters per second. Figure 8b illustrates the output gradient from the dam's body. As can be seen, gradient had its lowest value at the beginning of the dam and gradually increased towards the center of clay core. In other words, gradient was less than 0.04 in 140 meters from the dam and reached its maximum value, more than 0.32, in the center of clay core. After passing the core, the output gradient decreased again and reached to less than 0.02 at 430 meter distance from the dam. Thus, the presence of core had an important effect on decreasing the exit gradient.
From these figures, it can be concluded that the core is very necessary in the dam to lower the seepage line and decrease the total head and seepage through the dam. It can be noticed that the seepage line is raised when the core does not exist.
Another important parameter in earth fill dams is the head of water in body of the dam. For better illustration, the amount of investigated parameters are shown in the nodes of the different layers. The amount of total head in the body of the dam is illustrated in Figure 9a. As shown in this figure, in the nodes of the initial layer, it begins from 50 meters and reaches 40 meters. In the layer on the back of the clay core, it starts from 20 meters and reaches 0.
In order to assess the effects of the clay core more precisely, the amount of total head in clay core is demonstrated in Figure 9b. In the nodes of the initial layer, it starts from 40 meters and reaches 32 meters in the center. At the end, it starts from 27 meters and reaches 20 meters. As mentioned before, another important parameter is the pore water pressure in the body and core of the earth dam. The amount of pore water pressure is demonstrated in Figure 10a. These results show that in case of nodes in the initial layer, pore water pressure started from 450 kPa and reached -250 kPa in nodes close to clay core. Having passed the core, it gradually increased from about -500 kPa to zero. Figure 10b depicts the amount of pore water pressure in the nodes of the clay core. As it shown, in the initial nodes, the amount of pore water pressure was around 400 kPa and then decreased to the very low amount of -550 kPa in case of the middle nodes. For the final nodes, it increased again to about 200 kPa. This figure shows the substantial variations of the pore water pressure in the clay core of this earth dam.
The amount of overburden pressure in the nodes of clay core is shown in Figure 11. As can be seen in this figure, in the initial nodes, the value of the overburden pressure was about 40 meters and it decreased to -57 meters in the middle nodes. Then in the final nodes, it increased to 20 meters. By comparing Figures 10b and 11, one can see the accordance of the results of these two figures. In the initial nodes of the first layer of the clay core, water flux slightly increased to 0.3 m 3 /s and then drastically decreased to -1.19 m 3 /s. Finally, it gradually rose to zero in the final nodes of the core.
The amounts of water flow speed in the body and core of the dam is depicted in the Figures 13a & 13b, respectively. In the initial nodes of the primary layer of the dam's body, the flow speed started from zero and reached its maximum at the middle nodes. Then, in the final nodes it declined to zero again. In the initial nodes of the first layer of the clay core, the flow speed started from 0.025 m 3 /s and reached its maximum value (about 0.25 m 3 /s) in the distance of 250 meters. In the final nodes of the core it decreased to about 0.02 m 3 /s. This figure shows that the variations of the water flow speed in the clay core are substantially higher than that in the body (cover) of the earth dam.
The rising of water level at the upstream of an earth dam may cause seepage pressure to downstream and increase pore water pressure causing the soil shear strength to decrease. Whereas the drawdown of upstream water level results in the increase of pore water pressure in the dam body and the seepage pressure to the upstream [18].

Conclusion
Regarding the water crisis in the country, the role of dams in water storage, supplying the required drinking water, and prevention and control of flood is undeniable. On the other hand, given that dams are mainly built near the places with dense population, destruction of a dam is not only economically destructive, but also is a great threat for the people living downstream.

Figure 13. Water flow speed in a) the body of the dam and b) the clay core
In the current design of the dam, there are uncertainties in the construction of the dam that are not considered in the usual design methods. The purpose of this paper was to investigate the earth dam leakage with regard to uncertainty in the hydraulic parameters of the soil.
Investigating and obtaining hydraulic parameters of soil which requires mathematical and probabilistic calculations, using mathematical software which has high ability in this field, and the preliminary data used for simulation was prepared. Monte Carlo simulation with the limit equilibrium approach was the basis for discussing the stability of earth dams under probabilistic conditions. The results of the analysis showed that considering the uncertainty parameters in the construction and design of a dam will increase the safety of the dam against breaking. Due to the hydraulic conditions of the core soil, total head value, pore water pressure, and water flux in the area of the core varied. The results also indicated that uncertainty in the hydraulic parameters of the Alborz Dam is significant and therefore the risk in this dam is really important.
From this study, it can be found that the quantity of seepage increases considerably when the dam is without core, therefore, the core is very necessary to decrease the value of seepage through the earth dam. Generally speaking, the use of unlimited boundary elements allows the element to be expanded at points far away from the point of effect without extending the geometry of the model. It is suggested that in future studies the influence of other distribution variables as input probability variables, and correlation of the parameters in the Monte Carlo analysis could be considered.

Funding
This project was financially supported by Department of Civil Engineering, Tafresh University.

Conflicts of Interest
The authors declare no conflict of interest.