Computational Geosciences (2020) 24:17–35
https://doi.org/10.1007/s10596-019-09879-8

ORIGINAL PAPER

Phreatic seepage flow through an earth dam with an impeding strip
A. R. Kacimov 1

2& N. D. Yakimov & J. Šimůnek

3

Received: 11 September 2017 / Accepted: 20 August 2019 / Published online: 26 December 2019
# Springer Nature Switzerland AG 2019

Abstract
New mathematical models are developed and corresponding boundary value problems are analytically and numerically solved for Darcian flows in earth (rock)–filled dams, which have a vertical impermeable barrier on thedownstream slope. For saturated flow, a 2-D potential model considers a free boundary problem to Laplace’s
equation with a traveling-wave phreatic line generated by a linear drawup of a water level in the dam reservoir.
The barrier re-directs seepage from purely horizontal (a seepage face outlet) to purely vertical (a no-flow boundary).
An alternative model is also used for a hydraulic approximation of a 3-D steady flow when the barrier is only apartial obstruction to seepage. The Poisson equation is solved with respect to Strack’s potential, which predicts the
position of the phreatic surface and hydraulic gradient in the dam body. Simulations with HYDRUS, a FEM-codefor solving Richards’ PDE, i.e., saturated-unsaturated flows without free boundaries, are carried out for both 2-D and
3-D regimes in rectangular and hexagonal domains. The Barenblatt and Kalashnikov closed-form analytical solutions
in non-capillarity soils are compared with the HYDRUS results. Analytical and numerical solutions match well whensoil capillarity is minor. The found distributions of the Darcian velocity, the pore pressure, and total hydraulic heads
in the vicinity of the barrier corroborate serious concerns about a high risk to the structural stability of the dam due
to seepage. The modeling results are related to a “forensic” review of the recent collapse of the spillway of theOroville Dam, CA, USA.
Keywords 2-,3-D transient seepage with a phreatic surface . earth dams . analytical models . Richards’ equation . fields of porepressure head . water content and velocity . heaving and suffusion2010 Mathematics Subject Classification• 76S05Flows in porous media; filtration; seepage• 30E25Boundary value problems• 35K61Nonlinear initial-boundary value problems for nonlinear parabolic equations• 65M60Finite elements, Rayleigh-Ritz and Galerkin methods, finite
methods
Electronic supplementary material The online version of this article(https://doi.org/10.1007/s10596-019-09879-8) contains supplementarymaterial, which is available to authorized users.
* A. R. Kacimov
anvar@squ.edu.om; akacimov@gmail.com; https://www.squ.edu.om/agr/Academics/Departments/Soil-Water-andAgricultural-Engineering

1

Department of Soils, Water and Agricultural Engineering, Sultan
Qaboos University, POBox 34, Al-Khod, 123 Muscat, Sultanate ofOman

2

Department of Thermal Engineering, Kazan State PowerEngineering University, Kazan, Russia

3

Department of Environmental Sciences, University of California atRiverside, Riverside, CA, USA

N. D. Yakimov
nyakimov@inbox.ru
J. Šimůnek
jsimunek@ucr.edu

18Comput Geosci (2020) 24:17–35

1 Introduction
Barriers and other impermeable or low-permeable entities,
naturally occurring or engineered inside porous massifs, significantly alter the incident seepage flows [1]. In geotechnicalengineering of earth dams, several types of barriers are designed for reducing seepage: clay cores, sheet pilings, cutoff-walls, among others. A motivating example for this paper
is the Oroville Dam (California) crisis of 2017. On the eve of
this crisis, Yakimov and Kacimov [2] accentuated the necessity to study a potentially dangerous phenomenon of seepage
across cutoff-walls in porous foundations of dams. In this typeof Darcian flow, invisible but insidious hydraulic gradients
across the barrier can cause its mechanical suffusion and deterioration of counter-seepage functions. In February 2017,
the main spillway of the Oroville Dam failed. Figure 1 showsa photo, where points F1 and F2 will be explained later, and dis half-width of the spillway chute.
The Oroville Dam proper is an engineered earthen dam that
forms a reservoir known as Lake Oroville behind it. The crisis
involved the need to release water from the reservoir during aperiod of rising water levels and significant precipitation (see[3, 4]). The initial releases were directed to the main concretelined spillway constructed adjacent to the dam proper along
the top of a bounding geologic formation. Part one of the crisis
involved the appearance of a crater in the concrete spillway
during the discharge, the subsequent uplifting of a part of theconcrete slab, and the exposure of the geologic materials underneath to significant erosion from high velocity flows. Part
two of the crisis involved a redirection of releases to a second,
unlined emergency spillway adjacent to the main spillway,
across the same body of poor geologic materials, which led
to additional significant erosion there. Although the cause of
the erosion in both cases was triggered by the flows occurring
across the top surface of this formation, its cause (in bothcases) was largely attributed to an improper understandingof how the (natural) geologic materials within the formation(including, presumably, resident saturation levels) would

d
F1

F2

Fig. 1 Photo of a collapsed spillway of the Oroville dam

behave under spillway use. It is key to note that the failed
geologic materials were not part of the engineered dam properand that their intrinsic integrity under conditions of main spillway use had been previously questioned.
The event spurred a surge of debates in the media, concernsof NGOs/public, comments by academics (e.g., Bea [5–7]),
and official “forensic” investigations (see e.g. [3]), which involved professional geotechnical engineers, professional societies, and federal/state authorities (i.e., US Army Corps ofEngineers (USACoE), The Association of State Dam Safety
Officials, U.S. Society on Dams, Canadian Dam Association,Federal Energy Regulatory Commission, California
Department of Water Resources, among others). For example,
Bea et al. [8], Cobos-Roa and Bea [9], and Steedman andSharp [10], among others, reviewed the failure of the New
Orleans levees and pinpointed the importance of seepage erosion. Bea [5, 6], investigating a similar failure in 2017, raised aquestion: did a sudden rise of the reservoir level of theOroville dam induce seepage in the part of the dam embankment close to it spillway and were these major factors thathave contributed to its collapse?
As after any catastrophe of this kind, the spotlight is on the
durability of hydraulic structures, especially “old” ones. The
inventory of USACoE includes almost 90,000 of dams in theUSA. NRC [11] identified 74,119 earth (rock)-filled dams in
the National Inventory of Dams (NID) [67] (USACoE [68])and 22,500 km of levees in the National Levee Database;according to NRC, about 56% of dams in NID will be older
than 50 years by 2020. Seepage-triggered failures of “new”
dams also generate a number of studies on the internal erosionof earth-rock fillings of embankments (see, for example, [12,13]).
Allegedly anticipated effects of global warming suggestthat structures will be subjected more often to extreme hydrologic events, which are characterized by intermittent droughts
and torrential rains of ever-increasing amplitudes (Vahedifard
et al., [14], Vano et al. [15]). For example, Lake Oroville(behind the Oroville dam) was in the center of public attention
for several years during the most recent drought in California.
The dried-out banks and the bed of the lake and the embankment itself were shown in the media as an illustration of “climate change,” just before the February-17 deluge.
Reservoirs like Lake Oroville are operated to achieve two,often conflicting goals: (1) storage of water for use in the drysummer months and (2) flood prevention through the maintenance of an open storage volume, which often requires controlled releases even when not needed for water supply. Thus,
operationally, water levels are typically kept lower in anticipation of an upcoming rainy season. The effects of climatic
swings that are now being experienced suggest that the strategy for this kind of dual-use operation must be reevaluated.
These swings from very dry to flooded conditions of earthworks exacerbate the natural aging of dams, many of which

19

Comput Geosci (2020) 24:17–35

were designed decades ago for hydrologically relatively lessvariable regimes of operations. Now, old earth (rock)-filleddams confront the following phenomena:
&

&
&

&
&

rapid drawdowns and drawups of the water level in thedam reservoir (lake) and tailwater accompanied by fluctuations in the water content of dams’ filling from dry, desiccated conditions to wet high pore-pressure ones,
infiltration into the dam’s body and riverbanks caused bytorrential rains,surge waves and currents released through the dams’ culverts and over the spillways, among others.

Seepage-induced/facilitated erosion of soil in earth(rock)-filled dams depends on dimensions of the embankment, its zonation (texture of the soil), andreservoir-tailwater levels. For example, piping, suffusion, and heaving (uplift) of the dam soil and parts ofthe constructed structure are controlled by seepage characteristics, mostly the hydraulic gradients and pore pressure [3, 16–18]. The criteria of stability against erosion
are commonly developed and tested in simple column
experiments, i.e., in well-controlled 1-D saturated seepage flow fields (see, for example, [12, 19–23]). Since
seepage in real dams is transient, involves free andmoving boundaries [24] of phreatic surfaces, and unsaturated zones above them, seepage flow (hydraulic gradients and pore pressures) needs to be characterizedusing mathematically complex 2-D or 3-D fields.
The geotechnical mitigation (forestalling) of seepage
erosion is achieved via the construction of toe, chimney,blanket drains, and drainage galleries/tunnels inside thedam body and in the downstream berm, relief wells (vertical and horizontal), cutoff walls, sheet piling, grout curtains, cushioning of key embankment components bygraded filter envelopes and geotextile wrapping, among
others [25–35]. Some dam components (such as spillways)
are added to the embankment long after the original design/construction, and therefore their impact on seepage inthe dam may not be properly incorporated into the local/
overall analysis of the structure’s stability, as Troyanowski[36] assessed and warned for dams’ spillways.
The overall post-crisis opinion is that the collapse of thespillway was triggered by high-velocity overland flows occurring overexposed and un-engineered geologic surfaces. In thispaper, we address a different factor: seepage near an impermeable barrier on the downstream slope of an earth dam that
concurs with what France et al. [3] stated: “A detailed study ofthe seepage flow under the slab would be necessary in thiscase to determine the uplift pressure distribution.”
The configuration of the Oroville Dam crisis is somewhat
different than the idealized components of the problem addressed in this paper; in particular, we assumed the following:

&
&

&
&

The dam body is a rectangular parallelepiped that makes
possible conjugation of analytical and finite-element
modeling.
We ignore the herringbone drains under the slabs of thespillway [3].
The spillway is a vertical rather than a tilted barrier.
The barrier is impermeable to both seepage from the soil
to the atmosphere and for leaks from the surface water into
the soil subjacent to the spillway slabs; in other words, we
do not study a coupled surface-subsurface flow which
may have contributed to the crisis [3].
Soil is homogeneous and isotropic.
The lake stage in our models is also simplified as compared with the prototype dam.

France et al. [3] pointed out that no detailed analysis ofseepage (in the report, called “groundwater flow”) has been
done by their forensic team. The main questions relevant toseepage are:a) Can transient seepage with a phreatic surface, which propagates in a “pre-dried” embankment body, be a plausible
cause of the collapse of an earthwork?
b) Can new analytical and numerical models of “globalscale” seepage through earth dams add better insights towhat the classics [37, 38] posed as potential threats to
dam stability?
c) How to connect these large (dam) scale seepage models towhat is tested in column experiments with “core samples”
of the embankment porous materials, subject to high andvariable gradients (see, for example, [39])?
d) Can relatively fast external hydrological/hydraulicdrivers, such as a rapid rise of the reservoir level or intensive rain, induce singular (ruinous) consequences via relatively slow seepage flows within the components of aporous dam?

