Advances in Water Resources 145 (2020) 103732

Contents lists available at ScienceDirect

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

Seepage to ditches and topographic depressions in saturated andunsaturated soils
A.R. Kacimov a,∗, Yu.V. Obnosov b,1, J. Šimůnek caDepartment of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Muscat, Oman
Institute of Mathematics and Mechanics, Kazan Federal University, Kazan, Russiac
Department of Environmental Sciences, University of California Riverside, CA, United Statesba r t i c l e

i n f oKeywords:
Analytic and HYDRUS solutions for Darcian
2-D and axisymmetric ﬂows in saturated andunsaturated soils towards drainage ditches andtopographic depressions
Evaporation and seepage exﬁltration fromshallow groundwaterComplex potential and conformal mappings
Method of images with sinks and sources for
the Laplace equation and ADE
Boundary value problems involving seepage
faces on Earth and Mars
Isobars, isotachs, constant piezometric head,
and Kirchhoﬀ potential linesa b s t r a c t
An isobar generated by a line or point sink draining a conﬁned semi-inﬁnite aquifer is an analytic curve, to whicha steady 2-D plane or axisymmetric Darcian ﬂow converges. This sink may represent an excavation, ditch, or wadi
on Earth, or a channel on Mars. The strength of the sink controls the form of the ditch depression: for 2-D ﬂow,
the shape of the isobar varies from a zero-depth channel to a semicircle; for axisymmetric ﬂow, depressions as ﬂat
as a disk or as deep as a hemisphere are reconstructed. In the model of axisymmetric ﬂow, a ﬁctitious J.R. Philip’s
point sink is mirrored by an inﬁnite array of sinks and sources placed along a vertical line perpendicular to ahorizontal water table. A topographic depression is kept at constant capillary pressure (water content, Kirchhoﬀ
potential). None of these singularities belongs to the real ﬂow domain, evaporating unsaturated Gardnerian soil.
Saturated ﬂow towards a triangular, empty or partially-ﬁlled ditch is tackled by conformal mappings and thesolution of Riemann’s problem in a reference plane. The obtained seepage ﬂow rate is used as a right-hand sidein an ODE of a Cauchy problem, the solution of which gives the draw-up curves, i.e., the rise of the water levelin an initially empty trench. HYDRUS-2D computations for ﬂows in saturated and unsaturated soils match well
the analytical solutions. The modeling results are applied to assessments of real hydrological ﬂuxes on Earth andpaleo-reconstructions of Martian hydrology-geomorphology.

1. Introduction“I had to live in the desert before I could understand the full value of grassin a green ditch.”Ella Maillart
Analytical models of 2-D seepage towards drainage ditches andtrenches, constructed by civil, geotechnical, and agricultural engineers,
used the machinery of complex variables (Anderson, 2013, Aravin andNumerov, 1953, Bear, 1972, Kirkham and Powers, 1972, PolubarinovaKochina, 1962, 1977, hereafter abbreviated as PK-62,77, Skaggs et al.,1999, Strack, 1989, Vedernikov, 1939), in particular, by tackling freeboundaries of Darcian ﬂows, the so-called phreatic surfaces.
We recall (see, e.g., Radcliﬀe and Šimůnek, 2010) that transient, 3D, saturated-unsaturated ﬂows in porous media (when both water andsoil are incompressible) obey the Richards equation:
𝜕𝜃= ∇(𝐾(𝑝)∇ℎ)𝜕𝑡where 𝜃(t, x, y, Z) is the volumetric moisture content, ∇ is the nabla
operator (in Cartesian or cylindrical coordinates), K(p, x, y, Z) is thehydraulic conductivity function, h(t, x, y, Z) = p + Z is the total head,p is the pressure head, Z is a vertical coordinate, and p(𝜃) is a capillarypressure (water retention) function, ﬁxed for each soil. Eq. (0) involves→→

the Darsy law 𝑉 = −𝐾(𝑝)∇ℎ, where 𝑉 is the Darcian ﬂux vector, and theprinciple of mass conservation.
Boundary value problems (hereafter abbreviated as BVPs) are solvedfor Eq. (0) by specifying initial conditions, e.g., 𝜃(0, x, y, Z), as well as
imposing physically meaningful boundary conditions (e.g., Dirichlet’s,
Neumann’s). Only numerical codes like HYDRUS-3D (Šimůnek et al.,2016) tackle such problems for arbitrary 3D transient ﬂows. Eq. (0) is
a highly nonlinear parabolic PDE. For steady ﬂows, its LHS vanishes,
and the equation becomes elliptic. If the ﬂow is purely saturated andthe porous medium is homogeneous, then K(p)=Ks , where Ks is the sat-(0)Abbreviations: ADE =, advective dispersion equation; BVP =, boundary value problem; DF =, Dupuit-Forchheimer; LHS, RHS=, left hand side, right hand side; ODE=, Ordinary Diﬀerential Equation; OSDP =, Optimal Shape Design Problem; PK-62,77 =, Polubarinova-Kochina, P.Ya., 1962. Theory of Ground Water Movement.
Princeton University Press, Princeton. The second edition of the book (in Russian) was published in 1977, Nauka, Moscow.; VG =, van Genuchten.
∗
Corresponding author at: Department of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Al-Khod 123, PO Box 34, Sultanate of Oman.
E-mail addresses: anvar@squ.edu.om (A.R. Kacimov), yobnosov@kpfu.ru (Yu.V. Obnosov), jsimunek@ucr.edu (J. Šimůnek).
1
ORCID ID:orcid.org/0000-0001-9220-7989
https://doi.org/10.1016/j.advwatres.2020.103732
Received 5 May 2020; Received in revised form 21 July 2020; Accepted 11 August 2020
Available online 12 August 20200309-1708/© 2020 Elsevier Ltd. All rights reserved.

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

urated hydraulic conductivity, and Eq. (0) is reduced to a linear Laplace
equation. For a special class of soils with an exponential K(p) function(Gardner, 1958), the Kirchhoﬀ transform makes possible a linearizationsuch that Eq. (0) is reduced to a linear advective dispersion equation(ADE) with respect to a Kirchhoﬀ potential. The ADE can be analyticallysolved for solitary or systematically placed singularities (Kacimov, 2007,
Philip, 1968, 1971, Pullan, 1990, Raats, 1971, 1972). More details on
BVPs for the Laplace equation and the ADE are given in Sections 2-5 of
the MS.
The total head for saturated seepage in undeformable porous media
obeys the Laplace equation and, therefore, the mathematical commonality between ﬂows of ideal ﬂuids (Zhukovskii, 1948) and pore water
motion has been widely explored (PK-62,77, Strack,1989). The methodof images and the theory of BVPs (Henrici, 1993) have been engaged. In
this paper, we combine these traditional analytical and new numerical(HYDRUS-2D, Šimůnek et al., 2016) methods to problems of groundwater and soil water movement from shallow conﬁned and unconﬁned
aquifers towards ditches and topographic depressions, common in desert
landforms of arid-hyperarid climates of Oman and Mars.
A practical motivation of our work stems from a daunting problem of groundwater inundation caused by a rapid rise of the water table in perched unconﬁned and shallow conﬁned aquifers, detected in many urban areas of the world (see, e.g., Attard et al., 2016,
Barron et al., 2013, Chaudhary, 2012, Coda et al., 2019, Howard andIsraﬁlov, 2012, Jha et al., 2012, Kazemi, 2011,Lerner, 2003, Lerner andHarris, 2008, Medovar et al., 2018, Naik et al., 2008, Porse et al.,2016, Preene and Fisher, 2015, Quan et al., 2010, Schirmer et al.,2013, Vázquez-Suñé et al., 2005, Vogwill, 2016). Surprisingly, the arid
countries of Arabia – despite the hydrological mantra on water deﬁcit
and groundwater depletion – encounter the same problem of waterlogging of urban structures by shallow groundwater as the cities inhumid zones (e.g., Abu-Rizaiza et al., 1989, Abu-Rizaiza, 1999, AlRawas and Qamaruddin, 1998, Al-Sefry and Şen, 2006, Al-Senafy, 2011,
Al-Senafy et al., 2015, Alsharhan et al., 2020, Bob et al., 2016,
Kreibich and Thieken, 2008). It is noteworthy that urbanization in theSouthwest of the USA caused similar negative hydrological impacts (ampliﬁcation of ﬂash ﬂood intensities due to reduced evapotranspiration,see, e.g., Goudie, 2013).
In the metropolitan area of Muscat, at the campus of Sultan Qaboos
University (SQU), the groundwater inundation has become evident even
without any piezometry. In recent years, subsurface water seeped out innatural topographic depressions, man-made excavations, ditches, anddrainage channels called wadis. Large wet spots or ponded areas at SQU
are bio-marked by adjacent lush, wild vegetation. In the desert climateof Oman, this vegetation in some places pops up virtually as Maillart’s“grass in green ditches” (see the epigraph). The studied SQU ﬁeld is
located in a hyperarid zone of the Batinah Coast, Oman, Arabia (e.g.,
Alsharhan et al., 2001, 2020). The investigated wadi section was about
1 km long, with several excavated pits-trenches and natural depressions.
For instance, one of the excavations has the following coordinates: Northing: 619721, Easting: 2608729, 47.5 m above the level in the Sea of
Oman. Urbanization of the ﬁeld started in the mid-1980s, before which
the site was a desert with scattered shrubs and trees. Currently, a signiﬁcant part of the campus area is paved, which reduced natural evapotranspiration. Waterlogging and associated lush vegetation are geotechnically deleterious, putting foundations of the SQU buildings, roads,
subsurface cables, etc., in jeopardy, even when the salinity of the rising groundwater is small. Secondary salinization of the topsoil is also
widespread in Oman when the rise of the water table is from an aquiferof poor groundwater quality. Fieldwork was carried out in NovemberJune 2020 at the SQU campus: auger holes, ditches, and excavations
were constructed, and groundwater monitored, to assess soil moistureoozing ﬂuxes above shallow perched and conﬁned aquifers detected in
the study area.
The motivation of this work also stems from our intention to attract
the attention of Planet Mars hydrologists, who model seepage of Mar-

Advances in Water Resources 145 (2020) 103732

