Acta hydrotechnica 36/65 (2023), Ljubljana                                                                                       Open Access Journal 
ISSN 1581-0267                       Odprtodostopna revija 
95 
UDK/UDC: 532: 556.32(282) Prejeto/Received: 18.12.2023 
Izvirni znanstveni članek – Original scientific paper Sprejeto/Accepted: 23.02.2024 
DOI: 10.15292/acta.hydro.2023.06 Objavljeno na spletu/Published online: 18.03.2024 
TWO-DIMENSIONAL UNCONFINED SEEPAGE FLOW TOWARD A HIGHWAY CUT SLOPE  
DVODIMENZIONALNI TOK PODZEMNE VODE S PROSTO GLADINO PROTI 
AVTOCESTNEMU VKOPU 
Yebegaeshet T. Zerihun1 
1 David & James – Engineering and Environmental Consultancy, 204 Albion Road, Victoria 3350, Australia  
Abstract 
The unconfined gravity-flow system near a free-outflow boundary such as a highway cut slope was 
investigated by using a higher-order numerical model. Unlike the Dupuit–Forchheimer equation, which is 
applicable mainly to a hydraulic flow problem, the proposed model accounts for the effects of the vertical 
component of the flow for a full treatment of the problem of plane phreatic flow. The model equations were 
numerically solved by means of a finite-difference scheme. Their accuracy was then verified using the 
solutions of the full two-dimensional potential-flow method and rigorous mathematical approaches, revealing 
that for a face slope flatter than 70º, the differences between the solutions of the model and the earlier 
approaches for the relative seepage-face height were nearly negligible. The comparison results also 
demonstrated the substantial effects of the slope of the downstream face on this height and the seepage 
discharge. Furthermore, the accuracies of the model predictions for the phreatic-surface profile and the 
distributions of the piezometric head at different vertical sections are much better than the earlier method, 
which approximates the trapezoidal-shaped aquifer by its equivalent rectangular one. Such a satisfactory 
performance may be attributed to the model’s higher-order correction factor for the effects of the phreatic-
surface curvature and steep slope. 
Keywords: unconfined seepage flow, phreatic surface, artificial cut, seepage discharge, non-hydrostatic 
groundwater flow, seepage-face height. 
Izvleček 
Z numeričnim modelom višjega reda je bil raziskan gravitacijski tok pri prosti gladini v bližini meje prostega 
odtoka, kot je na primer vkop avtoceste. Za razliko od Dupuit–Forchheimerjeve enačbe, ki se uporablja 
predvsem za hidravlični tok, predlagani model upošteva vplive vertikalne komponente toka za celovito 
obravnavo problema ravninskega toka vode pri prosti gladini. Enačbe modela toka so bile numerično rešene s 
pomočjo diferenčne sheme. Natančnost rešitve je bila preverjena ob upoštevanju dvodimenzionalnega toka in 
s strogimi matematičnimi pristopi. Za naklone vkopa, ki je položnejši od 70º, so razlike v višini vode na 
precejnem robu med rešitvami tega modela in prejšnjimi pristopi skoraj zanemarljive. Primerjava rezultatov 
 
1 Stik / Correspondence: zyebegaeshet@gmail.com  
© Zerihun Y.; This is an open-access article distributed under the terms of the Creative Commons Attribution – NonCommercial – 
ShareAlike 4.0 Licence. 
© Zerihun Y.; Vsebina tega članka se sme uporabljati v skladu s pogoji licence Creative Commons Priznanje avtorstva – 
Nekomercialno – Deljenje pod enakimi pogoji 4.0.  
 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
96 
je prav tako pokazala na znatne vplive naklona vkopa na višino vode na precejnem robu in odtok. Poleg tega 
je natančnost modelne napovedi za potek gladine podzemne vode in porazdelitev hidravličnih višin na 
različnih navpičnih odsekih veliko boljša od prejšnje metode, ki vodonosnik trapezne oblike upošteva kot 
ekvivalentni pravokotni vodonosnik. Zadovoljivo rešitev je možno pripisati modelnemu korekcijskemu 
faktorju višjega reda za učinke ukrivljenosti gladine podzemne vode in strmega pobočja. 
Ključne besede: precejanje pri prosti gladini, gladina podzemne vode, vkop, odtok, nehidrostatsko precejanje, 
višina vode na precejnem robu. 
 
