Advances in Water Resources 164 (2022) 104182

Contents lists available at ScienceDirect

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

Seepage to staggered tunnels and subterranean cavities: Analytical and
HYDRUS modeling
A.R. Kacimov a, *, Yu V. Obnosov b, J. Šimůnek caDepartment of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Oman
N.I. Lobachevsky Institute of Mathematics and Mechanics, Kazan Federal University, Kazan, Russiac
Department of Environmental Sciences, University of California Riverside, CA, USAb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Subterranean holes
Complex potential
Darcian velocity
Isobars
Seepage force
HYDRUS modeling

Darcian, 2-D flows to subterranean holes are studied analytically (by the methods of complex analysis) andnumerically (by HYDRUS). For flow towards two circular or quasi-circular tunnels, reconstructed as isobars
generated by two sinks under a ponded homogeneous soil surface, the flow nets, the velocity vector fields, andRiesenkampf’s seepage force vectors are found. The position of the two tunnels is optimized using a criterion of
the total area of the empty space comprised by the isobars with the admissible seepage flow rate as a constraint
and the locus of the tunnels as a control variable. The case of a partially-filled tunnel, the contour of which iscomposed of two conterminous isobaric and equipotential arcs, is also analyzed. For a grouted tunnel, therefraction problem for two potential fields in two subdomains of the half-plane of the seepage flow domain issolved for an arbitrary contrast between the hydraulic conductivities of the liner and ambient soil. The tunnel
grouting is, generally, eccentric with respect to the tunnel contour. High hydraulic gradients in the grouting aredetected, which is a long-term seepage-harbinger for any liner materials used by geotechnical engineers.

(continued )List of abbreviations and notations
SymbolDescription

PK-62,77

Polubarinova-Kochina (1962). Theory of Ground
Water Movement.
van Genuchten
Ponded water depth
Volumetric soil water content
Residual water content
Saturated water content (porosity)
Saturated hydraulic conductivities
VG soil parameters
Pressure head
Total (piezometric) head =p+y in analytical and p+z in
HYDRUS models
Time
Cartesian coordinates in analytical solutions
Cartesian coordinates in HYDRUS
Specific discharge (Darcian velocity)
Intensities of sinks (seepage flow rates of tunnels)
Depth of the sink and centre of circular tunnel under
the ponded soil surface

VG
Н
θ
θr
θs
K, Kg
n, α
p
h
t
x,y
x,z
V
Q 1, Q 2
d, dc

Unit
m
[-]
[-]
[-]
m/day
[-], m− 1
m
m
day
m
m
m/day
m2/daym(continued on next column)

w=ϕ+iψ ,z=x+iyϕ, ψRComplex potential, complex physical plane varaibles inanalytical solutions
Velocity potential and stream function
Radius of circular tunnel

m2/day,mm2/daym

1. Introduction
Subterranean openings, as man-made structures, appeared at thedawn of human civilizations, viz.:• holes (quasi-cylindrical in shape and about 2 × 1 m in sizes) into
which coffins (Christian) or wrapped dead bodies (Islam) are placed
and backfilled, at a depth of 1− 2 m from the ground surface of
cemeteries;• metro, highway, sewage network tunnels at depths of 100− 101 m;• horizontal wells, in particular, qanat (falaj in Oman) tunnels tapping
unconfined alluvium aquifers, at depths of 101− 102 m;• Repositories of nuclear waste in deep rock at depths of 102− 103 m,* Corresponding author.
E-mail address: anvar@squ.edu.om (A.R. Kacimov).
https://doi.org/10.1016/j.advwatres.2022.104182
Received 6 December 2021; Received in revised form 7 March 2022; Accepted 26 March 2022
Available online 27 March 2022
0309-1708/© 2022 Elsevier Ltd. All rights reserved.

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182• civil defense bunkers of WWI-WWII in the urban infrastructure, at
depths of 100− 101 m;

He et al., 2021; Zhang et al., 2019, 2020b). Mathematically, boundary
value problems (BVPs) to the elliptic (e.g., Laplace’s, advective disper­
sion, and modified Helmholtz’s) or parabolic (e.g., nonlinear
Richards-Richardson and linear/nonlinear Boussinesq diffusion) PDEs
are formulated. The tunnel itself and the seepage-feeding boundary arecommonly assumed to be a horizontal equipotential (and therefore
isobaric) line or surface. A slightly curved water table with accretion iscommonly approximated by a horizontal isobar. The tunnel itself isassumed to be either a circular or quasi-circular isobar or equipotential
contour. When the seepage BVP is solved, i.e., the flow net, the scalar
field of pore pressure, and the vector-field of the Darcian velocity are
determined, the analysis of strains-stresses determines the stability ofthe tunnel’s slopes (the tunnel’s roof zone is most risky) and near-tunnel
skeleton for stability. All types of soils’ deformations, erosion, andcollapsing (suffusion, colmation, heaving, sliding, etc., see, e.g.,
Fernández and Alvarez, 1994; Ivars, 2006; Lee and Nam, 2001, 2004;
Lee et al., 2006; Pand and Dias, 2016; Verruit and Booker, 2000; Zhang
et al., 2014) should be avoided. If this geomechanical examination of thetunnel and its soil vicinity detects unstable conditions-regimes (seep­
age-induced deformations or soil’s degradation), then protection mea­sures (drainage and lining) are designed-constructed.
However, a rigorous mathematical analysis of seepage into clustersof tunnels and grouted tunnels is missing. Specifically, in this paper, we
present new analytical and numerical solutions to BVPs for seepage
flows towards:

among others. These openings (hereafter called tunnels) are of a regularshape and well-known loci-orientation with respect to the vector g ofgravitational acceleration.
Natural (uncontrolled) subsurface openings intersperse the soil androck continuum:• large gardens’ and crop fields’ areas have an intricate network ofholes made by subterranean mammals, e.g., shrewmice, pocket go­phers, moles, opossums, voles, and raccoons (10− 1-100 m);• holes made by social insects in their nests and adjacent soils (10− 310− 2 m);• wild desert trees and shrubs generate multiple holes made by plantsroots, which in arid deserts can span 60− 70 m deep;• catacombs of karst caves and cavities;
among others. These holes, unlike engineered entities, are irregular in
their location and orientation (see, e.g., photos in the chapter “Tunnel
erosion” in Zhu and Xu, 2021).
All types of the abovementioned openings (hereafter called tunnels)
are subject to seepage of ambient pore water, which commonly perco­lates down. The soil water interacting with the tunnels can originate
from rainfall, irrigation, and ponding of the ground surface by flash
floods, among others. Moreover, some tunnels are often purposely
constructed deep in the beds of shallow lakes, lagoons, or even seafloors(e.g., for masking the ICBM structure or storing nuclear wastes, e.g., theSwedish Final Repository near Forsmark, 60 m deep under the floor of
the Baltic Sea).
Seepage flows into tunnels (and in their vicinity) is a crucial phe­nomenon determining their very sustainability and operational effi­
ciency. For example, seepage into a tunnel of a London tube is a greatnuisance. Intensive seepage of the highest possible flow rate is the only
motivation to construct a qanat (falaj) hydraulic structure. Railroad,highway, urban, dam, and water transportation tunnels vary in theirtechnical specifications (see, e.g., Ashjari et al., 2019; British, 2014,
Chen et al., 2017, 2016; Hornero et al., 2021; Jin et al., 2016; Li et al.,2022, 2021; Preene and Chrimes, 2021; Wang et al., 2019; Wu et al.,2020; Xu et al., 2021; Xue et al., 2021) but most of them are subject tothe inrush of the ambient groundwater with ensued hazards to thestructures and ambient environment-ecosystems (e.g., subsidence of theground surface).
Hydrologists, hydrogeologists, soil physicists, and civil (geotech­nical, groundwater) and military engineers studied (numerically andanalytically) 2D, 3D, steady-state or transient seepage into a solitarytunnel or a cluster of 2-3 tunnels, either grouted or not, placed in ho­
mogeneous or heterogeneous, isotropic or anisotropic soils, which wereposited to be porous (double-porous or fractured) media (continua). The
skeleton was assumed to be either rigid or deformable and groundwater
either compressible or incompressible (see, e.g., Aghdam et al., 2021;
An et al., 2021; Butscher, 2012; Dassargues et al., 2021; El Tani, 2003,
2010; Farhadian and Katibeh, 2017; Fernandez and Moon, 2010a, b;
Font-Capó et al., 2011; Fujii and Kacimov, 1998; Gattinoni and Scesi,
2010; Gokdemir et al., 2021, 2022; Gong et al., 2018; Guo et al., 2021;
Hassani et al., 2018; Houben et al., 2022; Huang et al., 2021; Huangfu
et al., 2010; Kacimov, 2005, 2007a, 2009; Kolymbas and Wagner, 2007;
Lei, 1999; Li et al., 2018; Liu and Li, 2021; Maleki, 2018; Mityushev andAdler, 2019; Park et al., 2008a, b; Perrochet, 2005a, b; Perrochet and
Dematteis, 2007; Qin et al., 2020; Shahbazi et al., 2021; Su et al., 2017;
Sweetenham et al., 2017; Tang et al., 2018a, b; Wang et al., 2008; Wu
et al., 2021; Yang and Yeh, 2007; Yang et al., 2016, 2021; Ying et al.,2018; Wang et al., 2014; Zhang and Franklin, 1993; Zhang et al., 2020a).
While the Darcy law for the saturated or unsaturated pore water
motion is postulated, the hydraulic gradients close to tunnels’ contours
are often so high that non-Darcian corrections are necessary (see, e.g.,a) two empty tunnels,b) a tunnel partially filled with water,c) a grouted tunnel.
We answer the following main questions:• How much pore water seeps into the tunnels?
• Are the soils and liners around the tunnel able to withstand the porepressures, their gradients, hydraulic gradients, and seepage forces?
Moreover, from the ontological azimuth, we believe that mostmodern modelers in subsurface mechanics (except Strack, 2017, and his
group and a few others, like Tong et al., 2022) have abandoned thepowerful methods of the theory of holomorphic functions, developed,
mastered, and applied to the problems of agricultural horizontal
drainage by Vedernikov (1939) and Polubarinova-Kochina (1962, 1977)(hereafter abbreviated as PK-62,77). Indeed, most of theabove-reviewed analytical solutions to BVPs are obtained by themethods of images and Fourier series expansions applied to a harmonic
function (piezometric head). Our analytical and HYDRUS solutions
presented below give the spectrum of seepage characteristics simulta­
neously: the flow net (piezometric head and stream function), porepressure distributions, vector fields of the Darcian velocity (hydraulicgradients), and the Riesenkmapf’s seepage force function (see PK-77).
Thus, in this paper, we expand and – despite our focus on the subsur­
face (Gaya-Pluto’s realms) - astronomize the analysis of seepage totunnels (including ones, still invisible from the space, on Mars, butmathematically acting as equipotential and isobars (Kacimov et al.,2022a), engaging two synergetic modeling weapons: the theory of hol­omorphic functions and HYDRUS.
2. Analytical solution for empty and partially-filled tunnels
We introduce the gauge pressure of pore water, P, P= γ w p, where p is
the pressure head, γ w =ρ g, ρ is a constant water density, and g is themagnitude of gravity acceleration. The Richards-Richardson equation(see, e.g., Radcliffe and Šimůnek, 2010) governs the dynamics of porewater in both saturated and unsaturated forms:

2

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 1. Vertical cross-sections. Two empty circular tunnels, placed one above another, under a ponded soil surface (a). An empty tunnel under accreted phreatic
surface (b).

∂θ= ∇(K(p)∇h).
∂t

Zingg and Anagnostou, 2018), are ignored in our analysis.

(1)

In Eq. (1),∇ is the nabla operator in Cartesian coordinates (Oxy) (seeFig. 1a, which illustrates a vertical cross-section), θ(t, x, y) is a volu­
metric water content, k(p) is the hydraulic conductivity function, and h(t, x, y) = p + y is the total (piezometric) head. In fully saturated ho­
mogeneous isotropic soils, k(p) = K = const (the saturated hydraulicconductivity). Then for steady-state or transient flows without freeboundaries (e.g., phreatic surfaces), Eq. (1) is reduced to the Laplace
equationΔϕ(t, x, y) = Δψ (t, x, y) = Δp(t, x, y) = 0,

