HYDROLOGICAL SCIENCES JOURNAL
2021, VOL. 66, NO. 5, 795–808
https://doi.org/10.1080/02626667.2021.1890327

Water table rise in urban shallow aquifer with vertically-heterogeneous soils:
Girinskii’s potential revisited
A. R. Kacimova, A. Al-Maktoumi

a,b

and J. Šimůnek

ca

Department of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Muscat, Oman; bWater Research Centre, Sultan Qaboos
University, Muscat, Oman; cDepartment of Environmental Sciences, University of California Riverside, Riverside, California, USA
ABSTRACTARTICLE HISTORY

Seepage through an aquifer, the hydraulic conductivity of which varies vertically, is studied using the
Dupuit-Forchheimer approximation (Girinskii’s potential) and numerically by FDM-MODFLOW and FEMHYDRUS-2D. In urban water hydrology, the effect of compaction of the top stratum of an aquifer on theflow rate and the position of the water table (characterized by an integral quantity of the saturated/dryarea) is analysed. The area of the saturated zone is evaluated as a function of the conductivity ratio of thetwo strata, and their thicknesses. Conductivity varying linearly and exponentially is also analyticallystudied. Boundary-value problems are solved for a linear (nonlinear) ordinary differential equation if
the evaporation rate is constant (decreasing exponentially with depth). The MODFLOW and HYDRUS-2D
simulations give consistent fields of the piezometric head, pressure head, and Darcian velocity, as well as
the water table position, which may have a global minimum (watershed).

Received 2 April 2020
Accepted 26 June 2020

“From the earliest traceable cosmical changes down to the latest
results of civilization, we shall find that the transformation of thehomogeneous into the heterogeneous, is that in which Progressessentially consists.”
Spencer, H., 1875. Illustrations of Universal Progress: A Series ofDiscussions. New York: D. Appleton & Company.

1 Introduction
In urban hydrology, one concern regards alterations to thegroundwater system with anthropogenic impacts in generaland moisture fluxes, water table dynamics and piezometricsurface variations in sub-city soils, vadose zone and aquifers
in particular (see e.g. Lerner 2003, Vázquez-Suñé et al. 2005,
Lerner and Harris 2008, Naik et al. 2008, Quan et al. 2010,
Kazemi 2011, Chaudhary 2012, Howard and Israfilov 2012, Jhaet al. 2012, Barron et al. 2013, Schirmer et al. 2013, Preene andFisher 2015, Attard et al. 2016, McGrane 2016, Porse et al.
2016, Vogwill 2016, Medovar et al. 2018, Coda et al. 2019).
Howard (2017, p. 207) asked: Does urbanization of soils ele­
vate or depress the pristine water table? We recall Burghardt(1994, p. 207):
In many cases the groundwater-table drops in urban areas.
As a reason for this the reduction of water infiltration by sealing
was suggested (Brechtel 1982). But there exist also results from
Berlin, which show the opposite effect due to a diminishing of theevapotranspiration. An enhancement of the renewal of ground­water takes place (Renger et al. 1987).

EDITOR

A. Fiori
ASSOCIATE EDITOR

G. Chiogna
KEYWORDS

Dupuit-Forchheimer model;
unconfined aquifer;evaporation from a phreatic
surface; MODFLOW;HYDRUS-2D

Changes of “virgin” hydrological conditions/storage under
cities trigger numerous phenomena when the subsurface becomeswetter (damage to the waterlogged footings-foundations of build­ings, heaving of roads and other urban structures, increased flash
flood amplitudes, exfiltration of groundwater into local depres­sions and breeding of harmful insects, de-aeration of the rootzone of urban parks and gardens; salinization of the topsoil if
groundwater is saline, among others). A critical rise of shallow
water tables (also called “groundwater inundation”) in arid urbanzones is caused by several factors (leaks from breached subsurface
pipes and tanks, irrigation return from over-watered city parks/
gardens, paving of ground surface and ensued hampering of
evaporation, among others; see e.g. Abu-Rizaiza et al. 1989, AlRawas and Qamaruddin 1998, Abu-Rizaiza 1999, Al-Sefry andŞen 2006, Kreibich and Thieken 2008, Al-Senafy 2011, Al-Senafy
et al. 2015, Bob et al. 2016). Similarly, groundwater recharge inhumid regions can increase significantly when these areas are
growing, in particular due to a significant loss of evapotranspira­tion (Minnig et al. 2018).
In this paper, we analyse the impact of a special type ofheterogeneity, viz. vertically varying soil hydraulic properties,
on a shallow unconfined aquifer. In the Muscat Metropolitanarea, groundwater pumping stopped decades ago (when farms
vanished), concomitant with heavy irrigation of expanding cityparks and urban water supply coming from desalination
plants. In contrast, ornamental vegetation is watered by ter­tiary treated water provided to Muscat Municipality via cen­
tralized subsidies. Large areas of Big Muscat are now covered

CONTACT A. R. Kacimov
anvar@squ.edu.om ; akacimov@gmail.com
Department of Soils, Water and Agricultural Engineering, Sultan Qaboos University AlKhod 123, PO Box 34, Sultanate of Oman
Omani page: https://www.squ.edu.om/agr/Academics/Departments/Soil-Water-and-Agricultural-Engineering© 2021 IAHS

796

A. R. KACIMOV ET AL.

by urban structures (buildings, roads, car parks, etc.). Most
urban areas in the world, since the earliest epochs of trans­
forming civilizations (see the epigraph above), are subject to aheavy load of these structures that causes compaction of cities’
substrate and subsidence (i.e. topsoil’s bulk density and capil­
larity increase and its permeability decreases). We answer thefollowing question: Do these changes increase or decrease thenatural location of the phreatic surface?
The structure of our paper is as follows. In Section 2, werevisit analytical solutions of Polubarinova-Kochina (1962,1977) (hereafter abbreviated as PK-1962 and PK-77) to onedimensional (1D) Dupuit-Forchheimer (hereafter abbreviatedas DF) groundwater flows through a soil layer, the hydraulicconductivity of which varies with depth as follows:● piece-wise continuously (Subsection 2.1 – no evapora­tion; Subsection 2.2 – with evaporation from the watertable),● exponentially (Subsection 2.3, no evaporation), and● linearly (Subsection 2.4, no evaporation).
In Section 3, we apply MODFLOW to model flow in a twolayered aquifer – in particular, one subject to evaporation from ashallow water table. In Section 4, we use HYDRUS-2D to model
the same two-layered aquifer as in Subsection 2.1 and Section 3,
taking into account the capillarity of the two strata. In Subsection
4.1, we consider topsoil as an impermeable boundary, and in
Subsection 4.2, we assume a given flux and pressure head along
this boundary.

