Journal of Hydrology 598 (2021) 126413

Contents lists available at ScienceDirect

Journal of Hydrology
journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Analytical traveling-wave solutions and HYDRUS modeling of wet wedges
propagating into dry soils: Barenblatt’s regime for Boussinesq’s
equation generalized
A.R. Kacimov a, *, J. Šimůnek bab

Department of Soils, Water and Agricultural Engineering, Sultan Qaboos University, Oman
Department of Environmental Sciences, University of California Riverside, CA, USA

A R T I C L E I N F O

A B S T R A C T

This manuscript was handled by Corrado Cor­
radini, Editor-in-Chief, with the assistance ofStephen Worthington, Associate Editor

The classical Barenblatt solution of an initial-boundary value problem (IBVP) to the parabolic Boussinesq
equation, which gives a rectangular triangle of full saturation, propagating from a reservoir into an adjacentporous bank with a vertical slope, is shown to coincide with a solution of IBVP to the elliptic Laplace equation
with a phreatic surface along which both isobaricity and kinematic conditions are exactly met. For an arbitrarybank slope, a saturated wedge, which propagates (translates) into dry soil, is also explicitly found. The analyticalsolutions favorably compare with the results of HYDRUS-2D modeling, i.e., with the FEM solutions of the same
IBVPs to the Richards equation. Applications to geotechnical engineering of dykes subject to the impact of flash
floods are discussed by comparisons of phreatic lines, loci of the fronts, isobars, equipotential contours, vector
fields of Darcian velocity, isotachs, and streamlines in the three models. For example, it is shown that a rapid
drawup of the reservoir level induces hydraulic gradients, which may cause seepage-induced erosion of theporous medium, in particular, lessivage.

Keywords:
Dry dike subject to seepage from a reservoir
with a water level rising at a constant rate
Similarity solution
Transient complex potential
Translating phreatic surface
Isobars-streamlines-isotachs in HYDRUS

1. Introduction
The water level rise in dam reservoirs, ponds, rivers, lakes, oceans,and other surface water bodies with porous banks has hydrological,
geotechnical, and geopolitical implications. Transient infiltration into
these banks mitigates flash flood waves and smoothens peaks inhydrographs of watercourses (e.g., wadis in Oman and other arid
countries or rivers in humid countries; Denkers, 2021), albeit one usu­
ally does not differentiate “transmission losses” into vertical infiltration
through the bed and banks of unregulated ephemeral streams (see, e.g.,
Noorduijn et al., 2014, Shanafield et al., 2021). The river stage is aboundary condition for bank filtration or infiltration from temporarily
filled ponds in managed aquifer recharge (MAR) systems (Abdeldayem
et al., 2020, Alam et al., 2021, Maliva, 2020). Dikes, levees, and earth
dams, as well as natural soil slopes, when in contact with open water
bodies with rapidly rising water levels, are prone to large-scale collapses
and miniature (local) losses of structural stability of adjacent imbibedsoil massifs (see, e.g., Cheng et al., 2021; Denkers, 2021; Pasetto, 2014;

Yan et al., 2019). During the run-up tide stages, the rising seawater levelhas a cyclostationary hydroecological impact on coastal biosystems(McLachlan and Defeo, 2017).
In this paper, we use analytical solutions and HYDRUS-2D modeling(Šimůnek et al., 2016) to investigate transient seepage from an open
water body with a rising water level into an adjacent slope of a porous
earth body (hereafter called a reservoir and a dike, respectively). Steadyand transient seepage from reservoirs through porous dykes (levees,
earth filled embankments) have been studied both analytically andnumerically as a part of design, construction, operation, maintenance,
and forensic investigation in case of collapses of these earth works (see,e.g., Aravin and Numerov, 1953; Butera et al., 2020; Casagrande, 1937;
Cedergren, 1989; Chahar, 2004; Chapuis and Aubertin, 2002; Desai,
1972; Design, 1987; Fell et al., 2005; Freeze, 1971; Hansen and Rosh­anfekr, 2012; Havenith et al., 2015; Jafari et al., 2019; Jaswal et al.,2020; Jia et al., 2009; Jiao and Post, 2019; Kacimov and Brown, 2015;
Kacimov and Obnosov, 2012; Kacimov and Yakimov, 2001; Kacimov
et al., 2021; Kacimov et al., 2020; Kreuter, 1921; Liggett and Liu, 1979;* Corresponding author at: Department of Soils, Water and Agricultural Engineering, College of Agricultural and Marine Sciences, Sultan Qaboos University, AlKhod 123, PO Box 34, Oman.
E-mail addresses: anvar@squ.edu.om (A.R. Kacimov), jsimunek@ucr.edu (J. Šimůnek).
URL: https://www.squ.edu.om/agriculture/Academic-Department/Soils-Water-and-Agricultural-Engineering (A.R. Kacimov).
https://doi.org/10.1016/j.jhydrol.2021.126413
Received 5 February 2021; Received in revised form 17 April 2021; Accepted 29 April 2021
Available online 4 May 2021
0022-1694/© 2021 Elsevier B.V. All rights reserved.

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413

Fig. 1. a) Barenblatt’s rectangular triangle
as a right-ward propagating saturated zone;b) Barenblatt’s complex potential triangle
propagating left-ward. x0 and y0 are hori­zontal and vertical coordinates, respectively,z(t) = x + i*y is a complex traveling variable,
g is gravitational acceleration, r is the rate of
the water level rise, k is the saturated hy­
draulic conductivity, m is fillable porosity,and w(t) = φ(t) +iψ (t) is a transient complex
potential, where φ = − k h is the velocity
potential and ψ is the stream function.

