Abstract

A generalized model for well hydraulics in confined aquifers is presented. The lattice Boltzmann method (LBM)-based altered-velocity model is used to simulate well hydraulics in homogeneous and heterogeneous confined aquifers for different pumping conditions. LBM is applied to simulate well hydraulics in two different configurations of heterogeneous aquifer, where (1) a circular disk of material or (2) an infinite linear strip of material is embedded in a matrix of differing hydraulic properties. The effect of hydrogeologic boundaries on drawdown curve in confined aquifers is simulated using LBM. The LBM is further applied to simulate drawdown during harmonic pumping well test. The LBM data in lattice units are scaled to physical units using non-dimensional discharge rate and diffusion length. The results from LBM were found to be sensitive to the relaxation parameter and the RMS error for drawdown was found to be less than 1% for later time. This work will verify the potential of LBM to simulate well hydraulics in homogeneous and heterogeneous confined aquifers for constant, step, and harmonic pumping.

INTRODUCTION

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.

LATTICE BOLTZMANN METHODS

The LBM solves the dynamic interaction between hypothetical groups of particles represented by a particle distribution function (fj), which defines the occurrence of a group of particles with momentum at any location x and time t. The fj stream on discrete lattices is updated by a collision mechanism. The linearized discrete Boltzmann equation with Bhatnagar–Gross–Krook (BGK) type collision, as shown in Equation (1), is solved on discrete lattices (Qian et al. 1992): 
formula
(1)
ej is the microscopic velocity of particle groups, fjeq is equilibrium distribution function, and τ is a relaxation parameter that indicates the rate at which the system approaches equilibrium through a series of collision and streaming. The macroscopic properties, such as density (ρ) and macroscopic velocity (u), are calculated by taking zero and first order moment of fj, as (Qian et al. 1992) and . The equilibrium distribution function fjeq for the Navier–Stokes equation (NSE) solver is (Qian et al. 1992): 
formula
(2)

The weights wj 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).

LBM-BASED MACROSCOPIC FLOW MODEL

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).

A resistance force R can be related to the pressure drop across a porous medium and it could be a tensor to account for directional-dependent properties. Freed (1998) defined this approach as an averaging of steady state NSE over the local sub-volume in the Stokes regime. A LBM-based porous media model is an extension of the basic LBM that is obtained by altering the local macroscopic velocity during the collision step. Some new notation is introduced to define the change in macroscopic velocity: is the pre-collision velocity, u′ is the post-collision velocity, and is the mean centred velocity. An external force F is (Freed 1998): 
formula
(3)
which is equivalent to the term ρRq, and 
formula
(4)
The resistance field is related to the permeability as 
formula
(5)
Mean centred velocity is the actual macroscopic velocity, and is calculated as (Freed 1998): 
formula
(6)
where 
formula
(7)
and 
formula
(8)
The equations for the porous media model are (Freed 1998): 
formula
(9)
and 
formula
(10)
where u′ and are given by (6) and (7). The resulting macroscopic hydrodynamics will be governed by: 
formula
(11)
and momentum equations at the N-S scale 
formula
(12)

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.

Transient groundwater flow modeling in confined aquifer

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.

Figure 1

Schematic of pumping well and drawdown in a confined aquifer.

Figure 1

