# Analytically Supported Numerical Modeling of Horizontal and Radial Collector Wells

Publication: Journal of Hydrologic Engineering

Volume 26, Issue 12

## Abstract

Riverbank filtration (RBF) is widely used in drinking water production all over the world. The horizontal and radial collector wells are among the most important water production facilities in case of the RBF aquifers. In this study several types of modeling methods are used to analyze the hydraulic conditions of horizontal and radial collector wells. The unsteady and steady-state flow in homogeneous, infinite, semi-infinite, and square-shaped bounded aquifers are investigated. The effect of delayed gravity response is modeled. Three-dimensional analytical and numerical finite-difference (FD) solution techniques are applied during the simulations to provide more reliable calculation results. Comparative analyses are carried out to optimize the wellbore drawdown in FD simulations. This important value is controlled by the local head loss due to the radial flow in vertical plane normal to the axes of screens. The multimesh FD simulation software FLOW is used for infinite and semi-infinite aquifers, while the MODFLOW Revised Multi-Node Well (MNW2) package is applied to the square-shaped aquifers. The late time partitioning (LTP) method is introduced to estimate the late time drawdown in infinite aquifers using semianalytical steady-state drawdown equations considering circular recharge boundaries. The latter value is derived by the semianalytical CW software. The use of results of various mostly analytically based solutions confirmed or improved the accuracy of the approximate numerical solutions. Such analytical support is vital for proper evaluation of the operation of horizontal and radial collector wells.

## Introduction

Today, drinking water supply is a strategic task all over the world. One of the most important drinking water sources in Central Europe, especially in Hungary, is riverbank filtration (RBF), which is based on abstraction wells located close to rivers (Grischek et al. 2012). RBF offers various advantages, such as protection against chemicals and pathogens, buffering capacity during contaminant transport and extreme hydrologic events, and a pretreatment step during water treatment (Sandhu and Grischek 2012). As a result of the bank filtration (BF) process particles, biodegradable organic compounds, pathogens, and nitrate are completely or partly removed from the surface water (Dragon et al. 2018). At the end of the BF process the quality of raw water often meets the drinking water quality standards of different sites (Hiscock and Grischek 2002). Through RBF, many communities along riverbanks pump drinking water from alluvial aquifers using different types of wells (Ray et al. 2002). In addition to the quality issues of water production, optimization of the quantity of water production should be taken into account. Properly selected and designed water production facilities can contribute to achieving higher yield with one construction. Such facilities may include horizontal wells and radial collector (also referred to as Ranney or star) wells with horizontal arms, which can provide high flow rate. RBF-based drinking water production and horizontal collector wells are used by many cities all over the world such as Warsaw, London, Budapest, and Belgrade (Babac and Babac 2009; Hunt 2003). It is important to know the hydraulics of these complex systems in order to optimize safe and sustainable drinking water production. The aim of our work is to better understand the hydraulic conditions of radial collector wells and to calculate the water levels and flow rates using different modeling methods.

One of the main objectives of hydrogeological modeling is to approximate the hydrogeological observations with the help of different modeling methods (Szucs et al. 2006). These methods can help maintain sustainable water extraction (Buday et al. 2015), especially in drinking water supply.

In this study, the steady-state and unsteady flow to wells with horizontal screens is analyzed. Nonleaky (confined or unconfined) homogeneous and isotropic aquifers are considered. Linear systems with no wellbore storage are assumed in comparative analyses. The nonlinear processes in the screen pipes (axial turbulent friction loss, kinetic energy variation) and the turbulent screen orifice loss are briefly overviewed in the next section. The saturated thickness of the aquifer is assumed as a time-invariant parameter, that is, the thickness reduction effect by depletion of the free surface at the top is ignored. The thickness, hydraulic conductivity, specific storage, and specific yield constitute the hydrogeologic parameters. The number, depth, direction, drilling radius, length, and distance of screened intervals characterize the simulated wells. Zero, one, and four fixed-head boundaries are considered in the presented simulation examples.

In the three presented analyses, infinite, semi-infinite, and square-shaped bounded aquifers are considered. Three solutions (Renard and Dupuy 1991; Morozov 2018; Milojevic 1963) are used to verify the CW software (Szekely 2021), which is utilized in all three examples. Zhan and Zlotnik (2002) and Bakker et al. (2005) provide conceptual and quantitative support for comparative analyses for infinite and bounded aquifers. More methods are used in (1) comparative simulations; and (2) the Appendix demonstrating the technique for unsteady-state drawdown called late time partitioning (LTP).

## Overview of Methods Used in Verifications and Comparative Analyses

Szekely (2021) presented an analytically based method for collector wells in homogeneous, vertically anisotropic aquifers and the related CW software developed for steady-state flow conditions. This program is the base tool providing analytical support.

Fig. 1 shows the vertical section of a radial collector well with arms at two different levels. The two-level option is used in the example of the section “Unsteady Flow to the Two-Level Radial Collector Well in an Unconfined Aquifer.” The five-arm example in the section “Steady-State Flow to the Five-Arm Radial Collector Well in a Square-Shaped Aquifer” demonstrates the case where all the arms are drilled at the same level. Finally, the horizontal well of the section “Unsteady Flow to Horizontal Well in an Infinite Unconfined Aquifer” may be approximated as a two-arm radial collector well with no caisson in the center.