Masciopinto et al., 1994; Loáiciga, 2020; Mishra and Singh, 2005; Navaset al., 2016; Neuman and Witherspoon, 1970; Nelson-Skornyakov, 1949;
Nichiporovich, 1973; Nouridjanyan et al., 1988; Peter, 1982; Rehbinder,
1997; Soil Mechanics, 1979; Stark et al., 2017; Sun et al., 2017; Tayfur
et al., 2005; The International Levee Handbook, 2013; Yang et al., 2019;
Zhou et al., 2020). Mathematical models of Darcian flows through dikes
should answer the following questions related to the seepage process:modeling of a dyke, having an arbitrary reservoir slope, with anincreasing reservoir water level.
2. Barenblatt’s solution revisited for a dike with a vertical slope
In the 1950ies, Barenblatt (see Aravin and Numerov, 1953; Bare­
nblatt et al., 1984, 1990; Polubarinova-Kochina, 1962, 1977, 1996;
Kochina and Kochina, 1994, 19961) obtained a number of similarity
solutions to the Boussinesq equation, which mathematicians ascribed tothe class of “porous medium equations” (see, e.g., Samarskii et al., 1995;
Vázquez, 2007), which stems from the Polubarinova-Kochina legacy(see, e.g., Ockendon and Howison, 2002). The initial-boundary value
problem (IBVP) to this nonlinear diffusion PDE for seepage studied
below reads:
(
)
∂hs
∂
∂hs
m
hs
, 0⩽x0 ⩽xIf (t),
=k
∂x0
∂t
∂x0
(2.1)
hs (x0 , 0) = 0, 0⩽x0 < ∞,• What is the shape of the moving phreatic surface, and where is thefront’s locus (the tip of this surface)?
• What are the hydraulic gradients in the flow domain, especially intheir “critical zones”?
• Do these gradients exceed a limit of safety against seepage erosion?
To address these questions, we use three different models: theDupuit-Forchheimer model for purely saturated, vertically-averaged
flow with a front (the Boussinesq PDE), the nonlinear potential model
for purely saturated phreatic 2D seepage in a vertical plane (the Laplace
PDE with a moving free boundary), and the saturated–unsaturated flow
model in the same plane (Richards’ PDE without any fronts or freeboundaries). We used these three different models to capture all essen­tial processes (transiency of reservoir’s level as a driver of seepage,gravity, Darcian resistance of dike’s soil and its capillarity) for theconsidered surface water-pore water hydrological system.
In all three models, we assume that seepage is one-phase (watermoisture motion), Darcian, and transient, taking place in a homoge­
neous and isotropic soil. The soil has the saturated hydraulic conduc­tivity k and fillable porosity m. The latter is defined in PolubarinovaKochina (1962, hereafter abbreviated PK-62) and is related to the
HYDRUS soil characteristics in the van Genuchten (hereafter VG) model,
in which m = θs − θr , and where θ s and θ r are the saturated and residual
volumetric water contents, respectively. For transient phreatic (uncon­
fined) seepage in aquifers and earthworks, the main factor, which con­trols the change of pore water storage is the motion of the phreatic
surface (see, e.g., Barenblatt et al., 1984, and PK-62), which is caused byvariations of the boundary conditions. We also assume that both waterand porous skeleton are incompressible.
Although most studies of transient seepage through dykes focus on
the regimes with a drawdown of the reservoir level (most dangerous interms of a rapid collapse of the reservoir-facing slope of the dyke), below
the drawup scenarios (for which the threat of seepage-induced gradientsis subtly pernicious) are modeled. In all problems below, the horizontal
free water level JF (Fig. 1a) in the reservoir rises at a constant rate r, 0 <r = const [m/day]. The engineering novelty of our work is based on thehs (0, t) = r t
where hs(x0,t) is the water table (phreatic surface) height above a hor­izontal bedrock M1L2 (Fig. 1a). A static Cartesian coordinate system x0y0
originates at the triple point O of the contact of initially (t < 0) dry dikesoil, empty reservoir, and bedrock. The abscissa, xIf (t), of the tip, If, of
the wetting front (at which both hs and seepage flux are zero) indicatesthe maximum extent of imbibition into dike’s soil, owing to the rise of
the water level in the reservoir. In the parlance of mathematicians, theexpanding segment OIf is a “finite support” of a “wetting signal” (arectangular triangle Gz in Fig. 1a), which propagates inland with a“finite speed”2.
The Boussinesq equation eliminates (by the Dupuit-Forchheimer
approximation) the vertical coordinate y0, i.e., posits a constant1
In her latest, 1994 and 1996 books in Russian, Polubarinova-Kochina P. Ya.
changed her name to “Kochina P. Ya.”.
2
The first analytical solution to a nonlinear heat conduction problem (theBoussinesq equation belongs to this class) was obtained by Zel’dovich andKompaneets (1950). Finiteness of this speed was historically surprising to those
who worked with nonlinear diffusion and heat conduction. For soil physicists
working with wetting of absolutely dry soils, the “finite speed” was trivial. The
excitement with a mathematically discovered (in the 1950ies) ”finite speeds”
was strange, albeit the “localization” and “stopped fronts” in the so-called“blow-up” regimes (Samarskii et al., 1995) of reservoirs’ filling are indeed
counterintuitive.

2

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413

Fig. 2. a) Saturated tongue propagating into a dry dike with an arbitrary bank slope, b) upper part of the tongue, a traveling-wave wedge, c) Seepage from a reservoir
over a wedge-shaped impermeable bedrock.

piezometric head along vertical lines in Gz and purely horizontal
seepage along the bedrock (Fig. 1a). Barenblatt’s exact solution to IBVP
(2.1) (see, e.g., Barenblatt et al., 1990, PK-62) reads√̅̅̅̅̅̅̅̅̅̅̅̅
hs (t, x0 ) = r t −r m/k xo
(2.2)

The kinematic condition (see, e.g., Dicker and Babu, 1974; Eq.
(2.10), Chapter XV in PK-62; Eq. (3) of Strack in Wang et al., 2011) is:( )2 ( )2
∂ϕ
∂ϕ
∂ϕ
∂ϕ
+
+k
=0
(2.5)
m +
∂t
∂x
∂y
∂y

i.e., the phreatic surface FCIf in Fig. 1a is a segment of a straight tiltedline that caps Gz. This triangle expands into the dike as a traveling wave.
Barenblatt’s flow domain Gz(t) is bounded from the right by a right-ward
moving phreatic surface FCIf.
We show now that Barenblatt’s solution (2.2) is also exact in the full
2D transient potential model, in which the Darcian velocity vector→
V (t, x, y) is not apriori assumed to be horizontal but, generally speaking,has both horizontal, u(t, x, y), and vertical, v(t, x, y), components. To
illustrate this, we select a system of Cartesian traveling-wave co­ordinates xFy, which move vertically upward at the same rate r as thereservoir level, i.e., the traveling abscissa axis always coincides with thisascending ray JF (Fig. 1a).
→We recall that Darcy’s law states V = − k∇h, where the piezometrichead h(t, x, y) for incompressible media obeys Laplace’s equation:

intrinsic
nonlinearity
stems
from
the
term
( )2 ( )2∂ϕ∂ϕ|V| = ∂x + ∂y = u2 +v2 on the LHS of Eq. (2.5). Both Eqs. (2.4)∇2 h(t, x, y) = 0

An2

