Icarus 372 (2022) 114719

Contents lists available at ScienceDirect

Icarusjournal homepage: www.elsevier.com/locate/icarus

Seepage-evaporation controlled depletion of initially water-filled reservoirs
on Earth and Mars: Analytic versus HYDRUS modeling
A.R. Kacimov a, *, Yu.V. Obnosov b, J. Šimůnek caDepartment of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Oman
Institute of Mathematics and Mechanics, Kazan Federal University, Kazan, Russiac
Department of Environmental Sciences, University of California Riverside, CA, USAb

A R T I C L E I N F O

A B S T R A C T

Keywords:
Seepage from axisymmetric craters
Gardnerian unsaturated seepage
Hunt’s saturated seepage from paraboloidalcap
HYDRUS modeling

Analytical solutions are obtained for water extinction from an axisymmetric crater, filled at t < 0 and depleted byevaporation and transient infiltration into a Gardner or capillarity-free homogeneous soil during the time in­terval 0 ≤ t ≤ Te. The extinction time Te is found for crater beds, the shapes of which are shallow cones, spherical,
spheroidal, and paraboloidal caps. An instantaneous seepage flow rate, Q(t), is approximated by truncated twoterm formulae of Wooding for a zero-depth disk in Gardner's soil or Hunt for paraboloidal craters (soils with nocapillarity). The instantaneous evaporation losses are the product of a constant A-pan evaporation rate and theshrinking area of a flat horizontal disk of the free water, which dwindles in the crater. In HYDRUS simulations ofa van Genuchten soil, the Reservoir Boundary Condition is used for a falling water level in the ponded de­pressions. Cones and paraboloids are selected as craters, initially fully or partially filled with free water at t = 0,
and infiltrating until extinction. The results are presented as drawdown curves and – for shallow craters - attest agood match between analytical approximations and HYDRUS numerical simulations. Experiments with theextinction of water from small axisymmetric ponds in dune sand are also carried out. They allow blitz-evaluationof hydraulic parameters of the subjacent sand. Hydrological implications for commingling surface-subsurface(pore) water entities in terrestrial and Martian environments are discussed.

List of abbreviations and notations
SymbolGA
PK-62,
77
VG
В
θ
θr
θs
Ks
n, α

Description
Green-Ampt
Polubarinova-Kochina, P.Ya., 1962, 1977. Theory of
Ground Water Movement.
Van Genuchten
Ponded water depth in the crater
Volumetric soil water content
Residual water content
Saturated water content (porosity)
Saturated hydraulic conductivity
VG soil parameters

αP

Sorptive number
Pressure head
Total (piezometric) head
Time
Radial coordinate
Vertical coordinate
Volume of free water in the pond

p
h
t
r
z
V

Unit
–
–
–
cm
[− ]
[− ]
[− ]
cm/h
[− ],
cm− 1
cm− 1
cm
cm
h
cm
cm
cm3

“That last day does not bring extinction to us, but change of place.”
Marcus Tullius Cicero
1. Introduction
Water in all forms (surface water, groundwater, liquid soil moisture,ice, and vapor in both atmosphere, on the land surface, and inside thepore space) and all types of its transport in the hydrological cycle(runoff, seepage, infiltration, deep percolation, groundwater motion,
evaporation, precipitation) is of paramount importance for all humanactivities (agriculture, civil engineering, etc.). Recent irrefutable NASAdiscoveries of water on the Moon and Mars and other planets and theirsatellites in the Solar system determine the long-term strategies ofconquering Space. The hydrology, hydrogeology, hydrogeochemistry,* Corresponding author.
E-mail addresses: anvar@squ.edu.om (A.R. Kacimov), yobnosov@kpfu.ru (Yu.V. Obnosov), jsimunek@ucr.edu (J. Šimůnek).
https://doi.org/10.1016/j.icarus.2021.114719
Received 10 December 2020; Received in revised form 14 September 2021; Accepted 21 September 2021
Available online 28 September 2021
0019-1035/© 2021 Elsevier Inc. All rights reserved.

A.R. Kacimov et al.

Icarus 372 (2022) 114719

Let us introduce the pore pressure, P, P = ρ g p, where p [cm]1 is thepressure head, ρ [g/cm3] is liquid's (water, brine, oil, etc.) density, and g[cm/s2] is gravity acceleration. Obviously, ρ is planet-invariant, and g =981 and 371 cm/s2 on Earth and Mars, respectively. We use theterrestrial acceleration in all our computations below (rescaling to Marsor other gravitating entities is straightforward, modulo g).
In this paper, we consider axisymmetric groundwater and vadosezone motion, which obeys the Richards-Richardson equation:

