Haldia Institute of Technology Publishing International Journal of HIT Transaction on ECCN Int.J.HIT.TRANSC:ECCN.Vol.7: Issue 1A (2021) Page 25-61 ISSN: 0973-6875 Available Online at www.hithaldia.in/locate/ECCN All Rights Reserved ORIGINAL CONTRIBUTION Analysis of Seepage Beneath a Sheet Pile Wall in Homogeneous Soil by Finite Element method Ajit Kumar Paria Department of Civil Engineering, Haldia Institute of Technology, West Bengal, India ABSTRACT An analytical method was followed to find flow characteristics, like flow rate, exit gradient that are to be determined for seepage through a uniform and homogeneous soil. The finite element program, SEEP/W was introduced to analyze flow characteristics of a sheet pile wall a homogeneous soil. Extensive analyses were performed for different conditions, including soil layer thickness, soil hydraulic conductivity and sheet pile wall depth. The flow rate and exit gradient were determined for each condition analyzed. The analysis can help a designer to get solutions to the seepage problem. The results were depicted in simple charts. The method was verified by comparing the results with those obtained from SEEP/W program, and good agreement was noted. KEYWORDS: seepage, sheet pile, homogeneous soil, exit gradient, hydraulic conductivity. 1. INTRODUCTION A sheet pile wall consists of a series of sheet piles driven side by side into the ground, thus forming a continuous vertical wall for the purpose of retaining an earth bank. Sheet pile walls are retaining walls constructed to retain earth, water or any other fill materials. These walls are thinner in section as compared to masonry or concrete retaining walls. Certain configurations of sheet piling are termed as bulkheads or cofferdams. A sheet pile retaining wall of the water front is called the bulkhead. A bulkhead may be constructed in open water and backfilled or it may be constructed in natural ground and then the earth removed from its face by dredging. The former type is known as a fill bulkhead and the latter a dredged bulkhead. A cofferdam is a reasonably water-tight enclosure made of sheet pile walls, usually temporary, built around a working area for the purpose of excluding water during construction. Sheet pile walls are employed as bulkheads in piers, docks and harbours, and in sea walls, break waters, and other shore protection works. Seepage *Corresponding Address: The stability of earth structures and of natural deposits is dependent not only upon the static properties of the soil but also upon the forces produced by water as it seeps or flows through the pores of the soil. As an aid to his judgment in the design of earth structures or the stabilization of earth deposits, the engineer should be able to estimate through analyses the magnitude of seepage forces and pressures, and the quantities of water flowing through the soil. Because the loss of head as water flows through soils is a reaction against the soil grains, and because the strength of cohesionless soils is produced by pressure between the grains, water flowing downward through fine grained soils increases the strength of such soils. When water flows upward through the soil, gravity pressure between the grains is reduced by the seepage forces and the strength of the soil is reduced. 2. OBJECTIVE & SCOPE OF PRESENT STUDY From literature review it is noted that ample scope of works are to be carried out on the study P a g e | 25 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) of ‘seepage analysis beneath sheet pile’. Literature study also presents works related to seepage analysis below the base of earthen dam and rockfill dam have been done quite more. So keeping this in view the objective of the present investigation is to study the seepage analysis beneath the anchored bulkhead as an earth retaining and water front earth retaining structure. Following are the scopes of the work considered in the present study: 1. Design of cantilever earth retaining sheet pile wall in dry condition. 2. Design and seepage analysis below the cantilever earth retaining sheet pile wall considering water table at different levels using SEEP/W software 3. Design and seepage analysis of anchored bulkhead earth retaining sheet pile for different levels of water table using SEEP/W software 3. METHODOLOGY This chapter presents the methods, equations, procedures, and techniques used in the formulation and development of the seepage function. It is of value to be familiar with this information in order to analyze flow problems. An understanding of these concepts will be of great benefit in applying the software SEEP/W, resolving difficulties, and judging the acceptability of the results. 4. FLOW OF FLUID: NAVIER-STOKES EQUATIONS We shall be concerned only with the equations governing the motion of viscous, incompressible fluids. Preserving an analogy with twodimensional solids, (u) and (υ) become velocities in the (x) and (y) directions respectively and (ρ) is the mass density as before. Also, (Fx) and (Fy) are body forces in the appropriate directions. Conservation of mass leads to: ∂ρ ∂ ∂ ---- + ---- (ρu) + ---- (ρv) = 0 ∂t ∂x …………………….(4.1) ∂y But due to incompressibility this may be reduced to: ∂u ∂v ……………………………(4.2) ---- + ---- = 0 ∂x ∂y Conservation of momentum leads to: ∂u ∂u ∂u ∂σx ∂τxy ρ (---- + u ---- + v ----) = Fx + (----- + ------) ∂t ∂x ∂y ∂x ∂y ∂v ∂v ∂v ∂σy ∂τxy ρ (---- + u ---- + v ----) = Fy + (----- + ------) ∂t ∂x ∂y ∂y ……………….(4.3) ……………….(4.4) ∂x where (σx, σy, and τxy) are stress components as defined for solids. ISSN: 0973-6875 P a g e | 26 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Introducing the simplest constitutive parameters (u) (the molecular viscosity), (λ) (taken to be –⅔μ) and (p) (the fluid pressure) the following form of the stress equations is reached: ∂u ∂v ∂u σx = – p + λ (---- + ----) + 2μ ---∂x ∂y ∂x ∂u ∂v ∂v ……………………….(4.5) σy = – p + λ (---- + ----) + 2μ ---∂x ∂y ∂u ……………………….(4.6) ∂y ∂v τxy = τyx = μ (---- + ----) ∂y …………………………..(4.7) ∂x Combining equations (4.2) and (4.5), (4.6) and (4.7) a form of the “Navier Stokes” equations can be written: ∂u ∂u ∂u 1 1 ∂p 1 μ ∂ ∂u ∂v μ ∂ 2u ∂2u ---- + u ---- + v ---- = --- Fx – --- ---- + --- --- ---- (---- + ----) + --- (----- + -----) ∂t ∂x ∂y ρ ρ ∂x 3 ρ ∂x ∂x ∂y ρ ∂x ∂v ∂v ∂v 1 1 ∂p 1 μ ∂ ∂v μ ∂2v ∂u 2 …(4.8) ∂y 2 ∂2v ---- + u ---- + v ---- = --- Fy – --- ---- + --- --- ---- (---- + ----) + --- (----- + -----) ∂t ∂x ∂y ρ ρ ∂y 3 ρ ∂y ∂x ∂y ρ ∂x 2 …(4.9) ∂y 2 On introduction of the incompressibility condition, these can be further simplified to: ∂u ∂u ∂u 1 1 ∂p μ ∂2u ∂2u ---- + u ---- + v ---- = --- Fx – --- ---- + --- (----- + -----) ∂t ∂x ∂y ρ ρ ∂x ρ ∂x2 ∂y2 ∂v ∂v ∂v 1 1 ∂p μ ∂2v ∂2v ……(4.10) ---- + u ---- + v ---- = --- Fy – --- ---- + --- (----- + -----) ∂t ∂x ∂y ρ ρ ∂x ρ ∂x 2 ……(4.11) ∂y 2 For steady state conditions, the terms (∂u/∂t) and (∂v/∂t) can be dropped resulting in ‘coupled’ equations in the ‘primitive’ variables (u, v and p). The equations are also, in contrast to those of solid elasticity, nonlinear due to the presence of products like u(∂u/∂x). Ignoring body forces for the present, the steady state equations to be solved are: ∂u ∂u 1 ∂p μ ∂2u ∂2u u ---- + v ---- + --- ---- – --- (----- + -----) = 0 ∂x ∂y ρ ∂x ρ ∂x2 ∂y2 ∂v ∂v 1 ∂p μ ∂2v ∂2v ISSN: 0973-6875 ………………….…(4.12) P a g e | 27 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) u ---- + v ---- + --- ---- – --- (----- + -----) = 0 ∂x ∂y ρ ∂y ρ ∂x2 ………………….…(4.13) ∂y2 Proceeding as before, and for the moment assuming the same shape functions are applied to all variables (Smith and Griffiths, 1988): u = N u, v = N v, p = N p We have for a single element, treating the terms (u) and (v) as constants (ū = N uo) and ( purposes of integration: ∂N ∂N ū ---- u + 1 ∂N μ ∂2N ∂y ρ ∂x ρ ∂x ∂N ∂N 1 ∂N μ ∂2N ū ---- v + ∂x μ ∂2N ---- u + --- ---- p – --- ----- u – --- ----- u = 0 ∂x ρ ∂x …………………….(4.14) ρ ∂y 2 2 μ ∂2N ---- v + --- ---- p – --- ----- v – --- ----- v = 0 ∂y = N vo) for the ρ ∂x2 …………………….(4.15) ρ ∂y2 Multiplying by the weighting functions and integrating as usual yields: ∂N ∂N ∫∫ Nū ---- u dx dy + ∫∫ N ∂x 1 ∂N ---- u dx dy + --- ∫∫ N ---- p dx dy ∂y μ ρ ∂N ∂x μ 2 ∂2N – --- ∫∫ N ----- u dx dy – --- ∫∫ N ----- u dx dy = 0 ……………..(4.16) ρ ∂N ∂x2 ∂N ∫∫ Nū ---- v dx dy + ∫∫ N ∂x μ ρ 1 ∂y2 ∂N ---- v dx dy + --- ∫∫ N ---- p dx dy ∂y ρ ∂y ∂2N μ ∂2N – --- ∫∫ N ----- v dx dy – --- ∫∫ N ----- v dx dy = 0 ρ ISSN: 0973-6875 ∂x2 ρ ……………..(4.17) ∂y2 P a g e | 28 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Integrating products by parts where necessary and neglecting resulting contour integrals gives: ∂N ∂N ∫∫ Nū ---- dx dy u + ∫∫ N 1 ∂N μ ∂N ∂N ---- dx dy u + --- ∫∫ N ---- dx dy p + --- ∫∫ N ---- ---- dx dy u ∂x ∂y ρ ∂x μ ρ ∂x ∂x ∂N ∂N + --- ∫∫ N ---- ---- dx dy u = 0 ρ ∂N ∂N ∫∫ Nū ---- dx dy v + ∫∫ N ∂x …………..(4.18) ∂y ∂y 1 ∂N μ ∂N ∂N ---- dx dy v + --- ∫∫ N ---- dx dy p + --- ∫∫ N ---- ---- dx dy v ∂y ρ ∂y μ ρ ∂x ∂x ∂N ∂N + --- ∫∫ N ---- ---- dx dy v = 0 ρ …………..(4.19) ∂y ∂y The set of equations is completed by the continuity condition: ∂N ∂N ∫∫ N (---- u + ---- v) dx dy = 0 ∂x …………………………(4.20) ∂y Collecting terms in (u, p and v) respectively leads to an equilibrium equation (Taylor and Hughes, 1981): ………………………(4.21) Where, ∂N C11 = ∫∫ (Nū ---- + N ∂N μ ∂N ∂N ---- + --- ---- ---- + --- ---- ----) dx dy ∂x 1 μ ∂N ∂N ∂y ρ ∂x ∂x ρ ∂y ∂y ∂N C12 = ∫∫ --- N ---- dx dy ρ ∂x ISSN: 0973-6875 P a g e | 29 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) C13 = 0 ∂N C21 = ∫∫ N ---- dx dy ∂x C22 = 0 ∂N C23 = ∫∫ N ---- dx dy ∂y C31 = 0 1 ∂N C32 = ∫∫ --- N ---- dx dy ρ ∂y C33 = C11 Simplified Flow Equations In many practical instances, it may not be necessary to solve the complete coupled system described in the previous section. The pressure (p) can be eliminated from equations (4.10) and (4.11) and if vorticity (w) is defined as: ∂u ∂v w = ---- – ---∂y ………………………..(4.22) ∂x This results in a single equation: ∂w ∂w ∂w μ ∂2w ∂2w ---- + u ---- + v ---- = --- (----- + -----) ∂t ∂x ∂y ρ ∂x2 …………………………(4.23) ∂y2 Defining a stream function (ψ) such that ∂ψ u = ----- ……………………………………(4.24) ∂y ISSN: 0973-6875 P a g e | 30 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) ∂ψ v = – ---- ……………………………(4.25) ∂x An alternative coupled system involving (ψ) and (w) can be devised, given here for steady state conditions: ∂2ψ ∂2ψ …………………………(4.26) ----- + ----- = w ∂x2 ∂y2 μ ∂2w ∂2w ∂ψ ∂w ∂ψ ∂w --- (----- + -----) = ---- ---- – ---- ---ρ ∂x2 ∂y2 ∂y ∂x …………………………(4.27) ∂x ∂y This clearly has the advantage that only two unknowns are involved rather than the previous three. However, the solution of equations (4.26) and (4.27) is still a relatively complicated process and flow problems are sometimes solved via equation (4.23) alone, assuming that (u) and (v) can be approximated by some independent means. In this form, equation (4.23) is an example of the 'diffusion-convection' equation, the second order space derivatives corresponding to a 'diffusion' process and the first order ones to a 'convection' process. The equation arises in various areas of engineering, for example sediment transport and pollutant disposal (Smith, 1976 and 1979). If there is no convection, the resulting equation is of the type: ∂w μ ∂2w ∂2w ---- = --- (------ + ------) ∂t ρ ∂x2 ……………………………(4.28) ∂y2 which is the 'heat conduction' or 'diffusion' equation well known in many areas of engineering. A final simplification is a reduction to steady state conditions, in which case: ∂2w ∂2w ------ + ------ = 0 ISSN: 0973-6875 ……………………………(4.29) P a g e | 31 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) ∂x2 ∂y2 Leaving the familiar "Laplace" equation. In the following sections, finite element formulations of these simplified flow equations are described, in order of increasing complexity. Simplified Fluid Flow: Steady State The form of Laplace’s equation (4.29) which arises in geomechanics, for example concerning groundwater flow in an aquifer (Muskat, 1937), is: ∂2Ф ∂2Ф kx ------ + ky ------ = 0 ∂x2 ……………………………(4.30) ∂y2 where (Ф) is the fluid potential and (kx), (ky) are permeabilites in the (x) and (y) directions. The finite element discretisation process reduces the differential equation to a set of equilibrium type simultaneous equations of the form: KP. Ф = Q ………………………….(4.31) where (KP) is symmetrical and (Q) is a vector of net nodal inflows/outflows. The typical terms in the matrix (KP) are of the form (Smith and Griffiths, 1988): ∂Ni ∂Nj ∂Ni ∂Nj ∫∫ (kx ---- ----- + ky ----- -----) dx dy ∂x ∂x ………………………………(4.32) ∂y ∂y With the usual finite element discretisation: Ф = NФ ……………………………….(4.33) A convenient way of expressing the matrix (KP) in equation (4.31) is: KP = ∫∫ TT K T dx dy ……………………………(4.34) where the property matrix (K) is analogous to the stress strain matrix (D) in solid mechanics; thus: ISSN: 0973-6875 P a g e | 32 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) …………………………(4.35) assuming that the axes of the permeability tensor coincide with (x) and (y). The (T) matrix is similar to the (B) matrix of solid mechanics and is given by: …………………..(4.36) where: (Ni) are shape functions. Finally, it is worth noting that equation (4.34) can also be arrived at from energy considerations. The equivalent energy statement is that the integral: 1 ∂Ф 1 ∂Ф ∫∫ [--- kx (-----)2 + --- ky (----)2 ] dx dy 2 ∂x 2 ………………….(4.37) ∂y Simplified Fluid Flow with Advection If pollutants, sediments, tracers, etc., are transported by a laminar flow system, they are at the same time translated or 'advected' by the flow and diffused within the governing differential equation for the twodimensional case is (Smith et al., 1973): ∂2Ф ∂2Ф ∂Ф ∂Ф ∂Ф cx ----- + cy ------ – u ---- – v ---- = ----∂x2 ∂y2 ∂x ∂y ………………………(4.38) ∂t where (u) and (v) are the fluid velocity components in the (x) and (y) directions (compare equation (4.23)). The extra advection terms –u(∂Ф/∂x) and –v(∂Ф/∂y) compared with equation (4.38) lead to unsymmetric components of the "stiffness" matrix of the type: ISSN: 0973-6875 P a g e | 33 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) ∂Nj ∂Nj ∫∫ (–u Ni ----- – v Ni -----) dx dy ∂x ……………………..(4.39) ∂y which must be added to the symmetric diffusion components given in equation (4.32). When this has been done, equilibrium equations like equation (4.31) are regained. Mathematically, equation (4.38) is a differential equation which is not self-adjoint (Berg, 1962), due to the presence of the first order spatial derivatives from a finite element point. A second consequence of non-self-adjoint equations is that there is no energy formulation equivalent to equation (4.37). It is clearly a benefit of the Galerkin approach that it can be used for all types of equation and is not restricted to self-adjoint systems. Equation (4.38) can be rendered self-adjoint by using the transformation: ux vy Ф = exp (----) exp (----) 2cx …………………(4.40) 2cy Volumetric Water Content Functions Fundamental to the formulation of a general seepage analysis is an understanding of the relationship between pore-water pressure and water content. As water flows through soil, certain amounts of water are stored or retained within the soil structure. The amount of water stored or retained is a function of the pore-water pressure and the characteristics of the soil structure. For a seepage analysis, it is convenient to specify the stored portion of the water flow as a ratio of the total volume. This ratio is known as the volumetric water content. In equation form: Θ = Vw/V ……………………….(4.41) where: Θ = volumetric water content, Vw = volume of water, and V = total volume. ISSN: 0973-6875 P a g e | 34 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) The volumetric water content (Θ) is dependent on the pore-water pressure. Figure (4.1) illustrates this relationship, which is also known as the soil-water characteristic function, (Fredlund, and Rahardjo, 1993). When the degree of saturation is (100%), the volumetric water content is equivalent to the soil porosity, which is defined as the volume of voids divided by the total volume. Considering a completely saturated soil where the pore-water pressure is near zero, and the total external load on the soil remains constant, as the pore water pressure becomes more positive, the effective stress will decrease. This causes the soil to swell and results in an increase in its water content. As the porewater pressure becomes negative, the soil begins to desaturate and the water content decreases. Ultimately, the soil becomes completely desaturated, and the water content no longer changes with a further decrease in pore-water pressure. The slope of the soil-water characteristic curve (designated as mw) represents the rate of change in the amount of water retained by the soil in response to a change in pore-water pressure. When the pore-water pressure is positive, (mw) is equivalent to (mv), the coefficient of compressibility for one dimensional consolidation. The parameter (mw) is required in a transient seepage analysis. Figure (4.1): General Form of Volumetric Water Content Function (after Fredlund and Rahardju, 1993). Soil-water characteristic functions for fine-grained (clay) soils may be relatively flat, while for coarsegrained soils (sand), the function may be quite steep. Figure (4.2) presents actual volumetric water content ISSN: 0973-6875 P a g e | 35 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) curves obtained by Ho, (1979), for fine sand, silt, and clay. The variation in these curves demonstrates the effect that the soil properties have on the characteristic functions. Hydraulic Conductivity Functions Water in liquid form can be considered to flow along a web of interconnected but continuous conduits. Decreasing the water content has the effect of decreasing the size and number of conduits, thereby reducing the capacity to conduct water through the soil. Ultimately, when the soil is dry, the capacity to conduct water along continuous water-filled conduits disappears. When the soil is saturated, all of the available conduits are utilized, and the water hydraulic conductivity capacity is at a maximum. Figure (4.2): Actual Volumetric Water Content Functions for Fine Sand, Silt, and Clay (after Ho, 1979). The capacity of soil to conduct water can be viewed in terms of hydraulic conductivity (or the coefficient of permeability). In this context, the hydraulic conductivity is dependent on the water content. Since the water content is a function of pore-water pressure and the hydraulic conductivity is a function of water content, it follows that hydraulic conductivity is also a function of pore water pressure. Figure (4.3) presents a curve showing a typical relationship between hydraulic conductivity and pore-water pressure. ISSN: 0973-6875 P a g e | 36 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Techniques have been developed for predicting the hydraulic conductivity function from a soil-water characteristic function. Establishing the characteristic function is generally not as complicated as conducting laboratory tests to measure the conductivity function. Figure (4.3): A Hydraulic Conductivity Function. The program SEEP/W uses the Green and Corey, (1971), procedure for estimating the conductivity function from a soil-water characteristic function. Defining the hydraulic conductivity for negative pore-water pressure regions makes it possible to analyze problems involving unsaturated flow as well as saturated flow. Flow Law The program SEEP/W is formulated on the basis that the flow of water through both saturated and unsaturated soil follows Darcy's Law which states that: q=ki ……………………(4.42) where: q = specific discharge, ISSN: 0973-6875 P a g e | 37 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) k = hydraulic conductivity, and i = gradient of fluid head or potential. Darcy's Law is often written as: ………………………..(4.43) v=ki where (v) is known as the Darcian velocity. The actual average velocity at which water moves through the soil is the Darcian velocity divided by the porosity of the soil (Harr, 1962). Governing Equations The governing differential equation used in the formulation of flow in the Program SEEP/W is: ∂ ∂H ∂ ∂H ∂Θ ---- (kx -----) + ---- (ky -----) + Q = ----∂x ∂x ∂y ∂y …………………….(4.44) ∂t where: H = total head, kx = hydraulic conductivity in the x-direction, ky = hydraulic conductivity in the y-direction, Q = applied boundary flux, Θ = volumetric water content, and t = time. This equation states that the difference between the flow (flux) entering and leaving an elemental volume at a point in time is equal to the change in the volumetric water content. More fundamentally, it states that the sum of the rates of change of flows in the (x) and (y) directions plus the external applied flux is equal to the rate of change of the volumetric water content with respect to time. ISSN: 0973-6875 P a g e | 38 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Under steady-state conditions, the flux entering and leaving an elemental volume is the same at all times. The right side of the equation consequently vanishes and the equation reduces to: ∂ ∂H ∂ ∂H ---- (kx -----) + ---- (ky -----) + Q = 0 ∂x ∂x ∂y …………………….(4.45) ∂y Changes in volumetric water content are dependent on changes in the stress state and the properties of the soil. The stress state for both saturated and unsaturated conditions can be described by two state variables (Fredlund and Morgenstern, 1976 and Fredlund and Morgenstern, 1977). These stress state variables are (ua) and (ua - uw) where is the total stress, (ua) is the pore-air pressure, and (uw) is the pore-water pressure. The program SEEP/W is formulated for conditions of constant total stress; that is, there is no loading or unloading of the soil mass. The second assumption is that the pore-air pressure remains constant at atmospheric pressure during transient processes. This means that (σ – ua) remains constant and has no effect on the change in volumetric water content. Changes in volumetric water content are consequently dependent only on changes in the (ua – uw) stress state variable, and with (ua) remaining constant, the change in volumetric water content is a function only of pore-water pressure changes. A change in volumetric water content can be related to a change in pore water pressure by the equation: ∂Θ = mw ∂uw ………………..(4.46) where (mw) is the slope of the storage curve. The total hydraulic head is defined as : H = uw/γw + y ………………(4.47) where: uw = pore-water pressure, γw = unit weight of water, and ISSN: 0973-6875 P a g e | 39 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) y = elevation. Equation (4.47) can be rearranged as: uw = (H – y) ………………….(4.48) Substituting equation (4.48) into (4.46) gives the following equation: ∂Θ = mw γw ∂(H – y) ………………..(4.49) Now can be substituted into equation (4.44), leading to the following expression: ∂ ∂H ∂ ∂H ∂(H – y) ---- (kx -----) + ---- (ky -----) + Q = mw γw ---------∂x ∂x ∂y ∂y …………………….(4.50) ∂t Since the elevation is a constant, the derivative of (y) with respect to time disappears, leaving the following governing differential equation: ∂ ∂H ∂ ∂H ∂H ---- (kx -----) + ---- (ky -----) + Q = mw γw ----∂x ∂x ∂y ∂y …………………….(4.51) ∂t Finite Element Formulation of the Flow Problem Coordinate Systems The global coordinate system used in the formulation of SEEP/W is the first quadrant of a conventional (x, y) Cartesian system. The local coordinate system used in the formulation of element matrices is presented in Figure (4.4). In the case of quadrilateral elements, nodes 5, 6, 7, and 8 are secondary nodes. In the case of triangular elements, nodes 5, 6, and 7 are secondary nodes. The local and global coordinate systems are related by a set of interpolation functions. SEEP/W uses the same functions for relating the coordinate systems as for describing the variation of the field variable (head) within the element. The elements are consequently isoparametric elements. The (x) and (y) coordinates anywhere in the element are related to the local coordinates and to the (x, y) coordinates of the nodes by the following equations: ISSN: 0973-6875 P a g e | 40 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) x = <N> {X} ………………(4.52) y = <N> {Y} ………………..(4.53) where <N> is a vector of interpolating shape functions and {X} and {Y} are the global (x, y) coordinates of the element nodes. The interpolating functions are expressed in terms of local coordinates. Therefore, once a set of local coordinates (r, s) have been specified, the corresponding global coordinates can be determined by equation (4.52) and equation (4.53). Figure (4.4): Global and Local Coordinate Systems. Interpolating Functions The program SEEP/W uses a general set of interpolating functions presented by Bathe, (1982). These general functions are suitable for elements which have none, some, or all of the secondary nodes defined. This allows for considerable versatility in the types of elements that can be used. The interpolating functions in terms of local coordinates (r) and (s) for quadrilateral and triangular elements are given in Appendix (A), Tables (A.1) and (A.2), respectively. The functions represent a linear equation when the secondary nodes are missing and a quadratic (nonlinear) equation when the secondary nodes are included. Field Variable Model ISSN: 0973-6875 P a g e | 41 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) To formulate a finite element analysis, it is necessary to adopt a model for the distribution of the field variable within the element. Since the field variable in the seepage analysis is the total head (H), it is necessary to adopt a model for the distribution of (H) within the element. The program SEEP/W assumes that the head distribution within the element follows the adopted interpolating functions. This means that the head distribution is linear when the secondary nodes are missing, and the head distribution is nonlinear when the secondary nodes are present. In equation form, the head distribution model is: ……………………(4.54) h = <N> {H} where: h = head at any local coordinate, <N> = vector of interpolation function, and {H} = vector of heads at the nodes. Interpolation Function Derivatives The constitutive relationship for a seepage analysis is Darcy's Law. q=ki The gradient (i) is one of the key parameters required in the finite element formulation. The following presents the procedure used by SEEP/W to compute the gradient. From the adopted head distribution model, the head at any point within the element in terms of the nodal heads is: h = <N> {H} ……………………(4.55) The gradients in the (x) and (y) directions are ∂h ∂ <N> ix = ---- = --------- {H} ∂x ISSN: 0973-6875 ………………….(4.56) ∂x P a g e | 42 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) ∂h ∂ <N> ………………….(4.57) iy = ---- = --------- {H} ∂y ∂y The interpolating functions are written in terms of (r) and (s) and not in terms of (x) and (y). The derivatives must consequently be determined by the chain rule of differentiation, as follows: ∂h ∂ <N> ………………….(4.58) iy = ---- = --------- {H} ∂y ∂ <N> ∂y ∂ <N> ∂x ∂ <N> ∂y -------- = -------- ----- + --------- ---∂s ∂x ∂s ∂y ……………………(4.59) ∂s This can be written as: ……………………………(4.60) where [J] is the Jacobian matrix and is defined as: ………………………….(4.61) The derivative of the interpolation function with respect to (x) and (y) is called the (B) matrix and can be determined by inverting the Jacobian matrix and rewriting the equations as: ISSN: 0973-6875 P a g e | 43 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) ………………………(4.62) Recall from equations (4.52) and (4.53) that: x = <N> {X} y = <N> {Y} ∂x ∂ <N> ----- = -------- {X} ∂r ∂r ∂x ∂ <N> ----- = -------- {X} ∂s ∂s ∂y ∂ <N> ----- = -------- {Y} ∂r ∂r ∂y ∂ <N> ----- = -------- {Y} ∂s ∂s Substituting these values into equation (4.61), the Jacobian matrix becomes: ISSN: 0973-6875 P a g e | 44 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) ………………………(4.63) The derivatives of the interpolation functions with respect to (r) and (s) are required to compute the Jacobian matrix (equation 4.63) and to compute the flow gradients (equations 4.56, 4.57, and 4.62). The derivatives of the interpolation functions with respect to (r) and (s) used by SEEP/W for quadrilateral and triangular elements are given in Appendix (A), Tables (A.3) and (A.4), respectively. The following notation is used in the preceding tables: ∂Ni Ni, r = ---- ; ∂Ni Ni, s = ----- ∂r ∂s The Jacobian matrix is a (2 x 2) matrix …………..(4.64) The inverse of [J] is: ………………(4.65) The determinant of [J] is: det [J] = (J11 x J22 – J21 x J12) ………………….(4.66) Finite Element Equations The finite element equation that follows from applying the Galerkin method of weighted residual to the governing differential equation (equation (4.51)) is: ∫v ([B]T [C] [B]) dv {H} + ∫v (λ <N>T <N>) dv {H}, t = q ∫A (<N>T) dA …………..(4.67) where: ISSN: 0973-6875 P a g e | 45 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) [B] = gradient matrix, [C] = element hydraulic conductivity matrix, {H} = vector of nodal heads, λ = mw λw <N>T <N> = [M] = mass matrix {H}, t = ∂h/∂t = change in head with time q = unit flux across the side of an element, and <N> = vector of interpolating function. For a SEEP/W two-dimensional analysis, the thickness of the element is considered to be constant over the entire element. The finite element equation can consequently be written as: t1 ∫A ([B]T [C] [B]) dA {H} + t1 ∫A (λ <N>T <N>) dA {H}, t = q t1 ∫L (<N>T) dL …..(4.68) where (t1) is the element thickness. When (t1) is a constant, the integral over the volume becomes the integral over the area and the integral over the area becomes the integral over the length from corner node to corner node. In an axisymmetric analysis, the equivalent element thickness is the circumferential distance about the symmetric axis. The complete circumferential distance is (2) radian times (R). Since SEEP/W is formulated for one radian, the equivalent thickness is (R). The finite element equation for the axisymmetric case is: ∫A ([B]T [C] [B] R) dA {H} + ∫A (λ <N>T <N> R) dA {H}, t = q ∫L (<N>T R) dL …..(4.69) The radial distance (R) is not a constant as is the thickness (t1) for a two-dimensional analysis within an element; consequently, (R) is a variable in the integration. In abbreviated form, the finite element equation is: [K] {H} + [M] {H}, t = {Q} …………………..(4.70) where: ISSN: 0973-6875 P a g e | 46 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) [K] = element characteristic matrix = t1 ∫A ([B]T [C] [B]) dA, or = ∫A ([B]T [C] [B] R) dA [M] = mass matrix = t1 ∫A (λ <N>T <N>) dA, or = ∫A (λ <N>T <N> R) dA {Q} = applied flux vector = q t1 ∫L (<N>T) dL, or = q ∫L (<N>T R) dL Equation (4.70) is the general finite element equation for a transient seepage analysis. For a steady-state analysis, the head is not a function of time and consequently the term [M]{H},t vanishes, reducing the finite element equation to: [K] {H} = {Q} ……………………(4.71) Time Integration The finite element solution for a transient analysis is a function of time as indicated by the {H},t term in the finite element equation. The time integration can be performed by a finite difference approximation scheme. Writing the finite element equation in terms of finite differences leads to the following equation (Segerlind,1984): (w Δt [K] + [M]) {H1} = Δt ((1 – w) {Q0} + w{Q1}) + ([M] – (1 – w) Δt [K]) {H0} ….(4.72) where: Δt = time increment, w = a ratio between (0) and (1), ISSN: 0973-6875 P a g e | 47 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) {H1} = head at end of time increment, {H0} = head at start of time increment, {Q1} = nodal flux at end of time increment, {Q0} = nodal flux at start of time increment, [K] = element characteristic matrix, and [M] = element mass matrix. SEEP/W uses the Backward Difference Method, a method that sets (Θ) to (1.0). The transient finite equation for (Θ) equal to (1.0) is: (Δt [K] + [M]) {H1} = Δt {Q1} + [M] {H0} ………………..(4.73) As indicated by this equation, in order to solve for the new head at the end of the time increment, it is necessary to know the head at the start of the increment. Stated in general terms, the initial conditions must be known in order to perform a transient analysis. Numerical Integration SEEP/W uses Gauss numerical integration to form the element characteristic matrix [K] and the mass matrix [M]. The integrals are sampled at specifically defined points in the elements and then summed for all the points. The following integral (from equation (4.70)): t ∫A ([B]T [C] [B]) dA can be replaced by: n Σ [Bj]T [Cj] [Bj] det |Jj | W1j W2j ……………………….(4.74) j=1 where: j = integration point, ISSN: 0973-6875 P a g e | 48 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) n = number of integration points, det |Jj | = determinant of the Jacobian matrix, and W1j, W2j = weighting factors. The number of sample (integration) points required in an element depends on the number of nodes and the shape of the elements. One-point integration for a triangular element results in a constant gradient throughout the element. The number of integration points is denoted as the integration order. The appropriate integration order is a function of the presence of secondary nodes. When secondary nodes are present, the interpolating functions are nonlinear and consequently a higher integration order is required. It is also acceptable to use four-point integration for quadrilateral elements which have secondary nodes. This is called a reduced integration order (Bathe, 1982). Acceptable results can be obtained with reduced integration. For example, reduced integration is useful in saturated zones where the hydraulic gradient is low and the hydraulic conductivity is constant. Selective use of reduced integration can greatly reduce the required number of computations. Hydraulic Conductivity Matrix The general form of the SEEP/W hydraulic conductivity matrix is: …………………..(4.75) where: C11 = kx cos 2α + ky sin 2α C22 = kx sin 2α + ky cos 2α C12 = kx cos α sin α – ky sin α cos α C21 = C12 ISSN: 0973-6875 P a g e | 49 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) The parameters (kx), (ky), and are defined in Figure (4.5). Figure (4.5): Definition of Hydraulic Conductivity Matrix Parameters. When (α) is zero, [C] reduces to: …………….(4.76) The parametric (kx) is always determined from the hydraulic conductivity function. Parameter (ky) is then computed as (kx) multiplied by K-Ratio. In equation form, ky = kx . K-Ratio ……………………(4.77) Mass Matrix As first presented in equation (4.70), the element mass (or storage) matrix for a two-dimensional analysis is defined as: [M] = t1 ∫A (λ <N>T <N>) dA and for an axisymmetric analysis as: [M] = ∫A (λ <N>T <N> R) dA where: λ = mw. γw, mw= slope of storage curve, ISSN: 0973-6875 P a g e | 50 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) t1 = element thickness, R = radial distance from the symmetric axis to the integration point, and <N> = vector of interpolating functions. SEEP/W uses a lumped formulation, to establish the mass matrix. The details of this formulation are given by Segerlind, (1984). The mass matrix is formed by numerical integration. This has the advantage of making it possible for the parameter (mw) to vary throughout the element. SEEP/W obtains an (mw) value from the storage function for each integration point. SEEP/W computes (mw) from the slope of a straight line between the old and new pore-water pressures at a Gauss point, as illustrated in Figure (4.6). Figure (4.6): Computation of (mw). The slope of this straight line can be viewed as the average rate of change during one increment of time. This is considered to be a more realistic value than taking the derivative of the function at a specific point. An exception to this procedure is when the old and new pore-water pressures are nearly identical. In this case, SEEP/W computes (mw) by calculating the derivative of the function at the average of the old and new porewater pressures. Flux Boundary The nodal boundary flux (Q) for a two-dimensional analysis is defined as: ISSN: 0973-6875 P a g e | 51 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) q t1 ∫L (<N>T) dL ……………(4.78) for an axisymmetric analysis as: q ∫L (<N>T R) dL ……………(4.79) and for a plan view analysis as: q ∫L (<N>T R) dL ……………(4.80) where: {Q} = vector of nodal flux, q = unit flux across the side of an element, t1 = element thickness, A = area of the element, R = radial distance from the symmetric axis to the element corner nodes, and <N> = vector of interpolating function. Solutions to the integrals are dependent on the analysis type and on the presence of secondary nodes. For two-dimensional (i.e., vertical section) and axisymmetric analyses, the integrals are solved by close form solutions as illustrated in Figure (4.7). However, for plan view analysis, the contributing area of a node is computed by numerical integration in the same way as forming the mass matrix. ISSN: 0973-6875 P a g e | 52 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Figure (4.7): Contributing Areas for Vertical Two-Dimensional (Plane Strain) Elements with Width Equal to (1) Unit. Two types of flux boundaries may be specified in SEEP/W namely: a nodal flux boundary (Q) and a unit flux boundary (q). A nodal flux boundary (Q) can be specified directly on the boundary nodes. A unit flux boundary (q) must be specified along the boundary edges of the elements, except for a plan view analysis. It is required to identify the edges of the elements across which a (q) boundary should be applied. Based on this specific element edge information, the integration is performed and determines the applied flux (Q) at the nodes. Assembly and Solving of Global Equations The element matrix for each element in the discretized finite element mesh can be formed and assembled into a global system of simultaneous equations. The finite element solution requires the solving of the system of simultaneous equations. SEEP/W stores the coefficients of the global system equations using a compressed row storage scheme (Barrett et. al, 1994). This is a general scheme that makes no assumptions about the sparsity structure of the matrix. Instead of a full matrix with many zero elements, the compressed row storage scheme stores the matrix with (3) vectors: one for the non-zero elements, one for the column index and one for the row pointers. As a result, it provides significant saving in storage memory particularly in large finite element meshes. The equation solver can accommodate missing elements in the array. This feature makes it possible to add and delete elements from a finite element mesh without renumbering the nodes and elements. Unlike the Gauss elimination solver, sorting the node number vertically or horizontally will have no effect to the computational time and accuracy of the solution with this improved iterative solver. As a result of the iterative solver, the solution of the system of simultaneous equations is to a large degree dependent on the convergence tolerance and the number of iterations. The Gauss elimination skyline direct solver is also included in the program. The computational speed of a direct solver is depending on the bandwidth of the finite element mesh. The direct solver is fast for simple problems, but it can be slower than the iterative solver in complex problems with large bandwidths. Gradients and Velocities ISSN: 0973-6875 P a g e | 53 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Once the solution has converged and the nodal heads are known, SEEP/W computes the hydraulic gradients and Darcian flow velocities at each of the integration points within each element. The gradient at each Gauss or integration point is computed from the equation: ………………..(4.81) where: ix = gradient in (x) direction, iy = gradient in (y) direction, [B] = gradient matrix as defined in finite element equations, and {H} = vector of total head at the nodes. The Darcian velocities at each Gauss point are computed from the Equation: ……………………..(4.82) where: vx = velocity in (x) direction, vy = velocity in (y) direction, and [C] = hydraulic conductivity matrix. SEEP/W stores in an array the hydraulic conductivity at each Gauss point used in the formulation of the finite element equations. The same hydraulic conductivity values are later used to compute the velocities. The SEEP/W velocity is actually the specific discharge, which is the total flux (Q) divided by the full cross-sectional area (voids and solids alike); it is not the actual speed with which the water moves between the soils particles. The actual microscopic velocity is: v = q/n ISSN: 0973-6875 ………………..(4.83) P a g e | 54 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) where: v = average linear velocity, q = specific discharge, and n = porosity Flow Quantities SEEP/W has the ability to compute the seepage quantity that flows across a user-defined section. This quantity can be computed from the nodal heads and the coefficients of the finite element equation. For example, considering a mesh with only one element, as illustrated in Figure (4.8). The objective is to compute the total flow across a vertical section of the element. Equation (4.73) can be rewritten as: ΔH [K] {H} + [M] ---- = {Q} …………………..(4.84) Δt Figure (4.8): Illustration of a Flux Section. ISSN: 0973-6875 P a g e | 55 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) In a steady-state analysis, the storage term [M] ΔH/Δt becomes zero, and the equation can be reduced to: …………(4.85) [K] {H} = {Q} The global set of finite equations for one element is as follows: ……………………..(4.86) From Darcy’s Law (equation (4.42)), the total flow between two points is: ΔH Q = k A ---- ……………………………..(4.87) L The coefficients (c) are a representation of (k.A/l) in equation (4.86). Therefore, the flow from node (i) to node (j) is: Qij = Cij (Hi – Hj) ……………………..(4.88) In a transient analysis, because of material storage, the calculation of the total flow quantity must include the storage effect. The change in flow quantity due to the storage term can be expressed as: ………………….(4.89) where ΔH1, ΔH2, ΔH3, and ΔH4 are the changes of total head at the various nodes between the start and the end of a time step. In general, the average change of total head from node (i) to node (j) can be expressed as: ΔHij = (ΔHi + ΔHj)/2 ISSN: 0973-6875 ……………………………(4.90) P a g e | 56 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Therefore, the change in flow quantity from node (i) to node (j) due to a change in storage is: ΔHij ………………..(4.91) Qij = mij -----Δt The total flow quantity from node (i) to node (j) for a transient analysis then becomes: ΔHij Qij = cij (Hi – Hj) + mij ------ ……………….(4.92) Δt The total flow quantity through the flux section shown in Figure (4.8) is: Q = Q21 + Q24 + Q31 + Q34 …………………………(4.93) 5. RESULT Estimation of Seepage beneath Sheet Pile in Homogenous Soil Figure (5.6). Infinite elements are used on the left and right vertical boundaries. Figures (5.4) and (5.5) present the dimensions of the problem and the flow net solution. The permeability of the soil material is (1x10-4 m/sec). The boundary nodes of the upstream and downstream surfaces are designed as head boundaries with total head equals to (2.5 m) and (0.5 m) respectively. Default boundary conditions (no flow) are assumed for all other boundaries. Figure (5.7) presents the SEEP/W computed head distribution and flow vectors, while Figure (5.8) shows the counters of total head and the velocity flow. The value of seepage is (0.933 x 10-4), the pore pressure at point (A) is (23.11 kN/m2), and the exit gradient is (0.22194). Based on the flow net solution, the seepage under the wall is (1x10-4 m3/sec/m), the pore pressure at point (A) is (23.3 kN/m2), and the exit gradient is (0.23). The steady state seepage of the sheet pile problem is analyzed using SEEP/W. The finite element mesh used for the analysis is shown in ISSN: 0973-6875 P a g e | 57 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Figure (5.4): The Dimensions of sheet pile wall problem Figure (5.5): The Flow Net Solution of the Problem ISSN: 0973-6875 P a g e | 58 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Figure (5.6): The Finite Element below the sheet pile wall Figure (5.7): The Vectors of Flow. ISSN: 0973-6875 P a g e | 59 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) Figure (5.8): The Flow Net and Contours of Total Head as calculated by Finite Element Method. To this point, it can be concluded that the program SEEP/W can successfully simulate flow problems through porous media. 6. CONCLUDING REMARKS From the present study it can be concluded that, 1. SEEP/W program could successfully simulate the flow problems through porous media under the sheet pile. 2. The seepage analysis using SEEP/W provides the flow net diagram and water head at different points below the sheet pile wall. Amount of seepage loss can also be estimated using this software. 3. The amount of seepage is decreased when the embedment depth of sheet pile is increased. 4. The above information is useful for the design of water front sheet pile wall. REFERENCES [1] Chandra Satish and Prakash Anand, “A Study of 3-D Seepage below Hydraulic Structures with Cutoffs.” [2] Chen Qun and Zhang L. M., “Three-dimensional analysis of water infiltration into the Gouhou rockfill dam using saturated-unsaturated seepage theory” [3] Das. B. M., “Advance Soil Mechanics”. McGraw Hill International Edition, 1983. [4] Desai C. S., “Free Surface Seepage Through Foundation and Berm of Cofferdams.” [5] Desai M., and Li G. C,. “Stress and Seepage Analysis of Earth Dams.” [6] Feng Zheng-yi and Jonathan T. H. Wu, “The epsilon method: analysis of seepage beneath an impervious dam with sheet pile on a layered soil.” [7] Mondal. S., “Finite Element Analysis of Seepage Flow Under Cofferdam”. ME Thesis submitted to Bengal Engineering College (DU), Howrah 711 103. ISSN: 0973-6875 P a g e | 60 Ajit Kumar Paria./ Int.J.HIT.TRANSC:ECCN. Vol.7: Issue 1A (2021) [8] Rama Rao B. S. and Rao J. Visweswara, “Seepage into Sheet Pile Cofferdam”, Journal of the Hydraulic division, September 1973 [9] Reddy A. S., Sridharan K., Kulkarni V. G., “Seepage under a Dam on Sloping Soil with Discontinuities in Slopes. [10] Manna M C, Bhattacharya A K, Choudhury S, Maji S C, “Groundwater Flow Beneath a Sheetpile Analyzed using Six-noded Triangular Finite Elements.” [11] Tomlin Geoffrey R., “Seepage analysis through zoned anisotropic soils by computer”. [12] Suriyachat Phanuwat, “Flow Characteristics under Sheet Pile in Anisotropic Porous Soil”. [13] Schroeder W. L. and Roumillac Philippe, “Anchored bulkheads with sloping Dredge lines”. [14] Basak Punyabrata, “Non-Darcy flow and its implications to seepage problems”. [15] Griffiths D.V. and Fenton Gordon A., “Seepage beneath water retaining structures founded on spatially random soil”. [16] Griffiths, D. V., and Fenton, G. A., (1997), "Three-Dimensional Seepage through Spatially Random Soil", Journal of Geotechnical and Geoenvironmental Engineering, ASCE, Vol. 123, No. 2, February. [17] Hadi, A. A., (1978), "Seepage through Embankments into Adjacent Drains". M.Sc. Thesis, College of Engineering, University of Baghdad, Iraq. [18] Harr, M. E., (1962), "Groundwater and Seepage", McGraw-Hill Book Company, New York. ISSN: 0973-6875 P a g e | 61