2 Conceptualization of the flow problem
In our analysis, we consider Darcian, isothermic, single-phase,and incompressible water flow in a rigid porous skeleton (soilor rock-filling); saturated hydraulic conductivity k [m/s] andthe van Genuchten–Mualem (VGM) soil hydraulic characteristics [40] are constant.
A physical 3-D domain is symmetric with respect to thebarrier axis (Fig. 2a) and we study the right half of it. A system
of Cartesian coordinates x1x2y is centered at point A1. An
impermeable horizontal rectangle A1A2D3D1 bounds the flow
domain from below. The barrier is assumed to be an impermeable vertical rectangle D1F1F2D2 placed on a vertical soil

20Comput Geosci (2020) 24:17–35

Fig. 2 A schematic of 3-D phreatic flow through an earth dam(a), a vertical cross-section at
x1=a (b)

y

a
b

E2

H(t)

c

F1
B1

A2

x2

A1
f

kL

B2
B1

TailwaterJ
C1 J2 3

D2

wall

H0

C0 C2

C3

x1

D1
DL

Seepage face
B3

B1

a

3-D phreatic
surface

B3

C3
C2 D3
C0

AL

g

F2

Spillway

B3

E1
Lake

F3
d

g
C1

H(t)

Constant
head
line

H(t)

Wall
H0

b

slope D1F1F3D3 that is an idealization of a spillway chute in
Fig. 1.
We introduce the total hydraulic head h(t, x1, x2, y)=p(t, x1,
x2, y)+y, where p is the pressure head. The reservoir waterlevel H(t) is rising such that the vertical face A1B1B2B3A2(the dam reservoir boundary) is under a time-variable condition in terms of h. The tailwater level H0 is constant and,correspondingly, the total head h over the rectangle
D2C2C3D3 is constant both in time and space. Through thehorizontal rectangle E1E2F3F1 (y=c in Fig. 2a) which represents the ground surface (a dam crest), rain recharges the watertable.
The 3-D phreatic surface B1B2B3J3J2C1 (sketched in Fig.
2a) is affected by the barrier and outcrops at the tailwatervertical slope with a seepage face J3C3C2J2. A vertical crosssection x2=b in Fig. 2a depicts a phreatic line B3J3, with outflow through a seepage face segment C3J3 and a constant headsegment D3C3. Another phreatic line, B1J1 is shown in a vertical cross-section x2=0 (Fig. 2a) In this cross-section, thesegment C1D1 is an impermeable wall and line of symmetry,
i.e., the horizontal components of the Darcian velocity vector
vanish there.
In 3-D, the rectangle D2J2J3D3 in Fig. 2a is a vent foroutseeping water. The adjacent saturated soil is routinely
viewed as a zone of the potential risk of seepage-induced
erosion. However, despite Troyanowski’s [36] warning, soil
adjacent to the impermeable surface D1F1F2D2 (i.e., theunder-chute soil) is seldom scrutinized for stability.
The trace of the 3-D phreatic surface B1B2B3J3J2C1 in thevertical plane x1=a is shown in Fig. 2b as curve C1J2J3. The
phreatic surface rises due to two drivers: the rise of the reservoir level H(t) and accretion from the vadose zone. This creates the following threats to dam stability:

&

&

&

A1

D1

A2

J3
C3

H0

D3

A high pore pressure on the soil side of the barrier
D1F1F2D2 creates a force perpendicular to this surface thatcan cause uplift-heaving of the concrete blocks and theirjacking. Even small seepage-induced displacements of theblocks can be exacerbated by supercritical surface flow
(shown in Fig. 1a) on the “external” (with respect to theporous medium) surface of the barrier.
High vertical upward-oriented hydraulic gradients in thevicinity of point C1 in Fig. 2b can cause dislodging andtranslocation of soil particles, especially, of the fine soil
fraction under the spillway blocks; that reduces the blocksoil grip.
High lateral (in the x2-direction) seepage gradients due to aphreatic mound due to the high variation of the elevationof the points on the curve C1J2J3 in Fig. 2b can cause
suffusion of soil under the barrier, as well as piping in
the open (seepage face) zone of the tailwater. Clearly,
without a barrier (concrete spillway) C1J2J3 in Fig. 2b
has a much milder slope and less threat of lateral
suffusion.

Various analytical and numerical models have been used
for studies of seepage in dam bodies (e.g., [41–53]). Below wefollow the recommendation of Alvi [54] and use various
mathematical models for determining 2-D and 3-D functions
h and p, as well as hydraulic gradients within the dam.
In Section 3, we use the Dupuit-Forchheimer (DF) approximation and present a new analytical solution to the problem
of a steady saturated flow occluded by the barrier (Fig. 2a). In
Section 4, we revisit and further develop the Kacimov andYakimov [55] analytical solution to a 2-D transient problem
of a potential saturated flow with a phreatic surface rising in acapillarity-free soil. A vertical barrier occludes this rise. We

21Comput Geosci (2020) 24:17–35

solve a moving boundary value problem (BVP) and calculatebuild-up of the seepage-induced pore pressure and hydraulicgradients. These analytical results are compared with
HYDRUS [56] simulations for a saturated-unsaturated flow
in soil with capillarity. There are no moving or free boundaries
in HYDRUS. In Section 4, the HYDRUS model is used tosimulate flow in a rectangular domain representing the dam’s
vertical cross-sections close to the plane x2=0 (Fig. 2a). The
phreatic line in HYDRUS is a zero-pressure isobar, whichmoves with the family of all other isobars.
In Section 5, we consider the early stage of the deluge
when the earthwork is dry, and the drawup in the reservoir
generates a saturated “tongue,” which propagates inside ofthe dam. We consider a monotonic and non-monotonic
reservoir hydrograph. The former is modeled using the
Barenblatt analytical solution ([57]; referred to as PK77) for capillarity-free soil on an impermeable foundation.
For the latter, after a rapid rise of the reservoir level up toa peak stage, a relatively slow drawdown takes place (thisis the actual situation in the Oroville lake when after
February 2017 the reservoir level receded). We recur to
the DF approximation and consider a transient nonlinear
Boussinesq equation. Unlike in Section 3, we assume thatthe dam bed is a “leaky boundary”, i.e., that it is made ofa thin claylike stratum (aquitard shown as ALA1D1DL in
Fig. 2a), which transmits water down from the dam body.
The Kalashnikov [58] analytical solution predicts that thesaturated zone in a vertical cross-section is a right triangle, whose hypotenuse (a moving phreatic surface) tilts
with an initially increasing and then decreasing slope.
The HYDRUS simulations are compared with the
Kalashnikov solution at a time when the dynamics of porewater is controlled by the rise of the water level in the reservoir
rather than the capillarity of the dam body and substratum.
HYDRUS results are also compared with Barentblatt’s solution, which is an exact one for both the Boussinesq equationand the 2-D potential model with a nonlinear phreatic surface
boundary condition. HYDRUS reconstructs a perfectly prismatic saturated “tongue,” which propagates from the reservoir
into the dam at a constant rate of drawup.
In Section 6, we use HYDRUS to simulate a 3-D transient saturated-unsaturated flow in the entire box shownin Fig. 2a. These results (in the limit of constant H in Fig.
2a) are compared with the analytical solution in Section 6.

3 Analytical solution for vertically averaged
steady flow impeded by a strip-barrier
In this Section, we neglect capillarity and flow transiency. We
assume that the dam reservoir level is constant, i.e.,
H(t)=H1=const. We also adopt the DF approximation (seePK-77, [1]) that eliminates the vertical coordinate y in

Fig. 2a, i.e., the total head h(x1,x2) in the saturated zone coincides with the height hs(x1,x2) of a steady phreatic surface
above the horizontal impermeable base A1A2D3D1. The impact of rain and the vadose zone on the saturated DF flow is
accounted for through the recharge rate N=const>0 [m/s],
which is space- and time-independent over the rectangle
E1E2F3F1 (Fig. 3a). The analytical model below describes
how DF flow “winds” (horizontally) around the impermeablebarrier D1F1F2D2. Fig. 3a represents a planar view.
The water table position hs(x1, x2) and the discharge vector!
F S ðx1 ; x2 Þ obey the following BVP in terms of Strack’s potential (see, for example, [1, 59]):
!
F S ¼ ∇Φ;
∇2 Φðx1 ; x2 Þ ¼ −N

