EVALUATION OF SEEPAGE AND STABILITY OF DUHOK DAM

One of the main causes of the earth dam failure is the seepage. This seepage can cause weakening in the earth dam structure, followed by a sudden failure due to piping or sloughing. For this purpose a finite element method through a computer program, named SEEP2D, was used to determine the free surface seepage line, the quantity of seepage through the dam, the pore water pressure distribution, the total head measurements and the effect of anisotropy of the core materials of Duhok zoned earth dam. First, the accuracy of the program was tested via the data of experimental dam and the results showed an acceptable accuracy of the program. The effect of the ratio of the permeability in the horizontal direction to that in the vertical direction (Kx/Ky) on seepage was tested and results indicated an increase in seepage quantity as this ratio increased. The stability of Duhok zoned earth dam was analyzed using a slope stability computer program, named STABIL2.3. The program is verified through a dam example of known factor of safety (solved by hand calculations). The results of the verification indicated a good accuracy of the program. The slope stability analysis results showed that the factor of safety decreases with the increase of Kx/Ky ratio. The analysis of the results of this study showed that Duhok zoned earth dam is safe against the danger of piping and slope sloughing under the present operation levels. Also, the present study showed that the field piezometers readings of the dam are not accurate.


INTRODUCTION
Duhok dam constructed in 1987 is located on Rubar Duhok river 2km to the north of Duhok city center. Both main and coffer dams are earth fill type with central clay core followed by one layer of filter on both sides and supported by sand and gravel shell as shown in Fig. (1). One of the important stages in design of earth and rock-fill dam is the exact evaluation of seepage quantity, pore water pressure distribution used in the analysis of slope stability, total head measurements and hydraulic gradients in various sections of the dam to ensure stability and avoid endangering effects such as piping and slope sloughing [1].
The main objectives of this study are: 1. Theoretical seepage analysis for Duhok zoned earth dam using the finite element method through a computer program named (SEEP2D) at different sections. 2. Comparing the theoretical results with the field piezometer readings. 3. Investigating the effect of anisotropy of core materials. 4. Analyzing the slope stability of different dam sections through a computer program named (STABIL2.3) performing the conventional limit equilibrium method of slices, the simplified Bishop method and Fellenius method.

SEEPAGE FLOW THROUGH EARTH DAMS
Seepage flow of water through porous media depends on the soil media, type of flow, properties of liquid and hydraulic gradient. Seepage piping account for approximately 50% of all earth dam failures [2]. Many methods have been developed to solve seepage problems, these methods can be classified as analytical, experimental and approximate methods [3,4,5,6].
Mishra and Sing [7] analyzed a steady state 2D flow through a homogeneous levee with a horizontal toe drain resting on an impervious base using fragments method. Dupuit, Schaffernak, Pavlovsky, Cassagrande, and Polubarinova-kochina, developed approximate methods to calculate the seepage quantity and the free surface seepage line of earth dams [3,4,5]. Billstein et al. [8] used experimental models to determine discharge, pore water pressure, seepage face and free surface profile. Most of mentioned methods are valid for simple geometries, isotropy and homogeneous soil media, but with complex geometry, anisotropy, non homogeneous and nonlinear material behavior these methods are difficult to be used and may limit the solution accuracy. Therefore, during the last 30 years, considerable advances have been made in the theory and practice with the help of high speed digital computers. As a result, the numerical methods (finite element, finite difference and boundary element) and their applications are widely used for the solution of many engineering problems. The finite element method (FEM) is one of these numerical methods that can be used for accurate solution of many complex engineering problems. Over the years, the FEM seems to be the most powerful and versatile tool over other numerical methods for solving a wide variety of practical problems efficiently [9]. Muhammad [10] used the finite element method to study seepage through earth dams. Abo [11] used the FEM to determine the free surface position, the seepage quantity, the pore water pressure distribution, the total head measurements and the anisotropy of the material in AL-Adheem zoned earth dam. Thieu et al. [12] used 2-D and 3-D finite element method for modeling the seepage in a saturated/unsaturated soil system of the earth dams under steady state and transient seepage conditions. Bardet and Tobita [13] used a finite difference method for the solution of unconfined seepage with an unknown free surface.