tian groundwater into remotely scanned relief features of the Red Planet,
and the legacy of experts in irrigation and drainage engineering, who– during the 20th century - assembled a good arsenal of analytical andnumerical methods for investigations of the motion of pore water towards drains (Bear, 1972, Kirkham and Powers, 1972, Skaggs et al.,1999, 2012, Strack, 1989, Vedernikov, 1939, Zhukovskii, 1948).
Mars is a hyperarid planet where the remotely analyzed landforms(craters, dunes, endorheic lakes and playas, pits, alcoves, gullies,
recurring slope lineae, ﬂuvio-lacustrine basins, fans, shallow dust
horizons, streaks, amphitheater‐headed valleys, sapping valleys,etc.) and hydrological elements (springs, outseeps, deep fracture
conduits, seepage spots), which convey liquid water, brine, andMartian gases, are studied. The alleged paleo and current dynamics
of Martian pore ﬂuids are elaborated, for example, in Abotalib andHeggy (2019), Bhardwaj et al. (2019),Boatwright and Head (2019),
Edwards and Piqueux (2016), Goldspiel and Squyres (2011),
Grimm et al. (2014), Hobbs et al. (2014), Kereszturi et al.,2010, Kochel and Piper (1986), Luo et al. (2011), Malin andCarr (1999) Malin and Edgett (2000), Marra et al. (2014, 2015),
Michalski et al. (2013), Mukherjee et al. (2020), and Salese et al. (2019).
For hypothetical Martian aquifers (both unconﬁned and artesian),
Ks is assessed from proxy data (see, e.g., the posited “nominal” value inGoldspiel and Squyres (2011); Luo et al. (2011) evaluated the hydraulicconductivity from the solution of inverse problems of the hydrologicalcycle on Mars; Boatwright and Head (2019) assumed that the conductivity exponentially decreases with the Martian depth). We note that even
the Darcy law is sometimes not correctly formulated by Martian hydrologists (see, e.g., an erroneous Eq. (1) in Abotalib and Heggy, 2019).
Although purely unsaturated ﬂows are also hypothesized (Edwards andPiqueux, 2016,Grimm et al., 2014) to occur on Mars, to the best of ourknowledge, Martian hydrologists (e.g., Luo et al., 2011) have so far used
only ﬂow models for saturated soils.
Models of Darcian ﬂows towards drainage ditches commonly conceptualize the ditches as rectangles in vertical cross-sections having vertical slopes (Afruzi et al., 2014, Barua and Alam, 2013,
Chahar and Vadodaria, 2008, 2012, Gureghian and Youngs, 1975,
Sarmah and Tiwari, 2018, Youngs, 1975, 1990, 1994). However, natural ephemeral river channels, both on Earth and Mars, as well as
constructed drainage ditches and trenches, often have non-rectangularshapes (Grotzinger et al., 2014, Kocurek et al., 2020). Excavations(pits) and natural depressions are often axisymmetric and mild-sloped(see, e.g., Goldspiel and Squyres, 2011). Only in few studies havephreatic seepage and ﬂows in the unsaturated zone to non-rectangular(trapezoidal and curvilinear) draining entities been analytically examined (PK-77 reported solutions by Bazanov, p. 150 of her book,and Vedernikov, pp.181-182, see also Ilyinsky and Kacimov, 1992a,
Kacimov, 2005, 2006a, Kacimov and Obnosov, 2002). We note that although Luo et al. (2011) sketched trapezoidal draining channels in theirmodels of groundwater discharge into these channels, they posited the
Dupuit-Forchheimer (hereafter abbreviated as DF) model, which actually ignores the shape of channels. Marra et al. (2014, 2015) conceptualized, and studied experimentally in sandboxes, seepage into real trapezoidal channels, allegedly occurring on Mars, but unfortunately did notrefer to a plethora of theoretical studies of this type of seepage to draining channels on our own planet (see, e.g., Aravin and Numerov, 1953,
PK-62,77, Vedernikov, 1939).
This work is organized as follows. In Section 2, we develop a simple analytical solution for a 2-D (not DF!) conﬁned ﬂow to a curvilinear
ditch draining a saturated, conﬁned aquifer. Some natural wadi channels
in our study area are arcuate, like ones on Mars (Malin and Carr, 1999),
but as the ﬁrst approximation, we ignore the 3-D eﬀect. Our ditch contour is an isobar generated by a line sink. In Section 3, we obtain a similar solution for a point sink in an axisymmetric ﬂow towards a cratershaped isobar, the piezometric head being still a harmonic function. In
Section 4, we assemble an inﬁnite array of J.R. Philip’s sinks-sourcesof the advective dispersion equation (ADE) for the Kirchhoﬀ potential,

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

which models an axisymmetric, steady water movement from a horizontal water table towards a topographic depression in a partially saturated
so-called Gardner soil. In Section 5 and Appendix I, 2-D seepage towards
a triangular empty or partially-ﬁlled ditch is studied by the machineryof holomorphic functions. All analytical solutions of Sections 2-5 aretested against HYDRUS-2D numerical simulations. In Section 6, we outline applications to arid zone hydrology, and in particular, to Martianhydrology.
2. A line-sink-generated seepage face ditch draining a conﬁned
aquifer
Unlike Hobbs et al. (2014), who could not identify conﬁning layers,
in the form of bedrocks and caprocks, of hydrostratigraphic units in theMartian regolith massif, we - when making excavations on our SQU site
- indeed crashed a low-permeable caliche layer at a depth of 1.5-2 m.
Typical values of Ks for low-permeable (cemented) layers (also called
calcrete, petrocalcic horizons, see, e.g., Duniway et al., 2010), formed by
the precipitation of carbonates or gypsum due to intensive evaporationin semi-arid, arid and hyper-arid regions, are given in the literature:10−4 to 10−6 cm/s (Al-Yaqout and Townsend, 2004), 0.0002 to 0.008
m/day (Heilweil and Watt, 2011), 2.8∗ 10−7 m/s (Mujica and Bea, 2020),
among others. In most reported cases, these strata hydrologically cap thevadose zone rather than a perched or conﬁned aquifer in our situation.
We could not excavate deeper to the next low-permeable horizon,
and therefore, in Fig. 1a, the aquifer is not conﬁned from below. We
note that genesis of the caliche in the SQU soil proﬁle is similar towhat Michalski et al. (2013) argued about the impact of the upwelling
Martian groundwater on the sequence of cemented sediments of the
Red Planet. Speciﬁcally, a periodic “upwelling” (rise) of groundwaterin SQU shallow perched unconﬁned aquifers, and intensive evaporation

Advances in Water Resources 145 (2020) 103732

from the water table has – allegedly – formed the caliche and “selfconﬁned” this aquifer by a caprock (Fig. 1a). When we “punctured” thecaliche by the JCB scoop, groundwater gushed into the excavation andﬁlled it to the level of ℎ∗𝐴 (Fig. 1a) within 2 hours. We collected soilsediment samples, measured Ks using standard techniques (see, e.g., AlShukaily et al., 2020b), and ascribed the obtained value to the texturalclass of HYDRUS-2D sands.
Fig. 1a shows a vertical cross-section of an empty ditch B0 MC0 ,which drains a dry soil layer above a horizontal aquifuge and a conﬁned,
homogeneous and isotropic aquifer (having at saturation a hydraulicconductivity Ks ) beneath an originally impermeable layer (caprock)
ABCA. Point A in Fig. 1a is at inﬁnity on the Riemann sphere. If the
excavation depth H∗ , counted from ABCA (point M is the deepest pointof the ditch), is below the aquifuge, groundwater discharges into theditch BMC (symmetric and having the width 2 L) with the total ﬂow
rate 𝑄∗𝑑 (m2 /s).
We assume that ﬂow takes place only below the aquifuge.
Goldspiel and Squyres (2011) call this layer “aquiclude”, that is – inthe vernacular of groundwater hydrology (Strack, 1989) – a “Martian
misnomer,” because Goldspiel and Squyres (2011) model impermeable
rather than leaky conﬁning layers. At a point Ah (Fig. 1a), located a certain distance 𝐿∗𝐴 (𝐿∗𝐴 >> 𝐿)| from O, a remote piezometer shows an elevation ℎ∗𝐴 . For simplicity, we ignore groundwater exﬁltration upward
through the caliche in Fig. 1a, i.e., leaky layer (aquitard) ﬂow scenarios,
as considered by Kacimov and Obnosov (2008, 2019).
We introduce Cartesian coordinates Ox∗ y∗ and the complex physicalcoordinate z∗ =x∗ + i y∗ (Fig. 1b). In this section, we study the case of anempty ditch, i.e., the pressure head, p∗ , along BMC is zero. The emptiness of the ditch can be ensured by either a high Manning slope in the
direction perpendicular to the plane of Fig. 1a, so that all seeped water
rapidly ﬂows away as surface water (see, e.g., Al-Shukaili et al., 2020a),Fig. 1. A vertical cross-section of an empty ditch B0 MC0 a); the right half of the empty ditch constructed such that the bottom of the ditch, BMC, is a seepage faceof the ﬂow in the domain Gz , generated by a ﬁctive line sink at the origin O(x,y) b); the complex potential domain Gw for the whole ditch c).

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 2. a) Half-contours of sink-generated ditches for Qd = 0.2, 1, and 5.86 (curves 1, 2, and 3, respectively). b) The depth H of a sink-generated ditch and its area
A (curves 1 and 2, respectively) as functions of the seepage ﬂow rate Qd .

or by intensive evaporation, provided H∗ is relatively small and climate
is arid enough. We note that surface water ﬂow perpendicular to theplane of Fig. 1, i.e., along the ditch, can be easily tracked on Earth. In ourﬁeldwork, we detected what Malin and Edgett (2000) called a Martian“seepage-fed surface runoﬀ,” slow Hortonian motion downslope a ditch,with velocities of few mm/s. On Mars, the juxtaposition of groundwaterand surface water ﬂows, reconstructed from the contemporary Martian
landforms, is purely hypothetical (Malin and Carr, 1999).
As we have already pointed out in Section 1, a saturated 2-D ﬂow
from a homogeneous aquifer into a ditch (Fig. 1a) obeys Laplace’s equation:Δℎ∗ (𝑥∗ , 𝑦∗ ) = 0,(1)⃖⃖⃗∗ (𝑢∗ , 𝑣∗ ) = −𝐾𝑠 ∇ℎ∗ is the Darcian velocity vector, u∗ and v∗ arewhere 𝑉its horizontal and vertical components, h∗ is the total piezometric head,h∗ = p∗ - y∗ , where the pore water pressure is P∗ = 𝜌d gp∗ , g is gravity acceleration (9.8 m/s2 on Earth and3.7 m/s2 on Mars), and 𝜌d ispore ﬂuid’s (groundwater, brine, soil moisture) density. We emphasizethat Eq. (1) and all the results below do not assume a quasi-horizontal
pore water motion, as it is postulated, for instance, in the DF model of
Boatwright and Head (2019) and Luo et al. (2011).
We introduce the complex potential w∗ = 𝜑∗ + i𝜓 ∗ , 𝜑∗ = −Ks h∗ . The
stream function 𝜓 ∗ and a velocity potential 𝜑∗ are conjugate harmonic
functions. Due to symmetry, we consider only the right half of the physical domain Gz (Fig. 1b). The following dimensionless quantities are
introduced:(𝐻, 𝐿𝐴 , 𝑥, 𝑦, 𝑧, 𝑟, ℎ, ℎ𝐴 , 𝑝) = (𝐻 ∗ , 𝐿∗𝐴 , 𝑥∗ , 𝑦∗ , 𝑧∗ , 𝑟∗ , ℎ∗ , ℎ∗𝐴 , 𝑝∗ )∕𝐿, (𝑤, 𝑄𝑑 )= (𝑤∗ , 𝑄∗𝑑 )∕(𝐾𝑠 𝐿), (𝑉 , 𝑢, 𝑣) = (𝑉 ∗ , 𝑢∗ , 𝑣∗ )∕𝐾𝑠
where r(0 ≤ r < ∞) is the modulus of point z=x+i y in the aquifer(Fig. 1a).
The complex potential domain Gw is depicted in Fig. 1c. The piezometric head is assumed to be zero at points B and C. The image of the
BMC isobar in G is a speciﬁc curve whose shape is unknown. Far awayfrom the ditch, the piezometric head is constant (hA ) along a circle ofa radius LA in Gz (Fig.1a). The image of this circle is a segment Ah Av in
Gw (Fig. 1c).
Now we use the trick implemented by Zhukovskii (1948) in ﬂuid
mechanics of ideal ﬂuids (see recent applications in Belotserkovsky and
Lifanov, 1992), also employed in problems of subsurface mechanics of steady and transient Darcian ﬂows (Fujii and Kacimov, 1998,
Kacimov, 1992, 2000a, 2007, 2009, Kacimov et al., 2009, Kacimov andObnosov, 2002, PK-62,77, Strack, 1989, 2020, Vedernikov, 1939). We
place a line sink of intensity Q=2Qd at the origin of coordinates. The

complex potential of this sink is:
𝑄

𝑄

ℎ=

𝑄𝑑
ln𝑧,
𝜋

𝑤 = −i 2𝑑 − 𝜋𝑑 ln 𝑧,

𝑧 = 𝑟 exp [i𝜈],𝑄𝑄𝜓 = − 2𝑑 − 𝜋𝑑 v,

𝑟=√𝑥2 + 𝑦2 ,𝑄|𝑉 | = 𝜋𝑟𝑑 ,(2)

where 𝜈 is the argument of z. It is evident from Eq. (2) that the streamlines of this ﬂow are rays converging to point O in Fig. 1b. The isotachs
and equipotential lines are semicircles centered at the same point.
Now we restrict r in Eq. (2) to only a positive-pressure part of thehalf-plane in Fig.1a, i.e., we use Eq. (2) for the reconstruction of BMC
by plotting the curve, which is determined by the following explicit formula:√
𝑄𝑑 ln 𝑥2 + 𝑦2 = 𝜋𝑦, 0 ≤ 𝑥 ≤ 1, −𝐻 ≤ 𝑦 ≤ 0.
(3)
The depth H of the ditch is found from the solution of the equation
𝑄𝑑 Log|𝑦| = 𝜋𝑦,(4)