sedimentology, soil physics-mechanics of hyperarid deserts on the Earthserve as proving grounds to infer the physical properties of the liquid,
solid, and gaseous phases of these planets. This will facilitate theforthcoming space missions and elucidate the past geological epochs of,say, Mars in backward modeling-reconstructions (see, e.g., Davila et al.,2021, Groemer et al., 2020).
In this paper, we extend our analytical and numerical studies oftransient 2-D seepage from canals and “dry wells” (Al-Shukaili et al.,2019, 2020; Kacimov et al., 2020, 2021; Sasidharan et al., 2018) toaxisymmetric flows from non-cylindrical reservoirs with water levels
falling due to seepage into terrestrial/Martian soils and evaporation into
the atmosphere. Interaction (exchange) between surface and subsurface
waters has been thoroughly investigated in hydrology, limnology, soil
physics, geomorphology, and other Earth sciences (see, e.g., Sopho­
cleous, 2002, Woessner, 2020). Paleodynamics of depletion of the freewater storage in Martian lakes/reservoirs came to space scientists'
attention upon NASA indagation and the dazzling, albeit so far byremote sleuthing, discovery of deep groundwater (e.g., by a GPR tech­nique). Mathematical modelers also started to portray the topology of
Darcian seepage in Martian aquifers. In the current purview of the subdiscipline of Martian hydrology, a hypothesized dynamics of surface
water in lakes, channels, and other fluvial or lacustrine entities iscommingled with saturated Darcian flows (Baum and Wordsworth,
2020, Bhardwaj et al., 2019, Boatwright and Head, 2019, Fukushi et al.,2019, Goldspiel and Squyres, 2011, Grimm et al., 2014, 2017, Hobbs
et al., 2014, Kochel and Piper, 1986, Luo et al., 2011, Kereszturi et al.,2010, Malin and Carr, 1999, Malin and Edgett, 2000, Marra et al., 2014,2015, Michalski et al., 2013, Mukherjee et al., 2020, Salese et al., 2019).
The commonality of landforms on Mars and Earth, for instance, gullies
or catchments' drainage networks, and the terrestrial scrutiny of the
impact of Darcian flows on their evolution (see, e.g., Blank et al., 2009,
Conway et al., 2019 versus Petroff et al., 2013) stimulate further studies
of subsurface mechanics on the two planets.
Soil physics of infiltration, deep percolation, groundwater motion in
unconfined aquifers and the capillary fringe above the water tables,
evaporation, redistribution, root water uptake, runoff, phase changes ofpore water (including freezing-thawing and vapourizationcondensation), and heat and solute transport in unsaturated porous
media have been thoroughly studied, including validation of models via
field and laboratory experiments, in terrestrial systems (Radcliffe andŠimůnek, 2010; Šimůnek et al., 2018). Without claiming a comprehen­sive review of the whole scope of reconnoitered extra-terrestrial Darican
flows, it seems that the unsaturated and tension-saturated water (brine)
flows on Mars (including the vadose zone and capillary fringe) havebeen ignored in mathematical models, i.e., the porous Martian media
were posited capillarity-free. Only recently, qualitative statements
about the role of capillarity have emerged: “It is possible that a change insoil permeability might affect the local capillary rise and thus the surface
visibility” (Lark et al., 2020). Neglecting capillarity is reasonable for large
Martian craters and mega-channels. However, for the forthcoming
human missions to Mars and the realization of Tsiolkovsky-Musk's
projects (populating Mars and other planets/satellites), a prodigious
interest to apposite models of infiltration and saturated-unsaturated
seepage from small-size man-made excavations, temporarily filledwith water, is emerging. For example, human colonies on Mars may
dispose of sewage in Martian lagoons. The astro-physio-hydro-biologicalquestion will then arise: how fast do the impounded-trickled liquids
percolate-evaporate in the triad of Martian sediments-rocksatmosphere? Upon comparisons of the above referenced Darcian
models, used by Martian hydrologists, with the legacy of terrestrial
irrigation-drainage sciences and engineering (Kacimov et al., 2020), we
deem that several aspects of transient seepage still deserve to beunearthed.

∂θ= ∇(K(p)∇h )∂t(1)where ∇ is the nabla operator in cylindrical coordinates, θ(t, r, z) is a
volumetric water content [− ], K(p) [cm/h] is the hydraulic conductivity
function, h(t, r, z) =p + z is the total (piezometric) head [cm], z is avertical coordinate oriented upward. It is well-known (PolubarinovaKochina, 1962, 1977, hereafter abbreviated as PK-62,77, Strack, 2017)that the intrinsic permeability, κ = K μ/(ρ g) [cm2], where μ [cP] isfluid's dynamic viscosity, as well as the capillary pressure function p(θ )(modulo sign, equivalent to the water retention function), and relative
water permeability, κ(θ ), depend on the properties of the porous skel­
eton only (provided water and air, as two fluids occupying the porespace, are fixed). The two functions are invariants in the whole Solar
System (modulo g-rescaling).
We present analytical and numerical solutions to boundary value
problems for eq. (1), which is coupled with the depletion of the storageof hydrostatically modeled water in a surface pond (crater). This inter­action of surface-subsurface water brings about a new (Reservoir)
boundary condition for axisymmetric flows, implemented recently in
HYDRUS (Šimůnek et al., 2018), which accounts for a falling headproblem for reservoirs of different shapes. This condition was analyti­
cally involved in 1-D seepage flows (Philip, 1992) and 2-D flows in a
vertical plane (Al-Shukaili et al., 2020).
The results of this paper will outreach to soil hydrologists and water
resources managers working in arid zones, where topographical de­pressions on Earth's surface are temporarily filled by episodic rains or
irrigation events (e.g., sprinkling over crop fields). For example, the2018 Mekunu cyclone caused torrential rains in the dune fields of thehyperarid Empty Quarter in Oman and Saudi Arabia, with large inter­
dunal valleys converted into ephemeral lakes. For local Bedouins con­
fronting exiguous water sources, a vital question was: How long can thepond be utilized as an ephemeral oasis to provide drinking water forcamels? The longer, the better. At the instance when the extant surface
water vanishes, groundwater pumping has to start providing watersupply to the cattle. However, the infiltrated water may still percolate
through a vadose zone before reaching a deep aquifer from which
pumping is feasible. Environmentally opposite abject scenarios emerge
in illegal-clandestine dumping of toxic liquid waste in depressions ofdesert sand or wadis. To detect and catch the dumping criminals, theslower infiltration of the waste into the desert sand, the better. In theensuing litigations, sampling of the extant liquid from a pond andcomparing water quality in it with one in remote contaminatedgroundwater wells is the key evidence in forensic hydrology. For
agronomists, ecohydrologists, and irrigation engineers, depletion of
temporary static surface water bodies (owing to evaporation and infil­
tration) is also of utmost practical interest because a conjunctive use ofwater resources requires a smart alteration of operation of surface
pumps and well pumps (like during the last 2011–2017 Big Drought in
California). There is a fascinating analogy between soils-sedimentslandforms-hydrology on Mars and in hot deserts of Earth (see, e.g.,
recent projects by space scientists in Oman and Israel: https://oewf.
org/en/portfolio/amadee-20/;
https://oewf.org/en/portfolio/ama
dee-18/ and Groemer et al., 2020).

1

2

The units are indicated solely for matching HYDRUS modeling below.

A.R. Kacimov et al.

Icarus 372 (2022) 114719

2. Analytical solution for falling head Seepage from shallow
axisymmetric ponds into Gardner's and capillarity-free soilsdV= − E(t) − Q(t)dt

(3)

In eq. (3), E(t) > 0 [cm3/h] are evaporation losses, which are pro­portional to the open area of the pond. Consequently,

