## Abstract

The time-dependent flow of groundwater through an aquifer, generated from rest by pumping from a well close to a river, is calculated without making Dupuit's approximation. The governing equations reduce to the diffusion equation for the pressure head with mixed boundary conditions at the surface of the aquifer and at the base of the river. Using transform techniques, the problem is reduced to an infinite set of linear equations. The steady-state solution provides a guide to the numerical solution of the time-dependent problem. Results are presented for the flux from the river into the aquifer as well as for the flux from the aquifer into the river and out again. Explicit expressions for these fluxes are obtained in the case of rivers that are narrow compared with the aquifer depth. The steady-state flux is then much smaller than the transient flux. Results are significantly different from those obtained using the Dupuit approximation.

## 1. Introduction

Vast water withdrawals have dramatically changed local and regional water budgets of aquifers and streams. The fraction of the pumping rate supplied by depletion of an adjacent stream is important for water resource management, because its knowledge facilitates the adjudication of water rights for each well near the stream. It is now commonplace that aquifers play a dominant role in irrigation.

A detailed evaluation of stream depletion rates by numerical modelling requires an extensive database that includes climate and land use (Sophocleous 2005). Lacking this information, decisions are often based on an early analytical model by Jenkins (1968), which is regarded as the standard tool for water management and rights adjudication.

A common feature of subsequent computational modelling is the use of the Dupuit assumptions, namely, the physical parameters are constant in each zone and the vertical flow is negligible. The latter has led authors to construct two-dimensional semi-analytic solutions based on a plan view of the stream and aquifer. Hunt (1999), Zlotnik & Huang (1999) and Butler *et al*. (2001) introduced a streambed leakage parameter to replace the previously assumed fully penetrating stream. Butler *et al*. (2007) extended this idea to include leakage into the aquifer from a low-permeability underlying aquitard. A more detailed analysis of this flow is given by Zlotnik & Tartakovsky (2008).

A fully three-dimensional initial value problem, generated by the start of constant pumping from a point well, is formulated here and, at an early stage, simplified by integration in the streamwise direction. Thus, the contaminant capture aspects are discarded in favour of determining the flux ratio of the total depletions from aquifer and stream. The resulting two-dimensional calculation, based on an elevation view of the stream and aquifer, rejects the Dupuit assumption of negligible vertical flow and finds significant dependence on the well depth of the groundwater flow and associated flux ratio. The limit of zero vertical hydraulic conductivity is demonstrated to be mathematically singular.

## 2. Mathematical model

Consider inviscid fluid of uniform depth *B* in the presence of a shallow stream of width *W* that lies above the strip , with directed vertically downwards. The initial position of the overlying water table is assumed to be bounded by and the impermeable boundary is located at . A pumping well at the point , where and , operates at a rate *Q* after . In the aquifer with anisotropic hydraulic conductivity *K*_{h} and *K*_{v} in the horizontal and vertical directions, respectively, and specific storage coefficient *S*_{s}, the hydraulic head satisfies the groundwater flow equation (cf. Freeze & Cherry 1979)(2.1)subject to the initial condition,(2.2)and the boundary conditions,(2.3)(2.4)(2.5)(2.6)In (2.1) *S*_{s} is the volume of water that a unit volume of aquifer releases from storage due to a unit decline in hydraulic head . Equation (2.4) defines a spatially constant head interface whose roles include the determination of streambed/aquifer leakage. In (2.5), the macroscopic velocity of the moving water table at the ‘free surface’ is the specific yield *S*_{y} times the microscopic velocity in the aquifer. Note that, although *K*_{v}≪*K*_{h} in practice, the boundary-value problem does not allow the limit *K*_{v}→0 to be invoked. This feature crucially distinguishes this study from those cited in the introduction, all of which neglected *K*_{v} and considered only horizontal flow with stream and aquitard effects modelled by equivalent boundary conditions. (Fig. 1 of Zlotnik & Tartakovsky (2008) displays an apparently vertical *y*-axis that is actually horizontal.)