1. Introduction 
The analysis of unconfined saturated flow began 
with the work of Dupuit (1863), who developed a 
simple but useful model of groundwater flow by 
considering only the horizontal flow conditions and 
ignoring the processes occurring in the vertical 
direction. The Dupuit approach assumed that the 
horizontal length scale of a typical unconfined 
aquifer is much larger than the vertical length (Bear 
et al., 1968, p. 191). As a result, the energy losses 
associated with the horizontal flow process 
dominate the total losses within the saturated 
aquifer. The Dupuit approach was later generalized 
by Forchheimer (1886, 1914). It is obvious that the 
Dupuit–Forchheimer (DF) approach is applicable 
mainly to a hydraulic flow problem where the 
streamlines are nearly horizontal. In the case of an 
unsaturated-saturated gravity-flow system near the 
free-outflow boundary such as a stream bank or a 
highway cut slope, the vertical component of the 
flow cannot be ignored. For this flow situation, the 
DF approach does not predict the presence of a 
seepage-face boundary at the groundwater-surface 
water interface (see, e.g., Muskat, 1946, p. 362). 
The existence of a seepage face provides the 
necessary transition between the internal phreatic 
surface and the external equipotential boundaries. 
Several investigators (e.g. Hall, 1955; Vauclin et al., 
1979; Billstein, 1998; Chapman and Ong, 2006) 
utilized the physical model tests to better understand 
the characteristics of the unconfined flow near a 
seepage-face boundary. The results of such tests 
were also used to fine tune the performance of 
various types of numerical models. The 
development of a higher-order numerical model 
capable of simulating a seepage-flow scenario, 
where vertical flow component and seepage-face 
formation are significant, is the task of this study. 
In the past, many researchers proposed higher-order 
approaches to analyzing unconfined groundwater-
flow problems. Dagan (1967) was probably the first 
to develop second-order governing equations for 
unsteady groundwater flow by applying a procedure 
expanding in powers of a small parameter. His 
approach accounted for the effects of the vertical 
component of the flow. Based on the assumption of 
a quadratic variation of a piezometric (hydro-
potential) head across a vertical section, Kashef 
(1969) proposed simplified equations for a steady 
unconfined groundwater flow by analyzing the 
hydrodynamic forces within the flow medium. The 
results of the equations for discharge and seepage-
face height were compared with rigorous 
mathematical solutions, demonstrating improved 
accuracy over the DF approach. For a homogeneous 
and isotropic aquifer with a sloping downstream 
face, however, his method provides neither 
analytical nor numerical solutions to the phreatic-
surface profile and the hydro-potential distribution. 
Using a similar assumption for the profile of the 
piezometric head, Knight (2005) proposed 
approximate solutions to the unconfined seepage-
flow problems. Muleshkov and Banerjee (1987) 
employed conformal mapping techniques to obtain 
a solution to the problem of the two-dimensional 
(2D) steady-state seepage flow. However, the 
practical application of their analytical solution is 
restricted by the assumption of infinite flow 
domains. As an alternative, Chapman and Dressler 
(1984) proposed governing equations in a bed-fitted 
curvilinear coordinate system for a transient 
unconfined groundwater flow over curved 
impermeable bedrock. Using the assumption of a 
small bedrock curvature, Chapman and Ong (2006) 
developed a simplified version of the Chapman and 
Dressler flow profile equation for a non-hydrostatic 
shallow groundwater flow. Their equation can be 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
97 
retrieved from the Hilberts et al. (2004) equation by 
considering flow in a constant-width aquifer. Since 
the bedrock curvature is incorporated without 
resorting to higher-order terms in the shallowness 
expansion, the resulting governing equations of 
these studies do not separate the effects of the 
phreatic-surface curvature from the bedrock 
curvature and become identical to the DF equation 
for the case of mild-slope planar bedrock. The 
Picard iteration method, which was originally 
proposed by Matthew (1991) for open-channel 
flow, was employed by Castro-Orgaz and Hager 
(2014) to extend the DF approach. Using the power 
series expansion of the harmonic stream function 
(Boussinesq, 1871; Rayleigh, 1876) and Darcy’s 
law, Di Nucci (2018) developed a higher-order 
differential equation for a curved unconfined flow 
over horizontal impermeable bedrock. Additionally, 
the numerical solutions of this equation for 2D 
seepage flow through a rectangular-shaped dam 
were investigated, and the results of the phreatic-
surface profile showed satisfactory agreement with 
the numerical test data. Zerihun (2018) presented a 
higher-order approach for analyzing curvilinear 
groundwater flow in a constant-width unconfined 
aquifer. His method allowed for the effect of non-
hydrostatic pressure, thereby overcoming the 
limitations of the DF approach. More recently, this 
approach was extended to deal with unconfined 
flows in convergent- and divergent-type hillslope 
aquifers with non-uniform bedrock slopes (Zerihun, 
2020). Although those equations had relatively long 
and complicated non-linear coefficients, they 
yielded reasonably accurate solutions to the 
problems of groundwater flow with a phreatic 
surface. 
This study was undertaken to address some of the 
above deficiencies by developing an alternative 
numerical model for a 2D unconfined groundwater 
flow. The field of interest is the problem of a steady-
state phreatic flow through a homogeneous aquifer 
with a slanting side cut, which is commonly 
encountered along highways and railways (see 
Figure 1). For this type of flow problem, 
approximate analytical solutions based on a 
complex potential-flow theory (e.g. Polubarinova-
Kochina, 1962, pp. 301-308; Kacimov and 
Obnosov, 2012) and the results of the empirical 
methods (e.g. Kozeny, 1931; Casagrande, 1937; 
Fukuchi, 2018) are available. From the 
computational point of view, however, the 
evaluation of some of the resulting formulas is not 
straightforward. In comparison, the proposed model 
is simple to use and has a higher-order correction 
term for the effects of the phreatic-surface 
curvature. As a result, it allows the treatment of 
complex unconfined flow behaviors that may 
influence the predictions of seepage discharge, 
phreatic-surface profiles, and seepage-face height, 
and obviates the deficiencies of some of the existing 
models. The model’s capacity for analyzing the 
salient features of a non-hydrostatic seepage-flow 
problem is also investigated by comparing its 
numerical solutions with the results of the full 2D 
potential-flow method and rigorous mathematical 
approaches. The complete solutions of the 
unconfined flow problem considered here are 
relevant for assessing the long-term stability of the 
slope of a highway cut section. In the following 
sections, the derivation of the governing equations 
and the development of the numerical solution 
procedure based on a finite-difference 
approximation will be presented. 
 
Figure 1: A steeply sloping side cut along the 
curved section of a highway (Photograph courtesy 
of ROADEX Network, Finland). 
Slika 1: Strm vkop vzdolž zakrivljenega odseka 
avtoceste (fotografija z dovoljenjem ROADEX 
Network, Finska). 
 
 
 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
98 
2. Governing equations  
The flow problem under consideration is a 
groundwater flow in an unconfined aquifer over 
curved impermeable bedrock with a steep slope. 
Figure 2 shows a schematic of the groundwater-
flow problem in a Cartesian coordinate system, 
where x  is horizontal in the streamwise direction; 
y  is vertically upward; and z  is horizontal in the 
lateral direction. The following assumptions are 
made to derive the equations using the Matthew 
(1991) step-wise iterative approach: (1) the 
variations of flow parameters in the transverse 
direction are ignored; (2) the aquifer is isotropic and 
homogeneous; and (3) the capillary fringe effects 
are ignored. 
 
Figure 2: Definition sketch of curvilinear 
groundwater flow with a phreatic surface. 
Slika 2: Definicijska skica krivolinijskega toka 
podzemne vode z gladino podzemne vode. 
The first approximation corresponds to the lowest-
order equation of an unconfined groundwater flow 
in a sloping aquifer. For such a gradually-varied 
flow, the expression for the piezometric or hydraulic 
head is given by (Zerihun, 2018) 
 ( )
( )
,
1
,, Y
HH
tx +

