Advances in Water Resources 145 (2020) 103740

Contents lists available at ScienceDirect

Advances in Water Resources
journal homepage: www.elsevier.com/locate/advwatres

An approximate analytical solution for non-Darcian ﬂow in a conﬁned
aquifer with a single well circulation groundwater heat pump system
Kun Tu a,b,c,∗, Qiang Wu a,b, Jirka Simunek c, Ke Zhu a,b, Chaofan Chen e,f, Weiyu Zheng d,
Yifan Zeng a,b, Shengheng Xu a,baCollege of Geoscience and Surveying Engineering, China University of Mining and Technology, Beijing 100083, China
National Engineering Research Center of Coal Mine Water Hazard Controlling, Beijing 100083, Chinac
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAdSchool of Energy and Mining Engineering, China University of Mining and Technology, Beijing 100083, ChinaeHelmholtz Centre for Environmental Research (UFZ), Leipzig 04318, GermanyfApplied Environmental Systems Analysis, Dresden University of Technology, Dresden 01069, Germanyb

a r t i c l e

i n f oKeywords:
Izbash equation
Linearization method
Non-Darcian ﬂow
Geothermal energy
Groundwater heat pump systema b s t r a c t
An analytical model is proposed in this study to describe transient drawdown induced by non-Darcian ﬂow ina conﬁned aquifer with a single well circulation groundwater heat pump system. The Izbash equation and alinearization method are employed to describe non-Darcian ﬂow in the horizontal direction of a conﬁned aquifer
and to approximate the nonlinear term in the governing equation, respectively. By applying a combination of
the Laplace and Fourier cosine transforms, an approximate analytical solution in the Laplace domain is obtained,
which is numerically inverted to obtain transient drawdown in the time domain using the Stehfest algorithm
method. The results of the derived analytical solution for the special case of Darcian ﬂow (m = 1) correspondwell with the existing solution derived using Darcy’s law. The steady-state analytical solution in the time domain
is obtained by applying the Fourier cosine transform. Moreover, the sensitivity analysis is performed to investigate
the inﬂuence of selected parameters, such as the power index m, the radial hydraulic conductivity Kr , the aquifer
speciﬁc storage S, and the length of the sealed section d2 , on drawdown. The results show that each parameterhas its inﬂuence period on drawdown, and that drawdown is more sensitive to the power index m compared toother parameters.

1. Introduction
It is well known that Darcy’s law, which deﬁnes a linear relationshipbetween the ﬂuid ﬂux and the hydraulic gradient, has been commonlyused to address groundwater ﬂow problems in practical engineering applications. However, the ﬂuid ﬂow becomes non-Darcian when the ﬂuid
ﬂux is smaller or larger than a speciﬁc ﬂux interval (Wilkinson, 1955;
Slepicka, 1961; Bear, 1972; Zoorabadi et al., 2015; Zhou et al., 2015).
When the relationship between the ﬂuid ﬂux and the hydraulic gradient deviates from classical Darcy’s law, the ﬂow becomes non-Darcian.
Under these circumstances, researchers point out that Darcy’s law is
not suﬃcient to describe groundwater ﬂow (Miller and Low, 1963;
Bear, 1972; Barker and Herbert, 1977; Sen, 1985). A variety of equations
has been derived to characterize the relationship between the ﬂuid ﬂux
and the hydraulic gradient for non-Darcian ﬂow (Forchheimer, 1901;
Izbash, 1931; Muskat, 1938; Rose, 1951; Escande, 1953). Generally,
these equations can be divided into two categories: power functions and∗polynomial functions (Wen et al., 2013). The ﬁrst category assumes thatthe ﬂuid ﬂux is a power function of the hydraulic gradient, while thesecond category states that the ﬂuid ﬂux can be described as a secondorder polynomial function of the hydraulic gradient.
Among equations of these two categories, the Forchheimer andIzbash equations are widely applied to describe non-Darcian ﬂow (e.g.,Sen, 1989; Sen, 1990; Wen et al., 2008a,b; Moutsopoulos et al., 2009;
Mathias and Todman, 2010; Yeh and Chang, 2013; Sedghi-Asl et al.,2014; Chen et al., 2015; Houben, 2015). Until now, a variety of (approximate or semi-) analytical solutions for non-Darcian ﬂow toward
a fully or partially penetrating well has been derived for these two
equations. For instance, Sen (1989, 1990) presented an analytical solution, which takes into account the non-Darcian ﬂow described using the Forchheimer law toward an inﬁnitesimal well or a large diameter well. An approximate analytical solution considering the eﬀect
of non-Darcian ﬂow, which is deﬁned by the Forchheimer equation,
was derived by Wu (2002). Moutsopoulos and Tsihrintzis (2005) pre-Corresponding author at: College of Geoscience and Surveying Engineering, China University of Mining and Technology, Beijing 100083, China.
E-mail addresses: tukun2014@163.com (K. Tu), wuq@cumtb.edu.cn (Q. Wu).

https://doi.org/10.1016/j.advwatres.2020.103740
Received 31 January 2020; Received in revised form 21 July 2020; Accepted 14 August 2020
Available online 18 August 20200309-1708/© 2020 Elsevier Ltd. All rights reserved.

K. Tu, Q. Wu and J. Simunek et al.

Advances in Water Resources 145 (2020) 103740

Fig. 2. The schematic diagram of the mathematical model for a single-well circulation system.

Fig. 1. A schematic of a single well circulation groundwater heat pump system(following Wu et al., 2019).

sented an approximate analytical solution for nonlinear ﬂow throughporous media described using the Forchheimer equation. Another approximate solution for Forchheimer’s ﬂow toward a well is presented
by Mathias et al. (2008). Applying the Izbash equation, Sen (2000) provided an analytical solution for transient drawdown for non-Darcian
ﬂow to a fully penetrating well in a conﬁned aquifer. Later, based on the
Izbash equation, many eﬀorts have been made by Wen et al. (2008) todevelop analytical solutions describing non-Darcian ﬂow towards a fully
penetrating well. They presented several approximate analytical solutions in the Laplace domain, considering the eﬀect of wellbore storage
and a well radius in diﬀerent aquifer systems. Furthermore, for nonDarcian ﬂow toward a partially penetrating well, Wen et al. (2013) presented a semi-analytical solution in the Laplace domain while neglecting
the eﬀect of the well radius. They also investigated non-Darcian ﬂow toward a larger diameter, partially penetrating well in a conﬁned aquifer(Wen et al., 2014). An analytical solution describing non-Darcian ﬂow toa partially penetrating well in a conﬁned aquifer while considering theeﬀect of a ﬁnite-thickness skin was presented by Feng and Wen (2016).
In addition to analytical solutions, numerical solutions for nonDarcian ﬂow toward wells have also been developed. Wu (2002) employed the ﬁnite diﬀerence method to solve non-Darcian ﬂow describedusing the Forchheimer equation through a fractured reservoir. Based on
the assumption that the Forchheimer equation can describe non-Darcian
ﬂow, Mathias et al. (2008) investigated ﬂow toward a fully penetrating
well in a conﬁned aquifer using the numerical model based on the ﬁnite diﬀerence method. Wen et al. (2009) developed a similar numericalmodel using the ﬁnite diﬀerence method and the assumption that the
Izbash equation can be employed to describe the non-Darcian ﬂow. Finally, Mathias and Wen (2015) performed a numerical simulation of
non-Darcian ﬂow in leaky aquifers using the Forchheimer equation.
A special well structure, referred to as a single well circulation system (Fig. 1), is commonly applied in practical engineering applications.
In a single well circulation system, a single borehole is divided into twopartially penetrating wells by well packers (as shown in Figs. 1 and 2),
which are used to block water injected in the injection well to ﬂow
to the pumping well. The lower partially penetrating well is a pumping well with the pumping rate Q, and the upper partially penetrating