2 Analytical solutions for two-layered aquifers – DF
approximation
A vertical cross-section of the pristine (pre-compaction)
groundwater flow is shown in Fig. 1(a). We assume a horizontal
bedrock on the horizontal scale of a city, a (tens to hundreds ofmetres), and an unconfined aquifer with hydraulic conductivity
K0 subtended by a horizontal impermeable bedrock OD. As
commonly conceptualized at such a scale, we ignore capillarity.
City hydrologists usually monitor (with piezometers) a steadystate (or quasi-steady-state) groundwater motion over annual
periods. Two piezometers in Fig. 1(a) show the water levels h1
and h2. At a prototype site in Oman (the Sultan Qaboos
University campus), the water table depth under the groundsurface varies from about 10 cm to 2 m within a few hundred
metres (size a in Fig. 1) between observation points.
Polubarinova-Kochina (1962, p. 1977) solved a free boundary
problem analytically for a partial differential equation (theLaplace equation), which models a two-dimensional (2D) flowin Fig. 1(a) (the so-called “dam problem”), which involves aseepage face on the tailwater slope of the dam. If the slope of thephreatic surface is large, then a three-dimensional (3D) seepage
face emerges also in monitoring boreholes (piezometers), which
gauge 1D regional flows. As an interpose, we refer to figs 5.4, 5.6,
and 5.8 of Yeh et al. (2015), where unconfined flows in aquifers
are sketched analogously to fig. 217 of PK-62, which illustrates
seepage through an embankment.
We assume a horizontal bedrock OD and a flat groundsurface AE (Fig. 1(b,c)). We consider vertically averaged flow

Figure 1. Vertical cross-sections of (a) a homogeneous aquifer, (b) a two-layered aquifer with no evaporation from AE, and (c) a two-layered aquifer with evaporationfrom the ground surface.

HYDROLOGICAL SCIENCES JOURNAL

in the x-direction of Fig. 1(b), i.e. we employ the DF approx­
imation (valid for sufficiently small h1/a). Mathematically, aboundary value problem (BVP) for a linear DF-ODE ignores
the seepage face above point C (Fig. 1(b)). We examine three K(y) functions. The Charny formula (PK-62, PK-77) for ahomogeneous case (K0) and extensions of this formula toarbitrary K(y), even for flows subtended by an interface
above subjacent static saline water (Strack 2017, Toller andStrack 2019), give an exact value of the flow rate q through thesaturated zone (Fig. 1(a)). Remarkably, the DF approximation
results in the same value of q as the 2D theory, albeit thephreatic lines in the two models are different. In this paper,
we select the area of the dry zone above the phreatic surface(Fig. 1) as a criterion of groundwater inundation.

extremum of Aw(d) for the layering in Fig. 1(b). Philosophically,
we prove that the pinnacle of Spencer’s Progress (see the epi­
graph) is attained at a certain hierarchical organization of thecomposite in Fig. 1(b).
The Girinskii potential is:hðxÞF ð x Þ ¼ ò 0 ðyhðxÞÞK ð yÞdy(1)that for the two-strata aquifer in Fig. 1(b) yields:
�Ko 2
h ; 0 � y � ho ;
� 2 K1F ð x Þ ¼ K0 2
(2)
h
2hh
h0 Þ2 ; ho � y � h1 :
0
0
2
2 ðh
Without evaporation, F(x) obeys a linear ODE, which inte­
grates twice:d2 F¼ 0; ) F ¼ q x þ C
dx2

2.1 Two-strata unconfined aquifer with no evaporation
In this subsection, we revisit the Girinskii potential (see PK-62,
her chapter X, section 4, Example 21), which is used for solving
the DF flow problem depicted in Fig. 1(b). The total saturated
thickness and a commanding piezometric head of the aquifer are
h1, the thickness of the bottom layer of conductivity K0 is h0. The
downstream head, h2, obeys the double inequality h1 > h0 > h2. No
accretion (evaporation) takes place along the water table, which isa smooth curve, except at point B (the abscissa of this point is x0),
where the phreatic line ABC crosses a horizontal interface I1I2
between the two strata. The upper layer has conductivity K1.
Obviously, in the two limits of the thickness d of the top
layer, bounded by the double inequality 0 � d � h1 , viz. at
d = 0 and d = h1, the water table is the same because in

797

(3)

where the two constants of integration (q and C) are deter­
mined from the BVP:
� K1
K0 2
ðh1
h0 2h1 h0
2
2
�
h2
K0
K0 0 ) q ¼
2h1 h0 h20
22a

F ð0 Þ ¼ C ¼

h0 Þ2 ;F ða Þ ¼

K1
h22 þ ðh1
K0

�
h0 Þ

2

(4)
Physically, q in Equation (4) is the groundwater inflow rate
through AO in Fig. 1(b). Because evaporation, infiltration, andleakage through OD are ignored in Equation (3), q is the flowrate along the whole distance0 � x � a.
The water table is described by the following equation:

