On the coupled unsaturated–saturated flow process induced by vertical, horizontal, and slant wells in unconfined aquifers

Conventional models of pumping tests in unconfined aquifers often neglect the unsaturated flow process. This study concerns the coupled unsaturated–saturated flow process induced by vertical, horizontal, and slant wells positioned in an unconfined aquifer. A mathematical model is established with special consideration of the coupled unsaturated–saturated flow process and the well orientation. Groundwater flow in the saturated zone is described by a three-dimensional governing equation and a linearized threedimensional Richards’ equation in the unsaturated zone. A solution in the Laplace domain is derived by the Laplace– finite-Fourier-transform and the method of separation of variables, and the semi-analytical solutions are obtained using a numerical inverse Laplace method. The solution is verified by a finite-element numerical model. It is found that the effects of the unsaturated zone on the drawdown of a pumping test exist at any angle of inclination of the pumping well, and this impact is more significant in the case of a horizontal well. The effects of the unsaturated zone on the drawdown are independent of the length of the horizontal well screen. The vertical well leads to the largest water volume drained from the unsaturated zone (W ) during the early pumping time, and the effects of the well orientation on W values become insignificant at the later time. The screen length of the horizontal well does not affect W for the whole pumping period. The proposed solutions are useful for the parameter identification of pumping tests with a general well orientation (vertical, horizontal, and slant) in unconfined aquifers affected from above by the unsaturated flow process.