For the screened sections of radially located horizontal arms, the following geometric parameters are considered: number, direction, depth $d$, length $l$, drilling radius ${r}_{d}$, and distance of the midpoint of the screen to the center of the caisson $c$. The caisson is considered as part of the aquifer of thickness $b$; therefore, its diameter and impermeable wall have no hydraulic impact on the flow patterns around the well. The CW software is designed and tested for horizontal and radial collector wells under three types of constant-head boundaries generating recharge through (1) the vertical cylindrical surface (circular recharge boundary), (2) the vertical plane (RBF), or (3) the top of the aquifer (well under the river bed). This software is the main tool in accuracy analysis and analytical support of (confirming or improving) the numerical models.

The closed-form solution by Renard and Dupuy [1991, Eq. (5)] combines two analytical methods. It was used to confirm the CW simulated well drawdown for square-shaped or circular discharge area with the effect of vertical anisotropy (Szekely 2020). The comparative analysis showed a difference in the third digit only. As shown in the Appendix, this solution in combination with the LTP method provides a reasonable approximation for the late time unsteady drawdown in infinite aquifers.

The semianalytical solution by Morozov [2018, Figs. 2(a and b)] was used to verify the CW software for both the steady-state well drawdown and the flux distribution along the well screen. Circular discharge areas with four and eight radially positioned horizontal wells were assumed and the ratio of the flux variation to the uniform drawdown was calculated. For each well layout, three options with different aquifer thicknesses are evaluated. According to Szekely (2020), the 114 data by Morozov (personal communication) and the corresponding CW simulated values exhibit a small mean absolute relative deviation of 0.0159 and a high 0.9995 correlation coefficient.

Milojevic [1963, Eq. (9)] presented an empirical formula to define the relationship between the pumping rate and the steady-state drawdown in the arms of a single radial collector well. The well operates in the semi-infinite aquifer with constant-head boundary that simulates the effect of RBF. The referenced RBF solution for an isotropic aquifer is based on laboratory experiments and considers uniform head in the arms. The presented comparative analysis by Szekely (2020) is performed with the following parameters: aquifer thickness $b=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, drilling depth $d=5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, drilling radius ${r}_{d}=0.25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, length of fully screened arms $l=25\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, number of $\mathrm{arms}=4$, distance to the $\mathrm{river}=200\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. The six listed parameters are in the validity range of the empirical equation published by Milojevic (1963, Table 3). Because no caisson is simulated, the distance required by CW between the origin and the middle of the screened section of the lateral is set to the half-screen’s length of $c=l/2=12.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. The ratio of $Q/k=20$ is applied in the analysis where $Q$ and $k$ are the pumping rate (${\mathrm{m}}^{3}/\mathrm{day}$) and the hydraulic conductivity ($\mathrm{m}/\mathrm{day}$), respectively. The two methods show a low relative deviation of 0.01125. The comparative analysis was extended to wells with 8 and 12 arms. The higher but still low relative deviations of 0.01987 (eight arms) and 0.05219 (12 arms) are acceptable when comparing results obtained with experimental and theoretical analyses.

The CW software has also been successfully verified for nonlinear hydraulic processes. Following Cyr (2007), the steady-state flow around and in the horizontal well was analyzed as outlined next. A presumably infinite and homogeneous aquifer beneath a water reservoir is assumed. Rough, turbulent flow through the screen orifice and in the screen pipe is considered. The impact of the dynamic pressure by the momentum effect is also taken into account. The improved modeling option called Modified HGS-2 by Cyr (2007) and the relevant CW simulation show appropriate agreement of both the calculated heads and fluxes (Szekely 2020).

The finite-difference (FD) numerical groundwater flow modeling is used in all examples. FLOW software is used for simulating the unsteady and steady-state flow in large aquifers of infinite or semi-infinite extension. The MODFLOW software involving the Revised Multi-Node Well (MNW2) package (Konikow et al. 2009) is used to simulate the steady-state flow to the five-arm radial collector well in a square-shaped aquifer. FLOW uses the point-centered FD scheme, whereas MODFLOW employs the block-centered FD scheme.

The analytical methods and WT and WF software for vertical wells allow for unsteady drawdown analysis of single wells or well fields, respectively. These methods are applied to the replacing vertical wells and provide useful tools to confirm the results obtained for horizontal wells. The short, screened replacing vertical wells are positioned along the horizontal screens of the simulated collector well. They are located in the centers of segments or horizontal line sinks constituting the horizontal screens of the CW model. The well line operates at a discharge rate of the segments and generates independent unsteady drawdown data for analytical support.

The FW method combines the FLOW and WF concepts and simulates the drawdown in the wells of the replacing well line. The FW procedure uses the FLOW numerical modeling environment by calling the routine for well drawdown in the nodes.

The MODFLOW MNW2 and FW software use the option of calculating separate drawdown in the wellbore and for the associated nodes of the FD simulation mesh incorporating the relevant screen section. They allow for taking into account the parameters of screen loss or skin effect caused by change of hydraulic conductivity outside the screen.

## Comparative Analysis of Flow to Horizontal and Vertical Wells

In the examples presented in the following subsections, various layouts of the numerical models and pumping wells are applied. Table 1 summarizes the main features of the simulation models.

Example | Model area (${\mathrm{m}}^{2}$) | Number of meshes | Model layers | Laterals of wells | Flow conditions |
---|---|---|---|---|---|

Section “Horizontal Well” | ${10}^{4}\times {10}^{4}$ | 3 | 2–6–12 | 2 | Unsteady |

Section “4-Arms Collector Well” | ${10}^{4}\times {10}^{4}$ | 3 | 2–6–12 | 4 | Unsteady to steady state |

Section “5-Arms Collector Well” | $602\times 602$ | 1 | 12 | 5 | Steady state |

