A general analytical model for head response to oscillatory pumping in unconfined aquifers : effects of delayed gravity drainage and initial condition

Oscillatory pumping tests (OPTs) provide an alternative to constant-head and constant-rate pumping tests for determining aquifer hydraulic parameters when OPT data are analyzed based on an associated analytical model coupled with an optimization approach. There are a large number of analytical models presented for the analysis of the OPT. The combined effects of delayed gravity drainage (DGD) and the initial condition regarding the hydraulic head are commonly neglected in the existing models. This study aims to develop a new model for describing the hydraulic head fluctuation induced by the OPT in an unconfined aquifer. The model contains a groundwater flow equation with the initial condition of a static water table, Neumann boundary condition specified at the rim of a partially screened well, and a free surface equation describing water table motion with the DGD effect. The solution is derived using the Laplace, finite-integral, and Weber transforms. Sensitivity analysis is carried out for exploring head response to the change in each hydraulic parameter. Results suggest that the DGD reduces to instantaneous gravity drainage in predicting transient head fluctuation when the dimensionless parameter a1 = Syb/Kz exceeds 500 with empirical constant , specific yield Sy, aquifer thickness b, and vertical hydraulic conductivity Kz. The water table can be regarded as a no-flow boundary when a1 < 10−2 and P < 104 s, with P being the period of the oscillatory pumping rate. A pseudo-steady-state model without the initial condition causes a time-shift from the actual transient model in predicting simple harmonic motion of head fluctuation during a late pumping period. In addition, the present solution agrees well with head fluctuation data observed at the Savannah River site. Highlights. An analytical model of the hydraulic head due to oscillatory pumping in unconfined aquifers is presented. Head fluctuations affected by instantaneous and delayed gravity drainages are discussed. The effect of the initial condition on the phase of head fluctuation is analyzed. The present solution agrees well with head fluctuation data taken from field oscillatory pumping.