Introduction
In addition to conventional vertical wells, horizontal and slant pumping wells have been broadly used in the petroleum industry and in environmental and hydrological applications in recent decades.Horizontal and slant pumping wells are commonly installed in shallow aquifers to yield a large amount of groundwater (Bear, 1979) or to remove a large amount of contaminant (Sawyer and Lieuallen-Dulam, 1998).Horizontal and slant wells have some advantages over vertical wells (Yeh and Chang, 2013;Zhan and Zlotnik, 2002); e.g., horizontal and slant wells yield smaller drawdowns than the vertical wells with the same pumping rate per screen length.Horizontal and slant wells have long screen sections, which can extract a great volume of water in shallow or low permeability aquifers without generating significant drawdowns.Hantush and Papadopulos (1962) first investigated the problem of fluid flow to a horizontal well in the hydrologic sciences.Since then, this problem has not been of great concern in the hydrological science community because of the limitations in directional drilling techniques and high drilling costs.With significant advances in directional drilling technology over the last 20 years, the interest in horizontal and/or slant wells has been reignited.Until now, flow to horizontal and/or slant wells has been investigated in various aspects, including flow in confined aquifers (Cleveland, 1994;Zhan, 1999;Zhan et al., 2001;Kompani-Zare et al., 2005), unconfined aquifers (Huang et al., 2011(Huang et al., , 2016;;Rushton and Brassington, 2013;Zhan and Zlotnik, 2002;Mohamed and X. Liang et al.: On the coupled unsaturated-saturated flow process induced by wells in unconfined aquifers Figure 1.The schematic diagram of groundwater flow to a horizontal well (a) and a slant well (b) in an unsaturated-saturated system.Rushton, 2006;Kawecki and Al-Subaikhy, 2005), leaky confined aquifers (Zhan and Park, 2003;Sun and Zhan, 2006;Hunt, 2005), and fractured aquifers (Nie et al., 2012;Park and Zhan, 2003;Zhao et al., 2016).The readers can consult Yeh and Chang (2013) for a recent review of well hydraulics on various well types, including horizontal and slant wells.
As demonstrated in previous studies, horizontal and slant wells have significant advantages over vertical wells in unconfined aquifers, and thus they were largely used in unconfined aquifers for pumping or drainage purposes.However, none of above-mentioned studies considered the effects of unsaturated processes on groundwater flow to horizontal and slant wells in unconfined aquifers.In the case of flow to vertical wells in saturated zones, the effects of unsaturated processes have been investigated by several researchers (Kroszynski and Dagan, 1975;Mathias and Butler, 2006;Tartakovsky and Neuman, 2007;Mishra andNeuman, 2010, 2011).For example, Tartakovsky and Neuman (2007) considered axisymmetric unsaturated-saturated flow for a pumping test in an unconfined aquifer and employed one parameter that characterized both the water content and the hydraulic conductivity as functions of the pressure head, assuming an infinite thickness of the unsaturated zone.Mishra andNeuman (2010, 2011) extended the solution of Tartakovsky and Neuman (2007) using four parameters to represent the unsaturated zone properties and considered a finite thickness of the unsaturated zone (Mishra and Neuman, 2010) as well as the wellbore storage (Mishra and Neuman, 2011).The main results from the studies concerning vertical wells indicated that the unsaturated zone often had a major impact on the S-shaped drawdown curves.
The following question remains: are these conclusions for vertical wells also applicable to horizontal and slant wells when coupled unsaturated-saturated flow is of concern?Specifically, how important is the wellbore orientation on groundwater flow to a horizontal or slant well considering the coupled unsaturated-saturated flow process?In order to answer these questions, we establish a mathematical model for groundwater flow to a general well orientation (vertical, horizontal, and slant wells) considering the coupled unsaturated-saturated flow process.We incorporate a threedimensional linearized Richards' equation into a governing equation for groundwater flow in an unconfined aquifer.We employ the Laplace-finite-Fourier-transform and the method of separation of variables to solve the governing equations for coupled unsaturated-saturated flow.This paper is organized as follows: we first present the mathematical model and the solution in Sects. 2 and 3, respectively, then describe the results and discussion in Sect. 4. We summarize this study and draw conclusions in Sect. 5.

Mathematical model
The schematic diagrams of flow to horizontal and slant wells in an unsaturated-saturated system are represented in Fig. 1a and b, respectively.Similar to the conceptual model used by Zhan and Zlotnik (2002), the origin of the Cartesian coordinate is located at the bottom of the saturated zone with the z axis along the upward vertical direction and the x and y axes along the principal horizontal hydraulic conductivity directions.The horizontal and slant well screens are located in the saturated zone with a distance z w from the center point of the screen (0, 0, z w ) to the bottom of the saturated zone.The slant well has three inclined angles γ x , γ y , and γ z with the x, y, and z axes, respectively, with the three angles satisfying cos 2 (γ x ) + cos 2 (γ y ) + cos 2 (γ z ) = 1.The horizontal well is a specific case of the slant well when γ z = π/2.The saturated zone is assumed as an infinite lateral extent of unconfined aquifer with a slight compressibility; it is spatially uniform and anisotropic (Tartakovsky and Neuman, 2007).The saturated zone is below an initially horizontal water table at z = d, and the unsaturated zone is above z = d with an initial thickness b.
In order to solve the problem of groundwater flow to a horizontal or slant well, we first solve the governing equation of groundwater flow to a point sink.The mathematical model for groundwater flow to a point sink (x 0 , y 0 , z 0 ) in a homogeneously anisotropic saturated zone is given by where s is the drawdown (the change in hydraulic head from the initial level) in the saturated zone [L].K x , K y , and K z are the saturated principal hydraulic conductivities in the x, y, and z directions, respectively [L T −1 ]; Q is the pumping rate (positive for pumping and negative for injecting) [L 3 T −1 ]; δ( q ) is the Dirac delta function [L −1 ]; S S is the specific storage [L −1 ]; d is the saturated zone thickness [L]; and t is the time since the start of pumping [T].It is noteworthy that the aquifer is assumed to be homogenous and spatially uniform in this study.Despite the fact that a real-world aquifer is likely to be heterogeneous and/or non-uniform, there is evidence that a moderately heterogeneous aquifer may sometimes behave as an averaged "homogeneous" system for pumping-induced groundwater flow problems.This interesting phenomena may be due to the diffusive nature of groundwater flow, which can somewhat smooth out the effect of the heterogeneity to a certain degree (Pechstein et al., 2016;Zech and Attinger, 2016).Flow in the unsaturated zone induced by pumping in the unconfined aquifer is governed by the Richards' equation.Due to the nonlinear nature of the Richards' equation, it is difficult to analytically solve this equation except in some specific cases.Kroszynski and Dagan (1975) proposed a first-order linearized unsaturated flow equation by expanding the dependent variable in the Richards' equation as a powerfunction series when the pumping rate was less than Kd 2 , where K is the saturated hydraulic conductivity of a homogeneous medium.The readers can find the details of the linearized equation derivation in previous studies (Kroszynski and Dagan, 1975;Tartakovsky and Neuman, 2007).With such a linearized treatment, it becomes possible to analytically solve the equation of flow in the unsaturated zone.The linearized three-dimensional unsaturated flow equation is adopted in this study as follows: where u is the drawdown in the unsaturated zone [L].The functions k 0 (z) and C 0 (z) are the zero-order approximation of the relative hydraulic conductivity [dimensionless] and the soil moisture capacity [L −1 ] at the initial water content of θ 0 , respectively; k is the relative hydraulic conductivity 0 ≤ k ≤ 1; C(≥ 0) is the specific moisture capacity [L −1 ] and where κ > 0 is the constitutive exponent [L −1 ] and S y is the specific yield [dimensionless].As mentioned in the introduction, this two-parameter model was extended to a fourparameter model by Mishra andNeuman (2010, 2011).The four-parameter model may be closer to a realistic situation.However, a model with more parameters has disadvantages as well.First, it is more difficult to determine the values of those parameters precisely from a practical standpoint.Second, the predictive capability of a model with more parameters may not be better than that of a model with fewer parameters.For a discussion of this issue, one may consult Voss (2011a, b) and/or Bredehoeft (2005).In this study, we focus on the question of how important the wellbore orientation is for groundwater flow to a horizontal or slant well considering the coupled unsaturated-saturated flow process.To focus on answering this question, we prefer to use a simpler model with a balance that keeps the most important physical processes in the model but, at the same time, ignores the secondary effects.
Eq. (3b) shows that, at the water table (z = d), a smaller κ leads to a smaller C 0 (z) and a larger retention capacity (Kroszynski and Dagan, 1975;Tartakovsky and Neuman, 2007); i.e., water in the unsaturated zone becomes more difficult to drain.In this study, we assume the upper boundary of the unsaturated zone as a no-flow boundary condition in Eq. ( 2c) by neglecting the effects of both infiltration and evaporation during pumping.Because typical pumping tests usually cover much shorter periods of time relative to the duration of infiltration and evaporation processes, this assumption can hold for most field conditions, particularly for land with sparse vegetation where the influence of plant transpiration is limited as well.
The saturated and unsaturated flows are coupled at their interface by the continuities of pressure and vertical flux across the water table.Following linearization, they take the form The above linearized equations of Eq. ( 4a) and (4b) assume that the variation in the water table is minor in respect to the total saturated thickness.This assumption works better for horizontal wells and slant wells than for vertical wells, provided that the same pumping rate is used.This is because horizontal wells and slant wells will generate much less drawdown over laterally broader regions, while vertical wells tend to generate laterally more concentrated and much greater drawdown near the pumping wells (Zhan and Zlotnik, 2002).

Solution for a point sink
The solution to Eq. ( 1a) is obtained by the Laplace and finitecosine Fourier transform.The Laplace domain solution of Eq. ( 1a), subject to the initial conditions of Eq. ( 1b) and the boundary conditions of Eq. ( 1c) and (1d), is given as (Zhan and Zlotnik, 2002) where The subscript "D" denotes the dimensionless terms, and the definitions of all dimensionless variables are presented in the Supplement (Sect.S1).p is the Laplace transform parameter with respect to the dimensionless time, and the overbar denotes a variable in the Laplace domain; ω n is the nth eigenvalue of the Fourier transform, and it will be determined later; K 0 is the modified second-kind Bessel function of zero order; and r D = (x D , y D ) and r 0D = (x 0D , y 0D ) are the dimensionless radial vectors of the observation point and the sink point, respectively.The solution to Eq. ( 2a) is obtained by the Laplace transform and the method of separation of variables (Sect.S2).It is given as where where The eigenvalues of the finite-cosine Fourier transform ω n can be obtained by substituting Eqs. ( 5) and ( 7) into the continuities of a normal (vertical) flux equation (Eq.S6b).Details can be found in Sect.S3.On the basis of the method illustrated above, it is straightforward to obtain the Laplace domain solution s D for the case of the unconfined aquifer with a free water table boundary and without the unsaturated zone influence (Zhan and Zlotnik, 2002; abbreviated as the ZZ solution hereafter) as well as the groundwater flow to a horizontal well in a confined aquifer (Zhan et al., 2001; abbreviated as the ZWP solution hereafter).The solutions s D for these two special cases require different ω n values.For the free water table condition, the ω n is the root of ω n tan(ω n ) = p/σ (Zhan and Zlotnik, 2002).For the confined aquifer case, ω n = nπ/α z and n = 0, 1, 2, . . .(Zhan et al., 2001).

Solution for a slant pumping well
Due to the linearity of the mathematical models in Eqs. ( 1) and ( 2), the principle of superposition can be employed to extend the basic solutions of Eqs. ( 5) and ( 7).Thus, on the basis of the principle of superposition, the drawdown induced by a line sink in the saturated zone can be obtained by integrating the solution of Eqs. ( 5) and ( 7) along the well axis, provided that the pumping strength distribution along the well screen is known.The precise determination of the pumping strength distribution along a horizontal or slant well involves complex, coupled aquifer-pipe-flow (Chen et al., 2003) in which the flow inside the wellbore (pipe flow) can experience different stages of flow schemes from laminar, to transitional turbulent, to fully developed turbulent flow.Such complex coupled well-aquifer flow is beyond the scope of this study, and one may consult the recent studies of Blumenthal and Zhan (2016) or Wang and Zhan (2016) for more details.However, often one may adopt a first-order approximation of a uniform flux distribution to treat horizontal or slant wells, particularly when the well screen lengths are not extremely long (kilometers).Such an approximation has been justified by Zhan and Zlotnik (2002).In this study, a uniform flux distribution will be utilized for horizontal or slant wells hereafter to obtain the solutions.
The drawdown in the saturated and unsaturated zones due to a slant pumping well can be written as and where s ID and u ID are the Laplace transforms of s ID and u ID , respectively.They are defined in the same way as s D and u D in Eqs. ( 5) and ( 7), respectively.L D = α x L/d is the dimensionless length of the slant well screen (L); z wD = α z z w /d is the dimensionless elevation of the center of the pumping well screen; l is a dummy variable; and α x sin γ z cos γ y ) 2 .s ID and u ID will respectively reduce drawdown in the saturated and unsaturated zones due to a horizontal well when γ z = π/2.It is noteworthy that these solutions can be straightforwardly extended to situations of locationdependent pumping rates as long as the flux rate distribution along the wellbore is known a priori.To do so, one simply modifies Eqs. ( 9) and ( 10) using a location-dependent flux function inside the integration.
The drawdown in an observation (vertical) well located in the saturated zone that is screened from z l to z u (z u > z l ) can be calculated using the average of the point drawdown in Eq. ( 9) along the observation well screen (Zhan and Zlotnik, 2002): where s oD is the Laplace transform of s oD and s oD is defined in the same way as s D in Eq. ( 5); z uD = α z z u /d, z lD = α z z l /d.
It should be noted that our solutions do not account for the wellbore effects of the pumping and observation wells.Indeed, the wellbore effects have introduced additional complexity to the solutions that are already substantially more complex than the solutions excluding the unsaturated flow process.To avoid the influence of wellbore storage effects, we make the following proposal that could be implemented in future investigations of the coupled saturated-unsaturated flow process: using pack systems to insulate the screens of pumping and the observation wells; thus, wellbore storages will not be a concern.

