Boundary Integral Equation Method for Unsteady Two Dimensional Flow in Unconfined Aquifer

In this paper, we employed the use of Boundary Integral Equation Method to obtain numerical solutions of specific unconfined aquifer flow problems. Of the two formulations presented in this paper, that in which the piezometric head and its normal derivative are assumed to vary linearly with time over each time step has proved more accurate than that in which both piezometric head and its normal derivative remain constant at each node throughout each time step. Comparisons between Method 2 of section-3 and the analytical solutions have demonstrated the superior accuracy of the integral equation formulation.


Introduction
The study of groundwater flow is of great importance due to the rapid increasing for the demand of water. This has necessitated the need to explore new methods of exploiting water, especially water in aquifers, one of the main tasks in studies concerning aquifers flow is the determination of the distribution of the piezometric head and its normal derivative n φ ∂ ∂ . In order to find the distribution of piezometric head, φ , throughout an aquifer it is necessary to solve a partial differential equation subject to specified initial and boundary conditions, [1], [4]. It is desirable to calculate analytical solutions to these problems, but this is seldom possible due, often, to irregularities in the shape of aquifer boundaries. Consequently, approximate solutions must be calculated either by analogue methods, which are based on the similarity between groundwater flow and other physical systems obeying the same partial differential equation, or by numerical techniques [4], [7]. Many numerical techniques been employed to solve both steady and unsteady groundwater flow problems. Probably the simplest and most widely used of these is the finite element method and the finite difference method, [14], [3], [8].
Another more direct method of solving time-dependent problems is with integral equations, in which the time derivative is written in finite-difference form, [11], [12], [13].
In this way, the solution at one time can be calculated from the solution at a previous time in terms of an integral equation containing only boundary data and a surface integral over the entire solution domain [10], [16]. The principal advantage of integral equation techniques is that if only boundary data is used in the calculations the dimensionality of the problem is essentially reduced by one. Thus, the amount of data preparation is considerably less than that required by either the finite-difference or finite-element methods. Solutions within the flow region can readily be calculated from an integral equation solution, as it will be shown, and the advantage of reduced data preparation is still significant because only points at which the solution is required need to be specified.
The flow of groundwater is essentially a three dimensional problem, but in many real situations the horizontal dimension of an aquifer are much greater than the vertical dimensions. In such cases it is possible to accept the Dupuit approximation [5], [15], which assumes that flow is horizontal. Consequently, the problem is reduced to two dimensions in the horizontal plane. The relative scarcity of published work in groundwater literature concerning the use of the boundary element method (BEM) to solve the equation that describes the flow in unconfined aquifer, is mainly due to the numerical difficulties encountered in applying this method to resolve non-linearity. The accuracy of the BIEM depends on the size of nodes and on the type of the interpolation function [9].
In this paper, BIEM with the constant interpolation function and BIEM with the linear interpolation function were formulated and used for the solution of two dimensional, unsteady homogenous, unconfined aquifer flow (not including well). The results obtained by the two methods were investigated and a comparison between them and with the exact solution was held.