Introduction
Numerous attempts have been made by researchers to the study of the oscillatory pumping test (OPT) that is an alternative to constant-rate and constant-head pumping tests for determining aquifer hydraulic parameters (e.g., Le Vine et al., 2016;Christensen et al., 2017;Watlet et al., 2018).The concept of the OPT was first proposed by Kuo (1972) in petroleum literature.The process of the OPT contains extraction stages and injection stages.The pumping rate, in other words, varies periodically as a sinusoidal function of time.
Compared with traditional constant-rate pumping, the OPT in contaminated aquifers has the following advantages: (1) it is low cost because it does not dispose contaminated water from the well, (2) it has a reduced risk of treating contaminated fluid, (3) it has smaller contaminant movement, and (4) it has a stable signal that is easily distinguished from background disturbance such as the tide effect and varying river stage (e.g., Spane and Mackley, 2011).However, the disadvantages of the OPT include the need of an advanced apparatus producing periodic rate.Oscillatory hydraulic tomography adopts several oscillatory pumping wells with different frequencies (e.g., Yeh and Liu, 2000;Cardiff et al., 2013;Zhou et al., 2016;Muthuwatta et al., 2017).Aquifer heterogeneity can be mapped by analyzing multiple data collected from observation wells.Cardiff and Barrash (2011) reviewed articles associated with hydraulic tomography and classified them according to nine categories in a table.
Various groups of researchers have worked with analytical and numerical models for the OPT; each group has its own model and investigation.For example, Black and Kipp (1981) assumed the response of confined flow to the OPT to be simple harmonic motion (SHM) in the absence of the initial condition.Cardiff and Barrash (2014) built an optimization formulation strategy using the Black and Kipp analytical solution.Dagan and Rabinovich (2014) also assumed hydraulic head fluctuation to be SHM for the OPT at a partially screened well in unconfined aquifers.Cardiff et al. (2013) characterized aquifer heterogeneity using the finite element-based COMSOL software that adopts SHM hydraulic head variation for the OPT.In contrast, Rasmussen et al. (2003) found that hydraulic head response tends toward SHM at a late period of pumping time when considering the initial condition prior to the OPT.Bakhos et al. (2014) used the Rasmussen et al. (2003) analytical solution to quantify the time after which hydraulic head fluctuation can be regarded as SHM since the OPT began.As mentioned above, most of the models for the OPT assume hydraulic head fluctuation to be SHM without the initial condition, and all of them treat the pumping well as a line source with infinitesimal radius.
Field applications of the OPT for determining aquifer parameters have been conducted in recent years.Rasmussen et al. (2003) estimated aquifer hydraulic parameters based on 1 or 2 h period of the OPT at the Savannah River site.Maineult et al. (2008) observed spontaneous potential temporal variation in aquifer diffusivity at a study site in Bochum, Germany.Fokker et al. (2012Fokker et al. ( , 2013) ) presented spatial distributions of aquifer transmission and the storage coefficient derived from curve fitting based on a numerical model and field data from experiments at the southern city limits of Bochum, Germany. Rabinovich et al. (2015) estimated aquifer parameters of equivalent hydraulic conductivity, specific storage, and specific yield at the Boise Hydrogeophysical Research Site by curve fitting based on observation data and the Dagan and Rabinovich (2014) analytical solution.They conclude that the equivalent hydraulic parameters can represent the actual aquifer heterogeneity of the study site.
Although a large number of studies have been made in developing analytical models for the OPT, little is known about the combined effects of delayed gravity drainage (DGD), a finite-radius pumping well, and the initial condition prior to the OPT.An analytical solution to such a question will not only have important physical implications but will also shed light on OPT model development.This study builds an improved model describing hydraulic head fluctuation induced by the OPT in an unconfined aquifer.The model is composed of a typical flow equation with the initial condition of a static water table, an inner boundary condition specified at the rim of the partially screened well for incorporating the finite-radius effect, and a free surface equation describing the motion of the water table with the DGD effect.The analytical solution of the model is derived by the methods of the Laplace transform, finite-integral transform, and Weber transform.Based on the present solution, sensitivity analysis is performed to explore the hydraulic head in response to the change in each hydraulic parameter.The effects of DGD and instantaneous gravity drainage (IGD) on the head fluctuations are compared.The quantitative criterion for treating the well radius as infinitesimal is discussed.The effect of the initial condition on the phase of head fluctuation is investigated.In addition, curve fitting of the present solution to head fluctuation data recorded at the Savannah River site is presented.

Mathematical model
Consider an OPT in an unconfined aquifer, illustrated in Fig. 1.The aquifer is of unbound lateral extent with a finite thickness b.The radial distance from the centerline of the well is r; an elevation from the impermeable bottom of the aquifer is z.The well with an outer radius r w is screened from elevation z u to z l .The flow equation describing spatiotemporal head distribution in aquifers can be written as where D r = K r /S s and α = K z /K r , h(r, z, t) is the hydraulic head at location (r, z) and time t, K r , and K z are respectively the radial and vertical hydraulic conductivities, and S s is the specific storage.Considering that the water table is a reference datum where the elevation head is set to zero, the initial condition is expressed as The rim of the well bore is regarded as an inner boundary under the Neumann condition, expressed as where l = z u − z l is screen length, and Q and ω = 2π/P are respectively the amplitude and frequency of the oscillatory pumping rate (i.e., Q sin(ωt)) with a period P .Water table motion can be defined by Eq. (4a) for IGD (Neuman, 1972) and Eq.(4b) for DGD (Moench, 1995): where C y = K z /S y and κ = 1/ , with being an empirical constant and S y being the specific yield.Note that Eq. (4b) reduces to Eq. (4a) when κ → ∞ or = 0.The impervious aquifer bottom is under the no-flow condition: The hydraulic head far away from the pumping well remains constant, written as Define dimensionless variables and parameters as follows: where the over bar stands for a dimensionless symbol.Note that the magnitude of a 1 is related to the DGD effect (Moench, 1995) and γ is a dimensionless frequency parameter.With Eq. ( 7), the dimensionless forms of Eqs. ( 1)-( 6) become, respectively, h = 0 at t = 0, (9) Equations ( 8)-( 13) represent the transient DGD model when excluding Eq. ( 11a) and represent the transient IGD model when excluding Eq. (11b).