( �K
x¼

1
q

�
�
1
K0 h0 ðh h1 Þ ; 0 � x � xo ; h0 � h � h1 ;h0 Þ2 ðh h0 Þ2
2 ðh1
�
��
1
h0 Þ2 K0 h2 þ h20 2h1 h0 ;xo � x � a; h2 � h � h0
2q K1 ðh1homogeneous formations (e.g. earth dams, Fig. 1(a)) the shapeof the phreatic surface does not depend on the conductivity. We
can use the fully saturated area Aw(d) under the curve AC as aquantifier of waterlogging. Obviously, Aw(0) = Aw(h1), therefore,
from the Rolle’s theorem, there should be a certain value of d atwhich Aw attains at least one minimum or maximum. This isalso apparent from a geotechnical engineering practice: standard
anti-seepage measures in earth dams (clay cores, sheet pilings,
impermeable barriers inside the dam bodies, colmation, etc.)
reduce the seepage flow rate through the dam. However,
upstream of an engineered low-permeable zone (barrier), loca­
lized inside the dam, the phreatic surface is hoisted (as com­pared with a “pristine,” homogeneous soil), while downstreamof this “lens of altered permeability,” the water table drops.
Consequently, Aw may either increase or decrease as comparedwith the homogeneous case. In this subsection, we find such an
1

In PK-77 this example is dropped.

(5)

From Equation (5), x0 is evaluated at h = h0. The area of thesaturated zone under ABC and above the horizontal line y = h2
in Fig. 1(b) is expressed as:
h0

h1

Aw ¼ ò xu ðhÞdh þ ò xb ðhÞdh
h0

h2

(6)

where the integrands xu and xb are taken from the first andsecond parts of Equation (5).
We introduce the octad of dimensionless quantities:��� � � � �
h ; x ; h2 ; h0 ; a ¼ h;x;hh21;h0 ;a ; FK0¼F; kr ¼ K1 =K0 ; q� ¼ q=ðk0 h1 Þ;
h2
1A�w ¼ Aw =h21 and, for the sake of brevity, drop the * superscript.
We evaluate the integrals in Equation (6) and obtain:
Aw ðh0 ; a; kr ; h2 Þ ¼

a h30 ðkr 1Þ 3h0 ðkr 1Þ þ 2kr(7)
3 h20 ðkr 1Þ 2h0 ðkr 1Þ þ kr h22

798

A. R. KACIMOV ET AL.

It is obvious from Equation (7), as well as from Yakimov’s
comparison theorems (see e.g. Ilyinsky et al. 1998), that variationsof three independent variables in the tetrad, viz. a; kr ; and h2 ,generate monotonic variations of Aw. For example, if h2 increases
from 0 to 1 (i.e. h2 = h1) (provided a and kr are fixed), then Aw
increases. As mentioned in Section 1, the function Aw(h0) is lesstrivial. Specifically, the first derivative testðh 0
dAw ðh0 Þ
¼
dh0

1Þ

3 h22

3 h0 h22
3ð

h22

kr

�

2h0 ðkr

Equation (8) is a quartic equation with respect to h0. Twocomplex roots are discarded. We evaluated two real roots ofEquation (8) using the routine Solve of Wolfram’s (1991)
Mathematica.2 According to the Ferrari formula, the firstroot is h0 = 1, which is the limit of a homogeneous aquifer.
The second root (dropped here for the sake of brevity) gives a
nontrivial internal extremum.
Next, we differentiated Equation (8) and used the second
derivative test for Aw(h0). We found that the two real roots ofEquation (8) give two global extrema: end and internal ones.
The extrema are either minima or maxima, depending onwhether kr > 1 or kr < 1. At the smallest possible value of
h0 = h2, we have another end extremum (local).
Figure 2 (top panel) shows Aw(h0) for the dyad (h2,a) = (0.1,2) and kr = 0.01, 0.1, and 0.3 (curves 1–3, respectively). For curve
3, for example, the minimum and maximum are indicated byarrows on the third curve. These extrema are in comport withthe analytical values: the global internal minimum is attained at(Aw,h0) = (1.25,0.25), the end local minimum is at (Aw,h0) = (1.28,
0.1), and the end global maximum is (Aw,h0) = (1.35, 1) (here­
after we retain two significant digits in computed values). In Fig.
2 (bottom panel), the curves Aw(h0) are plotted for the dyad (h2,a) = (0.1, 2) and kr = 10, 6, and 2 (curves 1–3, respectively). A
dimensionless dry area, Ad, in Fig. 1(b) equals a(1 − h2) − Aw, i.e.
the extrema of Ad(h0) can be inferred from the analysis above,
modulo a linear transformation.

Usually, an equivalent (also called “effective” or “apparent”)
hydraulic conductivity of heterogeneous aquifers is defined for
confined flows (see chapter 6 in Yeh et al. 2015). We define thisquantity for an unconfined flow in Fig. 1(b) as:
Keq ¼

3 h20 ðkr1Þ þ h30 ðkr

1Þ þ h20 ðkr

1 Þ þ kr Þ

1Þ

q
h1

�
k r ðk r

2

(9)

h2

1ÞÞ

¼0

(8)

From Equation (4) we know q for a two-strata aquifer(Fig. 1), which we put into Equation (9). From the Dupuit
formula and Equation (9) we obtain Keq0 for a homoge­
neous aquifer. Then we obtain the ratio of equivalentconductivities:
λ¼

Keq
2 h0
¼
Keq0

h20

h22 þ kr ð1
1 h22

h0 Þ2

(10)

If little information is known about the vertical heterogeneity,then (considering a two-layered formation in Fig. 1(b)) thevalue of λ in Equation (10) assesses a relative error in theconceptualization of the aquifer as a homogeneous massif.
For example, for a dyad (h0, h2) = (0.7,0.5), the ratios λ are
0.88 and 12.88 for kr = 0.01 and 100, respectively. Therefore,
for a compacted (less permeable) top layer, the error in theconceptualized “homogenization” is much less than for thecase of a more permeable top stratum.

2.2 Two-strata aquifer with evaporation
In this subsection, we consider the case h1 > h2 > h0, i.e.
the DF water table is located above I1I2 in Fig. 1(c), andadd evaporation (of a constant intensity e) from AC. As
compared with fig. 5.9 of Yeh et al. (2015), the direction ofarrows (fluxes) is reversed. Although the curve h(x) is

Figure 2. Areas Aw of the saturated zone under ABC as a function of h0 for (top) (h2, a) = (0.1, 2) and kr = 0.01, 0.1, and 0.3 (curves 1–3, respectively) and (bottom) (h2,a) = (0.1, 2) and kr = 10, 6, and 2 (curves 1–3, respectively).

2

All programs, which we wrote in Wolfram’s Mathematica, Visual MODFLOW, and HYDRUS-2D, are open-access files and are available from the authors upon request.

HYDROLOGICAL SCIENCES JOURNAL

smooth, for sufficiently high values of e, a trough emerges
on the water table (Fig. 1(c)). The deepest point of thetrough is Rs, the abscissa and ordinate of which are xm andhm, respectively. Moreover, if e exceeds a certain thresholdvalue, then this trough point drops below I1I2. Then two
corner points emerge on the water table (for the sake ofbrevity we do not consider such intense evaporation). The
trough in Fig. 1(c) also implies a stagnation point Ps and awatershed RsPs (vertical segment), i.e. groundwater fluxes
from OA and DC horizontally infiltrate into the aquifer
from the left and right, respectively, and collide along RsPs
such that no cross-flow from OA to DC (as in Fig. 1(a,b))
takes place.
Physically, the evaporation intensity depends on the thicknessd(x) (Fig. 1(c)) of the vadose zone. PK-77 and Kacimov et al.
(2019a) considered a linearly and exponentially varying e[d(x)]that resulted in BVPs for linear and nonlinear ODEs. A constant e
is a reasonable approximation for shallow water tables, i.e. smalld(x).
The governing ODE (3) and its solution are transformed asfollows:
d2 F
¼
dx2

e; ) F ¼x2
e þq xþC
2(11)

where the Girinskii function is given by the second line in
Equation (2).
We add the constant e* = e/K0 to dimensional quantities ofSection 2.1, solve BVP for the parabolic function F(x) (11), andobtain:e a 1 h2þ½2h0 þ kr ð1 þ h2 2h0 Þ�;22a2q x e x2 ¼ 2h0 ð1 hÞ kr ðh h0 Þ2 þ kr ð1

q¼

799

(12)
h0 Þ2

Equation (12) manifests a DF hyperbola as the water table.
Obviously, for the ODE (11) with a “sink term,” q in Equation(12) is the groundwater flow rate only through the commanding
segment OA. This rate, as is evident from Equation (12),decreases with x, can attain zero (at RsPs in Fig. 1(c)), and caneven change the sign. This is illustrated in Fig. 3(a) where thewater tables h(x) are shown for the tetrad (h0, h2, e, a) = (0.8, 0.9,
0.1, 2) and four values of kr = 100, 10, 1, and 0.01 (curves 1–4,
correspondingly). We skip over a routine analysis of extrema ofh(x). The water table for the selected tetrad becomes nonmono­
tonic, even for a homogenous aquifer (curve 3). We report the
final minimum only for curve 4: (hm, xm) = (0.87, 1.4).
The dimensionless areas of the whole saturated zone, AW,
and of the dry zone, Ad, are defined as:a

AW ¼ ò hðxÞdx; AW ¼ a
0

AW

(13)

The integrand in the integral (13) is expressed by solving thequadratic equation in Equation (12). We substitute the physi­
cally meaningful root h(x) into Equation (13) and integrate.
For the sake of brevity, we skip over the explicit expression for
AW(h0, h2, kr, e, a). In Fig. 3(b), we plot the function AW(h0) for
the triad (h2, e, a) = (0.9, 0.1, 2) and three conductivity ratios
kr = 0.1, 0.5, and 10 (curves 1–3, respectively). The first two
curves depict increasing functions, and the third curve repre­
sents a decreasing function. Figure 3(c) shows AW(kr) for the

