A generalized model for pumping well hydraulics in confined aquifers

Theis (1935) developed an analytical model called the groundwater flow equation to estimate transient change in head due to pumping in a homogeneous, isotropic confined aquifer of infinite aerial extent. The groundwater flow equation is based on conservation of mass and Darcy's law with assumptions including laminar radial flow, fully compressible aquifer, and negligible well diameter. Cooper & Jacob (1946) proposed an approximate solution for the Theis equation valid for later time (with the dimensionless well function argument u < 0.02). The Cooper–Jacob model is most commonly used to estimate the drawdown after pumping in a confined aquifer as a forward solution and is also used as an inverse method for estimation of aquifer hydraulic properties such as storage coefficient, S (volume of water discharged per unit area per unit drawdown after pumping), and transmissivity, T (product of hydraulic conductivity and thickness of aquifer) using the drawdown data. Jacob (1947) described drawdown as a combination of aquifer and well head losses. Aquifer loss is a linear function of discharge rate but well loss is a non-linear function of discharge rate and hence accounts for the effect of unsteady, non-Darcian flow (often in the well itself) on drawdown. The exponent of non-linearity remains a matter of debate and is generally estimated using drawdown data. Hantush (1964) further improved the accuracy of the Theis model by considering the finite non-zero diameters of pumping wells. Cooper et al. (1965) refers to a well–aquifer system as an analog of a spring-mass system, after a series of field tests showed strong oscillations of the water column in pumping wells. Thiem (1906) had proposed a simple solution to estimate the hydraulics parameter, T, using steady-state water level in two observation wells near the pumping well.

Barker & Herbert (1992) divided well losses into several components: partial penetration, non-Darcian flow, loss in the gravel pack, loss through the screen slots, frictional loss in the well casing, and loss in the well screen. Van Tonder et al. (2001) proposed an empirical relationship between drawdown and discharge to fit the drawdown data obtained from a step drawdown test (SDT) in the Karoo formations of South Africa because the Cooper–Jacob model could not fit the observed data. This empirical model has six parameters and its well loss parameter is made a function of time, contrary to the Cooper–Jacob model where the parameters remain constant. Mathias et al. (2008) presented an approximate solution for drawdown by replacing the Darcy equation with the Forchheimer equation. Mathias & Todman (2010) verified this model against four different sets of field data and the model was able to fit the observed field data which were overestimated by the Cooper–Jacob method. The Forchheimer flow to well model is numerically solved using finite difference (Choi et al. 1997; Ewing & Lin 2001), finite elements (Ewing et al. 1999; Ewing & Lin 2001; Kolditz 2001), Boltzmann's transformation (Sen 1988), and Laplace transformation (Mathias et al. 2008). In all of these models, drawdown is estimated as the sum of aquifer loss and well loss, which are further parameterized as a function of well and aquifer properties. Edwards (2010) cited a relationship between fractional error in drawdown estimates and fractional error in parameters. The fractional error in drawdown was found to be most sensitive to transmissivity and flow rate. Hence, there is a need to develop a model with the least set of parameters to minimize error in modeling results.