Transient solution for unconfined aquifer
The Laplace transform and finite-integral transform are applied to solve Eqs. ( 8)-( 13) (Latinopoulos, 1985;Liang et al., 2017Liang et al., , 2018)).The transient solution can then be expressed as with where c 0 = ap 0 for IGD and is a 1 p 0 /(p 0 +a 2 ) for DGD, i is the imaginary unit, Im (-) is the imaginary part of a complex number, K 0 (-) and K 1 (-) are the modified Bessel functions of the second kind of zero order and one, respectively, and β n is the positive roots of the following equation: tan The method for finding the roots of β n is discussed in Sect.2.3.The detailed derivation of Eqs. ( 14a)-( 14k) is presented in the Supplement.The first term on the right-hand side (RHS) of Eq. ( 14a) exhibits exponential decay due to the initial condition since pumping began; the second term defines SHM with amplitude A t (r, z) and phase shift φ t (r, z) at a given point (r, z).The numerical results of the integrals in Eqs. ( 14b), (14e), and (14f) are obtained by the Mathematica NIntegrate function.

Calculation of β n
The eigenvalues β 1 , . . ., β n and the roots of Eq. ( 15) can be determined by applying the Mathematica function FindRoot based on Newton's method with reasonable initial guesses.The roots are located at the intersection of the curves plotted by the RHS and left-hand side (LHS) functions of β n in Eq. ( 15).The roots are very close to the vertical asymptotes of the periodical tangent function tan β n .When c 0 = ap 0 , the initial guess for each β n can be considered to be β 0,n + δ, where β 0,n = (2n − 1)π/2, n ∈(1, 2, . . .∞), and δ is a small positive value set to 10 −10 .When c 0 = a 1 p 0 /(p 0 + a 2 ), the initial guess is set to β 0,n −δ for a 2 −ζ ≤ 0. There is an additional vertical asymptote at The initial guess is therefore set to β 0,n + δ for β 0,n on the LHS of the asymptote and is set to β 0,n − δ for β 0,n on the RHS.

Pseudo-steady-state solution for unconfined aquifer
A pseudo-steady-state (PSS) solution h s accounts for SHM of head fluctuation at a late period of pumping time and satisfies the following form (Dagan and Rabinovich, 2014): where H (r, z) is a space function of r and z.This defines a PSS IGD model as Eqs.( 8)-( 13), excludes Eqs. ( 9) and (11b), and replaces sin(γ t) in Eq. ( 10) by e iγ t .Substituting Eq. ( 17) and ∂h s /∂t = Im(iγ H (r, z)e iγ t ) into the model results in The resultant model is independent of t, indicating that the analytical solution of H (r, z) is tractable.The Weber transform, defined in Eq. (B1) of the supporting material, may be considered to be a Hankel transform with a more general kernel function.It can be applied to diffusion-type problems with a radial-symmetric region from a finite distance to infinity.For groundwater flow problems, it can be used to develop the analytical solution for the flow equation with a Neumann boundary condition specified at the rim of a finite-radius well (e.g., Lin and Yeh, 2017;Povstenko, 2015).Taking the transform and the formula of e iγ t = cos(γ t) + i sin(γ t) to solve Eqs. ( 18)-( 22) yields the solution of h s , expressed as with the Bessel functions of the first kind of zero order, J 0 (-) and one J 1 (-), as well as the second kind of zero order, Y 0 (-) and one Y 1 (-): e (2+z l )λ w − e (2+z u )λ w + e (2+2z l +z u )λ w , (23l) where λ 2 w = (ξ 2 + iγ )/µ, σ = iγ a, H p = 2/(π µξ λ 2 w ) and D = 2(σ cosh λ w + λ w sinh λ w ), and Re (-) is the real part of a complex number.Again, one can refer to the supporting material for the derivation of the solution.Equation (23a) indicates SHM for the response of the hydraulic head at any point to oscillatory pumping.Note that Eq. (23f) reduces to H m ξ dξ for a fully screened well when z l = 0 and z u = 1.

