RESEARCH ARTICLE
10.1029/2020WR027529
Key Points:• A general mathematical model
describing groundwater ﬂow in a
conﬁned aquifer for a single‐wellcirculation system is developed• An analytical solution for drawdown
in the Laplace domain and thesteady‐state analytical solution in
the time domain are obtained• The characteristics of groundwaterﬂow and the effects of different
parameters on drawdown areinvestigatedAn Analytical Solution of Groundwater Flow in a
Conﬁned Aquifer With a Single‐WellCirculation System
Kun Tu1,2,3 , Qiang Wu1,2 , Jirka Simunek3, Chaofan Chen4,5
Shengheng Xu1,2, and Yang Wang6, Ke Zhu1,2, Yifan Zeng1,2,

1

College of Geoscience and Surveying Engineering, China University of Mining and Technology, Beijing, China, 2National
Engineering Research Center of Coal Mine Water Hazard Controlling, Beijing, China, 3Department of Environmental
Sciences, University of California, Riverside, CA, USA, 4Helmholtz Centre for Environmental Research (UFZ), Leipzig,
Germany, 5Applied Environmental Systems Analysis, Dresden University of Technology, Dresden, Germany, 6Central
Plains Institute of Intelligent Geology, Zhengzhou, China

Abstract In this study, a general analytical model for the single‐well circulation system is developed toCorrespondence to:Q. Wu,
wuq@cumtb.edu.cn

Citation:
Tu, K., Wu, Q., Simunek, J., Chen, C.,
Zhu, K., Zeng, Y., et al. (2020). An
analytical solution of groundwater ﬂow
in a conﬁned aquifer with a single‐wellcirculation system. Water Resources
Research, 56, e2020WR027529. https://doi.org/10.1029/2020WR027529
Received 20 MAR 2020
Accepted 9 JUN 2020
Accepted article online 12 JUN 2020

analyze transient drawdown in a conﬁned aquifer. The analytical solution of transient drawdown in the
Laplace domain, which is numerically inverted into the time domain using the Stehfest method, isderived by employing a combination of the Laplace and Fourier cosine transforms. The characteristics oftransient drawdown and the effects of different parameters related to the single‐well circulation systemon drawdown are investigated. Furthermore, the analytical solution under steady‐state conditions is
obtained using the Fourier cosine transform. The results show that steady drawdown contours are
symmetric around a horizontal midplane of an aquifer and vary tremendously with distance from the wellaxis. The contours of drawdown around the sealed section are dense, meaning that the hydraulic gradientin this area is relatively large. The sensitivity analysis, performed to evaluate the characteristics ofdrawdown to changes in each parameter, indicates that the radial hydraulic conductivity and the lengthof the sealed section have a large impact on the drawdown and that each parameter has its inﬂuence
period on the drawdown.