Schematic of pumping well and drawdown in a 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 (L3/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.

Model verification against the field data

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.

Table 1

List of parameters in LB units and physical units for Figure 2

Qreal Treal Sreal rreal QLBM TLBM SLBM rLBM 
(m3/d) (m2/d) – (m) – (lu3/ts) (lu2/ts)   (lu) 
2,929.7 1,367.2 2 × 10−4 122 0.33 0.0703 
2,929.7 1,367.2 2 × 10−4 244 0.33 0.0703 
Qreal Treal Sreal rreal QLBM TLBM SLBM rLBM 
(m3/d) (m2/d) – (m) – (lu3/ts) (lu2/ts)   (lu) 
2,929.7 1,367.2 2 × 10−4 122 0.33 0.0703 
2,929.7 1,367.2 2 × 10−4 244 0.33 0.0703 

Aquifer properties and pumping rate are scaled to LB units using non-dimensional numbers.

A negative flux (QLBM) 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 QLBM = 0.0703 lu3/ts is set at the geometric center (96, 96) of the domain and the observation well is rplu from the pumping well at (96 + rp, 96). The transmissivity (T) of the aquifer is 1 lu2/ts and the storage coefficient is equal to 1; hence, the hydraulic diffusivity is 1 lu2/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.

Figure 2

Drawdown curve in observation well, at distance rp from pumping well, from the LB model (solid line) and field experimental data (Lohman 1972) shown as open symbols are compared. The estimated hydraulics parameters T and S used in LBM are given in Lohman (1972).

Figure 2

Drawdown curve in observation well, at distance rp from pumping well, from the LB model (solid line) and field experimental data (Lohman 1972) shown as open symbols are compared. The estimated hydraulics parameters T and S used in LBM are given in Lohman (1972).

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).

Effect of hydrogeologic boundaries

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 (rr) and imaginary (ri) from the point of interest where drawdown is estimated.

Figure 3

Cross section of an idealized semi-infinite confined aquifer bounded by an impermeable barrier or a river boundary on the right (modified from Batu 1998).

Figure 3

Cross section of an idealized semi-infinite confined aquifer bounded by an impermeable barrier or a river boundary on the right (modified from Batu 1998).

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.

Table 2

Parameter values for barrier and river boundaries as shown in Figure 3

  Qreal T1real Sreal rreal QLBM TLBM  rLBM 
(m3/d) (m2/d) – (m) – (lu3/ts) (lu2/ts) SLBM (lu) 
Figure 3(a) 2,929 1,367 2 × 10−4 560 0.051 0.33 28 
Figure 3(b) 2,929 1,367 2 × 10−4 560 0.051 0.33 28 
  Qreal T1real Sreal rreal QLBM TLBM  rLBM 
(m3/d) (m2/d) – (m) – (lu3/ts) (lu2/ts) SLBM (lu) 
Figure 3(a) 2,929 1,367 2 × 10−4 560 0.051 0.33 28 
Figure 3(b) 2,929 1,367 2 × 10−4 560 0.051 0.33 28 

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.

Figure 4

Drawdown in an observation well located between the barrier wall and the pumping well as shown in Figure 3(a).

Figure 4

Drawdown in an observation well located between the barrier wall and the pumping well as shown in Figure 3(a).

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.

Figure 5

Drawdown observed in a confined aquifer which is bounded on the right side by a river boundary (Figure 3(b)). The drawdown in the observation well is compared against the analytical solution given by Hantush (1959).

Figure 5

Drawdown observed in a confined aquifer which is bounded on the right side by a river boundary (Figure 3(b)). The drawdown in the observation well is compared against the analytical solution given by Hantush (1959).

Error analysis

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 (rp) 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 (QLBM) 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.

Figure 6

Error (E) between LB model and Theis solution (analytical solution) for different τ at R = 2.

Figure 6

Error (E) between LB model and Theis solution (analytical solution) for different τ at R = 2.

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 QLBM = 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.

Figure 7

Error (E) for different R at τ = 1.

Figure 7

Error (E) for different R at τ = 1.

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.

Pumping well hydraulics in heterogeneous confined aquifer

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.

Figure 8

Configuration of heterogeneity in non-uniform aquifers. (a) A radial disk, of radii R and hydraulic properties T1 and S1, is embedded in an infinite matrix of differing hydraulic parameters (T2, S2). (b) A linear strip of width d and hydraulic parameter T1, and S1 is surrounded by aquifer materials of hydraulic properties T2, T3, and S2, S3. The dashed line represents the infinite extent of aquifer in that direction.

Figure 8