Figure 3. Two-layered aquifer with evaporation: (a) water tables h(x) for (h0, h2, e, a) = (0.8, 0.9, 0.1, 2) and kr = 100, 10, 1, and 0.01 (curves 1–4, respectively); (b) area ofthe saturated zone AW(h0) for (h2, e, a) = (0.9, 0.1, 2) and kr = 0.1, 0.5, and 10 (curves 1–3, respectively); and (c) areas of the saturated zone AW(h0) for (h2, e, a) = (0.9, 0.1,
2) and h0 = 0.1, 04, and 0.7 (curves 1–3, respectively).

800

A. R. KACIMOV ET AL.

same triad as in Fig. 3(b) and three values of h0 = 0.1, 04, and
0.7 (curves 1–3, respectively). Figure 3(c) manifests that thesaturated area is more sensitive to the variations of the thick­ness of the top layer d = h1 − h0 at kr < 1, i.e. just for the case ofcompaction (which is fairly trivial). We emphasize that one
must check that the water table (Equation (12)) lies above I1I2of the two strata, because otherwise a three-branched watertable has to be found, i.e. Equation (12) is not valid.

2.3 Aquifer with hydraulic conductivity varying
exponentially with depth
In this subsection, similarly to Ameli et al. (2016) and Toller
and Strack (2019), we tackle the case of K = K0e−r y, where r is
an arbitrary real constant. We add to the list of dimensionlessvariables of Subsection 2.1 the quantity r* = r h1 (and drop thesuperscript). For no evaporation, ODE (3) holds and integrates
as in Section 2.1, where Equation (1) is expressed as� hðxÞ �1�FðxÞ ¼ 2 1 e r hðxÞ ð1 þ r hðxÞÞ þe r hðxÞ
rr
0 � h � 1; 0 � x � a

�
1 ;
(14)

The constants of integration in Equation (3) are:Fð0Þ ¼ C ¼ r12 ð1 e r Þ 1r ;� rrhe ð1þrÞ e 2 ð1þrh2 Þ
q ¼ r1aþ e r h2 h2r�e r þ1h2(15)

The water table and the wetted area above the horizontal liney = h2 are:
�
1 e r ð1 þ r Þx¼
rq
Aw ¼ a2er

e r h ð1 þ r hÞ
þ e r hhrer h2 ð2 þ rð 2 þ rer ðh2 1ÞÞðh22rðer þ er h2 ð 1 þ rer ðh2 1ÞÞÞe r þ1

�
h

1 ÞÞ
(16)

Figure 4(a) illustrates the water tables calculated using Equation(16) for the dyad (h2, a) = (0.5, 2) and r = −10, 0.1, and 10 (curves
1–3, respectively). For too-steep water tables, like curve 1 in Fig. 4(a), the DF approximation is poor because it ignores a seepage
face, which physically emerges on DE (Fig. 1(b)).
The dependence of Aw on r is plotted in Fig. 4(b) for a = 2
and h2 = 0.2, 0.5, and 0.8 (curves 1–3, respectively).

Figure 4. An aquifer with exponential K(y): (a) water tables h(x) for h2 = 0.5, a = 2,
and r = −10, 0.1, and 10 (curves 1–3, respectively); and (b) areas Aw(r) of thesaturated zone for a = 2 and h2 = 0.2, 0.5, and 0.8 (curves 1–3, respectively).

�h22 þ b=3 1 h32
;
q¼
2a
�
1 h2 þ b=3 1 h3
x¼a
; 0 � x � a; h2 � h � 1; (17)
1 h22 þ b=3ð1 h32 Þ��ð1 h2 Þ 8 þ 4h2 þ b 3 þ 2h2 þ h22
Aw ¼ a12ð1 þ h2 Þ þ 4bð1 þ h2 þ h22 Þ
1

We used Equation (17) and plotted in Fig. 5 (top) the
qubic parabolas x(h) for the dyad (h2, a) = (0.1,2) andb = −0.99, 1, and 10 (curves 1–3, respectively). The areasof the wetted curvilinear triangle between the water table
and the horizontal line y = −h2 are shown in Fig. 5(bottom) as Aw(b) for h2 = 0.7, 0.8, and 0.9 (curves 1–
3, respectively).

2.4 Aquifer with hydraulic conductivity varying linearlywith depth

2.5 Evaporation rate varying with depth

In this subsection, we present calculations for Example 1 from
PK-62 (see her chapter X, section 4), i.e. K = K0 (1 + b y), whereb is a parameter quantifying the linear change of conductivity
with depth (we add a dimensionless parameter b* = b h1 tothose in Section 2.1 and drop the “*” superscript) and where
evaporation is neglected. Obviously, 1 � b � 1. We pre­sent only the final solution for the BVP:In this subsection, we consider the evaporation rate
decreasing with d (Fig. 1). Linear or nonlinear functionse(d) were studied (see e.g. Hu et al. 2008, Kats 1968, PK77). Without any loss of generality, here we select theexponential function (Kacimov et al. 2019a). We consider
the same two-layered aquifer as in Subsection 2.2 butwith h2 > h0. Then Equation (11) is transformed into

HYDROLOGICAL SCIENCES JOURNAL

801

Figure 5. Water table depths (top) and wetted areas above the tailwater horizon (bottom) for an aquifer with a linearly varying conductivity. Water table depths are for(h2, a) = (0.1,2) and b = −0.99, 1, and 10 (curves 1–3, respectively). Wetted areas are for a = 2 and h2 = 0.7, 0.8, and 0.9 (curves 1–3, respectively).

Figure 6. The water table in a rectangular composite having h1 = 10 m, a = 50 m,
h0 = 5 m, h2 = 6 m, e0 = 0.01 m/d, K0 = 7 m/d, and K1 = 5 m/d. Curves 1 and 2
correspond to f = 0 (a constant evaporation rate) and an exponentially decaying
rate with f = 0.5 1/m, respectively.

d2 F
¼
dx2

e0 exp½ f ðh1hðxÞÞ�

�

K1
ðh
2

K0 2
F ðx Þ ¼
h
2 0

2hh0

where (e0, f) is a dyad of empirical positive constants (theevaporation rate from the saturated soil surface and an expo­
nential decay factor). Equation (18) represents a BVP for anonlinear ODE. We used the routine NDSolve ofMathematica to solve it.
We used dimensional quantities because the right-hand side(RHS) in Equation (18) is an empirical evaporation function(Kats 1968, Hu et al. 2008). Specifically, we selected h1 = 10 m,
a = 50 m, h0 = 5 m, h2 = 6 m, e0 = 0.01 m/d, K0 = 7 m/d, andK1 = 5 m/d. Curves 1 and 2 in Fig. 6 present the water table for
f = 0 and 0.5 1/m, respectively. The two curves are almost
identical (e.g. at the middle point x = 25 m, the values of h are
8.14 and 8.17 m, respectively). We note that the “sink” term, i.e.
the RHS in Equation (18), can be any other nonlinear function.