Several models have been developed recently to estimate the hydraulic parameters of aquifers. Ackerer & Delay (2010) used the adaptive multi-scale triangulation inversion technique on limestone karstic aquifer well data to estimate the heterogeneity of the aquifer. There are two types of grids in this approach: the parameter grid and the computational grid. Parameters are defined as local values on the parameter grid and their values are superimposed onto the computational grid by linear interpolation. The parameter grid is refined (downscaling of parameters) independent of the computational grid to lower the objective function. Lin et al. (2010) proposed an adjoint state method hybrid with the modified Tabu method to simultaneously optimize anisotropic transmissivity and storage coefficient in homogeneous and heterogeneous aquifers. Wen et al. (2017) presented a semi-analytical model to predict drawdown during exponentially decayed rates of abstraction in confined aquifers and the two inflection points in the drawdown-time curves were used to predict the aquifer parameters. Soueid Ahmed et al. (2016) concluded after a series of numerical tests that the harmonic pumping test can capture the full range of heterogeneity in a confined aquifer. Le Coz et al. (2017) applied multiple-point statistics to recover the hydraulic properties of a karst aquifer in Poitiers, France. The nature of drawdown was found to be strongly dependent on karst features. Demir et al. (2017) applied a Bayesian approach to recover hydraulic conductivity based on drawdown data and its time derivative from multiple pumping sites. The global curve fitting method (GCFM) employs drawdown and recovery data from various pumping tests to recover the hydraulic parameters of a confined aquifer (Li & Qian 2013). Li et al. (2014) developed empirical functions to determine the optimum pumping duration for GCFM, which is based on a change in gradient of the drawdown curve and distance between the pumping well and the observation well. The duration of pumping will vary when there is more than one observation well at varied distances from the pumping well. Houben (2015) reviewed various assumptions regarding flow regimes in porous media, and governing equations for different flow regimes (laminar and inertial). Design criteria for pumping wells is discussed in terms of maximum permissible or critical entrance (or approach) velocity at the screen, erosion and transport of particles, critical Reynolds number, and critical radius to avoid non-linear head losses in the gravel pack and the aquifer. Well hydraulics is shown to be dependent on well geometry such as partial penetration, and diameter and length of screen and borehole.

Ginzburg & d'Humières (2007) applied the lattice Boltzmann method (LBM)-based anisotropic-advection-dispersion solver to simulate saturated flows in layered heterogeneous anisotropic aquifers. Anwar & Sukop (2009) solved the groundwater flow equation with LBM by using the analogy between heat equation and groundwater flow equation. The groundwater flow equation was also solved using the LBM-based macroscopic flow model where damping force in LBM is analytically related to the hydraulic parameters of a confined aquifer. Budinski et al. (2015) presented a solution for two-dimensional isotropic groundwater flow equations in phreatic aquifers using LBM on non-orthogonal grids. The proposed solution compared well against analytical and finite difference solutions.

We are presenting a lattice Boltzmann (LB) method-based model to directly simulate aquifer response to pumping based on discharge rate (Q) and hydraulic parameters: transmissivity (T) and storage coefficient (S). The underlying concept is similar to Cooper et al. (1965)'s spring-mass analogy. The drawdown will directly respond to dynamic interaction between pumping rate, aquifer damping (T), and compressibility (S). The drawdown in the LB model does not have any components such as formation loss or well loss. The LB model does not require the solution of parameterized empirical equations using complex mathematical and statistical approaches as adopted by several models cited above. This work is a continuation of Anwar & Sukop (2009), in which two variants of LBM were applied to simulate boundary-value and transient groundwater flow problems. This article presents some novel ideas and results, such as: (1) a scaling approach is proposed and verified to scale constant and variable pumping rate from LB units to physical units which is extremely for simulating real-world scenarios; (2) the LB model is applied to estimate drawdown in homogeneous and heterogeneous aquifers for different boundary conditions and pumping conditions (constant, step, and harmonic/periodic pumping); (3) the LB model is validated against experimental data for harmonic pumping rate, which is an extremely difficult process to simulate even with more complex models such as Darcy–Forchheimer (Song & Renner 2007). LBM does not use a physical system of units (e.g., SI units or FPS); therefore, it is critical to use a non-dimensional number for scaling LB data to physical units. This paper will verify the potential of LBM to simulate well hydraulics in homogeneous and heterogeneous aquifers for constant, step, and harmonic pumping under different hydrogeologic boundary conditions. The motivation for this work is to establish LBM as a generalized method for modeling transient groundwater hydraulics in confined aquifers. This is a 2-D model and well hydraulics are simulated for a horizontal 2-D plane in a confined aquifer. The well is assumed to be fully penetrating in the third dimension which is the saturated thickness of the confined aquifer. The model assumes homogeneous uniform property in the vertical direction; hence, the proposed 2-D model is not suitable for solving well hydraulics in a confined aquifer which is heterogeneous in the vertical plane. This paper is organized as follows: in the section below LBM is introduced, followed by a section discussing scaling approach and simulation results, and finally, conclusions are presented.