In an axial cross-section (Fig. 1a), a shallow axisymmetric reservoir isassumed to be filled with water. In an aerial view at t < 0, the ponded
surface represents a disk (Fig. 1b) of a radius R0 [cm], which shrinks
with time, such that the radius of a free water surface DA decreases as 0
< R(t) < R0. The extinction time, Te, when surface water ponding ceases,
i.e., R(Te) = 0, is a part of the solution.
Without any loss of generality, we suppose that there is no freeboard
in the pond at t < 0, i.e., the initial water level OA0 is at the same horizon
as the soil surface. This Section deals with four geometries of the pond's
beds: conical, spherical, spheroidal, and paraboloidal. The soil of thecrater bed is putatively homogeneous, isotropic, and either Gardnerianor capillarity-free. For the former, a dyad of soil hydraulic parameters is(Ks, αP), where Ks [cm/day] is the saturated hydraulic conductivity, andαP [1/cm] is the Philip sorptive number (see, e.g., Angulo-Jaramillo
et al., 2016). Then in eq. (1), K=Ks exp.[αP p], − ∞ < p ≤ 0. For thelatter (a capillarity-free system), only Ks is needed to model seepage
from the reservoir.
We also surmise that the soil extends indefinitely deep at z < 0, i.e.,
we follow Wooding (1968), Weir (1986), and Hunt (1973) in dealing
with a porous half-space (without any confining substrata), into which
seepage flow propagates. Strictly speaking, an unbounded porous halfspace is a reasonable approximation for small-size reservoirs. Indeed,at large scales (on Earth, as well as on other planets and their satellites in
the Solar System, e.g., Titan), the value of Ks decreases with depth (see,e.g., Horvath and Andrews-Hanna, 2021, Horvath et al., 2016, Reh­binder and Isaksson, 1998, Toller and Strack, 2019). Colmation of the
bed is ignored, i.e., the contour CA0 in Fig. 1a is a line along which thepiezometric head is spatially constant at any fixed instance t.

(4)E(t) = e π R2

where e [cm/h] is the evaporation rate (which we assume to be constantat t > 0). For assessment of seepage losses, Q(t) > 0 [cm3/h], we adopt
the Lembke method of successive steady states (PK-62,77), i.e., we as­sume that at any time t > 0 seepage from the pond attains a new steadystate determined by pond's radius R(t). Hydrological balances equations,analogous with ours eqns. (3)–(4), were used in Horvath and AndrewsHanna (2021) and Warren et al. (2021).
In Subsections 2.1–2.3, we approximate Q for our dwindling reser­voir (cone and spherical-spheroidal caps) by the Wooding (1968) for­mula, which is widely used in infiltrometry (see, e.g., Angulo-Jaramillo
et al., 2016):()
4
Q(t) = Ks πR2 (t) 1 +(5)παP R(t)In eq. (5), the pond's depth is assumed to be zero, and, therefore,
Wooding's purely unsaturated seepage is posited. Strictly speaking, forreal reservoirs, a saturated “bulb”, bounded by an “inverted water table”
under the bed of a non-zero depth B(t) should be taken into account (see,e.g., Sophocleous, 2002, Warrick and Zhang, 1987) and the corre­
sponding Q would be higher than that in eq. (5).2
We put eqns. (4) and (5) into eq. (3) and get an initial value problem
for a nonlinear 1st order ODE:
tanω

dR
4Ks
= − e − Ks −,dt
παP R(t)0 ≤ t ≤ Te

(6)

where Te is to be found.
We introduce dimensionless quantities: (R*, B*, z*, r*) = (R, B, z, r)αP,
t* = tKsαP, e* = e/Ks and drop “*” (for the sake of brevity). Then eq. (6)becomes:()
dR
4tanω = − e +1+(7), 0 ≤ t ≤ Te , R0 ≥ R ≥ 0.
dt
π R(t)

2.1. Conical pond
The initial cone's depth, B0 [cm] (Fig. 1a, solid line), obeys theinequality R0 ≫ B0, i.e., a bank slope ω = arctan(B0/R0) is small. A
straight line CAA0 in Fig. 1a, which makes the cone by rotation aroundOz, is z = -B0 + r tan ω.
The origin, O, of a cylindrical coordinate system (z,r,γ) in Fig. 1a is at
the center of the disk r < R0, z = 0; the angular (azimuthal) coordinate γ
vanishes from the analysis.
At t < 0, seepage from the crater was steady-state, constant values of
R0 and B0 were maintained by a constant channel (runoff) inflow into
the crater from its tributaries (Fig. 1b). At t = 0, the pond started toempty due to two draining factors: seepage and evaporation.
We note that Warren et al. (2021) studied a “reverse” process of
extravasation of groundwater through the bottom of a water-filled
Martian crater, albeit taking into account dynamics of axisymmetricsurface water flow hydraulics downslope crater's rim. We ignore evap­oration of liquid pore water into soil-atmosphere vapor, albeit on thetime scale of post-extinction dynamics of a subsurface hydrologicalsystem, evaporation-sublimation becomes an important component of
the Martian hydrological cycle (see, e.g., Baum and Wordsworth, 2020,
Fukushi et al., 2019, Grimm et al., 2017, Lark et al., 2020).
During the time span 0 ≤ t ≤ Te, water is lost from a ponded surface
such thatEq. (7) resembles the Green-Ampt (GA) eq. (PK-62,77). ODE (7)separates and integrates:()
4 log[(4 + π(1 + e)R)/(4 + π(1 + e)R0 )] R − R0t(R) = tanω−.
(8)π(1 + e)2
1+e
Eq. (8) can not be resolved explicitly with respect to R. For GA
infiltration, PK-62,77 elaborated asymptotics R(t), which she derived
from series expansions. We can apply a similar analysis to eq. (8).
Using the ParametricPlot routine of Wolfram's Mathematica
Wolfram (1991), in Fig. 2 (left panel) we plot the curves R(t) for tan ω =0.1, e = 0.1, and R0 = 1, 0.6, and 0.2 (curves 1, 2, 3, respectively). The
right panel demonstrates the extinction time Te as a function of R0 for tan
ω = 0.1 and e = 1.
2.2. Spherical crater
In this subsection, we consider ponds, the beds of which are sphericalcaps. The generating contour is depicted in Fig. 1a by a dashed-line CA0.
A water-filled spherical cap (half of it at an instance t is shaded inside
CA0OC in Fig. 1a) is made of a sphere of a radius Rs. We limit our analysis
with the case of a circular arc, rotated with respect to the Oz axis, but

Lim(R(t) , B(t) ) = (0, 0)t→Tei.e., R(Te) = 0, B(Te) = 0.
The surface water storage in the pond is characterized by the volume
V (cm3) of the cone://
V(t) = πR2 (t)B(t) 3 = πtanωR3 (t) 3(2)

2
Seepage from ponds of non-zero depth will be numerically modeled by
HYDRUS in Section 3.

The principle of mass conservation yields:
3

A.R. Kacimov et al.

Icarus 372 (2022) 114719

Fig. 1. a) Axial cross-section of a crater; b) aerial view of crater's disk on the Earth/Martian surface.

Fig. 2. Radius R(t) of the surface water disk for conical ponds having tan ω = 0.1, and R0 = 1, 0.6, and 0.2 (left panel curves 1, 2, 3, respectively), evaporation rate e= 0.1 (left panel); extinction time Te (R0) for a conical crater with tan ω = 0.1 and e = 1 (right panel).