Total volume drained from the unsaturated zone for a slant well
The dimensionless total volume drained from the unsaturated zone to the saturated zone (water flux across the water table) can be obtained by where W D is the Laplace transform of W D , W D = W 4πα 3 z Q , and W is the total volume drained from the unsaturated zone; φ = L D α z cos(γ z )/(2α x ).
It is difficult to obtain closed-form solutions by analytically inverting the Laplace transforms of Eqs. ( 5), ( 7), ( 9), (10), and ( 12), and thus a numerical inverse Laplace method is employed in this study.There are several numerical inverse Laplace methods, such as the Stehfest method (Stehfest, 1970), the Zakian method (Zakian, 1969), the Fourier series method (Dubner and Abate, 1968), the Talbot algorithm (Talbot, 1979), the Crump technique (Crump, 1976), and the de Hoog algorithm (de Hoog et al., 1982).Each method is best fitted for a particular type of problem (Hassanzadeh and Pooladi-Darvish, 2007).Chen (1985), Zhan et al. (2009a, b), and Wang and Zhan (2013) have successfully employed the Stehfest algorithm to obtain the solution in the time domain for similar problems to the one in this study.For references on different inverse Laplace methods, one can consult the reviews of Kuhlman (2013) and Wang and Zhan (2015).In this study, we use the Stehfest method to invert the Laplace solutions into the solutions in the time domain.In order to ensure the accuracy of the Stehfest method, several numerical exercises have been performed against the benchmark solutions for several special cases of the investigated problem.

