Asymptotic Expansions for Moench’s Integral Transform of Hydrology

Theis’ theory (1935), later improved by Hantush & Jacob (1955) and Moench (1971), is a technique designed to study the water level in aquifers. The key formula in this theory is a certain integral transform H[g](r, t) of the pumping function g that depends on the time t and the relative position r to the pumping point as well as on other physical parameters. Several analytic approximations of H[g](r, t) have been investigated in the literature that are valid and accurate in certain regions of r, t and the mentioned physical parameters. In this paper, the analysis of possible analytic approximations of H[g](r, t) is completed by investigating asymptotic expansions of H[g](r, t) in a region of the parameters that is of interest in practical situations, but that has not yet been investigated. Explicit and/or recursive algorithms for the computation of the coefficients of the expansions and estimates for the remainders are provided. Some numerical examples based on an actual physical experiment conducted by Layne-Western Company in 1953 illustrate the accuracy of the approximations.


Introduction
Pumping from a well, group of wells or any kind of water reserve has been the subject of a deep study by specialized geologists in the last century. One of the most important subjects of investigation in this area is the prediction of the water level in the aquifer due to pumping mechanisms and pumping techniques that have evolved over time due to continuous technological improvements. It would be an impossible mission to cite all the important contributions to this topic, most of which were published in specialized journals of geology. Without the aim of being exhaustive, a few of them that give an idea about the state of the art in the problem of predicting water levels in aquifers due to different pumping mechanisms are mentioned: [1][2][3][4][5][6][7][8][9]. In order to obtain a complete picture of the problem, a review of these references and references therein is suggested.
Until 1971, the mathematical models that describe the water level changes were based on Theis' theory [8] and the superimposition of the effects of each pumping well [7,9] for both leaky and nonleaky aquifers. The first precise mathematical model was proposed in [4], which is used to predict water level changes at a certain point of an aquifer when a variable pumping is acting from one or more wells. Their main contribution is the derivation of an integral representation of the drawdown of water in the aquifer in the form of a convolution integral, which is valid when the pumping is constant. Their theory is improved and generalized by [6] for an arbitrary pumping. The key formula that gives the response of the aquifer due to an arbitrary pumping is the convolution integral In this formula, H[g](r, t) represents the drawdown as a function of the distance r to the pumping point and of the time t elapsed since the beginning of the pumping, all of which are measured in natural unities. The parameter S is the (dimensionless) storage coefficient; T is the transmissivity coefficient (measured in units Length 2 /Time); P is the vertical permeability of the confining layer (measured in speed units Length/Time); m is the thickness of the confining layer (measured in length units Length); and g(t) is the pumping as a function of time (measured in units Length 3 /Time). The three parameters m, S and T are positive, and P is non-negative. Hantush and Jacob's formula [1,2] is the particular case of Moench's Formula (1) when the pumping is constant, g(t) = 1. On the other hand, as a difference with Theis' theory, Formula (1) assumes a negligible storage capacity of the confining layer and that the leakage rate through this layer is proportional to the drawdown in the aquifer. When P = 0, the aquifer is said leaky, and it is nonleaky when P = 0. Formula (1) gives a quite accurate description of the level changes in the aquifer. However, it has some limitations, since there are some physical situations, more or less relevant in practice, that are not taken into account in Moench's theory. They are described by Moench himself in [6]. One of them is, for example, the presence of boundaries in the aquifer that are not included in the theoretical analysis, although it could be handled by using image theory [6]. For a more detailed description of the limitations of Moench's theory, see [6].
A partially penetrating pumping well has an effect on the vertical components of the water flow. Thus, in the presence of partially penetration effects in a leaky confined aquifer, some terms must be added on the right-hand side of (1) [10,11]: terms that are a combination of elementary functions as well as Bessel functions. Other variations and alternative representations in terms of Bessel functions may be found in [5]. In the case of the presence of parallel unsteady-state flow, Vanderberg showed that the denominator in Formula (1) must be slightly modified u → u 3/2 or, in other words, that a constant pumping translates into g(t) = t −1/2 [12,13].
Therefore, from a mathematical point of view, the following integral transform of a function g for which the integral exists is of interest, with t > 0 or t = ∞ (in the case t = ∞, only the first integral above is considered).
The integral transform considered in Chapter 27 of [14] is the particular case t = ∞ of the above formula. Chapter 27 in [14] contains several asymptotic expansions in terms of Bessel functions as well as several illustrative examples. In addition, the particular case t = ∞ and g(u) = u p is a Krätzel integral, which was introduced in [15]. Apart from a constant factor, integral (1) is integral (2) with the identification x = P/(2m) and y = Sr 2 /(4T). Therefore, integral (2) is from now on called Moench's transform of the function g. Even for elementary pumping functions g(t), the right-hand side of (2) cannot be evaluated in terms of elementary functions. Many numerical techniques have been considered in the literature for the evaluation of the right-hand side of (2) and eventually implemented in computer codes and compared with experimental data obtained from several aquifer systems (see for example [2,6,9]). In this paper, analytic approximations of (2) are of interest. Several analytic approximations, convergent or asymptotic, have been proposed in the literature in several limits of the parameters x, y and t. According to the experiments described in [1][2][3][4][5][6][7][8][9] and from a practical point of view, it is of interest to evaluate S t [g](x, y) for a domain of the time variable t that is as large as possible: small and large values of the parameter x and moderate (positive) values of the parameter y. In addition, the most interesting pumping functions g(t) are power functions of the form g(t) = t ν , ν ∈ R.
Several asymptotic or convergent expansions, accurate for small values of x and y, have been derived in [9,16,17] for g(t) = t ν . A convergent expansion for large x and small y was derived in [18,19] for g(t) = 1. Asymptotic expansions for large t, unbounded values of y and x = 1, are given in [20]. A detailed summary of the known analytic approximations is given in [21]. To our knowledge, analytic approximations of S t [g](x, y) for large values of x, moderate values of y and valid for any value of the time variable t are not given in the literature. This region of the parameters is important in practical situations (see the experiment described in Section 4), and then, the objective of the present paper is the analytic approximation of S t [g](x, y) in this region. For the sake of generality, in principle, pumping functions g(t) are considered in a large functional space with special attention to the more interesting family g(t) = t ν . In the next section, the set of functions g(t) that we consider in this paper is specified.
In Section 2 an integral representation of S t [g](x, y) appropriate for the asymptotic analysis when xy is large is derived. Observe that when x is large and y bounded from below, the product xy is a large variable. In Section 3, the family of pumping functions g(u) = u ν , ν ∈ R is analyzed in more detail as they are more interesting from a practical point of view. An asymptotic expansion of S t [g](x, y) for large values of the product xy (large x with y bounded from below) and arbitrary values of t > 0 is obtained with explicit and recurrent expressions for the coefficients of the expansion. In Section 4, the accuracy of the approximations is analyzed by means of several numerical experiments based on an experiment carried out by the Layne-Western Company in Illinois in 1953. Some conclusions are pointed out in Section 5. A general pumping function g(u) is discussed in a separate appendix, deriving an asymptotic expansion of S t [g](x, y) for large values of xy and arbitrary values of t > 0. In this more general case, obtaining an analytic expression for the coefficients of the expansion is not possible in general, and the coefficients must be computed numerically or by using an algebraic manipulator.