3 MODFLOW modeling of an aquifer consisting of
two strata(18)2h0 Þ ; ho � y � h1

In this section, we numerically simulate the problem presented
in Fig. 1(b), using a finite difference method (FDM)

802

A. R. KACIMOV ET AL.

Figure 7. A vertical cross-section of a rectangular flow domain for the problem simulated by MODFLOW.

Figure 8. The dimensionless seepage flow rate q as a function of h0 for a = 5, h2 = 0.6, computed by MODFLOW and analytically (Equation 4) (curves 1 and 2,
respectively).

implemented in the code MODFLOW 2005 (Simcore Software2005). The model runs for steady-state seepage. We compare
the MODFLOW results with those obtained by analyticalsolutions derived in Section 2.
Figure 7 depicts a vertical cross-section of the MODFLOW
domain. The dyad of saturated conductivities is (K0,
K1) = (7.13, 0.25) m/d (selected from the HYDRUS soil catalog
for sand and loam, respectively). The effective porosity of bothlayers is ϴs = 0.43. The domain is gridded with 50 layers and250 columns, forming a total of 12 500 active cells, with cell
sizes of Δy = Δx = 0.2 m. The generalized conjugate gradientmethod with preconditioning, involving a “modified incom­
plete Cholesky” solver, was selected with a convergence criter­
ion of 0.001. The constant head on the left and rightboundaries have h1 = 10 m and h2 = 6 m, respectively (Fig.
7). The model was run for different thicknesses d of thecompacted top layer (Fig. 1(b)): 0, 0.4, 0.8, 1.6, 2.4, and 3.8 m.
In Fig. 8, Curve 1 presents the MODFLOW discharge q
converted into dimensionless quantities of Section 2.1 andplotted by the ListPlot of Mathematica using the six abovelisted values of h0 = h1 − d. The dimensional value of q in
MODFLOW for all simulated cases was extracted from thewater budget output file. Curve 2 in Fig. 8 is an analytical
function q(h0) plotted by the non-dimensional Equation (4).
Curves 1 and 2 are very close (the vertical scale in Fig. 8
exaggerates the differences). For instance, at h0 = 0.84, thenumerical and analytical values are 0.0601 and 0.0615, respec­tively, i.e. the difference is less than 3%.
The results of MODFLOW simulations for the case of
aquifer flow with evaporation are shown in Fig. 9(a–c).
Figure 9(a) illustrates the piezometric heads, and Fig. 9(b)plots the analytical and MODFLOW-computed water tables.
In Fig. 9(c), MODFLOW predicts the water table, which inter­
sects the horizon I1I2 (Fig. 1(c)) at two “corner points,” i.e. thetrough (Fig. 1(c)) descends into the substratum. This has beenachieved by lowering h2 from 9 (Fig. 9a) to 8.2 m. We note thatthe analytical solution for this case of the water table dipping
under I1I2 is not available.

HYDROLOGICAL SCIENCES JOURNAL

803

Figure 9. MODFLOW total heads, H(x,y), for (a) e = 0.0075 m/d, h0 = 8 m, a = 50 m, h1 = 10 m, h2 = 9 m, and (K0,K1) = (7.13, 0.25) m/d); (b) dimensionless water tablecalculated analytically by Equation (12) and using MODFLOW for e = 0.0075 m/d, h0 = 8 m, a = 50 m, h1 = 10 m, h2 = 9 m, and (K0, K1) = (7.13, 0.25) m/d; and (c)MODFLOW total heads, H(x,y), for e = 0.0075 m/d, h0 = 8 m, a = 50 m, h1 = 10 m, h2 = 8.2 m, (K0, K1) = (7.13, 0.25) m/d c).

Figure 10 presents the MODFLOW streamlines for the
scenario in Fig. 9(a). “Seeding points” (marked particles)
are released at depths of 3, 5, 7, and 9 m and distin­
guished by contrasting colours. The deepest water table

point calculated from the dimensionless value in Section
2.2 is 8.3 m, while MODFLOW gives 8.2 m. The x-coor­dinates of the minimum are 30.52 m and 30.6 m,
correspondingly.

804

A. R. KACIMOV ET AL.

Figure 10. Streamlines with “seeding points” (marked particles) released at
depths of 2, 4, 6, and 8 m for e = 0.0075 m/d, h0 = 8 m, a = 50 m, h1 = 10 m,
h2 = 9 m, and (K0,K-1) = (7.13, 0.25) m/d, computed by MODFLOW.

4 HYDRUS modeling of an aquifer consisting of two
strataIn this section, we use HYDRUS-2D (Šimůnek et al. 2016), whichtakes into account a saturated-unsaturated flow in the two-layered
formation of Fig. 1(b). A standard notation for the HYDRUS
ordinate axis is z. We will use non-standard (for HYDRUS)
notations: p for the pressure head and H = p − z for the total(piezometric) head.
4.1 No evaporation from AE
In this subsection, we use the same composite rectangle as in the
MODFLOW model above. Specifically, h1 = |OA| = 10 m,|OD| = 50 m, h0 = |OI1| = |OI2| = 7.5 m, h2 = 6 m (the origin of
HYDRUS coordinates is at point O and all notations are the sameas in Fig. 1). The upper layer and substratum are loam and sand,respectively, from the HYDRUS Soil Catalogue (Carsel andParrish 1988). The two pentads of soil hydraulic parameters forsand and loam are (Ks, ϴs, ϴr, α, n) = (7.13 m/d, 0.43, 0.045, 14.5
1/m, 2.68) and (0.25 m/d, 0.43, 0.078, 3.6 1/m, 1.56), respectively.
In the HYDRUS project, we select the following boundary
conditions: no flow along OD, and AE, H = 10 m and H = 6 malong OA and DC, respectively. The vertical segment CE is a
no-flow boundary. The regular finite element discretization
was used with a horizontal step of 20 cm and a vertical stepof 10 cm. The number of mesh entities is characterized by thetetrad (nodes, 1D elements, 2D elements, boundary meshnodes) = (25 351, 700, 50 000, 700). The simulation time is
100 d to reach the steady-state condition. The triad of iteration
criteria is (a maximum number of iterations, the water content
tolerance, the pressure head tolerance) = (10, 0.0001, 0.1 cm).
Steady-state seepage is established in about 17 d (the resultspresented below are for t = 100 d). While Fig. 11(a) shows thestreamlines, Fig. 11(b,c) present the fields of the vertical andhorizontal components, respectively, of the Darcian velocity
vector. From Fig. 11(b,c) one can infer that groundwater flow
is prevalently horizontal. In other words, the capillarity of the
topsoil layer does not lead to a significant violation of theassumptions of the DF model, which we used in Section 2.
Figure 11(d) illustrates the spatial field of the volumetric watercontent. From this coloured map, one can see a dry zone above
the water table, which dips under the line I1I2 (Fig. 1) at
x0 = 32.8 m. Remarkably, for x > x0 in Fig. 11(c), the loam
cap is wetter than the unsaturated zone in the sand substratum.
An intriguing question was also raised: Can the HYDRUS
flow rate q through the two layers in Fig. 11 be a nonmonotonic
function of d? We recall that in Kacimov et al. (2019b), nontrivial extrema of the flow rate were discovered for 1D verticalevaporation in a two-layered soil. In contrast, for basicallyhorizontal and prevalently saturated flow in Fig. 11, the
HYDRUS computations gave only a monotonic function q(d).