Equations (2.1) and (2.5) define the length scale *S*_{y}/*S*_{s} and the time scale . Accordingly, dimensionless variables, with origin in the centre plane of the river, are introduced by setting(2.7)Then the geometry indicated by (2.3)–(2.5) suggests that useful dimensionless parameters are(2.8)Thus, *w* is the physical river width/aquifer depth ratio while *D* and 2*a* are, respectively, the dimensionless depth and width, scaled according to the physical quantities in the diffusion equation and free-surface condition. In addition, by setting(2.9)the well location can be described by its fractional depth *Z* and horizontal displacement *X*(>0) river widths from the river. Note that (2.7) implies that the ratio of vertical and horizontal velocity scales is , which may not be negligible.

When (2.7) and (2.8) are substituted into the diffusion equation (2.1) and the boundary conditions (2.3)–(2.5), the dimensionless hydraulic head *h*(*x*, *y*, *z*, *t*) satisfies the dimensionless equation(2.10)and conditions(2.11)

(2.12)

(2.13)

(2.14)

The analysis is simplified when the dependence of *h* on *y* is of minor interest by then defining(2.15)whence (2.10) becomes(2.16)with subject to the same conditions as *h*. Thus, the two-dimensional analysis below regards the horizontal cross stream and vertical as the directions in which the flows of principal interest occur. The geometry is shown in figure 1.

The even and odd, with respect to the centre plane of the river, parts of the solution are considered separately by writing(2.17)

## 3. The evolution in time from rest

A similarity solution in the neighbourhood of (±*a*,0) shows, by enforcing the river and free-surface conditions (2.13) and (2.14), that the normal derivative has a square-root edge singularity at each point. Thus, condition (2.14) is satisfied by setting(3.1)where, with *T*_{m} denoting a Chebyshev polynomial,(3.2)and the unknown functions *α*_{n}(*t*), *β*_{n}(*t*) are to be determined by the later application of (2.13). This strategy is a standard technique for handling mixed boundary-value problems: the square root singularity is taken into account explicitly and orthogonality relations for Chebyshev polynomials are used. Constants *A*_{n}, *B*_{n} associated with the steady state, given by as *t*→∞, are introduced for convenience. Note that *V* is ultimately identified with the upward pressure gradient with , which, in the context of Darcy's Law, represents a flux of fluid into the aquifer. Authors commonly avoid this direction contrast by working with drawdown instead of pressure head.

Define even and odd Fourier transforms by(3.3)

Substitution of (3.3) into (2.16) and (3.1) then shows that, for each *k*(>0),(3.4)(3.5)

The time dependence is conveniently handled by use of the Laplace transform,(3.6)which, when applied to (3.4) and the conditions (2.12), (3.5), yields, after invoking (2.11),(3.7)(3.8)(3.9)Define, for each *k*≥0, the shift(3.10)

The solution pair of (3.7) that satisfies (3.8), (3.9) is given by(3.11)Evidently, are analytic functions of *q* and hence *p* except for simple poles at *p*=0, *q*=*q*_{0} given by(3.12)and given by(3.13)whose locations depend on *k*^{2}. It is convenient to include the pole at *p*=*q*_{0}−*k*^{2} in the sequence by setting .

Inversion of the Laplace transforms is according to the formula(3.14)where Re(*γ*)>0 so that all poles lie to the left of the integration path. The forcing terms in (3.11) yield a sum of readily evaluated residues but the series give an integral and a convolution involving the function *G*(*z*, *t*) whose Laplace transform is given by(3.15)The function *G*(*z*, *t*) solves the problemApplication of the inversion formula to (3.15) yields(3.16)However, the initial condition does not imply that the series vanishes at *t*=0. The impulsive forcing yields *G*(*z*, 0−)=0 butInversion of the Laplace transform pair (3.11) now yields(3.17)where(3.18)Note that *g* is symmetric in *z*_{0}, *z* and such that .

It is apparent that the inversion formula for contains a divergent integral, an expected feature because the solution of (2.16) contains the fundamental solution of the two-dimensional diffusion equation, , and hence a logarithmic singularity for all *t*>0. Divergent integrals can be avoided by subtracting out this singularity or, more simply, by focusing attention on the derivatives of . These, of course, are the functions of physical interest. Substitution of (3.17) into the inversion of (3.3) yields(3.19)The far-field behaviour is determined by noting that (3.12) implies when *k*≪1 while (3.13) ensures that the *m*≥1 contributions are exponentially small in time. An integration by parts facilitates the estimateOnly the *J*_{0} term in the series makes a significant contribution (Gradshteyn & Ryzhik 2000, §§6.693.1 and 6.693.2) and hence(3.20)for *x*≫1 at given *t*. Thus, there is no disturbance at *∞* at any finite *t*, consistent with assumed ∂*h*/∂*y*→0 as |*y*|→∞ in (2.15) and condition (2.6).