Preliminaries
For an appropriate analysis of S t [g](x, y), when xy is large, a suitable integral representation of the Moench transform (2) must be found. With this aim, a sequence of changes of the integration variable is introduced below. After this sequence of transformations, the integral (2) is written in the form of a combination of Laplace transforms, which is much more appropriate for an asymptotic analysis.

•
First change of variable. The dilatation u → s defined by means of the formula u = y x s is introduced in (2), This change of variable allows the identification of the phase function, f (s) := s + s −1 and the asymptotically relevant point s = 1.
• Second change of variable. The change of variable s → u defined by means of the formula u = 1 2 (s + s −1 ) − 1 is considered (see Figure 1a), which has two inverse functions, The functions s ± (u) are depicted in Figure 1b. Observe that s + (u)s − (u) = 1. The minimum value of u(s) is attained at s = 1 and its value is u = 0.  Figure 1). On the other hand, when s * > 1, the position of s * and s −1 * is interchanged in Figure 1, and then, the integration variable u in (3) runs from u = ∞ to u = 0 (with s(u) = s − (u)) and also from u = 0 to u = u * (with s(u) = s + (u)). Then, when s * < 1, only s − (u) is required in the computation. However, when s * > 1, then both, s + (u) and s − (u) are involved. It is clear from Figure 1 that where • Third change of variable. The two integrals in the second and third lines of (4) are suitable for an asymptotic analysis when xy is large by using Watson's lemma ( [14], Chap. 2); a brief and simple version of Watson's lemma is summarized in Appendix A. This is so because the maximum value of the exponent in the integration interval is attained at u = 0, which is the left end point of the integration interval, and only the factor √ u in the denominator of these integrals vanishes at this point, whereas the other factor √ 2 + u is analytic there. The situation is different for the integral in the first line of (4). Once again, the maximum value of the exponent in the integration interval is attained at u = u * , which is the left end point of the integration interval. However, now, although the whole denominator u(2 + u) is analytic there, when the variable t approaches the critical value y/x (s * approaches the critical value 1), the integration limit u * approaches the value zero, where u(2 + u) vanishes. Therefore, in the two last integrals, an expansion of the integrand at the asymptotically relevant point u = 0 is possible for any value of t, and Watson's lemma can be applied uniformly in the time variable t. On the other hand, in the first integral, this expansion is only possible when the time variable t is bounded away from the critical time It is assumed in this discussion that as a function of the variable u, the "composite" pumping function is analytic at u = u * for any value of t > 0. This assumption does not mean any restriction from a practical point of view, as any realistic pumping function meets this condition. Nevertheless, the specific requirements for the pumping function will be established below.
As an asymptotic expansion of the Moench transform uniformly valid in t ∈ [0, ∞) is sought, the above-mentioned problem needs to be solved. The solution is a third change of variable u → s with u = s 2 . In order to simplify the notation, the following functions are defined and the parameter Then, introducing the change of variable u = s 2 in (4), The asymptotic features of the three integrals have not changed, but now, the factor 1/ √ 2 + s 2 (in the three integrals) is analytic at the asymptotic relevant point s = 0 (in the two last integrals) and s = Λ (in the first one), even when t → y/x and Λ → 0. However, before proceeding further and apply Watson's lemma to these integrals, a more compact form of S t [g](x, y) that simplifies computations is considered, where we have introduced the step function χ(t) and the sign function χ: Either of the two integrals in (7) has the form of a Laplace transform, which is suitable for an asymptotic analysis when xy is large by using Watson's lemma ( [14], Chap. 2) for any value of t ≥ 0, whenever the asymptotic expansion of g( y/x α ± (s)) at s = 0 is uniformly valid for large x and bounded values of y. The two integrals in (7) shall be studied separately. However, for the sake of clarity, and because it is more common in practice [9], in the next section, the particular case of a pumping function g(u) of power type is considered: g(u) = u ν , and the study of a general pumping function g(u) is relegated to the Appendix B. As it will be seen below, in the case of a pumping function of power type, the coefficients of the asymptotic expansion of g( y/x α ± (s)) may be computed explicitly and do not depend on the parameters x, y (and then, it is trivially uniform in these variables), whereas in the case of a general pumping function, the asymptotic expansion of g( y/x α ± (s)) must be written in terms of arbitrary coefficients, and their dependence on the parameters x, y must be analyzed.