Integral Equation Techniques
The horizontal dimensions of many aquifer are orders of magnitude larger than vertical dimensions. In such cases, it is possible to make use of the Dupuit approximation, which assumes that velocities are horizontal and that piezometric heads are constant along a vertical line. Under these conditions, unconfined flows satisfy following partial differential equation; [2], [5] in which S is the storage coefficient or effective porosity, Q. is the flow rate in well i located at (x i , y i ), and the rate of recharge, R I , is positive for evaporation and negative for replenishment, φ is the piezometric head. The transmissivity, T, is given by [11]; in which Z O (x, y) is the elevation of the bottom boundary of the aquifer and B (x, y, t) is the saturated aquifer thickness (which depends on φ for unconfined flows).
If the Dupuit approximation is applied, confined flows satisfy Here S, is a function of the elasticity of the aquifer. The transmissivity, T, is given by Eqn. linearized by assuming that T is independent of changes in φ , [16], [2], then both confined and unconfined flows will obey a partial differential equation of the same form. Furthermore, if the aquifer is homogeneous and isotropic, both confined and unconfined flows are governed by: Where T and S are constants. Integral equation techniques was used in [11] to analyze problems with combinations of leaky, layered, confined, unconfined, and non isotropic aquifers under steady state conditions.
To illustrate the basic technique for unsteady flows, consider unconfined flow through a homogeneous, anisotropic, non-leaky region, R, that is bounded by a contour, Γ , and in which the principal axes of seepage are parallel to the x and y axes. The permeabilities for directions parallel to the x and y axes are given by k x and k y ' respectively [6], and if there is no recharge, Eqn. (1) can be written as In which * x x = and * * x y k y y k = If a backward-difference approximation is adopted for the time derivative as in [15], Eqn. (5) can be written as In which values of φ at time t are known and at time t t + ∆ are unknown In order to solve Eqn. (6) by integral equation techniques, it must be linearized [8], [14]. The following was suggested in [11]: In which, ( ) t t φ + ∆ ) is an average value of φ over the region. From Eqn. (6), it, can be shown that; Where: And r (P, Q) is the distance from a fixed point, P, on the boundary to another point, Q. The variable,θ , is the angle between the boundary tangents at P. The solution is obtained iteratively by assuming a value of ( ) t t φ + ∆ and then solving Eqn. (7) as in [7], [2]. The estimate of ( ) t t φ + ∆ is then revised and the process repeated until there is essentially no change in successive solutions. The obvious disadvantage of such a formulation is that because the final term on the right side of Eqn. (7) is non-zero, in general, integrations must be performed throughout the solution domain at each time step. Thus, although the method can be extended to solve problems' in leaky aquifers with recharge [12], it does not have the advantage of reducing the dimensionality of the problem. The use of time-dependent fundamental solutions to solve problems having the same form as Eqn. (4) was first presented and the transient heat conduction through an anisotropic medium was analyzed in [4]. The method has since been used by many, [2], [15] where the solutions to the following form of the heat conduction equation were presented; x y onΓ ∈ Γ (13) in which V is the temperature, C is the thermal diffusivity, and Eqn. (10) is a special form of Eqn. (4), Eqn. (11) is the initial condition which must be satisfied throughout the solution domain, R, and Eqns. (12) and (13) are the boundary conditions which must be satisfied on the respective portions of the boundary contour, in which r (P, Q) is the distance from a fixed point, P, within R, to another point, Q, and 0 t τ ≤ ≤ 0 [5], so (9) can be written as: If P is a point on the boundary contour; [ ] In which θ is the angle between the boundary tangents at p. The solution of equation (16) was calculated by assuming that V and V n ∂ ∂ remained constant over a sufficiently small time step and by performing the time integrations stepwise as in [4]. [16]. Using this assumption in Eqn. (15) and changing the order of integration gives; The time integrals on the right side of Eqn. (17) can be evaluated analytically, and it can be shown that; In which (r, n) is the angle between r (P, Q) and the unit normal at point Q.
The obvious disadvantage of using Eqn. (17) is that, as with the formulation as mentioned in [11], a surface integral must be calculated at the start of each time step, and thus negating, to some extent, the advantage of integral equations over the more traditional techniques. However, the methods do retain their accuracy, it was claimed in [11] and [16], that the grid spacing used to calculate surface integrals can be greater than that used in a finite-element solution [14]. In [4], both Time dependent fundamental solution, linear time interpolation functions in Eqn. (17) were used to solve axisymmetric problems, in an effort to improve accuracy. An alternative formulation for the solution of Eqn. (16) was presented in [11] by calculating the contour integrals, assuming that V and V n ∂ ∂ varied linearly along a straight line joining adjacent nodes and integrating the expression analytically around the boundary. The trapezium rule was used to calculate integrals in the time domain. It was concluded in [11], that the initial portions of the time integrations were critical and had to be calculated accurately to avoid the introduction of a consistent error into the solution. This formulation, also, required the evaluation of an integral over the solution domain. A solution of Eqn.
In which t l = 0 and t N = t. At first glance, it may appear that Eqn. (19), is not better than any of the previous formulations because of the surface integral that appears on the right side. However, because the governing partial differential equation, Eqn. (10), is linear with constant coefficients, the principal of superposition can be applied.
In which V o is the initial condition at t = 0, given by Eqn. (10), then; Which must be solved subject to the following boundary Although the amount of data storage and computational effort required is greater for this method than for the others in the field, the advantage of reduced data preparation is very significant. Eqn. (20) was used in [6] to analyze problems of thermal shock by assuming that V and remain constant