2.1. Two empty staggered tunnels with their centers-connecting segment
collimated with gravity
In this Subsection, we study analytically the case of two empty cir­cular tunnels staggered in such a manner that the line connecting theircenters, C1 (the upper tunnel) and C2 (the lower tunnel), is parallel withthe vector of gravity, g (see similar settings in Guo et al., 2021; Huang
et al., 2020). The emptiness of the tunnels is engineeringly ensured byperfect drainage of water dripping (exfiltrating) into them.
The distance between the tunnels’ centers is a, their radii are R1 andR2, and the depths of the centers, counted from the horizontal soil sur­face Ox, are d1 and d2. Water seeps in the lower half-plane of Fig. 1, withtwo circles cut of this half-plane, i.e., the flow domain Gz is mathe­
matically triply-connected, bounded by two circles and a straight line Ox(a degenerated circle).
In Fig. 1a, we show the case of a soil surface, ponded by a layer ofstatic water of a depth H. Obviously, p = H and φ = − KH along Ox, and p= 0 along the tunnels’ circumferences. We recall that vacuum drainage
has been tested in clayey soils (see, e.g., Grigor’ev, 1973). In vacuumed
tunnels, the boundary condition along the tunnels’ contours should be p=-pv, where pv is a positive constant determined by the vacuum-creating
air pump. Vice versa, if flooding-ponding of soil is fast, then the pore air
can be piston-displaced by imbibing water into the tunnel, and there­
fore, the excess positive pressure p = pv can be maintained in a“dead-end” subterranean opening like a dead root hole in Fig. 1b.
We follow PK-62 and select point I (infinity) as fiducial, i.e., assume
that at this point ψ = 0, w = − KH. On the Riemann sphere, the infinity
is imaged by a single point I, albeit this point triplicates in the physicalplane of Fig. 1a and our HYDRUS models below. We also follow PK-62
and Vedernikov (1939) and assume that the pore water is stagnant at→
infinity, i.e., V (∞) = 0, which hydrogeologically means that soil issubtended by an impermeable horizontal bedrock (no-flow boundary).

(2)