+

−
=

  (1a) 
 ,
Y
Yy
−
−
=

   (1b) 
where   refers to the phreatic-surface elevation; H  
is the saturated thickness of the aquifer; y  is the 
elevation of a point in the flow field; Y  is the 
bedrock elevation; ( ) /py +=  is the piezometric 
head or the potential function; p  is the pressure;   
is the unit weight of the fluid;   is the 
dimensionless vertical height; t  is the time; and 
( ) xxx YHY ++=
2
1 . 
Applying Darcy’s law to Equation (1a) results in 
 ,





+

−=


−= x
x Y
H
K
x
Ku

 (2) 
where u  is the velocity in the streamwise x -
direction, and K  is the average hydraulic 
conductivity of a saturated aquifer. In the above and 
succeeding equations, the subscript x  denotes 
partial differentiation with respect to the horizontal 
x -axis. Subsequent corrections can be applied to 
Equation (2) using the following Cauchy–Riemann 
conditions (Bear, 1988, p. 234): 
 ,
yx
Ku


−=


−=

  (3) 
 ,
xy
Kv


=


−=

  (4) 
where   is the stream function, and v  is the 
vertical velocity. For simplifying the algebraic 
manipulation, the vertical height above the bedrock 
( )Yyh −=  is used hereafter (see Figure 2). Based 
on this simple approach, Equation (3) is integrated 
with the aid of Equation (2) to give the following 
first-order approximation for the stream function:  
( ) ( ),,,, 1 txY
H
KHudhtx x
x  +





+

=−=   (5) 
where 1  is a function of integration. Imposing the 
boundary condition at the bedrock ( )0=  results in 
01 = . Differentiating Equation (5) with respect to 
x  yields an equation for the vertical velocity as 
follows: 
.





+

−





+

=


= x
x
xxx
xx Y
H
KYY
H
KH
x
v 

 (6) 
Substituting the above equation into Equation (4) 
and integrating the resulting expression with respect 
to h  gives the following higher-order 
approximation for the potential function: 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
99 
 
( )
( )
,
,
2
1
,,
2
2
2
K
tx
Y
H
HY
Y
HH
vdh
K
tx
x
x
x
xx
xx



−





+

+






+

−=−= 
 (7) 
where 2  is a function of integration. According to 
Matthew (1991), two approaches can be utilized to 
determine the expressions for this function. The first 
one is by imposing the dynamic boundary condition 
at the phreatic surface, ( ) 01 ==p , which results 
in the following equation for  : 
 
( )
( ) .1
2
1
,,
2
2






+

+−





+

−







 −






+

=
Y
H
Y
H
HY
Y
H
Htx
x
x
x
xx
xx



 (8) 
This higher-order equation accounts for the 
curvilinear nature of the phreatic flow in the vertical 
plane, and hence it can accurately describe the 
distribution of the piezometric head across the 
saturated depth. The second approach applies the 
boundary conditions of the stream function at the 
bedrock and phreatic surface. For this purpose, 
Equation (7) is first inserted into Equation (3) and 
then differentiated with respect to x . This results in 
a higher-order approximation for the horizontal 
velocity component as follows: 
 
.
2
2
22
2
2
x
Y
H
KY
Y
H
KHY
Y
H
KHY
Y
HKH
x
Ku
x
x
x
x
x
xx
xx
xx
x
xxx
xxx


+





+

+






+

−






+

−






+

=


−=





 (9) 
Integration of Equation (9) with respect to h  leads 
to the desired higher-order approximation for the 
stream function, i.e. 
( )
( ),,
2
6
,,
3
22
22
3
3
txH
x
Y
H
KHY
Y
HY
Y
H
YKH
Y
HKH
udhtx
x
x
x
x
xxx
xx
xx
x
xxx
xxx





+


+





+

−












+

+





+

+






+

−=−= 
 (10) 
where 3  is also a function of integration. The 
expressions for 3  and x 2  are obtained by 
applying the boundary conditions ( ) 00 ==  and 
( ) q−==1 , respectively, and are given by: 
 ,03 =   (11) 
,0
2
6
2
2
=+





+







−+












+

−





+

−=


H
q
Y
H
Y
HY
K
Y
H
YY
HH
KH
x
x
x
x
xx
xx
xx
xxxx
xxx
 
   (12) 
where q  is the discharge per unit width. 
Substituting the above equations into Equation (10) 
and simplifying the resulting expression yields the 
following equation: 
( ) ( )
( )
( ) .
2
6
,,
2
2
22
3
3



qY
HYKH
Y
H
YKH
Y
HKH
tx
x
xxx
xx
xx
x
xxx
xxx
−−





+

+
−





+

+
−





+

=
 (13) 
By making use of Equation (12) in Equation (9), an 
expression for the distribution of the horizontal 
velocity is obtained as follows: 
 
( )
( )
( ) .21
2
21
13
6
2
2
H
q
Y
HKHY
Y
H
KHY
Y
HKH
u
x
xxx
xx
xx
x
xxx
xxx
+−





+

+
−





+

+
−





+

=



 (14) 
Similar to the results of previous studies (e.g. 
Chapman and Dressler, 1984; Di Nucci, 2018; 
Zerihun, 2018, 2020), this equation includes terms 
that account for the effects of the vertical curvature 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
100 
of the streamline. Substituting Equation (8) in 
Equation (3) and equating the resulting expression 
to Equation (14) yields the following flow profile 
equation for an unconfined groundwater flow: 
( )
.0
2
3
2
2
=+





+

+






+

+−





+

−






+

+





+

KH
q
Y
H
Y
H
YYHY
HHY
Y
H
HHY
HH
x
x
x
x
xxxx
xxx
xx
xx
xxxx
xxx
 (15) 
The above equation implicitly incorporates not only 
the effects of non-uniform horizontal velocity and 
non-hydrostatic pressure distributions but also the 
effect of a steep bedrock slope for accurately 
modeling the problems of a groundwater flow with 
a phreatic surface. Thus, it differs from previous 
equations (e.g. Castro-Orgaz and Hager, 2014; Di 
Nucci, 2018) developed by the Boussinesq (1871) 
or Matthew (1991) approach. For a 2D unconfined 
saturated flow on a vertical plane, the depth-
averaged mass-conservation equation reads as 
(Bear, 1988, p. 376): 
 ,