Effect of unsaturated zone parameters
The main difference between the ZZ solution and the present solution is the upper boundary condition of the saturated zone.The ZZ solution considered the linearized free surface (kinematic) equation as the water table boundary that employed one parameter, i.e., a specific yield (S y ) to account for the gravity drainage after water table decline.The present solution represents a coupled water flow through both the unsaturated and saturated zones.The water table boundary is replaced by coupled interface conditions between the unsaturated and the saturated zones.Thus, the behavior of the drawdown in the saturated zone induced by the pumping wells will be affected by the unsaturated zone.To investigate the manner in which the dimensionless constitutive exponent κ D and the dimensionless unsaturated thickness b D impact the drawdown in the saturated zone induced by a horizontal pumping well, we plot the log-log graph of s ID vs. t D /r 2 D (the type curves) for different κ D and b D in Fig. 2a and b, respectively.We also compare our solution to the ZZ solution (unconfined aquifer) and the ZWP solution (confined aquifer).For convenience, we assume the horizontal well screen to be situated along the x direction, i.e., γ x = 0 and γ y = γ z = π/2.The other parameter values in Eq. ( 9) are σ = 1 × 10 −3 , L D = 1, γ = 0, α z = 1, x D = 0.5, y D = 0.05, z D = 0.8, and z wD = 0.5.The unsaturated flow has a significant impact on drawdown curves in the saturated zone when κ D is less than 10 (the unsaturated-saturated system has a large retention capacity, a small initial saturated thickness, and/or a relatively small vertical hydraulic conductivity).The impact of unsaturated flow decreases as κ D increases, becoming small or insignificant when κ D comes close to 1 × 10 3 .Our curve is almost the same as the curve of the ZZ solution when κ D = 1 × 10 3 (gray solid curve) and gradually deviates from the ZZ solution but approaches the ZWP solution as κ D decreases to 1 × 10 −5 (black solid curve).For a fixed initial saturated thickness when κ D is smaller, i.e., the unsaturated zone has a larger retention capacity and/or both the unsaturated and saturated zones have relatively smaller vertical hydraulic conductivity, water drainage from the unsaturated zone is impeded.This forces more water to be released from the compressible storage of the saturated zone, leading to a larger drawdown in the saturated zone.The opposite is true when κ D is larger.This is consistent with the findings in the vertical pumping well case (Tartakovsky and Neuman, 2007).In order to further investigate the effects of the unsaturated zone, Fig. 2c displays the drawdown curves in the unsaturated zone (u ID ) for different values of κ D (1 × 10 −5 , 1 × 10 −3 , 1 × 10 −1 , 1 × 10 1 , and 1 × 10 3 ) at z D = 1.5 with the other parameters the same as in Fig. 2a.As κ D increases, the retention capacity of the unsaturated zone decreases, and thus more water is released from the unsaturated storage.This leads to smaller drawdown in both the unsaturated and saturated zones.Figure 2d depicts the drawdown curves in the unsaturated zone for different values of b D (0.5, 1, 2, 10, and 100).As expected, the drawdown in the unsaturated zone decreases as b D increases because more water is stored in the unsaturated zone for larger b D .These results are consistent with the findings of Mishra andNeuman (2010, 2011).