H 1 2 −H 0 2¼ const for 0 ≤ x2 ≤ b; i:e:along A1 A2 ;2Φða; x2 Þ ¼ 0 for d ≤ x2 ≤ b; i:e:along D2 D3 ;∂Φðx1 ; bÞ ¼ 0 for 0 ≤ x1 ≤ a i:e:along A2 D3 ;∂x2∂Φðx1 ; 0Þ ¼ 0 for 0≤ x1 ≤ a i:e:alongA1 D1 ;∂x2∂Φðx2 ; aÞ ¼ 0 for 0 ≤ x2 ≤ d i:e:along D1 D2∂x1
Φð0; x2 Þ ¼ Φ1 ¼ k

ð1Þh2 −H 2where Φðx1 ; x2 Þ ¼ k s 2 0 [m3/s] is the discharge potential,∇ and ∇2are the gradient and Laplacian 2-D (i.e., in the plane!
x1,x2) operators. The vector F S has components Fx1 and Fx2along the corresponding axes (see Fig. 3a).
We decouple the solution to the BVP (1) in the following
manner:Φ ¼ ΦP þ ΦL

ð2Þ

where the first term satisfies the BVP for the Poisson equation:∇2 ΦP ðx1 ; x2 Þ ¼ −N
ΦP ð0; x2 Þ ¼ Φ1 for 0 ≤ x2 ≤ b; i:e:along A1 A2 ;ΦP ð0; x2 Þ ¼ 0 for 0≤ x2 ≤ b; i:e:along D1 D3 ;∂ΦP∂ΦPð x1 ; 0Þ ¼
ðx1 ; bÞ ¼ 0 for 0 ≤ x1 ≤ a; i:e:along A2 D3 andA1 D1
∂x2
∂x2

ð3Þ
The solution to the BVP (3) is trivial in 1-D:


x21
Na Φ1
−
ΦP ðx1 Þ ¼ −N þ
x1 þ Φ 1
2 a
2

ð4Þ

From eq. (4), it follows:
∂ΦP
Na Φ1
−
¼ −N x1 þ
2 a
∂x1

ð5Þ

The second term on the right-hand side (RHS) of eq. (2) is aharmonic function which satisfies the following BVP:

22

Comput Geosci (2020) 24:17–35

Fig. 3 The horizontal base 3-Dflow domain in Fig. 2a and therectangle in which Strack’s potential is found (a), the domain of
a fictitious flow and the boundary
value problem for the harmoniccomponent of Strack’s potential(b), the hodograph domain of afictitious flow (c), the complex
potential domain of a fictitious
flow (d), and an auxiliary halfplane (e)

A2

x2

x2

D3

Real
b

D2 stagnaon
point

L=0
Separatrix

d
D1 x1

A1

x2

Ficous A2
stagnaon
point

a

0,

L

D3

G1

GVf

L=0

L

x2

0,

L

D1

D2

u0

L

x1

D3

iQl

D3
A2

(W L=

L +i

G1 A2

A1

D1 uf
u0

(

i )

x1

c

0

L

iQ D2

G2

D2

Gz

A1

b

(Vf =u-i v)(z=x1+i x2)

G2

a

D2

-vf

Ql

L

G

L)

e

GWL

A1

D1

D2 D3 G1 A2

G1
L

d

A1

D1

∇2 ΦL ðx1 ; x2 Þ ¼ 0
ΦL ð0; x2 Þ ¼ 0 for 0≤ x2 ≤ b; i:e:along A1 A2
ΦL ð0; x2 Þ ¼ 0 for d≤ x2 ≤ b; i:e:along D2 D3 ;∂ΦL∂ΦLð x1 ; 0 Þ ¼
ðx1 ; bÞ ¼ 0 for 0 ≤ x1 ≤ a; i:e:along A2 D3 andA1 D1
∂x2
∂x2
∂ΦL
∂ΦP
Na Φ1
ða; x2 Þ ¼ −ða; x2 Þ ¼
þ
¼ u0 for d≤ x2 ≤ b; i:e:along D1 D2
2
∂x1
∂x1
a

ð6Þ
L
The last boundary condition for ∂Φ∂x1 ða; x2 Þ in eq. (6) fol-lows from eq. (5). For ΦL, eq. (6) thus states a mixed BVP(see the details of this type of BVPs in PK-77). We use thehodograph method to solve this BVP. Namely, we consider
ΦLas a potential of a fictitious pseudo-flow with a velocity!
V f ¼ ∇ΦL , a vector having a horizontal component uf anda vertical component vf. Then the complex function Vf=ufivf=dWL/dz is holomorphic in the complex plane z=x1+ix2(the corresponding domain Gz is shown in Fig. 3b). The
function WL = ΦL + iΨLis a complex potential of the pseudoflow. ΨL is a stream function of the pseudo-flow. It is harmonic and conjugated with ΦL via the Cauchy-Riemann conditions. The mirrored hodograph domain is a right angle GVf,
which is shown in Fig. 3c, u0 is a given (see eq. 6) positive
constant. We note that the point G1 in GVf is a stagnation
point of the fictitious flow; the real flow, which is quantified
by the full Strack potential Φ, has D1 as a stagnation point(compare Fig. 3a, b). The separatrix G1G2 divides thepseudo-flow into a part, which feeds D1D2 by the segment
D2D3, from one emanated from A1A2. The complex potential
domain, GWL, corresponding to Gz, is depicted in Fig. 3d.
GWL has a horizontal cut A2G1D3 along which ΨL = Ql thatis the pseudo-flow rate from A1A2 into D1G2. The total
pseudo-flow rate into D1D2 is Q. The pseudo-flow rate from
D2D3 into D2G2 is Qr=Q-Ql. It is noteworthy that the scalars

Q, Ql, and Qr are to be found. Both the real flow in Fig. 3a
and fictitious flow in Fig. 3b have a singularity at point D2
!
where both the Strack discharge vector F S ða; d Þ and the!
fictitious velocity V f ða; d Þ blow up.
We introduce an auxiliary variable ζ = ξ + iηand map
conformally the half-plane Gξ: ξ>0 (Fig. 3e) onto Gz by the
Shwartz-Christoffel function:dτζ
z ¼ ic1 ∫−1=λ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
1−λ2 τ 2 ð1−τ 2 Þ

ð7Þ

where the branch of the root is fixed in such a manner thatthe root is positive at ζ = ξ, η = 0, − 1 < ξ < 1. The mapping
constant is determined from the given width of Gz , z(-1)=a,c1=a/K’, where K’ is the complete elliptic integral of the firstpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
kind with the modulus 1−λ2 . We know the height of Gz in
Fig. 3b. Therefore, the relation z(1)=a+ib in eq. (7) gives:
b 2K
¼ 0
a
K

ð8Þ

where K is the complete elliptic integral of the first kind withthe modulus λ. From eq. (8) we find the value of λ by the
FindRoot and EllipticK routines of Wolfram’s [60]Mathematica.
From eq. (7) the affix δ of point D2 in Gζ (Fig. 3e) is:
d¼

a δ
dτ
∫ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
K 0 −1 1−λ2 τ 2 ð1−τ 2 Þð9ÞEquation (9) is a nonlinear equation with respect to δ,
which we again solve by the FindRoot routine of Wolfram’sMathematica.
Next, we map conformally Gζ onto the quadrant GVf (Fig.
3b) by the Shwartz-Christoffel function:

23

Comput Geosci (2020) 24:17–35

dW Lζ¼ c2 ∫−1 ðτ−δÞ−3=2 ðτ þ 1Þ−1=2 dτ þ u0dz
sﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2c2
1þζ
þ u0
¼−
ζ−δ
ð1 þ δ Þvﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
0
1
u
u
ζ
u

1
dτ
B
C
2
2
2
hsw ðξÞ ¼ u
t H 1 −H 0 @1− I Re ∫ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
A þ H 0;
2 2
−1=λ
1
ðτ−δÞð1−τ Þ 1−λ τ

Vf ¼

ð10Þ

x2 ð ξ Þ ¼

K þ F ðarcsinξ; λÞ;
K0−1 < ξ < δ

ð15Þ
where the branch of the root is selected positive along theξ−axis in Fig. 3e, in the interval ξ > δ (and therefore at ξ < −1). The mapping constant c2 will be found later. Obviously, in
the interval −1 < ξ < δ, the root in eq. (10) is imaginary (see theray D1D2 in Fig. 3c). From eq. (10), the singularity at point D2
is Vf ∼ (ζ − δ)−1/2that can be expressed via ζ(z) from eq. (7).
Obviously, z(ζ)is regular at point D2.
We assume that point A1 is fiducial, i.e., at this point WL = 0(the complex potential is defined up to an arbitrary complex
constant, see PK-77). ThenζW L ¼ ∫−1=λdW Ldzðτ Þ ðτ Þdτ
dτ
dz

ð11Þ

By using eqs. (10) and (7) in the integrand of eq. (11) we
get
W L ¼ u0 zðζ Þ−i

ζ
2ac2
dτ
∫ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ð12Þ
K 0 ð1 þ δÞ −1=λ ðτ−δÞð1−τ Þ1−λ2 τ 2 

Now we determine c2 from the condition at the point ξ = δ
where z=a+id and WL=iQ.
From the last line in eq. (6) we have a constant horizontalL
pseudo-velocity∂Φ∂x1 ¼ u0 along D1 D2 and therefore Q=u0*d.
Then, from eq. (12) we get
u0 a ¼δ2ac2dτ
∫ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
0
K ð1 þ δÞ −1=λ ðδ−τ Þð1−τ Þ1−λ2 τ 2 

ð13Þ

From eq. (13) we immediately express c2 and put it backinto eq. (12) which becomes
W L ¼ u0 zðζ Þ−i