well is an injection well for injecting water at the same rate Q. Thus, thesingle-well circulation system can be considered as a combination of two
partially penetrating wells. To the best of our knowledge, little research
work has been done on groundwater ﬂow in a single well circulation
system. Ni et al. (2011) presented an analytical solution describing thetransient drawdown of groundwater induced by the operation of a single
well circulation groundwater heat pump system. Subsequently, Tu et al.
(2019) and Tu et al. (2020) developed an analytical solution describing
the transient drawdown distribution and obtained a steady-state analytical solution in the time domain using the Laplace transform. It should
be pointed out that these studies on groundwater ﬂow in a single wellcirculation system are based on the assumption of Darcian ﬂow.
However, because of very high hydraulic gradients and enhanced
ﬂow velocities due to the convergence of ﬂow lines (Mathias et al., 2008;
Yeh and Chang, 2013), a considerable number of researchers assumed
that the ﬂow in coarse-grained and fractured media near the wells likely
becomes non-Darcian. The non-Darcian ﬂow may also occur due to highﬂuxes near the wells, especially when the pumping/injection rates are
relatively large (Wan et al., 2013; Wen et al., 2013; Wen et al., 2014;
Feng and Wen, 2016). However, the research on groundwater ﬂow in asingle well circulation system has so far been based on the assumptionof the validity of Darcy’s law, which may not be adequate to describe
the groundwater ﬂow in such a system.
Practical engineering implementations of single well circulation
groundwater heat pump systems require favorable aquifer conditions,which allow water ﬂow at relatively high velocities to continuously
provide a signiﬁcant amount of water. The hydraulic conductivity ofsuch aquifers should be at least 10−3 m/s or higher (Wu et al., 2015;
Rybach, 2015). Besides, relatively large pumping rates (usually about
50–100 m3 /h; Rybach, 2015) are designed to meet the enormous energy demand. Therefore, the ﬂow near the wells in a conﬁned aquifer
with a single well circulation system is particularly prone to become
non-Darcian as a result of high ﬂow velocities and hydraulic gradients.
To the best of our knowledge, most studies on non-Darcian ﬂow
have focused on systems with fully or partially penetrating wells (e.g.,Wu, 2002; Mathias et al., 2008; Wen et al., 2008; Yeh and Chang, 2013;
Wen et al., 2013; Mathias and Wen et al., 2015; Feng and Wen, 2016),
while no analytical model has so far considered non-Darcian ﬂow in asingle well circulation system. It is, therefore, necessary to develop a
new analytical model for investigating the eﬀects of non-Darcian ﬂow
in a conﬁned aquifer with a single well circulation system. This will

K. Tu, Q. Wu and J. Simunek et al.

Advances in Water Resources 145 (2020) 103740

help in improving our understanding of the complex groundwater ﬂow
regime in single well circulation systems, which can be very beneﬁcial
for future investigations of heat convection in such systems.
The aim of this study is to extend the analytical solution based on
Darcy’s law to non-Darcian ﬂow conditions and to investigate the effects of non-Darcian ﬂow in a conﬁned aquifer with a single well circulation system. This could have theoretical signiﬁcance for the applications of single-well circulation groundwater heat pump systems. Dueto relatively high pumping rates and large hydraulic conductivities insingle well circulation systems, the horizontal ﬂow velocity is usually
very high, and non-Darcian ﬂow near the wells likely occurs. The ﬂow
in the horizontal direction is, therefore, assumed to be non-Darcian. On
the other hand, the ﬂow velocity in the vertical direction is relatively
small compared with the horizontal ﬂow velocity. In order to make theanalytical model tractable, the ﬂow in the vertical direction is assumed
to be Darcian. Such assumptions have been often adopted when dealing
with non-Darcian ﬂow problems and when the ﬂow velocity in the vertical direction is considered (Wen et al., 2013; Wen et al., 2014; Feng andWen, 2016).
The Izbash equation is used to describe non-Darcian ﬂow, and a linearization approximation method proposed by Wen et al. (2008) is employed to solve the nonlinear term in the governing equation. This lattermethod has been shown to address this nonlinear problem eﬃciently.
The approximate analytical solution in the Laplace domain is obtainedby applying the Laplace and Fourier cosine transforms. The steady-statesolution in the time domain is also developed using the linearization approximation method and the Fourier cosine transform. Moreover, sensitivity analysis is performed to investigate the inﬂuence of diﬀerent
parameters on drawdown.

2.1. Model setup
Fig. 2 presents a schematic diagram that describes groundwater ﬂow
in a conﬁned aquifer with a single-well circulation system. Non-Darcian
ﬂow is likely to occur because of the high ﬂow velocity in the horizontal direction near the well, especially for a relatively large pumping
rate. In this case, non-Darcian ﬂow near the well needs to be considered. Due to a unique well structure, the ﬂow in the vertical direction,
which is assumed to be Darcian, cannot be neglected. Additional assumptions for this model are made: (1) the conﬁned aquifer is treatedas homogeneous and anisotropic with a uniform thickness d and inﬁnite
in the horizontal direction; (2) the underlying and overlying layers of
the conﬁned aquifer are both impermeable; (3) the speciﬁc discharge
qr in the horizontal direction is assumed to be constant along the screen(the same assumption as in Wen et al., 2013); (4) the pumping and injection rates Q (positive and negative values denote pumping and injection,
respectively) are constant; and (5) the well radius rw is inﬁnitesimallysmall and can be neglected.
Based on the above descriptions and assumptions, the governing
equation for the mathematical model can be written as follows:(1)

where t is time [T], rand z represent the radial and vertical coordinates,
respectively [L], S refers to the speciﬁc storage of an aquifer [L−1 ],qr and qz denote the horizontal and vertical ﬂuid ﬂuxes, respectively[LT−1 ], and s is groundwater drawdown [L].
Before the system starts functioning, groundwater drawdown is assumed to be zero, and the initial conditions for the system can be described as:𝑠(𝑟, 𝑧, 0) = 0𝑠(∞, 𝑧, 𝑡) = 0(2)

The boundary conditions at inﬁnity in the horizontal direction,r = ∞, and the bottom and top of a conﬁned aquifer, z = 0 and z = d,

(3)

and
𝜕𝑠(𝑟, 0, 𝑡)
=0
𝜕𝑧

(4)

and
𝜕𝑠(𝑟, 𝑑, 𝑡)
=0
𝜕𝑧

(5)

where d is the thickness of a conﬁned aquifer [L].
With the assumption that the wellbore storage is neglected, theboundary conditions near the well can be expressed as:()
0 ≤ 𝑧 ≤ 𝑑1()
𝑑 ≤ 𝑧 ≤ 𝑑1 + 𝑑2
( 1)
𝑑1 + 𝑑2 ≤ 𝑧 ≤ 𝑑⎧ −𝑄⎪ 2𝜋𝑑1
lim 𝑟𝑞𝑟 = ⎨0𝑟→0
⎪ 𝑄
⎩ 2𝜋𝑑3(6)

where the pumping and injection rate Q (positive for pumping and negative for injection) are the same and constant [L3 T−1 ]; d1 , d2 , and d3
denote the lengths of well screens for pumping, sealed, and injection
wells, respectively [L].
Izbash (1931) presented a power-law relationship between the hydraulic gradient and water ﬂux, which can be employed to describe thenon-Darcian ﬂow. In this model, the horizontal speciﬁc discharge described by the Izbash equation for two-dimensional ﬂow, same as in
Wen et al. (2013), can be given as:
𝑞𝑟 ||𝑞𝑟 ||2. Mathematical model𝜕 𝑞𝑟 𝑞𝑟 𝜕 𝑞𝑧
𝜕𝑠(𝑟, 𝑧, 𝑡)
+
+
=𝑆
𝜕𝑟
𝑟
𝜕𝑧
𝜕𝑡

are given as follows, respectively:

𝑚−1

= 𝐾𝑟

𝜕𝑠(𝑟, 𝑧, 𝑡)𝜕𝑟(7)

where Kr is the (quasi) (Wen et al. 2006) radial hydraulic conductivity,
and the non-Darcian exponent m in the Izbash’s equation is an empiricalcoeﬃcient that ranges between 1 and 2 and denotes the degree of deviation from linearity (Izbash, 1931; Bordier and Zimmer, 2000). When m< 1, the ﬂow is treated as pre-linear ﬂow, and when 1 < m ≤ 2, the ﬂow
becomes post-linear ﬂow (Izbash, 1931; Wen et al., 2006). It should bepointed out that Eq. (7) reduces to Darcy’s law when m = 1.
Due to various factors, such as shapes and arrangements of pores,
and complex ﬂow paths, the non-Darcian exponent m varies in space andtime (Bordier and Zimmer, 2000; Wen et al., 2006). It is thus quite difﬁcult to determine and quantify its value in practical engineering applications. Although some researchers (Soni et al., 1978; Watanable, 1982;
Chen et al., 2015) have attempted to determine the value of the nonDarcian exponent m, this issue has not yet been well resolved. In orderto make the proposed models tractable, the non-Darcian exponent m isoften assumed to be constant when dealing with non-Darcian ﬂow problems (e.g., Sen, 2000; Wen et al., 2006, 2008; Quinn et al., 2011, 2013;
Wen et al., 2013; Liu et al., 2016). However, adopting a constant nonDarcian exponent m in the Izbash’s equation imposes some limitations
when describing the non-Darcian ﬂow. To overcome these limitations
still remains an open issue and a very interesting topic for future investigations.
The ﬂow in the vertical direction is assumed to be Darcian, which is
described as:
𝑞𝑧 = 𝐾𝑧𝜕𝑠(𝑟, 𝑧, 𝑡)𝜕𝑧(8)where Kz represents the vertical hydraulic conductivity.
Since the ﬂow in the horizontal direction is against the directionof the r coordinate, the speciﬁc discharge qr is always negative. Thus,
Eq. (7) can be furtherly reduced as:
(

−𝑞𝑟

)𝑚

= 𝐾𝑟

𝜕𝑠(𝑟, 𝑧, 𝑡)𝜕𝑟(9)K. Tu, Q. Wu and J. Simunek et al.