SLOPE STABILITY
Every soil mass which has a slope at its end is subjected to shear stresses on internal surfaces or failure plans in the soil mass near the slope. This is due to the gravitational forces that try to pull down parts of the soil mass adjoining the slope. Several models and analytical techniques have been developed to determine the critical slip surface and the associated factor of safety such as method of slices. The factor of safety is: where, FS factor of safety, f failure shear strength of the soil, and shear stress of the soil. Fellenius (1936) (cited in [14]) developed the Ordinary Method of Slices (OMS) known as Felonious or Swedish method. It is assumed that the forces acting on sides of any slice are neglected. Bishop [15] developed a method based on the assumption that the interslice forces are horizontal. This method is called the Simplified Bishop Method (SBM) in which a circle slip surface is assumed. Janbu (1956) (cited in [16]) developed a simplified method that assumes zero interstice shear forces. Morgenstern and Price (1965) (cited in [17]) developed a general analysis in which all boundary and equilibrium conditions were satisfied and in which the failure surface might be of any shape, circular, non-circular or compound. Wright et al. [18] used the finite element method to determine the factor of safety and associated critical slip surface and compared with the simplified Bishop method. Martins et al. [19] used the rigid-plastic methods (Felonious and Bishop) and the elasto-plastic finite element method (FEM) to analyze the slope stability of embankment with and without surcharge. Hammah et al. [20] compared the finite element method with the limit equilibrium methods for slope stability analysis.

SEEPAGE THEORY
The governing partial differential flow equation for steady state condition two dimensional flow through anisotropic, homogeneous porous media is presented by Laplace equation as: where, y K , x K coefficient of permeability in (x,y) directions, H total head of water and equal to z p w , p pore water pressure, w unit weight of water, and

FINITE ELEMENT METHOD
The basic concept of the finite element method is to divide the problem region into subdomains (finite elements) connected at their common nodal points and that the unknown function of the field variable is defined approximately within each element. The approximate solution of each element expressed by continous function is as follows: where, noe total number of elements.
There are more different approaches to formulate the approximate solution of the problems. In the present investigation, the weighted residual method with Galerkin's criterion is used: (12) where, R the element residual, and j W the weighting function.
Various forms of weighting function sets can be used in practice, each leading to a different weighted residual approximation method. One of the best approximations is called the Galerkin's method, in which the weighting function is made equal to the shape function defining the approximation. Applying this principle in finite element method over the whole problem domain, Eq. (12) can be written as: where, e v the domain of element (e).

Al-Rafidain Engineering
Vol. 19 No.1 February 2011 Using weighted residual Galerkin's method, Eq. (17) can be written as: where,  (19) where, e k represents the element stiffness matrix and known as follows: Surface integral e i s in Eq. (20) is unknown. This condition doesn't cause any problem in the solution of the differential equation because it is not related to the nodes that exist on the free surface but to the reservoir boundary.
Eq. (20) doesn't apply in this part because the head is known; therefore, the element matrix can be simplified as: (21) where, i N interpolation or shape function.
The interpolation function for triangular element can be expressed as: Derivations of the above functions are shown in [21].

SLOPE STABILITY ANALYSIS
In the present study, two traditional methods for slope stability analysis of Duhok zoned earth dam are used. These methods are:

Fellenius (Swedish) Method
For steady state seepage condition the factor of safety using Fellenius method is:

Simplified Bishop Method
For steady state seepage condition the factor of safety using Bishop method is:

RESULTS AND DISCUSSIONS
The results obtained by the application of the finite element program, named (SEEP2D), to Duhok zoned earth dam are analyzed and discussed extensively. The results include the location of the free surface seepage line, the quantity of seepage through the dam, the pore water pressure distribution, the total head measurements and the effect of anisotropy of the core material. For the verification of the program, experimental data of homogeneous dam constructed by Kellogg (cited in Abo [11]) where used. Kellogg's dam was built with coarse sand of 38 cm height and 68 cm base length with upstream and downstream slopes of 1H:1.3V. The origin of the x and y coordinates of the points tested was located at upstream slope toe, for more details, see [21]. The results indicated an acceptable accuracy of the program as shown in Table (1). The results of the slope stability analysis are also presented and discussed using a computer program, named (STABIL2.3), which used the conventional slope stability methods of slices (simplified Bishop and Fellenius methods) for determining the minimum values of the factor of safety and the associated critical slip surfaces. The program is verified using three dam examples of known factors of safety (solved by hand calculations), for more details, see [21]. The results of the verification indicated a good accuracy of the program as shown in Table (2). The seepage and slope stability analysis of Duhok zoned earth dam shown in Fig. (1) is studied according to the material properties illustrated in Table (3) [22].
The seepage through Duhok earth dam, which is a non-homogeneous earth dam (zoned earth dam) consisting of the earth fill shell with central clay core is studied. According to "SEEP2D Primer" [23] a portion of downstream is omitted because the ratio of permeability between the shell and the core material is more than one hundred so as to avoid any confusion in the solution.

