Comparison of the fully penetrating well drawdown in leaky aquifers between finite and infinite radius of influence under steady-state pumping conditions

In the paper theoretical derivation of steady state groundwater well pumping from leaky aquifers with infinite and finite radius of influence are presented. Based on the extensive literature review following mainly Jacob and Hantush work equations were derived from the cylindrical Bessel partial differential equation and results expressed in the combination of modified Bessel functions of zero order of the first and the second kind (I0, K0). We have shown that equation for steady state well pumping in the infinite aquifer is infinite limit of Hantush integral. Mathematical characteristics of solutions for infinite and finite radius of well influence were combined in the way that they can be represented as relative and absolute differences of drawdowns of each model. In the case when available data do not allow us to make a decision on the type of the radius of influence of the pumping well, they can help us in the interpretation of various errors due to application of different analytical models of pumping test.


Introduction
The most important and reliable in situ methods for groundwater investigations are pumping tests. During them water is pumped from the well constructed in the soil or rock and groundwater drawdown in the surrounding of the pumping well is observed. Pumping tests are intended for determination of aquifer physical characteristics, aquifer water balance and groundwater chemical status. These are reasons why pumping tests and the theory which connected with them is central to the hydrogeological science. The theory of pumping test is well developed and very complex. Since the very beginning when Theis published first mathematical model for unsteady groundwater flow toward the well (Theis, 1935) many conceptual and mathematical models of pumping tests were developed (for the review see Batu, 1998;Lebbe, 1999;Yeh & Chang, 2013 and references there in).
The appearance of water in intergranular porosity is conceptualized with the aquifer model which usually consists of three main elements: saturated part, unsaturated part and hydrogeological barriers. Different combinations of these elements represent different hydrodynamical models of aquifers. In general, we are talking about two main hydrodynamical types of aquifers with transitions between them. The simplest one is confined aquifer where saturated part is confined between two impermeable hydrogeological barriers. The other main aquifer type is unconfined aquifer where groundwater level fluctuates depending on the recharge. In natural conditions aquifers are complex entities consisting of beds with very different geometries and hydraulic characteristics. Aquifers where several beds with different hydraulic characteristics and contrasts are present are often conceptualised as leaky aquifers. In real aquifers drawdown of groundwater level during the pumping test or full well operation is complex.
The very first study on leaky aquifers under steady-state groundwater flow was presented by De Glee (1930). Jacob (1946) extended the work on leaky aquifers by introducing transient effect of leakage. In his treatment the key assumption was that the vertical flow rate in the upper hydrogeological barrier defined as an aquitard is proportional to the drawdown distribution in the same bed. Later Hantush and Jacob (1955) and Hantush (1959) derived analytical solutions for unsteady-state groundwater flow in leaky aquifer for fully penetrating well of infinitely small diameter. In addition, Hantush (1960Hantush ( , 1964 and DeWiest (1965) assumed that the piezometric head in the aquitard overlying permeable part of the aquifer does not change during water withdrawal from the underlying pumped part. The validity of the assumption that Darcian groundwater flow in the permeable pumped part of the aquifer is horizontal and in the overlying aquitard is vertical was tested by Neuman and Witherspoon (1969). The errors introduced by this assumption are less than 5 % if the difference in permeability between confined bed and semiconfining beds is of at least two orders of magnitude. Herrea and Figueroa (1969) and Herrea (1970) presented a correspondence principle where only the storage in the confining layers was taken into consideration. Şen (2000) has widened leaky aquifer hydraulic theory to non-Darcian flow and latter this analysis was extended based on volumetric approach by Birpinar and Şen (2004). For recent review on leaky aquifer hydraulic see Yeh & Chang (2013).
During mathematical simulations of the drawdown the well radius of influence is very often represented as limiting factor in calculations. This parameter is difficult to determine in nature. It is also difficult to determine whether to use analytical model of finite or infinite radius of influence. Radius of influence is very often defined from empirical formulas such as Sichardt equations (Powers et al., 2007) or it's estimation is based on the expert judgement from the field study. From the theoretical and practical point of view it is interesting to observe differences in drawdown calculations between mathematical models which include finite or infinite radius of influence.
In the paper mathematical analysis of drawdown in leaky aquifer during pumping test under steady-state conditions is presented. The analysis for the pumping tests with fully penetrating well in leaky aquifers with finite and infinite radius of well influence is extended based on solutions of Jacob (1946) and Hantush (1960Hantush ( , 1964. The comparisons are represented based on the various ratios between the drawdown for each of the different radius types under assumption that all other physical characteristic and pumping rate are the same. Ratios between different drawdowns are interpreted with various types of differences that can be interpreted as error analysis. Theoretical concepts are illustrated with numerical simulation. Finally, theoretical and numerical results are discussed.

Conceptual model
If the aquifer is not perfectly confined with upper and lower impermeable hydrogeological barrier leakage to the central water yielding unit -confined bed -may occur through the underlain or overlain semiconfining layer or aquitard. Leaky aquifers being either single or part of multi-layered aquifer systems and the degree of leakage between beds may become significant depending on the thickness and hydraulic conductivity of the confined bed which gives the main part of the aquifer yielding water. During pumping water from the aquitard water is also extracted through the confining layer. The conceptual model of the leaky aquifer is shown in the fig. 1. In parallel to this model several other leaky aquifer conceptual models are available in the literature (Yeh & Chang, 2013) but are not taken in consideration.
The leakage rates from the semiconfining layer -aquitard may be significant depending on hydraulic gradients around the pumping well. In the mathematical model of the pumping test from the leaky aquifer the thickness of the saturated part b' and vertical hydraulic conductivity of the aquitard K z p are taken into the account. It is hypothesised that leakage of water from the aquitard is strictly vertical and that no storage in this bed is present. The latter condition means that change of piezometric potential in the aquifer is simultaneous to change in the confined part of the aquifer. This part has transmissivity T that is defined as T=Kb; a product of hydraulic conductivity K and the thickness b of the confined part of the aquifer. Storage coefficient S of this defines vertical elastic properties of the aquifer. As a consequence of pumping from the aquifer the drawdown s appears on starting piezometric head h 0 that is horizontal at any distance r from the well. The drawdown s at r at time t is defined as s(r,t)=h 0 -h(r,t). Radius of influence of groundwater pumping R is defined as the distance from the well where s(R,t)=0. Under steady-state pumping conditions at any time two models of radius of influence of groundwater pumping R can be defined. In the first mathematical model of the leaky aquifer R is finite and constant; R=const.
In the second mathematical model of the leaky aquifer R is infinite; R=∞.
Together with boundary and initial conditions already presented in the fig. 1 the following assumptions are applied in the mathematical model (Batu, 1998): the confined part of the aquifer is homogenous and isotropic, the extraction rate of the well Q is constant, the aquifer is horizontal, has constant thickness b and is overlain by an aquitard with constant vertical hydraulic conductivity K z p and constant thickness of the saturated part b', the well penetrates entirely confined part of the aquifer, the diameter of the well is infinitesimally small with no storage and the groundwater flow in the confined part of the aquifer is horizontal.

Governing equation
Basic governing equation of the leaky aquifer in the vertical plane of the x, y, z Cartesian coordinate system is defined as (for derivation see Miletić & Heinrich Miletić, 1981)  where is defined as leakage factor. In the cylindrical system of the coordinates r and φ the equation (1) can be rewritten as (Hantush, 1964) and in the case of homogenous and isotropic aquifer according to ϕ Equation (4) can be recognized as modified Bessel equation of zero order (Lebedeev, 1970). For parameters and other symbols see fig. 1.

Basic solution
In the solution of governing equation (4) following Hantush (1964) valid boundary and initial conditions are defined as: In (5) last condition is consequence of Darcy law. Final solution of the governing equation (Hantush, 1964) is where u is defined as Integral (6) is known also as Hantush integral. It is important in many fields of mathematical physics and hydrology (Harris, 1997(Harris, , 2001Prodanoff et al., 2006). Detailed derivation of basic solution of leaky aquifer partial differential equation is given elsewhere (Hantush 1964;Batu 1998;Lebbe, 1999).

Steady state-solution in infinite leaky aquifer
In real aquifers steady-state conditions are reached only after longer time t. If we suppose that t →∞ than the limit of u is defined as and at sufficiently large time t u is small enough that u≈0. Consequently, integration borders become from 0 to ∞ and (6) becomes Based on Gradshteyn and Ryzhik (1994;equation 3.471.9) From (9) and (10) also follows and K 0 is modified Bessel function of second kind of zero order. After short manipulation in (10) and (11) it can be shown that drawdown for the infinite leaky aquifer s I in steady-state conditions is Which is the same results as de Glee (1930) defined initially through another derivation of equations.

Steady state-solution in finite leaky aquifer
General solution of modified Bessel equation of zero order (4) is (Lebedev, 1972) where are C 1 ,C 2 -constants I 0 -modified Bessel function of first kind of zero order.
Boundary conditions are defined as (Jacob, 1946) Determination of C 1 and C 2 of (13) from (14) leads to the equation of drawdown s F in finite leaky aquifer. Constants are determined in the area where r ≤ R. Elaboration of constants is not simple and straightforward, it is based on derivatives of s and limit properties of K 0 (x) and I 0 (x) (2) " " + 1 − " = (4) functions. Derivation of C 1 and C 2 is first given in Jacob (1946) and thoroughly summarised and elaborated in Batu (1998) and Miletić and Heinrich-Miletić (1981). Presentation of this derivation is out of the paper's scope. After definition of constants it follows:

Definitions
In hydrogeology we are frequently encountering problem of choosing a proper aquifer conceptual model. Dealing with results of pumping test it is a question whether to use finite or infinite model of well radius of influence. Available geological data are often not detailed enough or some information are missing for choosing the proper conceptual model. In such situation for calculations several models are used and their results are compared with the actual field measurements.
In the engineering practice measurements and calculated values are often expressed together with certain errors, among them are relative difference ε r or absolute difference ε a . The expression of those help us to understand the reliability of predictions preformed based on measurements and differences among them in their application of the mathematical models. These concepts can be used also in the comparison between drawdown calculations in leaky aquifers with finite and infinite radius of well influence under steady state pumping conditions.
We can define following differences and quotients. If x is any quatitative measure absolute difference ε a is defined as the absolute difference between two measured values x 1 and x 2 If reference value x ref is present absolute difference ε a ' is defined as where x is any value from the model. The generalized relative difference ε r is defined as Alternatively, average relative difference ε r avg can be defined as If the reference value x ref is defined than relative difference ε r ' is Sometimes ε r ' is defined as (12) and (15) are bearing some similarities. They can be easily used for the comparison between the modelled drawdown s in the aquifer with finite radius of influence -s F with the aquifer with infinite radius of influence -s I under steady state conditions and the same pumping rate Q.

From mathematical point of view equations
After short manipulation it can be shown from (12) and (15) that From that point right hand part of the (22) in the brackets can be understood as factor which is correcting infinite aquifer drawdown s I to the finite aquifer drawdown s F . Based on this we can define correction factor c F and consequently c F is independent of Q and depends only on B and R which are geometrical and physical characteristics of the leaky aquifer. Due to aquifer physical characteristics relation s I ≥ c F is always present and due to the characteristics of modified Bessel functions K 0 (x) and I 0 (x) functional relation s I ≥ s F is always valid. Consequently, head in the leaky aquifer with the same hydraulic characteristics under the same pumping rate Q is higher in the aquifer with finite radius of influence than in the aquifer with infinite radius of influence. It can be illustrated that under some circumstance heads around the well in both cases can be nearly equal.
We can further elaborate relations by dividing equation (15) with (12) and gaining or From the properties of I 0 (x) and K 0 (x) in (25) and (26) (28) follows from (21).
Equation (28) is not just a mere mathematical expression. It explains relation between two drawdown curves. If we have two leaky aquifers; with infinite and finite radius of influence under the same pumping rate and the same hydraulic characteristics ε r '(s I ) explains relative difference between both drawdown curves. Depending on r in the interval 0< ε r '(s I )<1 ratio explains relative differences between both drawdown curves. If the ε r '(s I ) is close or equal to 1 the curves have the same spatial distribution, and if ε r '(s I )=0 drawdown curve of s I is beyond the radius of influence R of the finite leaky aquifer.
With the analogies of equations from (16) to (19) and according to the definition in (16) subtracting (15) from (12) following expression for ε a can be derived Absolute difference ε a is expressed in length units. Comparing to ε r '(s I ) in (28) which is independent on pumping rate Q absolute difference ε a depends on it. Relation between ε a and ε r '(s I ) is also obvious.
Generalized relative difference ε r from (18) and with the help of (12) and (15) can be defined as Generalized relative difference ε r is also dimensionless quantity. It can be applied for analysis when no preference to s I or s F are given. This is the case when we are not sure if model of finite or infinite radius of influence is valid and we want to keep both results. Consequently, average relative difference ε r avg followed from definitions above is defined as

Results and discussion
In the following chapter we are presenting numerical simulation results based on the previous mathematical theory. Simulations were performed with build in numerical functions of modified Bessel functions of the first kind I 0 and the second kind K 0 in Excel for Mac 16.16.1. and with the program for symbolic and numerical computation Mathematica for Mac version 11.3.0. Numerical results are discussed from the hydrogeological point of view.

Estimation of leakage factor B
Main physical parameter in simulations is leakage factor B defined in (2) which is combination of two other physical parameters and one variable which in our simulation can be considered as a constant. Those parameters are: transmissivity T of the confined unit and K z p which is vertical permeability of semi-confining layer while b' defines head in the later. Based on the expert judgement of T, K z p and b' we have estimated values of B. For simulations b' = 2 m was used. As expected in the real aquifers T was considered on the interval between 5⋅10 -2 m 2 /s and 10 -5 m 2 /s and K z p was considered in the interval from 10 -9 m/s to 10 -6 m/s. Calculated values of B are represented on double logaritmic scale in fig. 2.
In the range of applied T highest values of B are present at K z p = 10 -9 m/s. In this case B values are calculated in the interval between 173 m and 12,247 m. Lowest B values are present at K z p = 10 -6 m/s. In this case B values are calculated in the interval between 5.5 m and 387 m. Therefore, total simulated range of B is from 5.5 m to 12,247 m. From the simulated values we can see that B is influenced by the K z p . In leaky aqufers semi-confining layer with vertical permeability K z p = 10 -6 m/s due to the depression in confined layer it is highly unlikely that vertical flow will appear. Consequently, values of B in real aquifers tend to be in the higher part of the interval.
Estimated values of B can be considered also in the evaluation of ratio R/B which is important in presentation of simulation results on the relative scale r/R. Expected radius of influence R in real aquifers under the steady state conditions are in the range of 500 m to 20,000 m. Consequently, expected approximate range of R/B is in the interval from 0.005 to 20.

Simulation of relative difference ε r '(s I )
To illustrate behaviour of ε r '(s I ) given by equation (28) we have chosen aquifer with influence radius R of 5,000 m. Such radius of influence can be expected in many natural aquifers. Results of calculations are presented in the fig. 3 for leakage factors B from 50 m to 20,000 m. At relatively small values of B large part of the curve is flatter reflecting ε r '(s I )=0. At higher r curve sharply turns up to values near ε r '(s I )=1. With higher values of B curvature is becoming flatter and values of ε r '(s I ) are becoming to rise slowly. Curves below B=5,000 m which is the same value as chosen R are concave with higher B they become convex. For high B values and at lower r values ε r '(s I ) starts to rise quickly and then at middle values of r the curve flattens and become nearly linear.
By simple reasoning it can be shown from (28) that results for different radius of influence R can be presented on the relative scale r/R. Diagram presented in fig. 4 is valid for any R at the same ratio R/B between radius of influence and leakage factor. Shape of lines are the same as they are on the fig. 3 and therefore reflecting the same relations as they are in the diagram for exact radius of influence. The diagram in fig. 4 can be understood as scaled diagram. Similarly, as before at relatively small values of B large part of the curve is equal to ε r '(s I )=0 and then at the right side of the diagram the curve sharply turns up to values near ε r '(s I )=1. From the diagram we can observe that for values of R/B < 1 the curves are concave and for values R/B > 1 the curves turn to be convex.
Curves fluence becomes important with larger B values. Therefore, when semiconfining layer has several orders of magnitude lower vertical permeability than in confined layer it become important which analytical model for radius of influence is used. Differences become bigger close to the well comparing to higher B values where differences are important far away from the pumping well.

Simulation of generalized relative difference
Generalized relative differences are presented only for relative distance r/R and are given in the fig. 5. Shape of the lines for different ratios R/B are nearly similar to the lines in fig. 4 where relative difference is shown. They are more convex comparing to relative difference ε r '(s I ).

Comparison between relative differences
Lines in fig. 4 and fig. 5 have similar shape therefore one may ask question whether there is any difference between the lines. For the comparison we have calculated both differences for two different ratios R/B = 1 and R/B = 0.25 respectively. Results are shown on the fig. 6. In spite of the similar shapes of the curves differences among curves exist. From the equations (28) and (30) it can be shown that for the same physical parameters R and B relative difference ε r '(s I ) is always larger than ε r . In the theoretical part of the paper we are not representing derivatives of (28, 30) but it can be illustrated that ε r '(s I ) always approach value of 1 (right part of the curve) slowly than ε r . Behaviour of ε r is the consequence of its definition in (18) where comparing to (20) value is weighted by s I .

Conclusions
In spite of the fact that in hydrogeological quantification of aquifers and their groundwater flow numerical models are widely applied, development of analytical mathematical models is still important. Analytical approach to groundwater flow enables different and deeper insight into the relations between different geometrical elements in the aquifers and their conceptual models. Analytical models are important for the control of numerical results and are very often applied as a scoping calculation representing first step in the consideration of hydraulic conditions in the aquifer.
In the paper we have presented classical derivation of the head distribution in the leaky aquifer under steady state pumping conditions. We have shown that infinite limit of Hantush integral which represents solution of the non-steady state pumping conditions in the leaky aquifer is solution for the steady state conditions. We have shown that solutions for steady state conditions under finite and infinite pumping well radius of influence are mathematically similar and that based on this characteristics comparison between them can be performed. They can be represented as relative and absolute differences of drawdowns for each model. In the case when available data do not allow us to make a decision on the type of the radius of influence of the pumping well, they can help us in the interpretation of various errors due to application of different analytical models of pumping test. We have shown that at larger leakage factors B determination of radius of influence R for the large part of the aquifer is not important, they become important at larger factors B when contrast between permeabilities in the semi-confining unit and confined unit becomes larger. Under such condition differences in drawdown are important in the vicinity of the pumping well. Fig. 6. Comparison of simulation between ε r '(s I ) and ε r with the same ratio R/B plotted on the relative scale r/R. For further consideration similar relation for non-steady solutions of the leaky aquifer are also interesting.