and is equal to
𝐻 = −𝑄𝑑 𝑊 (−𝜋∕𝑄𝑑 )∕𝜋,(5)where W stands for the Lambert W-function (implemented as ProductLog in Wolfram’s (1991) Mathematica).
Fig. 2a shows the MC contours for sink-generated ditches (Fig. 1b)
plotted for Qd = 0.2, 1, and 5.86 (curves 1, 2, and 3, respectively). Curve
1 in Fig. 2b shows H(Qd ). From Eq. (3) it is obvious that for small Qd ,the ditch degenerates into a horizontal segment, while in another limit,lim 𝐻 = 1, BMC becomes a semicircle. The analytical solution (2)-(5)𝑄𝑑 →∞is similar to those obtained by Bazanov (see PK-62, 77) and Kacimov andObnosov (2002) for empty ditches reconstructed by the speciﬁcation ofboundaries corresponding to BMC in the hodograph or other auxiliary
planes (e.g., the plane of the Zhukovskii holomorphic function, which
is deﬁned as w–i z, see PK-62,77). We also recall the analytical andnumerical solutions of Dagan (1964) and Gjerde and Tyvand (1992),
who studied potential transient 2-D (not DF) ﬂows towards horizontal
drains (modeled by linear sinks and empty circles, correspondingly).
Curve 2 in Fig. 2b shows A(Qd ). The cross-sectional area, A, of BMC(in Fig. 1a) is evaluated by the NIntegrate routine of Mathematica as:
𝐴 = −21∫0𝑦(𝑥)d𝑥,(6)

Fig. 3 shows the results of modeling with HYDRUS-2D Šimůnek et al.,2016). We selected Qd = 0.2/𝜋 in the analytical solution (3). With thehelp of Eqs. (3) and ((5) we generated 10 points on MC, imported themfrom Mathematica into HYDRUS, and made a spline curve in dimensional
quantities used by HYDRUS. Also, HYDRUS uses z for the vertical coordinate. Fig. 3a shows the HYDRUS modeling domain bounded by an

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 3. The HYDRUS ﬂow domain for a half-ditch generated by a line sink a), isotachs b), and streamlines c).
Fig. 4. An axial cross-section of an axisymmetric crater draining a conﬁned aquifer.

isobaric curve MC, for which we selected L=100 cm. Soil is loam fromthe HYDRUS soil catalog. The horizontal and vertical segments CAh andMAv in Fig. 3 are no-ﬂow lines. We selected 𝐿∗𝐴 = 500 cm, i.e., we assumed a quarter-circle Ah Av with a radius of 500 cm to be a line with aconstant piezometric head (we recall that HYDRUS, unlike PK-62,77 andStrack, 1989, uses h as a notation for the pressure head), along which
we set ℎ∗𝐴 = 32.2 cm. This value was found at r∗ =500 cm for a selectedvalue of 𝑄∗𝑑 from the logarithmic variation of h∗ (x∗ ) along CA, counting from the ﬁducial point C where ℎ∗𝐶 = 0, see Eq. (2).Fig. 3b shows
the isotachs (also constant h∗ and 𝜑∗ contours), and Fig. 3c shows thestreamlines. The origin of the spatial coordinates is retained in Fig. 3 toindicate the position of the generating sink. HYDRUS results match very
well the analytical solution.

3. A point-sink-generated dimple draining a conﬁned aquifer
In this section, we use the same procedure as in Section 2 for thecase of axisymmetric ﬂow towards a depression BMC shown in Fig. 4
in an axial cross-section. We select a system of cylindrical coordinatesOz∗ 𝜌∗ 𝜅, wherez∗ is a vertical√√ coordinate (the same now as in HYDRUS∗∗22D),𝜌 = 𝑥 + 𝑦∗2 , 𝑟∗ = 𝜌∗2 + 𝑧∗2 and the angular coordinate 𝜅vanishes from the analysis due to axial symmetry. At the origin O, we
now place a point sink of intensity 𝑄∗𝑑 (m3 /s).
The piezometric head caused by a point sink is (PK-62,77):ℎ∗ =)𝑄∗𝑑 ( 11−,
2𝜋𝐾𝑠 𝑅 𝑟∗(7)

where R is the radius of the depression that coincides with the r∗ (and𝜌∗ ) coordinate of points B and C in Fig. 4.Eq. (7) sets up the total headequal to zero in these reference points. The depression of a depth H∗ isempty and, therefore, the pressure head 𝑝∗𝐵𝑀𝐶 = 0 along this curve in
Fig. 4 (a surface in 3-D). At a remote point Ah , located at a distance 𝑅∗𝐴(𝑅∗𝐴 >> 𝑅) from point O (Fig. 4), we have:
ℎ∗𝐴 =

𝑄∗𝑑

2𝜋𝐾𝑠 𝑅.

(8)

Point Ah is a surrogate inﬁnity of the inﬁnite point A on the Riemannsphere. The head of Eq. (8) is illustrated in a remote piezometer sketchedin Fig. 4. We recall (PK-62,77) that for ﬂow to a 3-D sink, the piezometrichead at inﬁnity is ﬁnite, unlike for a 2-D sink (Section 2), for which thehead at inﬁnity (point A) in Fig. 1b logarithmically blows up to inﬁnity.
We introduce dimensionless quantities: (H, z, 𝜌, r, h, p,
hA ,RA )=(𝐻 ∗ , 𝑧∗ , 𝜌∗ , 𝑟∗ , ℎ∗ , 𝑝∗ , ℎ∗𝐴 , 𝑅∗𝐴 ),𝑄𝑑 = 𝑄∗𝑑 ∕(2𝜋𝐾𝑠 𝑅2 ),(𝑉 , 𝑢, 𝑣) =(𝑉 ∗ , 𝑢∗𝜌 , 𝑣∗𝑍 )∕𝐾𝑠 , where (u,v) are now the radial and vertical velocity
components.
From Eq. (7), a zero-pressure isobar (seepage face) BMC is reconstructed, similarly to Eq. (3), by an equation:√1
𝜌=− 𝑧2 , 0 ≤ 𝜌 ≤ 1, −𝐻 ≤ 𝑧 ≤ 0.
(9)
(1 − 𝑧∕𝑄𝑑 )2
We note that some Martian hydrologists use the term “aquifer sapping face” and “crater wall” (Goldspiel and Squyres, 2011) for what iscalled a “seepage face” in terrestrial civil engineering (Strack, 1989).
Eq. (9) deﬁnes a quartic line. The depth H is found from Eq. (9),which at 𝜌 = 0 is reduced to a quadratic equation having the following

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 5. a) Half-contours of a sink-generated axisymmetric topographic depression draining a conﬁned aquifer for Qd = 0.2, 1, and 2 (curves 1, 2, and 3, respectively).
b) Depth and volume (curves 1 and 2, respectively) of an axisymmetric sink-generated depression in a conﬁned aquifer.

Fig. 6. HYDRUS isobars corresponding to Qd = 0.2 a), isotachs b), and streamlines c).

4. An isobaric depression reconstructed from Philip’s point-sinkarray and HYDRUS-2D simulations

physically meaningful solution:

𝐻=√
𝑄2𝑑 + 4𝑄𝑑 − 𝑄𝑑
2

.

(10)

We used Eqs. (9-10) and calculated the volume of the depression
BMC in Fig. 4 as a body of revolution:
𝑉 =𝜋

0

∫−𝐻