Figure 11. Streamlines (a), vertical (b) and horizontal (c) components of the Darcian velocity, and volumetric water contents (d) in the HYDRUS simulation for a twolayered domain 50 m long and 10 m tall, without evaporation.

HYDROLOGICAL SCIENCES JOURNAL

805

Figure 12. Water flow in a composite rectangle having |OA| = 20 cm, |OD| = 100 cm, |OI1| = |OI2| = 10 cm and a pressure head boundary condition along A1E with
p = −100 cm for t = 60 min: (a) HYDRUS velocity vectors field, (b) HYDRUS velocity distribution along OD, (c) HYDRUS water content field, (d) HYDRUS isobars
sandwiching the water table p = 0, and (e) the analytically calculated (Equation 12) water table in dimensionless variables corresponding to dimensional parameters in
Fig. 11(a).

4.2 Evaporation from AE
In this subsection, we use a composite rectangle different
from the rest of our paper. Specifically, h1 = |OA| = 20 cm,|OD| = 100 cm, h0 = |OI1| = |OI2| = 10 cm, h2 = 15 cm.
The upper layer and substratum are again loam and sand,respectively. Such a small size of the domain was selected
to amplify the effect of the soil capillarity, which is not
immediately visible on the much larger scale used in thearticle so far. It is well known (see e.g. experiments of
Silliman et al. 2002) that the impact of the capillary fringe
and the unsaturated zone on groundwater flow is not
scalable.
In the first HYDRUS project (Fig. 12), we select the follow­
ing boundary conditions: no flow along OD, H = 20 cm andH = 15 cm along OA and DB, respectively. Along AE, we
impose the pressure head boundary condition with
p = −100 cm.
In Fig. 12(a), we selected a horizontal segment AA1 (z = 20 cm,
0 < x < 8.05 cm), along which we imposed a no-flow condition.
This condition is physically justified by a very shallow water table