Pseudo-steady-state solution for confined aquifers
Applying the finite Fourier cosine transform to the model, Eqs. ( 18)-( 22) with S y = 0 (i.e., a = 0) for the confined condition, leads to an ordinary differential equation with two boundary conditions.In solving the boundary-value problem, the solution of h s for confined aquifers can be expressed as Eqs.( 23a)-( 23e) with H (r, z), defined as where λ 2 m = γ i + µ(mπ ) 2 .The derivation of Eq. ( 24) is also listed in the supporting material.For a fully screened well (i.e., z u = 1, z l = 0), the first term of the series (i.e., m = 0) remains and the others equal zero because of sin(z u mπ ) − sin(z l mπ ) = 0.The result is independent of dimensionless elevation z, indicating that the confined flow is only horizontal.

Special cases of the present solution
Table 1 classifies the present solution (i.e., Solution 1) and its special cases (i.e., Solutions 2 to 6) according to transient or PSS flow, unconfined or confined aquifer, and IGD or DGD.Each of the solutions (Solutions 1 to 6) reduce to a special case for a fully screened well.Existing analytical solutions can be regarded as special cases of the present solution, as discussed in Sect.3.4 (e.g., Black and Kipp, 1981;Rasmussen et al., 2003;Dagan and Rabinovich, 2014).

Sensitivity analysis
Sensitivity analysis evaluates hydraulic head variation in response to the change in each of K r , K z , S s , S y , ω, and ε.The normalized sensitivity coefficient can be defined as (Liou and Yeh, 1997) 14k) with the roots of Eq. ( 15) and c 0 = a 1 p 0 /(p 0 + a 2 ) for DGD.Solution 2 is the same as Solution 1 but has c 0 = ap 0 for IGD.Solution 3 equals Solution 1 with Eqs.(16a)-(16d) and β n = 0, π, 2π, . . ., nπ.Solution 4 is the component h SHM of Solution 1 for DGD.Solution 5 consists of Eqs. ( 23a)-(23m) for IGD.Solution 6 consists of Eqs. ( 23a)-( 23e), with H (r, z) defined by Eq. ( 24). a z u = 1 and z l = 0 for fully screened well.b The solution is independent of elevation.
where S i is the sensitivity coefficient of ith parameter, P i is the magnitude of the ith input parameter, and X represents the present solution in dimensional form.Equation ( 25) can be approximated as where P i , a small increment, is chosen as 10 −3 P i .

Results and discussion
The following sections demonstrate the response of the hydraulic head to oscillatory pumping using the present so-

Delayed gravity drainage
Previous analytical models for the OPT consider either confined flow (e.g., Rasmussen et al., 2003) or unconfined flow with IGD effect (e.g., Dagan and Rabinovich, 2014).Little attention has been paid to the consideration of the DGD effect.This section addresses the difference among these three models.Figure 2 shows the curve of the dimensionless amplitude A t at (r, z) = (1, 1) of Solution 1 versus the dimensionless parameter a 1 related to the DGD effect.The transient head fluctuations are plotted based on Solution 1, with a 1 = 10 −2 , 1, 10, 500; Solution 2 is for IGD, and Solution 3 is for confined flow.Define the relative error (RE) as where A t is the dimensionless amplitude predicted by Solution 2 for the case of a 1 = 500 or by Solution 3 for the case of a 1 = 10 −2 .The curves of the RE versus the period of the oscillatory pumping rate (i.e., P ) for these two cases are displayed.The range of P ≤ 10 5 s (1.16 days) contains most practical applications of the OPT.When 10 −2 ≤ a 1 ≤ 500, the A t gradually decreases with a 1 to the trough and then increases to the ultimate value of A t = 1.79 × 10 −2 .The DGD, in other words, causes an effect.When a 1 < 10 −2 , Solutions 1 and 3 agree on the predicted heads; the RE is below 1 % for P < 10 4 s (2.78 h), indicating that the unconfined aquifer with the DGD effect behaves like confined aquifer and that the water table can be regarded as a no-flow boundary when a 1 < 10 −2 and P < 10 4 s.When a 1 > 500, the head fluctuations predicted by both Solutions 1 and 2 are identical; the largest RE is about 0.45 %, indicating that the DGD effect is ignorable and that Eq. ( 4b) reduces to Eq. ( 4a) for the IGD condition.This conclusion is applicable for any magnitude of P in spite of P > 10 5 s.