𝜌2 (𝑧)d𝑧 =)√𝜋𝑄𝑑 ( 2
𝑄𝑑 + 6𝑄𝑑 + 6 − (4 + 𝑄𝑑 ) 𝑄𝑑 (4 + 𝑄𝑑 ) ,6(11)

Fig. 5a shows the isobars MC for Qd = 0.2, 1, 2 (curves 1-3, correspondingly). In Fig. 5b, the curves (10) and (11) are plotted. Obviously,
for lim 𝐻 = 1 and lim 𝑉 = 2𝜋∕3, BMC becomes a hemisphere in this
𝑄𝑑 →0

𝑄𝐷 →∞

limit, while for small Qd , we get a disk.
The results of HYDRUS-2D simulations (in dimensional quantities)
corresponding to the analytical case of Qd = 0.2 are presented in Fig. 6.
As in Section 2 Fig. 3), we considered a loamy soil. We selected theradius R of the HYDRUS axisymmetric depression to be 100 cm. Using
Eqs. (9) and ((10), we generated ten discrete points on MC and made a
HYDRUS spline (seepage face). The radius 𝑅∗𝐴 of a constant piezometrichead quarter-circle Ah Av was 500 cm. According to Eq. (8) the total headℎ∗𝐴 = 𝑝∗𝐴ℎ = 20 cm at point Ah and, therefore, along Ah Av , i.e., the HYDRUS pressure head is 𝑝∗𝐴𝑉 =520 cm at point Av . Fig. 6a presents steadystate pressure head isolines. Fig. 6b illustrated isotachs and Fig. 6c plots
the streamlines. Fig. 6bc perfectly ﬁt the analytical solution.
In the next section, we focus on an opposite case of steady ﬂow in apurely unsaturated Gardner soil, albeit the same trick will be used, viz.
we place a ﬁctitious sink “in the air” and use real isobars to constructtopographic craters (dimples).

Analytical solutions of Sections 2 and 3 ignored ﬂow in the unsaturated zone and the capillary fringe, which are important in the caseof small-scale drainage entities (e.g., Abit et al., 2008, Silliman et al.,2002). In this section, we incorporate into an analytical quasilinear
model all three physical factors, which control seepage: the gravitational
force (on Earth or Mars), Darcian resistance of the soil, and capillarity.
We use a constitutive relation for the phase conductivity:[]
𝑘𝑢𝑛 = 𝐾𝑠 exp 𝛼𝑝 𝑝∗ ,(12)where kun (p∗ ) for p∗ ≤ 0 is an unsaturated hydraulic conductivity, 𝛼 p(const > 0) is the sorptive number (1/m) (Gardner, 1958, see also Raats,1973). J.R.Philip (1968) pioneered utilizing Eq. (12) for modeling 2Dand 3D unsaturated ﬂows (see, e.g., Pullan, 1990, Raats et al., 2002).
The soil constant 𝛼 p is related to the VG-HYDRUS parameters 𝛼 andn via well-known empiric relations (Acharya et al., 2012); we used
Eq. (10) from Ghezzehei et al. (2007).
The evaporation arises from a steady upward ﬂow from a horizontal(shallow) water table W1 W0 W2 to a hot and dry soil surface BMC (Fig. 7)
where the pressure head is constant (𝑝∗ = −𝑝∗𝑠 ,𝑝∗𝑠 > 0). This surface has
a topographic depression with the deepest point M at a depth H∗ > 0
counted from the soil surface, which is horizontal and ﬂat, far awayfrom the topographic trough in Fig. 7, i.e., the water table is at a given
depth b∗ above a ﬂat soil surface (near remote points B and C). BMC
in Fig. 7 can also be related to the “sapping valleys” on Mars (see, e.g.,
Salese et al., 2019), albeit there is a diﬀerence in scales: the Martian
depth of the water table (𝑍𝑠∗ − 𝑍1∗ ) is 4-5 km, while at SQU, the shallow
groundwater table is only several tens of cm deep.

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 7. An axial cross-section of evaporative
ﬂow in an unsaturated soil from a horizontal
water table W1 W0 W2 to a topographically depressed isobar BMC.

The depression in Fig. 7 is axisymmetric, like in Section 3. The axisof symmetry coincides with the OZ∗ axis of a cylindrical system of coordinates OZ∗ 𝜌∗ (introduced similarly as in Section 3 above). To be consistent with Philip, we orient OZ∗ downward. Thus, the unsaturated soilin 3-D is sandwiched between the W1 W0 W2 plane and the soil surface
obtained by the revolution of BMC with respect to the axis OZ∗ in Fig. 7.
There are no physical singularities representing drains, subsurface emitters, or plant roots in this soil layer.
The contour BMC in Fig. 7 is obtained in the following manner(Philip, 1989, Obnosov and Kacimov, 2018). A ﬁctitious point sink ofintensity 𝑄∗𝑝 (m3 /s) is placed at the origin of coordinates, such that thewater table is beneath this sink at a depth of 𝑍1∗ > 0. A ﬁctitious horizontal isobaric plane F1 F0 F2 is placed at the level 𝑍 ∗ = −𝑍2∗ ,𝑍2∗ > 0,
as illustrated in Fig. 7. At this plane, the pressure head 𝑝∗ = −𝑝∗0 is
speciﬁed. Points B and C in Fig. 7 of an “intermediate isobar” between two sandwiching isobaric planes F1 F0 F2 and W1 W0 W2 are at thelevel of 𝑍𝑠∗ , which can be either positive or negative. Fictitious point
sources and sinks of intensity 𝑄∗𝑝 mirror each other with respect tothe planes F1 F0 F2 and W1 W0 W2 . Only the ﬁrst ﬁctitious source under
the water table, placed on the OZ∗ axis at a depth of 2𝑍1∗ , is shownin Fig. 7.The array of an inﬁnite number of such sinks-sources generates a family of isobars, some of which are real (two dashed lines𝑝∗ = −𝑝∗2 ,𝑝∗ = −𝑝∗1 ,𝑝∗2 < 𝑝∗1 < 𝑝∗𝑠 are sketched in Fig. 7), and others are ﬁctitious (two dotted lines 𝑝∗ = −𝑃2∗ ,𝑝∗ = −𝑃1∗ ,𝑝∗𝑠 < 𝑃2∗ < 𝑃1∗ are sketchedin Fig. 7). Mathematically, instead of BMC in Fig. 7, we can select another soil surface, for example, any of the two dashed isobars.
Practically, we proceed in the following manner. One goes to theﬁeld with a theta-probe and collects the moisture content, 𝜃 s , from thetopsoil of a real depression, having certain geometrical sizes b∗ , H∗ . If
𝜃 s (and hence 𝑝∗𝑠 ) is almost constant in the depression and on a ﬂat surface away from the depression, then the position (𝑍1∗ , 𝑍2∗ ) of the originof coordinates (sink’s locus) and the strength 𝑄∗𝑝 of Zhukovskii’s sinkssources array should be mathematically adjusted in such a manner thatthe physical (ﬁeld) depression becomes close to a mathematical isobar
BMC in Fig. 7. In what follows, we illustrate this algorithm.

Routinely (see Philip, 1968, Raats, 1970, Obnosov and Kacimov, 2018), we introduce the Kirchhoﬀ potential, Ω∗ , which is often
called the matric ﬂux potential:Ω∗ ( 𝜌 ∗ , 𝑍 ∗ ) =

𝑝∗

∫−∞

[
]
𝑘𝑢𝑛 (𝑢)d𝑢 = 𝑘𝑢𝑛 𝑝∗ (𝜌∗ , 𝑍 ∗ ) ∕𝛼 𝑃 ,Ω∗∞ = 𝐾𝑠 ∕𝛼 𝑃 ,(13)i.e., Ω∗∞ is the Kirchhoﬀ potential along W1 W0 W2 . Obviously, the lines
of a constant Kirchhoﬀ potential are also isobars (the values of corresponding pressure heads are found directly from Eq. (13)) and isowetness curves (the corresponding water content is found from Eq. (13) andthe soil-water retention function).
We introduce dimensionless quantities: (𝑧, 𝑏, 𝐻, 𝑧𝑠 , 𝑧1 , 𝑧2 , 𝜌, 𝑟, ℎ, 𝑝)=(𝑍 ∗ , 𝑏∗ , 𝐻 ∗ , 𝑍𝑠∗ , 𝑍1∗ , 𝑍2∗ , 𝜌∗ , 𝑟∗ , ℎ∗ , 𝑝∗ ) × 𝛼∕2,Ω=Ω ∗ / Ω∞ , 𝑄 𝑝 =∗∗∗∗𝑄𝑝 𝛼∕Ω∞ ,(𝑉 , 𝑢, 𝑣) = (𝑉 , 𝑢𝑝 , 𝑣𝑍 )∕𝐾𝑠 , for which the Richards Eq. (0) is
reduced to a linear ADE:
𝜕 Ω∗
𝜕 Ω∗
𝜕 Ω∗
ΔΩ∗ = 𝛼𝑃, 𝑢∗𝜌 = − ∗ , 𝑣∗𝑍 = 𝛼𝑃 Ω∗ −
,
(14)
∗
𝜕𝑍
𝜕𝜌
𝜕𝑍 ∗
where Δ is the Laplacian operator in (𝜌,Z), and 𝑢∗𝜌 and 𝑣∗𝑍 are the radial
and vertical velocity components. Then, we adapt Philip’s (1989) andObnosov and Kacimov (2018) solutions to the case of ﬂow in Fig. 7 andobtain a series:
[
[
]
[
]]
[
]
Ω0 exp 2(𝑧1 + 𝑧2 ) − exp 2(𝑧 + 𝑧2 ) + exp 2(𝑧 + 𝑧2 ) −1 𝑄𝑃 exp 𝑧
Ω=
−
[
]
8𝜋
exp 2(𝑧1 + 𝑧2 ) −1[]
exp(−𝑟) exp(−𝑟1 ) exp(−𝑟2 ) exp(−𝑟3 ) exp(−𝑟4 )
−
−
+
+
− …,
×
𝑟
𝑟1
𝑟2
𝑟3
𝑟4
√
√
√
𝑟 = 𝜌2 + 𝑧2 , 𝑟1 = 𝜌2 + (𝑧 − 2𝑧1 )2 , 𝑟2 = 𝜌2 + (𝑧 + 2𝑧2 )2 ,√√𝑟3 = 𝜌2 + (𝑧 + 2𝑧1 + 2𝑧2 )2 , 𝑟4 = 𝜌2 + (𝑧 − 2𝑧1 − 2𝑧2 )2 , …
(15)
where the constant Ω0 = exp [ − 𝛼 p z2 ],0 ≤ Ω0 ≤ 1 is the Kirchhoﬀ potential along F1 F0 F2 . At points B and C, which are located on a physicalisobar, whose remote «wings» are above the sink in Fig. 7, we have
Ω = Ωs = exp [ − 𝛼 p z2 ] for Ω0 ≤ Ωs ≤ 1.
The ﬁrst term on the RHS of Eq. (15) determines 1-D unsaturated
evaporative ﬂow from W1 W0 W2 to F1 F0 F2 in the layer z1 < z < −z2 .

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Fig. 8. Analytically computed isolines of Kirchhoﬀ’s potential for (z1 , z2 ,Ω0 ,QP )=(1, 1, 0.2, 1).

The terms in the square brackets of Eq. (15) are responsible for sinks
and sources, and consequently, for mirrored (imaged) with respect tothe plains z = z1 and z = −z2 . We emphasize that Philip (1989) andObnosov and Kacimov (2018) studied ﬂow for one real source representing a subsurface irrigation emitter, while all in Fig. 7 and Eq. (15) are
ﬁctitious. A trick similar to ours was used by Zhukovskii (1948) in his famous formula for aerodynamics of ideal ﬂuids. He combined mathematical singularities of a characteristic holomorphic function (two dipolesand vortex), the latter placed inside an airfoil (e.g., an impermeablecylinder). Zhukovskii also ignored a ﬁctitious ﬂow inside an airfoil (inthe vicinity of the vortex) and only considered an exterior of a mathematical separatrix-streamline (cylinder). We conﬁne real Darcian ﬂow
not by streamlines (for the sake of brevity, we have not even introduced
the Stokes stream function, see Obnosov and Kacimov, 2018), but by
two isobars (BMC and W1 W0 W2 in Fig. 7).
The results of computations based on Eq. (15), where the series wastruncated to 5 terms, are plotted by the ContourPlot routine of Mathematica in Fig. 8 for (z1 , z2 ,Ω0 , QP ) =(1, 1, 0.2, 1). As an example, inFig. 8, we selected an isobar BM as the real depressed isobaric soil surface maintained at the Kirchhoﬀ potential Ωs = 0.26. For this isobar, weused the FindRoot routine of Mathematica and evaluated zM = 0.389,
zB = zs = −0.177, i.e., the dimensionless depth of the depression in
Fig. 8 is H=0.566.
HYDRUS-2D simulations are presented in Fig. 9. We used silt loamwith the HYDRUS triad of the VG soil hydraulic parameters (𝛼, n, Ks )=(21/m, 1.41, 1.08 m/day). We followed Ghezzehei et al. (2007) andmade Gardners’s soil having 𝛼 P =1.3𝛼 n, i.e., 𝛼 P =3.67 1/m. In Fig. 9a,a HYDRUS axial section is shown for the case of the pressure head𝑝∗𝑀𝐶 = −0.82 m. The ﬂow domain, which corresponds to the analyticaldepression, has the HYDRUS coordinates zM = −0.13 m, zC = 0.23 m,
and zW2 = −0.55 m. The analytical isobar MC is converted to the HYDRUS spline curve similarly as in Sections 2 and 3. Fig. 9b and c show the
HYDRUS isobars and streamlines. As is evident from Fig. 9c, the ditch
funnels up the ﬂux of moisture from the water table that also evaporates
at the almost ﬂat part of the soil surface (near point C). In other words,
a simple ﬂow tube is realized in Fig. 9c.
In Fig. 9d, we assumed a wetter soil surface with 𝑝∗𝑀𝐶 = −0.63 m.
The origin of coordinates O is again, like in Figs. 8 and 9a, at the locusof the analytic imaginary sink. The geometry of the ﬂow domain has

Advances in Water Resources 145 (2020) 103732

changed as compared to Fig. 9a. The corresponding analytical depression MC has now the HYDRUS coordinates zM = −0.1 m,zC = 0.176 m,
and zW2 = −0.55 m. The isobars are shown in Fig. 9d. They seem uninteresting, like ones in Fig. 9b. Fig. 9e presents HYDRUS streamlines,
which are more interesting. Indeed, the ﬂow topology is diﬀerent from
one in Fig. 9c. Speciﬁcally, like in Kacimov and Youngs (2005), who
studied a seemingly trivial ﬂow domain having a nontrivial two-sheet
Riemann surface in the hodograph plane (see also Anderson, 2013). In
Fig. 9e, we see a separatrix Ss Sw Sd , which divides the unsaturated domain into three subdomains. In the ﬁrst one (counted from the left in
Fig. 9e), moisture ascends from segment W0 Sw of the water table tothe closest (deepest) section MSd of the ditch surface (like in Fig. 9c).
At the rest of the water table, segment Sw W2 , moisture descends (inﬁltrates) from Ss C on the soil surface (third subdomain). PointSw is astagnation point. Stagnation points and separatrix streamlines also appeared in Raats (1977) for steady ﬂows to an array of the parallel line
sinks in unsaturated Gardner soils. In HYDRUS, the numerical value ofthe velocity vector at this point (x=0.35 m) attains a sharp minimum(1.3 × 10−5 m/day).
The second, and most interesting, subdomain in Fig. 9e is boundedby the separatrix (streamline) and a curve segment Sd Hp Ss . Along thisisobar (we recall that 𝑝∗𝑀𝐶 = −0.63 m), moisture inﬁltrates into the soil
through Hp Ss and exﬁltrates back into the atmosphere through Sd Hp ,making a U-turn. Therefore, the whole ﬂow domain in Fig. 9e represents
a complex ﬂow tube similar to one in the potential ﬂow shown in Fig.
62 of PK-77. It is noteworthy that in potential ﬂows governed by the
Laplace equation, this type of separatrix (complex ﬂow tubes) emerges
in situations when the boundary of the ﬂow domain has three or moreequipotential components. Fig. 9c and e illustrate that in ADE-governed
unsaturated ﬂow, a complex topology can appear in a domain boundedby two lines (surfaces) of a constant Kirchhoﬀ potential.
A hinge point Hp in Fig. 9e is similar to ones in the classical Tothian (see, e.g., Kirkham, 1947, Toth, 2009) topology of purely saturated
ﬂows controlled by gravity, Darcian resistance, and topography of theisobaric ground surfaces (postulated by Toth), along which groundwaterinﬁltrates-exﬁltrates (but without taking into account unsaturated ﬂow
in the vadose zone). At the point Hp ,the Darcian velocity vector is tangential to the soil slope. We rephrase: if the velocity vector is decoupled
into the components normal and tangential to the local soil surface, thenat Hp , the normal component is zero. Along Sd Hp , this normal component is oriented from the soil to the air, whereas along Hp Ss , the normalcomponent points into the soil.
To summarize: in the analysis of ﬂow in an unsaturated soil represented in Fig. 9d and e, we amended the Tothian (2009) physical factors by capillarity of the soil. We then got a topological pattern shown(upside-down) at the forefront of the Freeze and Cherry (1979) and ina normal way in Figs. 2.25 – 2.26 of Radcliﬀe and Šimůnek (2010).
5. A triangular ditch in a conﬁned aquifer
In this section, we return to steady, and capillarity-free ﬂow in saturated soil. We consider 2-D potential ﬂow in a vertical plane (Fig. 10a)and use the notations of Section 2, i.e., the complex physical coordinate
is again z∗ =x∗ +i y∗ .
A caprock of a conﬁned aquifer is now breached by a triangularisosceles ditch B0 MC0 , having a bank slope 𝜋𝜔, 0 ≤ 𝜔 ≤ 1/2. The ﬂow
domain, Gz , is now the half-plane y∗ < 0 with a triangular indent BMC.
Ditch’s depth is H∗ = L∗ tan 𝜋𝜔.
Unlike in Section 2, we allow water to accumulate in the ditch up thehorizon B1 C1 at the level 𝐻𝑤∗ . The piezometric head at inﬁnity (point A
in Fig. 10ab) is inﬁnite but, as in Section 2, we replace this inﬁnity by a“surrogate inﬁnity,” viz. an equipotential line Av Ah (shown as a dashedsemicircle and a segment in Fig. 10a and b, respectively). Practically, ina remote piezometer located distance 𝐿∗𝐴 from the ditch axis, the total
head, ℎ∗𝐴 , counted from points B and C, is measured. The overall ﬂow
rate into the ditch is 𝑄∗𝑡 . The domain Gw of the complex potential w is

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 9. HYDRUS simulations of purely unsaturated ﬂow between a water table and an analytically found isobaric depression: a) ﬁnite element discretization of
the ﬂow domain and the boundary conditions, b) and d) pressure head contours, c) and e) streamlines. b) and c) ﬂow to a relatively dry soil surface, d) and e)
evaporation-inﬁltration at a relatively wet soil surface.