Effect of well orientation and well screen length
In this section, we first investigate the effect of the inclined angle of the slant well on the type curves.Figure 3 shows the comparison between the ZZ solution and our solution with κ D = 10 for three different angles of a slant well (γ z = π/4 and π/2) at two observation points (z D = 0.9 for Fig. 3a and z D = 0.1 for Fig. 3b) with the other parameters the same as in Fig. 2. Clearly, the smaller angle creates a larger drawdown at both observation points.For the horizontal well (γ z = π/2), the discrepancy between the ZZ solution and our solution is larger than that for the vertical well (γ z = 0) at the upper observation point (Fig. 3a).Such a discrepancy diminishes at the lower observation point (Fig. 3b).This reveals that the effects of the unsaturated zone on the drawdown exist at any angle of inclination of a slant well for the upper part of the aquifer, and this impact is more significant in the case of the horizontal well.The impact of the unsaturated zone decreases when the observation point moves downward, becoming further away from the unsaturated zone, as expected.
Here, we investigate the effect of the horizontal well screen length on the drawdown.Figure 4 illustrates the comparison between the ZZ solution and our solution for three different lengths of the well screen (L D = 0.1, 1, and 10) at two observation points with the other parameters the same as in Fig. 3.This indicates that the longer well screen leads to a smaller drawdown at both the upper and lower observation points.The discrepancy between the ZZ solution and our solution is identical for different well screen lengths.This reveals that the effects of the unsaturated zone on the drawdown are insensitive to the length of the horizontal well screen.
In order to clearly illustrate the drawdown pattern in the unsaturated-saturated system, the drawdown profiles for three different angles of a slant well (γ z =π/4 and π/2) at different dimensionless times (t D = 1 × 10 3 , 1 × 10 4 , and 1 × 10 5 ) are presented in Fig. 5 in vertical cross sections.The other parameter values in Eqs. ( 9) and ( 10) 05, z wD = 0.75, γ x = 0, and γ y = π/2.As time increases, the effect of pumping gradually propagates into the unsaturated zone (z D > 1).The vertical well leads to larger drawdown in the unsaturated zone than the slant and horizontal wells.The reason is that the vertical well screen is closer to the unsaturated zone.
The water flux across the water table (Eq.12) is the volume drained from the unsaturated zone to the saturated zone.It is somewhat related to the concept of specific yield when the coupled unsaturated-saturated zone flow process is simplified into a saturated zone flow process with the water table as a free upper boundary.Thus, Eq. ( 12) reflects the impact of the unsaturated zone on the water flow in the saturated zone.Figure 6 shows the changes in the dimensionless water flux across water table, W D , with t D from the ZZ solution and our solution at three angles of a slant well screen (γ z = π/4 and π/2) (Fig. 6a) and at three screen lengths of a horizontal well (L D = 0.1, 1.0, and 10) (Fig. 6b) with the other parameters the same as in Fig. 3.
For early times of pumping, W D increases with time, and at later times W D approaches an asymptotic value that is dependent on the unsaturated parameter κ D .W D decreases as κ D decreases.The small κ D reflects the large retention capacity of the unsaturated zone, and thus it impedes water from draining from the unsaturated zone during pumping.This results in more water released from the saturated zone storage and a larger drawdown in the saturated zone (Fig. 2a).The ZZ solution overestimates W D because it neglects the effects of unsaturated flow (Fig. 6a).The W D ∼ t D curves deviate from each other considerably at different angles of a slant well, particularly at the early time.One can see from Fig. 6a that the W D of the vertical well (γ z = 0) is the largest at early times, and the W D ∼ t D curves of the three angles eventually approach the same asymptotic value at later times.This means that the vertical well leads to the greatest water drainage from the unsaturated zone at early times, and the effects of the well orientation are insignificant as time increases.In contrast to the angle of a slant well, the screen length of a horizontal well appears to have almost no impact on W D for the whole pumping period (Fig. 6b).Similar to Fig. 6a, the magnitude of W D in Fig. 6b is only dependent on the unsaturated parameter κ D .

