Determination of Aquitard Storage from Pumping Tests in Leaky Aquifers

Abstract

The management of groundwater resources requires a thorough understanding of groundwater flow and storage. Aquitards, in particular, can store large amounts of groundwater, but are generally often overlooked or ignored in groundwater surveys. In this study, we show how a pumping test performed in a leaky aquifer can be used to estimate the storage coefficient of an overlying aquitard. A new analytical solution is presented for a constant-rate pumping test in a leaky aquifer taking into account the water released from storage in the aquitard. The utility of the method is demonstrated by a practical application using data from a pumping test conducted in a semi-confined sandy aquifer overlain by an aquitard consisting of a mixture of clay and sand. The results obtained show that the solution enables the accurate calculation of the storage properties of the aquitard.

1. Introduction

Effective management of groundwater resources requires a thorough knowledge of groundwater flow and storage [1]. Important hydrogeologic properties are the ability to transmit and store groundwater, which can vary over a wide range, depending on the type and stratigraphy of the soil layers. Aquifers can transmit significant quantities of water but usually have a low storage capacity, while aquicludes can have a large storage capacity but are virtually impervious. Aquitards show intermediate properties, as they are less permeable than aquifers but can often store large amounts of groundwater. In practice, aquitards can thus produce significant amounts of water from storage that are often much larger than from aquifers. However, research on groundwater resources usually focuses only on aquifers, while the importance of aquitards is often overlooked or ignored.

The storage capacity of soil strata is usually expressed by the specific storage, defined as the volume of water in a unit volume of porous medium released from storage resulting from a unit decline in the hydraulic head. The release in storage is due to the expansion of the water and compression of the soil. The total storage capacity of a saturated soil layer is expressed by the storage coefficient, also often referred to as the storativity (similar to transmissivity), defined as the specific storage coefficient multiplied by the thickness of the layer. An overview of the methods used to determine storage properties has been presented by van der Kamp [2] who recommended the use of in situ techniques, in particular pumping tests. Recent studies determine the specific storage of aquitards using field compaction data [3,4], laboratory testing of soil cores [5], order of magnitude estimates based on geological guidelines [6], or pore pressure responses to barometric pressure [7].

A pumping test is a field experiment in which a well is pumped at a constant flow rate and the drop in the groundwater level is measured in one or more surrounding piezometers or observation wells. Theoretical models are used to interpret the observed groundwater drawdown, i.e., decline in groundwater level, taking into account the type of aquifer and the boundary conditions of the groundwater system. The results are mainly used to estimate the hydraulic properties of the aquifer, usually the hydraulic conductivity and the storage coefficient.

In this study, aquitard storage is derived from a pumping test performed in a leaky aquifer, consisting of a semi-confined aquifer overlain by an aquitard and a phreatic aquifer. Analytical solutions for such systems were presented by Hantush and Jacob [8] who ignored aquitard storage, Hantush [9] who presented approximations for short periods of time and for long periods of time, Neuman and Witherspoon [10] who considered a three-layer system consisting of two leaky aquifers separated by an aquitard, and De Smedt [11] who improved the approximations of Hantush [9] and also presented an analytical solution valid for all times. Moench [12] and Feng and Zhan [13] presented semi-analytical solutions in the Laplace domain, considering additional interfering processes such as well bore storage and well skin, partially penetrating well, and aquifer anisotropy. These solutions are useful for testing hypotheses and the sensitivity of side effects, but are generally not suitable for pumping test analysis due to the large number of unknown parameters. Practical aspects of pumping tests in leaky aquifers can be found in pumping test analyses compendia, e.g., [14,15,16,17,18].

The goal of this study is to present a new analytical solution for analysis of a constant-rate pumping test performed in a leaky aquifer with the main objective of determining the aquitard storage coefficient. The usefulness of the method will be demonstrated by a practical application using data from a pumping test conducted in a semi-confined sandy aquifer overlain by an aquitard consisting of a mixture of clay and sand.

2. Materials and Methods

2.1. Theory