0
R
x
q
t
H
dhu
xt
H H
=


+


=









+


   (16) 
where   is the effective porosity of the aquifer, and 
R  is the rate of vertical recharge. For a curvilinear 
groundwater-flow problem, the variation of the 
phreatic-surface profile can be predicted by 
combining the above depth-averaged continuity 
equation with Equation (15). Additionally, the 
present approach gives complete solutions for the 
rate of seepage flow, height of the seepage surface, 
and hydro-potential distribution within the saturated 
region. In the case of phreatic flows with 
streamlines nearly parallel to the sloping planar 
bedrock, Equation (15) reduces to: 
 ,0=+





+
 K
q
Y
H
H x
x   (17) 
which is a Boussinesq-type equation developed by 
Childs (1971) and was adapted by Chapman (1980) 
for analyzing the phreatic-surface profiles. For 
horizontal bedrock, this equation degenerates to the 
DF equation for steady groundwater flow with no 
recharge. Because of the non-linear characteristics 
of Equations (15) and (16), closed-form solutions 
are not possible except in highly simplified cases. In 
this study, a numerical approach based on the finite-
difference approximations was employed to solve 
the equations. Using such a relatively simple 
numerical technique, the applicability of the 
proposed equations is systematically investigated 
for the test cases of non-hydrostatic seepage flows 
towards artificial cuts ( 1= ). Such transient types 
of flow problems may be approximately analyzed 
by a quasi-steady flow approach, ignoring the 
aquifer’s storage capacity and the time lag between 
the inflow and outflow peak discharges. 
 
3. Seepage through unconfined aquifers of 
finite base width  
Figure 3 displays a typical longitudinal profile of an 
artificial cut with horizontal impermeable bedrock. 
The cut bottom coinciding with the bedrock is taken 
as a datum for the measurements of elevation (
0=Y ). For a steady accretion-free flow, Equation 
(15) may be rewritten in the form: 
 .
32
32
K
q
H
HH
x
xx −=








+


 (18) 
Integrating the above expression with respect to x  
yields: 
 ,
23
4
23
+−=

+

x
K
qH
H
H
xx  (19) 
where 4  is an integration constant and can be 
determined from the specified flow depth at the 
upstream end section ( 0=j ) as a boundary 
condition. In practice, this section coincides with a 
groundwater divide or with a section where the 
horizontal Darcy velocity is almost negligible. 
Thus, the expression for 4  becomes: 
 .
23 0
2
0
0,
0
3
0
4

+

=
H
H
H
xx   (20) 
 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
101 
 
Figure 3: Definition sketch for unconfined seepage 
flow toward a highway cut slope. 
Slika 3: Definicijska skica za precejanje pri prosti 
gladini proti pobočju vkopa avtoceste.  
In the above equation, the subscript 0  denotes flow 
parameters at the upstream end section. As 
discussed before, a numerical approach is the only 
way of obtaining a solution to Equation (19) for 
such a curvilinear potential-flow problem. The 
spatial derivative term of this equation is discretized 
using the four-point finite-difference equation 
(Bickley, 1941), which results in the following 
expression after some mathematical manipulation: 
 
( )
( )
( ) ,
23
2
2
3
0
2
0
0,
0
3
0
2
112
3

+

+−=

++−

−+
H
H
H
xj
K
q
H
HHH
x
H
xx
j
j
jjj
j
j
 (21) 
where x  is the size of the step. With the aim of 
minimizing the numerical error due to spatial 
discretization, the step size is varied between 2 and 
4% of the horizontal length of the model domain. 
For computational nodes near the downstream end, 
the derivative term is discretized with the three-
point backward finite-difference approximation. 
Using the Newton–Raphson method with a 
numerical Jacobian matrix, the non-linear form of 
the algebraic equations is transformed into a linear 
one. The resulting implicit set of linear equations is 
then solved by means of the lower-upper (LU) 
decomposition method. A relative change in 
solution criterion with a tolerance of 610−  is used 
for examining the convergence of the numerical 
solutions. More details regarding the applied 
numerical scheme can be found in Zerihun (2018). 
3.1 Seepage discharge  
Using the Darcy law, a general expression for the 
rate of seepage flow toward a highway cut slope 
with horizontal impermeable bedrock can be 
obtained as follows (Hantush, 1962; Polubarinova-
Kochina, 1962, pp. 281-282; Kashef, 1969; 
Zerihun, 2018): 
 
( )
( ) .
2
2









−


−=


−==


H
dyy
x
K
dy
x
Kdyyuq
Y
YY




 (22) 