sketched in Fig. 10b. The boundary conditions for w are⎧ ∗∗⎪𝜙 = −𝐾 𝑦 along 𝐵 𝐵1 and 𝐶 𝐶1 ,⎪𝜙∗ = 𝐾𝐻𝑤∗ along 𝐵1 𝑀 and 𝑀 𝐶1 ,⎨𝜓 ∗ = 𝑄∗ ∕2, along 𝐴𝐵,𝑡⎪⎪𝜓 ∗ = −𝑄∗𝑡 ∕2, along 𝐴𝐶.
⎩(16)

To the dimensionless quantities introduced in Section 2, we add(𝐻𝑤′ 𝐿𝐴 , ℎ𝐴 ) = (𝐻𝑤∗ , 𝐿∗𝐴 , ℎ∗𝐴 )∕𝐿, 𝑄𝑡 = 𝑄∗𝑡 ∕(𝐾𝑠 𝐿).
We map the tetragon Gz onto a ha(a.lf-plane 𝜂 ≥ 0 of a reference
plane 𝜁 = 𝜉 + i𝜂 (Fig. 10c) with the correspondence of points M → 0,
A → ∞, B → 1, C → −1, B1 → 𝛽, C1 → −𝛽 . Appendix I presents the detailsof our solution. An alternative method for the analytical solution (seethe taxonomy of these methods in Aravin and Numerov, 1953, see also
Bear, 1972) involves a conformal mapping of Gz onto a mirror-image ofthe pentagon in the hodograph plane. Fig. 10d sketches such a half-planewith two semi-inﬁnite cuts for the case of an empty ditch.
Fig. 11 shows the function Qt (H) (curve 1) computed by Eq. (A9) for
hA = 3 and LA = 5. Similarly to Al-Shukaili et al. (2020b) and
Kacimov (1985), we introduce
the Morel-Seytoux dimensional seepage√
ﬂow rate 𝜇 = 𝑄∗𝑡 ∕(𝐾𝑠 𝐴∗ ), 𝐴∗ = 𝐻 ∗ 𝐿∗ .
Curve 2 in Fig. 11 plots the function 𝜇(H). The minimum of the function is 𝜇 ≈ 6.085 and H ≈ 1.177. This minimum can be interpreted as asolution to the following optimal shape design problem, OSDP (similar
to ones solved in Kacimov, 2005, 2006b):

OSDP 1. Determine an empty ditch, into which a minimal quantity𝑄∗𝑡 of water seeps, provided the ditch area A∗ and remote piezometricdata (ℎ∗𝐴 and 𝐿∗𝐴 ) are ﬁxed.
Curve 2 in Fig. 11 provides evidence that there is a unique andglobal solution to OSDP1. The minimum is “mild” (similar to ones in
Kacimov, 1985), i.e., triangular ditches having slopes close to the best
one will not deviate much from the best (minimal) seepage exﬁltration
rates. Problems similar to OSDP1 are common in civil engineering whena pit (e.g., for a building foundation) is excavated, and minimum exﬁltration is desired, or pit’s drainage measures are planned.
For comparisons, in Fig. 12, we present the results of HYDRUS simulations for the tetrad of parameters (H∗ , L∗ , 𝐿∗𝐴 ,ℎ∗𝐴 ) =(50 cm, 100 cm,
500 cm, 300 cm). Fig. 12a shows the isobars. Fig. 12b illustrates the isotachs. In the analytical solution, the Darcian velocity approaches inﬁnity
at points M and C. At the same time, it attains a minimal value, whichcorresponds to the tip of the cut in the hodograph domain (Fig. 10d).
HYDRUS shows the same. In Fig. 12c, we plot the distribution of thevelocities along the ditch side MC, where a minimum |V|=43.1 cm/day
is attained at sd =66.7 cm, where sd is the arc coordinate counted frompoint M in Fig. 10a. Curves like the one in Fig. 12c are useful for theassessment of erosional stability of slopes of ditches, since heaving andsuﬀusion are determined by exit values and directions of hydraulic gradient vectors ⃗𝑖. For example, PK-62,77 considers the condition |⃗𝑖| > 1 asa criterion of instability, according to which Fig. 12c illustrates that theditch side MC in Fig. 12 is absolutely unstable at any sd .

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 10. A vertical cross-section of seepage ﬂow towards a triangular ditch in a conﬁned aquifer a), a complex potential domain b), a reference plane c), and amirror-image in the hodograph plane d).

as in Fig. 11. In the regimes of partial ﬁlling, we can solve OSDPs similar
to OSDP1. Regional groundwater ﬂow can be easily taken into account,similarly to Ilyinsky and Kacimov (1992b). The analytical solution to
the steady problem can be readily extended to a transient regime of exﬁltration into a gradually-ﬁlling ditch. We assume that the value ℎ∗𝐴 in
Fig. 10a is high enough such that when we make a ditch, seepage into
it does not reduce much ℎ∗𝐴 , and the soil in Fig. 10a remains saturated(an aquifer remains conﬁned).
Let us assume, similarly as in HYDRUS, that during the transient
seepage phase, the porous skeleton and water are incompressible. If
evaporation from the ditch is ignored, then from the principle of massconservation, the following ODE follows:
d𝐴∗𝑠
d𝑡∗

Fig. 11. Dimensionless seepage ﬂow rates Qt (curve 1) and𝜇(curve 2) as functions of H for an empty ditch, hA = 3, and LA = 5.

The HYDRUS-computed Morel-Seytoux factor 𝜇 is equal to 6.97 forthe trench in Fig. 12, while the analytical solution gives 𝜇 =7.08. We
computed other triangular ditches in HYDRUS-2D, by using a discretevariation of the slope angle. We evaluated 𝑄∗𝑡 from the steady-state
limit of the cumulative ﬂux into the ditch from a “feeding semicircle”(Fig. 10a, dashed curve). Then we converted this dimensional 𝑄∗𝑡 into
the dimensionless 𝜇, which ﬁts well the analytical curve 2 in Fig. 11.
Collation of solutions for seepage into a triangular ditch andsink-generated ditch in Section 2 is done in the following manner.
For the case presented in Fig. 11, we use Eq. (2) and evaluate
Qd = 𝜋hA /LnLA = 3𝜋/Ln5 = 5.86. The corresponding isobar is shown
as curve 3 in Fig. 2. Using Eq. (6), we evaluated the area and the MorelSeytoux factor 𝜇 = 5.94. We surmise that this is a global optimum in theclass of arbitrary ditch shapes.
Next, we considered a general case of a partially-ﬁlled ditch(Fig. 10a). In Fig. 13, curves 1-3 show the graphs of the functions Qt (Hw )for three ditch slopes: 𝜋/6, 𝜋/4, and 𝜋/3, and the same hA = 3 and LA = 5[]= 𝑄∗𝑡 𝐻 ∗ (𝑡) ,(17)

where t∗ is dimensional time and 𝐴∗𝑠 (𝑡) is the area of the swelling triangle B1 C1 M in Fig. 10a (for comparisons see Al-Shukaili et al., 2020b,
where emptying of triangular trenches has been studied). Obviously, theinitial condition in Eq. (17) is 𝐴∗𝑠 (0) = 0, i.e., seepage-ﬁlling commencesfrom an empty ditch stage. We put 𝐴∗𝑠 (𝑡) = (𝐻 ∗ − 𝐻𝑤∗ (𝑡))2 cotan 𝜋𝜔 into
Eq. (17) and get the following Cauchy problem:() d𝐻𝑤∗ (𝑡)[]2 cotan 𝜋𝜔 𝐻𝑤∗ (𝑡) − 𝐻 ∗= 𝑄∗𝑡 𝐻𝑤∗ (𝑡) ,d𝑡∗𝐻𝑤∗ (0) = 0, 0 < 𝑡 < ∞.
(18)

We used the NDSolve routine of Mathematica to solve the nonlinear
1-st order ODE (18) numerically, with the RHS taken from Eq. (A9). For
this purpose, we interpolated 51 point values of 𝑄∗𝑡 [𝐻𝑤∗ ] (see curve 1)
in the steady-state problem. The results of computations are shown in
Fig. 13 in dimensionless quantities, to which we added a dimensionless time t = t∗ Ks /L. In Fig. 14, we plotted the functions Hw (t) obtained
from the solution of Eq. (18) for hA = 3 and LA = 5, and the same ditchslopes 𝜔 = 1/6, 1/4, and 1/3, (curves 1-3, correspondingly). As is evident from these draw-up curves, the ditches are ﬁlled out in ﬁnite timeintervals, tf =(0.115, 0.175, 0.24), respectively. At t > tf , seepage into
a ditch continues due to an artesian pressure in an aquifer, and exﬁltrating water has to be removed (e.g., pumped out) if we want to keep

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

Fig. 12. HYDRUS simulations for a triangular empty ditch having H∗ =50 cm, L∗ =100 cm, 𝐿∗𝐴 = 500 cm, and ℎ∗𝐴 =300 cm: isobars a), isotachs b), and Darcian velocities
along the ditch side MC c).

Fig. 13. The dimensionless seepage rate Qt into a triangular partially ﬁlled ditch
as a function of Hw for three slopes (𝜔=1/6, 1/4, and 1/3), hA =3, and LA =5.

Fig. 14. Dimensionless draw-up curves Hw (1, 2, and 3) as functions of dimensionless time t for diﬀerent values of 𝜔(1/6, 1/4, and 1.3,respectively).

points C and B in Fig. 10a at zero pressure heads (i.e., not ponded).
That is equivalent to the condition that the horizon BC remains static(see in Kacimov, 1991, a similar ﬂow regime in a horizontal plane).
The type of curves in Fig. 14 can be used for the determination of Ks(similarly to Al-Shukaili et al., 2020b and Kacimov, 2000b) that is also
useful for deep structure constructors when pumping tests in deep conﬁned aquifers are not feasible. If an aquifer is unconﬁned, then for asuddenly emptied ditch (or an auger hole, Kacimov 2000,b), tf = ∞A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732(Kacimov et al., 2019). Fig. 9e illustrates that routine measurements ofthe boundary condition on soil surfaces with depressions are not sufﬁcient for the determination of the directions and magnitudes of the
motion of water governed by the Richards equation.
7. Concluding remarks