Advances in Water Resources 145 (2020) 103740

2.2. Approximate analytical solution

2.3. An analytical solution at a steady state

The governing equation of the mathematical model is obtained bysubstituting Eqs. (8) and (9) into Eq. (1):When the operation time of the system is long enough, the changes indrawdown near wells become stable, and the ﬂow reaches steady-state
𝑟,𝑧,𝑡)
conditions. In this case 𝜕𝑠(𝜕𝑡= 0 and Eq. (12) becomes:

𝐾𝑟𝐾𝜕 2 𝑠(𝑟, 𝑧, 𝑡)𝜕𝑠(𝑟, 𝑧, 𝑡)𝜕 2 𝑠(𝑟, 𝑧, 𝑡)𝜕𝑠(𝑟, 𝑧, 𝑡)+ ( 𝑟)
+ 𝐾𝑧
=𝑆
(
)𝑚−1
2
𝑚
−1
𝜕𝑟
𝜕𝑡
𝜕𝑟
𝜕 𝑧2
𝑚 −𝑞𝑟
𝑟 −𝑞𝑟
(10)
It is noteworthy that Eq. (10) is a nonlinear equation due to the nonlinear term ( − qr )m − 1 . Solving such a nonlinear equation may not be
tractable using rigorous mathematical means. However, a linearization
approximation method proposed by Wen et al. (2008) may be used toaddress this problem. This method has been successfully applied to ﬁnd
analytical solutions for other non-Darcian problems (e.g., Wen et al.,2008; Wen et al., 2013; Feng and Wen, 2016). By applying the linearization approximation method, the nonlinear term can be approximately
described as:
(
)𝑚−1
(
)𝑚−1
𝑄
≈
(11)
−𝑞𝑟
2𝜋𝑟𝑑
(

𝑄
2𝜋𝑟𝑑

)𝑚−1

( √ ) 2 [(𝑠̄(𝑟, 𝑧, 𝑝) =

3−𝑚

𝛼
3−𝑚

∞
∑

4

)𝑚

√

𝑝𝑑 𝐾𝑟 𝛼Γ( √ ) 2 [(
+

𝑄
2𝜋𝑑1

(12)

𝜑

3−𝑚

3−𝑚

𝑄
2𝜋𝑑1

)𝑚

𝑑1 −

(

2
3−𝑚

𝑠𝑖𝑛

(

(
)

𝑛𝜋𝑑1
𝑑

√

𝑛=1

×𝐾 1−𝑚
3−𝑚

( √
2 𝜑
3−𝑚

3−𝑚

𝑟 2

)

𝑄
2𝜋𝑑3

(

𝑛𝜋𝑧
𝑑

𝑑3

]
1−𝑚

𝑟 2 𝐾 1−𝑚3−𝑚)𝑛𝜋𝑝𝐾𝑟 𝜑Γ
𝑐 𝑜𝑠

)𝑚

(
+

(

𝑄
2𝜋𝑑3

2
3−𝑚

)

)𝑚

𝑠𝑖𝑛

(

+

𝜔
3−𝑚

( √
2 𝛼
3−𝑚

𝑛𝜋 (𝑑1 +𝑑2 )
𝑑

3−𝑚

)

𝑟 2

3−𝑚

𝑄
2𝜋𝑑1

)𝑚

𝑛=1
1−𝑚

where 𝜔 =Eq. (12) is a linear partial diﬀerential equation, which can be solved
using the Laplace and Fourier cosine transforms. A detailed derivation of the solution of this partial diﬀerential equation can be seen in
Appendix A. An approximate analytical solution of non-Darcian ﬂow for
the drawdown in the Laplace domain is given as:2( √ ) 2 [(
∞ 4
∑

( √
2 𝜔

3−𝑚

𝜕 2 𝑠(𝑟, 𝑧, 𝑡)
𝜕 𝑧2

(

𝑄
2𝜋𝑟𝑑

)𝑚−1

𝜕 2 𝑠(𝑟, 𝑧)
=0
𝜕 𝑧2

(14)

For this linear partial diﬀerential equation, the analytical solution ata steady state in the time domain is obtained by applying the Fouriercosine transform (see Appendix B):
(
)𝑚
(
)𝑚
𝑄
𝑄
𝑑1 − 2𝜋𝑑
𝑑3
2𝜋𝑑1
3
𝑟1−𝑚
(𝑟, 𝑧) =()
𝑑 𝐾𝑟 Γ 3−2𝑚 (3 − 𝑚)

× 𝑟 2 𝐾 1−𝑚

Substituting Eq. (11) into Eq. (10), one can obtain:
𝜕 2 𝑠(𝑟, 𝑧, 𝑡) 𝑚 𝜕𝑠(𝑟, 𝑧, 𝑡) 𝑚𝐾𝑧
+
+
𝑟
𝜕𝑟
𝐾𝑟
𝜕 𝑟2
(
)𝑚−1
𝜕𝑠(𝑟, 𝑧, 𝑡)
𝑄𝑚𝑆≈
𝐾𝑟 2𝜋𝑟𝑑𝜕𝑡𝜕 2 𝑠(𝑟, 𝑧) 𝑚 𝜕𝑠(𝑟, 𝑧) 𝑚𝐾𝑧
+
+
𝑟 𝜕𝑟
𝐾𝑟
𝜕 𝑟2

3−𝑚

3−𝑚

𝑟 2

)

𝑠𝑖𝑛

(

𝑛𝜋𝑑1
𝑑

)

(
+

𝑄
2𝜋𝑑3

)𝑚

𝑠𝑖𝑛

)
√ (
𝑛𝜋𝐾𝑟 𝜔Γ 3−2𝑚

cos

(

𝑛𝜋𝑧
𝑑

(

𝑛𝜋 (𝑑1 +𝑑2 )
𝑑

)]

)
(15)

𝑚𝐾𝑧 𝑛2 𝜋 2 𝑄 𝑚−1() , and s is the drawdown in the time domain.
𝐾𝑟 𝑑 2 2𝜋𝑑2.4. Simpliﬁcation of the analytical solution
When m = 1, the ﬂow becomes Darcian ﬂow and the analytical solution of Eq. (13) is as follows:()()
𝑛𝜋 (𝑑1 +𝑑2 ) ⎤𝑛𝜋𝑑⎡
sin 𝑑 1sin⎥𝑑𝑄 ⎢𝑠(𝑟, 𝑧, 𝑝) =
+
⎥
2⎢
𝑑
𝑑
np
𝐾
𝜋
1
3
𝑟
⎢
⎥
𝑛=1
⎣
⎦
(√ )
(
)
𝑛𝜋𝑧
𝛽𝑟 cos
×𝐾0
𝑑
∞
∑

𝐾

(16)

where 𝛽 = 𝐾𝑧 𝑛𝑑𝜋2 + 𝑝𝑆
, and K0 (·) is the modiﬁed Bessel function of the
𝐾
2 2

𝑟

)]
1−𝑚

𝑟 2

𝑟

second kind with an order of zero. The Eq. (16) is the same as the analytical solution derived by Tu et al. (2020) for Darcian ﬂow in a conﬁned
aquifer with a single-well circulation system.
3. Results and discussion)(13)

𝑚𝐾𝑧 𝑛2 𝜋 2 +𝑚𝑝𝑆 𝑑 2 𝑄 𝑚−1
where 𝜑 =( 2𝜋𝑑 ) , 𝛼 = 𝑚𝑝𝑆( 𝑄 )𝑚−1 ; 𝑠̄ is the solution of
𝐾𝑟 2𝜋𝑑𝐾𝑟 𝑑 2

drawdown in the Laplace domain; Γ(·) denotes the gamma function; K𝜈 (·)
represents the modiﬁed Bessel function of the second kind with an orderof 𝜈; and n and p are the Fourier and Laplace variables, respectively.
Due to the presence of the Bessel function K𝜈 (·) term, it may not be
possible to obtain the analytical solution of Eq. (13) in the time domain.
However, several numerical inverse methods (e.g., Stehfest, 1970a,b;
Crump, 1976; De Hoog et al., 1982) can be employed to transform theanalytical solution in the Laplace domain to the time domain. In thiswork, an inverse Laplace transform is numerically conducted using theStehfest method, which has been shown to address such inverse problems eﬃciently and accurately (Wen et al., 2008).
A program has been developed in MATLAB to calculate the solution
for drawdown in the time domain based on the Stehfest method. The
number of terms N, which are used in the inversion, needs to be speciﬁed before performing it. Multiple numerical studies have shown thatthe solution gives accurate and stable results with N = 18. Similarly,
the inﬁnite series of Eq. (13) converges when n ≥ 20. Therefore, to increase computational eﬃciency, N = 18 and n = 20 are employed in thefollowing calculations.

3.1. Drawdown versus time as a function of the power index
The inﬂuence of the power index m (from 1 to 2) on drawdown observed at r = 5 m and z = 15 m is analyzed ﬁrst when Q = 60 m3 /h,
d = 40 m, d1 = 15 m, d2 = 10 m, d3 = 15 m, S = 0.0001 m−1 , Kr = 0.1(m/h)m , and Kz = 0.01 m/h. For the case of m = 1 (Darcian ﬂow), it
can be observed in Fig. 3 that the drawdown calculated using the new
analytical solution matches well the results of the analytical solution of
Ni et al. (2011). The results in Fig. 3 indicate that the power index inﬂuences drawdown during the entire pumping period. A larger value ofthe power index results in larger drawdown at the beginning and smaller
drawdown at later times (Fig. 3). A similar inﬂuence can also be found in
Wen et al. (2013). A larger power index in the Izbash equation indicatesmore signiﬁcant deviations from Darcian ﬂow and results in greater ﬂow
turbulence. At early times, the elastic release process of an aquifer may
be accelerated for a larger power index. Consequently, larger drawdown
is observed at early times for a larger power index. As pumping time increases, drawdowns for diﬀerent power indices approach a steady state
when t > 10 h, indicating that the process of elastic release of an aquifer
is approximately complete. Meanwhile, pumped water comes mainly
from the region relatively far away from the pumping well. Larger ﬂow
turbulence induced by a larger power index is likely to result in larger

K. Tu, Q. Wu and J. Simunek et al.

Fig. 3. The inﬂuence of the power index (m) on drawdown versus time.

Advances in Water Resources 145 (2020) 103740

Fig. 5. The inﬂuence of the aquifer speciﬁc storage (S) on drawdown versus
time.

a larger hydraulic conductivity at early times can be found. At late times,
when the ﬂow reaches a quasi-steady state, recharge occurs from areas
relatively far away from the well. Smaller hydraulic conductivities willthen result in a larger inﬂuence range of the depression cone and larger
drawdown. Additionally, it is worth noting that drawdown approaches
a steady state more quickly with increasing radial hydraulic conductivities (Fig. 4).
3.3. Drawdown versus time as a function of the aquifer speciﬁc storage

Fig. 4. The inﬂuence of the radial hydraulic conductivity (Kr ) on drawdown
versus time.

recharge. Therefore, it can be seen that drawdown is smaller for a larger
power index at late times. Moreover, Fig. 3 also indicates that drawdown
reaches a steady state more quickly when the power index is larger.
3.2. Drawdown versus time as a function of the radial hydraulicconductivity
Fig. 4 depicts temporal drawdown distribution curves observed at
r = 5 m and z = 15 m for diﬀerent radial hydraulic conductivities Krfrom 0.01 to 0.1 m/h. Other parameters are as follows: Q = 60 m3 /h,
m = 1.5, d = 40 m, d1 = 15 m, d2 = 10 m, d3 = 15 m, S = 0.0001 m−1 ,and Kz = 0.01 m/h. It is interesting to note that the inﬂuence of the radial hydraulic conductivity on drawdown is similar to that of the powerindex. As indicated in Fig. 4, drawdown increases with the radial hydraulic conductivity at the beginning, while smaller drawdown is obtained for larger radial hydraulic conductivities at late times. The reasons for this phenomenon can be explained as follows. At early times,
larger hydraulic conductivities induce faster groundwater ﬂow. Due tothe faster spreading of the suction zone in the depression cone, groundwater can be replenished promptly. Consequently, larger drawdown for

The curves in Fig. 5 show the temporal drawdown distributions observed at r = 5 m and z = 15 m calculated using the aquifer speciﬁcstorage S from 0.0001 to 0.001 m−1 . The other parameters used to calculate the curves in the ﬁgure are as follows: Q = 60 m3 /h, m = 1.5,
d = 40 m, d1 = 15 m, d2 = 10 m, d3 = 15 m, Kr = 0.1 (m/h)m , andKz = 0.01 m/h. It is demonstrated in Fig. 5 that drawdown is smaller
for a larger aquifer speciﬁc storage and vice versa at early times. A
larger aquifer speciﬁc storage indicates that the aquifer can release moregroundwater when other conditions remain the same. Moreover, at latetimes (t > 1 h), when the process of the storage release from the aquifer
is almost complete, drawdown curves merge for diﬀerent aquifer speciﬁc storages (Fig. 5). Drawdowns approach a steady state at late times,
and the inﬂuence of diﬀerent aquifer speciﬁc storages on drawdown isnegligible.
3.4. Drawdown versus time as a function of the size of the sealed section
Fig. 6 shows drawdowns versus time observed at r = 5 m and z = 15 m
for diﬀerent lengths of the sealed section d2 from 6 to 14 m. In thiscase, the thickness of the aquifer is kept constant, and the lengths ofwell screens for pumping and injection wells (d1 and d3 , respectively)
are correspondingly adjusted when the length of the sealed sectionchanges. Other parameters used in this case are: Q = 60 m3 /h, m = 1.5,
S = 0.0001 m−1 , Kr = 0.1 (m/h)m , and Kz = 0.01 m/h. The curves shownin Fig. 6 indicate that changes in drawdown are signiﬁcant when thelength of the sealed section is reduced from 14 to 8 m. Also, larger
lengths of the sealed section result in smaller drawdowns. Since theaquifer thickness and the pumping rate remain constant, the boreholewith a longer length of the sealed section has a shorter length of thepumping screen and a correspondingly larger pumping rate for a unit
length of the well screen. It is thus not surprising that larger drawdowns
are obtained for shorter lengths of the sealed section. However, the inﬂuence of the further reduction in the length of the sealed section (fromK. Tu, Q. Wu and J. Simunek et al.

Fig. 6. The inﬂuence of the sealed section (d2 ) on drawdown versus time.

Advances in Water Resources 145 (2020) 103740

Fig. 8. The inﬂuence of the power index (m) on drawdown versus distance.

Fig. 9. Drawdown contours at a steady state.

Fig. 7. The inﬂuence of the pumping time (t) on drawdown versus distance.

8 to 6 m) on drawdown is almost imperceptible. Also, diﬀerences indrawdown are smaller for smaller lengths of the sealed section (e.g.,
from 6 to 10 m) than for larger lengths of the sealed section. Moreover,when t > 1 h, drawdowns approach a steady state, and the curves areparallel.
3.5. Drawdown versus distance as a function of the pumping time
It is also important to analyze drawdown versus the radial distancefor diﬀerent pumping times (for t = 0.1, 1, 10 h, and at a steady state).
The other parameters are as follows: Q = 60 m3 /h, m = 1.5, d = 40 m,
d1 = 15 m, d2 = 10 m, d3 = 15 m, S = 0.0001 m−1 , Kr = 0.1 (m/h)m ,and Kz = 0.01 m/h. The steady-state analytical solution, Eq. (15), isemployed in this ﬁgure as a reference. It can be observed in Fig. 7 thatdrawdown curves merge for diﬀerent pumping times when the radial
distance r < 1 m, indicating that groundwater ﬂow in the area aroundthe pumping well reaches a steady state quickly. For the region further
away from the well (r > 1 m), drawdown curves start deviating fromeach other. This is because the ﬂow in this region is still unsteady andwill take longer to reach a steady state. Moreover, drawdown at t = 10 h
is very close to that at a steady state, which means that groundwater ﬂow

in a conﬁned aquifer for a single-well circulation system approaches asteady state relatively quickly.
3.6. Drawdown versus distance as a function of the power index
Temporal drawdown distribution curves versus a radial distance fordiﬀerent values of the power index m from 1 to 2 are shown in Fig. 8.
The other parameters used in this ﬁgure are the same as before, i.e.,
Q = 60 m3 /h, d = 40 m, d1 = 15 m, d2 = 10 m, d3 = 15 m, t = 10
h, S = 0.0001 m−1 , Kr = 0.1 (m/h)m , and Kz = 0.01 m/h. It can beobserved in Fig. 8 that a smaller value of the power index results ina larger drawdown for r < 100 m. However, at larger radial distances,drawdown for smaller values of the power index is smaller than for its
larger values. This is consistent with the results in Fig. 3, implying thatﬂow near the pumping well reached a steady state when the pumping
time t = 10 h, while a steady state was not reached in the area far awayfrom the pumping well (Fig. 8).
3.7. Drawdown contours at a steady state
Drawdown contours displayed in Fig. 9 are obtained using the analytical solution of Eq. (15) for steady-state conditions. The following
parameters are used in this scenario: Q = 60 m3 /h, m = 1.5, d = 40 m,
d1 = 15 m, d2 = 10 m, d3 = 15 m, S = 0.0001 m−1 , Kr = 0.1 (m/h)m ,K. Tu, Q. Wu and J. Simunek et al.

Advances in Water Resources 145 (2020) 103740

Table 1
Ranges of parameters used in the global sensitivity
analysis.
Uncertain parameter

Range

Power index m
Radial hydraulic conductivity Kr
Speciﬁc storage S
Length of the sealed section d2

1–2
0.01–0.1
0.00001–0.0001
6–14

and Kz = 0.01 m/h. As shown in Fig. 7, changes in drawdown are verylarge for small r (less than 0.5 m), causing contours to be too dense to
be displayed in Fig. 9. Therefore, the starting point on a horizontal axis
in Fig. 9 is at r = 0.5 m. Due to ﬂow symmetry, contours for steady-state
drawdown are symmetric along a horizontal line at z = 20 m (the middle of a conﬁned aquifer) where drawdown is equal to zero (Fig. 9) andvary tremendously along the well axis. When a single-well circulation
system functions properly, pressure diﬀerences between the injection
zone and the suction zone will yield relatively large hydraulic gradientsnear the wells. Thus, drawdown contours near the borehole are muchdenser than further away from the pumping or injection wells.
3.8. Sensitivity analysis
The global sensitivity analysis is an eﬀective method that can be employed to analyze the sensitivity of the model output to changes in each
input parameter (Saltelli et al., 2000). In particular, the Sobol’s indiceshave been widely employed as sensitivity measures in many applications involving hydrological (van Werkhoven et al., 2009; Ciriello and
Di Federico, 2012; Di Federico and Ciriello, 2012; Brunetti et al., 2016)
and environmental (Nossent et al., 2011; Pianosi et al., 2015) models.
This is because no assumptions of linearity or monotonicity are required
for the adopted interpretative model (Sobol’, 1993; Archer et al., 1997).
These indices can provide accurate information about the model output
variance related to a single parameter or associated with interactions ofmultiple parameters. A more detailed description of the Sobol’ method
can be found in Sobol’ (2001).
According to the deﬁnition given by Sobol’ (1993), the Sobol’ sensitivity indices can be expressed as the ratio of the partial model varianceto the total model variance:
First -order 𝑆𝑖 =𝑉𝑖
𝑉Second-order 𝑆𝑖𝑗 =

(17)
𝑉𝑖𝑗
𝑉

Fig. 10. The total sensitivity indices as a function of time for various parameters (the power index m, the aquifer speciﬁc storage S, the radial hydraulicconductivity Kr , and the length of the sealed section d2 ).

Fig. 10 presents temporal changes in total sensitivity indices for selected input parameters. It can be seen that two parameters show a signiﬁcant inﬂuence on the model output variance. These are the powerindex m and the radial hydraulic conductivity Kr . The total sensitivity
index for the power index m is initially small, then increases with time,reaches a peak value near t = 20 h, and then remains constant withtime. The radial hydraulic conductivity Kr shows a similar trend as thepower index m, except that it has a lower peak value when reaching a
steady state. The total sensitivity index for the speciﬁc storage S has initially a very high value, then sharply decreases until near t = 10 h, andthen gradually drops to almost zero after t> 40 h. This indicates thatthe speciﬁc storage S has a negligible inﬂuence on the model output.
Finally, the overall behavior of the total sensitivity index for the lengthof the sealed section d2 is remarkably similar to that of the speciﬁc storage S, except that the total sensitivity indices reach a constant minimum
value after large time, which is much larger compared with that for the
speciﬁc storage S.
4. Conclusions(18)

where Vi represents the partial variance of the model associated with
the ith parameter, Vij represents the partial variance of the model associated with the interaction of the ith and jth parameters, and V is thetotal variance. The ﬁrst-order, or principal sensitivity, indices Si , called
the “main eﬀect indices”, are used to describe contributions of a single
parameter to the total variance of the model output. The inﬂuence of
the interaction of two input parameters on the model output varianceis denoted as the second-order index, Sij . The total sensitivity indices,
𝑆𝑇𝑖 (Ciriello and Di Federico, 2012), represent the main eﬀect of a given
parameter (the ith parameter) and all its interactions with other parameters (up to the kth order, k denotes the number of input parameters),
which can be expressed as follows:∑
Total 𝑆𝑇𝑖 = 𝑆𝑖 +𝑆𝑖𝑗 + ⋯
(19)𝑗≠𝑖

The evaluation of sensitivity indices is performed for the following parameters: the power index m, the radial hydraulic conductivity
Kr , the aquifer speciﬁc storage S, and the length of the sealed section
d2 . Table 1 shows the range of all evaluated parameters. Other conditions used when conducting the global sensitivity analysis are as follows:
Q = 60 m3 /h, d = 40 m, d1 = d3 = (40-d2 )/2, and Kz = 0.01 m/h.

In this work, an analytical model has been developed to describe
non-Darcian ﬂow in a conﬁned aquifer with a sing-well circulation
system. Non-Darcian ﬂow in the horizontal direction is described using the Izbash equation, and a linearization approach presented by
Wen et al. (2008) is employed to approximate the nonlinear term in thegoverning equation. Then, using the Laplace and Fourier cosine transforms, the analytical solution in the Laplace domain is obtained, whichis numerically inverted into the time domain using the Stehfest method.
The main ﬁndings from this work can be drawn as follows:(1) Larger values of the power index m and the radial hydraulic conductivity Kr result in larger drawdown at early times and smaller
drawdown at late times. Groundwater ﬂow approaches a steadystate faster for larger values of the power index m and the radial
hydraulic conductivity Kr .
(2) Larger values of the aquifer speciﬁc storage S result in smaller
drawdown, while the inﬂuence of S on drawdown is negligible atlate times.
(3) Smaller lengths of the sealed section d2 result in larger drawdown. However, a further reduction in the length of the sealed
section d2 below a certain critical length leads to imperceptiblechanges in drawdown.

K. Tu, Q. Wu and J. Simunek et al.

(4) The contours for steady-state drawdown vary signiﬁcantly along
the well axis and are symmetric around a horizontal line in the
middle between the injection and suction sections (z = 20 m inour example) where drawdown is equal to zero.
(5) Each parameter has its inﬂuence period on drawdown. Drawdown is the most sensitive to the length of the sealed section d2
and the aquifer speciﬁc storage S at early times, and the powerindex m (in particular) and the radial hydraulic conductivity Kr
at late times.
Declaration of Competing Interest
We declare that we have no ﬁnancial and personal relationships withother people or organizations that can inappropriately inﬂuence ourwork, there is no professional or other personal interest of any nature orkind in any product, service and/or company that could be construed as
inﬂuencing the position presented in, or the review of, the manuscript
entitled, “An approximate analytical solution for non-Darcian ﬂow in a
conﬁned aquifer with a single well circulation groundwater heat pumpsystem”.
CRediT authorship contribution statement
Kun Tu: Conceptualization, Methodology, Writing - original draft,
Software. Qiang Wu: Supervision, Project administration, Funding acquisition. Jirka Simunek: Supervision, Resources, Writing - review &editing. Ke Zhu: Methodology, Software. Chaofan Chen: Conceptualization, Formal analysis. Weiyu Zheng: Software, Formal analysis.
Yifan Zeng: Conceptualization, Investigation. Shengheng Xu: Formal
analysis.
Acknowledgments
The authors gratefully acknowledge the ﬁnancial support from the
Key Program of the National Natural Science Foundation of China (Grant
No. 41430318), the National Scientiﬁc and Technical Support Programof China (Grant No. 2016YFC0801801), the National Natural Science
Foundation of China (Grant No. 41572222, 41602262, 41702261), andthe China Scholarship Council (No. 201806430053).
Appendix A. Derivation of the approximate analytical solution
The governing equation after applying the linearization method canbe given as:()𝑚−1 2𝜕 2 𝑠(𝑟, 𝑧, 𝑡) 𝑚 𝜕𝑠(𝑟, 𝑧, 𝑡) 𝑚𝐾𝑧
𝜕 𝑠(𝑟, 𝑧, 𝑡)
𝑄
+
+
2
𝑟
𝜕𝑟
𝐾𝑟 2𝜋𝑟𝑑
𝜕𝑟
𝜕 𝑧2
(
)𝑚−1
𝜕𝑠(𝑟, 𝑧, 𝑡)
𝑄
𝑚𝑆
=
(A1)
𝐾𝑟 2𝜋𝑟𝑑𝜕𝑡
The Laplace transform with respect to time t is applied to Eq. (A1) toget:()𝑚−1 2𝜕 2 𝑠̄(𝑟, 𝑧, 𝑝) 𝑚 𝜕 𝑠̄(𝑟, 𝑧, 𝑝) 𝑚𝐾𝑧𝜕 𝑠̄(𝑟, 𝑧, 𝑝)
𝑄
+
+
2
𝑟
𝜕𝑟
𝐾𝑟 2𝜋𝑟𝑑
𝜕𝑟
𝜕 𝑧2
(
)𝑚−1
𝑄
𝑚𝑆
=
𝑝𝑠̄(𝑟, 𝑧, 𝑝)(A2)
𝐾𝑟 2𝜋𝑟𝑑where p and 𝑠̄ are the Laplace variable and drawdown in the Laplacedomain, respectively.
Then, the Fourier cosine transform is used to the second-order partial
derivative with respect to coordinate z in Eq. (A2) to obtain:
[ 2
]
(
)
𝑑 2
𝜕 𝑠̄(𝑟, 𝑧, 𝑝)𝜕 𝑠̄(𝑟, 𝑧, 𝑝)𝑛𝜋𝑧𝑛2 𝜋 2𝑑𝑧 = −𝐹𝑐=𝑐 𝑜𝑠𝑠̂(𝑟, 𝑛, 𝑝)
(A3)
2
2
∫0
𝑑
𝜕𝑧
𝜕𝑧
𝑑2
where n (n = 0, 1, 2, 3, ⋅⋅⋅⋅⋅⋅) is the Fourier variable, and 𝑠̂(𝑟, 𝑛, 𝑝)
denotes the Fourier cosine transform of drawdown.

Advances in Water Resources 145 (2020) 103740

By substituting Eq. (A3) to Eq. (A2) we get:𝑚𝐾𝑧 𝑛2 𝜋 2 + 𝑚𝑝𝑆 𝑑 2
𝑑 2 𝑠̂(𝑟, 𝑛, 𝑝) 𝑚 𝑑 𝑠̂(𝑟, 𝑛, 𝑝)
+
=
2
𝑟
𝑑𝑟
𝑑𝑟
𝐾𝑟 𝑑 2

(

𝑄
2𝜋𝑟𝑑

)𝑚−1

𝑠̂(𝑟, 𝑛, 𝑝)(A4)

By deﬁning 𝜑 =𝑚𝐾𝑧 𝑛2 𝜋 2 +𝑚𝑝𝑆 𝑑 2 𝑄 𝑚−1( 2𝜋𝑑 ) , Eq. (A4) can be reduced to:
𝐾𝑟 𝑑 2𝑑 2 𝑠̂(𝑟, 𝑛, 𝑝) 𝑚 𝑑 𝑠̂(𝑟, 𝑛, 𝑝)+− 𝜑𝑟1−𝑚 𝑠̂(𝑟, 𝑛, 𝑝) = 0
𝑟
𝑑𝑟
𝑑 𝑟2

(A5)

Eq. (A5) is a linear diﬀerential equation of the second order. The
general solution of Eq. (A5) can be described as:
[
( √
)
( √
)]
2 𝜑 3−𝑚2 𝜑 3−𝑚1−𝑚𝑠̂(𝑟, 𝑛, 𝑝) = 𝑟 2 𝐶1 𝐼 1−𝑚
𝑟 2+ 𝐶2 𝐾 1−𝑚
𝑟 2
(A6)
3−𝑚
3−𝑚
3−𝑚
3−𝑚
where 𝐼 1−𝑚 (⋅) and 𝐾 1−𝑚 (⋅) are the modiﬁed Bessel functions of the ﬁrst
3−𝑚

3−𝑚

𝑚
and second kind with the order of 1−, respectively. The integration
3−𝑚constants C1 and C2 can be determined according to boundary conditions.
By applying the Laplace and Fourier cosine transforms to Eq. (3) one
can obtain:𝑠̂(∞, 𝑛, 𝑝) = 0(A7)Considering the property of the modiﬁed Bessel function of the ﬁrst
kind and substituting Eq. (A7) to Eq. (A6) one gets:
𝐶1 = 0(A8)

Eq. (A6) can then be rewritten as:( √)2 𝜑 3−𝑚1−𝑚𝑠̂(𝑟, 𝑛, 𝑝) = 𝐶2 𝑟 2 𝐾 1−𝑚
𝑟 2
3−𝑚
3−𝑚

(A9)

Combining the Izbash equation and the linearization method, theboundary conditions of Eq. (6) can be given as:
(
)𝑚
(
)
𝑄
⎧ − 1
0 ≤ 𝑧 ≤ 𝑑1
𝐾𝑟 2𝜋𝑟(𝑑1)𝜕𝑠(𝑟, 𝑧, 𝑡) ⎪
lim
=⎨
(A10)
0
𝑑 ≤ 𝑧 ≤ 𝑑1 + 𝑑2(
)𝑚 1 (
𝑟→0
𝜕𝑟
)
⎪ 1
𝑄
𝑑1 + 𝑑2 ≤ 𝑧 ≤ 𝑑⎩ 𝐾 2𝜋𝑟𝑑𝑟3

Applying the Laplace and Fourier cosine transforms to Eq. (A10), one
can obtain:
[(
)𝑚 (
)
𝑛𝜋𝑑1
𝑑 𝑠̂(𝑟, 𝑛, 𝑝)
𝑄
−𝑑
= lim
𝑠𝑖𝑛
𝑟→0
𝑟→0 𝑛𝜋𝑝𝐾𝑟
𝑑𝑟
2𝜋𝑟𝑑1
𝑑
) )]
(
)𝑚 ( (
𝑛𝜋 𝑑1 + 𝑑2
𝑄
+
𝑠𝑖𝑛
2𝜋𝑟𝑑3
𝑑
lim

(A11)

The modiﬁed Bessel function has the following properties(Spanier and Oldham, 1987):
𝑥𝑑 𝐾𝜈 (𝑥)+ 𝜐𝐾𝜈 (𝑥) = −𝑥𝐾𝜈−1 (𝑥)
𝑑𝑥

(A12)

and
𝐾 𝜈 (𝑥 ) = 𝐾 − 𝜈 (𝑥 )

(A13)

and
𝐾 𝜈 (𝑥 ) ≈Γ(𝜈) ( 𝑥 )−𝜈
, 𝑥 → 0, 𝜈 > 0
2
2

(A14)

From Eqs. (A9), (A12), (A13), and (A14), one has:( √ )− 23−𝑚()
𝜑𝑑 𝑠̂(𝑟, 𝑛, 𝑝)
1
1 √
2
lim
= − lim 𝐶2 𝑚 𝜑Γ𝑟→0
𝑑𝑟
2 𝑟→0 𝑟
3−𝑚
3−𝑚

(A15)

K. Tu, Q. Wu and J. Simunek et al.

Advances in Water Resources 145 (2020) 103740

The integration constantEqs. (A11) and (A15) as:2𝑑( √ ) 2 [(
3−𝑚

𝜑

𝑄
2𝜋𝑑1

3−𝑚

𝐶2 =

)𝑚

𝑠𝑖𝑛

C2

(

can

𝑛𝜋𝑑1
𝑑

)

(
+

then

𝑄
2𝜋𝑑3

)𝑚

be

𝑠𝑖𝑛

(

)
√ (
𝑛𝜋𝑝𝐾𝑟 𝜑Γ 3−2𝑚

obtained𝑛𝜋 (𝑑1 +𝑑2 )
𝑑

from

)]
(A16)

Substituting Eq. (A16) to Eq. (A9), one gets:2𝑑( √ ) 2 [(𝑠̂(𝑟, 𝑛, 𝑝) =

1−𝑚
2

×𝑟

3−𝑚

𝜑

3−𝑚

𝐾 1−𝑚
3−𝑚

𝑄
2𝜋𝑑1

( √
2 𝜑
3−𝑚

𝑟

)𝑚

𝑠𝑖𝑛

(

𝑛𝜋𝑑1
𝑑

)

(
+

𝑄
2𝜋𝑑3

)𝑚

𝑠𝑖𝑛

(

)
√ (
𝑛𝜋𝑝𝐾𝑟 𝜑Γ 3−2𝑚)𝑛𝜋 (𝑑1 +𝑑2 )
𝑑

)]

(A17)

When the inverse Fourier transform is employed to Eq. (17), one canobtain:𝑠̄(𝑟, 𝑧, 𝑝) =

(

1
𝑛𝜋𝑧
𝑠̂(𝑟, 0, 𝑝) +𝑠̂(𝑟, 𝑛, 𝑝)𝑐 𝑜𝑠
𝑑
𝑑 𝑛=1
𝑑

)

2
𝐶2 ′ =

3−𝑚

𝛼
3−𝑚

𝑄
2𝜋𝑑1

)𝑚

𝑑1 −

(

)
√ (
𝑝𝐾𝑟 𝛼Γ 3−2𝑚

𝑄
2𝜋𝑑3

)𝑚

𝑑3

( √ ) −2()
𝜔 3−𝑚𝑑 𝑠̂(𝑟, 𝑛)
1
1 √
2
lim
= − lim 𝐶4 𝑚 𝜔Γ𝑟→0
𝑑𝑟
2 𝑟→0 𝑟
3−𝑚
3−𝑚

2𝑑
𝐶4 =

𝑠̂(𝑟, 0, 𝑝) =

𝛼
3−𝑚

3−𝑚

𝑄
2𝜋𝑑1

)𝑚