au0 ζ
dτ
∫ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ; I 1
I 1 −1=λ ðτ−δÞð1−τ Þ1−λ2 τ 2 δdτ¼ ∫ qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ−1=λðδ−τ Þð1−τ Þ 1−λ2 τ 2

ð14Þ

We separate the real part in eq. (14) and get ΦL. Then, we
re-assemble eqs. (11) and (4) to get the final solution in termsof Strack’s potential (2).
We introduce dimensionless variables by normalizing all
geometrical sizes to a and the Darcian velocity to k.
The parametric equation of the phreatic curve C1J2 (Fig.
2b) is:where F is the incomplete elliptic integral of the first kind.
Similarly, from (7) and (14) we get parametric equations of
B1C1 as hss(ξ), x1(ξ), − 1/λ < ξ < − 1.
Suffusion, piping, and lessivage are controlled by the vec! !
tor of the hydraulic gradient (PK-77) J ¼ V =k. In the DF!
model, while J does not vary with y, it varies with x1 and x2 asfollows:!
J ¼ −∇hs ðx1 ; x2 Þð16Þ!
In eq. (16), J has components Jx1 and Jx2 in the x1 and x2
directions, respectively. For example, from eq. (15) along thebarrier Jx1 = 0 and the component of the magnitude of Jx2 is: d ½h ðξÞ  d ½x ðξÞ −1  sw
2J x2w ðξÞ ¼ 
 ; x2 ð ξ Þ

 dξ
dξ
¼

K þ F ðarcsinξ; λÞ;
K0−1 < ξ < δ

ð17Þ

In computations involving eqs. (15) and (17), we used theMathematica routines NIntegrate, EllipticF, andParametricPlot. Figure 4a, b, c shows hsw(x2), Jx2w(x2), andhss(x1) for dimensionless parameters b=4, d=3/15, H1=102/
50, H1=2/50, and N=0 (see also Figs. 2a and 3a for thecorresponding points). The choice of these quantities is connected to ones in numerical simulations of Section 6 below.
Figure 4a and b show a steep drop of the phreatic surface frompoint C1 to J2 and extremely high gradients, especially, close
to the rim of the barrier, D2J2. These gradients well exceed the!
safe limit of PK-77, which is j J j ¼ 1.
For assessments of heaving (uplifting) of the barrier blocks,
we introduce a dimensionless (normalized to ρga, where ρ isthe density of water) integral pore pressure acting on the barrier:d hs ð1;x2 ÞPi ¼ ∫0δ∫

0½hsw ðξÞ2 dx2 ðξÞ
dξ
dξ
2
−1

¼ ∫

½hs ð1; x2 Þ2
dx2
2
0

d

ðhs ð1; x2 Þ−yÞdydx2 ¼ ∫ð18Þ

i.e., an integral of the pore pressure over a curvilinear trapezium D1C1J2D2 (Fig. 2b). This trapezium is the “wettedarea” of the barrier, “submerged” under the phreatic surface.
The force (normalized to the specific weight of pore water and24Comput Geosci (2020) 24:17–35

Fig. 4 a The phreatic curve C1J2
for dimensionless parameters
b=4, d=3/15, H1=102/50, H1=2/
50, N=0. b The hydraulic gradientalong the barrier for
dimensionless parameters b=4,
d=3/15, H1=102/50, H1=2/50,
N=0. c The phreatic curve B1C1
for dimensionless parameters
b=4, d=3/15, H1=102/50, H1=2/
50, N=0a) due to the hydrostatic pressure isFi = Pi/Ar, where Ar is the“sub-phreatic” (wetted) area of D1C1J2D2. For the example
illustrated in Fig. 4, computations give Pi=0.14 and Fi=0.49.
Similarly, we can find the center of the seepage pressure exerted
on D1C1J2D2. If the barrier consists of blocks, then the seepage
force acting on a block, the point of its application, and the force
moments, which are needed for assessing the threat of jacking,can be easily found by integrations similar to eq. (18).

4 Transient 2-D flow in a vertical plane,
for a constant-rate drawup
In this section, we consider a transient element of 3-D flow in
Fig. 2a, i.e., seepage in the vertical plane x2=0. Figure 5 a shows
the vertical cross-section x1A1y, in which we introduce Cartesian
coordinates XY with the origin at a moving point C1 where therising phreatic surface intersects the vertical barrier. Although
the real 3-D phreatic surface drops from C1 to the rim of thebarrier (point J2 in Fig. 2a), in this Section we assume that thebarrier half-width d is high enough so that we can assume flow
to be 2-D in the plane XY, i.e., close to the line D1C1 in Fig. 2a.

4.1 Analytical similarity solution for phreatic flow
impeded by a vertical barrier
The water level in the dam reservoir rises with time t [days] ata constant rate R [m/day]. We assume that this rise continues
for a sufficiently long time such that the phreatic line B1C1
propagates in a porous slab of a width a as a traveling wave,
i.e., without any change in its shape (see details of the corresponding similarity solution in [55]). This regime is realized if
the impermeable bed A1A2D3D1 is far beneath B1C1, i.e., H(t)in Fig. 2a is high enough. Correspondingly, the base A1D1 is

not even shown in Fig. 5a. Then, we can ignore A1A2D3D1
and study flow in a half-strip AB1C1, where A is an “infinity
point.” A rising curve B1C1 caps the domain Gz in Fig. 5a. The
difference in elevations of points B1 and C1, the depth, tr, of
this advancing “trough” on the phreatic line is a part of thesolution.
Similarly to Section 2, the stability of C1A against seepageinduced erosion is controlled by the distribution of the porepressure head p(X,Y,t) and the Darcian velocity!
V ¼ −k∇hðX ; Y ; t Þ, especially, in the vicinity of Y=0, i.e.,
close to the point C1 in Fig. 5a. In the following text, we workin a system of coordinates, which translates upwards with aspeed R.
In the analytical solution below, we neglect accretion on thephreatic surface and capillarity of the soil, i.e., assume a piston
type rise (or imbibition, see [61]) of B1C1. This means that thesoil beneath this line is saturated and above it is dry, with an
instantaneous saturation occurring when B1C1 passes throughthis point. The porosity of the soil is 휃s and the irreducible
moisture content is 휃r (in the analytic solutions below weassume 휃r=0).
Due to symmetry, we can mirror the physical strip GZ about
the C1Y axis that results in an “image reservoir” on the right ofthe barrier (Fig. 5a). Point Bm is a mirror image of point B1. We
introduce a complex physical coordinate Z=X+iY. The complex
coordinates of points B1 and Bm are ∓a + i tr, correspondingly.
Similarly to Kacimov and Yakimov [55], we introduce acomplex potential w = ϕ + iψ,ϕ = − kh and a function Vc=dw/
dZ, which is complex-conjugated with the complexified
Darcian velocity V=u+iv, i.e., Vc=u-iv where u and v are,correspondingly, the X and Y components of the Darcian ve!
locity vector V . The domain GVc is a circular digon (a luneshaped polygon) shown in Fig. 5b. The curve BmC1B1 is an25

Comput Geosci (2020) 24:17–35
Fig. 5 A rising phreatic surface inthe vertical cross-section x1A1y(Fig. 2a) with a baffling verticalbarrier (a), a hodograph domain(b), a domain in the complex potential plane (c), and a reference
half-plane (d)R*(z=X+i Y)--v

Y(V=u-i v)
Bm

B1

tr
Real
C1
reservoir
with aGZ
constant
FictitiousarateporousRealdrawupslab
porous
slab

B1
u
uB1

A
GVc

X
Fictitious
reservoir

R<1

C1 vC1R>1bC1Real wall, line at which
3-D seepage in Fig.2abifurcates

A

a

(w=

B1
A

Bm

(

)
A

Gw

i )

G
C1 B1

Bm

A

C1

d
Bm

arc of a circle of a radius (k + θsR)/2, centered at a point i ⋅ (k− θsR)/2and at an angle πγ with the u-axis. Darcian velocities
at points B1 and C1 are horizontal and vertical with magnitudes uB1 and vC1, respectively. We select point A as fiducial,i.e., wA=0. Then, the complex potential domain Gw, shown in
Fig. 5c, is a digon with an angle πv between the phreatic
surface and the coordinate axis.
We introduce dimensionless variables (z∗, X∗, Y∗, h∗, tr∗)= ( z , X , Y, h , t r ) / a , V*; u*B1 ; v*C1 ¼ ðV; uB1 ; vC1 Þ=k;
w*; ϕ* ; ψ* ¼ ðw; ϕ; ψÞ=ðk*aÞ, and R ∗ = qsR/k. For the sake of brevity, we drop the superscript “*”. Kacimov andYakimov [55] obtained the shape of B1C1 and found that
pﬃﬃﬃ
pﬃﬃﬃ
uB1 ¼ R; vC1 = R; cotðπν Þ ¼ R; γ = 1 − 2 ⋅ ν (see also
the phreatic surfaces in Fig. 5a of [55]). We note that if R>1
then the center of the circle B1C1Bm in the Vc-plane is below
the uc -axis and γ>1/2 (Fig. 5b, dashed line).
Kacimov and Yakimov [55] introduced an auxiliary variable z = x + i × h and considered a reference half-plane Gζ,− ∞ < ξ < ∞ , η > 0, ζ = ξ + iη (Fig.4.1d) onto which Gz has
been mapped by the function, which is obtained by integrationof the following function:
"
#
dZ
i
eνπi
e−νπi
¼
Ξ ðζ Þ ¼−
dζ
π ð1 þ ζ Þν ⋅ð1−ζ Þ1−ν ð1 þ ζ Þ1−ν ⋅ð1−ζ Þνð19ÞNow, we use the PK-77 method for determination of thefunction Ωðζ Þ ¼ dw
dζ , where w(ζ) is the complex potential as afunction of ζ:#
pﬃﬃﬃ "
i⋅ R
enpi
e−npi
Ω ðζ Þ ¼
þ
π
ð1 þ ζ Þν ⋅ð1−ζ Þ1−ν ð1 þ ζ Þ1−ν ⋅ð1−ζ Þν

