Semi-analytical 3D solution for assessing radial collector well pumping impacts on groundwater–surface water interaction

Large quantities of groundwater and induced surface water may be withdrawn by pumping wells. While vertical wells are most common, radial collector wells are sometimes used both for water withdrawal and for remediation of contaminated groundwater systems (e.g., Bakker et al. 2005; Yeh & Chang 2013). The radial collector well is typically composed of a group of horizontal wells which radially extend from a caisson. Modeling of the interaction between groundwater, surface water and radial collector wells can provide useful insights into the physics of this complex interaction.

Grid-based numerical models have been widely used to investigate transient flow and advective particle movement between regional groundwater, surface water and wells using finite element (Ameli & Abedini 2016) and finite difference methods (Haitjema et al. 2010; Rushton & Brassington 2013). However, flow to a radial collector well is usually a multi-scale problem that is challenging if grid-based numerical methods are used. To ensure a proper representation of each radial arm with an arbitrary orientation and a small diameter, these models require a high grid resolution, which may lead to computational inefficiency. For example, to properly obtain the drawdown-discharge relationship of a single horizontal well in a small homogeneous aquifer with a regular geometry (a domain of 150 × 480 × 24 m), Haitjema et al. (2010) used a MODFLOW model of 1,846,314 cells. The implementation of the boundary condition along the well screen is also challenging in standard numerical models (Bakker et al. 2005). Flow toward radial arms cannot be explicitly incorporated and is approximated using head dependent boundary cells (Patel et al. 2010) or the drain package (Kelson 2012) accompanied by an entry resistance (conductance factor) in MODFLOW. Additionally, to simulate particle transport, use of a non-uniform random walk particle tracking (RWPT) scheme or advection-dispersion-reaction equation within grid-based flow simulators can be computationally expensive, and can be subject to numerical dispersion and artificial oscillations in the vicinity of each radial arm (Zhan & Sun 2007; Starn et al. 2012).

Grid-free analytical methods are also used for the simulation of pumping impacts in geometrically simplified groundwater systems (e.g., Luther & Haitjema 2000; Chang & Yeh 2007). For example, the analytic element method (AEM) is an efficient analytical scheme to emulate the 3D impacts of arbitrarily oriented pumping wells. The basic idea behind AEM is the representation of pumping wells by analytic elements (e.g., point sink, line sink), where each element has an analytic solution which satisfies exactly the groundwater governing equation. A line sink with a variable strength distribution may represent each arbitrarily oriented well screen, which can realistically and efficiently emulate the non-uniform flow behavior in the vicinity of a long well screen (e.g., Luther & Haitjema 1999; Steward & Jin 2003); no horizontal or vertical grid discretization is typically required. AEM-based models therefore have been widely developed for the simulation of 2D and 3D flow toward vertical well (Luther & Haitjema 1999; Bakker 2010), horizontal well (Bakker & Strack 2003; Steward & Jin 2003; Haitjema et al. 2010) and radial collector wells (Luther & Haitjema 2000; Patel et al. 2010) in aquifers with homogeneous properties. Furthermore, analytical models including AEM can provide fast, precise and continuous particle tracking solutions of flow paths and transit time toward arbitrarily oriented pumping wells (Haitjema 1995; Basu et al. 2012; Zhou & Haitjema 2012). However, in analytical flow and particle tracking models, the surface water geometry and groundwater-surface water interaction are typically simplified. Indeed, these models cannot properly emulate irregular geometry and variable material property (e.g., layer stratification) of the aquifer and surface water bodies or complex 3D flow and particle movement among groundwater, radial collector well and surface water.

Grid-free semi-analytical methods, which benefit from the strength of both analytical and numerical schemes, can alternatively be used to address challenging groundwater-surface water interaction problems with or without pumping wells. For example, Ameli & Craig (2014) and Ameli et al. (2013) have relaxed the constraints of traditional Series Solutions analytical method by enhancing this scheme with a simple numerical least-square algorithm. The resulting semi-analytical free boundary groundwater-surface water interaction models have recently been used and tested to simulate 2D and 3D saturated-unsaturated flow in geometrically complex multi-layer unconfined aquifers with various patterns of vertical heterogeneity (Ameli et al. 2016a, 2016b). These approaches have also been extended to explore the hydrological controls on subsurface transit time distribution (TTD) (Ameli et al. 2016a), subsurface transport of sorbing contaminants (Ameli 2017) and weathering evolution in the critical zones (Ameli et al. 2017). In addition to naturally complex aquifer geometry and stratification, and free boundary conditions applied at water table surface, the geometry and properties of surface water bodies (e.g., lake, river and seepage faces) are appropriately handled in these semi-analytical series solution models.