and (2.5) must hold on an apriori unknown moving interface (FCIr in ourcase, Fig. 1a) between fully saturated (positive pressure) and dry (zeropressure) soil that is another source of nonlinearity. PK-62 discussed
both the full nonlinear model with Eq. (2.5) and its various lineariza­
tions. Nonlinearity makes the IBVP difficult (see Crank, 1984 for other
IBVPs to elliptic PDEs with moving boundaries).
For our complex potential, we select the tip of the in-dyke creeping
saturated wedge, If, as a fiducial point, i.e., a horizontally translating
point in which w = 0 at any time. Using the Cauchy-Riemann conditions,
we get for the Barenblatt wedge in Fig. 1a the following potential (2D)
solution:√̅̅̅̅̅̅̅̅̅̅̅
wB = − r√k ̅̅̅̅̅̅̅̅̅̅̅
t + z r k m,√̅̅̅̅̅̅̅̅̅̅̅
φB (t, x, y) = − r k t + x r k m, ψ B (t, x, y) √= ̅̅̅̅̅̅̅̅̅̅̅
y r k m,VB (z)= ComplexConjugate(dw/dz)√̅̅̅̅̅̅̅̅̅̅̅̅̅̅ = r k m√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
√̅̅̅̅̅̅̅̅̅̅̅
vB = 0, uB = r k m = const, xIf = r k/ m t, tanβ =k/(r m)(2.3)

We follow PK-62 and Strack (1989) and introduce a complex trav­
eling variable z(t) = x + i*y and a transient complex potential w(t, z) =φ(t, x, y) + iψ (t, x, y), where φ = − k h is the velocity potential. Thispotential is harmonic-conjugated with the stream function ψ , which also
obeys Laplace’s equation. Our complex potential w is a holomorphic
function. The complexified Darcian velocity, V(t, z) = u + i * v, is anantiholomorphic function of z, while V = dw/dz (an overbar indicatescomplex conjugation) is holomorphic. The pressure head, p(t, x, y) = h-y,
is another harmonic function, p = P/( ρ g), where P is the gauge porepressure, ρ is water density, and g is gravity acceleration.
Along any moving water table, two boundary conditions must be
satisfied (PK-62; Wang et al., 2011). The first, isobaricity (or zero gauge
pore pressure), condition is:ϕ+ k y = 0(2.6)

Solution (2.6) obeys both the governing PDE (Laplace’s equation forh, ϕ, ψ, and p) inside the “swelling” saturated triangle in Fig. 1a and allboundary conditions on the sides of this triangle, including conditions(2.4) and (2.5).
The transient complex potential domain for solution (2.6), Gw, isanother rectangular triangle depicted in Fig. 1b. Gw is obtained from Gz√̅̅̅̅̅̅̅̅̅̅̅̅
by a scaling factor - r k m and a translation –r k t to the left. Clearly, theimage of Gz in the hodograph plane degenerates into a single point
because the Darcian velocity in solution (2.6) does not depend on z and t.
The phreatic surface and Darcian velocity in solutions (2.6) and (2.2)
are identical. Kacimov et al. (2020) showed that Barenblatt’s solution

(2.4)
3

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413(2.2) is close to what HYDRUS-2D models, provided dike’s soil has arelatively low VG capillarity, and the size of the positive pressure tri­
angle OFIf in Fig. 1a is large enough (that is a common situation ingeotechnical engineering applications if zoned earth-filled embank­ments have reservoir-facing shoulders with a rip-rap or other course
porous cover and the main threat to dike’s stability occurs at the rela­
tively late time of reservoir filling for a relatively high rates r).
We notice that the Boussinesq PDE (and all its linearizations, see
Basha, 2013, PK-62) can model the dynamics of only “normal” phreaticlines (surfaces), for which hs(x) at any t is a function. In many applications
to seepage through dikes and for unconfined aquifers in their near-shore
zone, which are subject to transient drivers of the reservoir-tailwaterlevels and seawater tides-ebbs, this is not the case: neither hs(x) nor the
inversion of hs(x) are - at a certain t - functions in a rigorous mathematicalsense (see illustrative examples in Denkers, 2021; Pasetto, 2014). In otherwords, a part of the phreatic surface is located not above but below an“upper” positive pressure zone, which “overhangs” a dry soil zone, whichin its own term “overhangs” another segment of a “normal” phreatic
surface, which caps a “lower” positive pressure zone. Consequently, awellbore (piezometer) screened across the whole depth of such a homo­
geneous (!) dyke aquifer would detect a transient inflow through several
horizons of a porous medium that is usually attributed to stratification(heterogeneity) in hydraulic properties. Such fancy flow topology, withdry soil zones sandwiched between “zigzagged” or “fingered” saturated
zones, may seem weird to groundwater hydrologists indoctrinated by“standard models” and software packages like MODFLOW.

Fig. 3. Functions β(ω) for the VG (HYDRUS-default) loam having m = 0.352
and r/k = 0.04, 4.0, 32, and 128 (curves 1–4, respectively) calculated from
Eq. (3.1).

Similarly to Section 2 above, we introduce the piezometric (total)head h(t, x, y) and stream function ψ(t, x, y). Both satisfy the Laplace
equation inside I1BFI2. The streamlines (three are shown as boldarrowed segments in Fig. 2b) of Darcian flow from the reservoir into thesoil are straight and perpendicular to the slope along the whole segment
FBI1 (except, of course, two singular endpoints, F and I1). In the “simi­
larity zone” (Fig. 2b), the constant total head lines (two are dashed inFig. 2b) are parallel to the slope.
Although the whole problem is essentially transient due to the rise of
the reservoir water level, the traveling wave solution is surprisingly
simple. In this solution, the piezometric head (as well as the stream
function and the pressure head) exactly satisfies the governing (Lap­lace’s) PDE inside the expanding ”seepage wedge” I1BFI2. The solution
also meets the boundary conditions along its both rays, which bound thewedge, viz. along I1BF and I2CF (see Appendix I) :√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
sin2ω + sin2 2ω − 8m r/kcos2ω + 8m r/kr mtanββ = arctan
, V=
sinω
4m r/k
(3.1)