Synthetic pumping test
In order to further verify our solution and to explore the capability of our solution to interpret pumping test results in the unsaturated-saturated system, we have conducted a synthetic numerical simulation.The synthetic case considers a pumping test in an unconfined aquifer with a slant pumping well (γ z = π/4, γ x = 0, and γ y = π/2).The aquifer parameter values are as follows: the unconfined aquifer thickness d is 10 m; the unsaturated zone thickness b is 5 m; the horizontal conductivity is K x = K y = 0.06 m min −1 ; the vertical conductivity is K z = 0.5K x ; the specific storage is S S = 1 × 10 −4 m −1 ; and the specific yield is S y = 0.3.The unsaturated flow is described by Eqs. ( 2) and (3) with the constitutive exponent κ = 0.1 m −1 .The discharge rate of the pumping well is Q = 1 m 3 min −1 , the length of the pumping well screen L is 5 m, and the center of the well screen is located at (x = 0, y = 0, and z = 5 m).
The coupled Eqs. ( 1)-( 4) of the unsaturated-saturated system are numerically solved by COMSOL Multiphysics (COMSOL Multiphysics GmbH, Göttingen, Germany), a robust Galerkin finite-element software package that includes a partial differential equation (PDE) solver for modeling the types of governing equations in this study.Figure 7a shows the spatial discretization of our COMSOL model, in which tetrahedrons are used as elements for the three-dimensional model and the elements near both the pumping well and the unsaturated-saturated interface are refined.The number of tetrahedral elements is 328 358.The time step increases exponentially, and the total number of time steps is 100 with a total simulation time of 220 min.Figure 7b presents an example of the vertical profiles (the xz plane) of the drawdown in the unsaturated-saturated system at t = 210 min.Figure 7b indicates that the COMSOL model provides a good reproduction of the drawdown in the unsaturated-saturated system induced by a slant pumping well.
First, we verify our solutions by comparing the drawdowns in both the saturated and unsaturated zones with the numerical solution for the same aquifer parameter values.Figure 8a and b shows the drawdown curves in the saturated zone at the observation point (x = 0, y = 1 m, and z = 9 m) and the drawdown curves in the unsaturated zone at the observation point (x = 0, y = 1 m, and z = 11 m) using the numerical solution (triangles) and our solution (solid curves).These fig- ures indicate that, in general, our solution satisfactorily fits the numerical solution in both the saturated and unsaturated zones, although the agreement becomes less satisfactory (but acceptable) at later times.The sizes of the tetrahedral elements will affect the accuracy of the numerical solution, especially near the pumping well and the unsaturated-saturated interface.Although we refine the mesh at these places, the sizes of these elements may be insufficiently small to completely remove the numerical errors.Our numerical exercises show that a finer element discretization for this model leads to substantially greater computational cost, probably due to the three-dimensional nature of the model.
Second, we investigate the errors in using the ZWP and ZZ solutions to explain the drawdown curves in the unsaturated-saturated system induced by the slant pumping well.Figure 8a shows a least squares fit of the ZWP (dashed curves) and ZZ (dotted curves) solutions to the numerical solution, yielding parameter estimates K x = K y = 0.13 m min −1 and S S = 1.1 × 10 −2 m −1 (for the ZWP solution) and K x = K y = 0.03 m min −1 , S S = 2.3 × 10 −4 m −1 , and S y = 0.32 (for the ZZ solution).Clearly, the ZWP solution fails to fit the numerical solution entirely and significantly overestimates the horizontal hydraulic conductivity and the specific storage by 1 or 2 orders of magnitude because it is a confined aquifer solution.The ZZ solution dramatically deviates from the numerical solution at the early and intermediate times, and it agrees with the numerical solution at later times.The ZZ solution underestimates the horizontal hydraulic conductivity and overestimates the specific storage and the specific yield.
A major disadvantage of the two older models (the ZWP and ZZ models) is that they do not consider the unsaturated flow process, and thus they cannot be used to characterize the parameters of the unsaturated zone.The newer model developed in this study, however, is capable of characterizing the parameters of both the saturated and unsaturated zones.As far as we know, this represents a significant improvement over the older models.Furthermore, as the older models do not consider the unsaturated flow process proven to be important for producing the drawdown-time curves in the saturated zone, they often cannot satisfactorily reproduce the observed drawdown-time curves in the saturated zone in actual real-world aquifer pumping tests.The newer model has resolved this issue successfully because the conceptual model is closer to the physical reality of flow in an unsaturatedsaturated system.