The weights w_j for a D2Q9 model are 1/3 for j = 1, 2, 3, 4 (main Cartesian axes) and 1/12 for j = 5, 6, 7, 8 (diagonals). The weights for the rest of the particles (j = 0) are 4/9 (Qian et al. 1992). The LB model achieves a solution for NSE without directly solving NSE; this makes the LB model appealing for solving flow in complex porous media at the pore scale (Anwar et al. 2008).

Flow in porous media changes at multiple scales: pore, laboratory, and field scale. The dominant processes and macroscopic governing equations may also change with the change in scale (Kang et al. 2002). The standard BGK model in LBM (Equation (2)) is an NSE solver and solves the flow field in domains that are explicitly divided into fluid-filled pore or conduit space, and solids, which define the locations of no-slip boundaries. This approach becomes too computationally intensive when an entire aquifer is under study.

A common method employed to solve such problems is the use of macroscopic properties, such as permeability and hydraulic conductivity, which consider the overall effect of the porous medium on the flow. Thus, Darcy's law, which applies at macroscopic scale, does not depend directly on the scale of the problem and a Darcy-scale LBM avoids the scale problem. Solid particles in a porous medium offer resistance to flow; however, well-connected pores facilitate the flow through the porous medium. There are multiple ways to solve the Darcy flux at macroscopic scale using LBM (Gao & Sharma 1994; Spaid & Phelan 1997; Dardis & McCloskey 1998; Freed 1998; Walsh et al. 2009). The most common approach is to introduce an external force by dynamically changing the local velocity during the collision step (Spaid & Phelan 1997).

This type of flexibility in numerical frameworks is applied in modeling complex flow in karst aquifers (Anwar et al. 2008). The density (ρ) given by the LB model is equivalent to head in the confined aquifer. The pumping and observation wells are assumed to be fully penetrating, the flow is radial, and the confining layers are non-leaky.

The LBM-based altered-velocity flow model is applied to simulate drawdown in confined aquifers after the pumping test. The two hydrologic parameters, transmissivity (T) and storage coefficient (S) are assumed to control the nature and magnitude of drawdown during the pumping test. A typical schematic of pumping well and drawdown are shown in Figure 1. Since both the head and the density are scalar quantity, we are using density (ρ) from the LB model as equivalent to head in the confined aquifer.

Diffusion length (ξ), a non-dimensional number, is used to scale time in lattice units (Anwar & Sukop 2009), ⁠; D = T/S is hydraulic diffusivity, t is time (T), and L is a linear scale in the system – e.g., L is the distance between pumping well and the observation well or the grid size. can be used for L when multiple wells are available in the aquifer; however, grid size is a good choice for L when drawdown is observed in the pumping well or during a slug test. Diffusion length (ξ) is commonly used to make a non-dimensional type curve. Dimensionless discharge rate (κ) (Cunningham & Reinhard 2002) is used to scale pumping rate, where Q (L³/T) is pumping rate. The grid size is used for L in this article. The relationship between the hydraulic parameters T, S and the LBM parameters are explained in Anwar & Sukop (2009). Transmissivity (T) is associated with resistance field R, and the storage coefficient S is equal to the relaxation time τ. The LBM-based altered-velocity flow model is applied to simulate drawdown in homogeneous and heterogeneous confined aquifers.

The drawdowns observed in two observation wells at = 122 m and 244 m from a pumping well are simulated using LBM. The aquifer properties, pumping rate and LBM parameters are listed in Table 1. R =0.33 and τ = 1 are selected arbitrarily and grid size 1 lu = 30.5 m.

A negative flux (Q_LBM) is set at the sink node in the LB model to simulate a pumping well and the drawdown curve is observed at some distance from the pumping well. A confined aquifer in a 192 lu × 192 lu numerical domain is simulated, where all four boundaries are set to be periodic to satisfy the assumption of infinite areal extent. A point sink node of strength equal to Q_LBM = 0.0703 lu³/ts is set at the geometric center (96, 96) of the domain and the observation well is r_plu from the pumping well at (96 + r_p, 96). The transmissivity (T) of the aquifer is 1 lu²/ts and the storage coefficient is equal to 1; hence, the hydraulic diffusivity is 1 lu²/ts. Diffusion length (ξ) is used to scale time steps between physical units and lattice units. The comparison between model result and field data (Lohman 1972) is shown in Figure 2.