centered at point O. The radius of the base of the cap at t > 0 is 0 ≤ R(t)≤ R0 ≤ Rs. The case of the center of a generating circular arc, positioned
above or below the Or axis in Fig. 1a, is mathematically acceptable.
Physically, natural reservoirs having R > Rs are hardly possible in loose
unconsolidated Martian or terrestrial sediments/soils because the angles
of repose of these media are never obtuse. For example, in Fig. 1a, thecase of ω > 90o would be physically extravagant.
The volume of the body of revolution of a quarter of a circular arc CA(Fig. 1a) with respect to the Oz axis is:(√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ) (√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ )1
Vs (t) = π Rs −R2s − R2 (t) 2 2Rs + R2s − R2 (t) , 0 ≤ t ≤ Te
(9)
3

dVs
π R3
dR
= √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
dt
R2 − R2 dt

To the dimensionless quantities introduced in Section 2.1 we add Rs*= RsαP (and drop the asterisk).
Then, substituting eq. (10) in the LHS of eq. (3) and using eqns. (4)and (5), we come to the Cauchy problem for the ODE, which governs theextinction of surface water:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
R2s − R2dR
= − (π R(1 + e) + 4 ), R(0) = R0 .
(11)
dt
π R2

Similar to Subsection 2.1 and using eq. (9), the change of water
storage in the shrinking cap is:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
22
R2s − R20 −t(R) = Rs − R −1+e4π(1 + e)2

ArcTanR0√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
R2s − R20 − R R2s − R2
R2s − R20 − R2

(10)

s