Here we combine the benefits of AEM and the series solution method to develop a general new groundwater-surface water interaction model in the vicinity of radial collector wells. Based on superposition, the free boundary series solution model developed by Ameli & Craig (2014) for 3D groundwater-surface water interaction is augmented with a set of analytic elements (line sinks) used to represent a radial collector well. The coupled Series-AEM model is able to provide a ‘continuous’ map of velocity and dispersion tensors in the entire domain which facilitates an efficient simulation of advective-dispersive transport through implementation of a continuous non-uniform RWPT scheme. Our integrated flow and transport scheme exactly satisfies the governing equations of saturated flow, precisely meets boundary conditions and continuously tracks particles all the way from water table to surface water and to a radial collector well without discretization artifacts.

By applying our new integrated flow and transport model presented here to two hypothetical situations, we demonstrate that it can be a robust tool for simulating complex interaction between radial collector wells, groundwater, the water table and surface water bodies, including meandering rivers and lakes. The model can determine subsurface flow pathline distribution, source zone extent and the percentage of well water captured from different sources. Potential uses of the model include determining well head protection areas and designing contaminant remediation approaches such as river bank filtration (RBF) or pump and treat systems.

The above equation is the general steady-state Series-AEM solution to groundwater-surface water interaction flow in the vicinity of a radial collector well where (Equation A.1, see Appendix A) is the 3D series solution to steady-state groundwater-surface water interaction flow at the m^th layer of an unconfined aquifer and (Equation A.2) is the AEM solution representing the presence of radial collector well.

The computational domain has a length of L_x and L_y in x and y directions, and is subdivided into M layers each with hydraulic conductivity of K_m (Figure 1). Layers, indexed downward from 1 to M, are bounded by the above and below. No-flow boundary condition is applied along the bottom bedrock (at z = 0) and sides of the aquifer (x = 0 & x = L_x & y = 0 & y = L_y). The free boundary water table surface (a priori unknown interface), ⁠, is defined as the surface with zero pressure head calculated using an iterative scheme (Equation A.8). The predefined constant infiltration flux (R) is applied along the at each iteration. The predefined constant head boundary condition equals the water level stage at surface water body (S) is also applied along the predefined location of surface water bodies (e.g., lake, river). Radial horizontal arms of the same length of are located at the elevation of ⁠. The predefined pumping rate of Q is applied at the location of radial collector well using the AEM method. The uniform head of the radial collector well (H_w) is calculated as part of the solution. The mathematical formulation of both series solution and AEM portions of Equation (3) as well as the implementation of boundary conditions through simple numerical least-square algorithm are explained in Appendices A and B. The Appendices are available with the online version of this paper.