The LB solution shows a good match against the field data given in Lohman (1972). This demonstrates the ability of the LBM-based altered-velocity flow model to simulate aquifer response to pumping in a homogeneous confined aquifer. The scaling approaches also proved to be successful in scaling LBM data to physical units. Superposition of the Theis solution can be used to estimate drawdown for multiple pumping wells or SDT in a confined aquifer. See Appendix A and Appendix B for verification of the LB model for different pumping conditions (available with the online version of this paper).

An infinite areal extent of confined aquifer is an assumption inherent in the Theis solution, which is not always valid due to the presence of hydrogeologic boundaries in the confined aquifer. A hydrogeologic boundary can be the edge of the aquifer, a source of recharge to the aquifer such as a lake or stream, or an impermeable boundary that does not allow recharge of the aquifer. Boundaries are generally classified as recharge or barrier boundaries (Todd 1959; Ferris et al. 1962; Fetter 1994).

A validation of the LBM-based well hydraulics model against the analytical solution for drawdown in an observation well near the hydrogeologic boundary is presented. Ferris et al. (1962) presented a theoretical foundation of the effect of one straight impermeable boundary on the drawdown. Figure 3(a) shows a cross section of an idealized confined aquifer with an impermeable barrier boundary on the right. There will be no groundwater flow across the barrier boundary, which serves as a hydrogeologic boundary condition. Image well theory is most commonly applied to quantify the effect of hydrogeologic boundaries on the drawdown (s) during pumping. An image well is assumed to be present on the other side of the hydrologic boundary and assumed to be located exactly the same distance from the hydrologic boundary as the real well. The image well acts like a recharge well when a stream is present as a hydrogeologic boundary; or the image well will act like a discharge well when an impermeable barrier is present as a hydrogeologic boundary. The principle of superposition due to real and imaginary wells is applied to estimate drawdown at any point in the confined aquifer. The radial distance (r) in the Theis solution will be the distance of real (r_r) and imaginary (r_i) from the point of interest where drawdown is estimated.

LB simulation is designed to solve the effect of barrier boundary on drawdown in an observation well during pumping from a confined aquifer. An impermeable boundary is set up in the model using no-slip condition instead of using an imaginary pumping well in the simulation. The parameters are given in Table 2 for hydrologic barrier boundary as shown in Figure 3(a). The numerical domain is 192 lu × 192 lu and bounded on the right side, whereas the other three boundaries are periodic. The real pumping well is located 28 lu from the barrier boundary and the observation well is located 9 lu from the barrier boundary. The scaling factor for spatial dimension (L) is 1 lu = 20 m for both boundary value problems.

The drawdown in the observation well is shown in Figure 4. The observation well is located between the real pumping well and the barrier wall. The drawdown from LBM and analytical solution shows a good match as shown in Figure 4.

Another type of hydrogeologic boundary (Figure 3(b)) is where the right side of the confined aquifer is bounded by a river. A constant head boundary condition (Zou & He 1997) is applied in LBM on the right side of the numerical domain to simulate the river as a hydrogeologic boundary. The numerical domain is 192 lu × 192 lu and other model parameters are shown in Table 2. The real pumping well is located 28 lu from the river and the observation well is located 9 lu from the river. Hantush (1959) presented an analytical solution for drawdown in an observation well in a confined aquifer which is bounded by one recharge boundary (river).

The comparison between the analytical solution and the LBM result for drawdown in an observation well is shown in Figure 5. This verifies the ability of LBM to handle boundary value problems in well hydraulics for confined aquifers.