One can consider a fully penetrating pumping well in a leaky aquifer as shown in Figure 1. The aquifer is semi-confined because it is underlain by an impervious base and overlain by an aquitard and a phreatic aquifer with a water table that is assumed to remain constant during pumping. Thus, the groundwater flow to the well comes from three sources: release from storage in the aquifer, release from storage in the aquitard, and leakage from the phreatic aquifer through the aquitard.

Assuming horizontal flow in the aquifer and vertical flow in the aquitard, the resulting drawdowns in the aquifer and in the aquitard are given by the following flow equations [8,9]:

$$S\frac{\partial s\left(r,t\right)}{\partial t}=\frac{T}{r}\frac{\partial }{\partial r}\left[r\frac{\partial s\left(r,t\right)}{\partial r}\right]+K{\left.\frac{\partial s^{\prime }\left(r,z,t\right)}{\partial z}\right|}_{z=0}$$

(1)

$$S_s\frac{\partial s^{\prime }\left(r,z,t\right)}{\partial t}=K\frac{{\partial }^2{s}^{\prime }\left(r,z,t\right)}{\partial z^2},$$

(2)

where $s$ and $s’$ [L] are the drawdowns (change in hydraulic head) in the aquifer and in the aquitard, respectively, $S$[-] is the aquifer storage coefficient, $T$ [L² T⁻¹] is the aquifer transmissivity, $S_s$ [L⁻¹] is the specific storage coefficient of the aquitard, $K$ [L T⁻¹] is the hydraulic conductivity of the aquitard, $r$ [L] is the radial distance measured from the center of the well, $z$ [L] is the elevation, and $t$ [T] is the time since pumping started. The initial condition is

$$s\left(r,0\right)=s^{\prime }\left(r,z,0\right)=0.$$

(3)

The boundary condition for the drawdown in the aquifer at the well, assuming an infinitesimal well radius, is

$$\underset{r\to 0}{\lim }\left[r\frac{\partial s\left(r,t\right)}{\partial r}\right]=-\frac{Q}{2\pi T},$$

(4)

where $Q$ [L³ T⁻¹] is the pumping rate, while infinitely far away from the well, the drawdown remains zero.

$$s\left(\infty ,t\right)=0.$$

(5)

The boundary condition for drawdown in the aquitard at the interface with the aquifer, $z=0$, is

$$s’\left(r,0,t\right)=s\left(r,t\right),$$

(6)

while at the top of the aquitard, $z=D$, the drawdown remains zero.

$$s’\left(r,D,t\right)=0,$$

(7)

where $D$ [L] is the thickness of the aquitard.

Hantush [9] derived a solution for the drawdown in the aquifer by means of the Laplace transform, as follows

$$\tilde{s}\left(r,p\right)=\mathcal{L}\left\{s\left(r,t\right);t\to p\right\}=\frac{Q}{2\pi Tp}{K}_0\left[r\sqrt{\frac{S}{T}\left(p+\frac{\sqrt{p{S}_{s}K}}{S}\coth \sqrt{\frac{p{S}_s{D}^2}{K}}\right)}\right],$$

(8)

where $\tilde{s}$ is the Laplace transformed drawdown in the aquifer, $\mathcal{L}$ is the Laplace transform operator, $p$ is the Laplace transform variable, and $K_0$ is the modified Bessel function of the second kind of order zero. This equation can be more easily expressed by introducing the aquitard storage coefficient $S^{\prime }=S_{s}D$ [-] and the aquitard conductance, generally referred to as the leakage coefficient, $C=K/D$ [T⁻¹], as

$$\tilde{s}\left(r,p\right)=\frac{Q}{2\pi Tp}{K}_0\left[r\sqrt{\frac{S}{T}\left(p+\frac{\sqrt{p{S}^{\prime }C}}{S}\cot \mathrm{h}\sqrt{\frac{p{S}^{’}}{C}}\right)}\right],$$

(9)

For the derivation of the analytical solution, we make use of the convolution theorem and the iterated Laplace transform [19] (p. 228), given as

$$\mathcal{L}\left\{{{\displaystyle \int }}_0^{t}f\left(\tau ,t-\tau \right)d\tau ;t\to p\right\}=\mathcal{L}\left\{\mathcal{L}\left\{f\left(t_1,t_2\right);t_1\to p\right\};t_2\to p\right\},$$