Derivatives of the discharge potential (Equation (3)) with respect to x, y and z provide a continuous field of Darcy fluxes in the x ((⁠⁠), y ((⁠⁠) and z ((⁠⁠) directions throughout the entire saturated zone. Continuous fields of mean pore water velocity (⁠⁠, and ⁠) subsequently are obtained by dividing Darcy fluxes by porosity. These velocities can then be used in non-uniform RWPT scheme to track particles and determine their transit times toward surface water bodies and/or radial collector well as explained in Appendix D.

This section describes an example used to demonstrate the numerical efficiency of the Series-AEM approach for simulation of the interactions between groundwater, a surface water body and a radial collector well in a geometrically complex two-layer unconfined aquifer. In this example, the radial collector well is composed of two arms of the same length of and diameter of 0.5 m located at an elevation of ⁠, in close vicinity of a surface water body (Figure 2). A predefined uniform head of is considered along the surface water body. The hydrological and hydrogeological parameters used are: ⁠, ⁠, R (infiltration rate) = 10⁻² m/d and Q (pumping rate) = 10,000 m³/d.

The flow solution and water table location (Figure 2) were obtained using an iterative scheme to locate the water table surface (Equation A.8) and the least-square method to minimize the errors in each iteration (Appendix B). The solutions were obtained using and control points per interface (the top of modeled domain and layer interface), and control points along each radial arm (see Appendix B). A small number of series terms of J = N = 30 (Equation A.1), and line segments of (Equation A.2) along each arm is used. To satisfy the no-flow condition along the four sides of the domain for the AEM portion of the solution an additional 17 image wells (as shown in Figure A.1) are considered. An additional 18 image wells are also situated below the flat bedrock to satisfy no-flow condition along bottom bedrock for the AEM portion of the solution (Equation A.4). A relaxation factor of is also considered in the iteration scheme (Equation A.8).

As stated in Appendix A, the Series-AEM solution developed here satisfies the governing equation exactly. The no-flow condition along the bottom bedrock by the AEM portion of the solution was exactly satisfied using 18 additional image wells. The solution and the constraint on total inflow into the radial collector well (Equation A.6) were met exactly. However, using least squares there are numerical errors in the implementation of boundary and continuity conditions along the top surface, bottom surface (only series solution portion), layer interface and radial arms. The definition of least-square equations used to construct the system of equations, and the equations used to calculate the normalized numerical errors along each interface (i.e. ⁠) are described in Appendix B and Appendix C, respectively. Figure 3(a) shows the largest normalized flux error across the water table surface occurs directly at the intersection of surface water and aquifer with a maximum of 3%. This error can be attributed to Gibbs phenomenon (Gibbs 1899) at the abrupt transition from Neumann to Dirichlet boundary conditions, where the exact solution is discontinuous. The mean absolute normalized flux error along this interface is 0.04%. At the remaining error evaluation points along the top interface which are in direct contact with the surface water body, the maximum normalized head error is on the order of 0.01% (not shown here). The normalized least-square errors across the layer interface and bottom bedrock are explained in the caption of Figure 3. The head uniformity condition along two radial arms is met with a maximum normalized head error of (not shown here).

In the original example, an initial water table elevation equal to H_r (water level stage at surface water body) was used within the iterative scheme; 60 iterations were required for the convergence of pressure head along water table to zero (Appendix A). In a separate analysis, to ensure that the solution is robust, seven different initial water table elevations (from z = 10 to z = 70 m) were used; the water table surface for all seven cases converged to the original surface as shown in Figure 2(b), implying that the least-square scheme calculated the global minimum of the system.

One hundred and twenty uniformly spaced particle-capture points located at six positions along the perimeter of each arm of the well are used to generate pathlines toward the radial collector well using back tracking of each captured particle (only 40 pathlines are shown in Figure 2(a)). These pathlines can also be used to approximate the percentage of well waters originating from the surface water body source ⁠, by assigning a flow quantity to each pathline. This flow is equal to one sixth of the calculated strength of the line sink segment corresponding to the pathline termination location. From the subsurface pathlines in Figure 2(a), is estimated as 45%.

In the second example, we simulate flow, advective-dispersive particle movement and particle transit time between groundwater, a meandering river and a radial collector well in a homogenous unconfined aquifer with an average vertical thickness of 50 m. This example assesses the impacts of pumping rate and radial collector well location on the percentages of well water captured from different sources, source zone extent and the TTD of groundwater discharged into the radial collector well and the river. The land surface topography used for this example (Figure 4) is taken from the Grand River watershed in southern Ontario. The aquifer is assumed to be homogeneous with ⁠, (longitudinal dispersivity of the porous medium; Appendix D), (transverse dispersivity of the porous medium; Appendix D) and ⁠. Along the meandering river, a uniform surface water hydraulic gradient of 0.21 m/km is assumed with a minimum surface water elevation of 50.18 m at downriver (x = 1,200 m) and maximum head of 50.42 m at upriver (x = 0) of the river.

All the solution parameters (e.g., J, N, ⁠) and the layout of image wells are the same as in the first example. The solution for this hypothetical example was obtained successfully, with least-square errors along evaluation surfaces similar to those in the first example.

In Figure 4, colored lines are pathlines toward radial collector wells originating from river and groundwater. The color scales depict the flux distribution (m/d) across the river bed: for a and a′, Q = 5,000 m³/d; for b and b′, Q = 10,000 m³/d; and for c and c′, Q = 15,000 m³/d. Two different configurations, A (x = 400, y = 500) and B (x = 1,030, y = 350) are considered for the caisson of the radial collector well, which consists of two arms (in x+ and y+ directions) with a length of ⁠, diameter of 0.50 m and an elevation of ⁠. The grey lines are contours of topographic surface elevations. The vertical thickness of the domain is on average 50 m. Configuration B is located in the inner section of a bend while configuration A is further away from the bending branches of the river with an arm parallel and the other arm perpendicular to the neighboring straight branch of the river.

Figure 4(a) and 4(a′) show that with the same pumping rate of Q = 5,000 m³/d, the source zone extent of the radial collector wells for the two configurations is significantly different; for configuration B (inside a bend of the river), the source zone extent is highly controlled by the river shape. The depleted flux distribution also demonstrates that, for configuration B, higher fluxes of river water with a maximum of 0.09m/d is depleted from a small part of the river while this maximum value for configuration A (beside a straight portion of the river) is 0.02m/d. As can be expected the percentage of captured well water originated from river (S_r) for position B (70%) is considerably larger than configuration A (11%). For larger pumping rates, a similar behavior between two configurations was observed. For Q = 10,000 m³/d (Figure 4(b) and 4(b′)) the maximum depleted flux (induced by pumping) at the river bed is 0.18m/d and 0.07m/d, with S_r equal to 80% and 39% for pumping at configurations B and A, respectively. For Q = 15,000 m³/d (Figure 4(c) and 4(c′)), the water is depleted at the river bed with a maximum of 0.27m/d and 0.125m/d, and S_r is equal to 80% and 45% for configurations B and A, respectively. Figure 4 also suggests that by increasing the pumping rate, the lateral source zone extent does not change for configuration B while for configuration A this area increases considerably. This can also be seen from the variation of (S_r) by increasing the pumping rate; while S_r for configuration A increases considerably (from 11% to 45%) when pumping rate is tripled (from Q = 5,000 m³/d to Q = 15,000 m³/d), S_r for configuration B increases only from 70% to 80%. As pumping rate increases from Q = 10,000 m³/d to Q = 15,000 m³/d, S_r for configuration B remains constant.

Transit time (or groundwater age) of water particles (the elapsed time that particles spend traveling through subsurface; Appendix D) and transit-time distribution (the probability density function of transit times; Appendix D) of water particles discharged into the river and/or radial collector are also influenced by pumping rate and the location of radial collector well.

Figure 5(a) and 5(a′) show that for position A of the collector well, as pumping rate increases from 5,000 to 15,000 m³/d, the mean groundwater age and the proportion of early and late arrival waters discharged into the river considerably decreases. The Gamma shape parameter of the probability density function of transit times increases from 0.78 to 1.03 which implies a decrease in the age variability of water discharged into the river. For position B, however, as pumping rate increases the Gamma shape parameter only slightly increases from 0.78 to 0.84 with a smaller decrease in the percentages of early and late arrival of the water particles discharged into the river, and a smaller decrease in mean groundwater age of the water particles discharged into the river compared with position A.

Positions A and B of the collector well are shown in Figure 4. To calculate the transit times of the water particles discharged into the river (top panel), a back tracking scheme from 440 uniformly spaced particle release points (located along the river) toward the phreatic surface was used for all cases. To calculate the transit times of the water particles discharged into the radial collector well (b and b′), a back tracking scheme from 120 uniformly spaced particle release points located along the radial collector well screens was used for all cases. In the bottom panel, the TTDs are dimensionless TTD which were normalized with respect to mean groundwater age (⁠⁠), and denote the ensemble of particle transit times from both groundwater and river sources discharged into the radial collector well.

Figure 5(b) and 5(b′) show that for position A of the radial collector well, the percentage of both early and late arrival waters (relative to mean groundwater age) captured by the radial collector well increases smoothly as pumping rate increases from 5,000 m³/d to 15,000 m³/d. The ensemble of mean groundwater age from both groundwater and river sources discharged to the radial collector well decreases from 1,585 to 847 days. For position B, again, the percentages of both early and late arrival waters (relative to mean groundwater age) captured by the radial collector well increases as pumping rates increases; however, here there is a sharper decay in early transit times compared to position A. This can be intimately tied to the small source zone extent of the radial collector well located in position B where a high percentage of water particles in the vicinity of the radial collector well (relatively) rapidly captured by the well. The ensemble mean age of water particles captured by the pumping well is significantly younger for the radial collector well located at position B compared to position A.

Results suggest that the semi-analytical series-AEM model developed here is able to precisely address 3D steady flow and advective-dispersive transport between the radial collector well, regional groundwater and surface water body despite the complicating presence of a free boundary, infiltration and natural geometry and stratification. The groundwater governing equation was satisfied exactly, and boundary conditions were implemented within an acceptable range of error. This precise flow and transport scheme was then used to assess the impact of pumping well rate and placement on: (1) groundwater transit-time distributions discharged into the surface water body and captured by the well; (2) local source zone extent; and (3) the percentage of well waters captured from different sources (groundwater or surface water). Results suggest that the impact of pumping well rate on the aforementioned measures depends mainly on the shape of meandering river in the vicinity of pumping well.

The general flow and transport model developed in this paper has the ability to precisely discern the sources of the collector well waters and the percentage of captured well waters from a river and/or groundwater. In addition, the grid-free semi-analytical approach presented here can be easily used for various layouts, orientations (e.g., vertical), elevations and length of radial arms; thus it is potentially a useful tool for the engineering design of water remediation and filtration, and groundwater protection systems.

Results show that increasing the pumping rate does not lead to the same consequences for radial collector wells located at different positions with respect to the river. When the well is surrounded by a meandering branch of the river (position B, Figure 4), increasing the pumping rate can only slightly increase the percentage sourced from the river and the source zone extent, but significantly decreases the mean transit time of river water captured by the radial collector well (Figure 5(b′)). These findings are important for the design of RBF systems (radial collector wells installed in the vicinity of surface water bodies to withdraw naturally filtered surface water), where the effectiveness of a given radial collector well in a RBF system is assessed by the ability of the well to capture more river water than groundwater but with sufficient transit time so as to filter out undesirable constituents (Moore et al. 2012).

Our results also show that when in the well is surrounded by a meandering branch of the river a high percentage of early arrival (relative to mean age) water particles captured by the well and this behavior is enhanced as pumping rate increases (Figure 5(b′)). This together with significantly shorter mean groundwater ages (for all pumping rates) and smaller source zone extent of the well located at position B suggest that a high percentage of water particles ‘rapidly’ captured from the vicinity of the radial collector well. For example, for the hypothetical well–river interaction example solved here, 95% of water particles are captured in almost two years by the radial collector well at position B when the pumping rate (Q) is 5,000 m³/d but only 3% are captured by the well in position A. This suggests that more stringent protective policies are required in the vicinity of the pumping wells located in the inner section of the bend of meandering rivers (e.g., position B), as a high percentage of water and potential contaminant at the surface (both in surface water body and/or along land surface) is rapidly captured by the well which gives only a short time to detect the potential groundwater concern. This finding can inform the protection practices (e.g., Well Head Protection Area systems) required around pumping wells, which are used to supply municipal drinking water.

A general semi-analytical Series-AEM model for the simulation of 3D flow and advective-dispersive transport between a radial collector well, groundwater and geometrically complex surface water bodies was developed and assessed. For two hypothetical examples studied here, this grid-free model precisely simulated 3D groundwater–surface water interaction in the vicinity of radial collector wells; mass balance was satisfied exactly over the entire domain except along layer interfaces and boundaries where boundary and continuity conditions were met with an acceptable range of error. This precise semi-analytical groundwater–surface water interaction model showed that the impacts of pumping well rate on the TTDs of water particles discharged into a river and an adjacent radial collector well, the local well source zone extent and the percentage of well waters captured from different sources significantly depend on the river shape in the vicinity of the pumping well. The results demonstrate the value of the model in its practical application to problems of water remediation (e.g., RBF system design) and groundwater protection (e.g., Well Head Protection Area determination).

We thank four anonymous reviewers of the Hydrology Research journal for their constructive comments. This work was funded in part from an Early Researcher Award (#ER10-07-108) to the second author from the Ontario Ministry of Research and Innovation.