### Unsteady Flow to Horizontal Well in an Infinite Unconfined Aquifer

Zhan and Zlotnik (2002, Section 5.2) analyzed the drawdown caused by horizontal, slanted, and vertical wells in a $b=30\text{-}\mathrm{m}$-thick infinite isotropic water table aquifer exhibiting hydraulic conductivity $k=0.0001\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{s}$, specific storage ${S}_{s}=0.00002\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{-1}$, and specific yield ${S}_{y}=0.2$. The $L=2l=20\text{-}\mathrm{m}$-long horizontal well was simulated as a line sink of uniform strength positioned at $y=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ and $z=15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ between points $x=\u201310\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ and $x=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, where $z$ marks the elevation above the bottom of the aquifer. A pumping rate of $Q=0.02\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{3}/\mathrm{s}$ was applied in calculations. By considering uniform flux distribution along the well screen, the unsteady drawdown in Observation Point A at $x=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, $y=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, and $z=15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ is calculated. The authors rotated the line sink around the middle point into the vertical $xz$-plane and compared the drawdown at angles 0°, 45°, and 90° of inclination. Based on Fig. 7 in Zhan and Zlotnik (2002), the referenced authors concluded that within the selected geometric layout the inclination has little effect and “the drawdown at the rotational axis of this slanted well is insensitive to this angle.” The inclination option of 90° was confirmed through comparison with results by Neuman (1974) or Moench (1997) for vertical wells. This conclusion is confirmed in the forthcoming analysis for the case of wells with uniform wellbore drawdown.

In this section we focus on the horizontal and vertical wells only and extend the simulation options in several ways as follows:

•

The line sinks of uniform strength by Zhan and Zlotnik (2002) are changed to wells of finite drilling radius ${r}_{d}=0.1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ operated under uniform wellbore drawdown condition.

•

The CW-LTP method (see the Appendix) is used as the late time reference tool.

•

The approximate closed-form solution by Renard and Dupuy (1991) for the wellbore is also used as an alternate analytical core ${s}_{1}$ of the LTP simulation.

•

The FLOW software (Szekely 1998) with the multimesh technique is used for numerical FD modeling of the drawdown evolution in and around the horizontal well. The three applied meshes provide sufficiently high vertical and lateral resolution near to the well and low resolution at large distance.

•

The analytical method based on the numerical Laplace inversion and the relevant WT software (Szekely 2021) are used for simulating the unsteady flow to the 20-m-long vertical well.

•

The WF software (Szekely 2021) is utilized to simulate the unsteady flow to the line of short screen vertical wells positioned along the axis of the horizontal well. This approximate replacing well line method provides an alternate analytical solution to the drawdown effect of the horizontal well.

•

The FW method (see the section “Overview of Methods Used in Verifications and Comparative Analyses”) provides additional information based on the replacing well line concept.

•

In WF and FW modeling, the pumping rates of the replacing wells are obtained from the flux distribution along the axis of horizontal well generated by the steady-state CW calculation.

Table 2 summarizes drawdown data calculated at the end of the 1,000-day-long simulation period. Point A is located at $x=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, $y=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$, and depth $z=15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. Next, we provide more information on the multimesh FLOW simulation method.

Calculation site | Models of horizontal well | Models with replacing vertical wells | ||||
---|---|---|---|---|---|---|

CW-LTP | Renard-Dupuy-LTP | FLOW | WT | WF | FW | |

Pumping well | 11.985 | 12.331 | 12.020 | 12.248 | 12.208 | 12.020 |

Point A | 5.460 | — | 5.421 | 5.529 | 5.478 | 5.425 |

The horizontal well operates in the innermost high spatial resolution Mesh III of $100\text{\hspace{0.17em}}\mathrm{m}\times 100\text{\hspace{0.17em}}\mathrm{m}$ area and is modeled with high horizontal transmissivities of ${10}^{7}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{2}/\mathrm{day}$ introduced between nodes along the screen (the hydraulic conductivity of blocks controlling the radial and spherical inflow to the wellbore are unchanged). The water extraction is applied to the series of nodes of the horizontal well exhibiting practically uniform drawdown.

The low-resolution Mesh I of $10\times 10\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{km}}^{2}$ area with no-flow lateral boundaries is used to simulate the response of the infinite aquifer under unsteady flow condition over 1,000 days. Because of the large model area, the simulation resulted in 0.003-m maximum drawdown on the boundary at a distance of 5 km. This causes negligible drawdown impact in the well; thus, the model may properly emulate the infinite flow domain.

The intermediate Mesh II of $1\times 1\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{km}}^{2}$ area and medium resolution provides smooth transition of drawdown and fluxes between Meshes I and III.

In this example, the nodal drawdown of Mesh III is considered as the wellbore drawdown. By using 1-m thickness for the tapped model layer, the late (1,000-day) wellbore drawdown data obtained with the CW-LTP and FLOW methods showed sufficient difference. Gradual reduction of the thickness of the tapped model layer and simultaneous increment of the thickness of the two associated model layers in Mesh III improved the fit between the two model outputs. At 0.344-m thickness of the pumped model Layer 7 the FLOW simulation showed the 12.020-m wellbore drawdown fitting to the result of the CW-LTP modeling. The P curve in Fig. 2 depicts the time drawdown data of the numerical FLOW (solid line) and the analytical WF (dots) simulations. The fit between the two data sets provides independent verification for the FLOW simulation.

