Theoretical Study of the Effect of Permeability Tensor upon Drainage of Soils

the influence of permeability tensor upon drainage of anisotropic soils under ponded water and steady recharge (rainfall) is theoretically investigated. Tensorial permeability has led to the formulation of mixed type partial differential equations. Since there is no analytical solution to this problem, the formulation is therefore solved numerically by the method of finite elements. The finite element formulation is implemented into a computer model which can be applied to any problem of seepage under steady state conditions. Two different example problems representing two different flow conditions under full anisotropy have been studied. Results of the model for the isotropic case were checked against exact mathematical solutions derived analytically for isotropic soils and found to be accurate which indicates that using this model for anisotropic soils is safe and sound.


1-INTRODUCTION
A traditional method of reclaiming salt affected soils involves ponding water on the field to leach the salts from the soil trough subsurface tile drainage system. Subsurface drainage provides a method to remove the excess water from a waterlogged area. Kirkham (1949) developed an exact mathematical solution for the flow of ponded water into drain tubes embedded in isotropic soil. Toksoz and Kirkham (1971), presented an analytical solution; for steady drainage of layered soils. Theoretical formulas for the height of all points of an arched shape phreatic surface were derived for layered soil drained by tubes. Hinesly and Kirkham (1966) solved the problem of simultaneous upward and downward seepage to tube drains in a homogeneous soil analytically. Najamii, Kirkham and Dougal (1987) presented a potential theory solution, of steady state seepage, for equally spaced horizontal tubes in a stratified two layer soil replenished by upward and downward recharge. As well as the influence of anisotropy of soils upon drain spacing, they assumed; that the principal axes; of anisotropy coincides with the global Cartesian coordinates. Bazarra et.al. (1986) attempted to study the influence of artesian as well as of anisotropy of soils upon drain spacing; they assumed that the principal axes of anisotropy coincide with the global Cartesian coordinates. Miles and Kitmitto (1989), derived a new drain formula, which they suggest its applicability to anisotropic and layered soils. Fipps, G. and Skaggs, R.W. (1991) presented a numerical solution to the Richards equation for combined saturated and unsaturated flow towards subsurface pipe drains. They used their numerical solutions to illustrate various aspects of the drainage process and, to evaluate some simple analytical expressions for calculating drain flow rates and water table elevations. They also presented an approximate solution for predicting transient drainage events. The predictions of the approximate solution were found to be in agreement with the numerical solution for elliptic water tables. Sharma, H.C. et.al (1991) studied analytically and experimentally the problem of ditch drainage in layered soils. By using Girinsky's potential theory they presented a mathematical solution for subsurface drainage in layered soils with different boundary conditions. They conducted experimental work using Hele-Shaw model to verify the results of their analytical solution. Experimental results expressed in terms of drain discharges and water table heights were found to be in good agreement with those obtained from the analytical solution. By using potential theory Barua and Tiwari (1996) derived an analytical solution to the problem of seepage of steady rainfall or irrigation recharge into homogeneous and anisotropic soil medium drained by equally spaced ditch drains resting on (a) an impermeable barrier, (b) a gravel substratum.
The validity, of their proposed solutions, expressed in terms of midway water table heights, were tested against those of the upper and lower bounds, of the water table heights produced by Youngs (1965) and found to be very reliable. In the study of Barua and Tiwari, soil anisotropy was simply defined as the permeability's (Kx, Ky) in the horizontal and vertical directions which were taken to be coinciding with the principal axes of the flow domain. Boger (1998) developed an analytical solution for the problem of transient threedimensional unsaturated, anisotropic seepage of layered slopes with and without the effects of capillary forces. In his study Boger showed the influence of the permeability tensor and concluded that seepage flow depends on the anisotropy of the soil. It further depends on the slope of the principal axis of the permeability tensor which has led to a new description of lateral flow mechanism in slopes.
The primary objective of this study is to investigate the influence of tensorial permeability upon drainage of soils. When the permeability of the soil is described as a full tensor there is no certain drain spacing or drain discharge formula that can be applied for the calculations of spacing's and discharges. Therefore, it is intended to simulate numerically a given set of drainage design parameters (s, r, P, d, and e) and check their influence upon the height of the phreatic surface above drain level. In modeling flow to horizontal subsurface conduits (drains), it is necessary to obtain a proper procedure for representing geometrically the conduit. Vimoke and Taylor (1962) represented the subsurface drain as a single point in their electrical resistance network studies. From the wave guide analogy, they calculated adjustment coefficients (Cd) to adjust the values of the resistors connected to the grid point representing the drain. The values where found to be function of the effective drain radius (r) and the size (s) of the square mesh surrounding the drain center. In this study the drain is represented as a single point in the finite element grids and the Vimoke and Taylor adjustment factor (Cd) is divided by two to obtain solutions that agreed well with the analytical solutions of Kirkham (1949) and of Toksoz and Kirkham (1971) for isotropic soils. Siyal, Skaggs and Van Genuchten (2010) used Software HYDRUS 2D/3D to analyze the influence of partial ponding on reclamation of salt affected soils. Their study showed that areas above ponded tile drains were more quickly leached than areas between the ponded drains due to more rapid infiltration and shorter travel distances along such areas. Results of the computer model of Siyal et.al were very consistent with the theoretical studies of Youngs and Harrison (2000) and the experimental results of the sand tank model of Mirjat et.al (2008).