Summary and conclusions
The coupled unsaturated-saturated flow process induced by vertical, horizontal, and slant pumping wells is investigated in this study.A mathematical model for such a coupled unsaturated-saturated flow process is presented.The flow in the saturated zone is described by a three-dimensional governing equation, and the flow in the unsaturated zone is described by a three-dimensional Richards' equation.The unsaturated zone properties are represented by the Gardner (1958) exponential relationships.The Laplace domain solutions are derived using the Laplace transform and the method of separation of variables, and the time domain solutions are obtained using the Stehfest method (Stehfest, 1970).The solution is compared with the solutions proposed by Zhan et al. (2001) (confined aquifer; the ZWP solution) and Zhan and Zlotnik (2002) (unconfined aquifer; the ZZ solution) and verified using the finite-element numerical solution.The conclusions of this study can be summarized as follows: 1.The unsaturated flow has a significant impact on drawdown in unconfined aquifers induced by the horizontal pumping well when the dimensionless constitutive exponent κ D is less than 10 (the large retention capacity of the unsaturated zone, the small initial saturated thickness, and/or the small vertical hydraulic conductivity).
2. For the small dimensionless unsaturated thickness b D (= 0.001), the drawdown curves approach the ZWP solution.For the large unsaturated thickness b D (= 100), the drawdown curves do not approach the ZZ solution because the impact of the unsaturated flow becomes significant at a fixed κ D of 0.1.
3. The effects of the unsaturated zone on the drawdown exist at any angle of inclination of a slant well, and this impact is more significant in the case of the horizontal well.The effects of the unsaturated zone on the drawdown are insensitive to the length of the horizontal well screen.
4. For the early time of pumping, the water volume drained from the unsaturated zone (W ) to the saturated zone increases with time.With progresing time, W approaches an asymptotic value that is dependent on the unsaturated parameter κ D .The vertical well leads to the largest W value during the early time of pumping, and the effects of the well orientation become insignificant at the later time.The screen length of the horizontal well does not affect W for the whole pumping period.
5. By comparing it with the synthetic pumping test data generated by the finite-element numerical model of COMSOL, one can see that our solution provides a good reproduction of the drawdown curves in both the saturated and unsaturated zones, while both the ZWP and ZZ solutions fail to fit the drawdown curves.They either underestimate or overestimate the horizontal hydraulic conductivity, the specific storage, and the specific yield.

Data availability
The computer program codes and data used in this study can be accessed by contacting the corresponding author directly.
The Supplement related to this article is available online at doi:10.5194/hess-21-1251-2017-supplement.

Figure 2 .
Figure 2. (a) The log-log plot of s ID against t D /r 2 D for different values of the dimensionless unsaturated parameter κ D , the ZWP solution (confined aquifer), and the ZZ solution (unconfined aquifer).(b) The log-log plot of s ID against t D /r 2 D for different values of the dimensionless unsaturated thickness b D , the ZWP solution (confined aquifer), and the ZZ solution (unconfined aquifer).(c) The log-log plot of u ID against t D /r 2 D for different values of the dimensionless unsaturated parameter κ D and (d) the log-log plot of u ID against t D /r 2 D for different values of the dimensionless unsaturated thickness b D .