The straight tangent line of the late time $s\u2013\mathrm{log}t$ plot for the pumping well represents the CW-LTP solution. The plot confirms that the combined CW-LTP method can be safely used at elapsed time of 1 day and later. The pseudo-steady-state drawdown between 0.001 and 0.1 days is the result of the delayed gravity response (Neuman 1974) and marks the transition from the early time confined to the late time unconfined flow conditions. This feature also appears in Fig. 7 by Zhan and Zlotnik (2002).

In Fig. 2, Curve A with dots exhibits the drawdown evolution in Observation Point A. The solid line shows the drawdown of the horizontal well obtained with the numerical FLOW modeling, whereas the dots mark the WT simulated analytical drawdown caused by the vertical well of a 20-m-long screen.

Fig. 3 shows the late (1,000-day) drawdown contour lines at the depth of the horizontal well in the $100\times 100\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{2}$ inner Mesh III obtained with the FLOW modeling option.

Five simulation models are presented, allowing for calculating the drawdown in Observation Point A away from the well screen(s). Models FLOW, WT, WF, and FW exhibit sufficient difference from the base CW-LTP method regarding the extraction scheme, vertical discretization, and skin parameters. These various options are validated by fitting the late wellbore drawdown to the base version. The last row of Table 2 shows the drawdown at Point A. The drawdown data of observation piezometers or wells depend little on the model selection and show noticeable coherence. For this reason, they supply vital information to estimate the aquifer parameters, especially in case of considerable local (skin, friction) losses around or in the pumping well.

Following Zhan et al. [2001 Eqs. (30) and (37)], the method by Rosa and Carvalho (1989) is used for additional and independent verification of the drawdown in the pumping well at $t=\mathrm{1,000}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{days}$. Applying the averaged specific storage $({S}_{y}+{S}_{s}\xb7b)\xb7{b}^{-1}$ to the equivalent and presumably confined aquifer the calculation yields 12.348-m well drawdown, which is close to the relevant data in Table 2. This coherence confirms the six methods used in the presented comparative analysis.

### Unsteady Flow to the Two-Level Radial Collector Well in an Unconfined Aquifer

In many areas, the horizontal arms of radial collector wells are drilled at several elevations above the bottom of the caisson. The following modeling example assumes a collector well with four arms of equal $l=10\text{-}\mathrm{m}$ length and ${r}_{d}=0.1\text{-}\mathrm{m}$ drilling radius. The arms are drilled from the 2-m-diameter caisson constructed at location $x=y=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. The thickness and hydraulic parameters of the unconfined aquifer of the previous analysis are used in the calculation. The well is also pumped at the same rate of $0.02\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{3}/\mathrm{s}$. The two upper arms of the east–west direction are drilled at a depth of 12.27 m below the initial water table. The lower arms of the north–south direction are positioned at a depth of 17.73 m.

Simulation options of infinite and semi-infinite aquifers are evaluated and the late time results are summarized in Table 3.

Calculation site | Infinite aquifer | Semi-infinite aquifer | |||
---|---|---|---|---|---|

CW-LTP | WF | FLOW | CW | FLOW | |

Caisson | 8.671 | 8.701 | 8.671 | 5.662 | 5.689 |

Point A | 4.350 | 4.351 | 4.321 | 1.547 | 1.544 |

The flow in the infinite aquifer is analyzed first. The CW software is used for calculating the drawdown in the caisson and in the 15-m-deep Piezometer A located at $x=y=20\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. A circular recharge boundary of radius ${R}^{*}=300\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ is assumed in calculations. The CW simulation yielded 6.828- and 2.508-m steady-state drawdown in the caisson and the piezometer, respectively. The unsteady drawdown at elapsed time of 1,000 days was calculated with the combined CW-LTP method. The calculation resulted in 8.671-m drawdown for the caisson and 4.350-m drawdown for the piezometer.

The analytical modeling with the WF software was performed next. In this analysis, each arm is split into five horizontal sections that are replaced with vertical wells. The curves with dots P (pumping well) and A (piezometer) in Fig. 4 show the simulated drawdown obtained for the option of infinite aquifer. The dots over period 0.0001–1,000 days display the analytical drawdown data provided with the WF software. At time of 1,000 days, the WF method resulted in close drawdown values of 8.701 m (caisson) and 4.351 m (piezometer). The results obtained with the well line approximation method fit well to the result of the numerical FLOW simulation for both the pumping well (Curve P) and the piezometer (Curve A).

The multimesh FLOW software is applied for the numerical drawdown simulation. The three-mesh layout of the previous model study is used. In this case for Mesh III, we keep the 0.01-m-thick top layer to emulate the water table but apply 11 model layers of 2.7272-m uniform thickness for the 30-m-thick saturated part of the aquifer. In Piezometer A at elapsed time 1,000 days, the FD simulated drawdown of 4.321 m showed reasonable fit to the CW-LTP calculated value. By contrast, sufficient underestimation was detected in the caisson. The earlier method of decreasing the thickness of the tapped model layer proved to be unsuccessful; therefore, we applied an alternate solution to this problem.

The head correction is based on introducing the extra drawdown caused by radial flow in the vertical plane normal to the axes of the horizontal arms. The nodes of the numerical FLOW model allocated along the arms are given high-conductivity lateral interconnections to simulate the impact of negligible head loss in the screen pipes. The caisson is modeled as a vertical well of drilling radius ${r}_{d}=1\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ open to the formation in the two discharged model layers. It operates at the presumably uniform head in the arms. The effect of head loss between the nodes and the associated arm sections is emulated as the assumed skin loss at the face of the caisson. The virtual skin parameters are estimated via fitting the drawdown data in the caisson obtained with the CW-LTP ad FLOW modeling. For Mesh III, we obtained 7.343-m drawdown in the arms and the 8.671-m target drawdown in the caisson. The solid lines in Fig. 4 at Curve P depict the FLOW simulated unsteady drawdown in the caisson.