𝑑1 −

(

)
√ (
𝑝𝐾𝑟 𝛼Γ 3−2𝑚

𝑄
2𝜋𝑑3

)𝑚

𝑑3

]
𝑟

1−𝑚
2

×𝑟

𝐾 1−𝑚
3−𝑚

( √
2 𝛼
3−𝑚

𝑟

3−𝑚
2

)

Finally, the analytical solution in the Laplace domain can be obtainedby substituting Eqs. (A17) and (A20) to Eq. (A18).
Appendix B. Derivation of the steady-state analytical solution
Similarly, the Fourier cosine transform is used to the second-orderpartial derivative with respect to coordinate z in Eq. (14) to get:
[ 2
]
(
)
𝑑 2
𝜕 𝑠(𝑟, 𝑧)𝜕 𝑠(𝑟, 𝑧)𝑛2 𝜋 2
𝑛𝜋𝑧
𝐹𝑐
=
𝑐 𝑜𝑠
𝑠̂(𝑟, 𝑛)(B1)
𝑑𝑧 = −
∫0
𝑑
𝜕 𝑧2
𝜕 𝑧2
𝑑2
where 𝑠̂ is drawdown in the Fourier domain.
By applying the Fourier cosine transform to Eq. (14) and combining
it with Eq. (B1), Eq. (14) can then be reduced to:
𝑑 2 𝑠̂(𝑟, 𝑛) 𝑚 𝑑 𝑠̂(𝑟, 𝑛)+− 𝜔𝑟1−𝑚 𝑠̂(𝑟, 𝑛) = 0
𝑑𝑟𝑟 𝑑𝑟