Solution of Homogenous Unsteady Unconfined Aquifer Flow Equation without Well
Homogenous regions without well can be modeled by [5]: in which θ is the interior angle between the boundary tangents at (x, y) and a 2 is the ratio of transmissivity to

Method (1) -BEIM with Constant Interpolation
The simplest way to integrate Eqn. (24) numerically is to divide the time interval (0, t) into N time steps and to assume In which G is as defined as 1 2 2 And F is as defined as  25) is calculated numerically as in [5], and the second integral is numerically approximated by Eqns. (19) and (20). The final term on the left hand side of Eqn. (24) was approximated as in [3], while the integral on the right hand side can be calculated by the numerical schemes using the integral formulations given in [5].

Method 2-BIEM with Linear Interpolation
A second, way to integrate Eqn. (24), is to divide the time interval (0, t) into N time steps as before, but to assume that both φ and n φ ∂ ∂ vary linearly with time through out each time step. Following the integration we obtain: By applying this equation to each node successively and using a quadrature formula such as the trapezium rule, a set of simultaneous equations is obtained which can be solved according to the boundary conditions for φ and n φ ∂ ∂ at each node for r=0, 0 r ≠ ,.
The solution for φ at any point (x, y) contained inside the region R, at time t = t N is given by Once φ and n φ ∂ ∂ have been determined at all boundary nodes by the methods given in. above, φ at any internal point can be calculated directly by numerically integrating the right hand side of Eqn. (29). The integration is greatly simplified because there are no singularities in the integrand.
The assumption that φ and n φ ∂ ∂ vary linearly between successive time steps improve accuracy too much.

Comparisons Between Method-1 and Method-2
The results obtained in the above section will be compared with known analytical solutions. Results can be expressed most simply by introducing the following dimensionless variables; In which L is a characteristic length in the problem. Although solutions can be obtained for a number of boundary contour shapes, one of the most commonly used will be the square because analytical solutions for flows in this region is readily calculated. As the boundary conditions on either side of a sharp corner may be different and because there is a discontinuity in geometry there, nodes are placed adjacent to, rather than exactly at the corners, as shown in the figure, that the distance between a sharp corner and the adjacent nodes should be approximately one quarter of the usual nodal spacing. The contribution that these corner nodes make to the numerical integrals are calculated by assuming that values of h and ' n φ ∂ ∂ on the boundary contour between such a node and the sharp corner are constant and equal to the values at the node. x'=1, are compared with those from the analytical solution.
Once again all the results are in excellent agreement.

Conclusion
The integral equation techniques presented here have been found to be satisfactory for solving unsteady, two dimensional unconfined aquifer flow problems. The use of time-dependent fundamental solutions has led to accurate integral equation formulations for solving the unsteady Dupuit equation in the homogeneous regions. Of the two formulations presented in this paper, that in which the piezometric head and its normal derivative are assumed to vary linearly with time over each time step has proved more accurate. Comparisons between Method 2 of section-3 and the analytical solutions have demonstrated the superior accuracy of the integral equation formulation.