Figure
Figure 2a presents the drawdown curves in the saturated zone for different values of κ D (1 × 10 −5 , 1 × 10 −3 , 1 × 10 −1 , 1 × 10 1 , and 1 × 10 3 ) with a fixed dimensionless thickness of the unsaturated zone b D of 0.5.The dimensionless constitutive exponent is κ D = κd/α z = κdK 1/3 D , where K D is the anisotropic ratio between the vertical hydraulic conductivity and the horizontal hydraulic conductivity.The unsaturated flow has a significant impact on drawdown curves in the saturated zone when κ D is less than 10 (the unsaturated-saturated system has a large retention capacity, a small initial saturated thickness, and/or a relatively small vertical hydraulic conductivity).The impact of unsaturated flow decreases as κ D increases, becoming small or insignificant when κ D comes close to 1 × 10 3 .Our curve is almost the same as the curve of the ZZ solution when κ D = 1 × 10 3 (gray solid curve) and gradually deviates from the ZZ solution but approaches the ZWP solution as κ D decreases to 1 × 10 −5 (black solid curve).For a fixed initial saturated thickness when κ D is smaller, i.e., the unsaturated zone has a larger retention capacity and/or both the unsaturated and saturated zones have relatively smaller vertical hydraulic conductivity, water drainage from the unsaturated zone is impeded.This forces more water to be released from the compressible storage of the saturated zone, leading to a larger drawdown in the saturated zone.The opposite is true when κ D is larger.This is consistent with the findings in Fig.2aalso shows that the drawdowns have typical "S" pattern curves when κ D ≥ 0.1.At early times, all curves are approximately identical due to the response of the confined storage and the minor effects of the upper boundary of the saturated zone; at intermediate times, the drawdowns of the ZZ solution and our solution increase slower than those of the ZWP solution due to the response of the additional storage of the upper boundary of the saturated zone.At later times, the increasing drawdown rates of the ZZ solution and our solution are nearly the same as those of the ZWP solution due to the combined effects of both storage mechanisms.The unsaturated zone controls the effects of the additional storage and the upper boundary of the saturated zone on drawdown curves.There are physical differences between the ZZ solution and our solution.The ZZ solution uses the storage factor S y (specific yield) at the upper boundary of the saturated zone.Such a storage factor at the upper boundary is greater than the actual storage capacity of the unsaturated zone when the unsaturated parameter κ D ≤ 10, leading to a slower water level decline for the ZZ solution.Such an effect will become insignificant over a long pumping time.Similar to κ D , the dimensionless unsaturated thickness b D also affects the drawdown behavior of the saturated zone, as shown in Fig.2bfor different values of b D (0.001, 0.01, 1, 10, and 100) with a fixed κ D = 0.1 and the same parameters used

Figure 3 .
Figure 3.The log-log plot of s ID against t D /r 2 D for different angles of the well screen and a comparison with the ZZ solution for (a) dimensionless piezometer location (0, 0.05, and 0.9) and (b) dimensionless piezometer location (0, 0.05, and 0.1.)

Figure 4 .Figure 5 .
Figure 4.The log-log plot of s ID against t D /r 2 D for different dimensionless lengths of the horizontal well screen and a comparison with the ZZ solution for (a) dimensionless piezometer location (0, 0.05, and 0.9) and (b) dimensionless piezometer location (0, 0.05, and 0.1.)

Figure 6 .
Figure 6.The log-log plot of W D against t D for different values of the dimensionless unsaturated parameter κ D and the ZZ solution with (a) three angles of the slant well screen (γ z = 0, π/4, and π/2) and (b) three dimensionless lengths of the horizontal well screen (L D = 0.1, 1.0, and 10).

Figure 7 .
Figure 7. (a) The grid mesh of the unsaturated-saturated system used in the Galerkin finite-element COMSOL Multiphysics program and (b) the vertical profiles (xz planes) of the drawdown in the unsaturated-saturated system at t = 210 min for the synthetic case.

Figure 8 .
Figure 8.(a) A comparison of the synthetic drawdown in the saturated zone generated from a numerical solution with fitted analytical solutions using the the ZZ solution, the ZWP solution, and our solution.(b) A comparison of the synthetic drawdown in the unsaturated zone generated from a numerical solution with our solution.
proach the solution of the unconfined aquifer with the linearized free water table boundary (the ZZ solution).