𝑚𝐾𝑧 𝑛2 𝜋 2 𝑄 𝑚−1() . The general solution of Eq. (B2) is ob𝐾𝑟 𝑑 2 2𝜋𝑑in which 𝜔 =tained as:𝑠̂(𝑟, 𝑛) = 𝑟

1−𝑚
2

(B2)

[
𝐶3 𝐼 1−𝑚
3−𝑚

( √
2 𝜔
3−𝑚

𝑟

3−𝑚
2

)
+ 𝐶4 𝐾 1−𝑚
3−𝑚

( √
2 𝜔
3−𝑚

𝑟

3−𝑚
2

)]
(B3)

where C3 and C4 are integration constants; 𝐼 1−𝑚 (⋅) and 𝐾 1−𝑚 (⋅) denote
3−𝑚

3−𝑚

the modiﬁed Bessel functions of the ﬁrst and second kind with the order𝑚of 1−, respectively. Using the boundary condition Eq. (3), the constant3−𝑚
C3 can be obtained by applying the Fourier transform and properties of

)𝑚

𝑠𝑖𝑛

(

𝑛𝜋𝑑1
𝑑

)

(
+

𝑄
2𝜋𝑑3

)𝑚

𝑠𝑖𝑛

(

)
√ (
𝑛𝜋𝐾𝑟 𝜔Γ 3−2𝑚𝑛𝜋 (𝑑1 +𝑑2 )
𝑑

)]
(B7)