where Δ is the Laplacian operator, (φ(t, x, y), ψ(t, x, y)) is the dyad of→
harmonic functions, viz. ϕ(t,x,y) = − Kh = ∇ V is the specific dischargepotential (also called “velocity potential” by PK), h(t, x, y) = p + y is the→total (piezometric) head, V (t, x, y) is the vector of specific discharge,and ψ is the stream function, conjugated with φ via the Cauchy-Riemannconditions (see, e.g., PK-62, 77, and Strack, 2017). We introduce acomplex physical coordinate z = x+i y, complex potential w = φ + iψ,and complexified Darcian velocity V=u +i v, where u(t,x,y) and v(t,x,y)
are the horizontal and vertical components of the specific discharge
vector. The function w(z) is holomorphic, and V(z) is antiholomorphic(PK-62). We also introduce the Riesenkampf force per unit volume of asaturated porous medium (see PK-77 for more details), which is the→→vector FR (t, x, y) = γw I − (1 − θs )(γs − γ w ) j = ∇Ξ where Ξ(t, x, y) = − p− γ0y, γ0 = θsγw + (1 − θs)γs, and γw,γs, and θs are the specific weights ofwater, solid particles, which make the porous medium, and porosity,
respectively. The force potential, Ξ(t, x, y), is also a harmonic function.
This force was evaluated in Kacimov and Obnosov (2016) for problemsof subsurface irrigation where seepage from buried emitters can cause
the heaving of the soil surface above them. In tunnels, which are underconstruction, 3-D effects of seepage near an advancing (not lined and not
grouted) tip of the tunnel (see, e.g., Chen et al., 2022; Hong et al., 2010;3

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 2. A vertical cross-section with two empty tunnels (the left panel) and a complex potential domain (the right panel). The physical flow domain, Gz, is the fourth
quadrant of the z-plane. Seepage is symmetric with respect to the Oz axis. In Fig. 2 (right panel), the complex potential domain, Gw (a half-strip with a horizontal cut
E2SE1) corresponds to Gz. A dashed horizontal segment SrS in Gw corresponds to a curved branch of the water divide (separatrix) SSr in Gz. The solid curves A1B1 andA2B2 in Gw are the images of tunnels’ right half-contours in Gz.

In other words, if there are no tunnels, then there is no seepage. A deep
drainage layer or a deep water table can be alternative infinity condi­
tions (Kacimov et al., 2022b, PK-62).
Instead of a ponded soil surface, we can also consider the case oftunnels placed under a phreatic surface, which is formed by recharge ofintensity e from the vadose zone (Fig. 1b, where only one tunnel isdepicted) through which the rain or irrigation water infiltrates. In thecase of a deep tunnel and free deep drainage to infinity (the “Free
Drainage” boundary condition in HYDRUS, see Simunek et al., 2018), a
steady-state phreatic line (a dashed curve in Fig. 1b) is only mildlydepressed above a seeping tunnel. Hence this curve can be approximatedby a straight horizontal line, and mathematically, we arrive at the case offlow in Fig. 1a (see, e.g., Butscher, 2012). If e and the tunnel sizes are
small, then unsaturated flow takes place near what in soil physics iscommonly interpreted as a subsurface cavity, which can act as animpermeable barrier (Philip, 1998; Kacimov, 2007a; Kacimov et al.,2022b). Yucca Mountain deep subsurface structures in Nevada (where eis small) were also investigated as potentially impermeable subterra­nean openings (see, e.g., Philip, 1998).
The topology of seepage flow in Fig. 1a is characterized by a fourbranched separatrix and a stagnation point S. The two separatrix
branches, the curved streamlines SlS and SrS, start from the soil surface
and end at point S (located between the tunnels). The remaining
branches (straight streamlines) SB1 and SA2 start at S and end at thebottom point B1 of the upper tunnel and apex A2 of the lower tunnel,
respectively. The segment SlSr of the ponded surface feeds the upper
tunnel with the flow rate Q1. The lower tunnel is fed by the rays ISl and
ISr with the flow rate Q2. Shi et al. (2020) studied topologically similar

seepage regimes to vertical wells, i.e., ignoring gravity and zero-pressure
boundaries.
For a single tunnel (drain), Vedernikov (1939) and later PK-62 used
the method of images and showed that a strictly circular draincompletely filled with water (i.e., a constant piezometric head drain) canbe modeled by a line sink placed at a certain depth d. Such a sink gen­erates a family of equipotential lines-circles, each centered at a certain
depth larger than d. Two families of streamlines (semi-circles) areorthogonal to the equipotential circles. Fujii and Kacimov (1998)(hereafter FAK-98) used the Vedernikov solution and showed that theisobar p = 0 generated by a line sink, albeit not exactly circular, can beconsidered as a “quasi-circle” of a vertical diameter dv and a horizontal
diameter dh. Fig. 2 sketches these isobars as solid-line ovals. Strictly
circular tunnels’ contours (modeled in Section 4 below) deviate fromthese ovals in the same sense as the Vedernikov (1939) (see also PK-62,77) equipotentials generated by line sinks deviate from strict circles,
which physically model the lines of contact between gravel packs ofhorizontal drains and the ambient saturated soil. If deviations of isobars
in our paper or Vedernikov’s equipotentials from circles is insignificant,then empty tunnels or drains can be approximated by line sinks. In otherwords, if – at a given sink’s depth d – its strength Q is not too high, thenthe isobar is almost circular, i.e., dv ≈ dh ≈ 2R, and the depth of the sinkis almost the same as the depth of the center of the circular tunnel.
FAK-98 also optimized the placement of an empty quasi-circular tunnel
under the ponded surface, viz. they proved that an empty tunnel of agiven size (diameter or cross-sectional area) seeps at a minimal flow rate
Qm if placed at a depth dm. Such an optimum is intuitively clear andfollows from the Rolle theorem. Indeed, if in Fig. 1, for a given solitary4A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 3. Isobars p = 0 and 0.4, and streamlines ψ = 0, ±0.5, ±1, ±1.5, ±1.8, ±2 for d1 = 0.7, d2 = 1.2, q1 = Q1/2 = 1, and q2 = Q2/2 = 1.25 (a); isobars p = 0,±0.5for d1 = 0.7, d2 =1.2, q1 = Q1/2 = 1, and q2 = Q2/2 = 1.9988 (b); vicinity of the isobar p = 0, flow net in the right half of Fig. 3b (c).

circular tunnel d1 → R1, then obviously Q1 → ∞, while at d1 → ∞, theflow rate also blows up, Q2 → ∞.
The method of superposition of singularities (real and fictitious) for
modeling of saturated and unsaturated flows is implemented (seeAppendix A), similarly to Kacimov (2007b), Kacimov and Al-Maktoumi(2019), Kacimov and Obnosov (2016, 2018), and Kacimov et al. (2018a,b, 2019, 2020). Since p = − (φ/K + y), obtaining the contours of tunnels’

isobars from Eq. (A2) in Appendix A is reduced to the solution of thefollowing implicit equation:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
Q1x2 + (y − d1 )2 Q2x2 + (y − d2 )2F(x, y) =+ ln− H + y = 0.
(3)
ln
2
2
2π
2
π
x + (y + d1 )x2 + (y + d2 )2
We introduce dimensionless quantities:()/
h* , p* , x* , y* , d1* , d2* , dS* , a* , dv*1 , dv*2 , dh*1 , dh*2 = (h, p, x, y, d1 , d2 , dS , a, dv1 , dv2 , dh1 , dh2 ) H,
)
/
/
(
/
/
Ar* = Ar H 2 , Q*1 , Q*2 , ϕ* , ψ * , w* = (Q1 , Q2 , ϕ, ψ , w) (KH), V * = V K, Fr* = Fr γ w

5

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 4. Vertical size dv1 of the upper tunnel made by the sink as a function of d2 for q1=1, d1=1, and q2=2 (curve 1) and q2=1 (curves 2) (a). Vertical size dv1 of theupper tunnel made by the sink as a function of d2 for q1 = 1, d1 = 1, and q2 = 2 (curve 1) and q2 = 1 (curves 2) (a), and for q1 = 1, d2 = 5, and q2 = 2 (curve 1) and q2= 1 (curve 2) (b). Vertical size of the lower tunnel dv2 as a function of the depth of the generating sink d2 for q1 = 1, d1 =1, and q2 = 2 (curve 1) and q2 = 1.5 (curve 2)(c), and as a function of the depth of the upper generating sink d1 for q1 = 1, d2 = 2, and q2 = 2.2 (curve 1) and q2 = 2 (curve 2) (d).

and regular increments of 0.26 and 0.35 for the equipotential lines andstreamlines, which are shown by dashed and solid lines, respectively. At
the point In indicated in Fig. 3c, the streamline ψ = -2.1 and equipo­
tential line ϕ =1.04 cross the tunnel contour. At this point, it is evident
that the streamline is not orthogonal to the line p = 0.
With a further increase of q2 (keeping other parameters in thisexample constant), one can not get two disconnected isobaric contours p= 0, i.e., two staggered tunnels in Fig. 1a can not even be approximatedby two mathematical sinks. Still, the equation p(0,x) is solvable but with
only two roots. If q > q2cr, a single p = 0 subsurface opening of a fancy
dumbbell shape can be modeled by two sinks. Two or more log singu­
larities could be used to model an elongated isobaric shape that mightlook like the cross-section of a typical qanat tunnel (Ebrahimi et al.,2021), similarly to what is done by the superposition of a line source anda line sink when impermeable “Rankine bodies” are obtained asstreamlines-separatrices (see, e.g., Kacimov and Obnosov, 1997; War­rick and Fennemore, 1995).
Similarly to variation of q2 up to its upper bound q2c, we varied q1
and d1 in the range 0 < d1 < d1cr and d2 in the range d2cr < d2 < ∞. These
results are skipped over for the sake of brevity.
We used the FindRoot routine of Mathematica to solve Eq. (3), in
which we put x = 0, with respect to y. We find yA1, yB1, yA2, yB2 andvertical “diameters” of the tunnels, dv1 and dv2. Fig. 4a shows dv1(d2) for
q1 = 1, d1 = 1, and q2 = 2 (curve 1) and q2 = 1 (curve 2). Fig. 4b shows
dv1(d1) for q1 = 1, d2 = 5, and q2 = 2 (curve 1) and q2 = 1 (curve 2). In
Fig. 4c, the curves dv2(d2) are plotted for q1 = 1, d1 = 1, and q2 = 2(curve 1) and q2 = 1.5 (curve 2). In Fig. 4d, the curves dv2(d1) are plottedfor q1 = 1, d2 = 2, and q2 = 2.2 (curve 1) and q2 = 2 (curve 2).
The horizontal sizes of the tunnels are determined slightly differently
from FAK-98. Specifically, we have to find both the ordinates andabscissae dh1/2 and dh2/2 of M1 and M2, i.e., the tunnel’s contour points
most “protruding” to the right in Fig. 2. For this purpose, we used
Mathematica to evaluate the function F1(x,y) = dx/dy along the isobars p= const. We skip over the lengthy explicit expression for F1. For the

Fig. 5. The horizontal sizes of both tunnels, dh1 and dh2, as functions of thedepth of the upper sink (curves 1 and 2, respectively) for q1 = 1, q2 = 1, and d2= 2.

and, for the sake of brevity, drop “*”s from them.
We used the routine ContourPlot of Wolfram’s Mathematica (1991)
and in Fig. 3a plotted the dimensionless isobars p = 0, and 0.4, for d1 =1,
d2 = 2, q1 = Q1/2 = 1, q2 = Q2/2 = 2, as well as the streamlines1. The
two contours p = 0, are almost circular and, therefore, can serve as an
excellent approximation of circular tunnels in Fig. 1a. Obviously, thestreamlines in Fig. 3a are not orthogonal to isobars. The segment A1O(between the surface and top of the upper tunnel) is a branch cut.
Fig. 3b shows the isobar p = 0 and p = ±0.5 for the same d1, d2, q1 as
in Fig. 3a but for a relatively large value of q2cr = 1.9988. In this criticalcase, the isobar p = 0 is made of two “coalescing” contours. In Fig. 3c,the vicinity of the two “osculating” components of this tunnel is zoomed.
The flow net is shown for the right half of the flow domain, with equal1
All our Mathematica programs are open codes and can be shared with thereferees/readers.

6A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 6. Area of both tunnels Ar(d1) for q1 = 1, d2 = 2, and q2 = 1, 0.8, and 0.6 (curves 1-3, respectively) (a), and for q1 = 1, a = 1, and q2 = 2 and 1 (curves 1-2,
respectively) (b). Magnitude of dimensionless Riesenkampf’s resultant force |Fr| along the ordinate axis for d1 = 1, d2 = 2, q1 = 1, q2 = 2, γs/γw = 2.65, and θs = 0.4,
segment OE1 (c), segment E1E2 (d), segment E1I (e).

isobar p = 0, the points M1 and M2 have the derivative dx/dy = 0.
Consequently, for a given d1, d2, q1, q2, we solved the system of two
equations:F(x, y) = 0,

F1 (x, y) = 0

d1max = 0.66 for the first curve.
There is a limitation for the quasi-elliptic approximation of tunnels’
shapes, which we used. The isobar p = 0, at the values of tetrad’s pa­
rameters close to the critical, can resemble axial sections of a teardrop(see Fig. 3b) or beetroot, which is weird in the application to tunnels.
The ordinate of the stagnation point yS =-dS is found as the solution(4)

with respect to x and y. The FindRoot routine of Mathematica solved thesystem (4) with an initial guess point close to (0, -d1) for the upper
tunnel (point M1) and another initial guess point close to (0, -d2) for the
lower tunnel (point M2). The roots xM1 and xM2, computed by Mathe­
matica, were eventually doubled to get dh1 and dh2 in Fig. 2.
Fig. 5 shows the horizontal sizes of both tunnels, dh1(d1) and dh2(d1)(curves 1 and 2, respectively) for q1 = 1, q2 = 1, and d2 = 2.
The dimensionless area, Ar = π(dv1dh1 + dv2dh2)/4, of both tunnels
was evaluated by approximating their shapes as ellipses having semiaxes (xM1, dv1/2) and (xM2, dv2/2). Fig. 6a illustrates the dependenceAr(d1) for q1 = 1, d2 = 2, and q2 = 1, 0.8, and 0.6 (curves 1-3, respec­tively). The curves have a global maximum, local minimum, and two
end extrema. For example, curve 2 has a global (end) minimum at d1 =0 and a global maximum at d1 = 0.425 that comports with the optimi­
zation in FAK-98. The local minimum and the local end maximum are,however, surprising. In Fig. 6b, Ar(d1) is plotted for q1 = 1and q2 = 2 and
1 (curves 1-2, correspondingly) with the distance a between tunnels’
centers fixed. Again, analogously to FAK-98, a single maximum isattained at a certain depth of the dyad of tunnels, e.g., Armax= 0.051 atof the equation V(0,y)=0. The complexified specific recharge is V(z) =dw/dz, where the overbar means complex conjugation (PK-62). Thenfrom Eq. (A1) we get:V(z) =Q1

d1

π z2 + d12

+

Q2

d2

π z2 + d22.

(5)

From Eq. (5) at x = 0, we get an equation with respect to y for
detection of point S:v(y) =

Q1

d1

π y2 − d12

+

Q2

d2

π y2 − d22

= 0.

(6)

Unlike the previous transcendent equations above, which we
tinkered by FindRoot, Eq. (6) is a quadratic equation that has two
explicitly found real roots. We selected a physically meaningful value(negative):√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
d1 d22 q1 + d2 d21 q2√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ .
yS = −(7)
d1 q1 + d2 q2
7

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 7. A partially filled solitary tunnel for d = 1 and q = 5, streamlines with an increment of 1 (labeled for the right half of the flow domain), isobar p = 0 (dashedcontour), and three equipotential lines φ = 0.7, 1.4, 1.6(a); right halves of compound tunnels, the upper part of which is made of a segment p = 0, and the bottom partis made of segments φ = 0.7 and 1.4 (b). Partially-filled tunnel AT1BT2A with the water level T1T2, seepage face T2AT1, and a constant piezometric level arc T1BT2
reconstructed for sink’s.

The trick now is to combine the isobar p = 0, generated by a single
sink of intensity Q = 2q and placed at the depth d, with certain equi­potential lines φ = φ1, φ2, φ3, φ4,.... These lines are circles centeredat points z = − de1,− de2,− de3,− de4,... and have the radii R1e, R2e, R3e, R4e,
… (see Vedernikov, 1939). Obviously, there is a double-inequality lim­
itation φm < φi < φM, i.e., not all equipotential circles − 1 < φi < ∞generated by a sink and determined by parameters Q, d intersect thecontour p = 0 generated by this sink.
Fig. 7a and b shows the case of q = 5 and d = 1. The dashed contour is
the isobar p = 0. A complete flow net is presented in Fig. 7a withstreamlines labeled for the right half of the flow domain. The isobar p=0
is plotted as a dashed line. Three equipotential lines ϕ=0.7, 1.4, and 1.6
are also plotted. Fig. 7a illustrates that the streamlines (Vedernikov’s,
1939 circles) are not orthogonal to the quasi-circular isobaric tunnel
contour, while they are orthogonal to the equipotentials (also Vederni­
kov’s circles).
In Fig. 7b, we sketched the right half of Fig. 7a. A segment of theequipotential line ϕ =0.7, together with the upper part of the isobar
p=0, makes a tunnel contour, which is filled with water up to the point
T1. This point is an intersection of the two curves, ϕ =0.7 and p=0. Thus,
the contour AT1 BT seeps as a line sink. The horizontal segment T1T2
represents the free surface of static water in the tunnel. Another partiallyfilled tunnel of the same q is made by a segment of the equipotential line
ϕ =1.4, which together with the upper part of the same isobar p=0
makes another tunnel contour. This tunnel is filled with water up to
point N1, which is an intersection of the two curves, ϕ =1.4 and p=0.
Thus, the contour AN1BN in Fig. 7b seeps as a line sink. The horizontal
segment N1N2 represents the free surface of static water in the second
tunnel.
The futile equipotential ϕ =1.6 is shown in Fig. 7a but not in Fig. 7b,
because this equipotential contour does not intersect with the isobar p =0, i.e. this equipotential line is physically useless. We used the Mathe­
matica routine FindRoot to solve a system of two equations: φ(x, y) = φi,p(x, y) = 0 with respect to x and y. For example, for φ = φN = 1.4, we got
xN1 = 0.36 and yN1 = -1.4. The maximum φM for this sink is 1.517 (theequipotential circle touches the isobar “internally”, i.e., there is almost
no static water in the tunnel) and the minimum φm for this sink is 0.419(the equipotential circle touches the isobar “externally”, i.e., the tunnel
is almost full of water).

Fig. 8. Vertical cross-section of a grouted tunnel.

For the parameters d1 = 1, d2 = 2, q1 = 1, and q2 = 2, this gives dS =1.27. Next, we put the root (7) back into Eq. (A2) where we also set x =0 and get the potential at point S. For the parameters in Fig. 3a, this gives
ϕS = 0.63 that - in the complex potential plane - determines the locus of
the tip of the cut Gw (Fig. 2).
The distribution of the magnitude of the dimensionless Rie­
senkampf’s resultant force |Fr| along the y-axis is plotted for d1 = 1, d2 =2, q1 = 1, q2 = 2, the ratios of specific gravities of water and solid
particles γs/ γ w = 2.65, and porosity θs = 0.4. Fig. 6c, d, and e show themodulus of the Riesenkampf force along the rays OE1, E1E2, and E2I,
respectively. The maximum force occurs along the perimeters of the twotunnels at four points, which are found by solving the equation p(0,y)=0.
The roots of these equations are yr =-0.97, -1.03, -1.93, -2.06, dependingon the initial guess point along the ordinate axis. The dashed rays
intersecting the curves in Fig. 6c–e are called the tunnel perimeters on
the graphs. For these six parameters, the force attains a minimum valueof zero at two points on the symmetry line: y = -1.152 and y = -2.624.
These minima are found by the FindMinimum routine of Mathematica.
2.2. Single partially-filled tunnel
In this Subsection, we consider a single tunnel partially filled withstatic-free water. We use Vedernikov’s (1939) solution, a special case of
Q1 = Q, Q2 = 0, d1 = d, see Eq. (A1). The source of this water inside thetunnel can be either seepage from a saturated ambient soil or convey­
ance of water from somewhere in the direction perpendicular to the Oxyplane.
8

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 9. The grouted tunnel with the centre of the liner circle d = 1.1 m and a radius R = 0.5 m. Sink’s depth dc =1 m and sink’s strength Q = 5 m2/day in (a) and (b).
The streamlines ψ1 and ψ2 = 0, ±0.5, ±1, ±1.5, ±2, ±2.5 m2/day for conductivity ratios kr = 0.3 (a) and 0.1 (b). Sink’s depth d =1.25 m and d = 0.6 (kr = 0.1) in (c)and (d). Isotachs |V1,2| = const [m/day], for R = 0.5 m, d = 1.3 m, dc = 1 m, Q = 2 m2/day, and kr = 0.01 (e).

9

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

3. Analytical solution for saturated seepage to grouted tunnel

function everywhere except a simple pole at z = - i d. Therefore, its
complex conjugate function V2 (z) can be represented as follows:

In this Section, we present an explicit and rigorous analytical solu­
tion to 2-D steady saturated seepage towards a solitary grouted tunnel
placed under a ponded surface Ox (Fig. 8). The soil hydraulic conduc­tivity is now K, and the hydraulic conductivity of the grouting (alsomodeled as a Darcian porous medium) is Kg. These quantities are arbi­
trary, but unlike the problems of seepage to horizontal drains or verticalwells, for which gravel packs (filters) are made of materials with Kg>>K(see, e.g., Kacimov and Obnosov, 2018) in wells with skin-effect(Obnosov et al., 2010) and for grouted or lined tunnels, an opposite
inequality holds, K>>Kg, i.e., seepage into the tunnel is inhibited.
The two porous media are separated by a circumference of a radius R(in dimensional quantities), centered at point z = -i dc. An empty tunnel
is modeled by an isobar p = 0 (Fig. 8). This contour is generated by a line
sink of intensity Q = 2q. The sink is placed at a depth d (eccentric with
the interface between the two media). The seepage domain Gz is made oftwo subdomains: S1 = {z : |z + idc | > R, Imz < 0}, and S2 = {z: |z + idc|
< R}. Obviously, 0 < dc − R < d < dc + R.
The solution trick is the same as in Section 2: we reconstruct anisobar p = 0 and, if it is almost circular, then this line is interpreted as acontour of an empty tunnel. By varying d, we get an almost annular or
eccentric grouting zone of conductivity Kg, while the interior of thedashed contour in Fig. 8 is discarded. Similarly, Zhukovsky placed avortex in the interior of his airfoils (Zhukovsky, 1949) and discarded a
subdomain of his flow domain, which is demarcated by a separatrix
obtained by superposition of the vortex and dipole at infinity (ambient
unidirectional flow of an ideal fluid).
In the notations of this Section, the solution is written in terms of twocomplex potentials, and Darcy’s laws with two contrasting hydraulic→conductivities: w1 = φ1 + iψ1, φ1 = − Kh1, V 1 = ∇ϕ1 and w2 = φ2 + iψ2,→
φ2 = − Kgh2, V 2 = ∇ϕ2 , p1 = − φ1/K − y, p2 = − φ2/Kg − y, defined in theexterior and interior of the circle in Gz of Fig. 8a. For further compari­sons with HYDRUS below, we will keep dimensional quantities. Alongthe interface between S1 and S2, the refraction conditions hold://
ϕ1 K = ϕ2 Kg , ψ 1 = ψ 2 .
(8)

V2 (z) =− ∞ < x < ∞,(10)

where V20 (z) is holomorphic in S2. The characteristic function V20 iswritten as:) ∞(Q 1 − Λ2 ∑
− QΛ
Λn
.
(11)
V20 (z) =+2π
2π (z − ia− 1 )
z
−
ian
n=0
The following notations are introduced in Eq. (11):/()/Λ = (kr − 1) (kr + 1), an = T u2n a , T(z) = s(1 + z) (1 − z),//√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
s = d2c − R2 , a = (d − s) (d + s), u2 = (dc − s) (dc + s),(12)where kr=Kg/K is a dimensionless conductivity ratio. Points ± is are
symmetrical to each other about the real axis and the circle ∂S1; an → swhen n tends to infinity. The characteristic function V1 in S1 is:()∞
Q(1 − Λ) ∑
1
1
.
(13)
V1 (z) =Λn−2πz − ian z + ian
n=0
The complex potentials in S1 and S2 are obtained by integration of
Eqs. (10), (11), and (13) with respect to z. Integration of Eq. (13) is done
from z = 0 to z. Point z = 0 is fiducial, i.e., at this point w1(0) = − 1. Thenwe get:w1 (z) =∞Q(1 − Λ) ∑ian − zΛn log− 1.
2π
ian + zn=0(14)

In Eq. (14), the branches of logarithmic functions are fixed in the zplane with a branch cut from − ian to ian through infinity; the branches
vanish at the point z = 0.
Integration of Eqs. (10) and (11) is done from the point Ag in Fig. 8,
viz. from z = − i(dc − R) to z. At point Ag, the refraction conditions (seePK-62) have to be met: Imw1 ( − i(dc − R)) = Imw2 ( − i(dc − R)) = 0,
kr Rew1 ( − i(dc − R)) = Rew2 ( − i(dc − R)). Then we obtain:
(
(
)
)
Q
z + idQΛz − ia− 1log
log
w2 (z) = −−+i(R − dc − a− 1 )
2πi(d + R − dc )
2π
(15)
) ∞
(
(
)
Q 1 − Λ2 ∑z − ianΛn log+ kr w1 ( − i(dc − R)).
2πi(R − dc − an )n=0

The derivatives of the complex potentials are V1 = dw1 /dz and V2 =dw2 /dz. Functions V1 and V2 should obey the conditions of the abscissa
axis acting as an equipotential line and static pore water at infinity:
ImV1 (x) = 0,

− Q
+ V20 (z),2π(z + id)

V1 (∞) = 0.

The refraction boundary condition, written in a complex form,should be met along the circular interface between S1 and S2, viz.:))
2Kg V1 (τ) = (K + Kg V2 (τ) + (Kg − K R2 (τ + idc )− 2 v2 (τ), |τ + idc | = R,In Eq. (15), the branches of all logarithmic functions vanish at the
point − i(dc − R). All branches are fixed in the z-plane with a cut fromthe points ian up to infinity, except for the first term on the RHS of Eq.
(15). For this logarithm, the branch cut extends from the point z = − iddown to infinity.
We note that analytical solutions for seepage to an empty grouted
tunnel reported and used by Fernandez and Moon (2010a), Li et al.
(2018), Qin et al. (2020), Wang et al. (2008), Wu et al. (2021), Yanget al. (2016), and Ying et al. (2018) do not take into account (and do noteven articulate) the phenomenon of refraction on the interface between
soil and grouting. Specifically, seepage to either a sink (andsink-generated zero pressure quasi-circular isobar in our solution above)or a circular isobar with an annular liner “sheath” is always 2-D, in both
S1 and S2 subdomains of the half-plane of Fig. 8. This is evident from
Eqs. (11) and (13) (or Eqs. (14) and (15)), which represent the super­
position of an infinite number of sinks and sources for each of whichflow is 2-D.
In Fig. 9, we show the grouting circumference of a radius R = 0.5 m(dimensional quantities are used for further comparisons with HYDRUS
results in Section 4). In Fig. 9 a and b, computations were carried out ford = 1.1 m, dc =1 m, and Q = 5 m2/day. The series are truncated by N =50 terms. The streamlines ψ1 and ψ2 = 0,
± 0.5,± 1,± 1.5,(9)
where τ is a complex variable along the circumference in Fig. 8, which we recall - is the line of mechanical contact between the soil andgrouting medium. The refraction condition in this form is common in thetheory of holomorphic functions (see, e.g., the review by Zverovich andLitvinchuk, 1968, and the book by Litvinchuk, 2000, p.135, Eq. (I)). Eq.
(9) incorporates two refraction conditions (8), usually formulated in realrather than complex variables, viz. the continuity of the pressure(piezometric) head and the normal component of the two flux vectors
across the interface (see, e.g., PK-62, and Strack, 2017).
Our problem is closely related to the refraction problem for two
arbitrary circles with conductivity contrasting with that of the plane. An
exact analytical solution of the latter problem was obtained by Obnosov(2011) for a general case of arbitrary singularities (not only
sinks-sources). We solved the problem (8) and (9) using the same
method as in Obnosov (2011). The details of derivations are presented inAppendix B.
In S2, the complexified Darcian velocity is an antiholomorphic
10

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 10. The transport domain with two tunnels (the upper one is half empty, the lower one is empty), and its discretization using triangular finite elements in
HYDRUS-2D (left), and the domain detail (right).

have to coincide up to the third digit. Obviously, the isotachs |V2| insidethe dashed contour in Fig. 9c have to be discarded in the context of thetunnel problem. The computed isotachs clearly show that thicker”
grouting of the tunnel in its “roof” zone helps in reducing the Darcian
velocities there, which, however, remain dangerously high in the vi­
cinity of point Bg (the tunnel floor where the grouting is “thin”). We
recall that PK-77 formulated a criterion-inequality of seepage-induced
instability in Darcian flows: the magnitude of the hydraulic gradientin an unconsolidated soil obeys |V1 /K| > 1 (in our dimensionlessquantities). If the gradient is less than 1 (e.g., between the contour andthe grouting-soil interface in Fig. 9e where 1 > |V1 /K | = 0.8), then thesoil may still be unstable to internal erosion (suffusion).
In the degenerate case of kr = 1 (no grouting, Section 2), for the set ofparameters in Fig. 9c, i.e., a single tunnel generated by a sink at d = 1.3
m and Q = 2 m2/day, we have a much smaller size Dv = 0.0038 m, i.e.,
grouting does comprise a large tunnel’s area of the same Q. On the other
hand, if we fix Dv = 0.92 m, then a tunnel without grouting will have amuch larger Q = 8.2 m2/day, as compared with the grouted case in
Fig. 9c, i.e., kr = 0.01 is an efficient counter-seepage measure. It isnoteworthy that - unlike the case of “formation damage” and skin effectin aged wells (see, e.g., Houben and Treskatis, 2007), which is detri­mental to wells’ productivity, suffusion, and translocation of finest
particulates of a porous medium (lessivage in the vernacular of soil
pedology) towards the grouting circle in Figs. 8 and 9 is a positive
phenomenon. Indeed, if the near-grouting soil gets extra clogging
accrued by a long-lasting seepage, this will only bolster the main missionof the anti-seepage aim of grouting. By shifting the sink up and downwith respect to the center of the grouting circle in Fig. 8 (which can also
be an engineering design parameter), i.e., by the effects of eccentricboundaries, we can fortify the tunnel from seepage where it is critical.
Optimization of the liner thickness and of the locus of the tunnel itself
with purely seepage-related criteria for a rigid porous skeleton of S1 andS2 can be a preamble of a geomechanical optimization (see, e.g., Li andZhou, 2015) that can be carried out similarly to Section 2 (we skip overthis analysis for the sake of brevity).

± 2,± 2.5 m2/day (thin lines) are found by the separation of theimaginary part in Eqs. (14) and (15). The tunnel contour p = 0 (labeled
dashed line) is determined using φ1 and φ2 obtained by the separation of
the real part in Eqs. (14) and (15). Panels a) and b) in Fig. 9a and b
illustrate the cases of kr = 0.3 and 0.1, respectively. Obviously, thestreamlines inside the dashed curve, i.e., in the air zone of the tunnel, areunphysical.
The reconstructed “quasi-circular” tunnel contour has the verticalsize, Dv, which is evaluated as the difference in the ordinates of theuppermost point At on the isobar p = 0 and its lowermost point Bt. We
used again the Mathematica routine FindRoot to solve the equation p2(0,y)=0 and determined the loci of these points. For panels a) and b) in
Fig. 9, the roots, i.e., the y-coordinates of At and Bt, are -0.689 m and-1.41 m for a) and -0.577 m and -1.72 m for b). Consequently, the cor­
responding tunnels’ “diameters” Dv are 0.72 m and 0.87 m.
One has to be cautious with the above-described procedure ofreconstruction of the p = 0 isobar. Specifically, if Q in the examples ofFig. 9 increases above a certain critical value, then a negative pressurezone emerges in the soil beneath point Bg in Fig. 8. In this case, the wholesolution (14) and (15), albeit mathematically (formally) correct, can not
be applied to the problem of seepage to a grouted tunnel, since the basic
assumption in the conceptual model of this Section is that the soil andgrouting seep under non-negative pore pressures.
To illustrate this, we placed the sink in Fig. 9c and d as: d =1.25 andd = 0.6, respectively. In Fig. 9c, the isobar p = 0 is slightly below point
Bg. In Fig. 9d, with a sink “hoisted” with respect to the center of thegrouting circle, the isobar p = 0 is “bulging” above point Ag in Fig. 8, i.e.,
protruding into the aquifer. Both cases in Fig. 9c and d can be considered
as a “failed” grouting, albeit mathematically dc and d in Fig. 8, and ouranalytical solution can be arbitrary.
Fig. 9e presents Mathematica-plotted (with the help of the Con­
tourPlot routine) isotachs |V1,2| = const, for d = 1.3 m, dc = 1 m, Q = 2
m2/day, and kr = 0.01. The circular interface between the soil andgrouting is a circle of a radius R = 0.5 m. The grouting-air interface isplotted as a dashed quasi-circle, which has Dv = 0.92 m.
In this case, the contrast in hydraulic conductivities of two porous
subdomains of Gz is so high that N = 50 is not enough in Eqs. (10), (11),and (13). We increased summation in the series up to N = 150 terms toachieve the following convergence criterion: at points Ag and Bg, thevelocities (V1 and V2) from the soil and grouting sides of the interface