For this two-level radial collector well, the hydraulic effect of the straight-line recharge boundary was additionally evaluated. The constant-head or zero drawdown boundary is introduced along the vertical plane $x=\u201350\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. This boundary of the north–south direction extends infinitely along $y$ and simulates the hydraulic effect of infiltration from the long river. The latter flows and provides fixed-head boundary conditions at the western boundary of the $100\times 100\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{2}$ model area of Mesh III shown in Fig. 5. It provides recharge to the well through the whole aquifer thickness. This recharge of this semi-infinite aquifer increases with time and converges to the stabilized value sustaining the steady-state flow conditions at $t\to \infty $.

In this modified simulation example, we keep the base parameters (size, resolution) of the previous multimesh modeling. The western sides of all three square-shaped model areas are positioned along the fixed-head boundary and the midpoints of those sides share the same location of $x=\u201350\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ and $y=0\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$.

Solid lines ${\mathrm{P}}^{*}$ and ${\mathrm{A}}^{*}$ in Fig. 4 exhibit the FLOW simulated drawdown–time curves for the pumping well and observation point. The good fit between the data for the CW and FLOW models at $t=\mathrm{1,000}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{days}$ confirms the negligible effect of the constant-head boundary on the short time drawdown evolution.

Dots at $t=\mathrm{1,000}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{days}$ mark the CW calculated steady-state drawdown for the pumping well, 5.662 m, and at Point A, 1.547 m. The FLOW simulated data (5.689 and 1.544 m) fit well to these values. The plot proves that in the presented riverbank environment, the drawdown close to the river stabilizes after a few tens of days.

Fig. 5 shows the drawdown contour lines at $t=\mathrm{1,000}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{days}$ in Model III at the depths of the arms for the option of recharge boundary at $x=\u201350\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$.

### Steady-State Flow to the Five-Arm Radial Collector Well in a Square-Shaped Aquifer

Bakker et al. (2005) presented the results of simulation of steady-state drawdown in and around a five-arm radial collector well. Multilayer and three-dimensional (3D) analytic element (AEM) methods are used in calculations. Radial flow is considered beyond some distance and 3D flow is simulated near the collector well. The high-conductivity isotropic ($k=150\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}/\mathrm{day}$) aquifer of $b=24\text{-}\mathrm{m}$ thickness is pumped at a rate of $Q=\mathrm{60,000}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{3}/\mathrm{day}$. The $L=60\text{-}\mathrm{m}$-long arms of ${r}_{d}=0.15\text{-}\mathrm{m}$ drilling radius are located at elevation between $z=2.85\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ and $z=3.15\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ above the impervious bottom. A caisson of 3-m radius is considered in calculations. Bakker et al. (2005) used a multilayer model that includes four layers below and seven layers above the screened interval.

In the modeling process, MODFLOW is used, which is the most widely used program in the world for simulating groundwater flow (Szucs et al. 2013). In this modeling system, the MODLFLOW MNW2 package for horizontal and slanted wells (Konikow et al. 2009) can be operated. Although this package was created to simulate horizontal and slanted wells, it can be a useful option to determine the hydraulic conditions of collector wells with horizontal arms. To know the behavior of the MNW2 package, hydraulic conditions calculated by MNW2 are compared with the results of Bakker et al. (2005) and the CW simulations.

The MODFLOW MNW2 simulation model (Konikow et al. 2009) covers a $602\times 602\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{2}$ square-shaped area (Fig. 6) and uses the vertical discretization applied by Bakker et al. (2005). By choosing an area of this size, the center of the collector well is located on the middle of the model area, and it provides symmetry during the simulation. The FD mesh has high lateral resolution ($\mathrm{\Delta}x=\mathrm{\Delta}y=2\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$) to achieve sufficient precision for providing appropriate fit to the selected alternate simulation models. The simulation aimed at calculating the head difference between the top of the aquifer and the caisson. The early test simulations resulted in underestimation of the head difference in comparison with the multilayer (0.64 m) and 3D AEM (0.70 m) solutions by Bakker et al. (2005) and the CW (0.66 m) calculated value. The deviation between the MODFLOW MNW2 and the alternate modeling results was eliminated by means of manual adjustment of the cell-to-well conductance (CWC) parameter of the MNW2 package. The final drawdown contour lines are shown in Fig. 6 where the vertical head difference is set to 0.66 m.

According to Fig. 6, the flow domain may be split into three different subdomains. The inner zone is controlled with the arms of the radial collector well. This 3D flow domain is marked with green contour lines that converge to circles at increasing distance from the center. The blue contour lines of the outer zone are affected by the square-shaped boundary. The intermediate red circular contour lines in the transition zone identify radial flow patterns. The latter feature allows for comparing with the alternate solutions corresponding to radial flow at distance.