3. A constant-rate rise of a water level for a dike with anarbitrary slope
In this Section, we generalize Barenblatt’s solution for a vertical bankto an arbitrary bank slope.
3.1. Analytical solution for saturated seepage with a straight tilted watertable propagating in soils without capillarity
In this Subsection, we study a large size dike (Fig. 2a). Its slope
I1BFU, adjacent to the reservoir, is slanted at an angle ω,0 ≤ ω ≤ π/2, tothe horizontal dike bed M2M1, which can be impermeable or highly
permeable (these two cases will be examined later by HYDRUS). We
study a transient seepage zone in the vicinity of the “climbing” tip F of
the wetting front (phreatic surface) FI2M2. In other words, we ignore
what happens near the segment M2M1 where perturbation of the flow by
the dike bed makes seepage essentially two-dimensional (2D) and sim­
ilarity conditions break. We also assume that the dike is broad enough, i.
e., its right boundary, schematically depicted by a dashed line L1L2 in
Fig. 2a, does not affect seepage near the slope. In Section 3.2, we willmodel by HYDRUS-2D dikes of a finite horizontal width.
In this Subsection, we ignore capillarity, which, we repeat, is
acceptable in relatively coarse earth works such as blitz-constructed (e.
g., by emergency bulldozing) embankments made of local topsoil (often
sand or local alluvium) for exigency protection against flash floods. In
gravitational earth dams, the dam shoulders (shells), buttressing the clay
core, are also commonly coarse (Kacimov and Obnosov, 2012).
Thus, we obtain a similarity solution for an analytic element zoomed
in Fig. 2b. Here, at the dike slope I1BF, we select a water-submerged
point B and a system of traveling Cartesian coordinates xFy. Point B is
the tip of the propagating surface and pore water at the time instance t =-r T, where T is an arbitrary time instance. At the snapshot time t = 0,
depicted in Fig. 2b, the tip F is the point where air, reservoir water, anddry and saturated soils make four conterminous wedges.
The ponded water level above B is H(T) = r T. Along the phreatic
surface FCM2, the pressure is atmospheric (zero). In the “similarity
zone” of Fig. 2b, FCI2, this surface is a straight line at an angle β, 0 < β <π/2, with the slope M1I1BF. The β value is to be found from the solution
below. The wedge I1BFI2 is a fully saturated flow domain.

where V is the magnitude of the Darcian velocity vector (constant in
the whole saturated wedge). In the limit of a vertical slope, ω = π/2, Eq.
(3.1) is reduced to the Barenblatt similarity solution. Fig. 3 shows
functions β(ω) calculated using Eq. (3.1) for the VG loam from the
HYDRUS Soil Catalogue having m = 0.352 (0.43–0.078) and different
ratios of r/k.
PK-62 considered the inequality in the hydraulic gradient, Ig = V/k >1, as a criterion of seepage-induced instability of the porous skeleton,
which is a key design criterion for earthen dams (De Mello, 1977). Eq.
(3.1) allows one to immediately check this criterion in the similarity
zone of Fig. 2b and, hence, along the whole slope M1F in Fig. 2a, because
Darcian velocity along M1I1 is higher than along I1F.
Capillarity can be added to the above used analytical solutions in thesame manner as in Chen and Young (2006) by imposing the condition p= -pc (a constant pressure head along I1F, a capillary fringe boundary in
the Green-Ampt model).
3.2. HYDRUS FEM modelingIn this subsection, we use HYDRUS-2D and compare this FEM code’s
results with the analytical solution above. There is no free surface in the
HYDRUS domain, and the Richards equation is solved in the entire flow
domain (Šimůnek et al., 2016). This domain is finite, i.e., infinite wedges
assumed in the analytical solutions have to be bounded, i.e., the domain
has to be large enough so that the external boundaries do not affect thesolution. We recall that semi-infinite domains, i.e., a quadrant and ahalf-plane, were used in the analytical solutions for seepage flows in
4

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413

instance. Far above M2 in the saturated tongue, these vectors have analmost constant magnitude (0.62 m/day) and the same direction,
showing a slight downward dip. For comparisons, Eqs. (2.6) and (3.1)
give 0.84 m/day and 0.85 m/day, respectively. The hydraulic gradientsin all three models are much higher than 1. Moreover, HYDRUS (Fig. 5b)gives even higher values of the hydraulic gradient along the slope close
to point F as one moves from below. In the vicinity of this point, HYD­
RUS velocity has a blip (maximum at F) due to increased inflow in bothhorizontal and upward direction into the initially dry dyke and is zeroalong FU.
We emphasize that the similarity solution (3.1) is not valid suffi­
ciently close to the bed of the dyke in Fig. 2a. However, (3.1) gives exact
V for a rather extravagant geometry shown in Fig. 2c. Here, an imper­
meable boundary is made of two rays: BM0 and BC. The ray BU bounds
from above a “tilted quadrant”, UBC, of initially dry soil. The reservoir,
the wedge M0BU, is initially empty. The water level there starts risingwith a constant rate r at t = 0. As a result, a positive-pressure rectangulartriangle FBIf propagates into the “tilted quadrant” in exactly the samemanner as the Barenblatt triangle in Fig. 1a. Seepage is 1D along thetilted bedrock BC in Fig. 2c (see Youngs, 1974, his Fig. 3 for an analo­
gous steady-state flow). Now we apply a variational theorem (Ilyinsky
et al., 1998) to the Laplace-equation governed seepage with a moving
phreatic surface (no capillarity). Specifically, let us “press-in” (rotate
counter-clockwise) the tilted ray BC of Fig. 2c into the horizontal ray
M1L2 in Fig. 2a. According to the variational theorem, this deformation
reduces the Darcian velocity along M1L2 (Fig. 2a) compared to its con­stant value along BC (Fig. 2c). This reduction is significant near point M1(B) and minor near point F. Consequently, the instantaneous discharges
obey the inequality:∫F∫F
QM1 F (t) =VM1 F (t, s)ds < QBF (t) =VBF (t)ds = VBF (t)|BF|(3.2)

Fig. 4. HYDRUS-2D trapezoidal flow domain.