(10)

and invert Equation (9) in two steps, writing it as

$$\tilde{s}\left(r,p_1,p_2\right)=\frac{Q}{2\pi T{p}_2}{K}_0\left[r\sqrt{\frac{S}{T}\left(p_1+\frac{C}{S}\sqrt{\frac{p_2{S}^{\prime }}{C}}\cot \mathrm{h}\sqrt{\frac{p_2{S}^{\prime }}{C}}\right)}\right].$$

(11)

First, we invert the equation with respect to $p_1\to t_1$

$${\mathcal{L}}^{–1}\left\{\tilde{s};p_1\to t_1\right\}=\frac{Q}{4\pi T{p}_2{t}_1}\exp \left(\frac{-S{r}^2}{4T{t}_1}\right)\exp \left(\frac{-C{t}_1}{S}\sqrt{\frac{p_2{S}^{\prime }}{C}}\cot \mathrm{h}\sqrt{\frac{p_2{S}^{\prime }}{C}}\right),$$

(12)

and next with respect to $p_2\to t_2$

$$\begin{gathered} {\mathcal{L}}^{–1}\left\{{\mathcal{L}}^{–1}\left\{\tilde{s};p_1\to t_1\right\};p_2\to t_2\right\}= \\ \frac{Q}{4\pi T{t}_1}\exp \left(\frac{-S{r}^2}{4T{t}_1}\right){\mathcal{L}}^{–1}\left\{\frac{1}{p_2}\exp \left(\frac{-C{t}_1}{S}\sqrt{\frac{p_2{S}^{\prime }}{C}}\cot \mathrm{h}\sqrt{\frac{p_2{S}^{\prime }}{C}}\right);p_2\to t_2\right\}, \end{gathered}$$

(13)

which, given Equation (10), leads to the following solution:

$$s\left(r,t\right)=\frac{Q}{4\pi T}{{\displaystyle \int }}_0^t\exp \left(\frac{-S{r}^2}{4T\tau }\right)F\left(\frac{C\tau }{S},\frac{C\left(t-\tau \right)}{S^{\prime }}\right)\frac{\mathrm{d}\tau }{\tau },$$

(14)

where the F-function is given by

$$F\left(x,y\right)={\mathcal{L}}^{–1}\left\{p^{–1}\exp \left(-x\sqrt{p} \cot \mathrm{h}\sqrt{p}\right);p\to y\right\},$$

(15)

De Smedt [11] presented approximations for small and large values of $y$ and an exact analytical solution for all values of $y$ by deriving the inverse Laplace transformation using the Bromwich integral. However, a simpler analytical solution can be obtained using a theorem for the Laplace transform inversion by means of the Fourier cosine transform [20].

First, Equation (15) is written as

$$F\left(x,y\right)={{\displaystyle \int }}_0^y{\mathcal{L}}^{–1}\left\{\exp \left(-x\sqrt{p} \cot \mathrm{h}\sqrt{p}\right);p\to y\right\}\mathrm{d}y.$$

(16)

Next, the Laplace transform inversion theorem using the Fourier cosine transform [20] is applied, resulting in

$$F\left(x,y\right)={{\displaystyle \int }}_0^y\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }\mathcal{R}\left[\exp \left(-x\sqrt{i\omega }\coth \sqrt{i\omega }\right)\right]\cos \left(\omega y\right)\mathrm{d}\omega \mathrm{d}y,$$

(17)

where $\mathcal{R}$[ ] denotes the real part of the expression inside the brackets, i is the imaginary unit, and $\omega$ is a dummy variable. Changing the order of integration yields

$$F\left(x,y\right)=\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }\mathcal{R}\left[\exp \left(-x\sqrt{i\omega }\coth \sqrt{i\omega }\right)\right]\sin \left(\omega y\right){\omega }^{-1}\mathrm{d}\omega ,$$

(18)

and evaluating the real part of the complex exponential function results in