The MODFLOW MNW2 modeling yielded discharge rates of arms between 10,963 and $\mathrm{12,167}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{3}/\mathrm{day}$, which contradicts to the radial symmetry of the flow in the transition zone. By contrast, the CW simulation at circular recharge boundary of radius ${R}^{*}=100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ concluded with the equal $\mathrm{12,000}\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{3}/\mathrm{day}$ discharge rate for all five arms. For comparative analysis of inflow distribution along the arms, this option also needs to be applied to the numerical model. This condition, referred to as analytical support, may be achieved by replacing the five-arm radial collector well with five separately pumped horizontal wells without direct hydraulic connection in the caisson. In this case, the inflow distribution of the eastern well along the $x$-axis is expected to be coherent with that of the CW simulation. Fig. 7 shows the CW output (dots) and the MODFLOW MNW2 data (line) in the centers of the 2-m-long sections constituting the 60-m-long arm. The negligible difference between the two solutions confirms the reliability of numerical MODFLOW MNW2 modeling of inflow data along the axes of screen pipes.

## Conclusions

In our work, we dealt with horizontal and radial collector wells, which are one of the most important well types in RBF systems. We aimed to better understand the hydraulic conditions of this kind of important water production facility. In this study, we used various methods for modeling the hydraulic conditions of horizontal and radial collector wells. The following six points summarize our main conclusions:

•

Unlike the vertical wells, the horizontal and radial collector wells at present do not have a widely recognized, safe, multioptional modeling environment and software tools. Therefore, application of several suitable techniques is recommended to derive reasonable model output in both the steady-state and unsteady (transient) flow analyses.

•

The analytical CW software is suitable for calculating the 3D steady-state drawdown in and around the horizontal or radial collector wells by considering homogeneous aquifers and constant-head recharge from circular or straight-line lateral boundaries or from the top. The combined application called the CW-LTP method is used to calculate the long-term unsteady drawdown considering the option of infinite aquifer. These data are accepted as the analytical target values to be reached in model fitting procedures.

•

By introducing the concept of replacing the well line method, the analytical WF and the numerical FW software are used for unsteady simulation of a horizontal well in infinite aquifers. Two more analytical methods (Renard-Dupuy-LTP and WT) are also involved in the verification procedure of solutions for horizontal well.

•

The multimesh FLOW software was used for numerical simulation of the unsteady drawdown variation caused by the horizontal well and the four-arm radial collector well in large, presumably infinite or semi-infinite aquifers. Analytical support is involved to optimize the thickness of the pumped model layer (horizontal well) and the skin parameter (collector well).

•

The MODFLOW MNW2 modeling was successfully applied for simulating the steady-state flow to the five-arm radial collector well. The simulated vertical head difference between the top model layer and the caisson was fitted to the data obtained with three independent calculation methods by means of manual adjustment of the CWC parameter.

•

The MODFLOW MNW2 modeling also yielded inflow distribution along arms, which is coherent with the result obtained by means of the CW software.

In future field studies, the two presented numerical simulation methods can be used for (1) heterogeneous aquifers with irregular, curve-linear boundaries; (2) water table aquifers with drawdown-dependent saturated thickness; and (3) wells with screens operated under the riverbed common in riverbank filtration schemes.

## Appendix. Late Time Drawdown Partitioning Method

Hantush (1964, Fig. 25) analyzed and visualized the unsteady drawdown around radial collector wells of arms’ length $l$ in infinite water table aquifers. He concluded that at distance more than $5l$ “the equation of drawdown for a collector well of at least two laterals on a line can be given by the Theis formula” (Hantush 1964). This means that the drawdown contour lines become practically circular at and beyond the distance ${R}_{c}=5\xb7l$ from the center of the well. For fully penetrating vertical fracture or horizontal trench, the same is demonstrated graphically by Van Tonder et al. (2002).

For single vertical wells operating in infinite confined or unconfined aquifers, the $s\u2013\mathrm{log}t$ plot at late time converges to the straight line and the Theis solution may be approximated by $Q/(4\xb7\pi \xb7T)\xb7\mathrm{ln}[(2.2458\xb7T\xb7t)/({r}^{2}\xb7S)]$, referred to as Jacob’s formula (Cooper and Jacob 1946). Here symbols $s$, $Q$, $T$, $S$, $r$, and $t$ denote the drawdown, discharge rate, transmissivity, storativity, distance, and time, respectively. An arbitrary recharge boundary $R$ is introduced by adding and subtracting the term $\mathrm{ln}(R)$. In this way the previously referenced single term unsteady drawdown formula can be partitioned or decomposed into the steady-state analytical core ${s}_{1}$ and the unsteady part as follows:where

$$s={s}_{1}+\frac{Q}{4\xb7\pi \xb7T}\xb7\mathrm{ln}\frac{2.2458\xb7T\xb7t}{{R}^{2}\xb7S}\phantom{\rule[-0.0ex]{1em}{0.0ex}}\mathrm{at}\phantom{\rule[-0.0ex]{1em}{0.0ex}}r\le R$$

(1)

$${s}_{1}=\frac{Q}{2\xb7\pi \xb7T}\xb7\mathrm{ln}\frac{R}{r}$$

(2)

Eq. (2) is valid at $r<R$. The term ${s}_{1}$ is the Dupuit equation at assumed constant-head circular recharge boundary of radius $R$, while the second term on the right-hand side of Eq. (1) defines the time-variant drawdown at $R$. Similar LTP drawdown approximation may also be introduced for vertical fractures, horizontal drains, or collector wells. It requires replacement of the Dupuit equation with the relevant steady-state solutions $s(x,y,z,R)$ for the selected abstraction scheme considering circular recharge boundary with an influence radius $R$. The following comparative analysis demonstrates and confirms this procedure.

Gringarten et al. (1974) developed analytical methods for two-dimensional (2D) unsteady drawdown in a confined homogeneous aquifer tapped by a fully penetrating vertical fracture of length $2l$. The options of uniform flux and infinite fracture conductivity or uniform head are considered. Table 1 in Gringarten et al. (1974) provides the dimensionless drawdownin the fracture of infinite conductivity. The dimensionless time is defined as