4. HYDRUS modelingIn this Section, we use the finite element model (FEM) HYDRUS-2D(Radcliffe and Šimůnek, 2010; Šimůnek et al., 2018) with the following
intentions:
11A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 11. Streamlines (a) in the transport domain displayed in Fig. 10, and the flow velocity field (b) and velocity vectors (c) in a zoomed domain around thetwo tunnels.

4.1. Tunnels with partial filling and their center-connecting line
collimated with gravity
In this Subsection, seepage is symmetric with respect to the line
connecting the tunnels’ centers, and therefore in HYDRUS, only the righthalf of the flow domain is modeled. In Fig. 10, we assume that the lower
tunnel is empty, while the upper one is partially filled with trickled
exfiltrated water, which accumulates at tunnel’s bottom (e.g., due topoor drainage).
We assume that both tunnels have a radius R = 1 m. The upper tunnel
is half-filled with static-free water, and C1M1 is a horizontal air-waterinterface. The tunnel centers are placed at z = -10 and -15 m. In HYD­
RUS, the flow domains are always finite, and the following rectangle
sizes in Fig. 10 are used: Lh = 20 m and Lv = 30 m. These values are large
enough (as compared with R) to mimic the infinity in Fig. 1a. Instead of asingle point I on the Riemann sphere, we have points I1, I2, and I3 in
Fig. 10. The segment OI1 is ponded, i.e., a constant pressure head p = 1 m
boundary condition is imposed there. There is no flow across the seg­ments I1I2, I2I3, I3B2, A2B1, and A1O. In the upper tunnel, p decreases
hydrostatically from ht = R/2 to 0 along the quarter circle B1M1. Hereand in Section 4.2, we consider the VG loam with (θr,θs,α, n, K) = (0.078,
0.43, 3.6 1/m, 1.56, 0.25 m/day).
A topological peculiarity of flow in this Subsection is the shape andposition of the separatrix. Unlike the case of two empty tunnels in Sec­tion 2, the stagnation point St is now located on B1M1, i.e., the constantpiezometric head segment of the bottom of tunnel 1. The branch SrSt1
divides all infiltrated water into two parts. Ponded water on the left ofpoint Sr seeps into the upper tunnel, while segment SrI1 feeds the lower
tunnel. The separatrix branch St1St2 demarcates a zone of seepage fromthe upper tunnel to the lower tunnel. At point St1, the separatrix kinks if
the flow is governed by the Laplace equation (Section 2). From thehodograph analysis, Kacimov and Youngs (2005) showed that the angle
between the separatrix branches at St1 is θ = 900, i.e., this streamline is
not orthogonal to the equipotential line (B1M1 in Fig. 10).
We skip over the details of the finite element meshing and the iter­
ative solution of the Richardson-Richards equation2. We solved theabove-stated boundary value problem (BVP) for this PDE with theFig. 12. Vertical cross-section with two arbitrary circular tunnels.