ð20Þ

c
Ω
Consequently, for V c ¼ dw
dZ ¼ Z from eqs. (19)-(20), we
obtain:V c ðζ Þ ¼

pﬃﬃﬃ ð1 þ ζ Þγ ⋅eνπi þ ð1−ζ Þγ ⋅e−νπiRð1 þ ζ Þγ ⋅eνπi −ð1−ζ Þγ ⋅e−νπið21ÞThen, integration of eqs. (19)-(20) yields:
1 1
1
tr ¼ Im ∫0 Ξ ðτ Þ dτ ¼ ∫0 Ξ ðτ Þ dτ þ 1
i
1 i⋅η
i⋅η
Y ðηÞ ¼ Im ∫0 Ξ ðτ Þ dτ ¼ ∫0 Ξ ðτ Þ dτ; 0 ≤ η < ∞; 0 ≥ Y > −∞ii⋅ηhðηÞ ¼ −tr þ ∫0 Ωðτ Þ dτ; 0≤ η < ∞; −tr < h < 0vðηÞ ¼ i⋅V c ði⋅ηÞ 0 ≤ η < ∞; R < vy < 0ð22ÞFrom eq. (22), the distribution of the pore pressure along
AC1 is obtained as a parametric equation:
Y ¼ Y ðηÞ; pðηÞ ¼ hðηÞ−Y ðηÞ þ tr;0≤η < ∞

ð23ÞSimilarly, the distribution of the velocity along AC1 is determined from eq. (22) by the parametric equation:
Y ¼ Y ðηÞ; v ¼ vðηÞ;0≤η < ∞

ð24ÞFigure 6 presents computations based on eqs. (22)-(24),which we carried out in Wolfram’s Mathematica for three
dimensionless drawup rates, R=4, 1, and 0.25. The corresponding curves 1v-3v are the Darcian vertically-oriented velocities v(-Y) along C1A in Fig. 5a; the curves 1h-3h show thevertical distribution of the total hydraulic head h(-Y) along thebarrier and the curve 1p is the distribution of the pore pressure
for R=4. These results illustrate that at dimensionless depths of

26

Fig. 6 Darcian velocities v (curves 1v-3v are for dimensionless drawup
rates R=4, 1, 0.25, respectively), the total hydraulic head h (curves 1h-3h
are for R=4, 1, 0.25, respectively) and the pore pressure p (for R=4) along
AC1 in Fig. 5aY ∼ − 3, water close to the barrier is almost stagnant, the totalhead is almost the same as in the reservoir, and the pore pressure is almost hydrostatic. However, closer to the rising phreatic surface (point C1 in Fig. 5a), the hydraulic gradients are
high, i.e., in this zone the soil adjacent to the barrier is temporary, when the traveling phreatic line B1C1 passes through,susceptible to erosion.

4.2 HYDRUS solution in a rectangle with an impeding
vertical barrier
For comparisons, in this Subsection, we return to dimensional
quantities and use HYDRUS [56] to simulate a transient 2-D
saturated-unsaturated flow in a rectangular domain with a vertical barrier, similar to C1A in Fig. 5a. We consider drawup for
0<t<7 days with a rate of R =50 cm/day. Below we use the
HYDRUS notations and coordinates (XH, ZH) with the originof coordinates at point A1 (Fig. 2a). The rectangle A1E1F1D1of sizes a=100 cm and c=450 cm in XH and ZH directions,
respectively, is shown in Fig. 7.
First, we selected loam as the dam material, for which weretrieved the physical properties from the HYDRUS Soil
Catalogue, i.e., k=25 cm/day, 휃s=0.43, 휃r=0.078, and theVedernikov and Bouwer [40] constants 훼=0.036 cm-1 andComput Geosci (2020) 24:17–35

n=1.56. In the analytical solution, the corresponding dimensionless drawup rate is 0.71.
The transport domain was discretized into finite elements with the following parameters and/or statistics: atargeted size of finite elements of 2.5 cm, 11,771 nodes,
440 1-D and 23,100 2-D finite elements. We assumed
that prior to the reservoir filling and a seepage event (att=0) the whole rectangle was dry with an initial pressure
head of pd=-10,000 cm. The boundary conditions wereselected as follows. No flow was assumed along thethree sides of the rectangle, A 1D1, D1F1, and F1E1.
Along the “expanding” boundary segment A1B1, waterlevel in the reservoir rises and a time-variable boundary
condition with a linear increase of the pressure headp(0, Z, t) with time for Z varying from point A 1(fixed) to point B1 (moving) was specified. In otherwords, the total head h along A1B1 was constant at
any fixed time instance and was increasing linearly as
a function of time. The “shrinking” vertical segment
B1E1 was assumed to be impermeable because the pressure head there was negative. This boundary condition
is consistent with one in the analytical model ofVedernikov and Bouwer [63].
Figures 7 a, b, c, d, and e show colored plots of p(X, Z, t) fortime instances t=1, 2, 3, 4, and 7 days. The similarity regime,
examined in the analytical solution (Fig. 5), is indeed
established after about 3 days. From Fig. 7a and a boundaryline chart p(X, 0, 1) generated by HYDRUS we obtain that at
t=1 day the tip, WB, of the phreatic surface is at the locus
XWn=76 cm. The isotachs V=const are shown in Figs. 7f–g
at time t=7 days. These isolines demarcate a “tongue” of highhydraulic gradients (their magnitude exceeds 2 near point B1)
close to the reservoir boundary. It is geotechnically worrisome
that, for example, at t=7 days (Fig. 7f), a zone of gradients ofmagnitude higher than 0.1 extends more than 100 cm beneaththe phreatic surface. Fig. 7g shows the contour lines of horizontal components, u, of the Darcian velocity vector plotted at
t=7 days. The magnitude of velocity at point C1 is VC1≈16 cm/
day. The analytical solution from Subsection 4.1 gives thecorresponding value of VC1 =21.5 cm/day when 휃s=0.43 isused in the definition of the linear average velocity and VC1=17.6 cm/day when the porosity is replaced by an “effective”
value, 휃e =s - 휃r =0.43-0.078, i.e., VC1=e*R. The depth of thetrough (Fig. 5a) obtained by HYDRUS is tr≈39 cm as compared with those obtained by the analytical solution, viz. tr =Fig. 7 The results of the HYDRUS computations for a constant rate of
drawup in a rectangle made of loam (k=25 cm/day, 휃s=0.43, 휃r=0.078,-1
훼=0.036 cm , and n=1.56; R=50 cm/day, the initial pressure head pd=10000 cm) (a-g) and sand (k=0.495 cm/min, 휃s=0.43, 휃r=0.045,-1
훼=0.145 cm and n=2.68; R=25 cm/min, the initial pressure head pd=100 cm) (h-i). The pressure heads are displayed in (a-e) and (h); water
fluxes in (f), (g), and (i).

Comput Geosci (2020) 24:17–35

27

28

Comput Geosci (2020) 24:17–35

49.72 cmand tr = 42.19 cm, depending on whether 휃s or 휃eis used in the definition of velocity.

Kalashnikov's hydrograph
120
100

The difference between the analytical and HYDRUS results is
reduced when the two models physically converge, i.e., when
the impact of soil’s capillarity is minimized. In this subsection,
we consider the same dam geometry as in Fig. 7a-g, but withthe dam material being sand (휃s=0.43, 휃r=0.045, and the van
Genuchten [40] constants 훼=0.145 cm-1 and n=2.68, k=0.495
cm/min). In this case, capillarity plays a relatively minor rolein the Richards equation. We assume that the reservoir level
rises with a rate of R =25 cm/min, the initial pressure head in
the rectangle is pd=-100 cm (this is very dry for sand).
Barenblatt (see PK-77 for full details) obtained anotheranalytical similarity solution for the problem of a “linear
drawup” of the reservoir water level with the propagation ofa wetting front into a dry dam without capillarity at relatively
small values of t, when the “saturated tongue” does not reachthe confining (right) barrier (D1F1 in Fig. 2a) at the other endof the domain. In Barenblatt’s solution, the dimensionalpﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
Darcian velocity, ua ¼ k R θs , is purely horizontal and constant in the entire flow domain. The time-expanding saturated
domain is exactly a right triangle, the hypotenuse of which is a“tilted” phreatic surface with a constant slope. The locus of thefront, i.e., the intersection of the phreatic line with the impermeable base of the dam (D1A1 in Fig. 2a) is
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
X Wa ðt Þ ¼ t k R=θs . It is noteworthy that Barenblatt’s“expanding triangle” is an exact solution to both the Laplace
equation in a potential transient 2-D flow model and to theBoussinesq 1-D equation, which is based on the DF
approximation.
HYDRUS isobars and isotachs for the selected example are shown in Fig. 7h and i, respectively. The positionof the front after 10 min of flooding is XWn≈56 cm (Fig.
7h) according to HYDRUS, and X wa =53.9 cm and59.6 cm according to the Barenblatt analytical solution,
depending on the choice of porosity or “effective” porosity, respectively. At the same time, the magnitude of velocity, although not exactly constant, in the bulk of therectangular trough is un ≈ 2 cm/min according to the numerical solution (Fig. 7i) and ua=2.3 cm/min and 2.1 cm/
min according to the Barenblatt analytical solution, depending on whether θsor θe is used, respectively. These
comparisons cross-verify that the phreatic surface model
based on the analytical solution and the numerical solution of the Richards equation provide similar results. This
substantiates a combined utilization of HYDRUS and analytical models.

80

hr(cm)