Fig. 1a and 2b, respectively. We select a HYDRUS trapezoid shown in
Fig. 4. The VG soil hydraulic functions and parameters are taken fromthe HYDRUS Soil Catalogue. We recall that the total head h = p + Z. The
notations h and p are different from the HYDRUS manual, where p stands
for the pore pressure head (modulo matrix potential).
We use X and Z3 for the horizontal and vertical coordinates (fixed) in
our HYDRUS simulations and corresponding figures, with the origin atpoint M1.
We use most of the default HYDRUS options, such as iteration criteria
and time step controls. The initial (i.e., at t = 0, t is the HYDRUS
simulation time, 0 < t < 10 days) pressure head p = -100 cm throughout
the modeled 2D domain (i.e., we assume an initially relatively dry dikesoil).
The HYDRUS domain size in Fig. 4 is indicated by the X and Z co­
ordinates of points U and L1. The reservoir bank slope is 10:1. The dyke
is made of the VG loam with k = 0.25 m/day, θs = 0.43, θr = 0.078, andthe VG constants α = 3.6 1/m and n = 1.56. The spatial discretization of
the transport domain with an element size of 0.5 m involves 23,171
nodes and 581 1D and 45,759 2D finite elements. Water starts rising in
the dike reservoir from the trapezoid’s impermeable bottom in Fig. 4
with the rate of r = 8 m/day. We assume that the horizontal dike crest
UL1 and the vertical tailwater slope L1L2 are impermeable, i.e., noevaporation occurs from these segments. This condition can be changed
to the atmospheric one in HYDRUS. However, even in hyperarid Oman,evaporation is much smaller than the effect of the rapid drawup of thereservoir level, and the zero flux condition is thus a good approximation.
Along M1FU, a time-variable boundary condition with a gradually
increasing water level is imposed for the pressure head. The segment FU
is assumed to be impermeable, using one of the HYDRUS special
boundary options imposing no flux on the section of a boundary with anegative pressure head.
First, we assume that M1M2L2 is impermeable (as in Section 2). In
Fig. 5a, isobaric contour lines are plotted for t = 10 days. Four contours
are highlighted. The contour line p = 0 (θ = 0.43) is sandwiched be­
tween the lines p = 4 m and p = -1 m and indicates the locus of the
propagating phreatic surface. All isobars, located sufficiently far fromM2, are almost parallel to the phreatic surface. The abscissa of M2 is 33.2
m. For comparisons, the Barenblatt solution (2.6) gives 23.84 m andfrom Eq. (3.1) we get xM2 = yF sin β/[sinω sin(β + ω)] = 23.6 m. Further
propagation of the phreatic surface in the HYDRUS simulation comparedto the analytical solutions is due to the capillary forces considered by thenumerical model and neglected by the analytical solutions. This
discrepancy will be larger for fine-textured soils than for coarse-textured
soils and reflects the relative importance of gravitational and capillary
forces characterized by the soil capillary length.
Fig. 5b illustrates simulated vectors of Darcian velocity at the same

M1

B

The HYDRUS simulations (Fig. 5b) qualitatively concur with thevariational theorem. The inequality (3.2) is useful for coarse soils (minor
capillarity) because it gives an upper bound for QM1 F , which requires a
numerical solution to the Laplace equation with a moving phreatic
surface via an elementary calculated QBF. It is noteworthy that, for theexample in Fig. 5b, the whole slope M1F is prone to seepage erosion(according to the PK-62 criterion) induced by the drawup. Hence, fine
fractions of the loam can be agitated, entrained by seepage, and trans­
located into the dyke interior (the phenomenon of lessivage).
When the dyke is made of the VG sand, i.e., from a much more
permeable (k = 7.13 m/day) porous medium with smaller capillarity (α= 14.5 1/m and n = 2.68), then the width of the HYDRUS trapezoid has
to be made larger than for loam (see Fig. 5c). This forestalls the impact of
the tailwater (rightmost) boundary on the encroaching saturated
tongue. At t = 10 days, HYDRUS gives XM2 = 123 m, and the velocity
magnitude ranges 4.7 m/day < V < 5.1 m/day along the ponded
reservoir slope M1F in the range 50 m < Z < 80 m. Assuming sand’s
fillable porosity m = 0.43–0.045, the Barenblatt solution (2.6) gives XM2= 121.72 m, V = u = 4.68 m/day, while Eq. (3.1) yields XM2 = 113.4 m,
V = 5.05 m/day. As discussed above, since the extent of capillarity issignificantly lower in the sand than in the loam, the numerical andanalytical results’ correspondence is much higher in the sand. It isnoteworthy that the Barenblatt solution for a vertical slope is a good
approximation for a sufficiently steep non-vertical reservoir slope, pro­vided the rate of drawup in the reservoir is large enough.
Next, we considered the effect of the boundary condition along
M1M2, i.e., along the bottom of the “saturated tongue” in Fig. 2a. Insteadof a no-flow condition, we applied the HYDRUS free drainage (a unit
vertical hydraulic gradient) condition along M1M2 for the same loamy
trapezoid as shown in Fig. 4 and the same r. The results are shown in
Fig. 5d, e. At the end of the simulation, the left flux QM1 F (10) is about 48
m2/day for the case with the impermeable bottom and 53 m2/day for thepermeable bottom, while the free drainage flux in this latter case is3
Do not confuse the complex variable z(t, x, y) in Section 2 with the verticalcoordinate Z in HYDRUS.

5

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413

Fig. 5. Pressure heads (a, c, d), the vector field of Darcian velocities (b), and streamlines (e) at t = 10 days for r = 8 m/day in a trapezoidal dyke made of the VG loam(a, b, d, e) or sand (c). The bed of the dyke is either impermeable (a–c) or permeable (d, e) with a free drainage boundary condition.

about 7 m2/d. Thus, only about 46 m2/d contributes to the change instorage (very similar to 48 m2/d for the impermeable bottom), which
explains a relatively small difference in the front’s advancement. For thecase of an impermeable bed of the dyke, the water storage added by
infiltration of reservoir water is defined as:∫tS(t) =QM1 F (u)du(3.3)

S(t) = r3/2

√̅̅̅̅̅̅̅̅̅ 2
k m t /2

(3.4)

One can see from Fig. 5d, e that the boundary condition at M1M2 hasa relatively small impact on the “traveling wave” zone near point F. For
example, at t = 10 day, the abscissa of the HYDRUS tip M2 is XM2 = 33.2
m and 28.2 m for the impermeable and permeable beds of the dyke,
respectively. Along the ponded slope FM1 in the range 50 m < Z < 80 m,
the HYDRUS velocity magnitude ranges 0.63 m/day < V < 0.633 m/day
for either boundary condition along M1L2. However, the boundary’s

0

For the Barenblatt solution (2.6), Eq. (3.3) gives:
6A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413

Fig. 6. Pressure heads, p, at t = 10 days for r = 4 m/day and an impermeable bed of a triangular dyke made of the VG sand (a) and loam (b).