$$F\left(x,y\right)=\frac{2}{\pi }{{\displaystyle \int }}_0^{\infty }\exp \left(-x{f}_1\right)\cos \left(x{f}_2\right)\sin \left(\omega y\right){\omega }^{-1}\mathrm{d}\omega ,$$

(19)

$$f_{1,2}=\sqrt{\frac{\omega }{2}}\left[\frac{\sin \mathrm{h}\left(\sqrt{2\omega }\right) \pm \sin \left(\sqrt{2\omega }\right)}{\cos \mathrm{h}\left(\sqrt{2\omega }\right)-\cos \left(\sqrt{2\omega }\right)}\right],$$

(20)

where the plus on the right corresponds to $f_1$ and the minus to $f_2$.

Thus, a Fourier type integral is obtained, which can easily be evaluated numerically due to the fast decaying exponential term in the integrand. Equations (14), (19) and (20) complete the analytical solution.

Note that when storage in the aquitard is ignored ($S^{\prime }=0)$, the F-function in Equation (14) becomes $F\left(C\tau /S,\infty \right)=\exp \left(-C\tau /S\right)$, and Equation (14) reduces to

$$s\left(r,t\right)=\frac{Q}{4\pi T}{{\displaystyle \int }}_0^t\exp \left(\frac{-S{r}^2}{4T\tau }-\frac{C\tau }{S}\right)\frac{\mathrm{d}\tau }{\tau },$$

(21)

which is the classical solution of Hantush and Jacob ignoring storage in the aquitard [8].

2.2. Test Case

To demonstrate the usefulness of the methodology, it is applied to a practical pumping test case, conducted at an industrial site near Oudenaarde, Belgium (Figure 2). The fully penetrating well is situated in a leaky aquifer composed of Pleistocene sand, located at a depth of 16 m to 27 m below ground level. The aquifer is underlain by a very thick clay substrate, which can be considered impervious, and overlain by a mixture of Pleistocene and Holocene sand and clay, which can be considered semi-pervious and therefore acts as an aquitard. The water table is at a depth of about 1.5 m below ground level. A series of ponds and ditches (Figure 2) enable us to stabilize the position of the water table during the pumping test.

The pumping tests lasted for about two and a half days. The pumping rate was monitored by an in-line volumetric flow meter, read at varying time intervals. The pumping rate varied between 360 and 400 m³/d and averaged 380 m³/d. Groundwater drawdowns were observed at six locations, specifically in three fully penetrating observation wells at 14 m, 45 m and 85 m distance from the pumping well and in three piezometers with screens of 2 m length located in the middle of the aquifer at 13 m, 86 m and 262 m distance from the pumping well. Water levels in the observation wells and piezometers were measured manually with portable water level sensors, initially at one-minute intervals for observations near the pumping well and later gradually expanded to larger intervals.

The observed drawdowns versus distance and time (Figure 3 and Figure 4) consist of 175 data points. The observed drawdowns range from 1 mm to 2.26 m and exhibit the typical behavior of a pumping test; they decrease with the radial distance from the pumping well and gradually increase with the duration of the pumping test. A gap in the data at about 0.09 days is due to a small drop in the pumping rate, while later larger gaps in the data correspond to overnight breaks.

3. Results

The analytical solution given by Equations (14), (19) and (20) is fitted to the drawdown observations by minimizing the squared differences between predicted and observed drawdowns using a non-linear least squares optimization algorithm, which enables the estimation of four parameters: transmissivity T and storage coefficient S of the aquifer, and leakage coefficient C and storage coefficient S’ of the aquitard. The convolution integral in Equation (1) is computed numerically using the QDAGS routine of the IMSL library [22] and the Fourier integral of Equation (6) with the QDAWF routine of the same library; both routines allow us to obtain results with absolute and relative errors of less than 10⁻⁸ in double precision.