1. Introduction
Among different types of groundwater heat pump systems, the single‐well circulation groundwater heatpump system is different from conventional groundwater heat pump systems—for example, a well doublet system (Banks, 2009; Galgaro & Cultrera, 2013) or a standing column well system (Abu‐Nadaet al., 2008; Deng et al., 2005—especially in terms of the well structure. In a single‐well circulation system(Figure 1), a single borehole is divided into two sections by well packers (the gray section in Figures 1 and2), which are placed in the middle of a single borehole to block injected water in the injection section of
the well from ﬂowing to the pumping section of the well. The single‐well circulation systems can be
regarded as a combination of two partially penetrating wells. The lower one is a pumping well, fromwhich groundwater is pumped at a rate Q, and the upper one is an injection well for injecting water at
the same rate Q.

©2020. American Geophysical Union.
All Rights Reserved.

TU ET AL.

A single‐well circulation groundwater heat pump system is a new technology that emerged in 2001 and sincethen has been employed to provide heating and cooling for buildings in China (Xu & Rybach, 2003). This
unique groundwater heat pump system has gained tremendous interest and has been increasingly used during the last two decades, especially in northern China. Although the system has many practical applications,the theoretical investigation of this system is somewhat limited. Xu and Rybach (2005) presented the development and several successful engineering implementations of this original and innovative technology interms of system designs, energy savings, and environmental protection. Wu et al. (2015) compared thesingle‐well circulation system with other conventional groundwater heat pump systems in terms of the well
conﬁguration, requirements on hydrogeological and thermal geological conditions, environmental impacts,
regulatory requirements, and so forth. Rybach (2015) presented engineering applications in China and othercountries and summarized operational principles of the single‐well circulation groundwater heat pump
1 of 14

Water Resources Research

10.1029/2020WR027529

Figure 1. An overview of a single‐well circulation system for shallow geothermal energy extraction (Wu et al., 2019).

system. Finally, Zeng et al. (2017) provided a detailed discussion on energy efﬁciency, cost‐effectiveness,
operation, and maintenance of this system using a case study in Nebraska.
The heat exchange in a single‐well circulation groundwater heat pump system is a result of combined groundwater ﬂow and heat transport. Heat convection, in particular, is inﬂuenced by ﬂuid ﬂow. Therefore, a goodunderstanding of groundwater ﬂow is an especially vital and essential step
before one can start studying the heat exchange in a single‐well circulation
groundwater heat pump system. Up to now, most of the research on
groundwater ﬂow focused on fully or partially penetrating wells.

Figure 2. The schematic diagram of a single‐well circulation system.

TU ET AL.

The topic of drawdown around a fully penetrating well has been discussed
in many works. For instance, Theis (1935) was the ﬁrst who presented ananalytical solution describing groundwater ﬂow toward a fully penetrating well in a conﬁned aquifer of an inﬁnite extent. Later, Papadopulos
and Cooper (1967) proposed a mathematical model describing groundwater ﬂow in a conﬁned aquifer while considering the inﬂuence of thewell radius and wellbore storage. Chen (1984) derived, based on several
existing solutions, an analytical solution describing groundwater ﬂow ina ﬁnite conﬁned aquifer with a zero drawdown condition at the exterior
boundary. Novakowski (1989) obtained a Laplace domain solution ofdrawdown in a composite conﬁned aquifer while considering the effectsof the well radius and wellbore storage. Subsequently, Yeh et al. (2003)
presented an analytical solution describing the head distribution for a slug
test in a two‐zone conﬁned aquifer system. Wen and Zhan (2008)
2 of 14

Water Resources Research

10.1029/2020WR027529

Table 1Comparison of Analytical Models for Partially and Fully Penetrating Wells
Well radius

Groundwater ﬂow

Skin effectZero

Darcian ﬂow

No

A conﬁned aquifer with an inﬁnite areal extent

Theis (1935)
Chen (1984)
Wen and Zhan (2008)
Lin et al. (2016)

Two partially penetrating wells
in one wellbore
A fully penetrating well
A fully penetrating well
A fully penetrating well
A fully penetrating well

Zero
Zero
Zero/Finite
FiniteDarcian ﬂow
Darcian ﬂow
non‐Darcian ﬂow
Darcian ﬂow

No
No
No
Yes

Hantush (1961)
Chiu et al. (2007)

A partially penetrating well
A partially penetrating well

Zero
FiniteDarcian ﬂow
Darcian ﬂow

NoYes

Wen et al. (2013)
Feng and Wen (2016)

A partially penetrating well
A partially penetrating well

Zero
Finitenon‐Darcian ﬂow
non‐Darcian ﬂowNoYes

A conﬁned aquifer with an inﬁnite areal extent
A ﬁnite conﬁned aquifer
A conﬁned aquifer with an inﬁnite areal extent
A two‐zone conﬁned aquifer with a ﬁnite areal
extentA conﬁned aquifer with an inﬁnite areal extent
A two‐zone conﬁned aquifer with an inﬁnite
areal extent
A conﬁned aquifer with an inﬁnite areal extent
A two‐zone conﬁned aquifer with an inﬁnite
areal extent

Model
This study

Well structure

Type of aquifer

proposed a general analytical model describing non‐Darcian groundwater ﬂow in a conﬁned aquifer by
employing a combination of the power law function and a linearization procedure. Wang et al. (2012)
presented an analytical solution describing groundwater ﬂow while considering a constant‐ﬂux pumping
along a ﬁnite radius well in the skin zone. Lin et al. (2016) proposed an analytical model describing
groundwater ﬂow in a composite conﬁned aquifer of a ﬁnite extent induced by a pumping test whileconsidering the effect of the well skin and a Robin‐type outer boundary.
There are also a lot of studies available in the literature regarding groundwater ﬂow to a partially penetrating
well (Table 1). The ﬁrst mathematical model for a partially penetrating well was presented by Hantush(1957), who derived an analytical solution by applying the Laplace and Fourier cosine transforms.
Since then, partially penetrating well systems have been investigated by many other researchers (e.g.,
Ataie‐Ashtiani et al., 2012; Chang & Chen, 2002; Chang & Yeh, 2009; Chen et al., 2010; Chiu et al., 2007;
Yang & Yeh, 2005; Zlotnik et al., 1998). For instance, steady‐state ﬂow in an unconﬁned aquifer to one or
multiple partially penetrating wells was studied by Luther and Haitjema (1999). Yang and Yeh (2005) presented a semianalytical solution describing groundwater ﬂow toward a partially penetrating well inducedby the constant head test in an inﬁnite extent aquifer. Taking into account the effect of a ﬁnite‐thickness
skin, Chiu et al. (2007) presented an analytical solution for a drawdown of groundwater ﬂow to a partiallypenetrating well. Barua and Bora (2010) presented a steady‐state solution for groundwater ﬂow inducedby pumping in a partially penetrating well in a ﬁnite conﬁned aquifer associated with the condition of theconstant head at the outer boundary. Ataie‐Ashtiani et al. (2012) derived an analytical solution for the capture zone of a partially penetrating well with an inﬁnitesimally small well radius and a ﬁxed pumping rate Q.
By employing the Izbash equation and a linearization procedure, Wen et al. (2013) presented an approximate analytical solution, which can describe the non‐Darcian ﬂow toward a partially penetrating well in a
conﬁned aquifer. Feng and Wen (2016) proposed an analytical model describing non‐Darcian ﬂow towarda partially penetrating well in a conﬁned aquifer, which considered the effects of the well skin. Table 1 summarizes and compares the available, above‐mentioned analytical models.
However, compared with the studies mentioned above of groundwater ﬂow in an aquifer with fully or partially penetrating wells, the research on a particular well structure of a single‐well circulation system is quite
limited. Only a few studies (Ni et al., 2011; Sorensen & Reffstrup, 1992; Tu et al., 2019) dealing with groundwater ﬂow in an aquifer with a single‐well circulation system can be found in the literature. Sorensen andReffstrup (1992) were the ﬁrst to present a simpliﬁed mathematical model for the analysis of groundwaterﬂow and heat transport for a new type of a groundwater heat pump system, which can extract and injectgroundwater with one well borehole. Ni et al. (2011) developed a mathematical model describing groundwater ﬂow induced by a single‐well system in a leaky conﬁned aquifer. Based on the superposition principle
and by adopting the analytical solution of a partially penetrating well derived by Hantush (1961), Ni
et al. (2011) obtained their solution by summing up the analytical solutions for these two partially penetrating wells. Although it is convenient to obtain analytical solutions using the superposition principle, thisapproach inevitably introduces problems such as complicated well conﬁguration functions. The Hantush
TU ET AL.

3 of 14

Water Resources Research

10.1029/2020WR027529

and Jacob well function (Hantush, 1961) makes the analytical solution more complicated and increases its
instability and computational cost. Subsequently, Tu et al. (2019) derived a closed‐form, steady‐state analytical solution to investigate groundwater ﬂow in a conﬁned aquifer by employing the separation variable
method and the Laplace transform. However, their solution still includes the complicated Hantush andJacob well function and the well conﬁguration function. Therefore, to overcome the above‐mentioned
problems, we aim at rigorously deriving the analytical solution based on the governing equation for asingle‐well circulation system by employing a combination of the Laplace and Fourier cosine transforms.
The main objective of this work is to develop a general analytical mathematical model describing groundwater ﬂow in a conﬁned aquifer for a single‐well circulation system with a unique well structure with twopartially penetrating wells in one single wellbore. The analytical solution for drawdown in the Laplacedomain is obtained by employing the Laplace and Fourier cosine transforms. The effects of the hydraulicconductivity, the length of the sealed section, and other parameters on the drawdown distribution were analyzed, which to our knowledge have not been investigated by other researchers. The steady‐state solution in
the time domain is also derived. The analytical solution of transient drawdown is used in the sensitivity analysis to investigate the effects of different parameters on drawdown and to identify which of these parametersmust be carefully considered when designing a single‐well circulation system.

2. Mathematical Model
In order to investigate groundwater ﬂow in a single‐well circulation system, a mathematical model wasdeveloped, as shown in Figure 2. The assumptions of this mathematical model are as follows: (1) A conﬁned
aquifer is assumed to be homogeneous, anisotropic, and inﬁnite in the horizontal direction; (2) underlying
and overlying rocks are considered to be impermeable, homogeneous, anisotropic, and of uniform thickness;(3) groundwater ﬂow in a conﬁned aquifer follows Darcy's law; (4) the storage and release of groundwater
due to rising and decreasing hydraulic heads are instantaneous; (5) pumping and injection rates Q (positive
for pumping and negative for injection) are constant; and (6) the well radius is inﬁnitesimally small and canbe regarded as 0.
2.1. Governing Equation
Base on the above‐listed assumptions, the equation describing groundwater ﬂow in a conﬁned aquifer can be
written as follows:
Kr

 2

∂ sðr; z; t Þ 1∂sðr; z; t Þ∂2 sðr; z; t Þ∂sðr; z; t Þ
þ
K
;
þ
¼S
z
∂r 2
r
∂r
∂z2
∂t

(1)

where r is the radial coordinate (L); t denotes time (T); s is groundwater drawdown (L); Kr and Kz represent the horizontal and vertical hydraulic conductivities, respectively (LT−1); and S refers to the speciﬁcstorage of an aquifer (L−1).
2.2. Initial Conditions
Initial drawdown in a conﬁned aquifer at time t = 0 is assumed to be constant and equal to 0:
sðr; z; 0Þ ¼ 0:

(2)

2.3. Boundary Conditions
The boundary conditions in inﬁnity in the radial direction and at the top and bottom of a conﬁned aquifer,r = ∞, z = d, and z = 0, respectively, are given as follows:
sð∞; z; t Þ ¼ 0(3)∂sðr; d; t Þ
¼0
∂z

(4)

and

and

TU ET AL.

4 of 14

Water Resources Research

10.1029/2020WR027529∂sðr; 0; t Þ
¼ 0;
∂z

(5)

where z is the vertical coordinate (L) and d denotes the thickness of a conﬁned aquifer (L). The boundary
conditions at the rim of the wellbore are as follows:
8 −Q
>
>
>
> 2πK r d1
∂s <0
limr→0 r ¼
∂r >
>
>
Q
>
:
2πK r d3

ð0 ≤ z ≤ d1 Þðd1 ≤ z ≤ d1 þ d2 Þ ;(6)

ðd1 þ d2 ≤ z ≤ dÞ

where d1, d2, and d3 are the lengths of pumping, sealed, and injection sections, respectively (L), and Q
represents pumping and injection rates (L3T−1), which are the same and constant in time.
2.4. A General (Transient) Analytical Solution
The mathematical model of a single‐well circulation system is fully described using the governing partial differential equation (Equation 1), initial conditions (Equation 2), and boundary conditions (Equations 3–6).
This system of equations can be solved by applying the Laplace and Fourier cosine transforms, which convert partial differential equations to ordinary differential equations. A detailed derivation of the analyticalsolution is given in Appendix A. The analytical solution of drawdown in the Laplace domain is expressed as

3


Nπd1
Nπ ðd1 þ d2 Þ


sin
sin
7 pﬃﬃﬃ 
Q 6
Nπz
d
d
∞
6
7
sðr; z; pÞ ¼ ∑N¼1
þ
5K 0 δr cos d ;NpK r π 2 4
d3
d1
2

(7)

K z N 2 π 2 Sp
þ , K0(·) denotes zeroth‐order modiﬁed Bessel functions of the second kind, and pKrK r d2and N (N = 1, 2, 3,…) are the Laplace and Fourier transform variables, respectively.
in which δ ¼

It is currently not possible to obtain the analytical solution of Equation 7 in the real‐time domain due tothe Bessel function. However, Equation 7 can be numerically inverted to the real‐time domain using several numerical inversion methods (e.g., Crump, 1976; De Hoog et al., 1982; Stehfest, 1970a, 1970b). In thisstudy, the Stehfest method is utilized to invert the Laplace domain solution numerically and to calculate
drawdown in the real‐time domain. The detailed introduction of the Stehfest method can be found in
Stehfest (1970a, 1970b).
A MATLAB program has been developed to calculate drawdown in the Laplace domain, which is theninverted numerically using the Stehfest method. This inversion method requires as input the number ofterms used in the inversion, n, and the number of terms of the inﬁnite series, N. Numerical tests indicate thatthe choice of n = 20 and N ≥ 10 gives accurate and stable results. To reduce the computational cost, n = 20
and N = 10 are used in all calculations.
2.5. A Steady‐State Analytical Solution
When the pumping and injection times are large, namely, t → ∞, the groundwater ﬂow will reach aquasi‐steady state. In this case, drawdown around partially penetrating wells stops changing with time.
Therefore, we have∂sðr; z; t Þ
¼ 0:
∂t

(8)

 2

∂ sðr; zÞ 1 ∂sðr; zÞ∂2 sðr; zÞ
þ Kz
þ
¼ 0:
2
∂r
r ∂r
∂z2

(9)

Then, Equation 1 can be reduced toKr

TU ET AL.

5 of 14

Water Resources Research

10.1029/2020WR027529

Applying the Fourier cosine transform to the second‐order partial derivative with respect to the z coordinate, the partial differential equation 9 issolved (see Appendix B) to give

3

2 
Nπd1
Nπ ðd1 þ d2 Þ


sin
sin
7 pﬃﬃﬃﬃ′ 
Q 6
Nπz
d
d
∞
7K 0
6
δ r cos
sðr; zÞ ¼ ∑N¼1
þ
;
5
4
2
d3
d1
NK r πd(10)where δ′ ¼

Figure 3. Comparison of the newly derived analytical solution fordrawdown (Equation 7) with the result of Ni et al. (2011) for parameters
given in the ﬁgure.

K z Nπ 2.
Kr d

According to Darcy's Law, steady‐state groundwater ﬂuxes can beobtained from the analytical solution for steady‐state drawdown(Equation 10) as follows:8∂s>>< νr ¼ K r
∂r
:
(11)
∂s
>
>
: νz ¼ K z∂z

dK 0 ðx Þ¼ − K 1 ðx Þ,dx
where x is the variable of the modiﬁed Bessel function. Thus, Equation 11 can be written as

According to the properties of the modiﬁed Bessel function, the mathematical relation


3

8
2 
Nπd1
Nπ ðd1 þ d2 Þ>rﬃﬃﬃﬃﬃﬃ sin
>


sin
>
>
7 pﬃﬃﬃﬃ′ 
Nπz
d
d
>
∞ −Q K z 6
>
6
7K 1
>
ν
¼
∑
þ
δ
r
cos
r
>
N¼1
5
>
πd K r 4
d
d3
d1
>
>
<
:

3

2 
>
Nπd1
Nπ ðd1 þ d2 Þ
>
>


sin
sin
>
>
7 pﬃﬃﬃﬃ′ 
Nπz
>
d
d
∞ −QK z 6
>
7K 0
>
þ
δ r sin> νz ¼ ∑N¼1 πdK 6
4
5
>
d
d3
d1
>
r
:

(12)

3. Results and Discussion
3.1. Veriﬁcation of the Analytical Solution
First, the newly derived analytical solution for transient drawdown is veriﬁed by comparing its results with the solution of Ni et al. (2011). In theirwork, a single‐well circulation system is regarded as a combination of one
partially penetrating well for injection and another one for pumping. The
analytical solution for a partially penetrating well was developed by
Hantush (1961). Based on the superposition principle, Ni et al. (2011)
obtained the analytical solution by summing up the analytical solutions
for these two partially penetrating wells in a single‐well circulation system. These analytical solutions (i.e., Equation 7 and Ni et al., 2011) for asingle‐well circulation system are used to calculate transient groundwater
drawdown for a system with the following parameters: Q = 60 m3/hr,
Kr = 0.1 m/hr, Kz = 0.01 m/hr, S = 0.0001 m−1, d = 40 m, d1 = 15 m,
d2 = 10 m, and d3 = 15 m, at r = 5 m and z = 15 m. Figure 3 shows thatthe results of the newly derived analytical solution correspond very wellwith the solution of Ni et al. (2011).
3.2. Steady‐State Drawdown
Figure 4. Contours of steady‐state drawdown and steady‐state ﬂux vectors
in a single‐well circulation system.

TU ET AL.

Figure 4 depicts contours for steady‐state drawdown and steady‐state ﬂux
vectors (Equations 10 and 12, respectively) in a single‐well circulation
6 of 14

Water Resources Research

10.1029/2020WR027529

system. The following parameters (the same as in Figure 3) were used inthis scenario: Q = 60 m3/hr, Kr = 0.1 m/hr, Kz = 0.01 m/hr,
S = 0.0001 m−1, d = 40 m, d1 = 15 m, d2 = 10 m, and d3 = 15 m. It can
be clearly observed that contours for steady‐state drawdown are symmetric around a horizontal midplane of a conﬁned aquifer and vary tremendously with distance from the well axis. The contours of drawdown
around the sealed section are dense, indicating that the hydraulic gradientin this area is relatively large. Additionally, due to ﬂow symmetry, absolute values of drawdown are very small around the midplane of a conﬁned
aquifer; drawdown is 0 at locations with z = 20 m (the middle of a conﬁned aquifer).
3.3. Radial Hydraulic Conductivity
Figure 5 shows the effects of different values of the radial hydraulic conductivity (Kr = 0.01, 0.02, 0.05, and 0.1 m/hr) on drawdown with time.
Figure 5. The effect of the radial hydraulic conductivity (Kr) on drawdown
Other parameters remain the same as in previous examples; that is,with time (for parameters given in the ﬁgure).
Q = 60 m3/hr, Kz = 0.01 m/hr, S = 0.0001 m−1, d = 40 m, d1 = 15 m,
d2 = 10 m, and d3 = 15 m, at r = 5 m and z = 15 m. It can be seen thatdrawdown increases with the radial hydraulic conductivity (Kr) at early times, while it decreases as Kr
increases at late times. Larger hydraulic conductivities will cause groundwater to react more quickly tothe start of water pumping/injection from/into the aquifer at early times, with the depression cone in thesuction zone spreading faster, while groundwater can be replenished in a timely manner. Drawdown is thus
larger at early times for larger hydraulic conductivities. In contrast, at late times, an increase in drawdown
for larger hydraulic conductivities is not as large as for smaller hydraulic conductivities. Moreover, it isworth noting that drawdown converges toward a steady state more quickly for larger radial hydraulic conductivities, as indicated in Figure 5.
3.4. Sealed Section
Figure 6 shows the impact of the length of the sealed section (d2) on groundwater drawdown in a single‐wellcirculation system. The sealed section is located in the middle of the conﬁned aquifer and between thepumping and injection sections of the well (see Figure 2). The same parameters as in the previous section
are used to obtain results displayed in Figure 6, that is, Q = 60 m3/hr, Kz = 0.01 m/hr, S = 0.0001 m−1,
d = 40 m, d1 = 15 m, and d3 = 15 m (i.e., d1 = d3), at r = 5 m and z = 15 m. Besides, the following lengthsof the sealed section are used: d2 = 6, 8, 10, 12, and 14 m, while d1 and d3 were correspondingly adjusted tokeep d constant. One can observe in Figure 6a that shorter lengths of the sealed section result in larger drawdown at early times, while the length of the sealed section has only a little effect on drawdown at late times.
This is due to the fact that shorter lengths of the sealed section require longer lengths of the pumping andinjection sections, which require smaller pumping and injection rates for a unit well length (note that the
total pumping and/or injection rates are constant). Thus, shorter lengths of the sealed section produce larger

Figure 6. The effect of the length of the sealed section (d2) when d1 = d3 (a) and various ratios κ (=d1/d3) for a constantd2 (b) on drawdown with time.

TU ET AL.

7 of 14

Water Resources Research

10.1029/2020WR027529

drawdowns at early times. Moreover, drawdowns approach a steady state(i.e., curves are horizontal and parallel) when t > 1 hr, which is different
from the previous case (the radial hydraulic conductivity). It is worth noting that the effect of further decreasing the length of the sealed section
from 8 to 6 m on drawdown becomes negligible. Differences between
drawdowns are larger for larger lengths of the sealed section (e.g., from10 to 14 m) than for smaller lengths of the sealed section.
The effect of the ratio κ of the lengths of the pumping and injection sections (κ = d1/d3, for d1 ≠ d3 and constant d and d2) is shown in
Figure 6b. As shown in Figure 6b, larger κ ratios result in smaller drawdowns. Drawdowns converge asymptotically toward a steady state for all
scenarios. In addition, it is worth noting that changes in drawdown forsmaller κ are larger than those for larger κ (from 1.0 to 1.8).
3.5. Aquifer Speciﬁc Storage
Figure 7. The effect of the aquifer speciﬁc storage on drawdown(S) with time.

Figure 7 demonstrates how drawdown changes with time for different
values of the aquifer speciﬁc storage (S = 0.0001, 0.0002, 0.0005, and
0.001 m−1) and for the same other conditions as before, that is,3
Q = 60 m /hr, Kr = 0.1 m/hr, Kz = 0.01 m/hr, d = 40 m, d1 = 15 m, d2 = 10 m, and d3 = 15 m, at r = 5 m
and z = 15 m. It can be observed that drawdown decreases with an increase in the speciﬁc storage (S) at earlytimes. In contrast, all drawdown curves for different values of the speciﬁc storage converge asymptotically tothe same value at late times. When other conditions remain the same, an aquifer with large speciﬁc storage
releases more water than an aquifer with small speciﬁc storage. Therefore, as shown in Figure 7, drawdown
decreases with the speciﬁc storage at early times. Since the process of water release from the aquifer storageis almost complete at late times, all drawdown curves eventually approach the same asymptotic value after
about 10 hr of operation and remain stabilized at late times. As a result, the impact of the speciﬁc storage on
drawdown is insigniﬁcant at late times.

3.6. Pumping Time
It is also interesting to analyze how drawdown varies with distance for different pumping times (t = 0.1, 1,
and 10 hr, and steady state). The other parameters in this scenario are again the same as before; that is,
Q = 60 m3/hr, Kr = 0.1 m/hr, Kz = 0.01 m/hr, S = 0.0001 m−1, d = 40 m, d1 = 15 m, d2 = 10 m, andd3 = 15 m at r = 5 m and z = 15 m. The analytical solution Equation 10 for the steady‐state drawdown (seesection 2.4) is utilized as a reference in this section. Figure 8 shows that drawdown changes only slightly withtime for relatively small distances from the well, indicating that groundwater ﬂow quickly reaches aquasi‐steady state in the area close to partially penetrating wells. Additionally, Figure 8 also shows thatthe curves for t = 10 hr and the steady state almost coincide, indicatingthat the entire ﬂow ﬁeld in the single‐well circulation system reaches
steady state relatively quickly.
3.7. Sensitivity Analysis
To estimate the effects of various parameters on drawdown, the sensitivity
analysis has been carried out for the speciﬁc storage S, the length of thesealed section d2, and the radial and vertical hydraulic conductivities, Kr
and Kz, respectively. Using the deﬁnition of the normalized parameter
sensitivity proposed by Kabala (2001), one can obtain

Ym;n ¼ Rn

Figure 8. The effect of the pumping time (t) on drawdown with distance.

TU ET AL.

∂X m
;
∂Rn

(11)

where Ym,n denotes the normalized sensitivity of the mth parameter Rn
at the nth time and Xm is a function of Rn and n. A ﬁnite‐difference
equation proposed by Yeh (1986) (also used by Simunek & van

8 of 14

Water Resources Research

10.1029/2020WR027529

Genuchten, 1996) can be used to approximate the partial derivativein Equation 11:∂X m X m ðRn þ ΔRn Þ − X m ðRn Þ
¼
:
ΔRn
∂Rn

(12)

The small increment ΔRn of 0.01 × Rn is used in this work, similarly
as in Yeh (1986), Simunek and van Genuchten (1996), Huang and
Yeh (2007), and Wen et al. (2013), with the latter two carrying out
the sensitivity analysis of parameters for non‐Darcian ﬂow problems
in a conﬁned aquifer. Other parameters are as before: Q = 60 m3/hr,
Kr = 0.1 m/hr, Kz = 0.01 m/hr, S = 0.0001 m−1, d = 40 m, d1 = 15 m,
d2 = 10 m, and d3 = 15 m, at r = 5 m and z = 15 m. Figure 8 shows thenormalized sensitivities of drawdown with respect to the change of
selected parameters. The sensitivity analysis indicates that theseparameters do not have any impact on drawdown at early times
Figure 9. The sensitivity analysis of drawdown to selected parameters (theand that each parameter has its own inﬂuence period on drawdown.
speciﬁc storage S, the length of the sealed section d2, and the radial and
It also indicates that drawdown is very sensitive to the change in thevertical hydraulic conductivities, Kr and Kz, respectively).
radial hydraulic conductivity. It is worth noting that the radial
hydraulic conductivity Kr has a negative inﬂuence on drawdown in
the period of 3 × 10−3 hr < t < 2 × 10−1 hr. The normalized sensitivity
then starts increasing until about 4 hr when it stabilizes. The sensitivity of drawdown with respect to thelengths of the sealed section d2 starts at about 3 × 10−3 hr, reaches the highest value at about 0.5 hr, and thenremains stable. The sensitivity of drawdown associated with the speciﬁc storage S starts increasing at about
3 × 10−3 hr, reaches the peak value at about 0.1 hr, and then gradually decreases until about 15 hr. Finally,
the inﬂuence of the vertical hydraulic conductivity on drawdown is small and negligible before 0.1 hr, thenkeeps increasing until about 10 hr, and then remains constant.
In practical engineering applications, the design of a single‐well circulation groundwater heat pump needs
to consider many factors, such as local hydrogeological conditions, thermal geological conditions, and the
conﬁguration of the well structure, in order to optimize the efﬁciency of the system and minimize its costs.
As indicated by the sensitivity analysis in Figure 9, it is necessary to investigate hydrogeological conditions
using pumping tests, in situ permeability tests, and geophysical surveys to obtain data related to hydraulicconductivities. The conﬁguration of the well structure in real engineering applications then needs to be
adjusted according to aquifer properties. As shown in Figure 6, shorter lengths of the sealed section result
in larger drawdown, which indicates that the thermal breakthrough will occur sooner and the system willbe less efﬁcient. Although longer lengths of the sealed section could delay the thermal breakthrough,the cost of drilling the well and the electricity consumption for pumping and injection will increase.
Therefore, factors such as the efﬁciency of the system, initial costs, and energy consumption should be comprehensively considered when designing a proper length of the sealed section for a single‐well circulation
groundwater heat pump system.

4. Conclusions
A general analytical mathematical model was developed in this study to investigate the characteristics ofgroundwater ﬂow in a single‐well circulation system. The governing partial differential equation is solved
simultaneously with initial and boundary conditions by applying the Laplace and Fourier cosine transforms.
The analytical solution for transient drawdown in the Laplace domain is obtained and then inverted numerically using the Stehfest method. The steady‐state solution for drawdown in the time domain for a single‐wellcirculation system is also developed. Additionally, the effects of different parameters on drawdown are also
analyzed. The main research conclusions can be drawn as follows:
1. The contours for steady‐state drawdown are symmetric around the horizontal midline of the aquiferand vary tremendously with distance from the well axis (in the case of d1 = d3). The drawdown contours around the sealed section are dense, indicating that the hydraulic gradient in this area is relatively large.

TU ET AL.

9 of 14

Water Resources Research

10.1029/2020WR027529

2. The length of the sealed section d2 has a signiﬁcant inﬂuence on drawdown in the case of d1 = d3. In thed1case of d1 ≠ d3 and a constant length of the sealed section d2, larger ratios of κ (κ ¼ ) result in smallerd3
drawdowns, while changes in drawdown become insigniﬁcant for further increases in κ values.
3. The spatial distribution of drawdown in a single‐well circulation system is sensitive to the radial hydraulic conductivity Kr and the length of the sealed section d2. Each parameter has its own inﬂuence period on
drawdown.
4. Drawdown is not sensitive to parameters S, d2, Kr, and Kz at early times, while it is very sensitive to theseparameters at late times, particularly to the radial hydraulic conductivity, Kr.

Appendix A: Derivation of the Analytical Solution of the Mathematical Model
The governing equation in the Laplace domain is obtained by using the Laplace transform with respect totime t:∂2 sðr; z; pÞ 1 ∂sðr; z; pÞ K z ∂2 sðr; z; pÞSþþ¼
psðr; z; pÞ;
∂r 2
r
∂r
∂z2
Kr
Kr

(A1)

where s is the Laplace transform of drawdown and p is the Laplace transform variable. The Fourier cosine
transform is then applied to the term in Equation A1 with the second‐order partial derivative with respect
to the z coordinate:
Fc



2
d ∂ sðr; z; pÞ∂2 sðr; z; pÞNπz
N 2 π2¼∫
dz ¼ − 2 bs ðr; N; pÞ:
cos
0
2
2
∂z
∂z
d
d

(A2)

Equation A1 is then rewritten by substituting Equation A2 as follows:d2bs ðr; N; pÞ 1 dbs ðr; N; pÞ K z
N 2 π2Spbs ðr; N; pÞ ¼
bs ðr; N; pÞ;
þ
þ
−
r
Kr
dr
Kr
dr 2
d2

(A3)

in which bs denotes the Fourier cosine transform of drawdown and N (N = 1,2,3,……) is the Fourier transform variable. After simplifying Equation A3, one obtains

d2bs ðr; N; pÞ 1 dbs ðr; N; pÞ
K z N 2 π 2 Sp
bs ðr; N; pÞ ¼ 0:
−
þ
þ
r
Kr
dr
K r d2dr 2
Equation A4 can be rewritten as follows by deﬁning δ ¼(A4)

K z N 2 π 2 Spþ :Kr
K r d2d2bs ðr; N; pÞ 1 dbs ðr; N; pÞþ− δbs ðr; N; pÞ ¼ 0:
r
dr
dr 2

(A5)

Equation A5 is a linear second‐order differential equation, which has the general solution as follows:bs ðr; N; pÞ ¼ C1 I 0

pﬃﬃﬃ pﬃﬃﬃ δr þ C 2 K 0 δ r ;(A6)

where C1 and C2 are integration constants, which can be determined using boundary conditions, and I0(x)and K0(x) denote the ﬁrst and second kind modiﬁed Bessel functions of zeroth order with an argument x,
respectively.
Applying the Laplace transform to Equations 3–5, the boundary conditions in the Laplace domain can begiven as follows, respectively:

TU ET AL.

∂sðr; 0; pÞ
¼ 0;
∂z

(A7)

sð∞; z; pÞ ¼ 0;

(A8)10 of 14

Water Resources Research

10.1029/2020WR027529∂sðr; D; pÞ
¼ 0:
∂z

(A9)

Using the Fourier cosine transform to the z coordinate in Equations A7–A9, we can obtainbs ð∞; N; pÞ ¼ 0:

(A10)

Taking into account the property of modiﬁed Bessel functions of the ﬁrst kind, substituting Equation A10
into Equation A6, and having C1 = 0, Equation A6 then can be rewritten asbs ðr; N; pÞ ¼ C2 K 0

pﬃﬃﬃ δr :Applying the Laplace and Fourier cosine transforms to Equation 6, one can obtain, respectively,
8
−Q 1
>
>
ð0 ≤ z ≤ d1 Þ>>>
< 2πK r d1 p
∂sðr; z; pÞlimr→0 r
0¼
ðd1 ≤ z ≤ d1 þ d2 Þ
>
∂r
>
>
Q
1
>
>
ðd1 þ d2 ≤ z ≤ dÞ:
2πK r d3 p
and





dbs ðr:N; pÞ
−dQ
Nπd1
dQ
Nπ ðd1 þ d2 Þ
¼
:
sin
sin
−
r→0
dr
2NpK r d1 π 2
2NpK r d3 π 2
d
d

lim r

(A11)

(A12)

(A13)

From Equations A11 and A13, we can obtain the integration constant C2:

 3

2 
Nπd1
Nπ ðd1 þ d2 Þ
sin
sin
7
dQ 6
d
d
6
7:
C2 ¼
þ
5
2NpK r π 2 4
d3
d1

(A14)

Substituting Equation A14 into Equation A11, we get

 3

2 
Nπd1
Nπ ðd1 þ d2 Þsinsin
7 pﬃﬃﬃ 
dQ 6
d
d
6
7K 0 δr :bs ðr; N; pÞ ¼
þ5
2NpK r π 2 4
d3
d1

(A15)

K z N 2 π 2 SpSpSp
þin Equation A4 becomes. Here, we deﬁne ω ¼ . Thus, theKrKrKr
K r d2governing equation for this case can be written as

When N = 0, the term

d2bs ðr; 0; pÞ 1 dbs ðr; 0; pÞ− ωbs ðr; 0; pÞ ¼ 0:
þ
r
dr
dr 2

(A16)

A similar approach is applied to solve Equation A16 and boundary conditions. Finally, the general solution
can be obtained as follows:bs ðr; 0; pÞ ¼ C2 ′ K 0

pﬃﬃﬃﬃ 
ωr :

(A17)

Applying the Fourier cosine transform to Equation A12, one has
limr→0 r





d1
d
dbs ðr; 0; pÞ
−Q 1
Nπz
Q 1
Nπz
cos
dz
cos
dz
¼ ∫0
þ∫d1 þd2¼ 0:
2πK r d1 p
d
2πK r d3 p
d
dr
N¼0
N¼0

(A18)

Based on Equation A17, one can obtain

TU ET AL.

11 of 14

Water Resources Researchlimr→0 r

10.1029/2020WR027529

pﬃﬃﬃﬃ pﬃﬃﬃﬃ 
dbs ðr; 0; pÞ¼ − limr→0 C2 ′ ωrK 1 ωr ¼ − C2 ′ :
dr

(A19)

Then, combining Equations A18 and A19, the integration constant C2′ in the case of N = 0 can be obtained as
C2 ′ ¼ 0:(A20)

Then, substituting Equation A20 to Equation A17, one hasbs ðr; 0; pÞ ¼ 0:

(A21)

Using the inverse Fourier cosine transform to Equation A15, one can obtain


1
2 ∞
Nπz
:
sðr; z; pÞ ¼ bs ðr; 0; pÞ þ ∑N¼1bs ðr; N; pÞcos
d
d
d

(A22)

Finally, by substituting Equations A15 and A21 to Equation A22, we can obtain an analytical solution in the
Laplace domain.

Appendix B: Derivation of the Steady‐State Solution
Applying the Fourier cosine transform to the term of the second‐order partial derivative of the z coordinate
in Equation 9, one can obtain
Fc



2
d ∂ sðr; z Þ∂2 sðr; zÞNπz
N 2 π2
bsðr; N Þ:
¼
∫
dz
¼
−
cos
0
∂z2
∂z2
d
d2

(B1)

Substituting Equation B1 to Equation 9, one getsd2bsðr; N Þ 1 dbsðr; N Þ K z
N 2 π2
bsðr; N Þ ¼ 0:
þ
þ
−
r
dr
Kr
dr 2
d2
By deﬁning δ′ ¼

(B2)

K z N 2 π2
, Equation B2 can be simpliﬁed as follows:
K r d2d2bsðr; N Þ 1 dbsðr; N Þ− δ′bsðr; N Þ ¼ 0:
þ
r
dr
dr 2

(B3)

The general solution of Equation B3 can be given as follows:bsðr; N Þ ¼ C3 I 0

pﬃﬃﬃﬃ pﬃﬃﬃﬃ δ′ r þ C 4 K 0
δ′ r ;(B4)

where C3 and C4 are integration constants. Using the boundary condition of Equation 3 and considering
the property of the ﬁrst kind of modiﬁed Bessel functions, one can have C3 = 0, then Equation B4 becomesbsðr; N Þ ¼ C 4 K 0

pﬃﬃﬃﬃ δ′ r :(B5)

Applying the Fourier cosine transform to Equation 6, one gets
limr→0 rdbsðr; N Þ
−dQ
Nπd1
dQ
Nπ ðd1 þ d2 Þ
¼
:
sin
sin
−
dr
2NK r d1 π 2
2NK r d3 π 2
d
d

(B6)

By integrating Equations B5 and B6, the constant C4 can be obtained:
TU ET AL.

12 of 14

Water Resources Research

10.1029/2020WR027529


 3


Nπd1
Nπ ðd1 þ d2 Þ
sin
sin
7
dQ 6
d
d
6
7:
C4 ¼
þ
4
5
2
2NK r π
d3
d1
2

(B7)

Substituting Equation B7 into Equation B5, one gets

 3

2 
Nπd1
Nπ ðd1 þ d2 Þsinsin
7 pﬃﬃﬃﬃ′ 
dQ 6
d
d
6
7K 0
bsðr; N Þ ¼
þ
δr :
4
5
2
2NK r π
d3
d1

(B8)

Using the inverse Fourier cosine transform to Equation B8, one can obtain


1
2 ∞
Nπz
:
sðr; zÞ ¼ bsðr; 0Þ þ ∑N¼1bsðr; N Þcos
d
d
d

(B9)

Finally, substituting Equation B8 to Equation B9, an analytical solution for steady‐state conditions can beobtained (Equation 10).

Data Availability Statement
There are no data‐sharing issues because any data used in this manuscript have been provided. Computercodes used in this manuscript are available upon request from the corresponding author.

Acknowledgments
The authors gratefully acknowledge theﬁnancial support from the Key Programof the National Natural Science
Foundation of China (Grant 41430318),
the National Scientiﬁc and TechnicalSupport Program of China (Grant
2016YFC0801801), the National
Natural Science Foundation of China(Grants 41572222, 41602262, and41702261), and the China ScholarshipCouncil (CSC 201806430053).

TU ET AL.

References
Abu‐Nada, E., Akash, B., Al‐Hinti, I., Al‐Sarkhi, A., Nijmeh, S., Ibrahim, A., & Shishan, A. (2008). Modeling of a geothermal standing
column well. International Journal of Energy Research, 32(4), 306–317.
Ataie‐Ashtiani, B., Shafei, B., Rashidian‐Dezfouli, H., & Mohamadzadeh, M. (2012). Capture zone of a partially penetrating well with skineffects in conﬁned aquifers. Transport in Porous Media, 92(2), 437–457.
Banks, D. (2009). Thermogeological assessment of open‐loop well‐doublet schemes: A review and synthesis of analytical approaches.
Hydrogeology Journal, 17(5), 1149–1155.
Barua, G., & Bora, S. N. (2010). Hydraulics of a partially penetrating well with skin zone in a conﬁned aquifer. Advances in Water Resources,33(12), 1575–1587.
Chang, C. C., & Chen, C. S. (2002). An integral transform approach for a mixed boundary problem involving a ﬂowing partially penetrating
well with inﬁnitesimal well skin. Water Resources Research, 38(6), 1071. https://doi.org/10.1029/2001WR001091
Chang, Y. C., & Yeh, H. D. (2009). New solutions to the constant‐head test performed at a partially penetrating well. Journal of Hydrology,369(1), 90–97.
Chen, C. S. (1984). A reinvestigation of the analytical solution for drawdown distributions in a ﬁnite conﬁned aquifer. Water Resources
Research, 20(10), 1466–1468.
Chen, S. J., Wu, L. C., & Liu, C. W. (2010). Analysis of contaminant transport toward a partially penetrating extraction well in an anisotropic aquifer. Hydrological Processes, 24(15), 2125–2136.
Chiu, P. Y., Yeh, H. D., & Yang, S. Y. (2007). A new solution for a partially penetrating constant‐rate pumping well with a ﬁnite‐thickness
skin. International Journal for Numerical and Analytical Methods in Geomechanics, 31(15), 1659–1674.
Crump, K. S. (1976). Numerical inversion of Laplace transforms using a Fourier‐series approximation. Journal of the ACM (JACM), 23(1),89–96.
De Hoog, F. R., Knight, J. H., & Stokes, A. N. (1982). An improved method for numerical inversion of Laplace transforms. SIAM Journal on
Scientiﬁc and Statistical Computing, 3(3), 357–366.
Deng, Z., Rees, S. J., & Spitler, J. D. (2005). A model for annual simulation of standing column well ground heat exchangers. Hvac&R
Research, 11(4), 637–655.
Feng, Q., & Wen, Z. (2016). Non‐Darcian ﬂow to a partially penetrating well in a conﬁned aquifer with a ﬁnite‐thickness skin. Hydrogeology
Journal, 24(5), 1287–1296.
Galgaro, A., & Cultrera, M. (2013). Thermal short circuit on groundwater heat pump. Applied Thermal Engineering, 57(1–2),
107–115.
Hantush, M. S. (1957). Nonsteady ﬂow to a well partially penetrating an inﬁnite leaky aquifer. Proc. Iraqi Sci. Soc., 1(10), 10–19.
Hantush, M. S. (1961). Drawdown around a partially penetrating well. Journal of the Hydraulics Division, 87(4), 83–98.
Huang, Y. C., & Yeh, H. D. (2007). The use of sensitivity analysis in on‐line aquifer parameter estimation. Journal of Hydrology, 335(3–4),
406–418.
Kabala, Z. J. (2001). Sensitivity analysis of a pumping test on a well with wellbore storage and skin. Advances in Water Resources, 24(5),483–504.
Lin, Y. C., Yang, S. Y., Fen, C. S., & Yeh, H. D. (2016). A general analytical model for pumping tests in radial ﬁnite two‐zone conﬁned
aquifers with Robin‐type outer boundary. Journal of Hydrology, 540, 1162–1175.
Luther, K., & Haitjema, H. M. (1999). An analytic element solution to unconﬁned ﬂow near partially penetrating wells. Journal of
Hydrology, 226(3–4), 197–203.

13 of 14

Water Resources Research

10.1029/2020WR027529

Ni, L., Li, H. R., Jiang, Y. Q., Yao, Y., & Ma, Z. L. (2011). A model of groundwater seepage and heat transfer for single‐well ground source
heat pump systems. Applied Thermal Engineering, 31(14–15), 2622–2630.
Novakowski, K. S. (1989). A composite analytical model for analysis of pumping tests affected by well bore storage and ﬁnite thickness skin.
Water Resources Research, 25(9), 1937–1946.
Papadopulos, I. S., & Cooper, H. H. Jr. (1967). Drawdown in a well of large diameter. Water Resources Research, 3(1), 241–244.
Rybach, L. (2015). Innovative energy‐related use of shallow and deep groundwaters—Examples from China and Switzerland. Central
European Geology, 58(1–2), 100–113.
Simunek, J., & van Genuchten, M. T. (1996). Estimating unsaturated soil hydraulic properties from tension disc inﬁltrometer data bynumerical inversion. Water Resources Research, 32(9), 2683–2696.
Sorensen, S.N., & Reffstrup, J. (1992). Prediction of long‐term operational conditions for single‐well groundwater heat pump plants. In:
Proceedings of the 27th Intersociety Energy Conversion Engineering Conference, vol. 4, 109–104, 114.
Stehfest, H. (1970a). Algorithm 368: Numerical inversion of Laplace transforms [D5]. Communications of the ACM, 13(1), 47–49.
Stehfest, H. (1970b). Remark on algorithm 368: Numerical inversion of Laplace transforms [D5]. Communications of the ACM, 13(10),624–625.
Theis, C. V. (1935). The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using
groundwater storage. Eos, Transactions American Geophysical Union, 16(2), 519–524.
Tu, K., Wu, Q., & Sun, H. Z. (2019). A mathematical model and thermal performance analysis of single‐well circulation (SWC) coupled
ground source heat pump (GSHP) systems. Applied Thermal Engineering, 147, 473–481.
Wang, C. T., Yeh, H. D., & Tsai, C. S. (2012). Transient drawdown solution for a constant pumping test in ﬁnite two‐zone conﬁned aquifers.
Hydrology and Earth System Sciences, 16(2), 441–449.
Wen, Z., Liu, K., & Chen, X. L. (2013). Approximate analytical solution for non‐Darcian ﬂow toward a partially penetrating well in a
conﬁned aquifer. Journal of Hydrology, 498, 124–131.
Wen, Z., Huang, G. H, & Zhan, H. B. (2008). An analytical solution for non‐Darcian ﬂow in a conﬁned aquifer using the power law
function. Advances in Water Resources, 31 (1), 44–55.
Wu, Q., Tu, K., Sun, H. Z., & Chen, C. F. (2019). Investigation on the sustainability and efﬁciency of single‐well circulation (SWC)
groundwater heat pump systems. Renewable Energy, 130, 656–666.
Wu, Q., Xu, S. H., Zhou, W. F., & LaMoreaux, J. (2015). Hydrogeology and design of groundwater heat pump systems. Environmental Earth
Sciences, 73(7), 3683–3695.
Xu, S. H., & Rybach, L. (2003). Utilization of shallow resources performance of direct use system in Beijng. Geothermal Resource Council
Transactions, 27, 115–118.
Xu, S. H., & Rybach, L. (2005). Development and application of a new, powerful groundwater heat pump system for space heating andcooling. In Proceedings World Geothermal Congress, pp. 1–5.
Yang, S. Y., & Yeh, H. D. (2005). Laplace‐domain solutions for radial two‐zone ﬂow equations under the conditions of constant‐head andpartially penetrating well. Journal of Hydraulic Engineering, 131(3), 209–216.
Yeh, H. D., Yang, S. Y., & Peng, H. Y. (2003). A new closed‐form solution for a radial two‐layer drawdown equation for groundwater underconstant‐ﬂux pumping in a ﬁnite‐radius well. Advances in Water Resources, 26(7), 747–757.
Yeh, W. W. G. (1986). Review of parameter identiﬁcation procedures in groundwater hydrology: The inverse problem. Water Resources
Research, 22(2), 95–108.
Zeng, Y. F., Zhou, W. F., & LaMoreaux, J. (2017). Single‐well circulation systems for geothermal energy transfer. Environmental Earth
Sciences, 76(7), 296.
Zlotnik, V., Chen, X. H., & Sun, B. (1998). Semi analytical evaluation of three‐dimensional velocity near a partially penetrating well in an
unconﬁned aquifer. Groundwater, 36(3), 514–519.

TU ET AL.

14 of 14