a) cross-compare the FEM results with those obtained analytically inthe previous Sections;b) incorporate various boundary conditions (e.g., tunnels partially fil­led with water rather than empty ones);c) model clusters of tunnels, centers of which are not on the samevertical line.
Generalizations b)-c) are not amenable to analytical considerations
for circular tunnels.
Without any loss of generality, we assume that the two tunnels are
perfectly circular in their vertical cross-sections (Fig. 1a). In all our
computations, we consider the van Genuchten’s (VG) soils (Šimůnek
et al., 2018), characterized by parameters θr,θs,α, n, K, where θris theirreducible volumetric water content, and α, n represent the VG shapeparameters. We changed the notation for the vertical Cartesian coordi­
nate to z, which is more common in soil hydrology.

2
All our HYDRUS-2D projects are publicly available and can be shared withthe referees/readers.

12A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 13. Seepage to two arbitrarily located empty tunnels in a VG loam for H = 1 m, Lh = 20 m, Lv = 30 m, R1 = 1.5 m, R2 = 2 m, C1 (7 m, -10 m), and C2 (13 m, -13m). Streamlines (a), Darican velocities (b).

following initial condition: p(z) decreasing hydrostatically from 31 m to
1 m when z increased from -30 m to 0. The steady-state flow velocity
field and streamlines are shown in Fig. 11.
Fig. 11a illustrates the HYDRUS streamlines, and Fig. 11b zooms thevicinity of the tunnels where the magnitude of Darcian velocity is colormapped. The point St1n with zero velocity is located between the FE node
having coordinates x = 0.59 m and z = -10.81 m and an outflow velocity
V = 0.0826 cm/d and the FE node having coordinates x = 0.81 m and z= -10.59 m and an inflow velocity V = 0.095 cm/d. Along the arc A1St1n,
the velocity decreases from the apex (tunnel’s roof) to the “quasi-stag­nation” point. Along St1nB1, the velocity increases again to 0.21 m/day.
For the lower tunnel, HYDRUS detects a red zone of high velocities near
the exfiltrating seepage face. Remarkably, the velocity magnitude in­creases from the apex A2, where V = 0.58 m/day, to B2, where V = 0.84
m/day. Despite the seepage plume from its ponded bottom, the upper
tunnel actually hydraulically “shades” its lower “brother,” which in­tercepts this plume. Indeed, HYDRUS shows that the bottom of thesecond tunnel is more susceptible to seepage erosion than the roof,
provided the criterion of seepage-induced instability is based on themagnitude of the hydraulic gradient, as PK-77 conceptualized in her
condition |∇h| > 1of instability of saturated porous massifs.
We used the Boundary Fluxes option of HYDRUS and retrieved theinfiltration and exfiltration fluxes as follows. From the ponded soil
surface, a standard infiltration flux qOI1 = 3.4285 m2/day. The exfil­
tration flux through A1M1St is qin1 = 0.836 m2/day. This quantity is asum of what comes through both the seepage face A1M1 and a constanttotal head arc M1St. The “hidden” infiltration flux through B1St to thelower tunnel is qout1 = 0.128 m2/day. The upper tunnel is hydrologicallya “flow-through” entity. The total flow rate into the second tunnel is asum of qout1 and qsin2 = 2.414 m2/day, which is the quantity seeping
through St2B2 directly from the ponded soil surface. The second tunnel isa “gaining” hydrological actor.
In HYDRUS, any degree of filling ht (not only R/2 as in the example
above) in any tunnel can be easily modeled. Moreover, at a specific value
ht in Fig. 10, the upper tunnel acts as a “passive” entity, for which

qin1=qout1, i.e., all water seeped into the tunnel from above is oozed back“downstream.”
For natural holes made by decayed roots of dead plants, such a“passive” regime is the only one possible because, unlike gainingloosing-flow through lacustrine hydrology entities (Ilyinsky and Kaci­mov, 1992; Kacimov, 2000, 2007b), from which water evaporates, ortunnels (usually equipped with pumps for dewatering of dripped water),macropores in soils do not have mechanisms to uptake the in-hole water(except macropores under ant domes where social insects consume
trickled water).
4.2. Arbitrarily staggered tunnels
In this Subsection, we consider the case of tunnels of radii R1 and R2
whose centers C1 and C2 are arbitrarily located at depths dc1 and dc2(Fig. 12). There is no symmetry in this problem, and we have to consider
the whole flow domain. Without any loss of generality with respect totunnels’ boundary conditions, we consider two empty unlined tunnels,
the centers of which are located distance j apart in the horizontal
direction.
The HYDRUS Cartesian coordinate system I4xz is selected such thatboth tunnels are located sufficiently far from the sides I1I2, I2I3, and I3I4of the simulation rectangle having sizes Lh and Lv and selected large
enough to avoid the “wall effect” (Fig. 12).
Fig. 13a presents the streamlines for seepage from I4I1 (ponded to thedepth of H =1 m) towards this tunnel couple in a VG loam for Lh = 20 m,
Lv = 30 m, R1 = 1.5 m, and R2 = 2 m, and the coordinates of tunnels’
centers (7 m, -10 m) and (13 m, -13 m). In Fig. 13b, the magnitude of the
Darcian velocity is color-plotted. Remarkably for both tunnels, themaximum Darcian velocity shifted from point A1 to the right and frompoint A2 to the left for the upper and lower tunnels, respectively. Thiscan be attributed to the mutual hydrodynamic shading of this tunnel
dyad. However, these maxima (about 0.35 m/day along both tunnels’
contours) are dangerously high because the gradients at this point are
0.35/0.25 >1, which is above the PK limit of erosional instability of13

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 14. Streamlines (top), pressure head distribution along OI3 (b), and flux distribution along OI1 (c) for d = -4 m, dc = -3.5 m, H = 1 m, Ri = 1 m, Re = 2, Kg = 0.1
m/day, and K = 1 m/day.

soils. A stagnation point S (numerically a zone rather than a point) is
indicated beneath the upper tunnel. The total flow rate into the tunnel is
Q1 + Q2 = 5.21 m2/day and Q2 = 3.46 m2/day.

-2.62 m, yB1 = -4.81 m and, therefore, Dv = 2.19 m, i.e., slightly higher
than that of the HYDRUS tunnel circle. The analytical velocity at theorigin of coordinates is vO = -0.23 m/day, while HYDRUS gives vO = -0.2
m/day. The analytical pressure heads at the interface points at theordinate axis are pA g= 2.16 m and pBg = 5.04 m, while HYDRUS results
are pAg = 2.03 m and pBg = 4.4 m. If we place the sink deeper in theanalytical model at d = 4.5 m, then (for the same Q) yA1 = -2.98 m, yB1 =-4.99 m, and therefore Dv is almost the same as the HYDRUS circle. The
value of vO is now exactly the same as one in HYDRUS, pAg = 2.22 m and
pBg = 4.76 m. Moreover, the maximum of velocity in Fig. 14c is the sameas in the analytical solution. Thus, HYDRUS results and this solution
match well. Consequently, HYDRUS can be used as a design tool when
the sizes of a circular tunnel and liner “sheath” are fixed, and one wants
to determine the seepage flow rate. The analytical solution can be usedin an “inverse” mode: if Q is fixed and the tunnel position with respect to
OI1 and the grouting circle are flexible (provided that a posteriori, thequasi-circular shape of the p = 0 isobar is checked, and seepage isverified to be at positive pressures everywhere in the porous “depressionzone” close to the tunnel bottom).

4.3. Grouted tunnels modeled by HYDRUS
In this Subsection, we model in HYDRUS a solitary empty tunnel of aradius Ri grouted by an eccentric annular liner, a “sheath” of a radius Re
such that the depths of the circles’ centers are d = -4 m and dc = -3.5 m,
respectively. The ponding depth H = 1 m.
The ambient soil and grouting material are made of hypotheticalsoils, in which the hydraulic conductivity, Kg = 0.1 m/day, which is 1/
10 of that of the soil, K = 1 m/day. We selected Lh = 7.5 m, Lv = 10 m, Ri= 1 m, and Re = 2 m, i.e., the thickness of the grouting is 1.5 m at theapex, while it is only 0.5 m at the bottom. The HYDRUS-simulated flowrate into the tunnel is Q/2 = 1.962 m2/day.
Fig. 14a presents the streamlines with a palette of the magnitude ofvelocity. Obviously, the grouting near the bottom of the tunnel is atextremely high hydraulic gradients.
A pressure head distribution along the vertical ordinate axis is shownin Fig. 14b. The grouting, especially the segment B1Bg (see also Fig. 8), is
subject to extremely high pressure head gradients. The |v(x)| plotted inFig. 14c is needed for controlling the horizontal size of the HYDRUS
domain, viz. Lh should be large enough to avoid the “wall effect” with thecriterion v → 0. at point I1.
For comparisons, we modeled the same tunnel and grouting using
the analytical formulae of Section 3, where we use Q = 3.924 m3/day
and the depth of the sink d = 4 m. A quasi-circular isobar p = 0 has yA1 =4.4. Phreatic flow with accretion
In this Subsection, we consider a grouted tunnel placed under a soil
surface with a given infiltration flux e, which is sufficiently small not togenerate surface ponding (Fig. 1b). A vadose zone is formed above thephreatic line p = 0. A phreatic fringe forms above this line and below theline p = -pca where pca is a positive constant, which characterizes thisfringe as a tension-saturated zone in the Vedernikov-Bouwer model (see14

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Fig. 15. HYDRUS-modeled seepage through an accreted water table and a grouted tunnel of Fig. 14. Isobaric contours sandwiched as -0.12< p < 0.2 for the
infiltration intensity e = 0.08 m/day from the soil surface (a); streamlines for e = 0.08 m/day (b), isobars for e = 0.06.

PK-62,77 and Kacimov et al., 2019). The same HYDRUS flow domain
and soils as in Fig. 14 are modeled, with the following modification of
the boundary condition: along the soil surface OI1, we now do not
specify a constant ponding depth but rather a constant infiltration rate e= μ K, where μ <1 is a constant determined by the rain intensity.
Obviously, μ should be small enough for such a condition to be physi­
cally realistic. If μ is above a certain threshold value, ponding of ahorizontal ground surface (or runoff) takes place.
Fig. 15a shows isobars as contour lines plotted in the range -0.12< p< 0.2 for infiltration factor μ = 0.08.The water table p = 0 is a “standard”
curve with an inflection point. The contour p=0.2 is a fancy “doubleogive”. The wiggling curve p = -0.12 can be considered as a kind of
capillary fringe boundary. The phreatic line (as well as the capillaryfringe) “caves” towards the tunnel, qualitatively similar to the shapes ofphreatic lines in seepage through homogeneous soils in Ilyinsky andKacimov (1992) and Kacimov and Obnosov (2021). In Fig. 15b, thestreamlines are shown for μ = 0.08. This Figure illustrates that imme­
diately under the soil surface, infiltration is vertical, but closer to thewater table, the streamlines deflect to the right. The shape of somestreamlines is quite whimsical, reflecting the impact of a draining tunnel
and a refracting effect of the grouting. In Fig. 15c, isobars are plotted for
μ = 0.06. In this case, the infiltration rate is relatively low such that thephreatic line does not “overhang” the tunnel. An unsaturated “gap” ispositioned above the tunnel apex. Above the water table, a purely un­
saturated (but 2-D) seepage occurs. If μ > 0.09, ponding of the soil
surface occurs in the right upper part of the HYDRUS domain.
If Kg is increased (keeping μ constant), then “caving” of the watertable becomes more pronounced such that it “cusps” towards the tunnel.

5. Concluding remarks
The case of a single quasi-circular tunnel modeled by a line sink,placed under a horizontal ponded soil surface, is amenable to analyticaltreatment both if the soil is homogeneous (Vedernikov, 1939) or
piece-wise homogeneous (Fujii and Kacimov, 1998, PK-77, and Section
3 above). The constraints and simplifications in the models used are: arigid porous matrix, incompressible, one-phase, and fully steady-state
saturated flows, which –despite potentially high hydraulic gradientsnear the modeled tunnels – are assumed to be Darcian. Transience can betaken into account via boundary conditions, provided no free (phreatic)
surfaces emerge.”
If the tunnel is strictly circular, then the theory of holomorphic
functions still gives an analytical solution for isobaric contours in ho­
mogeneous soils (Mityushev and Adler, 2019). Moreover, a tunnel of any(noncircular) shape can be conformally mapped onto a concentric cir­cular annulus (see the Bieberbach theorem, 1953), as Verruijt andBooker (2000) demonstrated for a circular tunnel. After that, theDirichlet problem for a harmonic velocity potential can be explicitly
solved by the Villa formula in terms of Schwarz’ integrals (see, e.g.,
Abubakirov and Aksent’ev, 2015; Akhiezer, 1970; Bieberbach, 1953).
Tong et al. (2022) demonstrated this procedure for the case of a rect­
angular tunnel. In the future, we plan to utilize this technique of solving
a mixed, rather than Dirichlet, boundary value problem for a strictly
circular but partially-filled tunnel.
For the case of two circular tunnels (Fig. 1a), the problem seems
prohibitively complex for analytical treatment. A triply-connected
seepage flow domain in Fig. 1a still can be conformally mapped onto
15

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