2-THEORY
The components of the discharge (q) and the hydraulic gradient (grad h) in the global coordinate system (x,y) are related to those in the local coordinate system (x',y') through the following transformations Wang and Anderson (1982): (1) Where the rotation matrix [R] is defined as In the local coordinate system the generalized form of Darcy's law is: (4) Substituting equations (1) and (2) into (4) and multiplying across by -1 gives: Equation (5) is the matrix form of Darcy's law in the global coordinate system after introducing the following identification: [K] = -1 (6) Equation (6) represents the mathematical definition of full anisotropy which is due to the rotation of the permeability tensor form local to global coordinates. The continuity equation for steady-state conditions in two dimensions is: Substituting equation (5) into (7) gives: ( ) is the angle of inclination of the major xaxis. Equation (8) represents the governing mixed partial differential equation of groundwater flow in the global coordinate system.

3-FINITE ELEMENT FORMULATION
The Galerkin technique is used to determine approximate solutions of equation (8) (2) is subdivided into triangular elements with denser elements around the drain as required by the principals of electrical analogy. For each finite element mesh, it is necessary to find the solution of the following equation: Trial solutions are chosen of the form given in equation (10): (10) Where (i =1, 2, 3... n) is a particular set of spatial functions known as shape functions that satisfy the boundary conditions. These functions are defined piecewise element by element. In the Galerkin method, the weighting functions are set equal to the shape functions which define the trial solutions. On the other hand ( ) are the undetermined coefficients of the (n) nodal points. Hence the expression in equation (9) is set orthogonal to all the functions of the system ( ). This is expressed as: (11) If the shape functions ( are selected in such a way that they have a value of unity at node (i) and zero at all other nodes in the domain (Ω), then the coefficients ( ) ( i = 1, 2, ... , n) are equal to the required functions ( at the (n) nodes. Substituting equations (7); (9) and (10) into equation (11) (13) and (14) are done by Gaussian Quadrature. The system of linear equations is solved by Guass-Seidel iteration.

4-EXAMPLE PROBLEMS
The finite element model has been used to solve two different flow problems under full anisotropy as shown in figures (1) and (2). Figure (1) represents the flow of ponded water; case (a) whereas figure (2) represents the problem of seepage under steady recharge; case (b). In the first example, the influence of full anisotropy upon drain rates of flow is studied for different ratios of drain depths and for a given ratio of anisotropy. In the second example, the variation of phreatic levels midway between two drain pipes is studied under the conditions of full anisotropy and constant recharge for different ratios of permeability and drain depths. The boundary conditions of each case are as follows:

5-NUMERICAL RESULTS
The validity of the finite element model of figure (1) is checked with the analytical solution of Kirkham (1949) for the isotropic conditions and presented graphically in figure (3) as the black line whereas the analytical solution of Kirkham (1949) is represented by the black line of figure (3). The figure shows the variation of the drain discharges for the conditions of b = 0.6 m, r = 0.01 m, 2s = 1.2 m and ratio of anisotropy = 1.The behavior of the finite element solution is expected to be similar to Kirkham's analytical solution. Both solutions verified the influence of drain position on drain discharge rates. Drain discharges were found to increase when the drains were raised above the impervious base of the tank. Up to a certain height, drain flow rates started to decrease due to reduction in the available hydraulic head difference between the ponded soil surface and the drain center. The maximum deviation in the calculated (Q) values by the FE model as shown in the figure are not more than 6% of those of the analytical solution. The reason for this slight deviation is attributed to the density of the elements chosen around the drain pipe. Figure (4)  ). This is mathematically interpreted in terms of the influence of the (cos sin ) term of eq. (8) which is maximum at (45°) and minimum at any other angle. The figures also indicate that water table heights of isotropic soils are higher than those of anisotropic soils, a reason which is simply justified by letting ( ) in equation (8) and dropping the cross product term ( ) cos sin thus leaving equation (8) in a form identical to Laplace's equation since, for any desired ( ) the term ( + ) of equation (8) is always one.

7-CONCLUSIONS
For a given ( ) the higher the (RA) value the smaller will be the value of (H). However the figures show a sort of increasing fluctuations in the values of (H) with an increase in (RA) values. The reasons behind such fluctuations are due to the effect of the cross product of the permeability tensor ( ). The effect of the cross product of this tensor vanishes at ( ) = 0, 90 and 180 degrees of rotation and is distinguishable at any other value of ( ). The conclusions regarding the behavior of tensorial permeability and its influence upon DRAINAGE of soils are considered helpful in deriving an analytical solution to such complicated problem. x, y = global coordinate axes.
[R] = transformation (rotation) matrix. Nj = shape function. Q = drain discharge. P = recharge (rainfall) rate. H = water table height at mid spacing of drains. r = drain radius. b =depth of flow region. d = height of drain above datum level (x-axis). s = semi-spacing of drains.