$${s}_{D}({t}_{D})=\frac{s}{\frac{Q}{2\xb7\pi \xb7T}}$$

(3)

$${t}_{D}=\frac{T\xb7t}{S\xb7{l}^{2}}$$

(4)

As the time increases, (1) the gradient stabilizes at $r<R$, (2) the $s\u2013\mathrm{log}t$ graph becomes linear, and (3) the flow patterns converge to the radial flow beyond the distance ${R}_{c}$. At this stage, the ${s}_{D}$ function for the fracture can be approximated with the LTP method.

Citrini (1951) analyzed the 2D steady-state flow to the radial system of symmetrically located $N$ drains fully penetrating the confined aquifer. The drains of zero width and uniform length $l$ extend from the origin and the drawdown are given as (Citrini 1951; Babac and Babac 2009)

$${s}_{1}=\frac{Q}{2\xb7\pi \xb7T}\xb7\mathrm{ln}\frac{R}{{r}_{e}}$$

(5)

The equivalent well radius ${r}_{e}$ is to be calculated as

$${r}_{e}=l\xb7{0.25}^{1/N}$$

(6)

In the case of single drain $N=2$, the flow patterns are identical to the vertical fracture of infinite conductivity. Another formula was presented by Dvoretskii and Yarmakhov (1998, p. 131). Both solutions yield identical results. Dvoretskii and Yarmakhov (1998) also published a closed-form analytical solution for the steady-state-specific (per unit length) flux along the axes of the drains (Dvoretskii and Yarmakhov 1998).

Renard and Dupuy (1991) presented the closed-form formula for the drawdown in a fully penetrating horizontal well positioned in the center of the symmetric 2D flow domain in the form of

$${s}_{1}=\frac{Q}{2\xb7\pi \xb7T}\xb7\mathrm{ln}(\frac{{R}^{*}}{l}+{({\left(\frac{{R}^{*}}{l}\right)}^{2}-1)}^{1/2})$$

(7)

Following Giger (1987), ${R}^{*}$ denotes either the radius of constant-head circular recharge boundary or the ${(A/\pi )}^{1/2}$ parameter, where $A$ is the area of the presumably square-shaped flow domain with zero drawdown at the boundary.

Table 4 compares the analytical (Gringarten et al. 1974, Table 1), the LTP modified (Citrini 1951; Renard and Dupuy 1991) solutions at the following parameters: $T=0.001\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{2}/\mathrm{s}$, $S=0.00001$, $Q=0.001\text{\hspace{0.17em}\hspace{0.17em}}{\mathrm{m}}^{3}/\mathrm{s}$, and $l=10\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$. A value of $R=100\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{m}$ is selected for calculations; this assumed radius of influence yielded the following steady-state ${s}_{1}$ data: 4.7677 m (Citrini model) and 4.7637 m (Renard and Dupuy model). The small deviation from the analytical data verifies the applicability of the LTP transformation for this flow problem.

Dimensionless time, ${t}_{D}$ | Dimensionless drawdown, ${s}_{D}$ | ||
---|---|---|---|

Analytical | Citrini-LTP | Renard-Dupuy-LTP | |

5 | 1.9312 | 1.9023 | 1.8998 |

50 | 3.0634 | 3.0536 | 3.0511 |

500 | 4.2127 | 4.2049 | 4.2024 |

5,000 | 5.3638 | 5.3562 | 5.3537 |

## Data Availability Statement

Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.

## Acknowledgments

The research was carried out within the 2018-1.2.1-NKP-2018-00011 “Drinking Water: Multidisciplinary Assessment of Secure Supply from the Source to the Consumers” project. We would like to thank the reviewers for their useful comments to improve our manuscript.

## References

Babac, D., and P. Babac. 2009. “Wells with horizontal drains.” In

*Theory, practice, calculation examples*. Belgrade, Serbia: Balby International.Bakker, M., V. A. Kelson, and K. H. Luther. 2005. “Multilayer analytic element modeling of radial collector wells.”

*Ground Water*43 (6): 926–934.Buday, T., P. Szűcs, M. Kozák, Z. Püspöki, R. W. McIntosh, E. Bódi, B. Bálint, and K. Bulátkó. 2015. “Sustainability aspects of thermal water production in the region of Hajdúszoboszló-Debrecen, Hungary.”

*Environ. Earth Sci.*73 (7): 1–11. https://doi.org/10.1007/s12665-014-3983-1.Citrini, D. 1951.

*Drenaggi a raggiera e sistemi di pozzi in una falda artesiana; l’ energia elettrica*. Milano, Italy: Istituto di idraulica e costruzioni idrauliche del Politecnico di Milano.Cooper, H. H., and C. E. Jacob. 1946. “A generalized graphical method for evaluating formation constants and summarizing well field history.”

*Am. Geophys. Union Trans.*27 (4): 526–534. https://doi.org/10.1029/TR027i004p00526.Cyr, M. D. 2007. “An aquifer-well coupled model. A refined implementation of wellbore boundary conditions in three-dimensional, heterogeneous formations.” M.S. thesis, Dept. of Civil Engineering, Queen’s Univ. Kingston.

Dragon, K., J. Górski, and M. Skotnicki. 2018. “Removal of natural organic matter and organic micropollutants during riverbank filtration in Krajkowo, Poland.”

*Water*10 (10): 1457. https://doi.org/10.3390/w10101457.Dvoretskii, P. I., and I. G. Yarmakhov. 1998.

