Approximate analysis of three-dimensional groundwater flow toward a radial collector well in a finite-extent unconfined aquifer

This study develops a three-dimensional (3-D) mathematical model for describing transient hydraulic head distributions due to pumping at a radial collector well (RCW) in a rectangular confined or unconfined aquifer bounded by two parallel streams and no-flow boundaries. The streams with low-permeability streambeds fully penetrate the aquifer. The governing equation with a point-sink term is employed. A first-order free surface equation delineating the water table decline induced by the well is considered. Robin boundary conditions are adopted to describe fluxes across the streambeds. The head solution for the point sink is derived by applying the methods of finite integral transform and Laplace transform. The head solution for a RCW is obtained by integrating the point-sink solution along the laterals of the RCW and then dividing the integration result by the sum of lateral lengths. On the basis of Darcy’s law and head distributions along the streams, the solution for the stream depletion rate (SDR) can also be developed. With the aid of the head and SDR solutions, the sensitivity analysis can then be performed to explore the response of the hydraulic head to the change in a specific parameter such as the horizontal and vertical hydraulic conductivities, streambed permeability, specific storage, specific yield, lateral length, and well depth. Spatial head distributions subject to the anisotropy of aquifer hydraulic conductivities are analyzed. A quantitative criterion is provided to identify whether groundwater flow at a specific region is 3-D or 2-D without the vertical component. In addition, another criterion is also given to allow for the neglect of vertical flow effect on SDR. Conventional 2-D flow models can be used to provide accurate head and SDR predictions if satisfying these two criteria.