a canonical domain (Bieberbach, 1953), e.g., a half-plane with two slits(the so-called Zhukovsky drains, see PK-62,77). However, to the best of
our knowledge, there is no explicit integral formula, similar to the Villaone, which would solve the Dirichlet problem in any triply-connecteddomain. Consequently, using numerical models, such as HYDRUS, isthe only solid method to examine 2-D seepage towards two, three, andmore circular tunnels.
The results obtained above illustrate that HYDRUS and analyticalsolutions match well. The method of images (Vedenrikov, 1939, PK-62)
is used for circular interfaces between soil and air (a flat ponded groundsurface) and soil and grouting. These interfaces are considered as“mirrors” with respect to which lines-sinks (and more general hydro­
dynamic singularities) are placed in such a manner that the physical, e.
g., the Dirichlet, refraction, Robin, and other types of boundary condi­
tions are exactly satisfied. HYDRUS is versatile and can model any
domain boundaries and interfaces for steady or transient, saturated andunsaturated seepage in soils with an arbitrary type of layered hetero­geneities. The dyad of HYDRUS and the machinery of holomorphic
functions allows one to model (behold) the intricate details of fluxes,
pore pressures, and forces in the vicinity of invisible subterranean
entities.

acquisition, Investigation, Methodology, Project administration, Re­
sources, Software, Visualization, Writing – original draft, Writing – re­view & editing. Yu V. Obnosov: Conceptualization, Formal analysis,
Funding acquisition, Investigation, Methodology, Resources, Software,
Visualization, Writing – original draft, Writing – review & editing. J.
Šimůnek: Conceptualization, Formal analysis, Investigation, Method­
ology, Resources, Software, Visualization, Writing – original draft,
Writing – review & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments
This work was supported by the grants DR\RG\17 and IG/AGR/
SWAE/22/02, Sultan Qaboos University, Oman, and by the Kazan
Federal University Strategic Academic Leadership Program ("PRIORITY2030"). Helpful critiques and comments by two anonymous referees areappreciated.

CRediT authorship contribution statement
A.R. Kacimov: Conceptualization, Formal analysis, Funding

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.advwatres.2022.104182.
Appendix A
The method of superposition of line sinks is widely used in analytical studies of Darcian flows in well hydraulics (see, e.g., Vuković and Soro, 1997).
Since the seminal Zhukovsky (1949) ideas, potential flows of an ideal fluid around an airfoil or their cascades are modeled by combining two sin­gularities: a real dipole (which is equivalent to a unidirectional flow at infinity) and a fictitious vortex (placed inside the airfoil). The method simplifies
the analytical treatment of multiple connected domains by reducing them to a simple-connected one but with compensated singularities. Below, we
follow this Zhukovsky-laid avenue and extend the FAK-98 analysis to multiple (staggered) tunnels, starting from a dyad of line sinks shown in Fig. 2.
The flow domain Gz is now a simply-connected domain, viz. the half-plane z < 0. The characteristic function w(z) has two logarithmic singularities: oneat point E1 (z = -d1) and another at point E2 (z = -d2). The function V(z) is antiholomorphic in the half-plane, except points z = -d1, z = -d2, z = -dS, and z= ∞. Superposition of Vedernikov’s (1939) complex potentials for two-line sinks with two fictitious line sources yields:w=

Q1 z − id1 Q2 z − id2
ln
+ ln
− kH.
2π z + id1 2π z + id2(A1)

We select the branches of function (A1) as follows. Due to the symmetry of flow and the flow domain Gz, we consider only its right half (the rightbottom quarter-plane in Fig. 2). In further HYDRUS modeling, this will save computation time. The corresponding complex potential domain Gw is a
half-strip of a height (Q1+Q2)/2 with a cut E2SE1 (Fig. 2). Obviously, from Eq. (A1), ψ = 0 along the streamline IE2. The curved streamline SrS bi­
furcates at point S. Along the cut in Gw and along the straight segments SE1 and SE2 in Gz, we have ψ = − Q2/2. Along the straight segments OE1 in Gz,
the stream function is ψ = − (Q1 + Q2)/2. The value of potential K hS (where hS=const <H) at the stagnation point is a part of the solution.
The FAK-98 trick is: the isobars A1B1 and A2B2 in Gz are determined by solving the equation p(x,y) = 0. In Gw (Fig. 2), the images of these two
isobars are shown as two “thick” curves. The loci of points A1, B1, A2, and B2 are computed by solving the equation p(0,y) = 0.
We will show below that this equation can have four or two roots. In a special case of three roots, the contours of two conterminous empty tunnels“coalesce,” In the case of four roots, the corresponding tetrad (yA1, yB1, yA2, yB2) is used for obtaining the vertical “diameters” of the two tunnels as yA1 yB1 and yA2 - yB2. If the reconstructed curves A1B1 and A2B2 in Gz are close to half-circles, then the two generated isobars p = 0 approximate well thegenuinely-circular contours of the tunnels in Fig. 1a.
Thus, we separate the real part in Eq. (A1) to obtain the potential:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
Q1x2 + (y − d1 )2 Q2x2 + (y − d2 )2
ϕ=
+
− KH.
(A2)
ln
ln
2
2
2π
2π
x + (y + d1 )x2 + (y + d2 )2Similarly, we separate the imaginary part in Eq. (A1) to get the stream function.

16

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Appendix B
In this Appendix, we provide full details of the solution to the problem of seepage in a two-component composite made of an aquifer half-plane of
one conductivity and a circular inclusion of a contrasting conductivity where a line sink, which generates a quasi-circular isobar as a tunnel contour, isplaced. Specifically, we prove the validity of Eqs. (11) and (13) in Section 3.
The boundary conditions along the real axis and the circular interface between S1 and S2 are
ReV1 (x) = 0, − ∞ < x < ∞, V1 (∞) = 0,
V1 (t) = AV2 (t) + BR2 (t + idc )− 2 V2 (t), |t + idc | = R,

(B1)

here A = (K + Kg)/2Kg, B = (Kg − K)/2Kg, and there is a line sink of intensity Q = 2q at the point –id, i.e.
V2 (z) = −qπ(z + id)

(B2)

+ V20 (z),where the function V20(z)is holomorphic in the domain S2.
We consider the Mőbius transformations of the z-plane and of the ζ-plane onto each other by the functions
ζ = T − 1 (z) =z + is,z − isz = T(ζ) = is

ζ+1
,
ζ− 1

(B3)

where the points ± is are symmetric about both the real x-axis and circumference ∂S2, i.e.
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
s = d2c − R2 .

(B4)
− 1

The first function (B2) maps conformally the lower half of the z-plane onto the unit disk of the ζ-plane and T
introduce the notations
T − 1 ( − id) =

d− s
= a,
d+s

T(ζ) + idc = i(dc + s)

u2 =dc − s,dc + sT(ζ) + ia = i(d + s)

: Sk → Ωk,k = 1, 2 (see Fig. B1). We

ζ− a= T2 (ζ),
ζ− 1

(B5)

ζ − u2= i(dc + s)T1 (ζ).
ζ− 1We introduce new unknown functions Vk* (ζ) = Vk (T(ζ)), k =1,2 and omit the superscript * below. Then, using conditions (B1) and notations (B2)–(B5), we obtain
ReV1 (τ) = 0, |τ| = 1,
/ V1 (1) = 0,
V1 (τ) = AV2 (τ) − Bu2 T1 (τ)2 V2 (τ), |τ| = u,

(B6)

and
V2 (ζ) = −q

πT2 (ζ)(B7)+ V20 (ζ).

According to the first boundary condition (B6), the function V1(ζ) can be analytically continued into the ring u < |ζ| < 1/u via the symmetry

principle. Therefore, the continued function (we keep the same notation for it) satisfies the condition V1 (ζ) = − V1 (1/ζ) and can be represented as(B8)

V1 (ζ) = V1+ (ζ) − V1+ (1/ζ),
where the function V1+ (ζ)is holomorphic in the circle|ζ| < 1/u. Without any loss of generality, we can assume that V1+ (1) = 0.
Next, we consider the following function

Fig. B1. Physical (left panel) and reference (right panel) planes for the problem of seepage to grouted tunnels.
17

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

⎧
Bqu2
⎪
⎪
⎪ V1+ (ζ) − AV20 (ζ) −( / ), |ζ| < u,⎪⎨
πT1 (ζ)2 T2 u2 ζF(ζ) =
⎪
⎪
( / )
Aq
u2
⎪
⎪
V20 u2 ζ , |ζ| > u.
+ V1+ (1/ζ) − B⎩−
π T2 (ζ)
T1 (ζ)2(B9)

It satisfies the continuity condition F+(τ) = F− (τ) for |τ| = u according to the second equality (B6) and due to representations (B5) has a simple poleat the point ζ = u2. Additionally, the function F(ζ) is bounded at infinity and F(1) = 0. In comport with Liouville’s theorem, we have F(ζ) = C/T1(ζ),
where C is a constant to be determined. Clearly
C = limT1 (ζ)F(ζ) = −
ζ→1

Bqu2

πT1 (1)2 T2 (u2 )

=− i

Bq

π (dc + s).

Then from Eq. (B9) we obtain (simple algebra is omitted):
V1+ (ζ) − AV20 (ζ) −−Aq

πT2 (1/ζ)iBq(dc − d)(ζ − 1)
= 0,π(d + s)(aζ − 1)

+ V1+ (ζ) −

u2

B

(

2

2

)

|ζ| < u

V20 u ζ = i

T1 (1/ζ)

(B10)

/
Bq

, |ζ| < 1π(dc + s)T1 (1/ζ)

u.

Exclusion of the function V1+ (ζ) from the system (B10) leads to the relation()
Λu2 (ζ − 1)2 ( 2 )qΛ(dc − d)(ζ − 1)(ζ − 1)Λ(ζ − 1)
V20 (ζ) =−V20 u ζ − i
−
2
−
2
2
2
π (d + s)(dc − s)(au ζ − 1) (d + s)(aζ − 1) (dc + s)(u ζ − 1)(u ζ − 1)
Let
(B11)

V20 (ζ) = (ζ − 1)V21 (ζ).
Then we get the following functional equation with respect to V21(ζ)()
q
V21 (ζ) = Λu2 L1 (ζ)V21 u2 ζ − i
f0 (ζ),π(ζ − 1)(B12)

where,(B13)

fj (ζ) = c1 Mj− 1 (ζ) − c2 Mj (ζ) − c3 Lj+1 (ζ).
Consequently,
Lj (ζ) =ζ− 1,
u2j ζ − 1

Mj (ζ) =ζ− 1,
au2j ζ − 1c1 =Λ(dc − d),(d + s)(dc − s)

c2 =

1
,
d+s

c3 =

Λ
.
dc + sFrom the definitions (B13), (B14) follows:(())
Lj (ζ)L1 (u2j ζ )= Lj+1 (ζ), Lj (ζ)M( − 1 u)2j ζ = Mj− 1 (ζ),
2j
2j
Lj (ζ)M0 u ζ = Mj (ζ), Lj (ζ)f0 u ζ = fj (ζ).

(B14)

(B15)