Effect of finite radius of pumping well
Existing analytical models for the OPT mostly treated the pumping well as a line source with an infinitesimal radius (e.g., Rasmussen et al., 2003;Dagan and Rabinovich, 2014).The finite difference scheme for the model also treats the well as a nodal point by neglecting the radius.This will lead to significant error when a well has a radius ranging from 0.5 to 2 m (Yeh and Chang, 2013).This section discusses the relative error in predicted amplitude, defined as where A t and A D&R are the dimensionless amplitudes at r = 1 (i.e., r = r w ) predicted by IGD Solution 2 and the Dagan and Rabinovich (2014) solution, respectively.Note that their solution assumes an infinitesimal radius of a pumping well and has a typo in that the term e −D w +1 − e −D w should read e β(−D w +1) − e −βD w (see their Eq.25). Figure 3 demonstrates the RE for different values of radius r w .The RE increases with r w as expected.For case 1 of r w = 0.1 m, both solutions agree well in the entire domain of 1 ≤ r ≤ ∞, indicating that a pumping well with r w ≤ 0.1 m can be regarded as a line source.For the extreme case 2 of r w = 1 m or case 3 of r w = 2 m, the Dagan and Rabinovich solution underestimates the dimensionless amplitude for 1 ≤ r ≤ 6 and agrees  (c) displays the curves of S i of h exp and A t at water table (i.e., z = 1) versus S y and ε.
to the present solution for r > 6.The REs for these two cases exceed 10 %.The effect of finite radius should therefore be considered in OPT models, especially when observed hydraulic head data are taken close to the well bore of a largediameter well.

Sensitivity analysis
The temporal distributions of the normalized sensitivity coefficient S i , defined as Eq. ( 26), with X = h exp as in Solution 1, are displayed in Fig. 4a for the response of exponential decay to the change in each of the six parameters, K r , K z , S s , S y , ω, and ε.The exponential decay is very sensitive to variation in each of the parameters K r , K z , S s , and ω because of |S i | > 0. Precisely, a positive perturbation in S s produces an increase in the magnitude of h exp , while that in K r or K z causes a decrease.In addition, a positive perturbation in ω yields an increase in h exp before t = 1 s and a decrease after that time.It is worth noting that the S i for S y or ε is very close to zero over the entire period of time, indicating that h exp is insensitive to the change in S y or ε and that the subtle change of gravity drainage has no influence on the exponential decay.In contrast, the spatial distributions of S i associated with the amplitude A t are shown in predicted head below the water table and slight variation at the water table.Solutions 1 and 2 agree with the h predictions because the head at z = 0.5 under the water table is insensitive to the change in S y or ε, as discussed in Sect.3.3.It is worth noting that the solution of Dagan and Rabinovich (2014) for PSS flow has a time-shift from the h SHM of Solution 2. This indicates the phase of their solution (i.e., 1.50 rad -1 rad approximates 57.30 • ) should be replaced by the phase of Solution 2 (i.e., φ t = 1.64 rad) so that their solution fits the h SHM of Solution 2 exactly.Figure 6 displays head fluctuations predicted by transient Solution 3 expressed as h = h exp +h SHM , and PSS Solution 6 as h s = A s cos(γ t − φ s ) for the partially screened pumping well in Fig. 6a and the full screen in Fig. 6b.The Rasmussen et al. (2003) solution for transient flow predicts the same h as Solution 3. The Black and Kipp (1981) results for PPS flow also predict a h SHM close to the prediction of Solution 3. The phase of Solution 6 (i.e., φ s = 1.50 rad for Fig. 6a and  1.33 rad for Fig. 6b) can be replaced by the phase of Solution 3 (i.e., φ t = 1.64 rad for Fig. 6a and 1.81 rad Fig. 6b) so that the h SHM prediction of Solution 3 is identical to the h s prediction of Solution 6.As concluded, excluding the ini-tial condition with Eq. ( 17) for a PSS model leads to a timeshift from the SHM of the head fluctuation predicted by the associated transient model, while the transient and PSS models give the same SHM amplitude.