Fig. 15. A draw-up curve computed by HYDRUS for 𝜔 =1/6, H∗ =50 cm, andL∗ =100 cm.

A transient analytical solution can be compared with a numericalone, which involves a new reservoir boundary condition of HYDRUS(Šimůnek et al., 2018). This new HYDRUS option allows for the consideration of variable storage of a furrow channel, or wellbore (e.g.,
Bristow et al., 2020, Sasidharan et al. 2018, 2019, 2020), the contoursof which are subject to the condition of a constant piezometric headon the submerged part and a seepage face on the empty part of thedraining entity. HYDRUS then considers ﬂow into or out of the reservoir
through its interface with the subsurface transport domain to depend on
the prevailing conditions in the transport domain, such as the positionof the groundwater table or piezometric heads, and external ﬂuxes, suchas pumping, or injection, or evaporation, or precipitation. The furrow
reservoir (Šimůnek et al., 2018) can be adapted for our triangular ditch.
The problem considered in Fig. 12 was rerun under transient conditions,starting with an empty ditch with dimensional sizes of H∗ =50 cm andL∗ =100 cm, corresponding to the case illustrated by curve 1 in Fig. 14.
The HYDRUS draw-up curve 𝐻𝑤∗ (𝑡∗ ) is shown in Fig. 15.
There is a perfect match between the results in Fig. 15 and curve
1 in Fig. 14 from the analytical solution in dimensionless variables.
We also compared the transient ﬂuxes into the half-ditch, i.e., theanalytical Qt /2 from Fig. 11 (curve 1) and one computed by HYDRUS.
For example, at t∗ =0.1day, the HYDRUS ﬂux is about 6075 cm2 /day,whereas the analytical ﬂux 𝑄∗𝑡 ∕2 is 6063 cm2 /day.
6. Comparison of terrestrial and martian evaporating spots
As an example, let us consider depression 1 (Qd =0.2) in Fig. 5. Assume R∗ =100 cm, H∗ =35 cm, and a loamy soil having Ks =25 cm/day.
Let the water evaporate from the disk 𝜋R∗ 2 with a constant intensity e.
In order to keep the crater empty, this value has to be e = 0.2 × 25/𝜋 ≈1.6 cm/day, which is - according to the Penman-Monteith equation andﬁeld measurements of evaporation - a too high value even in the Empty
Quarters of Oman.
If one wants to reconstruct the desertiﬁcation of Mars or intercept
ascending moisture ﬂuxes from the deep water table, such as, for example, in the process of Mars colonization, then the topography likethe one in Fig. 9 should be taken into account similarly as what wehave on Earth. Hummocky terrains or even earth dams (Kacimov andBrown, 2015) generate complex 2D and 3D Toth-type Darcian ﬂuxes.
Speciﬁcally, ﬂat landforms with a homogeneous vadose zone and a shallow water table are easy to examine for water storage and ﬂuxes: it
suﬃces to measure the depth of the water table and the pressure head(moisture content) of the topsoil to infer what is the value and directionof the ﬂux. If the vadose zone is heterogeneous, the ﬂux magnitudes canattain nontrivial extrema, even of a simplest type of a two-layered soil

The mainstream groundwater hydrology focuses on MAR (Managed
Aquifer Recharge), a method of replenishing aquifers, especially in arid
environments. We are putting forward the opposite concept of MAD(Managed Aquifer Discharge) of shallow aquifers, which have to be intelligently drained to mitigate the damage inﬂicted by groundwater inundation. The analysis in this paper will be useful in planned drainage
of water-logged built-up areas and further studies in dryland ecohydrology (Camporeale et al., 2019). Maillart’s desert ditches of Central Asia
(similarly as ones in humid Holland), can be viewed as a technique oflowering a shallow water table or piezometric surface. The shape ofthese ditches can be mathematically approximated by isobars made byhydrodynamic sinks, comparable with Zhukovskii’s method of a smart
composition of singularities in aerodynamics. The Darcian ﬂows governed by the Laplace equation in saturated media or ADE for steady
ﬂows in unsaturated media can be viewed as created by intelligently
located singularities. The resulting isobars model the outlets in seepage to empty curvilinear ditches and depressions. Seepage to triangular ditches, both empty and partially ﬁlled, is modeled by conformal
mappings and the BVP method. The transient stage of ﬁlling, i.e., thevariable storage in ditches, is examined by solving a Cauchy problem
for a nonlinear ODE with type curves showing how rapidly the ditch is
ﬁlled up. HYDRUS-2D simulations well agree with analytical modeling.
Analytical solutions and HYDRUS reconstruct jointly isobars, isohumes,
isotachs, and streamlines - the theoretical arsenal of soil physicists andcivil engineers working with Earth and Mars soil and water systems.
The proper design of ditches, pits, and other artiﬁcial excavations tobe operated in shallow water table environments, as well as an understanding of subsurface hydrology of natural wadis and depressions, canincorporate the ﬁndings of this paper.
Declaration of Competing Interest
The authors declare that they have no known competing ﬁnancial
interests or personal relationships that could have appeared to inﬂuence
the work reported in this paper.
Acknowledgments
This work was supported bySultan Qaboos University via thegrant “Rise of Water-table and Its Mitigation at SQU Campus”, grant
DR\RG\17 and it was carried out as part of the development programof the Scientiﬁc and Educational Mathematical Center of the Volga Federal District, agreement No. 075-02-2020-1478. Helpful comments by
two anonymous referees are highly appreciated.
Supplementary materials
Supplementary material associated with this article can be found, inthe online version, at doi:10.1016/j.advwatres.2020.103732.
Appendix
The Schwarz-Christoﬀel integral maps conformally the reference
half-plane on Gz :𝑧(𝜁 ) + i𝐻 = −𝑐1 (1 − i𝐻)= −(1 − i𝐻)

𝜁

∫0

−𝜔 2𝜔

(1 − 𝜏 2 )𝜏

B𝜁 2 (1∕2 + 𝜔, 1 − 𝜔)B(1∕2 + 𝜔, 1 − 𝜔)

d𝜏
,

(A1)

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek

Advances in Water Resources 145 (2020) 103732

where B and B𝜁 2 are the complete and incomplete beta-functions, respectively (Abramowitz and Stegun, 1969, formula 6.6.1). The positive constant 𝑐1 = 2∕B(1∕2 + 𝜔, 1 − 𝜔) in (A1) is determined from the condition
̄ = −𝑧(𝜁 )z(1) = −1. Note that z(𝜁 ) obeys the symmetry condition 𝑧(−𝜁)as it should be.
From Eq. (A1), along two seepage faces (BB1 and CC1 ), we have:()
B𝜉 2 (1∕2 + 𝜔, 1 − 𝜔)𝑦(𝜉) = 𝐻− 1 , 𝛽 ≤ |𝜉| ≤ 1,
(A2)
B(1∕2 + 𝜔, 1 − 𝜔)
where − 𝛽 is the aﬃx of point C1 . We determine 𝛽 from the relationy( − 𝛽) = −Hw , that gives:
B𝛽 2 (1∕2 + 𝜔, 1 − 𝜔)B(1∕2 + 𝜔, 1 − 𝜔)=

𝐻 − 𝐻𝑤.
𝐻(A3)

Next, we introduce a holomorphic function:
i𝑄𝑡(A4)
arcsin 𝜁 .
𝜋
The function (A4), in accordance with (16) and properties of thefunction arcsin, satisﬁes the conditions:
𝑊 (𝜁 ) = 𝑤(𝜁 ) −⎧Re 𝑊 (𝜉) = −𝑦(𝜉), 𝜉 ∈ (−1, −𝛽) ∪ (𝛽, 1),⎪
𝜉 ∈ (−𝛽, 𝛽),⎨Re 𝑊 (𝜉) = 𝐻𝑤 ,⎪Im 𝑊 (𝜉) = 0,𝜉 ∈ (−∞, −1) ∪ (1, ∞).
⎩Here y(𝜉) is determined in Eq. (A2). Then the function 𝑊0 (𝜁 ) =√𝑊 (𝜁 )∕ 1 − 𝜁 2 obeys the following conditions along the 𝜉-axis:√⎧Re𝑊 (𝜉) = −𝑦(𝜉)∕ 1 − 𝜉 2 , 𝜉 ∈ (−1, −𝛽) ∪ (𝛽, 1),
0
√
⎪
(A5)
𝜉 ∈ (−𝛽, 𝛽),⎨Re𝑊0 (𝜉) = 𝐻𝑤 ∕ 1 − 𝜉 2 ,⎪Re𝑊 (𝜉) = 0,
𝜉
∈
(−∞,
−1)
∪
(1
,
∞)
.
0
⎩
Eqs. (A5) comprise a Riemann (Schwartz) BVP (Henrici, 1993, PK62,77). The function W(𝜁 ) is ﬁnite at two transition points 𝜉 = ±1,
continuous at points ± 𝛽, and has a logarithmic singularity at inﬁnity. Hence, the function W0 (𝜁 ) has integrable singularities at points
𝜉 = ±1 and vanishes at inﬁnity. Therefore, the unique solution of the
problem (A5) gives the Schwartz operator for the upper half-plane(Henrici, 1993, PK-62,77). Thus, after simple algebra, we get
(
)
√
𝛽
1
2𝜁 1 − 𝜁 2𝑦(𝜏)d𝜏𝜏d𝑊 (𝜁 ) =
𝐻𝑤
−
,
√
√
∫0
∫𝛽
𝜋i
1 − 𝜏 2 (𝜏 2 − 𝜁 2 )1 − 𝜏 2 (𝜏 2 − 𝜁 2 )(A6)√
where, for the sake of deﬁniteness, the branch of the radical 1 − 𝜁 2
is ﬁxed positive for 𝜁 = 𝜉 ∈ ( − 1, 1) in the upper half-plane. Thisbranch and, therefore, the function (A6) satisfy the symmetry condition
𝑊 (−𝜁) = 𝑊 (𝜁 ). The integrals in (A6) are evaluated by the SokhotskiPlemely formula: the ﬁrst for 𝜁 → 𝜉 ∈ (0, 𝛽), and the second for 𝜁 → 𝜉∈ (𝛽, 1) (Henrici, 1993).
We illustrate the computations for the case of an empty ditch(Hw =H). We use the piezometric data at point Ah (Fig.10a), i.e.,𝜑A = −hA First, from Eq. (A1) along the ray CA we get:√𝑥(𝜉) = 1 + 2 1 + 𝐻 2−𝜉∫1− ∞ ≤ 𝜉 < −1.

−𝜔 2𝜔

(𝜏 2 − 1)𝜏d𝜏∕B(1∕2 + 𝜔, 1 − 𝜔),(A7)

From Eq. (A5) we determine the aﬃx − 𝛾 of point Ah in the reference
plane (Fig. 10c):
𝐿𝐴 = 𝑥(−𝛾).

(A8)We use the FindRoot routine of Mathematica to solve Eq. (A8) for 𝛾.
We put the √found root of Eq. (A8) into Eq. (A6), and at 𝜙𝐴 = 𝜙(−𝛾) =−𝑄𝑡 log(𝛾 + 𝛾 2 − 1)∕𝜋 + Re𝑊 (−𝛾) arrive at the equality:
(
)
√
𝛽
1
2 𝛾 𝛾 2 −1𝑦(𝜏)d𝜏
𝜏
d
𝑄𝑡 =
𝐻𝑤
−
√
√
√
∫0
∫𝛽
log(𝛾 + 𝛾 2 −1)
1−𝜏 2 (𝛾 2 −𝜏 2 )1−𝜏 2 (𝛾 2 −𝜏 2 )
𝜋𝜙𝐴
−
.
(A9)
√
log(𝛾 + 𝛾 2 −1)

