Boundary Integral Equation Method for Unsteady Two Dimensional Flow in Unconfined Aquifer
Azhari Ahmad Abdalla
Department of Mathematics, Prep-Year, University of Hail, Hail, Saudi Arabia
To cite this article:
Azhari Ahmad Abdalla. Boundary Integral Equation Method for Unsteady Two Dimensional Flow in Unconfined Aquifer. Pure and Applied Mathematics Journal. Vol. 5, No. 1, 2016, pp. 15-22. doi: 10.11648/j.pamj.20160501.13
Abstract: 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.
Keywords: BIEM, Unconfined Aquifer, Dupuit Assumption, Groundwater Flow
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 . 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, , . 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 , . 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, , , .
Another more direct method of solving time-dependent problems is with integral equations, in which the time derivative is written in finite-difference form, , , . 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 , . 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 , , 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 .
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.
2. 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; , 
in which S is the storage coefficient or effective porosity, Q. is the flow rate in well i located at (xi, yi), and the rate of recharge, RI, is positive for evaporation and negative for replenishment, is the piezometric head. The transmissivity, T, is given by ;
in which ZO (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.(2), but the aquifer thickness, B, is constant at any point. Consequently, Eqn. (1) is nonlinear in , while Eqn. (3) is linear. However, if Eqn. (1) is linearized by assuming that T is independent of changes in , , , 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  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 kx and ky' respectively , and if there is no recharge, Eqn. (1) can be written as
In which and If a backward-difference approximation is adopted for the time derivative as in , Eqn. (5) can be written as
In which,) is an average value of over the region. From Eqn. (6), it, can be shown that;
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 and then solving Eqn. (7) as in , . The estimate of 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 , 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 .
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,. A fundamental solution, u, of Eqn. (10) is given in  as:
in which r (P, Q) is the distance from a fixed point, P, within R, to another point, Q, and 0 , 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 remained constant over a sufficiently small time step and by performing the time integrations stepwise as in . . 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 , 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  and , that the grid spacing used to calculate surface integrals can be greater than that used in a finite-element solution . In , 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  by calculating the contour integrals, assuming that V and 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 , 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. (16) may be obtained by dividing the time interval into N time steps and reversing the order of integration to give;
In which tl = 0 and tN = 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 Vo is the initial condition at t = 0, given by Eqn. (10), then;
Which must be solved subject to the following boundary conditions Vi
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  to analyze problems of thermal shock by assuming that V and remain constant throughout each time step. Such a solution can be used as a starting point for the solution of the linearized groundwater flow equation (4), for when R' = 0. A second formulation, in which and vary linearly with time, will be used in an equation analogous to Eqn. (19) and the results compared with those obtained when and are assumed constant over each time step.
3. Solution of Homogenous Unsteady Unconfined Aquifer Flow Equation without Well
Homogenous regions without well can be modeled by :
in which is the interior angle between the boundary tangents at (x, y) and a2 is the ratio of transmissivity to storage coefficient i.e.
An algebraic approximation to Eqn. (24), can be written of M nodes, which gives M equations containing 2m unknowns ( and at each node). The remaining M equations required for a complete solution are obtained from the boundary conditions. To solve equation (24) two methods was used.
3.1. 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 that both h an remain constant at each node throughout each time step. The order of integration in Eqn (24) can be changed to give;
In which G is as defined as
And F is as defined as
The right hand side of Eqn. (25) contains only known information from previous time steps, while the left hand side contains the unknown values of and at time t = tN. Eqn. (25) is applied successively to each boundary node, and this together with the boundary conditions, gives a set of simultaneous equations which can be solved by Gauss Elimination to give the solution for and at each boundary node. The first integral on the left hand side of Eqn. (25) is calculated numerically as in , 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 , while the integral on the right hand side can be calculated by the numerical schemes using the integral formulations given in .
3.2. 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 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 at each node for r=0, ,.
The solution for at any point (x, y) contained inside the region R, at time t = tN is given by
Once and 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 vary linearly between successive time steps improve accuracy too much.
4. 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 on the boundary contour between such a node and the sharp corner are constant and equal to the values at the node.
Suppose that three sides of the square region shown in fig. 1 above are impermeable and that =0 everywhere within the region at t'=0, as given in  the exact solution is;
in which .
Numerically calculated values of , both on the boundary contour and at selected points within the flow region itself, are compared with the values obtained from the analytical solution of Eqn. (31).
Because at x' =l. 0 is infinite at t' =0 I (Eq. (32).), there are inherent problems in using Method 2 of section 3.2, in which h' and are assumed to vary linearly with time over each time step. When calculating values of and , at the end of the first time step, values of h' and at t'=0 are required, as described by Eqn. (24).
Thus, increasing the number of steps, the computational cost rises rapidly. It should be noted that the computational time required by Method 1 is approximately eighty percent of that red by Method 2. when same time steps are used for both sets of calculations. Consequently, the large instability in at X' =0, and the minor instability in , that occur near x'=l.0, are to be expected initially. However an excellent feature of this method is that despite the initial instability, the numerical solution closely approximates the analytical solution after only a few time steps. Also, numerically calculated values of within the flow region are in good agreement with Eqn. (31), even at the end of the first time step. The problem associated with the infinite value of at t'= 0 does not arise with Method 1 of section 3.1 because values of h' and an at t'=0 are not required for the solution at the end of the first time step, as stated in Eqn. (24). However, because Method 1 assumes that and are constant throughout each time step, time steps must be closely spaced to ensure accurate results. Immediately following the instantaneous rise in h' time steps are closely spaced, but as time proceeds and the step length increases, there is a loss in accuracy. Although the results can be improved by decreasing the time step length, the computational cost becomes prohibitive. It should be noted that the computational time required by Method 1 is approximately eighty percent of that required by Method 2. The instability that occurs with Method 2 following instantaneous changes in hi is not really considered a problem for groundwater flows because changes in are never instantaneous in real situations. If the instantaneous rise at x'=l.0 is replaced with;
In which E is a positive constant, very rapid changes can be modelled by choosing a sufficiently large value of E. More gradual changes in , are modelled by substituting a smaller value of E in the right side of Eqn. (33) The exact solution to this problem, calculated from Eqn. (31) using the Duhamel super position integral, 
in which. The right side of Eqn (35) is equal to zero at t' =0. for all values of except for . Consequently there is no discontinuity in 3h' and the problems experienced with Method 2 at t'=0 are eliminated. An instantaneous rise in , at x'=l.0 is closely approximated by substituting a values of in the right side of Eqn. (33). The numerical and analytical values of on the boundary contour, and the values at points within the flow region are compared with those from the analytical solution. The solutions for at x'=l.0 are compared also with the analytical solution values. There is now excellent agreement between the analytical solutions and those obtained using Method 2 in section 3.2, and there is no evidence of any numerical instability. There is virtually no difference between the solutions obtained using Method 1 for an instantaneous rise in h' and for h' This is to be expected because the two analytical solutions are almost identical except for different values of at x'=l.0 at t'=0, which do not appear in the numerical calculations. The accuracy is limited by the length of the time step as outlined previously. More gradual changes in h' at x'=l.0 can be modeled by substituting a value of in the right side of Eqn. (33). because , varies more slowly than in the previous examples, larger time steps were chosen. Solutions for h' around the boundary contour and at selected points within the flow region are and solutions for at x'=l.0 were compared with the analytical solutions. The solutions obtained using Method 2 are in excellent agreement with the analytical solutions, while those obtained using Method 1 are somewhat less accurate due to the large time steps used. The three cases examined so far all have steady-state solutions which are approached as t', becomes infinite. However, if the boundary condition at x'=l.0 is a periodic function of time, flow within the square region is oscillatory. For example, when the boundary condition at X'= 1.0 is;
the exact solution calculated from Eqn. (31) as in  is;
in which . Because of the relatively large time steps selected, numerical calculations have been attempted only with Method 2, the more accurate of the schemes presented in section 2. Also, since no discontinuity in or occurs at t' =0 (see Eqns. (37 and 38)), Method 2 will be stable. Values of h' on the boundary contour and at selected points within the flow region, and values of at x'=1, are compared with those from the analytical solution. Once again all the results are in excellent agreement.
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.