Configuration of heterogeneity in non-uniform aquifers. (a) A radial disk, of radii R and hydraulic properties T1 and S1, is embedded in an infinite matrix of differing hydraulic parameters (T2, S2). (b) A linear strip of width d and hydraulic parameter T1, and S1 is surrounded by aquifer materials of hydraulic properties T2, T3, and S2, S3. The dashed line represents the infinite extent of aquifer in that direction.

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−4. The ratio between transmissivity (T) of circular patch and surrounding zone is γ (T1 = γT2; R2 = γR1). 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 rp = 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.

Table 3

Parameter values for simulation of pumping well hydraulics in non-uniform confined aquifers as shown in Figure 8(a)

  Qreal T1 (rp)real R1 QLBM (T1)LBM (rp)LBM τ 
(m3/d) (m2/d) (m) – (lu3/ts) (lu2/ts) (lu) – 
Figure 8(a) 200 666.66 120 0.1 0.1 3.33 12 10 
Figure 8(a) 45 1,500 120 0.01 0.1 33.33 12 100 
Figure 8(a) 300 20,000 120 0.001 0.1 333.33 12 1,000 
  Qreal T1 (rp)real R1 QLBM (T1)LBM (rp)LBM τ 
(m3/d) (m2/d) (m) – (lu3/ts) (lu2/ts) (lu) – 
Figure 8(a) 200 666.66 120 0.1 0.1 3.33 12 10 
Figure 8(a) 45 1,500 120 0.01 0.1 33.33 12 100 
Figure 8(a) 300 20,000 120 0.001 0.1 333.33 12 1,000 

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.

Figure 9

Drawdown at 120m from the pumping well in a non-uniform aquifer as shown in Figure 8(a). The simulation parameters are shown in Table 3.

Figure 9

Drawdown at 120m from the pumping well in a non-uniform aquifer as shown in Figure 8(a). The simulation parameters are shown in Table 3.

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 T1 = T2=T3; 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.

Table 4

Parameter values for simulation of pumping well hydraulics in non-uniform confined aquifers as shown in Figure 8(b)

  Qreal T1 (rp)real R1 QLBM (T1)LBM (rp)LBM γ 
(m3/d) (m2/d) (m) – (lu3/ts) (lu2/ts) (lu) – 
Figure 8(b) 1,000 2,000 40 0.166 0.2 10 
  Qreal T1 (rp)real R1 QLBM (T1)LBM (rp)LBM γ 
(m3/d) (m2/d) (m) – (lu3/ts) (lu2/ts) (lu) – 
Figure 8(b) 1,000 2,000 40 0.166 0.2 10 
Figure 10

Drawdown in observation well at rp = 40m from the pumping well in non-uniform aquifer, as shown in Figure 8(b).

Figure 10

Drawdown in observation well at rp = 40m from the pumping well in non-uniform aquifer, as shown in Figure 8(b).

Periodic pumping test

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 m2. 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.

Figure 11

Well configuration for harmonic pumping well test (Becker & Guiltinan 2010).

Figure 11

Well configuration for harmonic pumping well test (Becker & Guiltinan 2010).

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 = Q0 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.

Table 5

Hydraulic parameters and LB parameters to simulate drawdown during harmonic well test as shown in Figure 11

Qreal Treal Sreal rreal R0 Q0LBM T1LBM SLBM rLBM 
(m3/s) (m2/s) – (m) – (lu3/ts) (lu2/ts)   (lu) 
0.0058 0.046 5 × 10−3 0.466 0.1 1.4 
Qreal Treal Sreal rreal R0 Q0LBM T1LBM SLBM rLBM 
(m3/s) (m2/s) – (m) – (lu3/ts) (lu2/ts)   (lu) 
0.0058 0.046 5 × 10−3 0.466 0.1 1.4 

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.

Figure 12

Drawdown from LB model for the harmonic pumping test for the well configuration as shown in Figure 11.

Figure 12

Drawdown from LB model for the harmonic pumping test for the well configuration as shown in Figure 11.

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.

CONCLUSION

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.

ACKNOWLEDGMENTS

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.

REFERENCES