impact is significantly larger on the velocity magnitude in the lower partof the inflow boundary.
Next, we selected a triangular dike with another, milder bank slope
1:1 (ω = π/4), r = 4 m/day, a permeable (free drainage) bed M1L2, andcompared VG’s loam and sand as porous media. Fig. 6a shows thepressure head contours for a dike made of the VG sand at t = 10 days andr = 4 m/day. The contour lines between 5 m and − 1 m are plotted tocapture the water table sandwiched in between. HYDRUS gives the front
position XM2 = 52 m. Along the ponded slope FM1 in the range 30 m < Z< 40 m, the HYDRUS velocity magnitude ranges 5.34 m/day < V < 7.36
m/day (decreasing to the dyke crest F). Eq. (3.1) yields a pretty close
value of XM2 = 51.43 m and the velocity V = 6.68 m/day.
Fig. 6b presents the same contour lines for a dike made of the VG
loam at the same time instance. As one can see, the less permeable soil
results in an “inverted water table” (Sophocleous, 2002). Obviously, thisinverted water table indicates a positive pressure wedge spreading
above dry soil (an “unstable” situation prone to vertical “fingering”). In
the analytical model (Subsection 3.1), the condition of such a seepage
regime is ω + β < π/2. For r/k = 32 and the VG loam, from Eq. (3.1) we
obtain that if ω < 1.281, the water table will overhang a dry soil beneaththe saturated wedge.
Let us return to a loamy trapezoid on an impermeable base (Fig. 5).
We placed four observation nodes: the first one at the bottom of thereservoir (X = 0, Z = 0), the second at the bottom of the dyke (22.5 m, 0),
and the third (22.5, 30.2) m and fourth (22.5, 75.2) just above thesecond and third. Fig. 7 shows p(t) for the first three observation nodes.
This figure shows that the water table and the ensued positive pressure
head zone arrive at observation nodes 2 and 3 at t = 6.7 days and 9.8
days, respectively. The front does not arrive at all (in 10 days of simu­
lations) at observation node 4.
Pressure tracking at specified HYDRUS points can be used for thedesign of screens’ length and depth of vertical (horizontal, directional)
wells to be drilled (dug) in porous banks of ephemeral streams (wadis, inOman). Indeed, local Omani farmers tap ephemeral perched aquifers
adjacent to wadis and reservoirs of recharge dams. The main questions
are: Where and how deep to drill (dig) a well? Shallower wells arecheaper but may not intercept the encroached (essentially “ephem­eral”4) saturated tongue. For example, a well screened to a depth ofobservation node 4 will be redundant (for the flash flood simulated inFig. 5b, c).

4
Obviously, a constant-rate water level rise in the reservoir will cease at a
certain time. After that, drawdown of the reservoir water level will commence.
The saturated triangle will slump and a “dry” observational point may remain
dry forever.

7

Journal of Hydrology 598 (2021) 126413

A.R. Kacimov and J. Šimůnek

Fig. 7. Pressure heads p(t) at three observation nodes [(0,0), (0,22.5), (22.5, 30.2)] in a loamy trapezoidal dyke (Fig. 5a, b).

We simulated other reservoir slopes ω, VG soils, drawup rates r, andother boundary conditions along M1M2.

2) VG = Van Genuchten
CRediT authorship contribution statement

4. Concluding remarks

A.R. Kacimov: Conceptualization, Methodology, Software, Formal
analysis, Resources, Writing - review & editing, Visualization, Project
administration, Funding acquisition. J. Šimůnek: Conceptualization,
Methodology, Software, Formal analysis, Resources, Writing - review &editing, Visualization.

Three models (the parabolic 1D Boussinesq equation, the elliptic 2D
Laplace equation with a moving phreatic surface, and the parabolic 2D
Richards equation), solved analytically and numerically (by HYDRUS2D), concur in their solutions of IBVPs if the water level in a reservoir
rises with a sufficiently high but constant rate and the imbibed soil is
initially dry. Specifically, the shape of the saturated zone advancing into
the dyke near the rising free water-soil point (tip F in Fig. 1a and 2) is a
traveling wedge, where Darcian velocities, and therefore the hydraulicgradients and threats of lessivage, are pretty close in all three models.
The intricacy of Darcian flows with fronts and phreatic surfaces“overhanging” dry soil zones should precaution modelers not to be
obfuscated by “standard hydrological concepts”.
There are Abbreviations

Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.
Acknowledgments

1) PK-62 = reference to Polubarinova-Kochina, P.Ya., 1962. Theory of
Ground Water Movement. Princeton University Press, Princeton +second edition of the book in Russian (published in 1977)

This work was supported by the Sultan Qaboos University via thegrant DR/RG/17. Helpful comments of three anonymous referees areappreciated.

Appendix I
In this Appendix, we derive Eq. (3.1).
We consider an arbitrary rectilinear streamline BC in Fig. 2b, which is far above the dike foundation M1M2 (Fig. 2a). Close to the foundation,
seepage is 2D, and I2M2 is not a straight line.
Obviously, from the rectangular triangle ABF, the hypotenuse |BF|=|AB|/sin ωπ. From the condition of a constant rate r, |AB|=r T.
We now select point F in Fig. 2b as fiducial, i.e., the total head h(t, x, y) is counted from F, viz.
h(0,0,0) = 0.
The pressure head p = h-y. At point B, pB = r T, yB = -r T, and hB = 0. From the right triangle AFB, we have:|BC| = |BF|tanβ =

rT
tanβ
sinω

(A1)

From the right triangle BCN:
yC = − r T − |BC|cosω = − r T −r Ttanβ
tanω

(A2)

8

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413

At point C on the phreatic line, pC = 0, and hence:
hC = yC = − r T −r Ttanβ
tanω

(A3)

→
Eqs. (A1), (A3), and Dacry’s law (we recall V = − k∇h = ∇φ) for the straight (but time-stretching) streamline BC yield:()
hB − hCsinω
|VBC | = k=k+ cosω
|BC|
tanβ

(A4)

From the conservation of mass (we recall that water and soil are posited to be incompressible), we have:
VC = m

d|BC|
tanβ
= mr
dT
sinω

(A5)

Eq. (A5) is a kinematic boundary condition for this 1D “titled” infiltration along BC. In the same form, this condition is used in the Green-Amptmodel of 1D vertical infiltration (see, e.g., PK-62).
We equate the right-hand sides of Eqs. (A4) and (A5) and get an equation:(A6)

m r/ktan2 β − tanβ sinωcosω − sin2 ω = 0

Eq. (A6) is a quadratic equation with respect to tan β. We retain only a positive root of this equation that, together with Eq. (A5), gives Eq. (3.1).

References