Introduction
The applications of a radial collector well (RCW) have received much attention in the aspects of water resource supply and groundwater remediation since rapid advances in drilling technology.An average yield for the well approximates 27 000 m 3 day −1 (Todd and Mays, 2005).As compared to vertical wells, RCWs require less operating cost, produce smaller drawdown, and have better efficiency of withdrawing water from thin aquifers.In addition, RCWs can extract water from an aquifer underlying obstacles such as buildings, but vertical wells cannot.Recently, Huang et al. (2012) reviewed semi-analytical and analytical solutions associated with RCWs.Since then, Yeh and Chang (2013) provided a valuable overview of articles associated with RCWs.
A variety of analytical models involving a horizontal well, a specific case of a RCW with a single lateral, in aquifers were developed (e.g., Park and Zhan, 2003;Hunt, 2005;Anderson, 2013).The flux along the well screen is commonly assumed to be uniform.The equation describing three-dimensional (3-D) flow is used.Kawecki (2000) developed analytical solutions of the hydraulic heads for the early linear flow perpendicular to a horizontal well and late pseudo-radial flow toward the middle of the well in confined aquifers.They also developed an approximate solution for unconfined aquifers on the basis of the head solution and an unconfined flow modification.The applicability of the approximate solution was later evaluated in comparison with a finite difference solution developed by Kawecki and Al-Subaikhy (2005).Zhan et al. (2001) presented an analytical solution for drawdown induced by a horizontal well in confined aquifers and compared the difference in the type curves based on the well and a vertical well.Zhan 2002) developed a semi-analytical solution of drawdown due to pumping from a nonvertical well in an unconfined aquifer accounting for the effect of instantaneous drainage or delayed yield when the free surface declines.They discussed the influences of the length, depth, and inclination of the well on temporal drawdown distributions.Park and Zhan (2002) developed a semi-analytical drawdown solution considering the effects of a finite diameter, the wellbore storage, and a skin zone around a horizontal well in anisotropic leaky aquifers.They found that those effects cause significant change in drawdown at an early pumping period.Zhan and Park (2003) provided a general semi-analytical solution for pumping-induced drawdown in a confined aquifer, an unconfined aquifer on a leaky bottom, or a leaky aquifer below a water reservoir.Temporal drawdown distributions subject to the aquitard storage effect were compared with those without that effect.Sun and Zhan (2006) derived a semianalytical solution of drawdown due to pumping at a horizontal well in a leaky aquifer.A transient 1-D flow equation describing the vertical flow across the aquitard was considered.The derived solution was used to evaluate the Zhan and Park (2003) solution, which assumed steady-state vertical flow in the aquitard.
Sophisticated numerical models involved in RCWs or horizontal wells were also reported.Steward (1999) applied the analytic element method to approximate 3-D steady-state flow induced by horizontal wells in contaminated aquifers.They discussed the relation between the pumping rate and the size of a polluted area.Chen et al. (2003) utilized the polygon finite difference method to deal with three kinds of seepagepipe flows including laminar, turbulent, and transitional flows within a finite-diameter horizontal well.A sandbox experiment was also carried out to verify the prediction made by the method.Mohamad and Rushton (2006) used MODFLOW to predict flows inside an aquifer, from the aquifer to a horizontal well, and within the well.The predicted head distributions were compared with field data measured in Sarawak, Malaysia.Su et al. (2007) used software TOUGH2 based on the integrated finite difference method to handle irregular configurations of several laterals of two RCWs installed beside the Russian River, Forestville, California, and analyzed pumping-induced unsaturated regions beneath the river.Lee et al. (2012) developed a finite element solution with triangle elements to assess whether the operation of a RCW near Nakdong River in South Korea can induce riverbank filtration.They concluded that the well can be used for sustainable water supply at the study site.In addition, Rushton and Brassington (2013a) extended Mohamad and Rushton (2006) study by enhancing the Darcy-Weisbach formula to describe frictional head lose inside a horizontal well.The spatial distributions of predicted flux along the well revealed that the flux at the pumping end is 4 times the magnitude of that at the far end.Later, Rushton and Brassington (2013b) applied the same model to a field experiment at the Seton Coast, northwest England.
Well pumping in aquifers near streams may cause groundwater-surface water interactions (e.g., Rodriguez et al., 2013;Chen et al., 2013;Zhou et al., 2013;Exner-Kittridge et al., 2014;Flipo et al., 2014;Unland et al., 2014).The stream depletion rate (SDR), commonly used to quantify stream water filtration into the adjacent aquifer, is defined as the ratio of the filtration rate to a pumping rate.The SDR ranges from zero to a certain value, which could be equal to or less than unity (Zlotnik, 2004).Tsou et al. (2010) developed an analytical solution of the SDR for a slanted well in confined aquifers adjacent to a stream treated as a constanthead boundary.They indicated that a horizontal well parallel to the stream induces the steady-state SDR of unity more quickly than a slanted well.Huang et al. (2011) developed an analytical SDR solution for a horizontal well in unconfined aquifers near a stream regarded as a constant-head boundary.Huang et al. (2012) provided an analytical solution for SDR induced by a RCW in unconfined aquifers adjacent to a stream with a low-permeability streambed under the Robin condition.The influence of the configuration of the laterals on temporal SDR and spatial drawdown distributions was analyzed.Recently, Huang et al. (2014) gave an exhaustive review on analytical and semi-analytical SDR solutions and classified these solutions into two categories.One group involved 2-D flow toward a fully-penetrating vertical well according to aquifer types and stream treatments.The other group included the solutions involving 3-D and quasi-3-D flows according to aquifer types, well types, and stream treatments.
At present, existing analytical solutions associated with flow toward a RCW in unconfined aquifers have involved laborious calculation (Huang et al., 2012) and predicted approximate results (Hantush and Papadopoulos, 1962).The Huang et al. (2012) solution involves numerical integration of a triple integral in predicting the hydraulic head and a quintuple integral in predicting SDR.The integrand is expressed in terms of an infinite series expanded by roots of nonlinear equations.The integration variables are related to those roots.The application of their solution is therefore limited to those who are familiar with numerical methods.In addition, the accuracy of the Hantush and Papadopoulos (1962) solution is limited to some parts of a pumping period; that is, it gives accurate drawdown predictions at early and late times but divergent ones at middle times.
The objective of this study is to present new analytical solutions of the head and SDR, which overcome the abovementioned limitations, for 3-D flow toward a RCW.A mathematical model is built to describe 3-D spatiotemporal hydraulic head distributions in a rectangular unconfined aquifer bounded by two parallel streams and by the no-flow stratums in the other two sides.The flux across the well screen is assumed to be uniform along each of the laterals.The assumption is valid for a short lateral within 150 m verified by agreement on drawdown observed in field experiments and predicted by existing analytical solutions (Huang et  2011, 2012).The streams fully penetrate the aquifer and connect the aquifer with low-permeability streambeds.The model for the aquifer system with two parallel streams can be used to determine the fraction of water filtration from the streams and solve the associated water right problem (Sun and Zhan, 2007).The transient 3-D groundwater flow equation with a point-sink term is considered.The first-order free surface equation is used to describe water table decline due to pumping.Robin boundary conditions are adopted to describe fluxes across the streambeds.The head solution for a point sink is derived by the methods of Laplace transform and finite integral transform.The analytical head solution for a RCW is then obtained by integrating the point-sink solution along the well and dividing the integration result by the total lateral length.The RCW head solution is expressed in terms of a triple series expanded by eigenvalues, which can be obtained by a numerical algorithm such as Newton's method.
On the basis of Darcy's law and the RCW head solution, the SDR solution can then be obtained in terms of a double series with fast convergence.With the aid of both solutions, the sensitivity analysis is performed to investigate the response of the hydraulic head to the change in each of aquifer's parameters.Spatial head distributions subject to the anisotropy of aquifer hydraulic conductivities are analyzed.The influences of the vertical flow and well depth on temporal SDR distributions are investigated.Moreover, temporal SDR distributions induced by a RCW and a fully penetrating vertical well in confined aquifers are also compared.A quantitative criterion is provided to identify whether groundwater flow at a specific region is 3-D or 2-D without the vertical compo-nent.In addition, another criterion is also given to judge the suitability of neglecting the vertical flow effect on SDR.