Application of the present solution to field experiment
Rasmussen et al. (2003) conducted field OPTs in a threelayered aquifer system containing one surficial aquifer, the Barnwell-McBean aquifer in the middle, and the deepest Gordon aquifer at the Savannah River site.Two clay layers dividing these three aquifers may be regarded as impervious strata.For the OPT at the surficial aquifer, the formation has an average thickness of 6.25 m near the test site.
The fully screened pumping well has a 7.6 cm outer radius.
The pumping rate can be approximated as Q sin(ωt), with Q = 4.16×10 −4 m 3 s −1 and ω = 2π h −1 .The distance from the pumping well to the observation well 101D is 6 m and 11.5 m to well 102D.The screen lengths from the aquifer bottom for well 101D and from the water table for well 102D are 3 m.For the OPT at the Barnwell-McBean aquifer, the formation mainly consists of sand and fine-grained material.The pumping well has an outer radius of 7.6 cm and pumping rate of Q sin(ωt), with Q = 1.19 × 10 −3 m 3 s −1 and ω = π h −1 .The observation well 201C is 6 m from the pumping well.The data of time-varying hydraulic heads at the observation wells (i.e., 101D, 102D, and 201C) are plotted in Fig.   (Wolfram, 1991).Note that a robust Gauss-Newton algorithm (Qin et al., 2018a, b) or an automatic optimization method "SCE-UA" (Duan et al., 1992;Wang et al., 2017) provides an alternative to the parameter estimation.Solutions 4 and 5 are used to predict depth-averaged head, expressed as (z u − z l ) −1 e j , where e j is the difference between predicted and observed hydraulic heads and M is the number of observation data (Yeh, 1987).The estimated parameters and associated SEE and ME are displayed in Table 2.The estimates of T , S, and D r given in Rasmussen et al. (2003) are also presented.The result shows that the estimated S y is very small, and the estimated T and S of Solution 3, Solution 6, or the Rasmussen et al. (2003) solution for confined flow are close to those of Solution 4 or 5 for unconfined flow, indicating that the unconfined flow induced by the OPT in the surficial aquifer is negligibly small.Little gravity drainage due to the DGD effect appears with a 1 = 20 for wells 101D and 102D, as discussed in Sect.3.1.Rasmussen et al. (2003) also revealed the confined behavior of the OPT-induced flow in the surficial aquifer.The estimated S y is 1 order of magnitude less than the lower limit of the typical range of 0.01-0.3(Freeze and Cherry, 1979), which accords with the findings of Rasmussen et al. (2003) and Rabinovich et al. (2015).Such a fact might be attributed to the problem of the moisture exchange limited by capillary fringe between the zones below and above the water table .Several laboratory research outcomes have confirmed that an estimate of S y at a short period of the OPT is much smaller than that determined by constant-rate pumping test (e.g., Cartwright et al., 2003Cartwright et al., , 2005)).In addition, the difference in T , S, or D r estimated by Solution 6 and that estimated by the Rasmussen et al. (2003) solution may be attributed to the fact that their solution assumes isotropic hydraulic conductivity (i.e., K r = K z ).In contrast, the transient Solution 3 gives smaller SEEs than the PSS Solution 6 or the Rasmussen et al. (2003) solution for the Barnwell-McBean aquifer and better fits the observed data at the early pumping periods as shown in Fig. 7. From those discussed above, we may conclude that the present solution is applicable to the real-world OPT.