REFERENCES
Anwar
S.
&
Sukop
M. C.
2009
Regional scale transient groundwater flow modeling using Lattice Boltzmann methods
.
Computers and Mathematics with Applications
58
,
1015
1023
.
Batu
V.
1998
Aquifer hydraulics: a comprehensive guide to hydrogeologic data analysis. Wiley-Interscience, New York
.
Becker
M. W.
&
Guiltinan
E.
2010
Cross-hole periodic hydraulic testing of inter-well connectivity
.
Proceedings, Thirty-Fifth Workshop on Geothermal Reservoir Engineering Stanford University
35
,
292
297
.
Bixel
H.
&
Larkin
B.
1963
Effect of linear discontinuities on pressure build-up and drawdown behavior
.
Journal of Petroleum Technology
15
,
885
895
.
Boonstra
J.
&
Boehmer
W.
1986
Analysis of data from aquifer and well tests in intrusive dikes
.
Journal of Hydrology
88
,
301
317
.
Budinski
L.
,
Fabian
J.
&
Stipić
M.
2015
Lattice Boltzmann method for groundwater flow in non-orthogonal structured lattices
.
Computers & Mathematics with Applications
70
,
2601
2615
.
Butler
J. J.
&
Liu
W. Z.
1991
Pumping tests in non-uniform aquifers – the linear strip case
.
Journal of Hydrology
128
,
69
99
.
Carter
R.
1966
Pressure behavior of a limited circular composite reservoir
.
Old SPE Journal
6
,
328
334
.
Choi
E.
,
Cheema
T.
&
Islam
M.
1997
A new dual-porosity/dual-permeability model with non-Darcian flow through fractures
.
Journal of Petroleum Science and Engineering
17
,
331
344
.
Clarke
D.
1987
Microcomputer Programs for Groundwater Studies
.
Elsevier
,
Amsterdam
,
The Netherlands
.
Cooper
H. H.
&
Jacob
C. E.
1946
A Generalized Graphical Method of Evaluating Formation Constants and Summarizing Well-Field History
.
US Department of the Interior, Geological Survey, Water Resources Division, Ground Water Branch
.
Cooper
H. H.
,
Bredehoeft
J. D.
,
Papadopulos
I. S.
&
Bennett
R. R.
1965
The response of well-aquifer systems to seismic waves
.
Journal of Geophysical Research
70
,
3915
3926
.
Demir
M. T.
,
Copty
N. K.
,
Trinchero
P.
&
Sanchez-Vila
X.
2017
Bayesian estimation of the transmissivity spatial structure from pumping test data
.
Advances in Water Resources
104
,
174
182
.
Ewing
R. E.
&
Lin
Y.
2001
A mathematical analysis for numerical well models for non-Darcy flows
.
Applied Numerical Mathematics
39
,
17
30
.
Ewing
R. E.
,
Lazarov
R. D.
,
Lyons
S. L.
,
Papavassiliou
D. V.
,
Pasciak
J.
&
Qin
G.
1999
Numerical well model for non-Darcy flow through isotropic porous media
.
Computational Geosciences
3
,
185
204
.
Ferris
J. G.
,
Knowles
D.
,
Brown
R.
&
Stallman
R. W.
1962
Theory of Aquifer Tests
.
Water Supply Paper 1536-E
.
US Geological Survey
,
Reston, VA
,
USA
.
Fetter
C. W.
1994
Applied Hydrogeology
.
Prentice Hall
,
Upper Saddle River, NJ
,
USA
.
Freed
D. M.
1998
Lattice-Boltzmann method for macroscopic porous media modeling
.
International Journal of Modern Physics C-Physics and Computer
9
,
1491
1504
.
Gao
Y.
&
Sharma
M.
1994
A LGA model for fluid flow in heterogeneous porous media
.
Transport in Porous Media
17
,
1
17
.
Hantush
M. S.
1959
Nonsteady flow to flowing wells in leaky aquifers
.
Journal of Geophysical Research
64
,
1043
1052
.
Hantush
M. S.
1964
Hydraulics of wells
.
Advances in Hydroscience
1
,
281
432
.
Jacob
C.
1947
Drawdown test to determine effective radius of artesian well
.
Transactions of the American Society of Civil Engineers
112
,
1047
1064
.
Kang
Q.
,
Zhang
D.
&
Chen
S.
2002
Unified lattice Boltzmann method for flow in multiscale porous media
.
Physical Review E
66
,
056307
.
Kolditz
O.
2001
Non-linear flow in fractured rock
.
International Journal of Numerical Methods for Heat & Fluid Flow
11
,
547
575
.
Li
W.
,
Englert
A.
,
Cirpka
O. A.
,
Vanderborght
J.
&
Vereecken
H.
2007
Two-dimensional characterization of hydraulic heterogeneity by multiple pumping tests
.
Water Resources Research
43
.
doi:10.1029/2006WR005333
.
Lin
H. T.
,
Tan
Y. C.
,
Chen
C. H.
,
Yu
H. L.
,
Wu
S. C.
&
Ke
K. Y.
2010
Estimation of effective hydrogeological parameters in heterogeneous and anisotropic aquifers
.
Journal of Hydrology
389
,
57
68
.
Lohman
S. W.
1972
Ground-water Hydraulics
.
Professional paper 708
.
US Geological Survey
,
Reston, VA
,
USA
.
Loucks
T.
&
Guerrero
E.
1961
Pressure drop in a composite reservoir
.
Old SPE Journal
1
,
170
176
.
Mathias
S. A.
&
Todman
L. C.
2010
Step-drawdown tests and the Forchheimer equation
.
Water Resources Research
46
,
1
9
.
doi:10.1029/2009WR008635
.
Mathias
S. A.
,
Butler
A. P.
&
Zhan
H.
2008
Approximate solutions for Forchheimer flow to a well
.
Journal of Hydraulic Engineering
134
,
1318
1325
.
Maximov
V. A.
1962
Unsteady state flow of compressible fluid to wells in a heterogeneous medium
.
Journal of Applied Mechanics and Techical Physics
3
,
109
112
.
Odeh
A.
1969
Flow test analysis for a well with radial discontinuity
.
Journal of Petroleum Technology
21
,
207
210
.
Qian
Y.
,
d'Humieres
D.
&
Lallemand
P.
1992
Lattice BGK models for Navier-Stokes equation
.
Europhysics Letters
17
,
479
484
.
Song
I.
&
Renner
J.
2007
Analysis of oscillatory fluid flow through rock samples
.
Geophysical Journal International
170
,
195
204
.
Soueid Ahmed
A.
,
Jardani
A.
,
Revil
A.
&
Dupont
J.
2016
Joint inversion of hydraulic head and self-potential data associated with harmonic pumping tests
.
Water Resources Research
52
,
6769
6791
.
Spaid
M. A. A.
&
Phelan
F. R.
Jr
1997
Lattice Boltzmann methods for modeling microscale flow in fibrous porous media
.
Physics of Fluids
9
,
2468
2474
.
Thiem
G.
1906
Hydrologische Methoden
.
JM Gebhardt's Verlag
,
Leipzig
,
Germany
.
Todd
D. K.
1959
Groundwater Hydrology
.
Chapman & Hall
,
London
,
UK
.
Van Tonder
G.
,
Botha
J.
,
Chiang
W. H.
,
Kunstmann
H.
&
Xu
Y.
2001
Estimation of the sustainable yields of boreholes in fractured rock formations
.
Journal of Hydrology
241
,
70
90
.
Walsh
S. D. C.
,
Burwinkle
H.
&
Saar
M. O.
2009
A new partial-bounceback lattice-Boltzmann method for fluid flow through heterogeneous media
.
Computers & Geosciences
35
,
1186
1193
.
Wen
Z.
,
Zhan
H.
,
Wang
Q.
,
Liang
X.
,
Ma
T.
&
Chen
C.
2017
Well hydraulics in pumping tests with exponentially decayed rates of abstraction in confined aquifers
.
Journal of Hydrology
548
,
40
45
.

Supplementary data