Jafari, Navid H., Cadigan, Jack A., Stark, Timothy D., Woodward, Mark L., 2019.
Phreatic surface migration through an unsaturated levee embankment. J. Geotech.
Geoenviron. Eng. 145 (11), 05019010. https://doi.org/10.1061/(ASCE)GT.19435606.0002139.
Jaswal, Mamta, Sinha, Rabindra Kumar, Sen, Phalguni, 2020. Delineation of phreatic
surface in soil type slope—A comparative study using physical and numericalmodeling. J. Min. Sci. 56 (3), 494–504.
Jia, G.W., Zhan, Tony L.T., Chen, Y.M., Fredlund, D.G., 2009. Performance of a largescale slope model subjected to rising and lowering water levels. Eng. Geol. 106 (1-2),
92–103.
Jiao, J., Post, V., 2019. Coastal Hydrogeology. Cambridge University Press.
Kacimov, A.R., Al-Maktoumi, A., Obnosov, Yu.V., 2021. Seepage through earth dam withclay core and toe drain: the Casagrande-Numerov analytical legacy revisited. in press
ISH J. Hydraul. Eng., 1–9. https://doi.org/10.1080/09715010.2019.1633694.
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, Anvar, Obnosov, Yurii, 2012. Analytical solutions for seepage near materialboundaries in dam cores: the Davison-Kalinin problems revisited. Appl. Math.
Model. 36 (3), 1286–1301. https://doi.org/10.1016/j.apm.2011.07.088.
Kacimov, A.R., Yakimov, N.D., 2001. Moving phreatic surface in a porous slab: An
analytical solution. J. Eng. Math. 2001 (40), 399–411.
Kacimov, A.R., Yakimov, N.D., Šimůnek, J., 2020. Phreatic seepage flow through anearth dam with an impeding strip. Comput. Geosci. 24 (1), 17–35.
Kochina, P.Ya., Kochina, N.N., 1994. Hydromechanics of Groundwater and Problems of
Irrigation. Fizmatlit, Moscow, p. 240 (in Russian).
Kochina, P.Ya., Kochina, N.N., 1996. Problems of Motion With Free Boundaries in
Subsurface Hydrodynamics. Published by the Editorial Board of the Journal“Uspekhi Fizicheskikh Nauk”, Moscow, 176 p. (in Russian).
Kreuter, F., 1921. Der Flussbau. Handbuch der Ingenieurwissenschaften 3(6).
Engelmann, Leipzig, Germany (in German).
Liggett, James A., Liu, Philip L.-F., 1979. Unsteady interzonal free surface flow in porous
media. Water Resour. Res. 15 (2), 240–246.
Loáiciga, Hugo A., 2020. Seepage face in steady-state groundwater flow between two
water bodies. J. Hydrol. Eng. ASCE 25 (9), 06020005. https://doi.org/10.1061/(ASCE)HE.1943-5584.0001997.
Maliva, R.G., 2020. Anthropogenic Aquifer Recharge. Springer, Cham.
Masciopinto, C., Passarella, G., Vurro, M., Castellano, L., 1994. Numerical simulations forthe evaluation of the free surface history in porous media. Comparison between two
different approaches. Adv. Eng. Softw. 21 (3), 149–157.
McLachlan, A., Defeo, O., 2017. The Ecology of Sandy Shores, 3-rd ed. Academic Press.
Mishra, G.C., Singh, A.K., 2005. Seepage through a levee. Int. J. Geomech. ASCE 5 (1),74–79.
Navas, Pedro, López-Querol, Susana, Yu, Rena C., Li, Bo, 2016. B-bar based algorithm
applied to meshfree numerical schemes to solve unconfined seepage problems
through porous media. Int. J. Numer. Anal. Meth. Geomecha. 40 (6), 962–984.
Neuman, Shlomo P., Witherspoon, Paul A., 1970. Finite element method of analyzing
steady seepage with a free surface. Water Resour. Res. 6 (3), 889–897.
Nelson-Skornyakov, F.B. (1949). Seepage in Homogeneous Medium. Soviet Science,
Moscow (in Russian).
Nichiporovich, A.A., (1973). Dams of Local Materials. Strojizdat, Moscow (in Russian).
Noorduijn, Saskia L., Shanafield, Margaret, Trigg, Mark A., Harrington, Glenn A.,
Cook, Peter G., Peeters, L., 2014. Estimating seepage flux from ephemeral stream
channels using surface water and groundwater level data. Water Resour. Res. 50 (2),1474–1489.
Nouridjanyan, S.Sh., Sarkisian, V.S., Gabayan, G.S., 1988. Filtration from reservoirs
during water level fluctuations. Proc. Acad. Sci. Armenian SSR, Earth Sci. XLI (2),64–71 in Russian.
Ockendon, J.R., Howison, S.D., 2002. Kochina and Hele-Shaw in modern mathematics,
natural science and industry. J. Appl. Math. Mech. 66 (3), 505–512.

Abdeldayem, O.M., Eldaghar, O., Mostafa, K., Habashy, M., Hassan, A.A., Mahmoud, H.,
Morsy, K., Abdelrady, A., Peters, R.W., 2020. Mitigation plan and water harvestingof flashflood in arid rural communities using modelling approach: A case study in
Afouna village, Egypt. Water 12 (9), 2565.
Alam, S., Borthakur, A., Ravi, S., Gebremichael, M., Mohanty, S., 2021. Managed aquifer
recharge implementation criteria to achieve water sustainability. Sci. Total Enviro.
144992.
Aravin, V.I., Numerov, S.N., 1953. Theory of Fluid Flow in Undeformable Porous Media.
Gostekhizdat, Moscow (in Russian). English Translation: Israel Program for Scientific
Translation, Jerusalem, 1965.
Basha, H.A., 2013. Traveling wave solution of the Boussinesq equation for groundwaterflow in horizontal aquifers. Water Resour. Res. 49 (3), 1668–1679.
Barenblatt, G.I., Entov, V.M., Ryzhik, V.M., 1984. Motion of Liquids and Gases in Natural
Formations. Nedra, Moscow (in Russian).
Barenblatt, G.I., Entov, V.M., Ryzhik, V.M., 1990. Theory of Fluid Flows Through NaturalRocks. Kluwer, Dordrecht.
Butera, Ilaria, Climaci, Marco, Tanda, Maria Giovanna, 2020. Numerical analysis ofphreatic levels in river embankments due to flood events. J. Hydrol. 590, 125382.
https://doi.org/10.1016/j.jhydrol.2020.125382.
Casagrande, A., 1937. Seepage through dams. J. New England Water Works 51, 295–336.
Cedergren, H.R., 1989. Seepage, Drainage and Flow Nets. Wiley, New York.
Chahar, Bhagu R., 2004. Determination of length of a horizontal drain in homogeneous
earth dams. J. Irrigation Drainage Eng., ASCE 130 (6), 530–536.
Chapuis, Robert P., Aubertin, Michel, 2002. A simplified method to estimate saturated
and unsaturated seepage through dikes under steady-state conditions. Can. Geotech.
J. 38 (6), 1321–1328.
Chen, L., Young, M.H., 2006. Green-Ampt infiltration model for sloping surfaces. Water
Resour. Res. 42, W07420. https://doi.org/10.1029/2005WR004468.
Cheng, Lin, Zhang, Anan, Cao, Buliang, Yang, Jie, Hu, Liangcai, Li, Yanlong, 2021. Anexperimental study on monitoring the phreatic line of an embankment dam based on
temperature detection by OFDR. Opt. Fiber Technol. 63, 102510. https://doi.org/
10.1016/j.yofte.2021.102510.
Crank, J., 1984. Free and Moving Boundary Problems. Clarendon Press, Oxford.
De Mello, V.F.B., 1977. Reflections on design decisions of practical significance toembankment dams. Géotechnique 27 (3), 281–355.
Denkers, A.J., 2021. Analysis of the Phreatic Line of Various Types of Primary Flood
Defences: Is the Phreatic Line Schematisation too conservative? BSc thesis Universityof Twente, Holland.
Desai, Chandrakant S., 1972. Seepage analysis of earth banks under drawdown. J. Soil
Mech. Found. Eng. Div. ASCE 98 (11), 1143–1162.
Design of Small Dams, 1987. USBR-ISDI, Washington.
Dicker, D., Babu, D.K., 1974. Two dimensional seepage in layered soil-destabilizing
effects of flows with an unsteady free surface. Water Resour. Res. 10 (4), 801–809.
Fell, R., MacGregor, P., Stapledon, D., Bell, G., 2005. Geotechnical Engineering of Dams.
Balkema, Leiden.
Freeze, R. Allan, 1971. Influence of the unsaturated flow domain on seepage throughearth dams. Water Resour. Res. 7 (4), 929–941.
Hansen, David, Roshanfekr, Ali, 2012. Assessment of potential for seepage-induced
unraveling failure of flow-through rockfill dams. Int. J. Geomech. ASCE 12 (5),560–573.
Havenith, H.B., Torgoev, I., Torgoev, A., Strom, A., Yuan, Xu.Y., Fernandez-Steeger, T.,2015. The Kambarata 2 blast-fill dam, Kyrgyz Republic: blast event, geophysicalmonitoring and dam structure modelling. Geoenvironmental. Disasters 2 (11).
https://doi.org/10.1186/s40677-015-0021-x.
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 Dyn. 33 (2), 157–168.