Figure 1 .
Figure 1.Schematic diagram for oscillatory pumping test at a partially screened well of finite radius in an unconfined aquifer.

Figure 2 .
Figure 2. Influence of delayed gravity drainage on the dimensionless amplitude A t and transient head h at r = 1, z = 1, predicted by Solution 1 for different magnitudes of a 1 related to the influence.

Figure 3 .
Figure 3. Relative error (RE) of the dimensionless amplitudes A t at the rim of the pumping well (i.e., r = r w ), predicted by IGD Solution 2 and the Dagan and Rabinovich (2014) solution.The well radius is assumed to be infinitesimal in the Dagan and Rabinovich (2014) solution and finite in our solution.

Figure 5 .
Figure 5. Head fluctuations at r = 6, predicted by (a) DGD Solution 1 and (b) IGD Solution 2. Solutions 1 and 2 are expressed as h = h exp + h SHM for transient flow.IGD Solution 5, expressed as h s = A s cos(γ t − φ s ), accounts for PSS flow.

3. 4
Figure5demonstrates head fluctuations predicted by DGD Solution 1 and IGD Solution 2, expressed as h = h exp +h SHM for transient flow and expressed by the IGD solution as h s = A s cos(γ t − φ s ) for PSS flow.The transient head fluctuation starts from h = 0 at t = 0 and approaches the SHM predicted by h SHM when h exp ∼ = m after t = 0.5P (i.e., 6 × 10 4 ).Solutions 1 and 2 agree with the h predictions because the head at z = 0.5 under the water table is insensitive to the change in S y or ε, as discussed in Sect.3.3.It is worth noting that the solution ofDagan and Rabinovich (2014) for PSS flow has a time-shift from the h SHM of Solution 2. This indicates the phase of their solution (i.e., 1.50 rad -1 rad approximates 57.30 • ) should be replaced by the phase of Solution 2 (i.e., φ t = 1.64 rad) so that their solution fits the h SHM of Solution 2 exactly.Figure6displays head fluctuations predicted by transient Solution 3 expressed as h = h exp +h SHM , and PSS Solution 6 as h s = A s cos(γ t − φ s ) for the partially screened pumping well in Fig.6aand the full screen in Fig.6b.TheRasmussen et al. (2003) solution for transient flow predicts the same h as Solution 3. TheBlack and Kipp (1981) results for PPS flow also predict a h SHM close to the prediction of Solution 3. The phase of Solution 6 (i.e., φ s = 1.50 rad for Fig.6a and 1.33 rad for Fig.6b) can be replaced by the phase of Solution 3 (i.e., φ t = 1.64 rad for Fig.6a and 1.81 rad Fig.6b) so that the h SHM prediction of Solution 3 is identical to the h s prediction of Solution 6.As concluded, excluding the ini- Figure 6.Head fluctuations at r = 6 predicted by Solutions 3 and 6 for (a) partially screened pumping well and (b) fully screened pumping well.Solution 3 is expressed as h = h exp + h SHM for transient flow.Solution 6, expressed as h s = A s cos(γ t − φ s ), accounts for PSS flow.

Figure 7 .
Figure 7.Comparison of field observation data with head fluctuations predicted by the present solution.Solutions 3 and 6 represent transient and PSS confined flows, respectively.PSS Solutions 4 and 5 stand for DGD and IGD conditions, respectively.
, with the upper elevation z u and lower elevation z l of the finite screen of the observation well 101D or 102D at the surficial aquifer.Note that Solutions 3 and 6 are independent of elevation because of the fully screened pumping well.The standard error of estimate (SEE) is defined as SEE = 1 model but replaces Eqs.(11a) by (11b).Substituting Eq. (17) into the result yields a model that depends on t, indicating that the solution h s to the PSS DGD model is not tractable.

Table 1 .
The present solution and its special cases.

Table 2 .
Rasmussen et al. (2003)timated by the present solution and theRasmussen et al. (2003)solution for OPT data from the Savannah River site.