The sensitivity of the model result to LB parameters is presented in this section. There are two hydraulic parameters, T and S, which depend on LB parameters τ and R. The drawdown (s) in observation well at 4 lu (r_p) from the pumping well is simulated using LB model for different τ and compared against the Theis solution. The domain is 240 lu × 240 lu with uniform R = 2 assumed, and the pumping flux (Q_LBM) at the center of the domain is 0.051 lu/ts. The percent error between the Theis and LB solutions as shown in Figure 6 is estimated by ⁠. The percent error (E) is high at early time steps (early pumping time) and then error decays exponentially to significantly lower values. The scaling factor (ξ and κ) for simulation should be selected so that early time steps in LBM represent a small fraction of the entire pumping duration. When scaling factor is properly selected in the LB model, the non-physical early time behavior will not produce significant error in the solution. It is preferable to keep τ > 0.7 to maintain error lower than 1%. The error is higher at τ closer to 0.5 because that is a stability limit in the LB model, but the error at later time is still less than 1% for τ >0.7.

The LB model sensitivity towards another parameter R is shown in Figure 7. The numerical domain is 240 lu × 240 lu and the pumping well is at the center of the domain. The rate of pumping is Q_LBM = 0.051 lu/ts and the relaxation parameter is τ =1. The LB model result for drawdown (s) at 4 lu from the pumping well is compared against the Theis solution. The percent error (E) between LBM and Theis solution is shown in Figure 7.

As shown in Figure 7, it takes much longer for error to become negligible when R is small (high transmissivity) whereas at R =3, the error decays to less than 0.1% in less than 100 ts. Therefore, a proper scaling factor must be used to ensure that error is quickly reduced to a negligible scale.

The verifications presented above are all based on a simple Theis solution; hence, a more rigorous verification is required to demonstrate the ability of LBM to solve complex well hydraulics problems in non-uniform confined aquifers. The proposed LBM is tested to simulate drawdown in heterogeneous confined aquifers, which cannot be predicted by the Theis solution.

The most commonly applied approach for pumping well hydraulics is based on a groundwater flow equation with assumption of homogeneous, isotropy and infinite areal extent of the aquifer (conventional type-curve method). The capture zone in homogeneous and heterogeneous aquifers for transient and steady-state conditions is remarkably different depending on the degree of heterogeneity in the aquifer. Hydraulic tomography is most commonly used as a field method to resolve spatial variability in hydraulics parameters (T and S) by performing multiple pumping well tests. Inverse modeling using geostatistics is also applied to resolve spatial heterogeneity in the aquifer; however, the inverse modeling is complex and computationally burdensome (Li et al. 2007). Butler (1988) and Butler & Liu (1991) proposed a semi-analytical solution for drawdown in two different configurations of non-uniform aquifers, as shown in Figure 8.

Butler (1988) proposed a semi-analytical solution for drawdown in a non-uniform aquifer configuration where a circular ‘patch’ of material is embedded in a matrix of differing hydraulic properties (Figure 8(a)). The hydraulic properties in the two zones are uniform but different from each other. The pumping well is located at the center of the disk. Butler & Liu (1991) presented another configuration for non-uniform aquifers where an infinite linear strip of material is embedded in a material of differing hydraulic properties (Figure 8(b)). Similar configurations are used by several researchers in groundwater and petroleum engineering (Loucks & Guerrero 1961; Maximov 1962; Bixel & Larkin 1963; Carter 1966; Odeh 1969; Boonstra & Boehmer 1986; Clarke 1987) to understand well hydraulics in heterogeneous aquifers. Butler (1988) showed that it is the contrast between the components which is important, not the property variations within each component. The analytical solution for drawdown is a function of hydraulic properties in each component irrespective of the location of pumping and observation wells. The LB model is applied to simulate head distribution or capture zone in a heterogeneous confined aquifer, as shown in Figure 8(a). The aquifer (240 lu × 240 lu) is split into two zones with different values for T and a fixed S = 2 × 10⁻⁴. The ratio between transmissivity (T) of circular patch and surrounding zone is γ (T₁ = γT₂; R₂ = γR₁). The radius (r) of the central disc is 50 m. The pumping well is set at the center of the disk (Figure 8(a)) and drawdown is measured in an observation well located at r_p = 120 m from the pumping well; hence the pumping and the observation wells are located in different zones of transmissivity. The grid resolution is 1 lu = 10 m and the relaxation parameter (τ) is 1 for all simulations. The parameter values in physical units and LB units are shown in Table 3.