Application of the (SEEP 2D) Program on Duhok Dam
For the seepage analysis of Duhok zoned earth dam, four sections are studied using a finite element computer program (SEEP2D) they are: sections (P1-P2, P3-P4, P5-P6 and P7-P8) and for different reservoir water levels (613.12, 611.16, 610.56 and 609.96 m, a. s. l.) [21]. These sections are located at piezometer stations shown in Fig. (2). In this paper the results of section P1-P2 are illustrated and discussed. The results of sections (P3-P4, P5-P6 and P7-P8) are shown else where see [21]. Section P1-P2 is located at station (150) m from the beginning of the dam along the crest length. Fig. (3) illustrates the construction of the mesh for this section. Table (4) represents a comparison of the finite element readings with field piezometer readings and showed that the finite element readings are greater than the field readings. This means that the phreatic line determined by the program is higher than the field phreatic line. So as to reach the field phreatic line the Kx/Ky ratio must be changed. For this purpose the Kx/Ky ratios are taken as (0.1, 0.2, 1, 5, and 10). The results indicate that the phreatic line is higher by increasing Kx/Ky ratio and is lower by decreasing it which agrees with the results given by [6,10,11]. When Kx/Ky increases the horizontal velocity is higher than the vertical velocity, therefore, the phreatic line will be more flat. To lower the phreatic line, Kx/Ky ratio should be less than 1, and this means that the coeffecient of permeability in the vertical direction (Ky) must be higher than the coefficient of permeability in the horizontal direction (Kx), but this condition is not correct as the case in the field. Since, in the compaction process the soil is laid in horizontal layers and then compacted, therefore, Kx value must be higher than Ky value. Low field phreatic line may be attributed to one or more of the following reasons: (1) Inaccurate field piezometer readings. (2) Bad installation of the piezometers because they were installed after the dam construction.
(3) Effect of consolidation of material. (4) Cracks or deformations in the dam. However there is no field evidence on that. Seepage quantity of section P1-P2 is shown in Table (5), while the total seepage quantity of the whole dam is shown in Table (6) for Kx/Ky values of (1, 5 and 10).  (4) to (13) show the pore water pressure distribution and total head measurements for Kx/Ky = 0.1, 0.2, 1, 5 and 10 ratios and reservoir water level 613.12m (a. s. l.). Figs (14) and (15) represent the relationship between seepage quantity and Kx/Ky ratio and indicate that the seepage quantity increases with the increase of Kx/Ky ratio.

Application of the (STABIL2.3) Program on Duhok Dam
According to the seepage analysis of the four different cross sections of Duhok zoned earth dam, the slope stability of these sections are also studied using the computer program (STABIL2.3) [21], but in this paper the results of section P1-P2 are presented only.
The downstream slope stability of this section of the dam is studied using the computer program to determine the minimum factor of safety and associated critical slip surface. Two conditions are considered for the slope stability analysis, with and without seepage. In case of seepage consideration the steady state seepage analysis is used. So as to study the effect of the free surface line on stability of slopes, the (0.1, 1 and 10) Kx/Ky ratios are used for reservoir water level 613.12m (a. s. l.). According to the properties of materials of the dam in Table (3), the program is operated. Tables (7) and (8) show the results of stability analysis of this section while, Figs. (16) and (17) show the critical slip surfaces determined. The results indicated that this section is safe against slope failure in both cases (with and without seepage). The results also indicate that the factor of safety decreases with the increase of Kx/Ky ratio, this is due to the increase of pore water pressure and the free surface line will be higher with the increase of the Kx/Ky ratio, this agrees with the results concluded by [24].   1. The finite element program (SEEP2D) can be used to analyze the homogeneous and non-homogeneous (zoned) earth dams. Therefore, the program is applied on four different sections of the Duhok zoned earth dam so as to determine the location of the free surface seepage line, the quantity of seepage through the dam, the pore water pressure distribution and the total head measurements. The results from this program showed acceptable accuracy compared with the experimental results. 2. The Kx/Ky ratio has a significant effect on the location of the free surface line. The free surface is higher by increasing the ratio and it is lower by decreasing the ratio. 3. The quantity of seepage increases as the Kx/Ky ratio increases, therefore, the loss of water through the dam increases. 4. The total seepage quantities of Duhok dam for Kx/Ky ratios (1, 5, 10) are (987.7, 4066.8, 7901.7) m 3 /day, respectively, for reservoir water level 613.12 m (a. s. l.). 5. The Kx/Ky ratio has a direct effect on the stability of slopes. Factor of safety decreases with the increase of Kx/Ky ratio, this is due to the increase of pore water pressure and the high free surface seepage line. 6. The slope stability analysis results obtained from the application of a computer program (STABIL2.3) on Duhok dam showed that the dam is safe against the danger of slope failure due to slippage under the operation heads. 7. Results of the present study showed that field piezometer readings of the dam are not accurate probably due to bad installation of piezometers and personal errors.