𝜔
3−𝑚

1−𝑚
2

3−𝑚

𝐾 1−𝑚
3−𝑚

𝑄
2𝜋𝑑1

( √
2 𝜔
3−𝑚

)𝑚

𝑠𝑖𝑛

(

𝑛𝜋𝑑1
𝑑

)

(
+

𝑄
2𝜋𝑑3

)
√ (
𝑛𝜋𝐾𝑟 𝜔Γ 3−2𝑚
)

3−𝑚

𝑟 2

)𝑚

𝑠𝑖𝑛

(

𝑛𝜋 (𝑑1 +𝑑2 )
𝑑

)]

(B8)

Applying the inverse of Fourier cosine transform to Eq. (B8), one has:𝑠(𝑟, 𝑧) =
(A20)

𝑄
2𝜋𝑑1

( √ ) 2 [(𝑠̂(𝑟, 𝑛) =𝑟( √ ) 2 [(

3−𝑚

𝜔
3−𝑚

2𝑑

where 𝛼 = 𝑚𝑝𝑆( 2𝑄)𝑚−1 . Then, with Eq. (A19), one has:𝐾𝜋𝑑2( √ ) 2 [(Substituting Eq. (B7) to Eq. (B4), one can obtain:

]
(A19)

(B6)

Combining Eqs. (B5) and (B6), the constant C4 can be given as:(A18)

When n = 0, similar processes are applied to solve Eq. (A5) andboundary conditions, and ﬁnally the integration constant C2 ′ can be obtained as:( √ ) 2 [(Using the Fourier cosine transform to Eq. (6), one can obtain:
[(
)𝑚 (
)
𝑛𝜋𝑑1
𝑑 𝑠̂(𝑟, 𝑛)
𝑄
−𝑑
lim
= lim
𝑠𝑖𝑛
𝑟→0
𝑟→0 𝑛𝜋𝐾𝑟
𝑑𝑟
2𝜋𝑟𝑑1
𝑑
) )]
(
)𝑚 ( (
𝑛𝜋 𝑑1 + 𝑑2
𝑄
+
𝑠𝑖𝑛
(B5)
2𝜋𝑟𝑑3
𝑑
With Eq. (B4) and the properties of the modiﬁed Bessel function as described in Eqs. (A12), (A13), and (A14), one has:

3−𝑚
2

∞
2∑

the modiﬁed Bessel function, leading to C3 = 0. Then, Eq. (B3) can befurther written as:
( √
)
1−𝑚
2 𝜔 3−𝑚𝑠̂(𝑟, 𝑛) = 𝐶4 𝑟 2 𝐾 1−𝑚
𝑟 2
(B4)
3−𝑚
3−𝑚

∞
(
)
𝑛𝜋𝑧
1
2∑
𝑠̂(𝑟, 0) +𝑠̂(𝑟, 𝑛)cos
𝑑
𝑑 𝑛=1
𝑑

(B9)

When n = 0, Eq. (B2) with corresponding boundary conditions canbe solved using similar methods. Then, using properties of the modiﬁed Bessel function from Eq. (A12) to Eq. (A14), the term 𝑠̂(𝑟, 0) in
Eq. (B9) can be obtained as:
(
)𝑚
(
)𝑚
𝑄
𝑄
𝑑1 − 2𝜋𝑑𝑑3
2𝜋𝑑13𝑠̂(𝑟, 0) =
𝑟1−𝑚
(B10)
)
(
2
𝐾𝑟 Γ 3−𝑚 (3 − 𝑚)
Substituting Eqs. (B8) and (B10) to Eq. (B9), we can eventually obtain the steady-state analytical solution in the time domain.
References
Archer, G.E.B., Saltelli, A., Sobol, I.M., 1997. Sensitivity measures, ANOVA-like techniques
and the use of bootstrap. J. Statist. Comput. Simulation 58 (2), 99–120.
Barker, J.A., Herbert, R., 1977. Non-Darcy ﬂow and its implications to seepage problem.
J. Irrig. Drain. Eng.-ASCE 459–473.
Bear, J., 1972. Dynamics of Fluid in Porous Media. Elsevier Science, New York.
Bordier, Zimmer, 2000. Drainage equations and non-Darcian modelling in coarse porous
media or geosynthetic materials. Journal of Hydrology.
Brunetti, G., Šimůnek, J., Piro, P., 2016. A comprehensive numerical analysis of the hydraulic behavior of a permeable pavement. J. Hydrol. 540, 1146–1161.
Chen, Y.F., Hu, S.H., Hu, R., Zhou, C.B., 2015. Estimating hydraulic conductivity of fractured rocks from high-pressure packer tests with an Izbash’s law-based empiricalmodel. Water Resour. Res. 51 (4), 2096–2118.
Chen, Y.F., Liu, M.M., Hu, S.H., Zhou, C.B., 2015. Non-Darcy’s law-based analytical models
for data interpretation of high-pressure packer tests in fractured rocks. Eng. Geol. 199,
91–106.
Ciriello, V., Di Federico, V., 2012. Similarity solutions for ﬂow of non-Newtonian ﬂuids
in porous media revisited under parameter uncertainty. Adv. Water Res. 43, 38–51.
Crump, K.S., 1976. Numerical inversion of Laplace transforms using a Fourier-series approximation. J. ACM (JACM) 23 (1), 89–96.

K. Tu, Q. Wu and J. Simunek et al.
De Hoog, F.R., Knight, J.H., Stokes, A.N., 1982. An improved method for numerical inversion of Laplace transforms. SIAM J. Sci. Stat. Comput. 3 (3), 357–366.
Di Federico, V., Ciriello, V., 2012. Generalized solution for 1-D non-Newtonian ﬂow in a
porous domain due to an instantaneous mass injection. Transp. Porous Media 93 (1),63–77.
Escande, L., 1953. Experiments concerning the ﬁltration of water through rock mass. Proceeding of Minnesota International Hydraulics Convention.
Feng, Q., Wen, Z., 2016. Non-Darcian ﬂow to a partially penetrating well in a conﬁned
aquifer with a ﬁnite-thickness skin. Hydrogeol. J. 24 (5), 1287–1296.
Forchhelmer, P.H., 1901. Wasserbewegun durch boden. Zeltschrlft des Verelnes Deutscher
Ingenleure 49, 1736–1749.
Houben, G.J., 2015. Hydraulics of water wells-ﬂow laws and inﬂuence of geometry. Hydrogeol. J. 23 (8), 1633–1657.
Izbash, S.V., 1931. O ﬁltracii v kropnozernstom materiale. USSR, Leningrad in Russian.
Liu, M.M., Chen, Y.F., Hong, J.M., Zhou, C.B., 2016. A generalized non‐Darcian radial
ﬂow model for constant rate test. Water Resour. Res. 52 (12), 9325–9343.
Mathias, S.A., Butler, A.P., Zhan, H., 2008. Approximate solutions for Forchheimer ﬂow
to a well. J. Hydraul. Eng. 134 (9), 1318–1325.
Mathias, S.A., Todman, L.C., 2010. Step-drawdown tests and the Forchheimer equation.
Water Resour. Res. 46 (7).
Mathias, S.A., Butler, A.P., Zhan, H.B., 2008. Approximate solutions for Forchheimer ﬂow
to a well. J. Hydraul. Eng. 134 (9), 1318–1325.
Mathias, S.A., Wen, Z., 2015. Numerical simulation of Forchheimer ﬂow to a partiallypenetrating well with a mixed-type boundary condition. J. Hydrol. 524, 53–61.
Miller, R.J., Low, P.F., 1963. Threshold gradient for water ﬂow in clay systems. Soil Sci.
Soc. Am. J. 27 (6), 605–609.
Moutsopoulos, K.N., Papaspyros, I.N., Tsihrintzis, V.A., 2009. Experimental investigationof inertial ﬂow processes in porous media. J. Hydrol. 374, 242–254.
Moutsopoulos, K.N., Tsihrintzis, V.A., 2005. Approximate analytical solutions of the Forchheimer equation. J. Hydrol. 309 (1–4), 93–103.
Muskat, M., 1938. The ﬂow of homogeneous ﬂuids through porous media. Soil Sci. 46 (2),169.
Ni, L., Li, H., Jiang, Y., Yao, Y., Ma, Z., 2011. A model of groundwater seepage andheat transfer for single-well ground source heat pump systems. Appl. Therm. Eng.
31, 2622–2630.
Nossent, J., Elsen, P., Bauwens, W., 2011. Sobol’ sensitivity analysis of a complex environmental model. Environ. Model. Softw. 26 (12), 1515–1525.
Pianosi, F., Sarrazin, F., Wagener, T., 2015. A Matlab toolbox for global sensitivity analysis. Environ. Model. Softw. 70, 80–85.
Quinn, P.M., Cherry, J.A., Parker, B.L., 2011. Quantiﬁcation of non-Darcian ﬂow observed
during packer testing in fractured sedimentary rock. Water Resour. Res. 47 (9).
Quinn, P.M., Parker, B.L., Cherry, J.A., 2013. Validation of non-Darcian ﬂow eﬀects inslug tests conducted in fractured rock boreholes. J. Hydrol. 486, 505–518.
Rose, H.E., 1951. Fluid ﬂow through beds of granular material. In: Some Aspects of Fluid
Flow. Edward Arnold, London, pp. 136–163.
Rybach, L., 2015. Innovative energy-related use of shallow and deep groundwaters-examples from China and Switzerland. Cent. Eur. Geol. 58 (1–2), 100–113.
Saltelli, A., Chan, K., Scott, E.M., 2000. Sensitivity Analysis. Wiley, New York, USA.
Sedghi-Asl, M., Rahimi, H., Salehi, R., 2014. Non-Darcy ﬂow of water through a packed
column test. Transp. Porous Media 101 (2), 215–227.
Sen, Z., 1985. Volumetric approach to type curves in leaky aquifers. J. Hydraul. Eng. 111(3), 467–484.
Sen, Z., 1989. Nonlinear ﬂow toward wells. J. Hydraul. Eng. 115 (2), 193–209.
Şen, Z., 1990. Nonlinear radial ﬂow in conﬁned aquifers toward large‐diameter wells.
Water Resour. Res. 26 (5), 1103–1109.
Sen, Z., 2000. Non-Darcian groundwater ﬂow in leaky aquifers. Hydrol. Sci. J. 45 (4),595–606.

Advances in Water Resources 145 (2020) 103740
Slepicka, F., 1961. The laws of ﬁltration and limits of their validity. In: IAHR Proceedings
9th Convention, 38, pp. 383–394.
Sobol, I.M., 1993. Sensitivity estimates for nonlinear mathematical models. Math. Model.
Comput. Exp. 1 (4), 407–414.
Sobol, I.M., 2001. Global sensitivity indices for nonlinear mathematical models and their
Monte Carlo estimates. Math. Comput. Simul 55 (1–3), 271–280.
Soni, J.P., Islam, N., Basak, P., 1978. An experimental evaluation of non-Darcian ﬂow inporous media. J. Hydrol. 38 (3–4), 231–241.
Spanier, J., Oldham, K.B., 1987. An Atlas of Functions. Hemisphere, New York, USA.
Stehfest, H., 1970a. Algorithm 368: Numerical inversion of Laplace transforms [D5]. Commun. ACM 13 (1), 47–49.
Stehfest, H., 1970b. Remark on algorithm 368: numerical inversion of Laplace transforms[D5]. Commun. ACM 13 (10), 624–625.
Tu, K., Wu, Q., Simunek, J., Chen, C., Zhu, K., Zeng, Y., Xu, S., Wang, Y., 2020. An analytical solution of groundwater ﬂow in a conﬁned aquifer with a single well circulation
system. Water Resour. Res. 56 e2020WR027529.
Tu, K., Wu, Q., Sun, H., 2019. A mathematical model and thermal performance analysis
of single-well circulation (SWC) coupled ground source heat pump (GSHP) systems.
Appl. Therm. Eng. 147, 473–481.
Watanabe, H., 1982. Comment on Izbash’s equation. J. Hydrol. 58 (3–4), 389–397.
van Werkhoven, K., Wagener, T., Reed, P., Tang, Y., 2009. Sensitivity-guided reduction ofparametric dimensionality for multi-objective calibration of watershed models. Adv.
Water Res. 32 (8), 1154–1169.
Wan, J.W., Huang, K., Chen, C.X., 2013. Reassessing Darcy’Law on water ﬂow in porous
media. Earth Sci.-J. China Univ. Geosci. 38 (6), 1327–1330 in Chinese.
Wen, Z., Huang, G., Zhan, H., 2006. Non-Darcian ﬂow in a single conﬁned vertical fracture
toward a well. J. Hydrol. 330, 698–708.
Wen, Z., Huang, G., Zhan, H., 2008. An analytical solution for non-Darcian ﬂow in a
conﬁned aquifer using the power law function. Adv. Water Res. 31 (1), 44–55.
Wen, Z., Huang, G., Zhan, H., 2008b. Non-Darcian ﬂow to a well in an aquifer–aquitard
system. Adv. Water Res. 31 (12), 1754–1763.
Wen, Z., Huang, G., Zhan, H., 2009. A numerical solution for non-Darcian ﬂow to a wellin a conﬁned aquifer using the power law function. J. Hydrol. 364 (1–2), 99–106.
Wen, Z., Huang, G., Zhan, H., Li, J., 2008a. Two-region non-Darcian ﬂow toward a wellin a conﬁned aquifer. Adv. Water Res. 31 (5), 818–827.
Wen, Z., Liu, K., Chen, X., 2013. Approximate analytical solution for non-Darcian ﬂow
toward a partially penetrating well in a conﬁned aquifer. J. Hydrol. 498, 124–131.
Wen, Z., Liu, K., Zhan, H., 2014. Non-Darcian ﬂow toward a larger-diameter partiallypenetrating well in a conﬁned aquifer. Environ. Earth Sci. 72 (11), 4617–4625.
Wilkinson, J.K., 1955. Flow of water through rock ﬁll and its application to the design ofdams. N.Z. Eng. 10 (11), 382–387.
Wu, Q., Tu, K., Sun, H., Chen, C., 2019. Investigation on the sustainability and eﬃciencyof single-well circulation (SWC) groundwater heat pump systems. Renew. Energy 130,656–666.
Wu, Q., Xu, S.H., Zhou, W.F., LaMoreaux, J., 2015. Hydrogeology and design of groundwater heat pump systems. Environ. Earth Sci. 73 (7), 3683–3695.
Wu, Y.S., 2002. An approximate analytical solution for non‐Darcy ﬂow toward a well in
fractured media. Water Resour. Res. 38 (3) 5–1.
Yeh, H.D., Chang, Y.C., 2013. Recent advances in modeling of well hydraulics. Adv. Water
Res. 51, 27–51.
Zhou, J.Q., Hu, S.H., Fang, S., Chen, Y.F., Zhou, C.B., 2015. Nonlinear ﬂow behavior at low
Reynolds numbers through rough-walled fractures subjected to normal compressive
loading. Int. J. Rock Mech. Min. Sci. 80, 202–218.
Zoorabadi, M., Saydam, S., Timms, W., Hebblewhite, B., 2015. Nonlinear ﬂow behaviorof rough fractures having standard JRC proﬁles. Int. J. Rock Mech. Min. 76, 192–199.