The function sets and are determined by setting the expressions (3.19) equal to zero on the interval , according to (2.13). The time-dependent disjoint integral equations thus obtained are converted to sets of time-dependent linear equations by the use of suitable orthogonal functions. The odd and even Chebyshev polynomials yield, respectively,(3.21)In a Darcy's Law interpretation, the fraction of flux drawn from the river is given, according to (3.1) and (3.2), by(3.22)where *G* is given by (3.16).

## 4. Numerical formulation

### (a) Formulation in the Laplace variable

In the analysis above, the Laplace transform is inverted first in order to obtain solutions that are real valued and display the expected features of diffusion. But for computational purposes the Fourier transform is inverted first in order to take advantage of an efficient scheme for inverting Laplace transforms.

When, as in (3.21), the river condition (2.13) is enforced by use of Chebyshev polynomials, the inversion of (3.3) yieldsOn substitution of (3.11) and use of the definition (3.15), the Laplace transforms of these equations give the linear systems(4.1)to be solved for at points on the chosen contour in the complex *p* plane.

The steady-state problem can be extracted by considering the *p*^{−1} term in (4.1):(4.2)

The fraction of flux drawn from the near half of the river is(4.3)according to (3.1). Note that (2.13) implies that the last integrand is independent of *x*. The Laplace transform of *F* may be decomposed into four parts, of which two come from integrating (3.2):(4.4)The contribution *F*_{2} comes from the odd part of the solution and vanishes when integrating over the entire river. The remaining two parts, which do not appear in the steady state, are(4.5)Their expressions in terms of the expansion coefficients , obtained by substitution of (3.11), are(4.6)

### (b) The steady state

With the steady-state problem derived as the small *p* limit of (4.1), the limit coefficients *A*_{n}, *B*_{n} can be found by letting *t*→∞ in (3.21). Thus, on noting thataccording to (3.15), and reintroducing the Chebyshev polynomial representations of the Bessel functions in order to evaluate the *k* integral in closed form (Gradshteyn & Ryzhik 2000, §3.981.8),(4.7)(4.8)These can be shown to be the same as (4.2).

For *a*≪*D*, which is typical in applications, an estimate of each matrix element is readily deduced (Gradshteyn & Ryzhik 2000, §6.574.2) to be given byIn this approximation, the matrices of coefficients in (4.7) and (4.8) have the respective inversesMoreover, if *x*_{0}=*O*(*D*), the left-hand sides are approximated byandwhich both vanish if *m*>0.

The fraction of flux drawn from the river is readily deduced from (3.22) to have large time limit *A*_{0}. In the above small *a*/*D* approximation,(4.9)

(4.10)

(4.11)Of interest is the *O*(*a*/*D*) term in *B*_{0} that exhibits a flux drawn into the far side and out of the nearside of the river. This feature is further highlighted by noting, from (3.19) with the use of (3.13) and (3.18), thatThus, in the even solution, the sinks draw fluid from both the river and the far field while in the odd solution fluid flows from the source to the sink with some passing through the river.

### (c) Numerical approach

The numerical procedure inverts the Laplace transform numerically using Talbot's (1979) algorithm to solve (4.1) for the expansion coefficients at given values of *p*. These complex-valued matrix equations are truncated at *n*=*N* with 0≤*m*, *n*≤*N*. Subsequently, the coefficients are calculated at the same value of *p* by evaluating (4.6), with the sum truncated at *N*−1 rather than *N* because the final odd coefficient in the truncation is anomalous. This can be seen by considering the finite truncation of the odd infinite matrix in §4*b* and noting the extra entry that leads to inappropriately large .