Mathematical model
Consider a RCW in a rectangular unconfined aquifer bounded by two parallel streams and no-flow stratums as illustrated in Fig. 1.The symbols for variables and parameters are defined in Table 1.The origin of the Cartesian coordinate is located at the lower left corner.The aquifer domain falls in the range of 0 ≤ x ≤ w x , 0 ≤ y ≤ w y , and −H ≤ z ≤ 0.
The RCW consists of a caisson and several laterals, each of which extends with length L k and counterclockwise with angle θ k , where k ∈ 1, 2, . ..N and N is the number of laterals.
The caisson is located at (x 0 , y 0 ), and the surrounding laterals are at z = −z 0 .
First, a mathematical model describing 3-D flow toward a point sink in the aquifer is proposed.The equation describing 3-D hydraulic head distribution h(x, y, z, t) is expressed as where δ( ) is the Dirac delta function, the second term on the right-hand side (RHS) indicates the point sink, and Q is positive for pumping and negative for injection.The first term on the RHS of Eq. (1) depicts aquifer storage release based on the concept of effective stress proposed by Terzaghi (see for example, Bear, 1979;Charbeneau, 2000), which is valid under the assumption of constant total stress.By choosing the water table as a reference datum where the elevation head is set to zero, the initial condition can therefore be denoted as Note that Eq. ( 2) introduces a negative hydraulic head for pumping, and the absolute value of the head equals drawdown.
The aquifer boundaries at x = 0 and x = w x are considered to be impermeable and thus expressed as ∂h/∂x = 0 at x = 0 (3) and Streambed permeability is usually less than the adjacent aquifer formation.The fluxes across the streambeds can be described by Robin boundary conditions as Aquifer hydraulic conductivities in x, y, and z directions, respectively (K 1 , K 2 ) Hydraulic conductivities of streambeds 1 and 2, respectively The number of laterals Q Pumping rate of point sink or radial collector well p Laplace parameter Shortest horizontal distance between the far end of lateral and aquifer lateral boundary S s , S y Specific storage and specific yield, respectively t Time since pumping t (K y t)/ S s y 2 0 w x , w y Aquifer widths in x and y directions, respectively w x ,w y w x /y 0 ,w y /y 0 X n Equaling X m,n defined in Eq. ( 39) with α m = 0 X n,k Defined in Eq. (45) x, y, z Cartesian coordinate system x, y, z x/y 0 , y/y 0 , z/H x k Coordinate x of the far end of the kth lateral x 0 , y 0 , z 0 Location of center of RCW x 0 , y 0 , z 0 x 0 /y 0 , 1, z 0 /H x 0 ,y 0 ,z 0 Location of point sink x 0 , y 0 , z 0 x 0 /y 0 ,y 0 /y 0 ,z 0 /H α m m π/w x β n