Suitable starting values for the model parameters are obtained as follows. In the early stage of the pumping test, the aquifer behaves as if it is confined and the drawdown can be approximated by the Theis well flow equation [21]. When the observed data s versus t/r² are plotted on log-log paper, a Theis type-curve can be fitted as an upper bound as shown in Figure 3. The model parameters are obtained as explained in the literature [14,15], resulting in T ≈ 75 m²/d and S ≈ 3 × 10⁻⁴. In the late phase of the pumping test, the aquifer behaves as if it is semi-confined and the flow becomes stationary. The drawdown can thus be described by the steady-state Hantush well flow equation [8], which depends on a leakage factor $B=\sqrt{T/\mathrm{C}}$ related to the radius of influence of the pumping well. Considering the local scale of the pumping test, we can assume that B = 300 m and the leakage coefficient can be estimated as C = T/B² ≈ 8. × 10⁻⁴ d⁻¹. Lastly, the starting value for the aquitard storage coefficient is assumed to be the same as for the aquifer, so S’ ≈ 3 × 10⁻⁴.

With these initial estimates, the optimization algorithm converges with a residual standard error of 0.034 m. The optimized model parameters and fitting criteria are listed in

Table 1. Calibrated aquifer and aquitard parameters (mean estimates and 95% confidence intervals) with storage in the aquitard, Equation (14), and without storage in the aquitard, Equation (21); also given are model fitting criteria: degree of freedom (DF), residual sum of squares (RSS), residual standard error (RSE), Akaike information criterion (AIC), and Bayesian information criterion (BIC).

Parameter	Units	Equation (14)	Equation (21)
T	m2/d	71.6 ± 0.9	75.1 ± 1.8
S	10−4	2.73 ± 0.16	4.75 ± 0.29
C	10−3 d−1	1.96 ± 0.11	2.07 ± 0.22
S’	10−4	15.4 ± 1.6	0
DF	-	171	172
RSS	m2	0.203	0.862
RSE	m	0.034	0.071
AIC	-	−676	−425
BIC	-	−660	−412

. Figure 4 shows a comparison between the observed and fitted drawdown versus time. The agreement between the observed drawdown and model predictions appears to be very close for all observation wells. The differences are between −8 cm and 9 cm. The maximum deviation of 9 cm occurs in the observation well at a distance of 14 m from the pumping well early after the start of the pumping test.

All model parameters are estimated with high confidence, as evidenced by the small 95% confidence intervals. The estimates are also acceptable from a physical point of view. Taking into account a thickness of 11 m for the aquifer and 14.5 m for the aquitard, the hydraulic conductivity becomes 6.51 m/d for the aquifer and 0.028 m/d for the aquitard, which is within expectations for these types of media. Note that the contrast between the aquifer and aquitard hydraulic conductivity is more than two orders of magnitude, justifying the assumption of horizontal flow in the aquifer and vertical flow in the aquitard [10]. For the specific storage, we find 2.5 × 10⁻⁵ m⁻¹ for the aquifer and 1.1 × 10⁻⁴ m⁻¹ for the aquitard, which is also within expectations for these type of porous materials, e.g., [15]. The specific storage of the aquitard is about four times larger than that of the aquifer. Yeh and Huang [23] analyzed different data sets from pumping tests in leaky aquifers, and suggested that aquitard storage can be ignored only if the ratio of the aquitard storage to the aquifer storage is less than 10⁻³.

To emphasize the importance of the aquitard storage, the Hantush and Jacob well flow equation [8] that ignores aquitard storage, given by Equation (21), is also fitted to the observational data, by estimating only three parameters: transmissivity T and storage coefficient S of the aquifer and leakage coefficient C of the aquitard. Starting with the same initial values as before, the optimization algorithm now converges with a residual standard error of 0.071 m. The results are listed in the last column of Table 1 and the predicted drawdown values are compared to the observations as shown in Figure 4. The agreement between the observed drawdown and model predictions is considerably lower than that obtained with the four parameter model, as the differences now range from −14 cm to 20 cm. The maximum deviation of 20 cm occurs again early after the start of the pumping test in the observation well 14 m from the pumping well. The estimated parameter values are acceptable from a statistical and physical point of view, but comparison of the resulting residual standard error and comparison of the observed and predicted drawdowns (Figure 4) show that the model with aquitard storage is physically more realistic and fits the data better.

4. Discussion