We integrate the initial value problem (11) and obtain:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ) ⎤((4 + π (1 + e)R) 4R0 + π (1 + e)R2s + ρ R2s − R20
16
⎥
⎢
(
−
log⎣
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ) ⎦,πρ(1 + e)22
2
2(4 + π(1 + e)R0 ) 4R + π(1 + e)Rs + ρ Rs − R
⎡

4

(12)

A.R. Kacimov et al.

Icarus 372 (2022) 114719

Fig. 3. Shrinking radius of the disk R(t) of a ponded spherical cap for Rs = 1.5, e = 0.1, and R0 =1.5, 1, 0.5 (curves 1, 2, 3, respectively, left panel); extinction time Te(R0) for ponds made by a spherical cap having Rs = 1.5; surface water is subject to evaporation e =1 (right panel).

where ρ =√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
π2 R2s (1 + e)2 − 16. The inequality πRs(1 + e) ≥ 4 must hold

The mass-balance (Cauchy's problem) is:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
b b2 − R2 (t)dR=−((1 +e)πR(t)+4), R(0) = R0 .
dt
πcR2 (t)as a necessary solvability condition of eq. (11). This is a purely mathe­
matical limitation, which pops up from the deficiency of our model,
which in eq. (5) simplifies a real crater (in the approximation of its
infiltration rate) to a flat disk.
The results of computation by eq. (12) are presented in Fig. 3 as R(t)for e = 0.1 and Rs = 1.5, R0 = 1.5, 1, 0.5 (curves 1–3, respectively).
Obviously, the emptying of the hemispherical crater (curve 1) is thelongest. The shallowest spherical cap in Fig. 3 (curve 3) demonstrates
almost no relatively flat part of the R(t) curve. All three curves in Fig. 3
have an almost vertical segment at t close to Te. Graph on the right panelpresents Te (R0) for Rs = 1.5 and e = 1.

Analytical integration of ODE (15) yields:√̅̅̅̅̅̅̅̅̅̅̅̅̅̅ √̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ⎞⎛√̅̅̅̅̅̅̅̅̅̅̅̅̅̅ √̅̅̅̅̅̅̅̅̅̅̅̅̅̅
22b2 − R20
R0 b2 − R20 − R b2 − R2c⎝ b − R −
4
⎠−
t(R)=
ArcTan−b2 − R20 − R2bπ(1+e)2
1+e(⎡√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ) ⎤2ππ+ 4R0 + ρ b2 − R20
(4
+
(1
+
e)R
)
(1
+
e)b
⎥
⎢
16c
)⎥
(
log⎢
⎦,
πρb(1 + e)2 ⎣ (4 + π(1 + e)R0 ) π(1 + e)b2 + 4R + ρ√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
b 2 − R2
(16)

2.3. Spheroidal crater√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
where ρ = π2 b2 (1 + e)2 − 16. A necessary condition for the solvabilityof eq. (15) is πb(1 + e) ≥ 4.
Using eq. (15), in Fig. 4 (left panel) we plotted R(t) for b = 2, c = 1, e
= 0.1, and R0 = 2, 1.5, 1 (curves 1, 2, 3, respectively). Similar to Fig. 3,
the longest extinction and a pronounced relatively flat segment of the R(t) curve in Fig. 4 is for the deepest crater, viz. a hemispheroid (curve 1).
Fig. 4 (right panel) shows Te (R0) plotted for the same b, c, and e.

A spheroidal crater is obtained by the revolution of an ellipse r2/b2 +z /c2 = 1 about the Oz-axis in Fig. 1a. As in the case of a spherical cap in
Section 2.2, the center of the generating ellipse coincides with the originof coordinates, and the slope of the bank of the crater at t = 0 does notexceed 90o. The two semi-axes of this ellipse are b and c. The shrinking
radius DA in Fig. 1a is R(t) ≤ R0 ≤ b. Additional dimensionless quan­tities are: (b*, c*) = (b, c)αP. Obviously, the case of a spherical cap is aparticular case of spheroidal caps.
Similarly to Subsections 2.1–2.2 (the algebra is omitted), we get thefollowing expression for the surface water storage in the pond:)()(√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 2√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
πc
V(t) =b−b2 − R2 (t)(13)
2b + b2 − R2 (t) , 0 ≤ t ≤ Te .
3b2

2.4. Paraboloidal crater
Al-Shukaili et al. (2020) ignored the soil's capillarity in theiranalytical study of the extinction of ponded water from a triangularchannel. In this subsection, we consider an evanescing ponding in ashallow paraboloid and also ignore capillarity. Thus, a generating curve
CA0 in Fig. 1a is a parabola:][z = B0 (r/R0 )2 − 1 ,(17)

Its rate of change, expressed via the time derivative of Wooding's
disk, is:dV
π cR3 (t)dR= √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
dt
b b2 − R2 (t) dt

(15)

(14)

Fig. 4. Shrinking radius R(t) of a ponded spheroid for b = 2, c = 1, e = 0.1, and R0 = 2, 1.5, 1 (left panel, curves 1-3, respectively); extinction time for a spheroidal
pond at b = 2, c = 1, and e = 0.1 (right panel).
5

A.R. Kacimov et al.

Icarus 372 (2022) 114719

where ε = В0/R0 is Hunt's (1973) small parameter. The tip of the parabolais z = -B0. For the freeboard f(t) > 0 in Fig. 1a, we have f(0) = 0 and f(Te) =B0 .
The instantaneous volume of ponded water is− B∫
0 +B(t)− B∫
0 +B(t)r2 (z)dz= πR20

V(t)= π− B0

(1+z/B0 )dz=

πR20 B2 (t) πB0 R4 (t)

− B0

2B0

=

2R20

(18)

Eq. (30) from Hunt (1973) is used to calculate instantaneous seepage
losses from a shrinking paraboloidal pond. In this equation, we selectHunt's (1973) parameters b = c = 0 and a = 1 that in our notations
yields:[/ ]Q(t) = πKs R2 (t) 1 + 1.29B0 R(t) R20(19)
Note that at any instance 0 ≤ t ≤ Te, Hunt's seepage is bounded by aphreatic line AHEH (Fig. 1a), above which soil is instantaneously drained(no capillarity). In other words, the flow domain is a plume, which isbounded by a shrinking line AHC where point AH “slides” to point C as
the time approaches Te. The difference between Hunt's and Wooding's
seepage models is discussed by Philip (1990).
We put eqns. (18) and (19) into eq. (3), keep the same eq. (4) in its
RHS (viz. counting evaporation from the shrinking disk of the free water
surface of our paraboloidal cap), and get:dV 2πB0 R3 (t) dR
=
dt
dt
R20

(20)
Fig. 6. The transport domain for a conical crater with its dimensions, initialand boundary conditions, and finite element mesh.

The dimensionless quantities: (R*, B*, z*, r*) = (R, B, z, r)/R0, Q* = Q/
(KsR02), V* = V/R03, t* = tKs/R0, e* = e/Ks, ε = B0/R0 with omitted
superscripts “*”s are used below. Eq. (19) is valid for relatively small ε.
Substitution of eqns. (20), (19), and (4) into eq. (3) leads to the
Cauchy problem for the following ODE:dR
1 + e + 1.29εR
= −
,
dt
2εR

R(0) = 1

Integration of a GA-type eq. (21) results in:)(
2
1 + e 1 + e + 1.29εRt(R) =
1− R+
log
1.29
1.29ε
1 + e + 1.29εpresent the details for the former body of revolution.
An axial cross-section is shown in Fig. 6. Due to symmetry, only half ofthe flow domain is presented as a pentagon OCA0UEO. The origin of cy­lindrical coordinates Orz is at the bottom left corner of the pentagon. The
sizes of the pentagon, Lh and Lv, are selected large enough for soil water topercolate freely and not to be constrained laterally by an artificial cylin­
drical boundary (a line UE in Fig. 6). We selected Lh = 800 cm and Lv =1000 cm. The following sizes characterize the conical crater: R0 = B0 = 1m, i.e., the bank slope is 45o. The water table at t = 0 is at the top of thetransport domain (at the soil surface) (Fig. 1a), i.e., f(0) = 0.
Soil is made of the van Genuchten's (VG) non-hysteretic loam (see, e.
g., Šimůnek et al., 2018), for which the soil hydraulic properties are(θr, θs, α, n, Ks)= (0.078, 0.43, 0.036 1/cm, 1.56, 24.96 cm/day). In
HYDRUS, we selected the 2D-General type of geometry and 2D-Axisym­
metric Vertical Flow. The total simulation time was TH = 60 days (i.e.,large enough to be greater than Te). The initial condition at t = 0 was p(0,r,z) = − 1000 cm.
The transport domain was discretized into an unstructured finite
element mesh (Fig. 6) with a targeted size of finite elements of 32 cm andthe mesh refinement of 5 cm on the CA0 boundary, which resulted in

(21)

(22)

Obviously, in eq. (22), t(1) = 0 and t(0) = Te.
Fig. 5 illustrates the effect of evaporation: the left panel plots R(t)according to eq. (22) for ε = 0.1 and e = 0.1, 0.4, and 0.7. The right panelof Fig. 5 illustrates the effect of the paraboloid's depth on the extinction
time Te(ε) for e = 0.1.
3. HYDRUS modelingIn this Section, we use HYDRUS (2D/3D) (Šimůnek et al., 2018) for anumerical solution of an axisymmetric Richards' eq. (1) applied to thefollowing shapes of the craters: a cone and paraboloidal cap. Below weFig. 5. The radius of the surface water disk R(t) in a paraboloidal reservoir for ε = 0.1 and e = 0.1, 0.4, and 0.7 (left panel, curves 1-3, respectively); extinction timeTe(ε) in a paraboloidal reservoir for e = 0.1.
6

A.R. Kacimov et al.

Icarus 372 (2022) 114719

2374 finite element nodes, 171 1D elements, and 4575 2D elements(triangles). The default HYDRUS iteration criteria were used. The seg­ments OC and UE were no-flow boundaries, and OD was a free drainage
boundary. We assumed that the segment A0U was also a no-flow
boundary, although evaporation from this moist segment can bemodeled in HYDRUS. Along the A0C segment, a new HYDRUS Reservoir
Boundary Condition (see Šimůnek et al., 2018; Sasidharan et al., 2018)with a Furrow reservoir type was used. Evaporation from the open
water body AD was also neglected.
Fig. 7 (left panel) shows HYDRUS-computed water depth for theconical project depicted in Fig. 6. All water infiltrated in 1.35 days.
Fig. 7 (right panel) presents a drawdown curve for R0 = 100 cm, B0 = 10cm. The HYDRUS extinction time for this shallower cone is 0.25 days. To
compare these results with analytical ones (i.e., eq. (8)), we convert αP to
α by the Ghezzehei et al. (2007) formula αP = 1.3α n that gives for the VG
loam αP = 0.076 1/cm. For the crater in Fig. 7, eq. (8), converted todimensional quantities, gives Te = 2.7 days, while for the bank slope ω =0.1 the analytical value is Te = 0.27 days. Therefore, for a deep cone, theWooding analytical approximation (neglecting the reservoir depth) is
inadequate. For shallow cones, there is a good match between HYDRUS
results and eq. (8).
We notice a steep slope of the B(t) curves in Fig. 7 at t close to Te
when only a small quantity of surface water remains in the cone (AlShukaili et al., 2020 got similar results for triangular channels). This iscaused mainly by the convex shape of the V(B) (water storage-water
depth) function for cones and a much smaller increase in the water
volume versus water depth (dV/dB) when water depths go to zero,compared to other geometrical shapes (spheroids, paraboloids). In our
HYDRUS model, we placed an observation point C (Fig. 6) at the cone's
tip and got the post-extinction limb of the p(t) curve − ∞ < p(t) < 0. In
the range − pc < p(t) < 0, where pc is a soil-specific constant, capillarity
stabilizes the excavation, i.e., the marcescent sediments on the slopes of
the crater are temporarily upheld by moisture suction (too parched andtoo wet soils are less stable against liquefaction, and Aeolian and othertypes of erosion).
The paraboloidal cap was simulated for a much larger domain (Lh =30 m and Lv = 25 m). The origin of cylindrical coordinates D0rz is now at
the top left corner of the transport domain. The paraboloidal cap was
generated using a function Z = br2 for b = 0.08 1/m and the origin of theparaboloidal cap, i.e., point C in Fig. 6, shifted to z = Z – 8 m. In thisexample, the cap's initial water table was 3 m below the soil surface (f(0)= 3 m), i.e., the initial water depth B0 = 5 m. The width of the parab­oloidal cap at the soil surface is R0 = 10 m.
Similarly, as above, the soil is made of the VG non-hysteretic loam.

The total simulation time was TH = 50 days (i.e., large enough to begreater than Te). The initial condition at t = 0 was p(0,r,z) = − 10 m.
The transport domain was discretized into an unstructured finite
element mesh with a targeted size of finite elements of 0.1 m, a
stretching factor of 4 (larger elements in the horizontal direction), themesh refinement of 0.025 m on the CA0 boundary, the mesh refinementof 0.2, 0.25, and 0.2 m in the corner nodes O, E, and U, respectively,
which resulted in 9190 finite element nodes, 376 1D elements, and
18,002 2D elements (triangles). As above, evaporation from the soil
surface or the water table was neglected. The Reservoir Boundary
Condition (see Šimůnek et al., 2018, Sasidharan et al., 2018) with aParaboloidal reservoir type was used along the A0C segment.
Fig. 8 shows the HYDRUS-computed water depth for the parabo­
loidal cap. The extinction time is about 11 days. For comparisons, theanalytical solution, eq. (22), gives for this paraboloid Te = 14.16 days. In
this example, capillarity again increases Q(t) and decreases Te. What areother factors affecting seepage? In analytical models of seepage fromchannels and reservoirs (free boundary problems for the Laplace equa­tion, the Vedernikov-Bouwer model), the variational theorems (see, e.g.,
Ilyinsky et al., 1998) guarantee that if an isobaric drainage layer OE in
Fig. 6 rises, then Q(t) increases. Therefore, we surmise that the HYDRUS“free drainage” condition along OE in Fig. 6 also increases Q(t)Fig. 8. The water depth B as a function of t in a paraboloidal cap simulatedby HYDRUS.

Fig. 7. The water depth B as a function of t in a conical crater simulated by HYDRUS: R0 = B0 = 100 cm (left panel) and R0 =100 cm, B0 = 10 cm (right panel). Note
different scales on axes in the left and right figures.
7A.R. Kacimov et al.

Icarus 372 (2022) 114719

compared with Hunt's solution, for which this condition holds infinitely
deep under the reservoir. On the other hand, the variational theorems(see, e.g., Ilyinsky et al., 1998) also rigorously prove that any artificial
confinement of free-boundary seepage by an external no-flow surface, e.
g., a cylindrical surface r = Lh, 0 < z < Lv in Fig. 6, reduces Q(t)compared with both Hunt's laterally unbounded flow domains. HYDRUS
calculates Q(t) based on the Richards equation, and Wooding's Q would
decrease if any impermeable boundaries confine the porous flow
domain.
Fig. 9 represents the results of HYDRUS computations for the same
paraboloidal crater as in Fig. 8. The panels show a sequence of coloured
maps of θ at five timeinstances: t = 1, 5, 10, 20, and 50 days. One cantrack, for example, evolution of a fully-saturated zone (the darkest red“bubble” bounded by an isohume θ = 0.43). Specifically, at t = 20 days
this zone separates from the bottom of an already empty crater while at t

= 50 days this zone does not exist because capillarity has “devoured” thesaturated “bubble”. That means that a piezometer screened sufficiently
deep under the crater will not detect any soil “wetting” signal p > 0.
4. Quick field experiments for infiltrometry
In October–November 2020, a series of experiments were conductedat the Agricultural Experimental Station (Sultan Qaboos University,
Oman; N: 23◦ 35′ E: 58◦ 9′ ). Excavations were shoveled in dry homo­
geneous dune sand (see the attached photo-video gallery of the experi­ments, Appendix 1). Several freshwater pulses of a temperature of
32–34o C were instantaneously poured into dug craters to stabilize theirshapes and wet the subjacent sand, which is relatively coarse andpermeable. After that, three liters of water were poured into a crater,
and the radius R(t) and the extinction time, Tem, were measured.

Fig. 9. HYDRUS–computed snapshots of the volumetric moisture content in an axial cross section for a paraboloidal cap in Fig. 8; time instances in the panels are
t=1, 5, 10, 20, and 50 days (left->right, up->down).
8

A.R. Kacimov et al.

Icarus 372 (2022) 114719

We present here the results for an initially water-filled crater having
B0 = 6 cm and R0 = 17 cm. For this crater, Tem = 1 min 10 s. For such ashort period of extinction, evaporation from DA in Fig. 1a can beneglected.
Using eq. (22) for this crater shape, we get dimensionless time Te =0.27. Therefore, Kse = R0Te/Tem = 3.9 cm/min. Obviously, adjustment of
Ks to other temperature values is needed using, for instance, the Hazen
formula (PK-77).
Experimental errors are likely due to the crater's shape not being
constant during infiltration. The poured surface water gyrates in thecrater, scours the sand, and deposits it with small rills on the crater bed.
The crater's shape is not exactly paraboloidal, and it is not easy to release
the whole volume V0 into the empty crater. Therefore, the instances t =0 and Tem were not precisely measured. The experimentally found valueof Kse is an overestimate of the real Ks because the used Hunt solution
ignores capillarity.
There are, unfortunately, no analytical formulae for Q, even forsteady-state flows, for axisymmetric seepage from ponds with free sur­faces (capillary fringes or phreatic surfaces without capillarity). We canuse an analogy with capillarity-counting seepage from soil channels (theVedernikov-Bouwer model, see, e.g., Al-Shukaili et al., 2020, PK-62,77,
Mahato and Dey, 2020). Indeed, for a Riesenkampf zero-depth channel
with a width bR and sand having an absolute value of hc of the pressure
head on the capillary fringe, the curve Q/Q0 as a function of hc/b isshown in Fig. 122, page 180, in PK-77. If we assume bR = 34 cm and hc &10 − 20 cm, then Q/Q0 & 1.5 − 2.5, i.e., if steady seepage losses fromthis channel are used for evaluation of Ks, then Kse will be overestimatedby 50–150% for such a small crater. We surmise that for large ephemeral
ponds and Martian craters (see, e.g., Bam et al., 2020), recorded
emptying of these ponds will give Kse much closer to the actual Ks, ascompared with mini-ponds reported in this Section. For terrestrial
ponds, estimated Ks and Philip's sorptive number, αP, or the van Gen­uchten dyad of (α, n), can be compared with HYDRUS simulations (see,e.g., Reynolds, 2015). This is similarly done for standard disc infil­
trometers, which evaluate an allegedly steady-state unsaturated flowfrom a vertical burette through an isobaric discharge disk, with calcu­
lations of the dyad of (Ks,αP) according to the Wooding formula. We also
note that if the geometry of the crater and the whole flow domain are
fixed, then for a constant value of Ks, capillary soils under steady-stateflow conditions always have higher Q(t) than soils without capillarity,
whatever capillarity model (Gardner, van Genuchten, or VedernikovBouwer with free boundaries) is selected. However, suppose soil's
texture changes, i.e., both Ks and capillarity parameters vary. In thatcase, the impact on Q(t) (and therefore Te) is opposite: coarser soils havehigher Ks but milder capillarity.

are conjugated with seasonal evaporation, can be done for ephemeral
ponds like those in Bam et al. (2020).
Our analytical and numerical models can be extended to the case ofheterogeneous porous beds. In large-scale Martian studies of seepage(see, e.g., Horvath and Andrews-Hanna, 2021), the thickness of porous
aquifers subjacent to large craters was several kilometers, and the hy­
draulic properties of these aquifers varied at kilometer scales. PK-62
pointed out that seepage losses of small-size earth channels are signifi­
cantly reduced by colmation, which clogs few centimeters of the nearbed soil. Colmated craters/ponds can be easily modeled by HYDRUS,
which can easily tackle two-component porous composites.
Declaration of Competing Interest
None.
Acknowledgments
This work was supported by the grants DR\RG\17 and IG/CAMS/
SWAE/18/01, Sultan Qaboos University, Oman, and the development
program of the Scientific and Educational Mathematical Center of theVolga Federal District, Russia, grant No. 075-02-2021-1393.
Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.icarus.2021.114719.
References
Al-Shukaili, A., Al-Mayahi, A., Al-Maktoumi, A., Kacimov, A.R., 2020. Unlined trench asa falling head permeameter: analytic and HYDRUS2D modeling versus sandbox
experiment. J. Hydrol. 583, 124568. https://doi.org/10.1016/j.
jhydrol.2020.124568.
Al-Shukaili, A., Al-Busaidi, H., Kacimov, A.R., 2019. Experiments, analytical and
HYDRUS2D modeling of steady jet of quasi-normal surface flow in rectangularchannel coupled with vertical seepage: Vedernikov-Riesenkampf’s legacy revisited.
Adv. Water Resour. 136, 103503.
Angulo-Jaramillo, R., Bagarello, V., Iovino, M., Lassabatere, L., 2016. Infiltration
Measurements for Soil Hydraulic Characterization. Springer, Berlin.
Bam, E.K., Ireson, A.M., van der Kamp, G., Hendry, J.M., 2020. Ephemeral ponds: arethey the dominant source of depression-focused groundwater recharge? Water
Resour. Res. 56 (3), e2019WR026640.
Baum, M., Wordsworth, R., 2020. Groundwater flow to Gale crater in an episodicallywarm climate. J. Geophysical Res. Planets 125 (8), e2020JE006397.
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.
Blank, J.G., Green, S.J., Blake, D., Valley, J.W., Kita, N.T., Treiman, A., Dobson, P.F.,
2009. An alkaline spring system within the Del Puerto Ophiolite (California, USA): a
Mars analog site. Planet. Space Sci. 57 (5–6), 533–540.
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,
infiltration, and groundwater flow. Planet. Space Sci. 171, 17–33.
Conway, S.J., de Haas, T., Harrison, T.N., 2019. Martian gullies: a comprehensive review
of observations, mechanisms and insights from Earth analogues. Geol. Soc. Lond.,
Spec. Publ. 467 (1), 7–66.
Davila, A.F., Warren-Rhodes, K., DiRuggiero, J., 2021. The Atacama Desert: a window
into late Mars surface habitability? In: Soare, R., Conway, S., Williams, J.-P.,
Oehler, D. (Eds.), Mars Geological Enigmas From the Late Noachian Epoch to thePresent Day. Elsevier, pp. 333–356.
Fukushi, K., Sekine, Y., Sakuma, H., Morida, K., Wordsworth, R., 2019. Semiarid climate
and hyposaline lake on early Mars inferred from reconstructed water chemistry at
Gale. Nat. Commun. 10 (1), 1–11.
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.
Goldspiel, J.M., Squyres, S.W., 2011. Groundwater discharge and gully formation on
Martian slopes. Icarus 211 (1), 238–258.
Grimm, R.E., Harrison, K.P., Stillman, D.E., 2014. Water budgets of Martian recurring
slope lineae. Icarus 233, 316–327.
Grimm, R.E., Harrison, K.P., Stillman, D.E., Kirchoff, M.R., 2017. On the secularretention of ground water and ice on Mars. J. Geophysical Res. Planets 122 (1),
94–109.
Groemer, G., Gruber, S., Uebermasser, S., Soucek, A., Lalla, E.A., Lousada, J., Sams, S.,
Sejkora, N., Garnitschnig, S., Sattler, B., Such, P., 2020. The AMADEE-18 Mars
analog expedition in the Dhofar region of Oman. Astrobiology 20 (11), 1276–1286.

5. Concluding remarks
Wray (2021) wrote about the applicability of mathematical andphysical models developed for terrestrial conditions to other planets, inparticular, Mars: “…alternative models often prove relevant in other plan­etary environments even if not in that for which they were first proposed.”
Cicero (see the epigraph) pointed out the link between the last day ofextinction and the ensued change of the locus of an entity, which weinterpret as a human soul or a finite initial volume of water V0, storedhydrostatically at t = 0 in a reservoir. Such a pond can be a miniatureterrestrial dimple on the surface of a crop field where local runoff iscollected after an intense rainfall event (sprinkling irrigation), or a large
Martian crater/lake, which suddenly loses its tributaries (Martianchannels). Our analytical and HYDRUS modeling for saturatedunsaturated seepage from axisymmetric reservoirs predicts how fast
water levels there drop with time. These results will serve to further
optimize the shapes of the reservoirs via the solution of isoperimetricproblems involving dynamics criteria, similar to КасиМов andТартаковский (1993). Also, the analysis of stable isotope enrichment of
surface water and pore water recharging aquifers, when seepage losses
9A.R. Kacimov et al.

Icarus 372 (2022) 114719
Petroff, A.P., Devauchelle, O., Seybold, H., Rothman, D.H., 2013. Bifurcation dynamics
of natural drainage networks. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 371,
20120365, 2004.
Philip, J.R., 1990. How to avoid free boundary problems. In: Hoffman, K.H., Sprekels, J.
(Eds.), Free Boundary Problems: Theory and Applications. Longman, New York,
pp. 193–207.
Philip, J.R., 1992. Falling-head ponded infiltration with evaporation. J. Hydrol. 138(3–4), 591–598.
Polubarinova-Kochina, P.Ya., 1962. 1977. Theory of Ground Water Movement, 2nd ed.
Princeton Univ. Press, Princeton. Nauka, Moscow (in Russian).
Radcliffe, D., Šimůnek, J., 2010. Soil Physics with HYDRUS: Modeling and Applications.
CRC Press, Taylor & Francis Group, Boca Raton.
Rehbinder, G., Isaksson, A., 1998. Large-scale flow of groundwater in Swedish bedrock.
An analytical calculation. Adv. Water Resour. 21 (6), 509–522.
Reynolds, W.D., 2015. A generalized variable-head borehole permeameter analysis for
saturated, unsaturated, rigid or deformable porous media. Eng. Geol. 185, 10–19.
Salese, F., Pondrelli, M., Neeseman, A., Schmidt, G., Ori, G.G., 2019. Geological evidenceof planet-wide groundwater system on Mars. J. Geophysical Res. Planets 124 (2),374–395.
Sasidharan, S.A. Bradford, Š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.
Šimůnek, J., Šejna, M., van Genuchten, M.T., 2018. New features of version 3 of the
HYDRUS (2D/3D) computer software package. J. Hydrol. Hydromechanics 66 (2),
133–142. https://doi.org/10.1515/johh-2017-0050.
Sophocleous, M., 2002. Interactions between groundwater and surface water: the state of
the science. Hydrogeol. J. 10 (1), 52–67.
Strack, O.D.L., 2017. Analytical Groundwater Mechanics. Cambridge Univ, Press, NewYork.
Toller, E.A., Strack, O.D., 2019. Interface flow with vertically varying hydraulicconductivity. Water Resour. Res. 55 (11), 8514–8525.
Warrick, A.W., Zhang, R., 1987. Steady two-and three-dimensional flow from saturated
to unsaturated soil. Adv. Water Resour. 10 (2), 64–68.
Warren, A.O., Holo, S., Kite, E.S., Wilson, S.A., 2021. Overspilling small craters on a dry
Mars: Insights from breach erosion modeling. Earth Planet. Sci. Lett. 554, 116671.
Weir, G.J., 1986. Steady infiltration from large shallow ponds. Water Resour. Res. 22(10), 1462–1468.
Woessner, W.W., 2020. Groundwater-Surface Water Exchange. Groundwater Project,
Guelph.
Wolfram, S., 1991. Mathematica. A System for Doing Mathematics by Computer.
Addison-Wesley, Redwood City.
Wooding, R.A., 1968. Steady infiltration from a shallow circular pond. Water Resour.
Res. 4 (6), 1259–1273.
Wray, J.J., 2021. Contemporary liquid water on Mars? Annu. Rev. Earth Planet. Sci. 49,
141–171.
КасиМов, А.Р., Тартаковский, Д.М., 1993. Оценка вреМени Миграции трассера в потоке
грунтовых вод. Журнал Вычислительной МатеМатики и МатеМатической Физики 33(11), 1751–1759 (in Russian). Engl. translation: Kasimov, A.R., and Tartakovsky, D.
M., 1993. Estimation of tracer migration time in ground water flow. Comput. Math.
Math. Phys., 33(11), 1535–1541.

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.
Horvath, D.G., Andrews-Hanna, J.C., 2021. The hydrology and climate of Mars during
the sedimentary infilling of Gale crater. Earth Planet. Sci. Lett. 568, 117032.
Horvath, D.G., Andrews-Hanna, J.C., Newman, C.E., Mitchell, K.L., Stiles, B.W., 2016.
The influence of subsurface flow on lake formation and north polar lake distribution
on Titan. Icarus 277, 103–124.
Hunt, B.W., 1973. Seepage from shallow reservoir. J. Hydraulics Division ASCE 99 (1),23–30.
Ilyinsky, N.B., Kacimov, A.R., Yakimov, N.D., 1998. Analytical solutions of seepage
theory problems. Inverse methods, variational theorems, optimization and estimates(A review). Fluid Dynam. 33 (2), 157–168.
Kacimov, A., Obnosov, Yu.V., Šimůnek, J., 2020. Seepage to ditches and topographic
depressions in saturated and unsaturated soils. Adv. Water Resour. 145, 103732
https://doi.org/10.1016/j.advwatres.2020.103732.
Kacimov, A.R., Al-Shukaili, A., Al-Mayahi, A., Al-Maktoumi, A., 2021. Triangular ditch of
fastest infiltration into porous substratum. J. Irrigation Drainage ASCE 147 (7),06021002.
Kereszturi, A., Möhlmann, D., Berczi, S., Ganti, T., Horvath, A., Kuti, A., Sik, A.,
Szathmary, E., 2010. Indications of brine related local seepage phenomena on thenorthern hemisphere of Mars. Icarus 207 (1), 149–164.
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.
Lark, L.H., Huber, C., Head, J.W., 2020. Anomalous recurring slope lineae on Mars:
Implications for formation mechanisms. Icarus, p. 114129.
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.
Mahato, R.K., Dey, S., 2020. Hydraulics of seepage from trapezoidal channels.
J. Hydraulic Eng. ASCE (12), 146, 04020083.
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
runoff on Mars. Science 288 (5475), 2330–2335.
Marra, W.A., Braat, L., Baar, A.W., Kleinhans, M.G., 2014. Valley formation bygroundwater seepage, pressurized groundwater outbursts and crater-lake overflow
in flume 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 inexperiments and on Mars. Earth Surface Dynam. 3 (3), 389–408.
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.
Mukherjee, S., Singh, D., Singh, P., Roy, N., 2020. Morphological and morphometricanalysis of a topographic depression near Huygens basin, Mars: identification of aputative endorheic playa. Geomorphology 351, 106912.

10