The analytical solution for drawdown at 120 m from the pumping well in a heterogeneous aquifer (as shown in Figure 8(a)) is compared against the LBM solution, as shown in Figure 9. Figure 9 shows drawdown for three different γ. The analytical solution presented in Butler (1988) shows that drawdown at any location is a function of hydraulic properties in each zone; hence, it is not important that one of the wells (pumping or observation) must be located in the same or different component of the non-uniform aquifer.

Another LB simulation is set up in a 240 lu × 240 lu domain to simulate drawdown for an aquifer shown in Figure 8(b). This aquifer essentially has three regions, one in the middle and one on each side. The central infinite strip is 100 m wide and the strips on each side of the central strip are semi-infinite. The pumping well is located in the geometric center of this aquifer (middle of central strip) and the observation well is at 40 m from the pumping well. The hydraulic parameter is T₁ = T₂/γ=T₃; hence, the central strip is γ times more (when γ >1) transmissive than the surrounding region.

The drawdown in the observation well obtained from the model is compared against the analytical solution, as shown in Figure 10. The model data are in good agreement with the analytical solution. The pumping rate and hydraulic parameters for the simulation are given in Table 4. LBM is further verified against field data for drawdown after periodic pumping in a confined aquifer.

Guiltinan & Becker (2015) presented results for harmonic (periodic) hydraulic tests in fractured sandstone in New York, USA to investigate inter-well hydraulic connectivity for geothermal reservoir characterization. Harmonic well testing means periodic excitation of a reservoir. The areal extent of the field test was limited to 100 m². The head in the pumping well was varied in sinusoidal fashion. The pumping well is located in the center and four observation wells were located at a distance of 7 m (diagonally) from the pumping well. The schematic of the well configuration is shown in Figure 11.

The head in the central pumping well is varied by raising and lowering the slug in a sinusoidal fashion using a winch driven by a computer controlled stepping motor. Discharge rate could not be measured but the head could be controlled precisely; hence, the hydraulic parameters S and T could not be estimated. The head in the observation wells was monitored using pressure transducers read by a data logger. The LB model is used to simulate drawdown in well 304. The LB parameters and the hydraulics parameters are given in Table 5. A point sink node in the center of the numerical domain is set and the strength of sink varies sinusoidal Q = Q₀ sin (ωt), where ⁠, P is the period of harmonic pumping, and t is time. Guiltinan & Becker (2015) indicated hydraulic diffusivity decaying by power law during harmonic pumping; hence, we assume that transmissivity is reducing during the harmonic pumping following a power law, to make hydraulic diffusivity a power law function.

The number of cycles from Figure 12 was counted to determine the period of harmonic pumping (P) in well 304 and the non-dimensional number diffusion length (ξ) was used to convert the period into LB units.

The drawdown obtained from the LB model for the harmonic pump test is compared against the field data as shown in Figure 12. The LB data show sinusoidal variation in the head in well 304 during the harmonic pumping with the same period as observed by Guiltinan & Becker (2015). The LB model data are in good agreement with the field data for harmonic pumping.

The LBM has been successfully applied to simulate pumping well hydraulics in homogeneous and heterogeneous confined aquifers. The model is relatively simple to use and does not require parameters for aquifer loss and well loss. The model can simulate aquifer response to the pumping well without going through complex mathematical formulations like Fourier, Boltzmann, or Laplace transformation. The non-dimensional numbers are essential to transform LB data into physical units. The ability to simulate harmonic pumping and periodic drawdown indicates that the model may further be applied to simulate more complex early-time unsteady oscillatory drawdown observed in highly permeable aquifers near the pumping well.

SA is grateful to Dr J. J. Butler, Jr for sharing his computer program to generate analytical solutions published in 1988 and 1991. SA is grateful to Dr Becker for sharing the field data for harmonic pumping shown in Figure 12. SA is grateful to Nvidia GPU grant program for providing a GPU card to execute LBM simulations presented in this article.