References
Abit, S.M., Amoozegar, A., Vepraskas, M.J., Niewoehner, C.P., 2008. Solute transport inthe capillary fringe and shallow groundwater: Field evaluation. Vadose Zone J. 7 (3),890–898.
Abotalib, A.Z., Heggy, E., 2019. A deep groundwater origin for recurring slope lineae on
Mars. Nat. Geosci. 12 (4), 235.
Abramowitz, M., Stegun, I.A., 1969. Handbook of Mathematical Functions Dover, NewYork.
Abu-Rizaiza, O.S., 1999. Threats from groundwater table rise in urban areas in developing
countries. Water Int. 24 (1), 46–52.
Abu-Rizaiza, O.S., Sarikaya, H.Z., Khan, M.A., 1989. Urban groundwater rise control: casestudy. J. Irrigat. Drain. Eng. ASCE 115 (4), 588–607.
Acharya, S., Jawitz, J.W., Mylavarapu, R.S., 2012. Analytical expressions fordrainable and ﬁllable porosity of phreatic aquifers under vertical ﬂuxes
from evapotranspiration and recharge. Water Resour. Res. 48, W125126.
https://doi.org/10.1029/2012WR012043.
Afruzi, A., Nazemi, A.H., Sadraddini, A.A., 2014. Steady‐state subsurface drainage of
ponded ﬁelds by rectangular ditch drains. Irrigat. Drain. 63 (5), 668–681.
Al-Rawas, A.A., Qamaruddin, M., 1998. Construction problems of engineering structures
founded on expansive soils and rocks in northern Oman. Build. Environ. 33 (2–3),
159–171.
Al-Sefry, S.A., Şen, Z., 2006. Groundwater rise problem and risk evaluation in major cities
of arid lands–Jedddah case in Kingdom of Saudi Arabia. Water Resour. Manage. 20(1), 91–108.
Al-Senafy, M., 2011. Management of water table rise at Burgan oil ﬁeld, Kuwait. Emir.
Emirates J. Eng. Res. 16 (2), 27–38.
Al-Senafy, M., Hadi, K., Fadlelmawla, A., Al-Fahad, K., Al-Khalid, A., Bhandary, H., 2015.
Causes of groundwater rise at Al-Qurain residential area, Kuwait. Procedia Environ.
Sci. 25, 4–10.
Alsharhan, A.S., Rizk, Z.E., 2020. Water Resources and Integrated Management of theUnited Arab Emirates (Vol. 3). Springer Nature, Cham.
Alsharhan, A.S., Rizk, Z.A., Nairn, A.E.M., Bakhit, D.W., Alhajari, S.A. (Eds.), 2001. Hydrogeology of an Arid Region: the Arabian Gulf and Adjoining Areas. Elsevier.
Al-Shukaili, A., Al- Busaidi, H., Kacimov, A.R., 2020. Experiments, analytical and HYDRUS2D modeling of steady jet of quasi-normal surface ﬂow in rectangular channel
coupled with vertical seepage: Vedernikov-Riesenkampf’s legacy revisited. Adv. Water Resour. 136. https://doi.org/10.1016/j.advwatres.2019.103503.
Al-Shukaili, A., Al- Mayahi, A., Al-Maktoumi, A., Kacimov, A.R., 2020. Unlined trench
as a falling head permeameter: Analytic and HYDRUS2D modeling versus sandbox
experiment. J. Hydrol. 583, 124568.
Al-Yaqout, A, Townsend, F., 2004. Applicability of caliche (gatch) as liner/cover in arid
climate landﬁlls: Laboratory and ﬁeld pad testing of permeability. Pract. Period. Hazard. Toxic Radioact. Waste Manage. 8 (4), 238–246.
Anderson, E.I., 2013. Stable pumping rates for horizontal wells in bank ﬁltration systems.
Adv. Water Resour. 54, 57–66.
Aravin, V. I., and Numerov, S. N., 1953. Theory of Fluid Flow in Undeformable PorousMedia, Gostekhizdat, Moscow. Engl. Translation: 1965 by the Israel Program Scientiﬁc
Translations, Jerusalem, Israel.
Attard, G., Winiarski, T., Rossier, Y., Eisenlohr, L., 2016. Impact of underground structureson the ﬂow of urban groundwater. Hydrogeol. J. 24 (1), 5–19.
Barron, O.V., Donn, M.J., Barr, A.D., 2013. Urbanisation and shallow groundwater: Predicting changes in catchment hydrological responses. Water Resour. Manage. 27,
95–115.
Barua, G., Alam, W., 2013. An analytical solution for predicting transient seepage into
ditch drains from a ponded ﬁeld. Adv. Water Resour. 52, 78–92.
Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York Elsevier, NewYork.
Belotserkovsky, S.M., Lifanov, I.K., 1992. Method of Discrete Vortices. CRC Press see alsoБелоцерковский, С.М., Лифанов, И.К., 1985. Численные методы в сингулярных
интегральных уравнениях и их применение в аэродинамике, теории упругости,электродинамике. Москва, Наука.
Bhardwaj, A., Sam, L., Martín‐Torres, F.J., Zorzano, M.P., 2019. Are slope streaks indicative of global‐scale aqueous processes on contemporary Mars? Rev. Geophys. 57 (1),48–77.
Boatwright, B.D., Head, J.W., 2019. Simulating early Mars hydrology with the MARSSIM
landform evolution model: New insights from an integrated system of precipitation,
inﬁltration, and groundwater ﬂow. Planet. Space Sci. 171, 17–33.
Bob, M., Rahman, N., Elamin, A., Taher, S., 2016. Rising groundwater levels problem inurban areas: a case study from the Central Area of Madinah City, Saudi Arabia. Arab.
J. Sci. Eng. 41 (4), 1461–1472.
Bristow, K.L., Šimůnek, J., Helalia, S.A., Siyal, A.A., 2020. Numerical simulations of
the eﬀects furrow surface conditions and fertilizer locations have on plant nitrogen and water use in furrow irrigated systems. Agricult. Water Manage. 232, 11.
https://doi.org/10.1016/j.agwat.2020.106044, 106044.
Camporeale, C., Perona, P., Ridolﬁ, L., 2019. Hydrological and geomorphological signiﬁcance of riparian vegetation in drylands. In: Odorico, P.D, Porporato, A. (Eds.), Dryland Ecohydrology. Springer, Cham, pp. 239–275.
Chahar, B.R., Vadodaria, G.P., 2008. Drainage of ponded surface by an array of ditches.
J. Irrigat. Drain. Eng. ASCE 134 (6), 815–823.
Chahar, B.R., Vadodaria, G.P., 2012. Steady subsurface drainage of ponded surface by anarray of parallel ditches. J. Hydrol. Eng. ASCE 17 (8), 895–908.
Chaudhary, M.T.A., 2012. Implications of rising groundwater level on structural integrityof underground structures – investigations and retroﬁt of a large building complex.
Struct. Surv. 30 (2), 111–129.

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek
Coda, S., Confuorto, P., De Vita, P., Di Martire, D., Allocca, V., 2019. Uplift evidences
related to the recession of groundwater abstraction in a pyroclastic-alluvial aquifer ofSouthern Italy. Geosciences 9 (5), 1–15 215.
Dagan, G., 1964. Linearized solution of unsteady deep ﬂow toward an array of horizontal
drains. J. Geophys. Res. 69, 3361–3369.
Duniway, M.C., Herrick, J.E., Monger, H.C., 2010. Spatial and temporal variability ofplant-available water in calcium carbonate-cemented soils and consequences for arid
ecosystem resilience. Oecologia 163 (1), 215–226.
Edwards, C.S., Piqueux, S., 2016. The water content of recurring slope lineae on Mars.
Geophys. Res. Lett. 43 (17), 8912–8919.
Freeze, R.A., Cherry, J.A., 1979. Groundwater. Prentice-Hall Inc, Englewood Cliﬀs.
Fujii, N., Kacimov, A.R., 1998. Analytically computed rates of seepage ﬂow into drains
and cavities. Int. J. Numer. Anal. Method Geomech. 22, 277–301.
Gardner, W.R., 1958. Some steady-state solutions of the unsaturated moisture ﬂow equation with application to evaporation from a water table. Soil Sci. 85, 228–232.
Ghezzehei, T.A., Kneafsey, T.J., Su, G.W., 2007. Correspondence of the Gardner and van
Genuchten–Mualem relative permeability function parameters. Water Resour. Res.
40, W10417. https://doi.org/10.1029/2006WR005339, 2007.
Gjerde, K.M., Tyvand, P.A., 1992. Transient free‐surface groundwater ﬂow due to an arrayof circular drainage ditches. Water Resour. Res. 28 (11), 2963–2972.
Goldspiel, J.M., Squyres, S.W., 2011. Groundwater discharge and gully formation on Martian slopes. Icarus 211 (1), 238–258.
Goudie, A.S., 2013. Arid and Semi-arid Geomorphology. Cambridge University Press.
Grimm, R.E., Harrison, K.P., Stillman, D.E., 2014. Water budgets of Martian recurring
slope lineae. Icarus 233, 316–327.
Grotzinger, J.P., Sumner, D.Y., Kah, L.C., Stack, K., Gupta, S., Edgar, L., Milliken, R., 2014.
A habitable ﬂuvio-lacustrine environment at Yellowknife Bay, Gale Crater, Mars. Science 343 (6169).
Gureghian, A.B., Youngs, E.G., 1975. The calculation of steady-state water-table heights
in drained soils by means of the ﬁnite-element method. J. Hydrol. 27 (1–2), 15–32.
Heilweil, V.M., Watt, D.E., 2011. Trench inﬁltration for managed aquifer recharge to permeable bedrock. Hydrol. Process. 25 (1), 141–151.
Henrici, P., 1993. Applied and Computational Complex Analysis. Volume 3: Discrete
Fourier Analysis, Cauchy Integrals, Construction of Conformal Maps, Univalent Functions. Wiley, New York.
Hobbs, S.W., Paull, D.J., Clarke, J.D.A., 2014. A hydrological analysis of terrestrial andMartian gullies: Implications for liquid water on Mars. Geomorphology 226, 261–277.
Howard, K.W., Israﬁlov, R.G. (Eds.), 2012. Current Problems of Hydrogeology in Urban
Areas, Urban Agglomerates and Industrial Centres (NATO Science Series, Vol. 8).
Springer, Dordrecht.
Ilyinsky, N.B., Kacimov, A.R., 1992. Problems of seepage to empty ditch and drain. Water
Resour. Res. 28 (3), 871–877.
Ilyinsky, N.B., Kacimov, A.R., 1992. Analytical estimation of ground-water ﬂow around
cutoﬀ walls and into interceptor trenches. Ground Water 30, 901–907.
Jha, A.K., Bloch, R., Lamond, J., 2012. Cities and Flooding: a Guide to Integrated Urban
Flood Risk Management for the 21st Century. The World Bank, Washington, DC.
Kacimov, A., Obnosov, Yu.V., 2019. Modeling of 2-D Seepage from aquifer towards
stream via clogged bed: the Toth-Treﬀtz legacy conjugated. Adv. Water Resour. 129,
231–251.
Kacimov, A., Obnosov, Yu.V., Simunek, J., 2019. Minimizing evaporation by optimal layering of topsoil: revisiting Ovsinsky’s smart mulching-tillage technology via Gardner-Warrick’s unsaturated analytical model and HYDRUS. Water Resour. Res. 55 (5),3606–3618.
Kacimov, A.R., 1985. Optimization of the shape of a triangular unlined canal. Power Tech.
Eng. 19 (1), 41–43.
Kacimov, A.R., 2000. Circular isobaric cavity in descending unsaturated ﬂow. J. Irrigat.
Drain. Eng. ASCE 126 (3), 172–178.
Kacimov, A.R., Brown, G., 2015. A transient phreatic surface mound, evidenced by a stripof vegetation on an earth dam. Hydrol. Sci. J. 60 (2), 361–378.
Kacimov, A.R., Obnosov, Yu.V., 2008. Leaky-layer seepage: the Verigin function revisited.
J. Eng. Math. 62, 345–354.
Kacimov, A.R., 1991. Steady, two-dimensional ﬂow of ground water to a trench. J. Hydrol.
127, 71–83.
Kacimov, A.R., 1992. Minimum dimension of total saturation bubble around an isolated
source. Fluid Dyn. 27, 886–889.
Kacimov, A.R., 2000. Analytic solution for transient ﬂow into a hemispherical auger hole.
J. Hydrol. 228, 1–9 v.
Kacimov, A.R., 2006. Analytic element solutions for seepage towards topographic depressions. J. Hydrol. 318, 262–275.
Kacimov, A.R., 2006. Analytical solution and shape optimisation for groundwater ﬂow
through a leaky porous trough subjacent to an aquifer. Proc. R. Soc. Lond. A 462,
1409–1423.
Kacimov, A.R., 2007. Dipole-generated unsaturated ﬂow in Gardner soils. Vadose Zone J.
6, 168–174.
Kacimov, A.R., 2009. Minimal-seepage depth of isobaric cavity under ponded conditions.
J. Irrigat. Drain. Eng. ASCE 135 (1), 108–110.
Kacimov, A.R., Al-Jabri, S., Sherif, M.M., Al-Shidi, S., 2009. Slumping of groundwatermounds: revisiting the Polubarinova-Kochina theory and modeling by analytic element method. Hydrol. Sci. J. 54 (1), 174–188.
Kacimov, A.R., Obnosov, Yu.V., 2002. Analytical determination of seeping soil slopes of aconstant exit gradient. Zeitschrift fur angewandte Mathematik und Mechanik 82 (6),363–376.
Kacimov, A.R., Youngs, E.G., 2005. Steady-state water-table depressions caused by evaporation in lands overlying a water-bearing substratum. J. Hydrol. Eng. (ASCE) 10 (4),295–301.