The results obtained for the practical test case, given in Table 1 and Figure 4, show that the fit of the three-parameter model ignoring storage in the aquitard is less good than for the four-parameter model that includes aquitard storage, especially for short periods of time and large distances. Aquitard storage therefore has a profound effect on the drawdown in the leaky aquifer. This implies that estimates of the storage coefficient of an aquitard can be obtained from observations of a pumping test performed in the underlying aquifer, a technique that has been largely ignored in practice, probably due to the lack of useful models to analyze such data.

Although the estimates of the aquifer transmissivity and the aquitard leakage coefficient are almost the same for both models, the estimates for the aquifer storage coefficient are very different and the four-parameter model predicts a storage coefficient of the aquitard approximately twice that of the aquifer. It is clear that in the three-parameter model, which ignores aquitard storage, all storage effects are attributed to the aquifer. Hantush [9] showed that for large distances, the storage coefficient of the aquifer apparently increases by one third of the storage coefficient of the aquitard. In the present case, this would mean S = (2.73 + 15.4/3) × 10⁻⁴ = 7.84 × 10⁻⁴, but this is larger than the value obtained by fitting the three-parameter model. The reason may be due to the fact that all observations are used in the fitting procedure, rather than just the observations at large distances. Nevertheless, these results indicate that estimates of aquifer storage can be inaccurate when storage in the aquitard is ignored. Rushton [24] studied the effect of aquitard storage using a numerical model and also concluded that unreliable parameter values may be deduced in a pumping test analysis if aquitard storage is ignored.

The model fitting criteria also clearly indicate that the four-parameter model including aquitard storage is superior. The residual sum of squares (RSS) obtained with the four-parameter model is considerably lower, while the residual standard error (RSE), i.e., the square root of the residual sum of squares divided by the residual degrees of freedom, is also significantly lower than for the three-parameter model. It is therefore justified to consider aquitard storage as an additional model parameter because this allows a much better fit of the observations.

The Akaike Information Criterion (AIC) is a statistical method for evaluating how well a model fits the data, and is used to compare different models and determine which is the best [25,26]. AIC is calculated based on the number of independent variables used to build the model and the maximum likelihood obtained with the following model [26]:

$$\mathrm{AIC}=2k-n \ln \left(\widehat{L}\right),$$

(22)

where $k$ is the number of model parameters, including the standard deviation of the residuals (differences between observed and predicted values), $n$ is the number of fitted data, and $\widehat{L}$ is the maximized value of the likelihood function for the model. The best model is the model that has the lowest AIC score. Similarly, the Bayesian Information Criterion (BIC) is a statistical criterion for model selection closely related to AIC, also based on the goodness of fit and model parameters, but using a different penalty for the number of parameters [26].

$$\mathrm{B}\mathrm{IC}=k\ln \left(n\right)-n \ln \left(\widehat{L}\right).$$

(23)

The model with the lowest BIC is preferred. An easy way to interpret models and allow a quick comparison and ranking is to look at the differences in AIC or BIC scores. A model that scores more than 10 points higher than another model can be discarded [26]. In the present case, the differences between both models for the AIC values and BIC values are more than two hundred, indicating an overwhelming strength-of-evidence against ignoring aquitard storage.

5. Conclusions

The new analytical solution for a constant-rate pumping test in a leaky aquifer taking into account storage in the aquitard presented in this study extends the work of Hantush [9] and De Smedt [11] because it is simpler from a mathematical point of view and valid for all time periods. It is also an alternative for the Neuman and Witherspoon [10] solution because the current solution is less complex and easier to compute and to apply in practice.

A computational procedure for the analysis of a pumping test conducted in a leaky aquifer was presented that enabled us to estimate the aquifer and aquitard properties by minimizing the sum of the squared differences between observed and predicted drawdowns. The practical application shows how the hydraulic properties of the system can be determined and accurate estimates can be made of the transmissivity and storage coefficient of the aquifer and of the leakage coefficient and storage coefficient of the aquitard. The fit between the model results and the observed drawdowns clearly demonstrate that the solution including storage in the aquitard is superior to the approach ignoring aquitard storage.

Hence, the solution presented in this paper appears to be very useful and accurate in practice to determine aquifer storage properties from the analysis of constant-rate pumping tests in leaky aquifers.