Pure and Applied Mathematics Journal
Volume 5, Issue 1, February 2016, Pages: 15-22

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

Email address:

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

1. 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 . 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.

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; [2], [5]


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 [11];


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 , [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 kx and ky' respectively [6], 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 [15], Eqn. (5) can be written as


In which values of at time t are known and at time  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,) 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 [7], [2]. 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 [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;





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 [4] as:


in which r (P, Q) is the distance from a fixed point, P, within R, to another point, Q, and 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 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  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. (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 [6] 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 [5]:


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 [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].

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.

Fig. 1. Nodal Arrangement for Square Region.

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 [9] 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, [14]



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 [9] 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.

5. 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.


  1. ALI A Ameli (2014),"Semi-analytical methods for simulating the groundwater-surface water interface", in fulfillment of the thesis requirement for the degree of Doctor of Philosophy in Civil Engineering Waterloo, Ontario, Canada.
  2. Albrecht, J. and Collaty, L. (Eds.),(1980) "Numerical Treatment of Integral Equations", Birkhauser Verlag, Basel, 1980.
  3. Bear.Jacob,(1978).Dynamic of Fluid in Porous Media,Dover Books on Engineering,ISBN-13: 978-0486453552
  4. Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C., (1984),Boundary Elements Techniques. Springer-Verlag, Berlin, 1984.
  5. Dupuit, J. (1863). Estudes Théoriques et Pratiques sur le mouvement des Eaux dans les canaux découverts et à travers les terrains perméables (Second ed.). Paris: Dunod.
  6. Dubois, M. and Buysse, M.,(1980),"Transient Heat Transfer Analysis by the Boundary Integral Equation Method", in New Developments in Boundary Element Methods, Brebbia, Co (Ed.), C.M.L. Publications, Southampton, 1980.
  7. Goldberg, M.A. (Ed.),(1979), "Solution Methods for Integral Equations: Theory and Applications, Plenum Press, New York, 1979.
  8. James A,Ligett,Philip L-F.Liu,(1983) "The Boundary Integral Equation Method for Porous Media Flow",London,George Allen and UNWIN,1983.
  9. Lafe, O.E., Liggett, J.A. and Liu, P.L.F.(1981) "BIEM Solutions to Combinations of Leaky, Layered, Confined, Unconfined, Nonisotropic Aquifers ll, Water Resources Research, Vol. lJ, No.5, pp. 1431-1444, October, 1981.
  10. Lennon, G.P., Liu, P.L.F. and Liggett, J.A.(1979) "Boundary Integral Equation Solution to Axisymmetric Potential Flows; 1. Basic Formulation; 2. Recharge and Well Problems in Porous Media ll Water Resources Research, Vol. 15, No.5, pp. 1102-1115, October, 1979.
  11. Lennon, G.P., Liu, 1L.F. and Liggett, J.A.(1980) "Boundary Integral Solutions to Three-Dimensional Unconfined Darcy's Flow", Water Resources Research, Vol. 16, No.4, pp. 65 658, August, 1980.
  12. Massonnet, C.E. (1965) "Numerical Use of Integral Procedures in Stress Analysis", Zienkiewicz, O.C. and Holister, G.S. (Eds.), John Wiley and Sons, Inc., New York.
  13. M.N. Vu a,n, S.T. Nguyen a,c, M.H. Vu a,b,(2015)," Modeling of fluid flow through fractured porous media by a single boundary integral equation", Engineering Analysis with Boundary Elements 59 (2015) 166–171.
  14. Norrie, D.H. and de Vries, G.(1978) "An Introduction to Finite Element Analysis ll, Academic Press, New York, 1978.
  15. Phoolendra K. Mishra and Kristopher L. Kuhlman(2013),Unconfined Aquifer Flow Theory - from Dupuit to present,physics.flu.dyn
  16. Wrobel, L.C. and Brebbia, C.A.(1981), "A Formulation of the Boundary Element Method for Axisymmetric Transient Heat Conduction", International Journal of Heat and Mass Transfer, Vol. 24, No.5, pp. 843-850, 1981.

Article Tools
Follow on us
Science Publishing Group
NEW YORK, NY 10018
Tel: (001)347-688-8931