We replace ζ with u2ζ in Eq. (B12). The result of this substitution we put into the RHS of Eq. (B12) and using the relations (B15), we arrive at:
)2
(
(
)
( 2
)
q
V21 (ζ) = Λu2 L2 (ζ)V21 u4 ζ − i
Λu L1 (ζ)f1 (ζ) + f0 (ζ) .
π(ζ − 1)
We replace ζ with u4ζ in Eq. (B12) and substitute the result into the RHS of the last equation. And so on. At the n-th step we obtain)n(()
V21 (ζ) = Λu2 Ln (ζ)V21 u2n ζ − i

n− 1∑( 2 )jq
Λu Lj (ζ)fj (ζ).
π(ζ − 1) j=0

In the limit, when n → ∞, taking into account Eq. (B11), the last relation gives the required solution in the form of a series uniformly and absolutely
convergent in the domain Ω2:
V20 (ζ) = − i

∞
(
q∑

π j=0

)j
Λu2 Lj (ζ)fj (ζ).

(B16)

We rewrite this solution as
(
)
∞
∑
)
( 2 )j ( 2qV20 (ζ) = − iΛuΛu c1 Lj+1 (ζ)Mj (ζ) − c2 Lj (ζ)Mj (ζ) − c3 Lj (ζ)Lj+1 (ζ) .
c1 M− 1 (ζ) +π

j=0Returning now to the z-plane we use relations (B3)–(B5), (B14), (B16) and obtain
(
)
− iΛ
− iΛ2
1
1
c1 M− 1 (ζ(z)) =, Λu2 c1 Lj+1 (ζ)Mj (ζ) = 2j−,z − ia− 1
z − aj z − bj+1
u
(
)
(
)
− i
1
1
− iΛ
1
1
c2 Lj (ζ)Mj (ζ) = 2j
, c3 Lj (ζ)Lj+1 (ζ) = 2j,−−z − bj z − bj+1
u z − aj z − bju18

Advances in Water Resources 164 (2022) 104182

A.R. Kacimov et al.

where
aj = s

1 + au2j,
1 − au2jThus,
V20 (z) = −bj = s

1 + u2j,
1 − u2ja0 = d,b0 = ∞.

( (
) (
)
(
)))
∞
∑
Λ
1
1
1
1
1
1
+
Λj Λ2
−
−
−
−
− Λ
=
z − iaj z − ibj+1
z − iaj z − ibj
z − ibj z − ibj+1
π z − ia− 1 j=0
(

q

(

−

) ∞
(
(
))
∞
∑
q 1 − Λ2 ∑
Λ
1+Λ
1
Λ
qΛ
Λj
− (1 − Λ)
Λj
−
+
=−
+
z
−
ia
z
−
ib
z
−
ib
π z − ia− 1
π
π
z
−
iaj
(z
−
ia
)
j
j
j+1
− 1
j=0
j=0
q

or due to Eq. (B2)V2 (z) = −) ∞(q 1 − Λ2 ∑
q
qΛ
Λj
.
−
+
π(z + id) π(z − ia− 1 )
π
z
−
iaj
j=0

From the first equality (B10), (B16) and (B8) one can obtain (we omit simple algebra)()∞
q(1 − Λ) ∑11
V1 (z) =Λj−.
z − iaj z + iajπ
j=0Thus, the validity of Eqs. (11) and (13) in Section 3 of the manuscript is proved.

References

Gokdemir, C., Li, Y., Rubin, Y., Li, X., 2022. Stochastic modeling of groundwater
drawdown response induced by tunnel drainage. Eng. Geol. 297, 106529 https://doi.org/10.1016/j.enggeo.2022.106529.
Gokdemir, C., Rubin, Y., Li, X., Xu, H., 2021. A vulnerability assessment method toevaluate the impact of tunnel drainage on terrestrial vegetation under various
atmospheric and climatic conditions. Adv. Water Res. 147, 103796.
Gong, B., Jiang, Y., Okatsu, K., Wu, X., Teduka, J., Aoki, K., 2018. The seepage control of
the tunnel excavated in high-pressure water condition using multiple times grouting
method. Processes 6 (9), 159.
Grigor’ev, V.M., 1973. Vacuum Water Drawdown. Moscow, Strojizdat (in Russian).
Григорьев. ВакууМное ВодопониЖение, Москва, Стройиздат. В. М. 1973.
Guo, Y., Wang, H., Jiang, M., 2021. Efficient iterative analytical model for undergroundseepage around multiple tunnels in semi-infinite saturated media. J. Eng. Mech. 147(11), 04021101. ASCE.
Hassani, A.N., Farhadian, H., Katibeh, H., 2018. A comparative study on evaluation ofsteady-state groundwater inflow into a circular shallow tunnel. Tunnel.
Underground Space Technol. 73, 15–25.
He, J.G., Pan, Z.Q., Chen, Y.F., Hu, R., Yang, Z., 2021. Control of non-Darcian flow byconsolidation grouting in the surrounding rocks of a concrete-lined pressure tunnel.
In: IOP Conference Series: Earth and Environmental Science, 861. IOP Publishing,
072081.
Hong, E.S., Park, E.S., Shin, H.S., Kim, H.M., 2010. Effect of a front high hydraulicconductivity zone on hydrological behavior of subsea tunnels. KSCE J. Civil Eng. 14(5), 699–707.
Hornero, J., Manzano, M., Custodio, E., 2021. Deciphering the origin of groundwater
inflow into the Talave tunnel (SE Spain). Sci. Total Environ. 789, 147904.
Houben, G., Treskatis, C., 2007. Water Well Rehabilitation and Reconstruction. McGrawHill, New York.
Houben, G.J., Collins, S., Bakker, M., Daffner, A., Triuller, F., Kacimov, A., 2022. Review:horizontal, directionally drilled and radial collector wells. Hydrogeol. J. 30,
329–357. https://doi.org/10.1007/s10040-021-02425-w.
Huang, F., Wang, D., Feng, Y., Zhang, M., 2020. Prediction of the collapse region induced
by a concealed karst cave above a deep highway tunnel. Adv. Civil Eng. 1–14.
https://doi.org/10.1155/2020/8825262 article 8825262.
Huang, Z., Zhao, K., Li, X., Zhong, W., Wu, Y., 2021. Numerical characterization ofgroundwater flow and fracture-induced water inrush in tunnels. Tunnel.
Underground Space Technol. 116, 104119.
Huangfu, M., Wang, M.S., Tan, Z.S., Wang, X.Y., 2010. Analytical solutions for steadyseepage into an underwater circular tunnel. Tunnel. Underground Space Technol. 25(4), 391–396.
Ilyinsky, N.B., Kacimov, A.R., 1992. Analytical estimation of ground-water flow aroundcutoff walls and into interceptor trenches. Ground Water 30, 901–907.
Ivars, D.M., 2006. Water inflow into excavations in fractured rock—a three-dimensionalhydro-mechanical numerical study. Int. J. Rock Mech. Mining Sci. 43 (5), 705–725.
Jin, X., Li, Y., Luo, Y., Liu, H., 2016. Prediction of city tunnel water inflow and its
influence on overlain lakes in karst valley. Environ. Earth Sci. 75 (16), 1–15.
Kacimov, A., Obnosov, Y.V., 2016. Tension-saturated and unsaturated flows from line
sources in subsurface irrigation: Riesenkampf’s and Philip’s solutions revisited.
Water Resour. Res. 52, 1866–1880.
Kacimov, A., Obnosov, Y.V., 2018. Well bore clogging of a pumping well in hydraulicinterference with an ambient groundwater flow: the Strack-Kostitsina refraction
problem in an annular composite redux. Hydrogeol. J 26 (7), 2531–2541.

Abubakirov, N.R., Aksent’ev, L.A., 2015. Explicit forms for the Schwarz integrals in anannulus and their applications. Siberian Math. J. 56 (1), 9–20.
Aghdam, A.T., Salmasi, F., Abraham, J., Arvanaghi, H., 2021. Effect of drain pipes on
uplift force and exit hydraulic gradient and the design of gravity dams using thefinite element method. Geotech. Geol. Eng. 39 (5), 3383–3399.
Akhiezer, N.I., 1970. Elements of the Theory of Elliptic Functions. Moscow, Nauka (inRussian). English Translation: 1990. American Mathematical Society.
An, Y., Li, J., Yue, J., Zhou, J., 2021. Effect of seepage force on tunnel face stability using
limit analysis with SRM. Geotech. Geol. Eng. 39 (3), 1743–1751.
Ashjari, J., Soltani, F., Rezai, M., 2019. Prediction of groundwater seepage caused byunclogging of fractures and grout curtain dimensions changes via numerical doubleporosity model in the Karun IV River Basin (Iran). Environ. Earth Sci. 78 (3), 1–17.
Bieberbach, L., 1953. Conformal Mappings. Chelsea, New York.
British Tunnelling Society, 2004. Tunnel Lining Design Guide. Thomas Telford, London.
Butscher, C., 2012. Steady-state groundwater inflow into a circular tunnel. Tunnel.
Underground Space Technol. 32, 158–167.
Chen, J., Min, F., Wang, S., 2022. Construction Technology of Large Diameter
Underwater Shield Tunnel. Springer Nature.
Chen, K.L., Wu, H.N., Cheng, W.C., Zhang, Z., Chen, J., 2017. Geological characteristics
of strata in Chongqing, China, and mitigation of the environmental impacts oftunneling-induced geo-hazards. Environ. Earth Sci. 76 (1), 1–16.
Chen, Y.F., Hong, J.M., Zheng, H.K., Li, Y., Hu, R., Zhou, C.B., 2016. Evaluation ofgroundwater leakage into a drainage tunnel in Jinping-I arch dam foundation insouthwestern China: a case study. Rock Mech. Rock Eng. 49 (3), 961–979.
Dassargues, A., Batelaan, O., Anceau, A., 2021. The first potentiometric map. Ground
Water 59 (5), 772–779.
Ebrahimi, A., Mehraban, Y., Omidvarborna, H., Vakilinejad, A., Al-Sayigh, A.R.S., 2021.
Kariz (ancient aqueduct) system: a review on geoengineering and environmental
studies. Environ. Earth Sci. 80 (6), 1–13.
El Tani, M., 2003. Circular tunnel in a semi-infinite aquifer. Tunnel. Underground SpaceTechnol. 18 (1), 49–55.
El Tani, M., 2010. Helmholtz evolution of a semi-infinite aquifer drained by a circulartunnel. Tunnel. Underground Space Technol. 25 (1), 54–62.
Farhadian, H., Katibeh, H., 2017. New empirical model to evaluate groundwater flow
into circular tunnel using multiple regression analysis. Int. J. Mining Sci. Technol. 27(3), 415–421.
Fernández, G., Alvarez Jr., T.A., 1994. Seepage-induced effective stresses and waterpressures around pressure tunnels. J. Geotech. Eng. ASCE 120 (1), 108–128.
Fernandez, G., Moon, J., 2010a. Excavation-induced hydraulic conductivity reduction
around a tunnel–Part 1: guideline for estimate of ground water inflow rate. Tunnel.
Underground Space Technol. 25 (5), 560–566.
Fernandez, G., Moon, J., 2010b. Excavation-induced hydraulic conductivity reduction
around a tunnel–Part 2: verification of proposed method using numerical modeling.
Tunnel. Underground Space Technol. 25 (5), 567–574.
Font-Capó, J., Vázquez-Suñé, E., Carrera, J., Martí, D., Carbonell, R., Pérez-Estaun, A.,
2011. Groundwater inflow prediction in urban tunneling with a tunnel boringmachine (TBM). Eng. Geol. 121 (1-2), 46–54.
Fujii, N., Kacimov, A.R., 1998. Analytically computed rates of seepage flow into drains
and cavities. Int. J. Numer. Anal. Methods Geomech. 22, 277–301.
Gattinoni, P., Scesi, L., 2010. An empirical equation for tunnel inflow assessment:
application to sedimentary rock masses. Hydrogeol. J. 18 (8), 1797–1810.

19

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Kacimov, A., Obnosov, Y.V., 2021. Infiltration-induced phreatic surface flow to periodic
drains: Vedernikov-Engelund- Vasil’ev’s legacy revisited. Appl. Math. Modell. 91,
989–1003.
Kacimov, A., Obnosov, Y.V., Al-Maktoumi, A, 2018a. Dipolic flows relevant to aquifer
storage and recovery: Strack’s sink solution revisited. Transp. Porous Media 123 (1),
21–44.
Kacimov, A., Obnosov, Y.V., Or, D., 2019. Evaporation induced capillary siphoning
through hydraulically connected porous domains: the Vedernikov-Bouwer model
revisited. Transp. Porous Media 129 (1), 231–251.
Kacimov, A., Obnosov, Y.V., Simunek, J., 2018b. Steady flow from an array of subsurface
emitters: Kornev’s irrigation technology and Kidder’s free boundary problems
revisited. Transp. Porous Media 121 (3), 643–664.
Kacimov, A., Obnosov, Y.V., Šimůnek, J., 2020. Seepage to ditches and topographic
depressions in saturated and unsaturated soils. Adv. Water Res. 145, 103732 doi.
org/10.1016/j.advwatres.2020.103732.
Kacimov, A., Obnosov, Y.V., Šimůnek, J., 2022a. Seepage-evaporation controlled
depletion of initially water-filled reservoirs on Earth and Mars: analytic versus
HYDRUS modeling. Icarus 372, 114719.
Kacimov, A., Obnosov, Y.V., Šimůnek, J., 2022b. Subterranean holes in the arid vadosezone as perturbations of descending pore water fluxes: J. R. Philip’s legacy revisited.
In: Abstracts of the International Conference “Water Resources Management andSustainability: Solutions for Arid Regions”, Dubai, March, 22-24, 2022.
Kacimov, A.R., 2000. Three-dimensional groundwater flow to a lake: an explicit
analytical solution. J. Hydrol. 240, 80–89.
Kacimov, A.R., 2005. Mathematical modeling of subsurface flow towards qanat tunnels.
In: Proc. International Conference on Qanat. Kerman, Iran. Iranian Academic Center
for Education, Culture and Research, pp. 239–251. November 23-24, 2005.
Kacimov, A.R., 2007a. Dipole-generated unsaturated flow in Gardner soils. Vadose Zone
J. 6, 168–174.
Kacimov, A.R., 2007b. Three-dimensional groundwater flow to a shallow pond: an
explicit solution. J. Hydrol. 337, 200–206.
Kacimov, A.R., 2009. Minimal-seepage depth of isobaric cavity under ponded conditions.
J. Irrig. Drain. Eng. 135 (1), 108–110. ASCE.
Kacimov, A.R., Al-Maktoumi, A., 2019. Riesenkampf’s vortex solution revisited for 2-D
commingling of groundwater in a three-layered aquifer: vertical-inclined-horizontal
seepage in aquitard. Adv. Water Res. 123, 84–95.
Kacimov, A.R., Obnosov, Y.V., 1997. Analytical solutions to a problem of sink-source
flow in a porous medium. Arab Gulf J. Sci. Res. 15 (2), 325–351.
Kacimov, A.R., Youngs, E.G., 2005. Steady-state water-table depressions caused byevaporation in lands overlying a water-bearing substratum. J. Hydrol. Eng. 10 (4),295–301. ASCE.
Kolymbas, D., Wagner, P., 2007. Groundwater ingress to tunnels–the exact analyticalsolution. Tunnel. Underground Space Technol. 22 (1), 23–27.
Lee, I.M., Nam, S.W., 2001. The study of seepage forces acting on the tunnel lining andtunnel face in shallow tunnels. Tunnel. Underground Space Technol. 16 (1), 31–40.
Lee, I.M., Nam, S.W., 2004. Effect of tunnel advance rate on seepage forces acting on theunderwater tunnel face. Tunnel. Underground Space Technol. 19 (3), 273–281.
Lee, S.W., Jung, J.W., Nam, S.W., Lee, I.M., 2006. The influence of seepage forces onground reaction curve of circular opening. Tunnel. Underground Space Technol. 22(1), 28–38.
Lei, S., 1999. An analytical solution for steady flow into a tunnel. Groundwater 37 (1),23–26.
Li, G., Wang, C., Wang, G., Xiao, Z., Wu, X., Jiang, F., 2021. Effect of the blocking waterand limiting discharge and surrounding rock permeability on the stability of subsea
tunnel. Geotech. Geol. Eng. 39 (2), 1365–1380.
Li, P., Wang, F., Long, Y., Zhao, X., 2018. Investigation of steady water inflow into a
subsea grouted tunnel. Tunnel. Underground Space Technol. 80, 92–102.
Li, P.F., Zhou, X.J., 2015. Mechanical behavior and shape optimization of lining structure
for subsea tunnel excavated in weathered slot. China Ocean Eng. 29 (6), 875–890.
Li, Z., Chen, Z.Q., He, C., Ma, C.C., Duan, C.R., 2022. Seepage field distribution and water
inflow laws of tunnels in water-rich regions. J. Mountain Sci. 19 (2), 591–605.
Litvinchuk, G.S., 2000. Solvability Theory of Boundary Value Problems and Singular
Integral Equations with Shift. Springer Science & Business Media.
Liu, J., Li, X., 2021. Analytical solution for estimating groundwater inflow into lined
tunnels considering waterproofing and drainage systems. Bull. Eng. Geol. Environ.
80 (9), 6827–6839.
Maleki, M.R., 2018. Groundwater seepage rate (GSR); a new method for prediction ofgroundwater inflow into jointed rock tunnels. Tunnel. Underground Space Technol.
71, 505–517.
Mityushev, V.V., Adler, P.M., 2019. Exact steady-state solution to air pumping fromunderground partly covered by an impermeable tarp. Acta Mech. 230 (3),1129–1144.
Obnosov, Y., 2011. Three–phase eccentric annulus subjected to a potential field induced
by arbitrary singularities. Q. Appl. Math. 69 (4), 771–786.
Obnosov, Y.V., Kasimova, R.G., Al-Maktoumi, A., Kacimov, A.R., 2010. Can
heterogeneity of the near-wellbore rock cause extrema of the Darcian fluid inflow
rate from the formation (the Polubarinova-Kochina problem revisited)? Comput.
Geosci. 36, 1252–1260.
Pan, Q., Dias, D., 2016. The effect of pore water pressure on tunnel face stability. Int. J.
Numer. Anal. Methods Geomech. 40 (15), 2123–2136.
Park, K.H., Lee, J.G., Owatsiriwong, A., 2008a. Seepage force in a drained circulartunnel: an analytical approach. Can. Geotech. J. 45 (3), 432–436.

Park, K.H., Owatsiriwong, A., Lee, J.G., 2008b. Analytical solution for steady-state
groundwater inflow into a drained circular tunnel in a semi-infinite aquifer: a revisit.
Tunnel. Underground Space Technol. 23 (2), 206–209.
Perrochet, P., 2005a. A simple solution to tunnel or well discharge under constantdrawdown. Hydrogeol. J. 13 (5-6), 886–888.
Perrochet, P., 2005b. Confined flow into a tunnel during progressive drilling: ananalytical solution. Groundwater 43 (6), 943–946.
Perrochet, P., Dematteis, A., 2007. Modeling transient discharge into a tunnel drilled in aheterogeneous formation. Groundwater 45 (6), 786–790.
Philip, J.R., 1998. Seepage shedding by parabolic capillary barriers and cavities. Water
Resour. Res. 34 (11), 2827–2835.
Polubarinova-Kochina, P.Y., 1962. Theory of Ground Water Movement 1977, second ed.
Princeton Univ. Press, Princeton, Nauka, Moscow (in Russian).
Preene, M., Chrimes, M., 2021. Groundwater lowering for construction of the Kilsby
tunnel–geological and geotechnical interpretation. In: Proceedings of the Institutionof Civil Engineers-Engineering History and Heritage, pp. 1–14.
Qin, Z., Wang, Y., Song, Y., Dong, Q., 2020. The analysis on seepage field of grouted andshotcrete lined underwater tunnel. Math. Probl. Eng. 1–15. https://doi.org/
10.1155/2020/7319054, 7319054.
Radcliffe, D., Šimůnek, J., 2010. Soil Physics with HYDRUS: Modeling and Applications.
CRC Press, Taylor & Francis Group, Boca Raton.
Shahbazi, A., Chesnaux, R., Saeidi, A., 2021. A new combined analytical-numericalmethod for evaluating the inflow rate into a tunnel excavated in a fractured rockmass. Eng. Geol. 283, 106003.
Shi, L., Lu, C., Ye, Y., Xie, Y., Wu, J., 2020. Evaluation of the performance of multiplewell hydraulic barriers on enhancing groundwater extraction in a coastal aquifer.
Adv. Water Res. 144, 103704.
Šimůnek, J., Šejna, M., van Genuchten, M.T., 2018. New features of version 3 of the
HYDRUS (2D/3D) computer software package. J. Hydrol. Hydromech. 66 (2),
133–142.
Strack, O.D.L., 2017. Analytical Groundwater Mechanics. Cambridge Univ. Press, NewYork.
Su, K., Zhou, Y., Wu, H., Shi, C., Zhou, L., 2017. An analytical method for groundwater
inflow into a drained circular tunnel. Groundwater 55 (5), 712–721.
Sweetenham, M.G., Maxwell, R.M., Santi, P.M., 2017. Assessing the timing andmagnitude of precipitation-induced seepage into tunnels bored through fractured
rock. Tunnel. Underground Space Technol. 65, 62–75.
Tang, Y., Chan, D.H., Zhu, D.Z., 2018a. Analytical solution for steady-state groundwater
inflow into a circular tunnel in anisotropic soils. J. Eng. Mech. 144 (9), 06018003.
ASCE.
Tang, Y., Chan, D.H., Zhu, D.Z., Guo, S., 2018b. An analytical solution for steady seepage
into a defective pipe. Water Sci. Technol. 18 (3), 926–935.
Tong, L., Zhou, H., Shen, H., Liu, H., 2022. A complex variable solution for shallow
rectangular tunnel in semi-infinite plane. J. Eng. Mech. 148 (4), 04022012. , ASCE.
Vedernikov, V.V., 1939. Theory of seepage and its application in irrigation and drainage.
Government Construction Technology Publishing House, Moscow, (in Russian).
Verruijt, A., Booker, J.R., 2000. Complex variable solution of Mindlin’s problem of an
excavated tunnel. In: Developments in Theoretical Geomechanics. Balkema,
Rotterdam, pp. 3–22.
Vuković, M., Soro, A., 1997. Groundwater Dynamics: Steady Flow. Highlands Ranch.
Water Resources Publications.
Wang, G.F., Wu, Y.X., Lu, L., Li, G., Shen, J.S., 2019. Investigation of the geological andhydrogeological environment with relation to metro system construction in Jinan,
China. Bul. Eng. Geol. Environ. 78 (2), 1005–1024.
Wang, L., Yu, L., Guo, Z., Wang, Z., 2014. Seepage induced soil failure and its mitigation
during suction caisson installation in silt. J. Offshore Mech. Arctic Eng. 136 (1),011103.
Wang, X., Wang, M., Zhang, M., Ming, H., 2008. Theoretical and experimental study ofexternal water pressure on tunnel lining in controlled drainage under high waterlevel. Tunnel. Underground Space Technol. 23 (5), 552–560.
Warrick, A.W., Fennemore, G.G., 1995. Unsaturated water flow around obstructions
simulated by two-dimensional Rankine bodies. Adv. Water Res. 18 (6), 375–382.
Wolfram, S., 1991. Mathematica. A System for Doing Mathematics by Computer.
Addison-Wesley, Redwood City.
Wu, J., Zhou, Z., Zhuang, C., 2021. A combined analytical-numerical method forgroundwater inflow into circular tunnels in drained conditions. Hydrogeol. J. 29 (7),2529–2543.
Wu, K., Cui, S., Liu, Y., Zhang, Q., Zhao, J., Zhang, Z., Han, Y., 2020. Study on themechanism of grouting under different tunnel depth of cross passage. Geotech. Geol.
Eng. 38 (3), 3291–3305.
Xu, P., Chen, C., Zhang, A., 2021. Research on construction technology and formation
deformation of shield tunneling under river. Geotech. Geol. Eng. 39 (5), 3627–3638.
Xue, Y., Kong, F., Li, S., Qiu, D., Su, M., Li, Z., Zhou, B., 2021. Water and mud inrush
hazard in underground engineering: genesis, evolution and prevention. Tunnel.
Underground Space Technol. 114, 103987.
Yang, G., Wang, X., Wang, X., Cao, Y., 2016. Analyses of seepage problems in a subsea
tunnel considering effects of grouting and lining structure. Mar. Georesour.
Geotechnol. 34 (1), 65–70.
Yang, J.J., Cheng, W.C., Wang, S., 2021. Advanced tunneling techniques and information
modeling of underground infrastructure. In: Proceedings of the 6th GeoChina
International Conference on Civil & Transportation Infrastructures: FromEngineering to Smart & Green Life Cycle Solutions – Nanchang, China, 2021. Cham.
Springer.

20

A.R. Kacimov et al.

Advances in Water Resources 164 (2022) 104182

Yang, S.Y., Yeh, H.D., 2007. A closed-form solution for a confined flow into a tunnel
during progressive drilling in a multi-layer groundwater flow system. Geophys. Res.
Lett. 34 (7), L07405.
Ying, H.W., Zhu, C.W., Shen, H.W., Gong, X.N., 2018. Semi-analytical solution forgroundwater ingress into lined tunnel. Tunnel. Underground Space Technol. 76,
43–47.
Zhang, D., Fang, Q., Lou, H., 2014. Grouting techniques for the unfavorable geologicalconditions of Xiang’an subsea tunnel in China. J. Rock Mech. Geotech. Eng. 6 (5),438–446.
Zhang, L., Franklin, J.A., 1993. Prediction of water flow into rock tunnels: an analyticalsolution assuming an hydraulic conductivity gradient. Int. J. Rock Mech. Mining Sci.
Geomech. 30 (1), 37–46.
Zhang, L., Zhao, D., Wu, J., Yang, W., Wang, W., Xin, D., 2020a. Prediction of water
inflow in Tsingtao subsea tunnel based on the superposition principle. Tunnel.
Underground Space Technol. 97, 103243.

Zhang, W., Dai, B., Liu, Z., Zhou, C., 2019. On the non-Darcian seepage flow field arounda deeply buried tunnel after excavation. Bull. Eng. Geol. Environ. 78 (1), 311–323.
Zhang, Y., Zhang, D., Fang, Q., Xiong, L., Yu, L., Zhou, M., 2020b. Analytical solutions of
non-Darcy seepage of grouted subsea tunnels. Tunnel. Underground Space Technol.
96, 103182.
Zhu, T., Xu, X., 2021. Watershed Erosion Processes. Springer, Cham.
Zhukovsky, N.E., 1949. Collected Works in Seven Volumes. V. 2 “Hydrodynamics”.
GITTL, Moscow (in Russian).
Zingg, S., Anagnostou, G., 2018. Tunnel face stability and the effectiveness of advancedrainage measures in water-bearing ground of non-uniform permeability. Rock
Mech. Rock Eng. 51 (1), 187–202.
Zverovich, E.I., Litvinchuk, G.S., 1968. Boundary-value problems with a shift for analytic
functions and singular functional equations. Russ. Math. Surv. 23 (3), 67.

21