Asymptotic Approximation of (7) for a Pumping Function of Power Type
Setting g(u) = u ν in (7), the Moench transform of this pumping function reads where with α ± (s) given in (5). From [14], Watson's lemma can be applied to the two integrals above whenever the two functions h ± (s) are infinitely differentiable at s = 0, and their derivatives are uniformly bounded by a positive constant times an exponential function e bs 2 , with b < 2 √ xy. From the definition (5) of the functions α ± (s), it is clear that both functions h ± (s) are infinitely differentiable at s = 0, and their derivatives are power functions bounded by a positive constant times an exponential e bs 2 for any b > 0. In order to apply Watson's lemma, the asymptotic expansion of h ± (s) needs to be considered at s = 0 ( [14], Chap. 2), The coefficients c ± k may be computed by using a symbolic manipulator. They may also be computed by means of a recurrence relation. In order to derive this recurrence, the above expansion must be introduced into the following differential equation satisfied by h ± (s), After some simplifications, and equating the coefficients of equal powers of s, it can be checked that the coefficients c n satisfy the recurrence relation The first few coefficients of the expansion are From Watson's lemma ( [14], Chap. 2), an asymptotic expansion of (8) for large xy follows by introducing the expansion (10) into the two integrals on the right-hand side of (8) and interchanging the sum and integral. For the first integral, the asymptotic expansion is given by and for the second one, with The integrals in the asymptotic expansion (11) may be computed exactly in terms of elementary functions, For the values of z required in the above formula, we have that On the other hand, the computation of the integrals (13) in the asymptotic expansion (12) is more elaborated. Integrating by parts in the integral (13), it is straightforward to obtain the recurrence relation with where erf(z) is the well-known error function [23] that models fast but continuous transitions between two constant values. The solution to the recurrence relation (15) is, for n = 0, 1, 2, 3, . . ., where empty sums must be understood as zero.
Therefore, from (8), (11), (12), (14) and (16), the following asymptotic expansion of the Moench transform for pumping functions of power type is finally derived, where for n = 0, 1, 2, . . ., the sequence of functions Ψ n (x, y, t) is defined in the form where empty sums (for n = 0) must be understood as zero.
Observe from (18) that the functions Ψ n (x, y, t) depend on the variables x, y, t through the combination Λ 2 √ xy. Moreover, the right-hand side of (18) is a sum of constants, an error function (whose absolute value is bounded by 1) and terms of the form (Λ 2 √ xy) k e −2Λ 2 √ xy , all of them bounded for any value of Λ 2 √ xy ∈ [0, ∞). Therefore, the numerators Ψ n (x, y, t) in terms of the expansion (17) are bounded functions of Λ 2 √ xy, which means that the terms of the expansion (17) are of the order O(( In particular, the first-order approximation (obtained by considering only the first term n = 0 of the sum (17)) is given by the following formula: The accuracy of Formula (19) is illustrated in Figure 2. The argument of the error function on the right-hand side of (19) (and also in all the error functions in the successive terms of the expansion (17)) is 2Λ √ xy, where Λ is defined in (6). At the critical value t = t * = y/x, we have that Λ = 0, and then this error function, multiplied by the sign function, presents a fast transition from the value −1 to the value 1 as t crosses this critical value t * . Then, this function encodes the fast (but smooth) transition that Moench's transform experiments at the critical time t = t * , from a value close to the zero value (as S 0 [u ν ](x, y) = 0) to a value close to It is shown in Figure 2.

Numerical Experiments
In Figure 2, the accuracy of Formula (19) (particular case of Formula (17)) for certain selected values of the parameters x, y, ν is shown. In this section, one step forward is given, and the accuracy of the approximation (17) in a more realistic situation is shown. Consider the experimental data obtained from the experiment described in Part 3 of [9] that was conducted by Layne-Western Company on 2 July 1953 on a village well at Gridley, McLean County, Illinois (for more details, see Part 3 in [9]). For this experiment, the parameter values are T = 10,100 gpd/ft = 0.9378 feet 2 /min, S = 0.00002, r = 824 feet, P/m = 560/18 gpd/feet 3 = 2.8886 × 10 −3 min −1 , and a pumping function g(u) = 220 gpm = 29.414 feet 3 /min. These parameter values translate into the values (x, y) = (P/(2m), Sr 2 /(4T)) = (144.429, 3.6202) used in Moench's transform S t [g](x, y).

Conclusions
A complete asymptotic expansion of Moench's transform (2) for large values of the parameter x with y bounded from below (or large values of xy) that is uniformly valid in the time variable t ∈ [0, ∞) is derived. It is given in (17) for the particular case (and more interesting in practice) of a power-type pumping function. For a more general pumping function, a complete asymptotic expansion of the Moench transform is given in Appendix B in Formula (A6). In any case, the expansion is given in terms of elementary functions and an error function with fixed argument 2Λ √ xy. For a power-type pumping function, Moench's transform, which measures the water level in the aquifer, experiences a fast (but continuous) transition between a low level and a top limit level: more precisely, from the zero value attained at t = 0 (S 0 [u ν ](x, y) = 0) to the limit value attained at t = ∞, In the case of a more general pumping function g(u), the transition is similar but with the factor y x ν 2 replaced by 1 2 g y x . This transition occurs around the critical time t * = y/x, where the Moench transform takes half of its maximum value, S t * [g](x, y) = 1 2 S ∞ [g](x, y). This fast but smooth transition is encoded in the error function that appears in the expansions (17) or (A6) at any order n of the approximation, as all the terms Ψ n (x, y, z) on the right-hand side of (17) or (A6) contain the error function. The argument of this error function is (see (18) and (6)) that vanishes at t = t * = y/x, and precisely the error function experiences the abovementioned transition when its argument vanishes. Moreover, the larger x is, the faster the transition is, as small variations of the time t around the critical value t * are amplified by the factor x.
The approximations (17) or (A6) have been tested by means of numerical examples that consider realistic parameter values, corresponding to actual physical experiments (see Section 4). As it can be observed in Figures 2 and A1 and Tables 1-4, the accuracy of the approximations is remarkable.
We have analyzed Moench's transform for large values of the parameter x with y bounded from below and t ∈ [0, ∞). It would be interesting to find other approximations of Moench's transform valid in new regions of the parameters, such as small values of x and/or y. Moreover, it would be interesting to find approximations uniformly valid for x and/or y in [0, ∞). This is the subject of future research.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Watson's Lemma
Watson's lemma ( [14], Chap. 2) is one of the simplest yet most powerful results in the theory of asymptotic expansions of integrals. It can be directly applied to Laplacetype integrals such as the ones analyzed in this paper. Basically, it consists of termwise integration of a series expansion of a factor of the integrand. More precisely: Theorem A1 (Watson's lemma). Consider the integral where x is a large positive variable and g : [0, ∞) → R is a smooth enough function. We assume that the integral exists for a large enough x, and g has an asymptotic expansion at t = 0 of the form Then, the integral F(x) admits the following asymptotic expansion, for large x, A formal derivation of this formula follows by introducing the asymptotic expansion (A2) of g(t) into the integral (A1) and interchanging summation and integration.

Appendix B. Asymptotic Approximation of (7) for a General Pumping Function
In this section, a general pumping function g(u) in (7) is considered. Then, where with α ± (s) given in (5). From ( [14], Chap. 2), Watson's lemma may be applied to these integrals whenever the functions h ± (s) are infinitely differentiable at s = 0, and their derivatives are uniformly bounded by a constant times an exponential function e bs 2 , with b < 2 √ xy. Therefore, as α ± (0) = 1, the pumping function g(u) is required to be infinitely differentiable at u = t * = y/x; that is, the pumping function must be infinitely differentiable at the critical time where the fast transition between the low and top levels of the aquifer occurs. Taking into account the definition (5) of the functions α ± (s), when g(u) is infinitely differentiable at u = t * , the composite function g y x α ± (s) is infinitely differentiable at s = 0. The pumping function g(u) and its derivatives are also required to be bounded by an exponential function. In order to apply Watson's lemma ( [14], Chap. 2), the Taylor coefficients c ± k of the function h ± (s) at s = 0 need to be computed, As a difference with respect to the case of a pumping function of power type, in this more general case, a recurrence relation for the coefficients c ± k cannot be provided. Nevertheless, in practical applications, only the first few coefficients c ± k are needed, and they may be easily computed by hand or by means of a symbolic manipulator.
In the case of a general pumping function g(u), the asymptotic analysis of Formula (A6) is a little bit more elaborated than in the case of the power-type pumping function analyzed in Section 3. In that case, the function h ± (s) given in (9) does not depend on the variables x, y, and then, the coefficients c ± k of its asymptotic expansion at s = 0 do not depend on the variables x, y, either (see Formula (10)). Using this fact, it is shown in Section 3 that the functions Ψ n (x, y, t) in (17) are bounded functions (and then, in particular, they are of the order O(1) when xy → ∞).
However, now, the function h ± (s) given in (A4) depends on the variables x, y through the combination y/x, and then, the coefficients c ± k in (A5) of its asymptotic expansion at s = 0 also depend on the variables x, y through their quotient y/x. The precise value of these coefficients depends, of course, on the actual pumping function g(u). However, using Faá di Bruno's formula for the derivative of a composite function [24] and the Leibnitz formula, regardless of what pumping functions are being considered (whenever it is infinitely differentiable at u = t * ), these coefficients have the following form where σ n,k are combinations of the derivatives of order 0, 1, 2, . . . , n of the function α ± (s) evaluated at s = 0. The precise value of these numbers is irrelevant; it is only necessary to note that they are independent of the variables x, y. It is clear from (A7) that the coefficients c ± k are bounded functions of x, y whenever x is bounded from below, and the successive derivatives of the pumping function evaluated at the critical time t * are bounded. Therefore, the functions Ψ n (x, y, t) in (A6) are, as well as they were in (17), bounded functions of x, y whenever the pumping function and its successive derivatives, evaluated at the critical time t * , are bounded.
It can be deduced from this formula that for large xy, only the value of the pumping function g(u) at the critical time u = t * is relevant to compute the Moench transform S t [g](x, y). The accuracy of this formula is illustrated with the following example. Consider a pumping function of rational type, g(u) = 1 1 + u , that models an extraction that starts at a constant rate and decreases like u −1 for a long time. The approximation is given by (A6), where the functions Ψ n (x, y, t) are defined in (18)  Observe that the coefficients c ± n of the above expansion h ± (s) ∼ c ± 0 + c ± 1 s + c ± 2 s 2 + c ± 3 s 3 + · · · are of the form specified in (A7). They are products of the successive derivatives g (n) (u) = (−1) n n!(1 + u) −n−1 , n = 0, 1, 2, . . ., evaluated at u = t * = y/x, by polynomials in the variable y/x of degree n.