Advances in Water Resources 145 (2020) 103732
Kacimov, A.R, 2005. Seepage to a drainage ditch and optimization of its shape. J. Irrigat.
Drain. Eng. (ASCE) 132 (6), 619–622.
Kazemi, G.L., 2011. Impacts of urbanization on the groundwater resources in Shahrood,
Northeastern Iran: Comparison with other Iranian and Asian cities. Phys. Chem. Earth36, 150–159.
Kereszturi, A., Möhlmann, D., Berczi, S., Ganti, T., Horvath, A., Kuti, A., Sik, A., Szathmary, E., 2010. Indications of brine related local seepage phenomena on the northernhemisphere of Mars. Icarus 207 (1), 149–164.
Kirkham, D., 1947. Studies of hillslope seepage in the Iowan drift area. Proc. Soil Sci. Soc.
Am. 12, 73–80.
Kirkham, D., Powers, W.K., 1972. Advanced Soil Physics. Wiley, New York.
Kochel, R.C., Piper, J.F., 1986. Morphology of large valleys on Hawaii—Evidence forgroundwater sapping and comparisons with Martian valleys. J. Geophys. Res. 91(B13), E175–E192.
Kocurek, G., Westerman, R., Hern, C., Tatum, D., Rajapara, H.M., Singhvi, A.K., 2020. Aeolian dune accommodation space for Holocene Wadi Channel Avulsion Strata, Wahiba
Dune Field, Oman. Sediment. Geol. 399, 105612.
Kreibich, H., Thieken, A.H., 2008. Assessment of damage caused by high groundwater
inundation. Water Resour. Res. 44, W09409.
Lerner, D.N.ed, 2003. Urban Groundwater Pollution: IAH International Contributions toHydrogeology 24. Balkema, Lisse.
Lerner, D.N., Harris, B, 2008. The relationship between land use and groundwater resources and quality. Land Use Policy 26S (2009), S265–S273.
Luo, W., Grudzinski, B., Pederson, D., 2011. Estimating hydraulic conductivity for theMartian subsurface based on drainage patterns—A case study in the Mare Tyrrhenum
Quadrangle. Geomorphology 125 (3), 414–420.
Malin, M.C., Carr, M.H., 1999. Groundwater formation of Martian valleys. Nature 397(6720), 589–591.
Malin, M.C., Edgett, K.S., 2000. Evidence for recent groundwater seepage and surface
runoﬀ on Mars. Science 288 (5475), 2330–2335.
Marra, W.A., Braat, L., Baar, A.W., Kleinhans, M.G., 2014. Valley formation by groundwater seepage, pressurized groundwater outbursts and crater-lake overﬂow in ﬂume
experiments with implications for Mars. Icarus 232, 97–117.
Marra, W.A., McLelland, S.J., Parsons, D.R., Murphy, B.J., Hauber, E., Kleinhans, M.G.,
2015. Groundwater seepage landscapes from distant and local sources in experiments
and on Mars. Earth Surf. Dyn. 3 (3), 389–408.
Medovar, Y.A., Yushmanov, I.O., Bronskaya, E.E., 2018. Changes in the ecological condition of a territory at the construction of buildings with a deep foundation undervarious geological conditions in Moscow. Water Resour. 45 (2), 65–72.
Michalski, J.R., Cuadros, J., Niles, P.B., Parnell, J., Rogers, A.D., Wright, S.P., 2013.
Groundwater activity on Mars and implications for a deep biosphere. Nat. Geosci.
6 (2), 133–138.
Mujica, C.R., Bea, S.A., 2020. Estimations of rooting depths and sources of plant-available
water (PAW) in ﬂatland petrocalcic soils under diﬀerent land uses. Geoderma 361,
114019.
Mukherjee, S., Singh, D., Singh, P., Roy, N., 2020. Morphological and morphometric analysis of a topographic depression near Huygens basin, Mars: Identiﬁcation of a putativeendorheic playa. Geomorphology 351, 106912.
Naik, P.K., Tambe, J.A., Dehury, B.N., Tiwari, A.N., 2008. Impact of urbanization on the
groundwater regime in a fast growing city in central India. Environ. Monitor. Assess.
146 (1–3), 339–373.
Obnosov, Yu.V., Kacimov, A.R., 2018. Steady Darcian ﬂow in subsurface irrigation of
topsoil impeded by substratum: Kornev-Riesenkampf-Philip legacies revisited. Irrigat.
Drain. 67 (3), 374–391.
Philip, J.R., 1968. Steady inﬁltration from buried point sources and spherical cavities.
Water Resour. Res. 4, 1039–1047.
Philip, J.R., 1971. General theorem on steady inﬁltration from surface sources, with application to point and line sources. Soil Sci. Soc. Amer. Proc. 35, 867–871.
Philip, J.R., 1989. Multidimensional steady inﬁltration to a water table. Water Resour.
Res. 25, 109–116.
Polubarinova-Kochina, P.Ya., 1962. Theory of Ground Water Movement. Princeton University Press, Princeton The second edition of the book (in Russian) was published in
1977, Nauka, Moscow.
Porse, E., Glickfeld, M., Mertan, K., Pincetl, S., 2016. Pumping for the masses: evolutionof groundwater management in metropolitan Los Angeles. GeoJ. 81 (5), 793–809.
Preene, M., Fisher, S., 2015. Impacts from groundwater control in urban areas. In: Proceedings of the XVI ECSMGE, pp. 2847–2852.
Pullan, A.J., 1990. The quasilinear approximation for unsaturated porous media ﬂow.
Water Resour. Res. 26 (6), 1219–1234.
Quan, R.S., Liu, M., Lu, M., Zhang, L.J., Wang, J.J., Xu, S.Y., 2010. Waterlogging risk assessment based on land use/cover change: a case study in Pudong New Area, Shanghai. Environ. Earth Sci. 61 (6), 1113–1121.
Raats, P.A., Smiles, D., Warrick, A.W.eds., 2002. Environmental Mechanics: Water, Massand Energy Transfer in the Biosphere. The Philip Volume. American GeophysicalUnion, Washington, DC.
Raats, P.A.C., 1970. Steady inﬁltration from line sources and furrows. Soil Sci. Soc. Amer.
Proc. 34, 709–714.
Raats, P.A.C., 1971. Steady inﬁltration from point sources, cavities and basins. Soil Sci.
Soc. Amer. Proc. 35, 689–694.
Raats, P.A.C., 1972. Steady inﬁltration from sources at arbitrary depth. Soil Sci. Soc. Amer.
Proc. 36, 399–401.
Raats, P.A.C., 1977. Laterally conﬁned, steady ﬂows from sources and to sinks in unsaturated soils. Soil Sci. Soc. Amer. J. 41, 294–304.
Radcliﬀe, D.E., Šimůnek, J., 2010. Soil Physics with HYDRUS: Modeling and Applications.
CRC Press, Taylor & Francis.

A.R. Kacimov, Yu.V. Obnosov and J. Šimůnek
Salese, F., Pondrelli, M., Neeseman, A., Schmidt, G., Ori, G.G., 2019. Geological evidenceof planet‐wide groundwater system on Mars. J. Geophys. Res. 124 (2), 374–395.
Sarmah, R., Tiwari, S., 2018. A two-dimensional transient analytical solution for a ponded
ditch drainage system under the inﬂuence of source/sink. J. Hydrol. 558, 196–204.
Sasidharan, S., Bradford, S.A., Šimůnek, J., DeJong, B., Kraemer, S.R., 2018. Evaluating
drywells for stormwater management and enhanced aquifer recharge. Adv. Water
Resour. 116, 167–177. https://doi.org/10.1016/j.advwatres.2018.04.003.
Sasidharan, S., Bradford, S.A., Šimůnek, J., Kraemer, S.R., 2019. Drywell inﬁltration and hydraulic properties in heterogeneous soil proﬁles. J. Hydrol. 570.
https://doi.org/10.1016/j.jhydrol.2018.12.073, 598-561.
Sasidharan, S., Bradford, S.A., Šimůnek, J., Kraemer, S.R., 2020. Groundwater
recharge from drywells under constant head conditions. J. Hydrol. 583, 14.
https://doi.org/10.1016/j.jhydrol.2020.124569, 124569.
Schirmer, M., Leschik, S., Musolﬀ, A., 2013. Current research in urban hydrogeology–A
review. Adv. Water Resour. 51, 280–291.
Silliman, S.E., Berkowitz, B., Šimůnek, J, van Genuchten, M.Th., 2002. Fluid ﬂow andsolute migration within the capillary fringe. Ground Water 40 (1), 76–84.
Šimůnek, J., Šejna, M., van Genuchten, M.Th., 2018. New features of version 3 of the
HYDRUS (2D/3D) computer software package. J. Hydrol. Hydromech. 66 (2), 133–
142. https://doi.org/10.1515/johh-2017-0050.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016. Recent developments and applications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Skaggs, R.W., Youssef, M.A., Chescheir, G.M., 2012. DRAINMOD: Model use, calibration,and validation. Trans. ASABE 55 (4), 1509–1522.

Advances in Water Resources 145 (2020) 103732
Skaggs, R.W., Van Schilfgaarde, J., Bartels, J.M., Hatﬁeld, J.L., Volenec, J.J.,
Bigham, J.M.eds., 1999. Agricultural Drainage. American Society of Agronomy, Madison, WI.
Strack, O.D.L., 1989. Groundwater Mechanics. Prentice Hall, Englewood Cliﬀs.
Strack, O.D.L., 2020. Applications of Vector Analysis and Complex Variables in Engineering. Springer, Cham.
Tóth, J., 2009. Gravitational Systems of Groundwater Flow: Theory, Evaluation, Utilization. Cambridge University Press, New York.
Vázquez-Suñé, E., Sánchez-Vila, X., Carrera, J., 2005. Introductory review of speciﬁc factors inﬂuencing urban groundwater, an emerging branch of hydrogeology, with reference to Barcelona, Spain. Hydrogeol. J. 13 (3), 522–533.
Vedernikov, V.V., 1939. Theory of Seepage and Its Application in Irrigation and Drainage.
Gosstrojizdat, Moscow in Russian.
Vogwill, R. (Ed.), 2016. Solving the Groundwater Challenges of the 21st Century (Vol.
22). CRC Press, Boca Raton.
Wolfram, S., 1991. Mathematica. A System for Doing Mathematics by Computer. Addison-Wesley, Redwood City.
Youngs, E.G., 1975. The eﬀect of the depth of an impermeable barrier on water-tableheights in drained homogeneous soils. J. Hydrol. 24 (3–4), 283–290.
Youngs, E.G., 1990. An examination of computed steady-state water-table heights in unconﬁned aquifers: Dupuit-Forchheimer estimates and exact analytical results. J. Hydrol. 119 (1–4), 201–214.
Youngs, E.G., 1994. Seepage to ditches from a ponded surface. J. Hydrol. 161 (1–4),
145–154.
Zhukovskii, N.E., 1948. Collection of Works. Gostekhizdat, Moscow in Russian.