In the steady case, there is no Laplace transform to consider and hence just a matrix problem to solve. The matrix elements in (4.2) are computed using numerical quadrature. After taking §4*b* into account, the integrand decays exponentially and is easily computed by using the SLATEC routine DQAWF with no trigonometric prefactor and a requested absolute accuracy of 2×10^{−13}. The right-hand side can be computed using either (4.2) or (4.7) and (4.8). In the first case, the integrand decays exponentially unless *z*_{0}=0 in which case DQAWF suffices again with the cos *kx*_{0} and sin *kx*_{0} terms taken into account explicitly alongside the more complicated oscillatory behaviour for large *k* due to the *J*_{2m+1} term. In the second case, DQAG was used with the same accuracy. The two methods gave answers that agreed to at least 28 significant figures.

For the unsteady case, the right-hand side is computed in the same fashion but the integral in the matrix elements poses more difficulty because, while the product of Bessel functions is oscillatory to leading order, the following term in its expansion for large values of *k* is not. Hence, integration routines designed for both oscillatory and non-oscillatory integrands struggle. To deal with this, the matrix elements are expressed as(4.12)The number *c* must be taken so that *kc*>2*n* and *kc*>2*m*+1. Of the four Hankel function products, one decays exponentially in the upper half-plane, the second decays exponentially in the lower half-plane, while the last two decay algebraically in either plane. The first two are hence treated by deforming the integration contour to a straight line in the upper and lower half-planes, respectively. A standard integration routine over an infinite interval is now suitable for all these integrals. Additional care would be required if the poles of were crossed when deforming the contour but this does not occur.

## 5. Results

Assigning the completely arbitrary parameter values *a*=0.9, *x*_{0}=1.2, *D*=1.4 and *z*_{0}=0.3, the convergence of the steady solution with the truncation order is examined by reference to the steady-state coefficients {*A*_{n}} and {*B*_{n}} as functions of *N*. Figure 2*a* displays the values of |*A*_{n}| and |*B*_{n}| for *N*=5, 10 and 20. Exponential decay is seen, with machine accuracy coming in around *N*=15, except for the anomalous value of *B*_{N} mentioned above. Figure 2*b* shows how |*A*_{0}| differs, as *N* increases, from its value at *N*=20. Henceforth, *N*=20 for steady calculations and *N*=5 for time-dependent cases, which correspond to accuracies at approximately 10^{−6}. The restriction to this accuracy in the Laplace inverse transforms leads to rapid calculations.

Now use values appropriate to the physical situation. Given the wide variation in aquifer and river characteristics, representative values are chosen to make the non-dimensional values convenient. The aquifer depth is 60 m and the river width is 30 m. The well is 60 m from the riverbank and the pump is at a depth of 20 m. The aquifer properties are *K*_{h}=144 m d^{−1}, *K*_{v}=9 m d^{−1}, *S*_{y}=0.15 and *S*_{s}=0.0001 m^{−1}. The resulting non-dimensional parameters are(5.1)Note that *a*/*D*=0.0625≪1, corresponding to the narrow-river regime.

Figure 3 shows the accuracy of the steady-state *a*/*D*≪1 approximation for the physical values considered, allowing *z*_{0} to span the whole depth of the aquifer. The approximation is barely distinguishable from the exact solution. Inspection of its sign shows that the small net flux is into the aquifer for deep wells, and out of the aquifer for shallow wells. Integrating the steady-state limit of (2.10) shows that the total flux into the aquifer is 1. Of that, part is through the river and part is from infinity. Evidently, the fluxes through the river for these parameter values are at most 2 per cent of the flux from infinity. The magnitude of the flux depends mostly on the distance of the well from the river.

Figure 4 shows the dividing line between the two cases of flux into and out of the aquifer through the bed of the river in (*x*_{0}, *z*_{0}) space, using the small *a*/*D* result that can be obtained from (4.9). In the small *a*/*D* limit, the steady overturning circulation given by *B*_{0} is always positive as can be seen from (4.10). That is to say, there is always flow into the aquifer over the right-hand half of the river. This is understandable, since the antisymmetric part of the flow has a sink to the right of the river, drawing water in, and a source to the left of the river.