Roots of Eq. (19) φ n
Equaling φ m, n defined in Eq. ( 32) with Roots of Eqs. ( 40) and ( 41), respectively Equaling µ m,n,0 defined in Eq. ( 36) with α m = 0 ν n,i Equaling ν m,n,i defined in Eq. ( 37 The free surface equation describing the water table decline is written as Neuman (1972) indicated that the effect of the secondorder terms in Eq. ( 7) can generally be ignored in developing analytical solutions.Equation ( 7) is thus linearized by neglecting the quadratic terms, and the position of the water table is fixed at the initial condition (i.e., z = 0).The result is written as Notice that Eq. ( 8) is applicable when the conditions |h| /H ≤ 0.1 and |∂h/∂x| + |∂h/∂y| ≤ 0.01 are satisfied.These two conditions have been studied and verified by simulations in, for example, Nyholm et al. (2002), Goldscheider and Drew (2007) and Yeh et al. (2010).Nyholm et al. (2002) achieved agreement on drawdown measured in a field pumping test and predicted by MODFLOW, which models flow in the study site as confined behavior because of |h| /H ≤ 0.1 in the pumping well.Goldscheider and Drew (2007) revealed that pumping drawdown predicted by the analytical solution by Neuman (1972) based on Eq. ( 8) agrees well with that obtained in a field pumping test.In addition, Yeh et al. (2010) also achieved agreement on the hydraulic head predicted by their analytical solution based on Eq. ( 8), their finite difference solution based on Eq. ( 7) with ∂h/∂y = 0, and the Teo et al. ( 2003) solution derived by applying the perturbation technique to deal with Eq. ( 7) with ∂h/∂y = 0 when |h| /H = 0.1 and |∂h/∂x| = 0.01 (i.e., α = 0.1 and ϕ |∂/∂x| = 0.01 at x = 0 in Yeh et al., 2010, Fig. 5a).On the other hand, the bottom of the aquifer is considered as a no-flow boundary condition denoted as Define dimensionless variables as h = (K y H h)/Q, t = (K y t)/S s y 2 0 , x = x/y 0 , y = y/y 0 , z = z/H , x 0 = x 0 /y 0 , y 0 = y 0 /y 0 , z 0 = z 0 /H , w x = w x /y 0 , and w y = w y /y 0 , where the overbar denotes a dimensionless symbol, H is the initial aquifer thickness, and y 0 is a distance between stream 1 and the center of the RCW.On the basis of the definitions, Eq. ( 1) can be written as where κ x = K x /K y and κ z = K z y 2 0 /(K y H 2 ).Similarly, the initial and boundary conditions are expressed as ∂h/∂y − κ 1 h = 0 at y = 0, ( 14) where

Head solution for point sink
The model, Eqs. ( 10)-( 17), reduces to an ordinary differential equation (ODE) with two boundary conditions in terms of z after taking Laplace transform and finite integral transform.The former transform converts h(x, y, z, t) into ĥ(x, y, z, p), δ(x −x 0 )δ(y −y 0 )δ(z−z 0 ) in Eq. ( 10) into δ(x − x 0 ) δ(y − y 0 )δ(z − z 0 )/p, and ∂h/∂t in Eqs. ( 10) and ( 16) into p ĥ − h t=0 , where p is the Laplace parameter, and the second term, the initial condition in Eq. ( 11), equals zero (Kreyszig, 1999).The transformed model becomes a boundary value problem written as with boundary conditions ∂ ĥ/∂x = 0 at x = 0 and x = w x , ∂ ĥ/∂y − κ 1 ĥ = 0 at y = 0, ∂ ĥ/∂y + κ 2 ĥ = 0 at y = w y , ∂ ĥ/∂z = −pγ ĥ/κ z at z = 0, and ∂h/∂z = 0 at z = −1.We then apply finite integral transform to the problem.One can refer to Appendix A for its detailed definition.The transform converts ĥ(x, y, z, p) in the problem into h(α m β n zp), δ(x − x 0 ) δ(y − y 0 ) in Eq. ( 18) into cos(α m x 0 )K(y 0 ), and κ x ∂ 2 ĥ/∂x 2 + ∂ 2 ĥ/∂y 2 in Eq. ( 18) into − κ x α 2 m + β 2 n h, where (m, n) ∈ 1, 2, 3, . ..∞, α m = m π/w x , K y 0 is defined in Eq. (A2) with y = y 0 , and β n represents eigenvalues equaling the roots of the following equation as (Latinopoulos, 1985) tan The method to determine the roots is discussed in Sect.2.3.In turn, Eq. ( 18) becomes a second-order ODE de- with two boundary conditions denoted as and Equation ( 20) can be separated into two homogeneous ODEs as and where h a and h b , respectively, represent the heads above and below z = −z 0 , where the point sink is located.Two continuity requirements should be imposed at z = −z 0 .The first is the continuity of the hydraulic head denoted as The second describes the discontinuity of the flux due to point pumping represented by the Dirac delta function in Eq. ( 20).It can be derived by integrating Eq. ( 20) from Solving Eqs. ( 23) and ( 24) simultaneously with Eqs. ( 21), ( 22), (25), and (26) yields the Laplace-domain head solution as and where a, b, and c are arguments.Taking the inverse Laplace transform and finite integral transform to Eq. ( 28) results in Eq. ( 31).One is referred to Appendix B for the detailed derivation.A time-domain head solution for a point sink is therefore written as where φ n and X n equal φ m,n and X m,n with α m = 0, respectively, and the eigenvalues λ 0 and λ i are, respectively, the roots of the following equations tan The determination for those eigenvalues is introduced in the next section.Notice that the solution consists of a simple series expanded in β n , double series expanded in β n and λ i (or α m and β n ), and triple series expanded in α m , β n , and λ i .Application of Newton's method with proper initial guesses to determine the eigenvalues β n , λ 0 , and λ i has been proposed by Huang et al. (2014) and is briefly introduced herein.
The eigenvalues are situated at the intersection points of the left-hand side (LHS) and RHS functions of Eq. ( 19) for β n , Eq. ( 40) for λ 0 , and Eq. ( 41) for λ i .Hence, the initial guesses for β n are considered as β v −δ if β v > (κ 1 κ 2 ) 0.5 and as β v +δ if β v < (κ 1 κ 2 ) 0.5 , where β v = (2n − 1) π/(2w y ) and δ is a chosen small value such as 10 −8 for avoiding being right at the vertical asymptote.In addition, the guess for λ 0 can be formulated as where the RHS second term represents the location of the vertical asymptote derived by letting the denominator of the RHS function in Eq. ( 40) to be zero and solving λ 0 in the resultant equation.Moreover, the guessed value for λ i is (2 i − 1)π/2 + δ.

Head solution for radial collector well
The lateral of RCW is approximately represented by a line sink composed of a series of adjoining point sinks.The locations of these point sinks are expressed in terms of (x 0 + l cos θ , y 0 + l sin θ , z 0 ) where (x 0 , y 0 , z 0 ) = (x 0 /y 0 , 1, z 0 /H ) is the central of the lateral, and l is a variable to define different locations of the point sink.The solution of head h w (x, y, z, t) for a lateral can therefore be derived by substituting x 0 = x 0 + l cos θ, y 0 = 1 + l sin θ , and z 0 = z 0 into the point-sink solution, Eq. ( 30), then by integrating the resultant solution to l, and finally by dividing the integration result into the sum of lateral lengths.The derivation can be denoted as where L k = L k /y 0 is the kth dimensionless lateral length.Note that the integration variable l (i.e., x 0 and y 0 ) appears only in X n and X m,n in Eq. ( 31).The integral in Eq. ( 43) can thus be done analytically by integrating X n and X m,n with respect to l.After the integration, Eq. ( 43) can be expressed as where is defined by Eqs. ( 31)-( 38), and X n and X m,n in Eq. ( 31) are replaced, respectively, by and with where X = x 0 + L k cos θ k and Y = 1 + L k sin θ k .Notice that Eq. ( 45) is obtained by substituting α m = 0 into Eq.( 46).When θ k = 0 or π , Eq. ( 45) reduces to Eq. ( 49) by applying L'Hospital's rule.

SDR solution for radial collector well
On the basis of Darcy's law and the head solution for a RCW, the SDR from streams 1 and 2 can be defined, respectively, as Again, the double integrals in both equations can be done analytically.Notice that the series term of 2 ∞ m=1 φ m,n X m,n cos (α m x) in Eq. ( 31) disappears due to the consideration of Eqs. ( 3) and ( 4) and the integration with respect to x in Eqs. ( 50) and ( 51) when deriving the SDR solution.The SDR 1 and SDR 2 are therefore expressed in terms of double series and given below: [κ z λ 0 cosh (z 0 λ 0 ) + p 0 γ sinh(z 0 λ 0 )], where 36) with α m = 0; ν n,i = ν m,n,i in Eq. ( 37) with α m = 0; X n,k is defined in Eq. ( 45) for θ k = 0 or π and Eq. ( 49) for θ k = 0 or π; and λ 0 and λ i are the roots of Eqs. ( 40) and ( 41) with α m = 0.

Confined aquifer of infinite extent
The head solution introduced in Sect.2.6.1 is applicable to spatiotemporal head distributions in confined aquifers of infinite extent before the lateral boundary effect comes.Wang and Yeh (2008) indicated that the time can be quantified, in our notation, as t = R 2 S s /(16K y ) (i.e., t = R 2 /(16y 2 0 ) for dimensionless time), where R is the shortest distance between a RCW and aquifer lateral boundary.Prior to the time, the present head solution with N = 1 for a horizontal well in a confined aquifer gives very close results given in Zhan et al. (2001).

Unconfined aquifer of infinite extent
Prior to the beginning time mentioned in Sect.2.6.2, the absolute value calculated by the present head solution, Eqs. ( 44) with N = 1, represents drawdown induced by a horizontal well in unconfined aquifers of infinite extent.The calculated drawdown should be close to that of the solution from Zhan and Zlotnik (2002) for the case of the instantaneous drainage from water table decline.

Unconfined aquifer of semi-infinite extent
When κ 1 → ∞ (i.e., b 1 = 0), Eq. ( 14) reduces to the Dirichlet condition of h = 0 for stream 1 in the absence from a lowpermeability streambed, and Eq. ( 19) becomes tan β n w y = −β n /κ 2 .In addition, the boundary effect occurring at the other three sides of the aquifer can be neglected prior to the beginning time.Moreover, when N = 1 and θ 1 = 0, a RCW can be regarded as a horizontal well parallel to stream 1.Under these three conditions, the present head and SDR predictions are close to those in Huang et al. (2011), the head solution of which agrees well with measured data from a field experiment executed by Mohamed and Rushton (2006).On the other hand, before the time when the boundary effect occurs at x = 0, x = w x , and y = w y , the present head and SDR solutions for a RCW give close predictions to those in Huang et al. (2012), the head and SDR solutions of which agree well with observation data taken from two field experiments carried out by Schafer (2006) and Jasperse (2009), respectively.

Sensitivity analysis
The hydraulic parameters determined from field observed data are inevitably subject to measurement errors.Consequently, head predictions from the analytical model have uncertainty due to the propagation of measurement errors.Sensitivity analysis can be considered as a tool of exploring the response of the head to the change in a specific parameter (Zheng and Bennett, 2002).One may define the normalized sensitivity coefficient as where S i,t is the normalized sensitivity coefficient for the ith parameter at time t, and P i represents the magnitude of the ith parameter.Equation ( 63) can be approximated as where P i is an increment chosen as 10 −3 P i (Yeh et al., 2008).

Results and discussion
This section demonstrates head and SDR predictions and explores some physical insights regarding flow behavior.In Sect.3.1, equipotential lines are drawn to identify 3-D or 2-D flow without the vertical flow at a specific region.In Sect.3.2, the influence of anisotropy on spatial head and temporal SDR distributions is studied.In Sect.3.3, the sensitivity analysis is performed to investigate the response of the head to the change in each hydraulic parameter.In Sect.3.4, the effects of the vertical flow and well depth on temporal SDR distributions for confined and unconfined aquifers are investigated.For conciseness, we consider a RCW with two laterals with N = 2, L 1 = L 2 = 0.5, θ 1 = 0, and θ 2 = π .The well can be viewed as a horizontal well parallel to streams 1 and 2. The default values for the other dimensionless parameters are w x = w y = 2, γ = 100, x 0 = 1, y 0 = 1, z 0 = 0.5, κ x = κ z = 1, and κ 1 = κ 2 = 20.

Identification of 3-D or 2-D flow at observation point
Most existing models assume 2-D flow by neglecting the vertical flow for pumping at a horizontal well (e.g., Mohamed and Rushton, 2006;Haitjema et al., 2010).The head distributions predicted by those models are inaccurate if an observation point is close to the region where the vertical flow prevails.Figure 2 demonstrates the equipotential lines predicted by the present solution for a horizontal well in an unconfined aquifer for x 0 = 10, w x = w y = 20, and κ z = 0.1, 1, and 10.The well is located at 9.5 ≤ x ≤ 10.5, y = 1, and z = 0.5 as illustrated in the figure.The equipotential lines are based on steady-state head distributions plotted by Eq. ( 44) with y = 1 and t = 10 7 .When κ z = 0.1, in the range of 10 ≤ x ≤ 13.66, the contours of the hydraulic head are in a curved path, and the flow toward the well is thus slanted.Moreover, the range decreases to 10 ≤ x ≤ 11.5 when κ z = 1 and to 10 ≤ x ≤ 10.82 when κ z = 10.Beyond these ranges, the head contours are nearly vertical, and the flow is essentially horizontal.Define d = d/y 0 as a shortest dimensionless horizontal distance between the well and a nearest location of only horizontal flow.The d is therefore chosen as 3.16, 1, and 0.32 for the cases of κ z = 0.1, 1, and 10, respectively.Substituting (κ z , d) = (0.1, 3.16), (1, 1), and (10, 0.32) into κ z d 2 leads to about unity.We may therefore conclude that the vertical flow at an observation point is negligible if its location is beyond the range of d< √ 1/κ z (i.e., d < H K y /K z ) for thin aquifers, an observation point far from the well, and/or a small ratio of K y /K z .

Anisotropy analysis of hydraulic head and stream depletion rate
Previous articles have seldom analyzed flow behavior for anisotropic aquifers, i.e., κ x (K x /K y ) = 1.Head predictions based on the models, developed for isotropic aquifers, will be inaccurate if κ x = 1.Consider w x = w y = 2, t = 10 7 for steady-state head distributions, and a RCW with parison with the case of κ x = 1.In Fig. 3b, the contours exhibit smooth curves in the strip regions of 1 ≤ y ≤ 1.45 for the case of κ x = 10 and 1 ≤ y ≤ 1.2 for the case of κ x = 50.For the region of y ≥ 1.45, the predicted heads for both cases agree well, and all the contour lines are parallel, indicating that the flow is essentially unidirectional.Substituting (κ x , y) = (10, 1.45) and (50, 1.2) into κ x (y −1) 2 results in a value of about 2. Accordingly, we may draw the conclusion that plots from the inequality of κ x (y − 1) 2 ≤ 2 indicate the strip region for κ x being greater than 10.Some existing models assuming 2-D flow in a vertical plane with neglecting the flow component along a horizontal well give accurate head predictions beyond the region (e.g., Anderson, 2000Anderson, , 2003;;Kompani-Zare et al., 2005).
Aquifers with K y H ≥ 10 3 m 2 day −1 can efficiently produce plenty of water from a well.RCWs usually operate with Q ≤ 10 5 m 3 day −1 for field experiments (e.g., Schafer, 2006;Jasperse, 2009).We therefore define significant dimensionless head drop as h > 10 −5 (i.e., |h| > 1 mm).The anisotropy of κ x <1 produces the drop in the strip areas of 1 ≤ x ≤ 1.48 for the case of κ x = 10 −3 in Fig. 3c and 1 ≤ x ≤ 1.32 for the case of κ x = 10 −4 in Fig. 3d.Substituting (κ x , x) = (10 −3 , 1.48) and (10 −4 , 1.32) into (x − x 0 − L 1 ) 2 /κ x approximates 52.9.This result leads to the conclusion that the area can be determined by the inequalities of (x −x 0 −L 1 ) 2 ≤ 52.9κ x and (x −x 0 +L 2 ) 2 ≤ 52.9κ x for any value of κ x in the range κ x <1.For a RCW with irregular lateral configurations, the inequalities become (x −max x k ) 2 ≤ 52.9κ x and (x − minx k ) 2 ≤ 52.9κ x , where x k is coordinate x of the far end of the kth lateral.The conclusion applies in principle to reduction in grid points for numerical solutions based on finite difference methods or finite element methods.On the other hand, we have found that Eq. ( 52) or (53) with various κ x predicts the same temporal SDR distribution (not shown), indicating that the SDR is independent of κ x .

Sensitivity analysis of hydraulic head
Consider an unconfined aquifer of H = 20 m and w x = w y = 800 m with a RCW having two laterals of L 1 = L 2 = 50 m, θ 1 = 0, and θ 2 = π and two piezometers installed at point A of (400, 340, −10 m) and point B of (400, 80, −10 m) illustrated in Fig. 4. As discussed in Sect.3.1, the temporal head distribution at point A exhibits the unconfined behavior in Fig. 4a because of κ z d 2 <1 while at point B displays the confined one in Fig. 4b due to κ z d 2 >1.The sensitivity analysis is conducted with the aid of Eq. ( 64) to observe head responses at these two piezometers to the change in each of K x , K y , K z , S s , S y , K 1 , L 1 and z 0 .The temporal distribution curves of the normalized sensitivity coefficients for those eight parameters are shown in Fig. 4a for point A and 4b for point B when m, and z 0 = 10 m.The figure demonstrates that the hydraulic heads at both piezometers are most sensitive to the change in K y , second-most sensitive to the change in K x , and third-most sensitive to the change in S y , indicating that K y , K x , and S y are the most crucial factors in designing a pumping system.This figure also shows that the heads at point A is sensitive to the change in S s at the early period of 4 × 10 −3 day <t<10 −1 day but at point B is insensitive to the change over the entire period.In addition, the head at point A is sensitive to the changes in K z and z 0 due to 3-D flow (i.e., κ z d 2 <1) as discussed in Sect.3.1.In contrast, the head at point B is insensitive to the changes in K z and z 0 because the vertical flow diminishes (i.e., κ z d 2 >1).Moreover, the head at point A is sensitive to the change in L 1 but the head at point B is not because its location is far away from the well.Furthermore, the normalized sensitivity coefficient of K 1 for point A away from stream 1 approaches zero but for point B in the vicinity of stream 1 increases with time and finally maintains a certain value at the steady state.Regarding the sensitivity analysis of SDR, Huang et al. (2014) has performed the sensitivity analysis of normalized coefficients of SDR 1 to the changes in K y , K 1 , and S s for a confined aquifer and in K y , K z , K 1 , S s , and S y for an unconfined aquifer.

Effects of vertical flow and well depth on stream depletion rate
Huang et al. ( 2014) revealed that the effect of the vertical flow on SDR induced by a vertical well is dominated by the magnitude of the key factor κ z (i.e., K z y 2 0 /(K y H 2 )), where y 0 herein is a distance between stream 1 and the vertical well.They concluded that the effect is negligible when κ z ≥ 10 for a leaky aquifer.The factor should be replaced by κ z a 2 (i.e., K z a 2 /(K y H 2 )) where a is a shortest distance measured from stream 1 to the end of a lateral of a RCW, and a = a/y 0 = 1 in this study due to N = 2, θ 1 = 0, and θ 2 = π .We investigate SDR in response to various z 0 and κ z a for unconfined and confined aquifers.The temporal SDR 1 distributions predicted by Eq. ( 52) for stream 1 adjacent to an unconfined aquifer are shown in Fig. 5a for z 0 = 0.5 and κ z a 2 = 0.01, 0.1, 1, 10, 20, and 30 and Fig. 5b for κ z a 2 = 1 and 30 when z 0 = 0.1, 0.3, 0.5, 0.7, and 0.9.The curves of SDR 1 versus t are plotted in both panels by the present SDR solution for a confined aquifer.In Fig. 5a, the present solution for an unconfined aquifer predicts a close SDR 1 to that for the confined aquifer when κ z a 2 = 0.01, indicating that the vertical flow in the unconfined aquifer is ignorable.The SDR 1 for the unconfined aquifer with κ z a 2 = 30 behaves like that for a confined one, indicating the vertical flow can also be ignored.The SDR 1 is therefore independent of well depths z 0 when κ z a 2 = 30 as shown in Fig. 5b.We may therefore conclude that, under the condition of κ z a 2 ≤ 0.01 or κ z a 2 ≥ 30, a 2-D horizontal flow model can give good predictions in SDR 1 for unconfined aquifers.In contrast, SDR 1 increases with decreasing κ z a 2 when 0.01<κ z a 2 <30 in Fig. 5a, indicating that the vertical flow component induced by pumping in unconfined aquifers significantly affects SDR 1 .The effect of well depth z 0 on SDR 1 is also significant as shown in Fig. 5b when κ z a 2 = 1.Obviously, the vertical flow effect should be considered in a model when 0.01<κ z a 2 <30 for unconfined aquifers.
It is interesting to note that the SDR 1 or SDR 2 induced by two laterals (i.e., θ 1 = 0 and θ 2 = π ) parallel to the streams adjacent to a confined aquifer is independent of κ z a 2 and z 0 but depends on the aquifer width of w y .The temporal SDR distribution curves based on Eqs. ( 52) and ( 53) with γ = 0 for a confined aquifer with w y = 2, 4, 6, 10, and 20 are plotted in Fig. 6.The dimensionless distance between the well and stream 1 is set to unity (i.e., y 0 = 1) for each case.The SDR 1 predicted by the solution by Hunt (1999) based on a vertical well in a confined aquifer extending infinitely is considered.The present solution for each w y gives the same SDR 1 as the Hunt solution before the time when stream 2 contributes filtration water to the aquifer and influences the supply of SDR 1 .It is interesting to note that the sum of steady-state SDR 1 and SDR 2 is always unity for a fixed w y .The former and latter can be estimated by (w y − 1)/w y and 1/w y , respectively.Such a result corresponds with that in Sun and Zhan (2007), which investigates the distribution of steady-state SDR 1 and SDR 2 induced by a vertical well.

Concluding remarks
This study develops a new analytical model describing 3-D flow induced by a RCW in a rectangular confined or unconfined aquifer bounded by two parallel streams and noflow stratums in the other two sides.The flow equation in terms of the hydraulic head with a point sink term is employed.Both streams fully penetrate the aquifer and are un-  52) and ( 53) with γ = 0 for confined aquifers when w y = 2, 4, 6, 10, and 20.
der the Robin condition in the presence of low-permeability streambeds.A first-order free surface Eq. ( 8) describing the water table decline gives good predictions when the conditions |h| /H ≤ 0.1 and |∂h/∂x| + |∂h/∂y| ≤ 0.01 are satisfied.The flux across the well screen might be uniform on a lateral within 150 m.The head solution for the point sink is expressed in terms of a triple series derived by the methods of Laplace transform and finite integral transform.The head solution for a RCW is then obtained by integrating the point-sink solution along the laterals and dividing the integration result by the sum of lateral lengths.The integration can be done analytically due to the aquifer of finite extent with Eqs.(3)-( 6).On the basis of Darcy's law and the head solution, the SDR solution for two streams can also be acquired.The double integrals of defining the SDR in Eqs. ( 50) and (51) can also be done analytically due to considerations of Eqs. ( 3)-( 6).The sensitivity analysis is performed to explore the response of the head to the change in each of the hydraulic parameters and variables.New findings regarding the responses of flow and SDR to pumping at a RCW are summarized below.
Groundwater flow in a region based on d< √ 1/κ z is 3-D, and temporal head distributions exhibit the unconfined behavior.A mathematical model should consider 3-D flow when predicting the hydraulic head in the region.Beyond this region, groundwater flow is horizontal, and temporal head distributions display the confined behavior.A 2-D flow model can predict accurate hydraulic head.
The aquifer anisotropy of κ x >10 causes unidirectional flow in the strip region determined based on κ x (y − 1) 2 >2 for a horizontal well.Existing models assuming 2-D flow in a vertical plane with neglecting the flow component along the well give accurate head predictions in the region.
The aquifer anisotropy of κ x <1 produces significant change in the head (i.e., h > 10 −5 or |h| > 1 mm) in the strip area determined by (x − maxx k ) 2 ≤ 52.9κ x and (x − minx k ) 2 ≤ 52.9κ x for a RCW with irregular lateral configurations.
The hydraulic head in the whole domain is most sensitive to the change in K y , second-most sensitive to the change in K x , and third-most sensitive to the change in S y .They are thus the most crucial factors in designing a pumping system.
The hydraulic head is sensitive to changes in K z , S s , z 0 , and L k in the region of d< √ 1/κ z and is insensitive to the changes of them beyond the region.
The hydraulic head at observation points near stream 1 is sensitive to the change in K 1 but away from the stream is not.
The effect of the vertical flow on SDR is ignorable when κ z a 2 ≤ 0.01 or κ z a 2 ≥ 30 for unconfined aquifers.In contrast, neglecting the effect will underestimate SDR when 0.01<κ z a 2 <30.
For unconfined aquifers, SDR increases with dimensionless well depth z 0 when 0.01<κ z <30 and is independent of z 0 when κ z ≤ 0.01 or κ z ≥ 30.For confined aquifers, SDR is independent of z 0 and κ z .For both kinds of aquifers, the distribution curve of SDR versus t is independent of aquifer anisotropy κ x .
and Zlot-C.-S.Huang et al.: Approximate analysis of three-dimensional groundwater flow nik (

Figure 1 .
Figure 1.Schematic diagram of a radial collector well in a rectangular unconfined aquifer.
C.-S. Huang et al.: Approximate analysis of three-dimensional groundwater flow fined by for β n , λ 0 and λ i y = w y .

Table 1 .
Symbols used in the text and their definitions.