4.3 The Barenblatt’s analytical solution versus
HYDRUS simulations for the “saturated tongue”
without a barrier impedance

60
40
20
0
0

100

200

300

400

500

600

t (min)Fig. 8 Kalashnikov’s one-peak hydrograph

5 Transient 2-D flow, drawup-drawdownhydrograph
A typical hydrograph (see, for example, Fig. 8), especially forrecharge dams in Oman [64], is a one peak curve, which rapidly
rises from a zero-level (the so-called “dry dams” which have
empty reservoirs prior to flash floods) and relatively slowly recedes. Under these type of hydrological drivers, the stability ofthe dam’s soil requires a special examination [65]. The whole
porous body of the dam, rather than a soil column, should betested under rising-receding hydraulic gradients and pore pressures. In this Section, we extend the analysis from Section 4, i.e.,
we compare an analytical solution with HYDRUS simulations
for a single-maximum hydrograph and a two-layered flow
domain.

5.1 Kalashnikov’s analytical solution
We consider a vertical cross-section, a rectangle A1E1F1D1 (Fig.
9a) made of a homogeneous material with the hydraulic conductivity k and porosity 휃s. However, unlike as in Section 4, aporous rectangle A1E1F1D1 is underlain by a layer A1D1DLAL,
a low-permeable medium, an “aquitard” of the conductivity
k1<<k (Figs. 2a and 9a). The thickness f of this layer is small,and the pressure head p along ALDL is zero. In the analyticalsolution, we introduce a Cartesian coordinate system xy with theorigin A1. Both porous rectangles (Fig. 9a) are initially (t<0) dry.
After t=0, the left face of the domain, A1E1, is assigned thepressure head boundary condition corresponding to the reservoir
level hp(t), varying as a one-maximum function shown in Fig. 8.
A saturation “tongue” propagates into A1E1F1D1. In terms of the
DF model, the thickness hB(t,x) of the “tongue” A1UWB (Fig. 9a)obeys the following BVP to the Boussinesq equation:
Fig. 9 The Kalashnikov’s straight-line phreatic surface evolving in abvertical cross-section of a two-layered soil displayed together with the
HYDRUS FE mesh and boundary conditions (a), the HYDRUScomputed isobars at t=40 (b) and t=140 min (c), and the HYDRUScomputed isotachs at t=100 min (d)Comput Geosci (2020) 24:17–35

29

30


∂hB
∂
∂hB
hB
¼
khB
−k 1 ;∂x∂t∂x
fhB ðt; 0Þ ¼ hp ðt Þ; hB ð0; xÞ ¼ 0
θsComput Geosci (2020) 24:17–35

qðt; xÞ ¼ −khB

∂hB
;
∂x

0 < x < xWB ðtÞ

ð25Þ
where q is the seepage flow rate through a vertical crosssection x and xwB(t) is the locus of the wetting front, point WB,
where q=0. Flow is quasi-horizontal and quantified by the firstterm on the RHS of eq. (25) and indicated by a yellow arrow
in Fig. 9a. Leakage through the aquitard is quasi-vertical anddescribed by the second term on the RHS of eq. (25). The
corresponding leakage flow rate ql(t,x) is indicated by a redarrow in Fig. 9a. Equation (25) ignores capillarity and hysteresis: at any point x and at any time t as the phreatic surface hB
rises or drops, the pore water immediately fills or abandons
the corresponding interstitial space.
Kalashnikov [58] reported a simple (and well-forgotten by
geotechnical engineers) explicit analytical solution to the nonlinear degenerate parabolic PDE (25) for a special one-peakhydrograph:
hp ðτ Þ ¼ 4hpmax exp½−τ ð1−exp½−τ Þð26Þ

where hpmax is the maximum of the reservoir level and 휏=t k/
s is a dimensionless time. Kacimov [62] used this solution forassessing a potential impact of the reservoir on a subsurface
structure (e.g., a waste repository or a missile silo). Now, werevisit Kalashnikov’s solution to compare it with the
HYDRUS simulations.
We introduce the following dimensionless variables: hBd =hB/hpmax, xd = x/hpmax, fd = f/hpmax, td = tk/(θshpmax), and kr =k1/k, and drop the subscript “d”. The solution to the problem(25)-(26) is a phreatic surface:pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
hB ðt; xÞ ¼ 4exp½−tk r = f ð1−exp
r =f
ﬃﬃﬃﬃﬃﬃﬃﬃﬃ
ﬃ Þ−2exp½−tk r =f  k r = f x;
p½−tk
0 < x < xWB ¼ 2exp½−tk r = f  f =k r ; 0 < t < ∞ð27Þ
i.e., a segment of a straight line. The slope of this line isalways negative, i.e., pore water in the dam is devoured byleakage. There is no seepage back to the reservoir and therefore no threat of suffusion along A1U (Fig. 9a). Obviously, theaquitard zone adjacent to the A1AL line, where the verticalhydraulic gradient varies according to hB(t,0)/f, is most susceptible to translocation of fine particulates from the porous
medium of the dam body into the “aquitard.”
From Eq. (27), the peak of the hydrograph hpmax=1 isattained at t=f/kr*log(2); the maximum instantaneous flow rateof seepage from the reservoir into the dam qmax ðt; 0Þ ¼ 32=2
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
7 k r =f is attained at t=f/kr*log(3/2); the maximum volumeof pore water stored in the dam Vmax(t) =max[휃sh(t,0)xwB(t)/pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ
2]=16θs =27* f =k r is attained at t=f/kr*log(3); the maximum

propagation width of the wetting front into the dam xWBmax
pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ=2 f =k r is asymptotically attained at t = ∞ .

5.2 HYDRUS 2D simulations
For comparisons with the analytical solution (27), the transport domain simulated by HYDRUS is a large sand-maderectangle A1E1F1DL (0<XH<350 cm, 0<ZH<250 cm; thedam body) and the origin of the spatial coordinates is point
AL (Fig. 9a). The thin rectangle A1D1DLAL (an “aquitard”
consisting of a VGM loam) is 5 cm thick. Hysteresis in thesoil-water retention function is neglected in both rectangles.
Simulations were carried out for 500 min. A time-variable
pressure head boundary condition (26) with hpmax=100 cm isapproximated using a linear interpolation between points ateach 10 min (i.e., we used 50 linear interpolation lines toapproximate the hydrograph in Fig. 8). Parameters of the FE
discretization are as follows: 40x40 layers, 1,681 nodes, 160
1-D elements, and 3,200 2-D elements (Fig. 9a). Figs. 9bc
show the isobars around the phreatic surface p=0 at
t=40 min and t=140 min, respectively (the hydrograph peakis at t= 89 min). The colored map of isotachs is shown in Fig.
9d at t=100 min. Up to about t=200 min, the HYDRUS phreatic surface and the Kalashnikov solution (27) match well.
After that, capillarity of the soil in the rectangles prevails,
and the numerical and analytical solutions diverge.

6 3-D saturated-unsaturated flow simulatedby HYDRUS
In this Section, we consider a 3-D transient saturatedunsaturated flow in a rectangular loam-made parallelepipedof sizes a*b*c (Fig. 2a), where a=50 m, b=200 m, c=110 m,
and d=15 m. The real dams (see, for example, Fig. 1-49 for theOroville Dam in [66]) are zoned, with a core, sloping shoulders, drains, and other elements, which have spatially varying
and contrasting soil properties. The initial conditions for totalpressure heads and/or water contents, which reflect transientconditions of infiltration (determined by rainfall) and reservoir
stage variations, may be rather complex. Therefore, a homogeneous flow domain (Fig. 2a) and simple initial and boundary conditions adopted in our HYDRUS-3D [56] simulations
given below serve as a simple qualitative analysis and as an
asymptotic (at large t) comparison with the analytical solution
described in Section 3.
We count the total head from the dam base y=0. We
assume the following initial conditions: at t=0 the total head
h=H0=2 m in the parallelepiped 0<x1<a, 0<x2<b, 0<y<H0(hydrostatic conditions, i.e., no seepage at t<0). Above this“thin” parallelepiped, which is at an initial “background”

31

Comput Geosci (2020) 24:17–35

full saturation, the dam’s loam is at negative hydrostaticpressure heads.
We assume the following boundary conditions (see Fig.
2a):
1. The rectangles A1E1F1D1 (the plane of symmetry, thefront face of the flow domain), A1A2D3D1 (the dam base),
A 2B 3F 3D 3 (the back face of the flow domain) andD1F1F2D2 (the vertical barrier, a shaded area in Fig. 2a)are no-flow boundaries.
2. The rectangle D2C2C3D3 (the lower right face of the flow
domain) remains at a constant total head h=H0 =2 m (astatic tailwater level) at all times.
3. The rectangle C2F2F3C3 (the upper right face of the flow
domain) is a potential seepage face, which transmits water
to the atmosphere from the 3-D flow domain (an actual
seepage face) when the pore pressure on its outflow side
becomes positive and is impermeable otherwise.
4. The “expanding” rectangle A1A2B3B1 (the lower left faceof the flow domain) is the surface of a contact between therising and (eventually) static reservoir water level. We
assume that the total head initially varies linearly:h ¼ H 0 þ R* t at 0 < t < T rise

ð28Þ