The time-dependent fluxes are shown in figure 5 that displays an immediate rapid adjustment to a maximum flux followed by a slow decay to the steady-state value. One finds that, at least in the small *a*/*D* regime, the flux is positive, corresponding to flow into the aquifer, as expected. The overturning circulation can be negative for very small times, but rapidly becomes positive. The adjustment to the steady state is slow, and it is at these large times that the flux can reverse.

In the small *a*/*D* limit, the approach to the steady state takes a power law form, *t*^{−1/2} to be precise. The coefficient of the *t*^{−1/2} decay can be computed by expanding the integrals(5.2)(5.3)where and are the integrals in (4.2). The only terms that contribute have *m*=0 and *n*=0. Hence(5.4)

(5.5)

Expand the linear equations for small *p* and obtain the correction terms satisfying(5.6)Then from the loop version of Watson's lemma, the corrections to the large time limit are . The small *a*/*D* limit can be found in closed form(5.7)(5.8)The other values are negligible. This is not quite all: corrections to the flux terms (4.6) are also required. These do not require any matrix inversion. As it happens *F*_{4} decays faster than *t*^{−1/2}, while(5.9)This term dominates *F*_{1} in the small *a*/*D* limit, and the result for that limit is(5.10)The inverse square root decay obtained is shown in dots in figure 5 and the following figures.

One can compute the time, *t*_{x}, at which the flux changes from being out of the river to being into the river for shallow wells. Combining (4.9) and (5.10) shows that , so this happens at very large times.

The effect of varying the width of the river, *W*, and the horizontal hydraulic conductivity, *K*_{h}, is shown in figure 6. This affects the steady state since the non-dimensional width of the river scales such as , all else being equal. The steady-state flux can be negative or positive, since the scaled *x*_{0} depends on *K*_{h}, i.e. the well can change from being ‘shallow’ to ‘deep’ without changing the dimensional *z*_{0}. The maximum flux increases with the size of the river, and the decay in time is very similar for all curves. The steady-state overturning circulation also increases with *W*. All these curves are in the small *a*/*D* regime, so the approximations derived above can be used.

The effect of varying the position of the well, i.e. *Z* and *x*_{0}, is shown in figure 7. Once again the curves are qualitatively similar. Shallower wells drive more overturning circulation, which is intuitively reasonable.

## 6. Summary

The major novelty of the present approach is the explicit consideration of the vertical structure of the problem. One of the features thus identified is that the flux of fluid out of the aquifer has a non-trivial dependence on the depth of the well and can reverse as the well is sunk further. The scalings employed, being derived from the hydraulic conductivities, the kinematic porosity and the specific storage coefficient, are natural for the problem and capture the unsteady growth of the flux, followed by a slow power-law decay.

Results have been presented in the small *a*/*D* limit, which is relevant to many real situations. In this limit, explicit expressions for the fluxes of fluid through the system and into and out of the river can be obtained for the steady-state limit, as well as the power-law approach to this limit. The formulae (5.10) are explicit closed-form representations of this behaviour. For wide rivers and shallow aquifers, the dependence of each flux on the physical parameters is more complicated.

The Dupuit approximation eliminates the parameter *a*/*D*, which tends to infinity since the vertical structure is eliminated. The limit *K*_{v}→0 is a singular limit in this respect. The results and *B*_{0}∼(*a*/*D*) when the well is far from the river show a tiny extraction from or input to the stream and that the dominant effect is through flow. These results differ from those in the cited literature, because the model is different. The Dupuit approximation corresponds to *D*→0, and heuristically, in this limit, *A*_{0} and *B*_{0} will increase, up to the *O*(1) values obtained previously (Hunt 1999; Butler *et al*. 2001; Bakker & Andersen 2003).

## Acknowledgments

The authors are grateful to Daniel Tartakovsky for introducing them to this problem and to Vitaly Zlotnik for sending them useful parameter values.

## Footnotes

- Received May 16, 2008.
- Accepted August 19, 2008.

- © 2008 The Royal Society