position in this area (Fig. 1). Indeed, Riesenkampf (see PK-62)
used a potential model of 2D seepage from a channel of a small
depth and inserted into his analytical model a horizontal no-flow
segment, along which tension-saturated water spreads laterally
from a point where p = 0 (the channel’s rim is analogous to point
C in our Fig. 1(b,c) to a point on the upper boundary of thecapillary fringe where p = −pc (pc is a soil capillary constantinversely related to van Genuchten’s α).
Similarly, the vertical segment CE in Fig. 12a (Fig. 1(c)) is also
a no-flow line. Segments AA1 and CE correspond to the stream­lines in the Vedernikov-Bouwer steady-state model (see e.g.
Vedernikov 1949, Kacimov et al. 2019c). They separate segmentsof constant but contrasting piezometric heads – i.e. without theseno-flow segments, a non-physical “shortcut” would take place.
The size of finite elements is about 2 cm, and the number ofmesh entities is characterized by the tetrad (nodes, 1D elements,
2D elements, and boundary mesh nodes) = (7468, 507, 14,576,
358). The simulation time is 60 min. The triad of iteration criteria
is (a maximum number of iterations, the water content tolerance,
and the pressure head tolerance) = (10, 0.0001, 1 cm).

806

A. R. KACIMOV ET AL.

Figure 13. Water flow in a composite rectangle having |OA| = 20 cm, |OD| = 100 cm, |OI1| = |OI2| = 10 cm and a flux boundary condition along A1E with e = 0.01 cm/min
for t = 60 min: (a) HYDRUS isobars sandwiching the water table p = 0, (b) HYDRUS pressure head distribution p(x) at a horizontal cross-section z = 12.86 cm, and (c)MODFLOW water table.

A steady-state seepage is established in about 30 min. Below,
the results are presented for t = 60 min. Figure 12(a) shows thevector field of Darcian velocities. Figure 12(b) presents the
distribution of velocities along OD. Point Ps (x = 70 cm, z = 0)
in Fig. 12(a) corresponds to a stagnation point in Fig. 12(b). The
water content field is shown in Fig. 12(c). Figure 12(d) illustrates
the phreatic surface AB, which is shown as a “thin band” plottedin HYDRUS by selecting the pressure heads varying in thebounds: −1 cm < p < 1(c)m. Point Rs (Fig. 1(c)) has the coordi­
nates xm = 70 cm and z = hm = 12.2 cm. For comparison, ananalytical water table (Section 2) is shown in Fig. 12(e) (indimensionless quantities). The analytical trough (point Rs) ischaracterized, upon conversion to dimensional quantities, bycoordinates of xm = 61.54 cm and z = hm = 14.5 cm.
In another HYDRUS project of this subsection, illustrated
in Fig. 13(a,b), the following boundary conditions were used:
no flow along OD, H = 20 cm and H = 15 cm along OA andDB, respectively. Along AE, we imposed the atmosphericboundary condition with the evaporative flux e = 0.01 cm/
min (which involved the minimum allowed pressure head
hCritA = 1000 m). The vertical segment CE was a seepage
face. Figure 13(a) shows the water table sandwiched between
the isobars p = 0.1 cm and p = −0.1cm. A cross-section line
chart plotted between the points (0, 12.86 cm) and (100 cm,
12.86 cm) (i.e. along a horizontal line with z = 12.86 cm) is
shown in Fig. 13(b). From this curve, we evaluated the ordinateof point Rs as xm = 75.86 cm. For comparison, we did
MODFLOW simulations for the small rectangular compositeof the same sizes as in HYDRUS, with grid sizes dx = dz = 0.25 cm(80 layers and 400 columns). Note that the main difference
between HYDRUS and MODFLOW simulations is thatHYDRUS considers flow in the unsaturated zone and assumes
that evaporation occurs at the soil surface, while in
MODFLOW, evaporation occurs at the groundwater table.
Figure 13(c) represents the water table, which has hm = 12.82 cm
at xm = 75.3 cm. From comparing Fig. 13(a and c), we canconclude that capillarity has a relatively minor impact on thewater table.

5 Concluding remarks
Soils, the vadose zone, and aquifers have vertically varying
hydraulic (saturated conductivity) and capillary (vanGenuchten parameters) properties that are determined bynatural processes of soil formation and anthropogenic altera­tions of the subsurface, e.g. by urban development, that fallwithin Spencer’s paradigm. In this paper, we used analyticaland numerical models to answer the main question: how thevertical variation of hydraulic-capillary properties of an aquifer
and a superjacent unsaturated zone affects a shallow watertable and other Darcian flow characteristics (velocity, total
and pressure head, water content) as compared with what the
Dupuit parabola predicts for a homogeneous aquifer in a 1D
steady state flow. We used the Girinskii potential and analyti­
cally (quasi-analytically) solved ODEs, which gave the watertable in a two-layered aquifer and an aquifer with conductivity
linearly and exponentially varying with depth.
Both constant-rate evaporation from the water table and theevaporation rate decreasing with the depth of the water table
under the ground surface were tackled. In the case of strong
evaporation (a shallow water table or evaporation exacerbatedby transpiration of Prosopis trees growing in the study area),
the water table has a global minimum between two verticalsegments with constant piezometric heads. Also, a nontrivial
extremum of the area of the saturated zone (a curvilinear
trapezium bounded by the inlet and outlet segments, thewater table, and the bedrock) is found for a two-strata aquifer.
The MODFLOW and HYDRUS-2D modeling of the twolayered aquifer without and with evaporation from the watertable or the soil surface, respectively, gave results consistent
with the analytical model. Waterlogging caused by ground­water inundation in arid climates is a daunting problem forcity managers and civil engineers, who can employ mathema­tical models as tools for the following:(a) Prediction of saturated-unsaturated flow characteristics
at various h1 and h2 values in Fig. 1.

HYDROLOGICAL SCIENCES JOURNAL(b) Contemplation of engineering measures of managed
aquifer recharge for mitigation of waterlogging, e.g.
horizontal or vertical drains, hydraulic barriers against
groundwater, and bioengineering (e.g. cultivation ofsalt-tolerant plants like Australian reeds, aimed at
intensive transpiration), among others.

807

Funding
Support from the Sultan Qaboos University via the grant “Rise of Watertableand its Mitigation at SQU Campus”, Project IG/VC/WRC/21/1, andby the SQU/DPS grant DR/RG/17 is appreciated.

ORCID
Applicability of the analytical and numerical models in thispaper to real-world aquifers depends on the quality of the fielddata to be collected from boreholes of a network of monitoring
wells (depths h1 and h2, thicknesses d and h0 of the strata, theirconductivities (see Fig. 1), and capillary constants of the top
layer in case of the HYDRUS-based analysis). Empirical con­stants r and b for the hydraulic conductivity, varying continu­
ously with depth, are, to the best of our knowledge, notmeasured in the field or lab but are tackled as fitting constants(see e.g. Toller and Strack 2019) that require an intricate
upscaling of discretely observed values of hydraulic propertiesof core samples and piezometric heads collected from different
depths in differential piezometers.
An empirical coefficient f of depth-dependent evaporationfrom a non-horizontal water table is another fitting parameter ofthe DF model, which has never before been correlated with fieldmeasurements of the moisture content (or capillary pressure) of
the topsoil and HYDRUS-2D simulations – which, ideally,should involve atmospheric boundary conditions. If the size ain Fig. 1 is large enough (say, tens of kilometres), then theseboundary conditions should be collected from several weather
stations located along AE in Fig. 1. What is an optimistic aspectof a generally pessimistic view on the feasibility of collecting
such extensive hydrologic and hydraulic data on subsurface
characteristics of the aquifers? Our models in this paper – bothanalytical and numerical – are simple and can be easily realized
by an urban hydrologist. The results from our models illustrate
non-trivial extrema and, as back-of-an-envelope estimates, canguide practitioners in their programme of data collection andfield experiments.

Highlights● Numerically and analytically, the impact of vertical variation of
hydraulic-capillary properties in an unconfined aquifer and a super­jacent unsaturated zone on velocity, total and pressure head, and watercontent fields is evaluated.
● The water table can have a global minimum in the case of strong
evaporation.
● A nontrivial extremum of the area of the saturated zone (a curvilinear
trapezium bounded by the inlet and outlet segments, the water table,
and the bedrock) is found for a two-strata aquifer.

Acknowledgements
Support from Sultan Qaboos University via the grant “Rise of Water-table
and its Mitigation at SQU Campus” and by the SQU/DPS grant DR/RG/
17 is appreciated. Helpful comments by two anonymous referees are also
appreciated.

Disclosure statement
No potential conflict of interest was reported by the authors.

A. R. Kacimov
http://orcid.org/0000-0003-2543-3219
A. Al-Maktoumi
http://orcid.org/0000-0001-8766-5993
J. Šimůnek
http://orcid.org/0000-0002-0166-6563

References
Abu-Rizaiza, O.S., 1999. Threats from groundwater table rise in urban
areas in developing countries. Water International, 24 (1), 46–52.
doi:10.1080/02508069908692133
Abu-Rizaiza, O.S., Sarikaya, H.Z., and Khan, M.A., 1989. Urban ground­water rise control: case study. Journal of Irrigation and Drainage
Engineering ASCE, 115 (4), 588–607. doi:10.1061/(ASCE)0733-9437(1989)115:4(588)
Al-Rawas, A.A. and Qamaruddin, M., 1998. Construction problems ofengineering structures founded on expansive soils and rocks in north­ern Oman. Building and Environment, 33 (2–3), 159–171. doi:10.1016/
S0360-1323(97)00048-6
Al-Sefry, S.A. and Şen, Z., 2006. Groundwater rise problem and risk
evaluation in major cities of arid lands–Jedddah case in Kingdom of
Saudi Arabia. Water Resources Management, 20 (1), 91–108.
doi:10.1007/s11269-006-4636-2
Al-Senafy, M., 2011. Management of water table rise at Burgan oil field,
Kuwait. Emirates Journal for Engineering Research, 16 (2), 27–38.
Al-Senafy, M., et al., 2015. Causes of groundwater rise at Al-Qurain
residential area, Kuwait. Procedia Environmental Sciences, 25, 4–10.
doi:10.1016/j.proenv.2015.04.002
Ameli, A.A., McDonnell, J.J., and Bishop, K., 2016. The exponential
decline in saturated hydraulic conductivity with depth: a novel method
for exploring its effect on water flow paths and transit time distribu­
tion. Hydrological Processes, 30 (14), 2438–2450. doi:10.1002/
hyp.10777
Attard, G., et al., 2016. Impact of underground structures on the flow ofurban groundwater. Hydrogeology Journal, 24 (1), 5–19. doi:10.1007/
s10040-015-1317-3
Barron, O.V., Donn, M.J., and Barr, A.D., 2013. Urbanisation and shallow
groundwater: predicting changes in catchment hydrological responses.
Water Resources Management, 27 (1), 95–115. doi:10.1007/s11269-0120168-0
Bob, M., et al., 2016. Rising groundwater levels problem in urban areas: acase study from the Central Area of Madinah City, Saudi Arabia.
Arabian Journal for Science and Engineering, 41 (4), 1461–1472.
doi:10.1007/s13369-015-1976-3
Burghardt, W., 1994. Soils in urban and industrial environments. Zeitschrift
für Pflanzenernährung und Bodenkunde, 157 (3), 205–214. doi:10.1002/
jpln.19941570308
Carsel, R.F. and Parrish, R.S., 1988. Developing joint probability distribu­
tions of soil water retention characteristics. Water Resources Research,
24 (5), 755–769. doi:10.1029/WR024i005p00755
Chaudhary, M.T.A., 2012. Implications of rising groundwater level on
structural integrity of underground structures – investigations andretrofit of a large building complex. Structural Survey, 30 (2), 111–129. doi:10.1108/02630801211228725
Coda, S., et al., 2019. Uplift evidences related to the recession ofgroundwater abstraction in a pyroclastic-alluvial aquifer ofSouthern Italy. Geosciences, 9 (5), 215, 1–15. doi:10.3390/
geosciences9050215
Harbaugh, A.W. (2005) MODFLOW-2005, the U.S. Geological Survey
modular ground-water model—the ground-water flow process. U.S.
Geological Survey Techniques and Methods 6-A16. U.S. Geological
Survey, Reston. https://pubs.usgs.gov/tm/2005/tm6A16/,

808

A. R. KACIMOV ET AL.

Howard, J., 2017. Anthropogenic soils. Cham, Switzerland: Springer.
Howard, K.W. and Israfilov, R.G., eds., 2012. Current problems of hydro­
geology in urban areas, urban agglomerates and industrial centres(NATO Science Series, Vol. 8). Dordrecht: Springer.
Hu, S., et al., 2008. Theoretical analysis of the limiting rate of phreatic
evaporation for aeolian sandy soil in Taklimakan Desert. Chinese
Science Bulletin, 53 (2), 119–124.
Ilyinsky, N.B., Kacimov, A.R., and Yakimov, N.D., 1998. Analytical solu­
tions of seepage theory problems. Inverse methods, variational theo­rems, optimization and estimates (A review). Fluid Dynamics, 33 (2),
157–168. doi:10.1007/BF02698697
Jha, A.K., Bloch, R., and Lamond, J., 2012. Cities and flooding: a guide tointegrated urban flood risk management for the 21st century.
Washington, DC: The World Bank.
Kacimov, A., Obnosov, Y.V., and Simunek, J., 2019b. Minimal evapora­tion by optimal layering of topsoil: revisitng Ovsinsky’s smart mulch­
ing-tillage technology via Gardner-Warrick’s unsaturated analyticalmodel and HYDRUS. Water Resources Research, 55 (5), 3606–3618.
doi:10.1029/2018WR024025
Kacimov, A.R., Obnosov, Y., and Al-Maktoumi, A., 2019a. Dipolic MAR“bubble” inside confined brine formation or floating “lens” on top of
unconfined saline aquifer. In: E.F. Escalante, ed. Proceedings of theInternational Symposium on Managed Aquifer Recharge (ISMAR 10)“Managed Aquifer Recharge: Local solutions to the global water crisis”.
Madrid, Spain, 333–340.
Kacimov, A.R., Obnosov, Y.V., and Or, D., 2019c. Evaporation-induced
capillary siphoning through hydraulically connected porous domains:the Vedernikov–Bouwer model revisited. Transport in Porous Media,
129 (1), 231–251. doi:10.1007/s11242-019-01285-z
Kats, D.M., 1968. The regularities of ground water evaporation in irrigation
lands of the arid zone. Water in the Unsaturated Zone, 83, 822–827.
Kazemi, G.L., 2011. Impacts of urbanization on the groundwater
resources in Shahrood, Northeastern Iran: comparison with otherIranian and Asian cities. Physics and Chemistry of the Earth, 36, 150–
159. doi:10.1016/j.pce.2010.04.015
Kreibich, H. and Thieken, A.H., 2008. Assessment of damage caused byhigh groundwater inundation. Water Resources Research, 44, W09409.
doi:10.1029/2007WR006621
Lerner, D.N., ed., 2003. Urban groundwater pollution: IAH international
contributions to hydrogeology 24. Lisse: Balkema.
Lerner, D.N. and Harris, B., 2008. The relationship between land use andgroundwater. resources and quality. Land Use Policy, 26S (2009), S265–
S273.
McGrane, S.J., 2016. Impacts of urbanisation on hydrological and water
quality dynamics, and urban water management: a review. Hydrological
Sciences Journal, 61 (13), 2295–2311. doi:10.1080/02626667.2015.1128084
Medovar, Y.A., Yushmanov, I.O., and Bronskaya, E.E., 2018. Changes
in the ecological condition of a territory at the construction ofbuildings with a deep foundation under various geological condi­
tions in Moscow. Water Resources, 45 (2), 65–72. doi:10.1134/
S0097807818060258
Minnig, M., et al., 2018. Impact of urbanization on groundwater recharge
rates in Dübendorf, Switzerland. Journal of Hydrology, 563, 1135–1146.
doi:10.1016/j.jhydrol.2017.09.058
Naik, P.K., et al., 2008. Impact of urbanization on the groundwater regime
in a fast growing city in central India. Environmental Monitoring andAssessment, 146 (1–3), 339–373. doi:10.1007/s10661-007-0084-6
Polubarinova-Kochina, P.Y., 1962. Theory of ground water movement.
Second ed. Princeton: Princeton Univ. Press. 1977, Theory of
Ground Water Movement. Nauka, Moscow (in Russian).
Porse, E., et al., 2016. Pumping for the masses: evolution of groundwater
management in metropolitan Los Angeles. GeoJournal, 81 (5), 793–809.
doi:10.1007/s10708-015-9664-0
Preene, M. and Fisher, S., 2015. Impacts from groundwater control inurban areas. Proceedings of the XVI ECSMGE GeotechnicalEngineering for Infrastructure and Development, 2847–2852. https://www.icevirtuallibrary.com/isbn/9780727760678
Quan, R.S., et al., 2010. Waterlogging risk assessment based on land use/coverchange: a case study in Pudong New Area, Shanghai. Environmental Earth
Sciences, 61 (6), 1113–1121. doi:10.1007/s12665-009-0431-8
Schirmer, M., Leschik, S., and Musolff, A., 2013. Current research inurban hydrogeology–A review. Advances in Water Resources, 51,
280–291. doi:10.1016/j.advwatres.2012.06.015
Silliman, S.E., et al., 2002. Fluid flow and solute migration within thecapillary fringe. Ground Water, 40 (1), 76–84. doi:10.1111/j.17456584.2002.tb02493.x
Šimůnek, J., Van Genuchten, M.T., and Šejna, M., 2016. Recent develop­
ments and applications of the HYDRUScomputer software packages.
Vadose Zone Journal, 15 (7), 1–25. doi:10.2136/vzj2016.04.0033
Strack, O.D., 2017. Vertically integrated flow in stratified aquifers. Journalof Hydrology, 548, 794–800. doi:10.1016/j.jhydrol.2017.01.039
Toller, E.A. and Strack, O.D., 2019. Interface flow with vertically varying
hydraulic conductivity. Water Resources Research, 55 (11), 8514–8525.
doi:10.1029/2019WR024927
Vázquez-Suñé, E., Sánchez-Vila, X., and Carrera, J., 2005. Introductory reviewof specific factors influencing urban groundwater, an emerging branch ofhydrogeology, with reference to Barcelona, Spain. Hydrogeology Journal,
13 (3), 522–533. doi:10.1007/s10040-004-0360-2
Vedernikov, V.V., 1949. Seepage with a draining or water-bearing stra­tum. In Doklady AN SSSR, 69 (5), 619–622(in Russian).
Vogwill, R., ed., 2016. Solving the groundwater challenges of the 21st
century. Vol. 22. Boca Raton: CRC Press.
Wolfram, S., 1991. Mathematica. a system for doing mathematics bycomputer. Redwood City: Addison-Wesley.
Yeh, T.C., Khaleel, R., and Carroll, K.C., 2015. Flow through heterogeneous
geologic media. New York: Cambridge University Press.
,