where R*=1 m/day is the rate of the reservoir water level
rise and Trise=100 days is the duration of the rise. At t>Trise,
the water level in the reservoir remains constant and therefore
the boundary condition along a stabilized rectangle A1A2
B3B1 is h=102 m at t>100 days.
5. The rectangle B1E1E2B3 (the middle left face of the flow
domain), which “shrinks” at 0<t<Trise is impermeable.
6. The rectangle E1E2F3F1 (the upper face of the flow domain) is a no-flow zone.
HYDRUS [56] simulates saturated-unsaturated flow
governed by Richards’ equation in the 3-D transport domain,which is discretized into finite elements of triangular prisms.
We used default HYDRUS parameters, such as iteration
criteria, time-step controls, the VGM soil hydraulic functions,the built-in Soil Catalogue [67], etc.
We selected the “3-D simple” domain type with the following dimensions: 50x200x110 m in the XH , YH , and ZH
Cartesian coordinates of HYDRUS (Fig. 10a), with the origin
at point A 1 (see also Fig. 2a where the corresponding
coordinates were x1, x2, and y). Triangular prisms were selected for the FE discretization of the parallelepiped with 41, 27,
and 87 layers in the XH, YH, and ZH directions, respectively.
The total number of nodes was 96,309, and the finite element
mesh involved 132, 2,080, and 178,880 1-D, 2-D, and 3-D
elements, respectively. Seven observation points were placedat characteristic locations of the parallelepiped for monitoring
the dynamics of flow parameters.
We assumed the following initial (t=0 d) conditions: thetotal head (counted from the ZH=0 plane in Fig. 10a) was
H 0 =2 m in the entire transport domain (0<X H <50 m,
0<YH<200 m, 0<ZH<110 m) (i.e., hydrostatic conditions, withno seepage at t<0). The transport domain is thus initially fully
saturated for ZH<2 m and has negative hydrostatic pressure
heads for ZH>2 m. The simulation was carried out for 0<t<365
days. We recall that the boundary condition along the reservoir
face XH=0 was time variable. A snapshot of t=365 days was
selected for visualization of the results.
We assumed that the dam is made of a homogeneous, isotropic, incompressible loam with the following van
Genuchten soil hydraulic parameters selected from the
HYDRUS Soil Catalogue: the saturated hydraulic conductivity k1=0.25 m/day, porosity of 0.43, the irreducible volumetricwater content of 0.078, and the van Genuchten parameters of
the soil water retention curve 훼 = 3.6 1/m and n=1.56. We
assumed no rainfall on the ground surface ZH=110 m and noroot water uptake from the dam soil.
Simulations corroborate our concerns about the seepageinduced threat to the dam’s stability. For example, Fig. 10gshows that along the major area of the seepage face the magnitude of the hydraulic gradient exceeds 1. Moreover, close to
the edge of the barrier (the vertical line XH=50 m, YH=15 m)
extremely high gradients (exceeding 10) are apparent, similarly as obtained by the analytical solution given in Fig. 4b.

7 Concluding remarks
Factors around (actual and potential) failure of dams’
spillways can be assessed with both analytic and numerical tools, and the paper presented what we learn from thatanalysis. Our objective was to use computer algebra and a
FEM model to study the fields of the hydraulic gradients(or Darcian velocities), total hydraulic heads, and porepressures within a porous dam body, in particular, in thelocus of static and transient phreatic surfaces. We developed new analytical solutions, revisited the old ones, andsurveyed the methods as they pertain to the problems of
seepage in earth dams. We used and cross-compared theresults of computations in:&&&analytical models for quasi 3-D steady DF flow, potential
transient 2-D phreatic surface seepage governed by
Laplace’s equation in a dam body resting on an impermeable substratum;
an analytical model for 1-D Boussinesq flow over a leaky
bottom;a numerical FEM model for saturated-unsaturated flows
governed by Richards’ equation.

32

Comput Geosci (2020) 24:17–35

Fig. 10 The 3-D HYDRUS transport domain with finite element
discretization (a), the distribution of the pressure head p at t=100 d (b),the positive pressure heads at t=365 d at the shell of the transport domain(c), isobaric contours around the phreatic surface displayed on the shell of
the transport domain at the left (YH=0) and back (XH=0) planes (d) andright (YH=200 m) and front (XH=50 m) planes (e), contours of volumetricwater contents in the capillary fringe above the phreatic surface at thefront plane (XH=50 m) at t=365 d (f), isotachs at the front plane (XH=50m) at t=365 d (g), and isobars, isotachs, and velocity vectors in the crosssectional plane (ZH=2 m) at t=365 d (h)The numerical and analytical solutions presented in ourpaper, generally speaking, match well. They indicate that
under selected initial and boundary conditions, the dam
soil adjacent to an impermeable barrier is indeed under ahigh risk of erosion due to seepage. As parametric inputs
to our models, we assumed simple flow geometries (rectangles and box-shaped domains) and soil homogeneitywith respect to the saturated hydraulic conductivity andVGM capillary constants. The real dam body went
through a decade-long consolidation-differential

settlement that in terms of our models would have resulted in the heterogeneity of k, θs, θr, and VGM parameters,
the anisotropy of k, and perhaps, hysteresis (not only inwater retention curves but perhaps also in relative permeability functions). Consequently, more realistic geometries(for example, a truncated prism as the dam body, i.e., atrapezium rather than a rectangle in a vertical cross-section) and heterogeneity of soil (rock) filling can bemodeled in the future by HYDRUS. For this extended
modeling, as well as for validation of both analytical

33

Comput Geosci (2020) 24:17–35

Fig. 10 (continued)

and numerical solutions (and for falsification of Hegel’s
skepticism about historical lessons), in the application to
the particular case of the Oroville Dam, the following
information would be beneficial:&systematic (regular) historic records from piezometers/
moisture content data loggers with interpolated dynamic
fields of 3-D transient pressure heads in the interior of thedam;several infiltration tests on the dam slopes and crest for theassessment of hydraulic properties of the soil surface and,
consequently, the prediction of accretion to the phreatic
surface during rainfall events and of evapotranspiration
during dry spells;
several geophysical surveys (e.g., by OhmMapper, gravimeters, sap flow meters used for the plants emerging on
the downstream slope of the dam) with transects paralleland perpendicular to the dam crest for proxy-data on the
moisture content/pore pressure distributions, plant roots’
topology, and position of the phreatic surface in the dam
body, the following phenomena.

&&

The modeling-data acquisition/processing program for theOroville dam can serve as a prototype for many othermonitoring-maintenance programs applied to hydraulic infrastructure in the USA, and many other countries.

2.

3.

4.

5.

6.

7.

8.

9.

10.
Acknowledgments AK was supported by SQU, grant IG/AGR/SWAE/
14/02. Critique and recommendations by an Associate Editor and anonymous referee are highly appreciated.
Abbreviations BVP, boundary value problem; DF, DupuitForchheimer; PDE, partial differential equation; PK-77, reference to thebook [57]; RHS, right hand side; VGM, van Genuchten – Mualem model

11.

12.
13.

14.

References
1.

Strack, O.D.L.: Groundwater mechanics. Prentice-Hall, Inc.,
Englewood Cliffs (1989)

15.

Yakimov, N.D., Kacimov, A.R.: Darcian flow under/through aleaky cutoff-wall: Terzaghi-Anderson’s seepage problem revisited.
Int. J. Numer. Anal. Methods Geomech. 41, 1182–1195 (2017).
https://doi.org/10.1002/nag.2668
France, J.W., Alvi, I.A., Dickson, P.A., Falvey, H.T., Rigbey, S.J.
and Trojanowksi, J.: Independent forensic teamreport: oroville dam
spillway incident. url: https://damsafety.org/sites/default/files/files/
IndependentForensicTeamReportFinal01–05-18.pdf 2018
Accessed September 7, 2019
Koskinas, A., Tegos, A., Tsira, P., Dimitriadis, P., Iliopoulou, T.,
Papanicolaou, P., Koutsoyiannis, D., Williamson, T.: Insights into
the Oroville dam 2017 spillway incident. Geosciences. 9(1), 37(2019). https://doi.org/10.3390/geosciences9010037 Accessed
September 7, 2019
Bea, R. G., 2017a. Preliminary root causes analysis of failuresof the Oroville dam gated spillway. Center for Catastrophic
Risk Management, University of California Berkeley, April
17, 2017
Bea, R.G.: Legislative oversight testimony report: Oroville dam.
Center for Catastrophic Risk Management, University ofCalifornia Berkeley, May. 11, 2017 (2017b)
Vahedifard, F., AghaKouchak, A., Ragno, E., Shahrokhabadi, S.,
Mallakpour, I.: Lessons fromthe Oroville dam. Science. 355(6330),
1139–1140 (2017)
Bea, R., Cobos-Roa, D., Storesund, R.: Discussion of “overview of
New Orleans levee failures: lessons learned and their impact onnational levee design and assessment” by GL Sills, ND Vroman,
RE Wahl, and NT Schwanz. J. Geotech. Geoenviron. Eng. ASCE.
135(12), 1991–1994 (2009)
Cobos-Roa, D., Bea, R.: Three-dimensional seepage effects at three
New Orleans levee breaches during Hurricane Katrina. Electron. J.
Geotech. Eng. 13, 1–26 (2008)
Steedman, R.S., Sharp, M.K.: Physical modelling analysis of the
New Orleans levee breaches. Proc. Inst. Civil Eng. Geotech. Eng.
164(6), 353–372 (2011)
National Research Council: Dam and levee safety and community
resilience: a vision for future practice. National Academies Press(2012)
Chang, D.S., Zhang, L.M.: Extended internal stability criteria forsoils under seepage. Soils Found. 53(4), 569–583 (2013)
Zhang, L.M., Chen, Q.: Seepage failure mechanism of the Gouhou
rockfill dam during reservoir water infiltration. Soils Found. 46(5),557–568 (2006)
Vahedifard, F., Robinson, J.D., AghaKouchak, A.: Can protracted
drought undermine the structural integrity of California’s earthen
levees? J. Geotech. Geoenviron. Eng. ASCE. 142(6), 02516001(2016)
Vano, J.A., Miller, K., Dettinger, M.D., Cifelli, R., Curtis, D.,
Dufour, A., Olsen, J.R., Wilson, A.M.: Hydroclimatic extremes as

34

16.
17.
18.
19.

20.