As noted by Casagrande (1937) and Muskat (1946, 
p. 291), the terminal end of the phreatic-surface 
profile is tangent to the downstream face of the 
aquifer at the exit point B (see Figure 3). Using such 
a qualitative characteristic of the profile of the 
groundwater flow, the horizontal component of the 
phreatic-surface velocity at B was determined by 
applying the Darcy law. Utilizing this known 
velocity, Equation (22) was integrated between the 
vertical sections at 0=j  and nj =  to give an 
expression for the seepage flux through a 
trapezoidal-shaped aquifer (for details, see the 
Appendix). The resulting equation for the seepage 
discharge reads as:  
 ,
)cot(2
22
0








−
−
=
B
B
HL
HH
Kq  (23) 
 ,sin
3
2
1 2 −=   (24) 
where   is the slope of the downstream face; BH  
is the depth of the saturated aquifer at nj = ; and L  
is the base width of the aquifer. Equation (23) is 
structurally similar to the discharge equation 
presented by Kashef (1969) after analyzing the 
hydrodynamic forces within the flow medium. This 
equation is directly coupled with Equation (19) for 
a complete numerical solution of the unconfined 
seepage-flow problem. 
In the absence of tailwater, the maximum seepage 
discharge (Polubarinova-Kochina, 1962, p. 308; 
Kashef, 1969) can be obtained by differentiating 
Equation (23) with respect to BH  and equating the 
resulting expression to zero as follows: 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
102 
 .0cot2cot 20
2 =+−  HHLH mSmS  (25) 
For a physically realistic solution, the discriminant 
of the above quadratic equation must be positive. 
Furthermore, its analytical solution satisfies the 
necessary condition 0HH B   and can be expressed 
as a non-dimensional form as follows: 
 ,
1tantan
21
2
2
0









−−==





H
HmS
m  (26) 
where ( )LH /0=  is the relative upstream depth of 
the saturated aquifer; mSH  is the maximum height 
of the seepage face for a downstream dry-bed 
condition; and m  is the relative maximum 
seepage-face height. Using the above equation, an 
expression for 20H  as a function of L  and   is 
obtained. Inserting this expression into Equation 
(23) leads to the following equation for the 
maximum seepage discharge: 
 ,tan mSm KHq =   (27) 
where mq  is the maximum seepage discharge. 
Combining Equations (23) and (27) results in an 
expression for the relative seepage discharge as 
follows: 
 ,
tan
2
1 20


















−
−








==




mSm
R
H
H
q
q
 (28) 
where R  is the relative seepage discharge, and 
( )0/ HH B=  is the relative saturated depth of the 
aquifer at nj = . Using the definition of   in 
Equation (26), the resulting expression after 
simplification becomes: 
 ,
1tan1tan
2/1
2
2









−−==




B
mS
m
H
H
 (29) 
where m  is the ratio of the maximum seepage-face 
height to the downstream saturated depth in the 
presence of tailwater. As mentioned before, 
rigorous mathematical approaches have been 
previously employed to analyze unconfined seepage 
flows. For estimating the discharge, Pavlovsky 
(1931) (Polubarinova-Kochina, 1962, p. 306) 
proposed an approximate equation, which is not a 
function of the downstream saturated depth, as 
follows: 
 .
cot11 22
0


−+
=
KH
q   (30) 
For 90 , Mikhailov (1952) systematically 
simplified the complex equation developed by 
Falkovich (Polubarinova-Kochina and Falkovich, 
1951) and introduced a relatively simple equation as 
follows: 
 ,
6cot
4
 −+
= B
KH
q   (31) 
where   is a constant, and its value depends on the 
magnitude of cot . At the maximum seepage 
discharge, equating Equation (30) to Equation (31) 
and simplifying the resulting expression yields: 
 
( )
,
cot114
6cot
22
0 



 −+
−+
==



H
HmS
m  (32) 
which differs from Equation (26) by the additional 
parameter  . Equations (30) and (32) are employed 
here to assess the performance of the proposed 
equations. 
 
3.2 Effects of the downstream slope  
As a part of this study, the effects of the 
downstream-face slope on the salient features of the 
unconfined seepage flow are examined. Figure 4a 
compares the predictions of Eq. (26) with the results 
of Eq. (32) for the relative maximum seepage-face 
height. As can be seen, the results of the proposed 
equation agree reasonably well with the results of 
the previous equation. For 70 , some 
discrepancies can be seen between the results of the 
two equations, with a maximum absolute error of 
3.6%. Considering the approximate nature of the 
simplified version of the previous equation, the 
performance of the proposed equation is 
satisfactory. For a constant  , increasing the slope 
decreases the relative seepage-face height. This 
effect is significant for 45  and 10.0 . For 
a steeper slope, the variation of the relative seepage-
face height with   is more gradual. 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
103 
 
Figure 4: (a) Variation of the relative seepage-face height m  with the slope of the downstream face  ; and 
(b) normalized seepage discharge N  as a function of the relative upstream saturated depth  . The results 
of Equations (26) and (23) in (a) and (b), respectively, are indicated by lines of different styles and colors. 
Slika 4: (a) Sprememba relativne višine precejanja m  z naklonom dolvodne meje  ; in (b) normaliziranega 
odtoka N  kot funkcije relativne gorvodne zasičene globine  . Rezultati enačb (26) in (23) v (a) oziroma (b) 
so označeni s črtami različnih slogov in barv. 
 
The effect of the slope of the downstream face on 
the normalized seepage discharge ( )0/ KHqN =  is 
also examined, and the result is depicted in Figure 
4b. The results of Equation (23) are favorably 
compared with the results of Equation (30). The 
computed maximum absolute error does not exceed 
3.4% for numerical results with   values within the 
considered range. For 30.0 , the effect of the 
face slope on the normalized discharge is 
insignificant. For a relatively higher value, 
however, increasing the slope decreases the 
normalized discharge. This is due to the effect of   
on the geometric characteristics of the phreatic-
surface profile. The overall results revealed that the 
slope of the downstream face significantly affects 
both the relative seepage-face height and the 
normalized seepage discharge. 
Likewise, the effect of the downstream-face slope 
on the relative seepage discharge is assessed using 
Equation (28). As shown in Figure 5, the reduction 
in seepage discharge decreases as   increases. For 
a larger face slope, the relative discharge varies 
more gradually with m/1 . Consequently, the drop 
in discharge is relatively small. For an artificial cut 
section with a flatter downstream-face slope, the 
converse is true. Apparently, the presence of 
tailwater affects the height of the seepage face ( SH
) and thus the relative discharge. The seepage-face 
height decreases with increasing tailwater depth (
mSS HH  ). As 1→m , the seepage discharge 
approaches the maximum value that corresponds to 
the condition of no tailwater. These results are 
consistent with previous results reported by 
Polubarinova-Kochina (1962, p. 308) for 
unconfined seepage flow through a triangular-
shaped dam. 
 
3.3 Phreatic-surface profile and piezometric 
head distributions  
In this section, the model’s accuracy is investigated 
by simulating unconfined seepage flow through a 
trapezoidal-shaped dam with vertical upstream and 
slanted downstream faces ( 45= ). For this 
problem, the slope of the phreatic surface ( 0=xH ) 
and the saturated depth were prescribed as boundary 
conditions at the upstream end section. The 
numerical technique applied here employs an 
iterative solution procedure based on the initial 
guess of the phreatic-surface profile. Using 
Equation (23), the discharge is first computed by 
utilizing the estimated depth of the saturated aquifer 
at nj =  ( mSB HH  ). If the discrepancy between 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
104 
the BH  employed to compute q  and the newly 
predicted BH  is significant, then the iteration 
process is continued until the discrepancy becomes 
insignificant (approximately less than 5 mm). 
 
 
Figure 5: Variation of the relative seepage 
discharge R  with the downstream saturated depth 
ratio m/1  for various values of downstream-face 
slope. 
Slika 5: Sprememba relativnega odtoka R  z 
razmerjem dolvodne zasičene globine m/1  za 
različne vrednosti naklona dolvodne meje.  
 
The computational result of the present model was 
compared with the result of the finite-element 
numerical model (Lacy and Prevost, 1987) and the 
result of an approximate method that replaces the 
trapezoidal longitudinal section with an equivalent 
rectangular one (Kashef, 1969). According to this 
method, the locus of the phreatic surface is 
determined by: 
( ) ,
4
exp1
2
1
2
1
2/1
2
2
0
2
0
Y
q
Kx
K
q
KH
qx
H
x +



























 −
−





+








−
=  (33) 
 ,
2
22
0
e
w
B
HH
K
q −
=   (34) 
where wH  is the tailwater depth, and eB  is the 
width of the equivalent rectangular longitudinal 
section. The equation for the width of this section is 
 .
2
cotB
e
H
LB −=   (35) 
As shown in Figure 6a, the result of the present 
model correlated reasonably well with the result of 
the finite-element numerical model with a mean 
relative error of 2.9%. Compared to these results, 
the approximate method underestimated the 
phreatic-surface elevations in the flow region near 
the sloping downstream face ( 40.0/ 0 Hx ). 
Consequently, the exit point of the line of seepage 
and the curvature of the phreatic surface were not 
accurately predicted. The result of the test 
confirmed that the present model is capable of 
accurately simulating the flow pattern in the region 
where the streamline curvature and inclination are 
significant. 
The results of the proposed model for the 
piezometric head distributions at vertical sections 
18.0/ 0 =Hx , 38.0/ 0 =Hx , and 60.0/ 0 =Hx  
were also compared with the solutions of the Kashef 
(1969) approximate method, which describes the 
vertical profile of the piezometric head by:   
 ( ) ,2 bb H −+=   (36) 
where b  is the piezometric head at the bedrock and 
can be determined from the following equation: 
 .1
4
3
4
2
0








−+=
e
b
B
x
H
HH
  (37) 
As shown in Figure 6b-d, the piezometric head at 
the bedrock decreases as the distance from the 
upstream end section increases. This is due to the 
effect of the vertical component of the flow, which 
is significant at a section near the discharge face. 
Similar to its result for the phreatic-surface profile, 
the approximate method did not accurately simulate 
the distributions of the piezometric head. It 
underestimated the head at different levels, 
especially at a vertical section near the discharge 
face where a maximum discrepancy of 9.4% 
occurred. Compared to the results of the present 
model and the numerical data of Lacy and Prevost 
(1987), the conventional DF approach, which is 
based on an approximation of the hydrostatic 
pressure distribution, overestimated the 
distributions of the piezometric head throughout the 
computational domain, i.e., 0.1/ = .
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
105 
 
 
Figure 6: (a) Normalized phreatic-surface profile for curvilinear seepage flow ( 47.0=N ). The normalized 
piezometric head distributions at vertical sections 18.0/ 0 =Hx , 38.0/ 0 =Hx , and 60.0/ 0 =Hx  are shown 
in (b), (c), and (d), respectively. 
Slika 6: (a) Normaliziran profil gladine podzemne vode za ukrivljene tokovnice ( 47.0=N ). Porazdelitve 
normalizirane hidravlične višine na navpičnih odsekih 18.0/ 0 =Hx , 38.0/ 0 =Hx  in 60.0/ 0 =Hx  so 
prikazane v (b), (c) oziroma (d).  
 
4. Summary and conclusions 
The problem of a steady-state phreatic flow through 
a homogeneous and isotropic aquifer with a slanting 
downstream face was investigated using a non-
hydrostatic model. The proposed model 
incorporates a higher-order correction term, which 
accounts for the effects of the vertical component of 
the flow. Hence, it allows the treatment of complex 
unconfined flow behaviors that may influence the 
predictions of seepage discharge and seepage-face 
height, and overcomes the shortcomings of some of 
the existing models (e.g. Kashef, 1969; Chapman 
and Ong, 2006). A numerical approach based on the 
finite-difference approximations was employed for 
the solutions of the model equations. The model’s 
capacity for analyzing the important features of the 
seepage-flow problem was investigated by 
comparing its numerical results with the solutions 
of the full two-dimensional potential-flow method 
and rigorous mathematical approaches. 
The following conclusions are drawn from this 
study: 
1. The solutions of the proposed equations for 
the relative seepage-face height and the 
normalized seepage discharge qualitatively 
agreed with the results of the previous 
simplified equations. The result of the 
comparison is quite promisimg when 
considering the complex nature of the 
problem of seepage flow, where the vertical 
velocity plays a significant role. Overall, 
the investigation results confirmed that the 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
106 
slope of the downstream face significantly 
affects both the relative seepage-face height 
and the normalized seepage discharge. 
Furthermore, the analysis of the relative 
discharge attested that the seepage 
discharge is strongly influenced by the 
downstream saturated depth ratio. 
2. For an unconfined seepage flow through a 
trapezoidal-shaped aquifer, the present 
model accurately reproduced the profile of 
the phreatic surface. In comparison, the 
Kashef (1969) approximate method 
underestimated the elevations of the 
phreatic surface in the flow region near the 
discharge face ( 40.0/ 0 Hx ). The result 
of the comparison also revealed that the 
accuracy of the model predictions for the 
distributions of the piezometric head at 
different vertical sections is much better 
than the approximate method and the DF 
approach. 
The test results presented in this study demonstrated 
the model’s capability of yielding reasonably 
accurate solutions to the problems of the unconfined 
seepage flow. The good quality of the 
computational results may be attributed to the 
model’s higher-order correction factor for the 
effects of the phreatic-surface curvature and steep 
slope. It is believed that the proposed model can 
serve as a numerical tool for solving engineering 
problems related to the stability analysis of a 
highway cut section and the design of the associated 
side channel. 
 
Appendix: Discharge equation development 
Along the streamline at the phreatic surface, if s  is 
the local coordinate tangent to the phreatic surface 
and   is the angle between the streamline and the 
x -axis at B, then the expression for the Darcy 
velocity becomes (Harr, 1991, p. 68): 
 
,cossin
cos,



K
s
K
x
Ku Bp
=


−=


−=
 (38) 
where Bpu ,  is the horizontal component of the 
phreatic-surface velocity at B. For unconfined 
seepage flow through a trapezoidal-shaped aquifer 
with horizontal bedrock ( 0=== xxxxxx YYY ), a 
general expression for the horizontal velocity at the 
phreatic surface, pu , is obtained from Equations 
(14) and (15) as: 
 ( ),1+−= HHKHu xxxp   (39) 
where xH  is the slope of the phreatic surface, which 
is equal to tan−  at B. Equating Equation (38) to 
Equation (39) gives the following expression: 
 .
sin1cos 22
,
BB
Bxx
HH
H
 −
=
−
=  (40) 
Using the condition of a nearly hydrostatic flow at 
the upstream end section, Equation (22) is 
integrated between the vertical sections at 0=j  
and nj = . The resulting expression becomes: 
 ,
232
23
,
2
00 BBBxx HHHH
K
qL
−−=  (41) 
 ,cot0 BHLL −=   (42) 
where 0L  is the length of the projection of the 
phreatic-surface profile on the x -axis (see Figure 
3). Substituting Equation (40) into Equation (41) 
and then simplifying the resulting expression gives 
the discharge equation for seepage flow as follows: 
 ,
2 0
22
0







 −
=
L
HH
Kq B   (43) 
 .sin
3
2
1 2 −=   (44) 
 
Notation 
eB  width of the equivalent rectangular-shaped 
aquifer [L] 
g  gravitational acceleration [LT-2] 
h  vertical height above the bedrock [L] 
H  depth of the saturated aquifer [L] 
BH  depth of the saturated aquifer at nj =  [L] 
mSH  maximum height of the seepage face [L] 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
107 
SH  height of the seepage face [L] 
wH  tailwater depth [L] 
xH , xxH , and xxxH  derivatives of the saturated 
depth [-, L-1, L-2] 
0H  saturated depth at the upstream end section 
[L]  
K  hydraulic conductivity [LT-1] 
L  base width of the aquifer [L] 
0L  length of the projection of the phreatic-
surface profile on the x-axis [L] 
p  pressure [ML-1T-2] 
q  seepage discharge per unit width [L2T-1] 
mq  maximum seepage discharge [L
2T-1] 
R  rate of vertical recharge [LT-1] 
s  local coordinate tangent to the phreatic 
surface [L] 
t  time [T] 
u  horizontal velocity [LT-1] 
pu  horizontal velocity at the phreatic surface 
[LT-1] 
v  vertical velocity [LT-1] 
yx, , and z  Cartesian coordinates [L] 
Y  bedrock elevation [L] 
xY , xxY , and xxxY  derivatives of the bedrock 
profile [-, L-1, L-2] 
  relative saturated depth at nj =  [-] 
m  relative maximum seepage-face height [-] 
  slope of the downstream face [deg] 
  unit weight of the fluid [ML-2T-2] 
x  step size [L] 
  phreatic-surface elevation [L] 
  effective porosity [-] 
  a constant factor [-] 
  dimensionless vertical height [-] 
  relative upstream saturated depth [-] 
  piezometric head [L] 
b  piezometric head at the bedrock [L] 
N  normalized seepage discharge [-] 
R  relative seepage discharge [-] 
  stream function [L2T-1] 
,, 21   and 3  functions of integration 
[L2T-1] 
4  an integration constant [L
2]  
m  ratio of the maximum seepage-face height 
to the downstream saturated depth [-] 
 
References 
Bear, J., Zaslavsky, D., Irmay, S. (1968) Physical 
Principles of Water Percolation and Seepage. Arid Zone 
Research - XXIX, UNSECO: Paris, France. 
Bear, J. (1988) Dynamics of Fluid in Porous Media. 
Dover Publications Inc.: Mineola, NY, USA. 
Bickley, W. G. (1941) Formulae for numerical 
differentiation. Math. Gaz., 25(263), 19–27, 
https://doi.org/10.2307/3606475. 
Billstein, M. (1998) Development of a Numerical Model 
of Flow through Embankment Dams. Licentiate Thesis, 
Department of Environmental Engineering, Luleå 
University of Technology, Luleå, Sweden. 
Boussinesq, J. (1871) Théorie de l’intumescence liquide 
appelée onde solitaire ou de translation, se propageant 
dans un canal rectangulaire (Theory of liquid wave 
referred to as a solitary or a translation wave, propagating 
in a rectangular canal). Comptes Rendus, Académie des 
Sciences, 72, 755–759 [in French]. 
Casagrande, A. (1937) Seepage through dams. J. New 
England Water Works Assoc., 51(2), 131–172. 
Castro-Orgaz, O., Hager, W. H. (2014) One-dimensional 
modeling of curvilinear free-surface flow: Generalized 
Matthew theory. J. Hydraul. Res., 52(1), 14–23. 
https://doi.org/10.1080/00221686.2013.834853.  
Chapman, T. G. (1980) Modeling groundwater flow over 
sloping beds. Water Resour. Res., 16(6), 1114–1118. 
https://doi.org/10.1029/WR016i006p01114.  
Chapman, T. G., Dressler, R. F. (1984) Unsteady shallow 
groundwater flow over a curved impermeable boundary. 
Water Resour. Res. 20(10), 1427–1434. 
https://doi.org/10.1029/WR020i010p01427.  
Chapman, T. G., Ong, G. (2006) A new equation for 
shallow groundwater flow over a curved impermeable 
boundary: Numerical solutions and laboratory tests. 
Water Resour. Res., 42, W03427, 
https://doi.org/10.1029/2005WR004437. 
Childs, E. C. (1971) Drainage of groundwater resting on 
a sloping bed. Water Resour. Res., 7(5), 1256–1263. 
https://doi.org/10.1029/WR007i005p01256.  
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
108 
Dagan, G. (1967) Second-order theory of shallow free-
surface flow in porous media. Q. J. Mech. Appl. Math., 
20, 517–526. https://doi.org/10.1093/qjmam/20.4.517.  
Di Nucci, C. (2018) Unsteady free-surface flow in porous 
media: One-dimensional model equations including 
vertical effects and seepage face. C. R. Mecanique, 346, 
366–383. 
Dupuit, J. (1863) Études Théoriques et Pratiques sur le 
Mouvement des Eaux dans les Canaux Découverts et a 
Travers les Terrains Perméables (Theoretical and 
Practical Studies on Water Movement in Open Channels 
and across Permeable Terrains). deuxième édition; 
Dunod: Paris, France, 229–293 [in French]. 
Forchheimer, P. (1886) Über die ergieblgkeit von 
brunnenanlagen und sickerschlitzen (On the productivity 
of well systems and seepage slits). Zeitschnft des 
Archttekten und Ingcnieurvereins zu Hannover, 32(7), 
538–564 [in German].  
Forchheimer, P. (1914) Hydraulik (Hydraulics). 1st ed.; 
Druck und Verlag Von B. G. Teubner: Leipzig und 
Berlin, Germany, 420–467 [in German]. 
Fukuchi, T. (2018) New high-precision empirical 
methods for predicting the seepage discharges and free-
surface locations of earth dams validated by numerical 
analyses using the IFDM. Soils Found., 58, 427–445. 
Hall, H. P. (1955) An investigation of steady flow toward 
a gravity well. La Houille Blanche, 10(1-2), 8–35. 
Hantush, M. S. (1962) On the validity of the Dupuit–
Forchheimer well-discharge formula. J. Geophys. Res., 
67(6), 2417–2420. 
https://doi.org/10.1029/JZ067i006p02417.  
Harr, M. E. (1991) Groundwater and Seepage. Dover 
Publications Inc.: Mineola, NY, USA. 
Hilberts, A. G. J., von Loon, E. E., Troch, P. A., Paniconi, 
C. (2004) The hillslope-storage Boussinesq model for 
non-constant bedrock slope. J. Hydrol., 291, 160–173. 
https://doi.org/10.1016/j.jhydrol.2003.12.043.  
Kacimov, A., Obnosov, Y. (2012) Analytical solutions 
for seepage near material boundaries in dam cores: The 
Davison–Kalinin problems revisited. Appl. Math. 
Model., 36, 1286–1301. 
https://doi.org/10.1016/j.apm.2011.07.088.  
Kashef, A. I. (1969) Groundwater movement toward 
artificial cuts. Water Resour. Res., 5(5), 1032–1040. 
https://doi.org/10.1029/WR005i005p01032.  
Knight, J. H. (2005) Improving the Dupuit–Forchheimer 
groundwater free-surface approximation. Adv. Water 
Resour., 28(10), 1048–1056. 
https://doi.org/10.1016/j.advwatres.2005.04.014.  
Kozeny, J. (1931) Grundwasserbewegung bei freiem 
Spiegel, Fluss und Kanalversickerung (Groundwater 
movement at free level, River and Channel Infilitration). 
Wasserkraft und Wasserwirtsehaft, 26(3), 28–31 [in 
German]. 
Lacy, S. J., Prevost, J. H. (1987) Flow through porous 
media: A procedure for locating the free surface. Int. J. 
Numer. Anal. Meth. Geomech., 11, 585–601. 
https://doi.org/10.1002/nag.1610110605.  
Matthew, G. D. (1991) Higher order, one-dimensional 
equations of potential flow in open channels. Proc. Instn. 
Civ. Eng., 91(3), 187–201. 
https://doi.org/10.1680/iicep.1991.14974.  
Mikhailov, G. K. (1952) Concerning the Seepage in 
Trapezoidal Dams on Horizontal Bedrock. 
Gidrotekhnika i Melioratsiya, No. 1, USSR [in Russian]. 
Muleshkov, A. and Banerjee, S. (1987) Seepage towards 
vertical cuts. J. Geotech. Eng., 113(12), 1419–1431. 
https://doi.org/10.1061/(ASCE)0733-
9410(1987)113:12(1419).  
Muskat, M. (1946) The Flow of Homogeneous Fluids 
through Porous Media. J. W. Edwards, Inc.: Ann Arbor, 
MI, USA. 
Pavlovsky, N. N. (1931) On Seepage through Earth 
Dams. Instituta Gidrotekhniki i Melioratsii, Leningrad, 
USSR [in Russian]. 
Polubarinova-Kochina., P. Ya., Falkovich, S. B. (1951) 
Theory of filtration of liquids in porous media. In: 
Advances in Applied Mechanics, von Mises, R. and von 
Kármán, T. (eds.), Academic Press, New York, NY, 
USA, 2, 153–225. 
Polubarinova-Kochina, P. Ya. (1962) Theory of 
Groundwater Movement. Princeton University Press: 
Princeton, NJ, USA. 
Rayleigh, L. (1876) On waves. Phil. Mag., 1(4), 257–
279. Available online: 
http://archive.org/details/londonedinburgh511876lond 
(accessed on 10 November 2022). 
Vauclin, M., Khanji, D., Vachaud, G. (1979) 
Experimental and numerical study of a transient, two-
dimensional unsaturated-saturated water table recharge 
problem. Water Resour. Res., 15(5), 1089–1101. 
https://doi.org/10.1029/WR015i005p01089.  
Zerihun, Y. T. (2018) Extension of the Dupuit–
Forchheimer model for non-hydrostatic flows in 
Zerihun Y.: Two-dimensional unconfined seepage flow toward a highway cut slope – Dvodimenzionalni tok podzemne 
vode s prosto gladino proti avtocestnemu vkopu  
Acta hydrotechnica 36/65 (2023), 95–109, Ljubljana 
109 
unconfined aquifers. Fluids, 3(2),  
https://doi.org/10.3390/fluids3020042. 
Zerihun, Y. T. (2020) A numerical investigation of 
transient groundwater flows with a phreatic surface along 
complex hillslopes. Slovak J. Civ. Eng., 28(1), 11–19. 
https://doi.org/10.2478/sjce-2020-0002.