*Electromagnetic and hydrodynamic techniques to the probing and performance of oil/gas reservoirs*. Moscow: Nedra Publishing House.Giger, F. 1987. “Low permeability reservoir development using horizontal wells.” In Proc., SPE Joint Low-Permeability Reservoir Symp./Rocky Mountain Regional Meeting, 18–19. Richardson, TX: OnePetro.

Gringarten, A. C., H. R. Ramey, and R. Raghavan. 1974. “Unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture.”

*Soc. Pet. Eng. J.*14 (4): 347–360. https://doi.org/10.2118/4051-PA.Grischek, T., D. Schoenheinz, P. Eckert, and C. Ray. 2012. “Sustainability of riverbank filtration—Examples from Germany.” In

*Groundwater quality sustainability*. 213–227. Washington, DC: International Association of Hydrogeologist.Hantush, M. S. 1964. “Hydraulics of wells.” In

*Advances in hydroscience*. 282–432. New York: Academic.Hiscock, K. M., and T. Grischek. 2002. “Attenuation of groundwater pollution by bank filtration.”

*J. Hydrol.*266 (3): 139–144. https://doi.org/10.1016/S0022-1694(02)00158-0.Hunt, H. 2003. “American experience in installing horizontal collector wells.” In

*Riverbank filtration, Improving water-source quality*, 29–34. Washington, DC: Water Science and Technology Library.Konikow, L. F., G. Z. Hornberger, K. J. Halford, R. T. Hanson, and A. W. Harbaugh. 2009.

*Revised multi-node well (MNW2) package for MODFLOW ground-water flow model. Techniques and methods 6–A30*. Washington, DC: USGS.Milojevic, M. 1963. “Radial collector wells adjacent to the river bank.”

*J. Hydraul. Div.*89 (6): 133–151. https://doi.org/10.1061/JYCEAJ.0000950.Moench, A. F. 1997. “Flow to a well of finite diameter in a homogeneous, anisotropic water table aquifer.”

*Water Resour. Res.*33 (6): 1397–1407. https://doi.org/10.1029/97WR00651.Morozov, P. E. 2018. “Steady fluid flow to a radial system of horizontal wells.”

*J. Appl. Mech. Tech. Phys.*59 (2): 273–280. https://doi.org/10.1134/S0021894418020104.Neuman, S. P. 1974. “Effects of partial penetration on flow in unconfined aquifers considering delayed gravity response.”

*Water Resour. Res.*10 (2): 303–312. https://doi.org/10.1029/WR010i002p00303.Ray, C., T. Grischek, J. Schubert, J. Z. Wang, and T. F. Speth. 2002. “A perspective of riverbank filtration.”

*J. Am. Water Works Assoc.*94 (4): 149–160. https://doi.org/10.1002/j.1551-8833.2002.tb09459.x.Renard, G., and J. M. Dupuy. 1991. “Formation damage effects on horizontal-well flow efficiency.”

*J. Pet. Technol.*43 (7): 786–869. https://doi.org/10.2118/19414-PA.Rosa, A. J., and S. Carvalho. 1989. “A mathematical model for pressure evaluation in an infinite-conductivity horizontal well.”

*SPE Form. Eval.*4 (4): 559–566. https://doi.org/10.2118/15967-PA.Sandhu, C., and T. Grischek. 2012. “Riverbank filtration in India-using ecosystem services to safeguard human health.”

*Water Sci. Technol.*12 (6): 783–790.Szekely, F. 1998. “Windowed spatial zooming in finite difference ground water flow models.”

*Ground Water*36 (5): 718–721. https://doi.org/10.1111/j.1745-6584.1998.tb02188.x.Szekely, F. 2020.

*Manual of the software package WELLSD version 1.3*. Somerset, England: WellsD.Szekely, F. 2021.

*Integrated well flow modeling*. Cleveland: LAMBERT Academic.Szucs, P., F. Civan, and M. Virag. 2006. “Applicability of the most frequent value method in groundwater modeling.”

*Hydrogeol. J.*14 (1–2): 31–43. https://doi.org/10.1007/s10040-004-0426-1.Szucs, P., F. Szekely, and B. Zakanyi. 2013. “Comparison of analytical and numerical approaches for simulating groundwater flow to multi screen wells.”

*Carpathian J. Earth Environ. Sci.*8 (2): 69–76.Van Tonder, G., I. Bardenhagen, K. Riemann, J. van Bosch, P. Dzanga, and Y. Xu. 2002.

*Manual on pumping test analysis in fractured-rock aquifers*. Berlin: Springer.Zhan, H., L. V. Wang, and E. Park. 2001. “On the horizontal-well pumping tests in anisotropic confined aquifers.”

*J. Hydrol.*252 (1–4): 37–50. https://doi.org/10.1016/S0022-1694(01)00453-X.Zhan, H., and V. A. Zlotnik. 2002. “Groundwater flow to a horizontal or slanted well in an unconfined aquifer.”

*Water Resour. Res.*38 (7): 1–13. https://doi.org/10.1029/2001WR000401.## Information & Authors

### Information

#### Published In

#### Copyright

This work is made available under the terms of the Creative Commons Attribution 4.0 International license, https://creativecommons.org/licenses/by/4.0/.

#### History

**Received**: Feb 10, 2021

**Accepted**: Jul 26, 2021

**Published online**: Sep 29, 2021

**Published in print**: Dec 1, 2021

**Discussion open until**: Feb 28, 2022

### Authors

## Metrics & Citations

### Metrics

### Citations

#### Download citation

If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download.