9

A.R. Kacimov and J. Šimůnek

Journal of Hydrology 598 (2021) 126413
Tayfur, Gokmen, Swiatek, Dorota, Wita, Andrew, Singh, Vijay P., 2005. Case study:
Finite element method and artificial neural network models for flow through
Jeziorsko earthfill dam in Poland. J. Hydraul. Eng. ASCE 131 (6), 431–440.
The International Levee Handbook, 2013. CIRIA, London.
Vázquez, J.L., 2007. The Porous Medium Equation: Mathematical Theory. Oxford
University Press.
Wang, X.S., Neuman, S.P., Strack, O.D., Verruijt, A., Jamali, M., Seymour, B., Bear, J.,
Cheng, A.H.-D., 2011. Methods to derive the differential equation of the free surface
boundary. Groundwater 49 (2), 133–143.
Yan, Long, Xu, Weiya, Wang, Huanling, Wang, Rubin, Meng, Qingxiang, Yu, Jun,
Xie, Wei-Chau, 2019. Drainage controls on the Donglingxing landslide (China)
induced by rainfall and fluctuation in reservoir water levels. Landslides 16 (8),1583–1593.
Yang, Yongtao, Sun, Guanhua, Zheng, Hong, 2019. Modeling unconfined seepage flow in
soil-rock mixtures using the numerical manifold method. Eng. Anal. Boundary Elem.
108, 60–70.
Youngs, E.G., 1974. Seepage rates and the horizontal flow approximation. Water Resour.
Res. 10 (4), 874–876.
Zel’dovich, Y.B. and Kompaneets. A.S., (1950) On the theory of propagation of heat withthe conductivity depending upon the temperature. Ed.: Lukirsky P., Vavilov S.,
Kurchatov I. et al. “Collection Dedicated to the 70th Birthday of A.F. Joffe”. Akad.
Nauk SSSR, 61-71 (in Russian); Зeльдoвич, Я.Б. and КoМпaнeeц, A.C. (1950). Ктeopии pacпpocтpaнeния тeплa пpи тeплoпpoвoднocти, зaвиcящeй oт тeМпepaтypы.
Cбopник, пocвящeнный 70 лeтию aкaдeМикa A.Ф. Иoффe. Peд.: Лyкиpcкий, П.,Baвилoв, C., Кypчaтoв, И. и дp. Mocквa, AH CCCP, 61-71.
Zhou, X.P., Wei, X., Liu, C., Cheng, H., 2020. Three-dimensional stability analysis of bankslopes with reservoir drawdown based on rigorous limit equilibrium method. Int. J.
Geomech., ASCE 20 (12), 04020229.

Pasetto, I., 2014. Groundwater flow through the test dike constructed with dredged
materials. MSc thesis. Universita Delgi Studi Di Padova.
Peter, P., 1982. Canal and River Levees. Elsevier, Amsterdam.
Polubarinova-Kochina, P.Ya., 1962. Theory of Ground Water Movement. Princeton
University Press, Princeton. Second edition of the book in Russian is published in
1977, Nauka, Moscow.
Rehbinder, G., 1997. Relaxation of pore pressure in a slender core of a rockfill dam.
J. Hydraulic Res. 35 (2), 161–176.
Samarskii, A.A., Galaktionov, V.A., Kurdyumov, S.P., Mikhailov, A.P., 1995. Blow-up in
Quasilinear Parabolic Equations. Walter de Grutyer, Berlin.
Shanafield, Margaret, Bourke, Sarah A, Zimmer, Margaret A, Costigan, Katie H, 2021. An
overview of the hydrology of non-perennial rivers and streams. Wiley
Interdisciplinary Rev.: Water 8 (2). https://doi.org/10.1002/wat2.v8.210.1002/
wat2.1504.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Soil Mechanics Note No7. The Mechanics of Seepage Analysis, 1979. USDA.
Sophocleous, Marios, 2002. Interactions between groundwater and surface water: thestate of the science. Hydrogeology J. 10 (1), 52–67.
Stark, Timothy D., Jafari, Navid H., Zhindon, J. Sebastian Lopez, Baghdady, Ahmed,
2017. Unsaturated and transient seepage analysis of San Luis Dam. J. Geotech.
Geoenviron. Eng., ASCE 143 (2), 04016093. https://doi.org/10.1061/(ASCE)
GT.1943-5606.0001602.
Strack, O.D.L., 1989. Groundwater Mechanics. Prentice Hall, Englewood Cliffs.
Sun, Guanhua, Yang, Yongtao, Jiang, Wei, Zheng, Hong, 2017. Effects of an increase inreservoir drawdown rate on bank slope stability: a case study at the Three Gorges
Reservoir. China Eng. Geol. 221, 61–69.

10