21.
22.

23.
24.
25.
26.

27.
28.
29.

30.

31.

32.
33.
34.

35.
36.

37.
38.
39.

40.

41.

Comput Geosci (2020) 24:17–35
challenges for the water management community: lessons from
Oroville dam and hurricane Harvey. Bull. Am. Meteorol. Soc.
100(1), S9–S14 (2019)
Bonelli, S. (ed.): Erosion in geomechanics applied to dams andlevees, p. 17. Wiley, London (2013)
Istomina, V. S.: Soil stability against seepage. Gosstrojizdat,
Moscow (in Russian) (1957)
Zhilenkov, V.N.: Hydroisolation properties of soils making coresand liners in large Dams. Energiya, Leningrad (in Russian). (1968)
Bendahmane, F., Marot, D., Alexis, A.: Experimental parametricstudy of suffusion and backward erosion. J. Geotech. Geoenviron.
Eng. ASCE. 134(1), 57–67 (2008)
Ke, L., Takahashi, A.: Strength reduction of cohesionless soil due tointernal erosion induced by one-dimensional upward seepage flow.
Soils Found. 52(4), 698–711 (2012)
Rönnqvist, H., Viklander, P.: Effects of suffusion in embankment
dam filters. Dams Reservoirs. 26(3), 118–125 (2016)
Sail, Y., Marot, D., Sibille, L., Alexis, A.: Suffusion tests on cohesionless granular matter: experimental study. Eur. J. Environ. Civ.
Eng. 15(5), 799–817 (2011)
Skempton, A.W., Brogan, J.M.: Experiments on piping in sandy
gravels. Géotechnique. 44(3), 440–460 (1994)
Crank, J.: Free and moving boundary problems. Clarendon Press,
Oxford (1984)
Design of Small Dams: USBR-ISDI. Washington (1987)
Engemoen, W.: Methods to mitigate internal erosion risks inexisting embankment dams. Proc. “6th International Conference
on Scour and Erosion”, pp. 1567–1575. Paris (2012)
Fell, R., MacGregor, P., Stapledon, D., Bell, G.: GeotechnicalEngineering of Dams. Balkema, Leiden (2005)
FEMA (Federal Emergency Management Agency): Filters for embankment dams – best practices for design and construction. (2011)
Hansen, D., Roshanfekr, A.: Assessment of potential for seepage
induced unraveling failure of flow-through rockfill dams. Int. J.
Geomech. ASCE. 12(5), 560–573 (2012)
ICOLD: Internal erosion of existing dams, levees and dikes, andtheir foundations. Bulletin. Volume 1: Internal Erosion Processes
and Engineering Assessment. ICOLD, Paris (2013)
Moran, R., Toledo, M.A.: Research into protection of rockfill dams
from overtopping using rockfill downstream toes. Can J. Civ. Eng.
38, 1314–1326 (2011)
Nichiporovich, A.A.: Dams of Local Materials. Strojizdat, Moscow(1973) (in Russian)Peter, P.: Canal and River Levees. Elsevier, Amsterdam (1982)
Wan, C.F., Fell, R.: Assessing the potential of internal erosion andsuffusion in embankment dams and their foundations. J. Geotech.
Geoenviron. Eng. ASCE. 134(3), 410–407 (2008)
Wolff, T.F.: Performance of levee under seepage controls: A criticalreview (No. ERDC/GSL-TR-02-19). Vicksburg Lab (2002)
Troyanowski, J.: Dam safety: Evaluating spillway condition. http://
www.hydroworld.com/articles/hr/print/volume-27/issue-2/
technical-articles/dam-safety-evaluating-spillway-condition.html(2008) Accessed September 7, 2019
Casagrande, A.: Seepage through dams. J. N. Engl. Water Works.
51, 295–336 (1937)
Cedergren, H.R.: Seepage, drainage and flow nets. Wiley, NewYork (1989)
Rochim, A.,Marot, D., Sibille, L., Thao Le, V.: Effects of hydraulicloading history on suffusion susceptibility of cohesionless soils. J.
Geotech. Geoenviron. Eng. ASCE. 143(7), paper 04017025 (2017)
van Genuchten, M.T.: A closed-form equation for predicting thehydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J.
44, 892–898 (1980)
Chapius, R.P., Aubertin, M.: A simplified method to estimate saturated and unsaturated seepage through dikes under steady- state
conditions. Can. Geotech. J. 38, 1321–1328 (2002)

42.
43.

44.
45.

46.

47.

48.

49.
50.

51.
52.

53.

54.

55.
56.

57.

58.

59.
60.
61.
62.

63.

Desai, C.S.: Seepage analysis of earth banks under drawdown. J.
SoilMech. Found. Eng. Div. ASCE. 98(SM11), 1143–1162 (1972)
Federico, F., Calzoletti, F., Montanaro, A.: Steady state seepage
flow through zoned earth structures affected by permeability defects. In: Numerical Methods in Geotechnical Engineering, pp.
311–316. Taylor & Francis, London (2010)
Freeze, R.A.: Influence of the unsaturated flow domain on seepage
through earth dams. Water Resour. Res. 7(4), 929–941 (1971)
Ilyinsky, N.B., Kacimov, A.R., Yakimov, N.D.: Designing theshape of soil slopes stable during seepage in Californian hillsides.
In: Chandler, R.J. (ed.) Slope stability engineering. proc. of theinternational conf. on slope stability, pp. 67–70. Thomas Telford,
London (1991)
Ilyinsky, N.B., Kacimov, A.R., Yakimov, N.D.: Analytical solutions of seepage theory problems. Inverse methods, variational theorems, optimization and estimates (A review). Fluid Dyn. 33(2),157–168 (1998)
Kacimov, A.R., Obnosov, Y.V.: Analytical solutions for seepage
near material boundaries in dam cores: the Davison-Kalinin problems revisited. Appl. Math. Model. 36, 1286–1301 (2012). https://doi.org/10.1016/j.apm.2011.07.088
Lee, S.S., Yamashita, T.: Finite-difference method three dimensional model for seepage analysis through Fordyce dam. In Advances inGeosciences. Hydrological Science. 6, 171–180 (2007)
Liggett, J.A., Liu, P.L.-F.: Unsteady interzonal free surface flow inporous media. Water Resour. Res. 15(2), 240–246 (1979)
Neuman, S.P., Witherspoon, P.A.: Finite element method of analyzing steady seepage with a free surface. Water Resour. Res. 6(3),889–897 (1970)
Rehbinder, G.: Relaxation of pore pressure in a slender core of a
rockfill dam. J. Hydraul. Res. 35(2), 161–176 (1997)
Stark, T.D., Jafari, N.H., Zhindon, J.S.L., Baghdady, A.:
Unsaturated and transient seepage analysis of San Luis dam, p.
04016093. J. Geotech. Geoenviron. Eng. ASCE (2016)
Vandenboer, K., van Beek, V., Bezuijen, A.: 3D finite element
method (FEM) simulation of groundwater flow during backwarderosion piping. Front. Struct. Civ. Eng. 8(2), 160–166 (2014)
Alvi, I. A., 2013. Engineers need to get real, but can't: The role ofmodels. In Structures Congress 2013: Bridging Your Passion withYour Profession, 916–927
Kacimov, A.R., Yakimov, N.D.: Moving phreatic surface in a porous slab: an analytical solution. J. Eng. Math. 40, 399–411 (2001)Šimůnek, J., van Genuchten, M.T., Šejna, M.: Recent developments
and applications of the HYDRUS computer software packages.
Vadose. Zone. J. 15(7), (2016). https://doi.org/10.2136/vzj2016.
04.0033
Polubarinova-Kochina PY (1962, 1977) Theory of ground water
movement. Princeton University Press, Princeton. Second editionof the book in Russian is published in 1977, Nauka, Moscow
Kalashnikov, A.S.: Some problems of the qualitative theory of nonlinear degenerate second-order parabolic equations. Russ. Math.
Surv. 42, 135–176 (1987)
Haitjema, H.M.: Analytic element modelling of groundwater flow.
Academic Press, San Diego (1995)
Wolfram, S.: Mathematica. A system for doing mathematics bycomputer, Addison-Wesley, Redwood City (1991)
Blunt, M.J.: Multiphase flow in permeable media: a pore-scale perspective. Cambridge University Press, Cambridge (2017)
Kacimov, A.R.: Estimation and optimization of transient seepage
with free surface. J. Irrigation and Drainage Engrg. 119(6), 1014–
1025 (1993)
Kacimov, A., Obnosov, Y. V. and Or, D.: Evaporation induced
capillary siphoning through hydraulically connected porous domains: the Vedernikov-Bouwer model revisited. Transp Porous
Med (Springer) in press (2019) 129: 231. https://doi.org/10.1007/
s11242-019-01285-z

35Comput Geosci (2020) 24:17–35
64.

Kacimov, A.R., Brown, G.: A transient phreatic surface mound,evidenced by a strip of vegetation in an earth dam shoulder.
Hydrol. Sci. J. 60(2), 361–378 (2015)
65. Jia, G.W., Zhan, T.L.T., Chen, Y.M., Fredlund, D.G.: Performanceof a large-scale slope model subjected to rising and lowering waterlevels. Eng. Geol. 106(1–2), 92–103 (2009)
66. Jansen, R. B.: Dams and public safety. Available from U. S.
Government Printing Office, Washington D. C. 20402 Stock,(024–003). (1980)

67.

68.

Carsel, R.F., Parrish, R.S.: Developing joint probability distributions of soil water retention characteristics. Water Resour. Res. 24,
755–769 (1988)
USACoE: National Inventory of Dams. http://nid.usace.army.mil/
cm_apex/f?p=838:1:0::NO (2016)

Publisher’s note Springer Nature remains neutral with regard tojurisdictional claims in published maps and institutional affiliations.