For hydraulic routing through coarse rockfill dams, there is still debate on whether the inertia terms might be neglected as a result of the drag force generated by the rock materials. In this study, a one-dimensional unsteady model for flow-through rockfill dams is built. For this purpose, inertia terms of Saint–Venant equations are disregarded. A semi-implicit scheme adopted for linearizing the nonlinear friction term within the time integration satisfies the Courant–Friedrich–Lewy stability criterion. The most challenging issue in the modeling of flows through rockfill dams is the appropriate definition of boundary conditions at the dam's exit zone. In addition to the analysis of different exit boundary conditions proposed in the literature, a Neumann-type boundary condition suitable for the non-inertia wave equation is also employed to estimate the exit boundary condition. This procedure is basically in appreciation of the nonlinear behavior of the water surface closer to the exit boundary. Due to the existence of the sloping edges in the trapezoidal-shaped dam, an effective length is considered for the solution domain. Finally, the model is compared with observed data and a dynamic wave model. A very good match is observed, which builds confidence in the presented modeling approach.

  • The problem of flow-through rockfill dams is considered numerically.

  • The non-inertia wave equation has been used to frame the governing equations.

  • A semi-implicit scheme is adapted for the solution.

  • Both rectangular and trapezoidal dam configurations are studied.

  • The exit boundary condition has been considered properly in the solution of the governing equations.

Graphical Abstract

Graphical Abstract
Graphical Abstract
a

nonlinear flow constant [L−1T]

b

nonlinear flow constant [L−2T2]

B

dam width [L]

a

nonlinear flow constant [L−1T]

b

nonlinear flow constant [–]

c

dimensional coefficient [L−1T2]

d

particle diameter [L]

e

void ratio [–]

Fr

Froude number [–]

g

gravity acceleration [LT²]

h

water depth [L]

h0

water depth at the start point for the analytical solution [L]

hdown

downstream water depth [L]

hexit

exit height [L]

hup

upstream water depth [L]

hx

water depth at distance x for the analytical solution [L]

i

hydraulic gradient [–]

K

empirical correction coefficient [–]

L

dam length [L]

m

hydraulic mean radius [L]

p

porosity [–]

q

discharge per unit width [L2T−1]

r

rating curve constant [–]

s

rating curve constant [–]

Q

discharge [L3T−1]

re

constant representing the particle surface area efficiency [–]

Sf

friction slope [–]

S0

bed slope [–]

t

time [T]

u

bulk velocity [LT−1]

uexit

exit velocity [LT−1]

W

Wilkins' empirical constant [L0.5T−1]

x

space [L]

α

hydraulic conductivity of the rockfill dam [LT−1]

β

constant representing the degree of turbulence [–]

ɛ

threshold value for wet/dry treatment [L]

ν

kinematic viscosity [L2T−1]

θdown

angle of the downstream face of the dam [–]

θff

angle of the emergent seepage field [–]

θup

angle of the upstream face of the dam [–]

The construction of a concrete spillway in some circumstances could be unimportant, costly and inflexible for flood control since it is designed only for a specific design storm (Townsend et al. 1991). They also obstruct the natural river flow and diminish the river's auto-purification behavior (Talaee et al. 2011). One alternative option could be the use of rocks in hydraulic structures, such as gabions and rockfill dams (Hansen et al. 1995; Chanson 1996). Rockfill dams without the impermeable core are considered one of the most effective structural methods for the rapid control of floods (Samani et al. 2013). More flow passes through the dam compared to other dam classes (Hansen & Roshanfekr 2012). A concrete layer can also cover rockfill structures on the upstream face or accompanied by an in-built spillway dam (Chanson 1996). These kinds of barriers are environmentally friendly (Kells 1996) and identified as detention dams attenuating and delaying the inflow hydrograph during floods (Findikakis & Tu 1985). In rockfill dams, flow is deviated from Darcy's law due to its coarse rock materials (Greenly & Joy 1996) and characterized by the non-Darcy regime, namely Forchheimer law, under the Dupuit's assumptions (Bari & Hansen 2002).

In terms of engineering practice, the engineers' most significant concern is the stability of rockfill structures in the vicinity of the exit point (Chanson 1996; Hansen et al. 2005). Within this region, the topmost unstable rock materials are located at the seepage flow's emergence with high flow velocity and steep exit gradient (Hansen et al. 2005; Siddiqua et al. 2011). As a consequence of improper design due to incorrect estimation of the exit height, the initiation of a few rock materials' motions causes a domino-like movement affecting other stones (Hansen & Roshanfekr 2012). One option to study the flow-through rockfill structures adjacent to the exit zone properly could be two-dimensional (2D) models under the mesh-based or mesh-free approaches to account for the non-hydrostatic pressure behavior (Townsend et al. 1991; Samani et al. 2003; Larese et al. 2015; Peng et al. 2017). However, these models are costly in computational time and suffering from some difficulties. In the 2D mesh-based methods, the mesh should be created for the whole domain even for the non-fluid part to find the zero-pressure line representing the phreatic surface's location; moreover, mesh generation is tricky on the dam's sloping faces. Although mesh-free methods can solve these two problems due to its Lagrangian nature (Peng et al. 2017), the definition of inlet and outlet boundaries is problematic in an unsteady inflow condition (Ferrand et al. 2017). Furthermore, they may experience unphysical pressure fluctuations (Jandaghian & Shakibaeinia 2020).

One-dimensional (1D) models can apply as powerful tools to study the flow-through rockfill dam that is very similar to the open-channel flow under Dupuit's assumptions (Bari & Hansen 2002), except for regions closer to the outlet that equipotential lines are substantially curved (Bear & Cheng 2010; Hansen & Roshanfekr 2012). In this case, the friction slope is not attributed solely to the roughness of the channel bed but to the properties of the porous structure of the dam as well (Bari & Hansen 2002). In flood routing through rockfill dams, previous studies have been mainly focusing on the lumped system models as an effective technique for flood control (Townsend et al. 1991; Hansen et al. 1995; Greenly & Joy 1996; Bari & Hansen 2002; Samani et al. 2003, 2004; Sarkhosh et al. 2018). However, there are a few studies in the literature regarding flood routing using a distributed system, so-called hydraulic routing. Hosseini & Joy (2007) developed an implicit one-dimensional dynamic wave formulation, ROCKFLOW, to analyze unsteady nonlinear flow-through valley fills, and rock drains as coarse porous media. They utilized a modified version of Saint–Venant equations and adopted the Newton–Raphson iteration due to the nonlinear system of equations. Ostad & Shayannejad (2015) studied the impact of the rockfill dam as a detention structure on a flood wave along a river stretch via the HEC-RAS model. Their research showed that the peak of the flood hydrograph after reaching the rockfill dam reduced significantly, and the effect of the rockfill dam as a detention structure was well observed.

Although the dynamic wave models are the most reliable methods for simulating unsteady flow in open channels, in practice, it is preferred to develop uncomplicated models for solving particular problems with fewer computational difficulties and higher efficiency (Fernández-Pato & García-Navarro 2016). Tsai (2003) reported in gradually varied unsteady flow (GVUF) in entirely or roughly prismatic channels, the inertia terms (convection and local accelerations) are usually of the same order in magnitude, but opposite in signs. Therefore, inertia terms in the momentum can be neglected, forming non-inertia wave models, since the gravity wave dominates over the inertia force (Ponce 1990; Singh 1996). One- and two-dimensional forms of these simplified models have widely applied with good accuracy comparable to the dynamic wave models in flood routing in natural channels (Hromadka & Yen 1986; Cappelaere 1997; Kazezyılmaz-Alhan & Medina 2007; López-Barrera et al. 2012; Weill et al. 2014; Fernández-Pato & García-Navarro 2016). A comparison of the explicit and implicit schemes for solving the non-inertia wave equation was thoroughly described by Caviedes-Voullième et al. (2020). The main restriction of the explicit schemes is the application of low time steps to satisfy the stability criterion. In this regard, implicit or semi-implicit methods are proposed to solve this restriction (Gross et al. 1998). Semi-implicit methods have the advantage of linearizing the nonlinear problems during the time integration without the need for an iteration scheme, such as Picard or Newton–Raphson methods so that a scalar quantity (e.g. water depth in the shallow water equations (SWE) or concentration in advection-diffusion equations) is solved implicitly and then the velocity is updated explicitly or semi-implicitly (Casulli 1990; Gross et al. 1998).

In the present work, this problem is considered in the realm of a one-dimensional non-inertia wave through rockfill dams. The main issue of solving the momentum equation is the nonlinearity arising from the friction term. In this regard, a semi-implicit numerical scheme is proposed to linearize the equation's system within the time integration. This scheme satisfies the Courant–Friedrich–Lewy (CFL) stability criterion in non-inertia wave models since the gravity wave is dominant for most shallow water applications (Giraldo & Restelli 2010). This method leads to solving linear algebraic equations using the Thomas algorithm, reducing the computational complexities arising from the iterative methods, El-Amin et al. (2011, 2016). The finite volume method (FVM) as the most effective mesh-based maintaining conservation principle locally (Eymard et al. 2000) is used for the spatial discretization. The estimation of the exit height as the boundary condition imposed in 1D mathematical simulations of the flow-through rockfill dams is the most arguable issue in the hydraulic literature of rockfill structures since the Dupuit's assumptions may not be accurately applied as discussed. Different kinds of the downstream boundary conditions presented in the literature are thoroughly investigated to improve the accuracy of the model for both steady and unsteady-state conditions. Moreover, a new strategy of the phreatic surface computation neighboring the outlet is presented based on the Neumann-type boundary condition in harmony with the non-inertia wave equation's nature. An effective length is considered for the solution domain due to the existence of the sloping boundaries in the trapezoidal-shaped dam.

One-dimensional SWE and a proper numerical scheme with associated initial and boundary conditions are presented in the following sections to fulfill the study's computational objectives. These nonlinear hyperbolic partial differential equations are derived from the incompressible Navier–Stokes equations by integrating the continuity and momentum equations over the cross-sectional area (Chow 1959). Therefore, it is called cross-section averaged SWE (Chang et al. 2011).

Governing equations

For the simulation of unsteady flows through rockfill dams, Saint–Venant equations are adopted similarly to open channel flows but considering the porosity into both the continuity and momentum equations (Hosseini & Joy 2007). These equations for a rockfill dam with a rectangular cross-section are as follows, respectively:
(1)
(2)
where h is the water depth (L), u is the bulk velocity of the flow (LT−1), p is the porosity (dimensionless), x is the space (L), t is the time (T), Sf and S0 are the friction and bed slopes, respectively (dimensionless) and g is the gravitational acceleration (LT−2). In non-inertia models, the pressure, gravity and friction force terms play the most significant roles, and so the convective and local accelerations (the third and fourth terms on the right-hand side of Equation (2)) may be neglected. Therefore, we have:
(3)

In open channels, it is customary to use empirical relationships such as the Manning or Chezy equations to account for the resistance caused, essentially, by the confining walls (Chow et al. 1988). It is conceivable to estimate the friction slope from a non-Darcy flow equation, a nonlinear relationship between the hydraulic gradient and the flow velocity. In this case, the resistance due to form drag dominates the bed friction (Hansen 1992). Numerous scientists have proposed a variety of non-Darcy coefficients, either in forms of power or quadratic polynomial functions, according to ample experimental data and particle characteristics (e.g. Ergun 1952; Wilkins 1955; McCorquodale et al. 1978; Stephenson 1979; van Gent 1996).

In 1998, Li et al. introduced a general non-Darcy relationship for seepage flow-through rockfill structures by conducting a large number of experiments with a broad spectrum of rockfill parameters. Hansen et al. (1995) evaluated experimentally various non-Darcy equations presented in the literature regarding the degree of turbulence and stated that the quadratic polynomial equation is superior to the power form due to its adaptability to the flow regime. In the present research, the quadratic polynomial form of the non-Darcy equation is adopted since the degree of turbulence might vary due to the unsteadiness and so it is compatible with the non-inertia wave model. In this regard, Forchheimer introduced an equation of motion, which takes into consideration the nonlinear effects arising from drag forces. This equation takes the form (Bear 2013):
(4)
where i is the hydraulic gradient (dimensionless), u is the bulk velocity of the flow (LT−1), a and b are coefficients that are obtained experimentally through fitting exercise. They are characteristics of rockfill structure and fluid properties. The above equation is the standard Forchheimer equation, which can be adopted for flow-through rockfill dams under steady-state conditions (Bari & Hansen 2002). On the other hand, an additional term is required to incorporate the unsteady-state effects resulting in the extended Forchheimer equation with the following structure (Hassanizadeh & Gray 1987):
(5)
where c is a constant. Hosseini & Joy (2007) examined the effect of the additional term of the Forchheimer equation (the third term on the right-hand side of Equation (5)). They concluded that it is insignificant since the acceleration of the gradually varied unsteady flow-through coarse porous media is not very fast. Bari & Hansen (2002) declared that it is reasonable to substitute the Forchheimer equation for the friction slope. Hence, in this study, the standard Forchheimer equation is included as the friction term into the momentum equation as outlined below:
(6)

Two types of rockfill dam configurations are considered in this study, namely rectangular-type and trapezoidal-type rockfill dams. In the rectangular-type rockfill dams, the computational domain is fixed in size with respect to the water surface profile. In other words, the computational area spans the whole length of the dam (L in Figure 1) so that the emergence point of the flow exactly occurs at the end of the dam length. However, in trapezoidal-type rockfill dams, this is not the case, and the downstream extent of water surface changes with time, and so the computational domain differs from the dam length, as shown in Figure 2 by L2.

Figure 1

Rectangular rockfill dam: the head values and the properties are defined in the cells, whereas the velocity and harmonic head defined at the faces.

Figure 1

Rectangular rockfill dam: the head values and the properties are defined in the cells, whereas the velocity and harmonic head defined at the faces.

Close modal
Figure 2

Trapezoidal rockfill dam: the solution domain (L2) varies with time.

Figure 2

Trapezoidal rockfill dam: the solution domain (L2) varies with time.

Close modal

Numerical scheme

Equations (1) and (6) are the fundamental equations applied to the model of the non-inertia wave. In this research, a finite volume (FV) scheme is used for the spatial discretization of the governing equations. Therefore, the discretization of momentum and continuity equations is outlined in the following:

The momentum equation

As can be seen from Equation (6), the momentum equation is nonlinear since term u2 exists in the friction term. It is possible to utilize a two successive time steps k and k + 1 for linearizing the quadratic nonlinear friction term within the time integration (Casulli 1990; Casulli & Cheng 1992; Sun et al. 2012; Salama et al. 2013, 2014, 2016; Fuhrmann et al. 2014). Therefore, the linearized momentum equation will be:
(7)

The continuity equation

In non-inertia wave models, it would be possible to combine the continuity equation with the reduced momentum equation. By doing this, the equations are converted to a single partial differential equation with only the head as the unknown variable (Hunter et al. 2005; López-Barrera et al. 2012). Replacing the velocity obtained from Equation (7) into the continuity (Equation (1)) yields:
(8)
The above equation is spatially discretized using both upwinding (for the first-order derivative), central scheme (for the second-order derivative) and temporally integrated via Backward Euler as outlined below:
(9)
where subscripts j − 1 and j refer to the successive space steps of centered cells, as shown with points in Figure 1, whereas subscripts and refer to the values at cell faces. The first-order derivative is discretized using an upwinding scheme knowing that the direction of velocity is known, and the second-order derivative is discretized using a cell-centered scheme. Water depth at cell face is calculated using the harmonic mean, while other properties are calculated using the arithmetic mean (please note that, in this work, we considered a homogeneous porous domain and so there is no difference between properties at faces and the center of the cell). The harmonic mean is calculated as:
(10)
It is worth mentioning that it is possible to rewrite Equation (10), as shown below:
(11)
where the coefficients A, B and C are defined in the following forms:
(12)
As can be seen from the system of equations, all terms are linear and so there is no need for iteration at each time step since the nonlinear algebraic equation has been transformed into a linear system of equation. After calculating water depths, the velocities are updated explicitly as follows:
(13)

Initial and boundary conditions

This problem represents a class of partial differential equations that are customarily identified as Initial-Boundary Value Problems. The initial and boundary conditions for the problem under investigation are as follows.

Initial condition: the solution for h (x, t) and u (x, t) requires knowledge of the initial conditions, i.e., h (x, 0) and u (x, 0) throughout the dam, at the beginning of the calculations. In this regard, the water depths at all nodes are considered to be approximately zero (ɛ = 0.0001) at t = 0 for both the rectangular and trapezoidal rockfill dam. It is worth noting that the threshold value (ɛ) is adopted to treat wet/dry interface (Aureli et al. 2000; Defina 2000). Consequently, the initial velocities are obtained from Equation (7) as equal to S0/a.

Upstream boundary condition: the upstream boundary condition (UBC) based on the discharge Q (xin, t) is automatically obtained considering an inflow hydrograph, and then using the upstream rating curve (i.e., hup = rQs where r and s are the rating curve constants), the UBC with respect to the hydraulic head is determined (i.e., h (xin, t) = hup).

Downstream boundary condition:h (xexit, t), defined here as hexit, represents the downstream boundary condition (DBC). Open outflow boundary condition (OOBC) is one of the most problematic boundary conditions in computational fluid dynamics (Gresho 1991). Regarding rockfill structures, many of the proposed exit boundary conditions assume the critical depth, without clear justifications (such as Stephenson 1979). Moreover, surface profiles across the dam using this approach do not seem to match that found experimentally except with the use of a correction parameter that is usually obtained through a fitting exercise with the experimental data (Bari 1999; Hansen et al. 2005). In this research, it is assumed that the flow is free at the exit zone, and so the emergent seepage is not affected by the tailwater. In this case, the DBC varies according to the dam shape, as described in the following.

  1. DBC for rectangular dam: as discussed in the Introduction, Stephenson (1979) suggested that the flow condition for the emergence point of the rockfill dam with a negligible bed slope is critical. In this case, the exit height may be calculated as:
    (14)
    where q is the discharge per unit width (L2T−1). Hansen et al. (2005) stated that the exit height of the flow resembles the depth at the brink-end of a free over-fall, which was presented by Henderson (1966) according to the occurrence of an M2 profile (or Dupuit Parabola). The formula is outlined below:
    (15)
  2. 2.
    DBC for trapezoidal dam: Hansen et al. (2005) extracted a relationship from Wilkins's non-Darcy flow equation (1955) for estimating the exit height:
    (16)
    where θdown is the angle of the downstream face of the dam as illustrated in Figure 2, W is Wilkins' empirical constant, which is 5.243 for crushed stones (L0.5T−1), and m is the hydraulic mean radius (L) which is approximately one-tenth of the particle diameter (d/10) or more precisely is determined in the following:
    (17)
    where e is the void ratio (dimensionless), re is a constant (dimensionless) representing the particle surface area efficiency that is around 1.3 for crushed stones.
Hansen et al. (2005) state that the above equation is only valid when the angle of the emergent seepage flow (effective flow angle), as shown with θff in Figure 3(a), is approximately equal to θdown. It means that for low discharges, which lead to small exit heights, the effective flow angle is much lower than the downstream dam slope and, consequently, Equation (16) is no longer applicable (Figure 3(b)). Another approach for employing the exit height as the downstream boundary condition is to utilize the similarity with the critical depth as discussed for the rectangular dam but using a correction factor:
(18)
where K is an empirical correction coefficient. Hansen et al. (2005) reported that its value ranges from 3 to 5, which is subject to the angle of the downstream dam face. In 2018, Sarkhosh et al. (2018) generalized the Pavlovsky's solution to coarse porous media using an optimization technique under the non-Darcy regime. For this purpose, they modified the hydraulic gradient to take the nonlinear effects into account in addition to the correction of the hydraulic conductivity as outlined in the following:
(19)
where hdown is the downstream water depth as depicted in Figure 2 (determined from downstream rating curve via known discharge), α is the hydraulic conductivity of the rockfill dam (LT−1), and β is a constant representing the degree of turbulence. The magnitude of these coefficients is attributed to the particle diameter assuming a uniform distribution of rockfill materials as formulated below:
(20)
Figure 3

(a) Effective flow angle at high discharge (adapted from Hansen et al. (2005)). (b) Effective flow angle at low discharge (adapted from Hansen et al. (2005)).

Figure 3

(a) Effective flow angle at high discharge (adapted from Hansen et al. (2005)). (b) Effective flow angle at low discharge (adapted from Hansen et al. (2005)).

Close modal

Neumann-type boundary condition for the exit point

Although Hosseini & Joy (2007) applied the critical-flow condition to their unsteady-state model as the DBC for dams with both vertical and sloping downstream faces, all the above equations are proposed under the steady-state condition where the M2 profile (Dupuit parabola) near the exit zone is developed. However, it is questionable whether the above equations could be suitable for the DBC under the unsteady-state condition since the M2 profile might not be formed near the outlet at some transient times. Hunt (2005) and also Castro-Orgaz & Dey (2014) found a discrepancy between the steady and unsteady curvilinear phreatic flows subjected by flood waves through fine porous media.

On the other hand, unsteady-state experiments for turbulent flow-through coarse porous media confirmed this variation in streamline curvatures due to the unsteadiness of the flow (Hosseini & Joy 2007). Therefore, in this study, a new strategy for considering the exit height as the DBC is presented. Singh (1996) suggested that for a diffusion wave modeling through open channels, the Neumann-type boundary condition is suitable for OOBC with a unique and stable solution. As a result, in addition to the analysis of the above equations as DBCs for both steady and unsteady flows, a Neumann-type BC is also adopted to estimate the exit height as described in the Supplementary Material, Appendix A, that is:
(21)

Algorithm

Equation (11) is a system of equations that should be solved implicitly at each time step using the Thomas algorithm. After water depth calculations, velocities are obtained explicitly from the momentum equation (Equation (13)). The algorithm for hydraulic calculations is demonstrated in Figure 4. The spatial discretization for the trapezoidal dam is problematic since the length of the solution domain varies with time due to the existence of the sloping boundaries. To tackle this problem, the following strategy is presented:

  • 1.

    The length of the dam (L) is divided into three segments L1, L2 and L3 as illustrated in Figure 2 

  • 2.

    At first attempt, the length of the solution domain is considered to be L2 as an effective length computed: where is the horizontal distance from the point of seepage entrance to the dam heel and L3 is the length between the dam toe and the final cell face at n +1. At the initial time, L3 is calculated as:

  • 3.

    The hydraulic calculation at a specific time step is performed and then is considered equal to .

  • 4.

    Since differs from the seepage emergence point (), an adjusting length (ΔL) is accounted: .

  • 5.

    L2 and L3 are updated: .

  • 6.

    The steps 3–5 are repeated until gets less than a desirable value (equal to ).

In the procedure described above, the mesh is generated within the effective length L2. In other words, the solution domain is restricted between upstream and downstream boundaries. This method is an advantage that leads to lower computational time compared to considering the whole length of the dam L in the trapezoidal-shaped dam.

Figure 4

Model algorithm.

Figure 4

Model algorithm.

Close modal

A computer program is developed to achieve the calculation objectives of this study. Afterwards, the accuracy of the model is evaluated. In this regard, the modal validation and application are conducted using the analytical solution, as well as experimental data for steady and unsteady-state conditions to build confidence in the modeling approach.

Steady-state condition

  • 1.
    At first, the model was compared with the analytical solution introduced by Hansen (1992) for the steady-state condition as outlined in Equation (22):
    (22)
    where h0 is the water depth at the start point (the exit height), hx is the water depth at the distance x. The other variables and parameters have already been defined. A rectangular dam with property as follows was adopted: dam length L = 10 m, dam width B = 1 m, bed slope S0 = 0.0, porosity p = 0.45, a = 1 and b = 54.44. Three constant upstream values equaling 1–3 m were considered using Δx = 0.01 m and Δt = 0.50 s. As can be seen from Figure 5, the results of the present model (non-inertia model) have a very good agreement with the analytical solution.

A comparison between the DBC computed from the Neumann BC (Equation (21)) and Equation (15) (Hansen et al. 2005) with the analytical solution (Hansen 1992) was summarized in Table 1. The exit height and velocity obtained from Equation (21) revealed a very close value to the one calculated from Equation (15) as shown in Table 1, so that the mean absolute error (MAE) arising from these two approaches is less than 2%. It means that the exit height showed similarity with the brink-end of a free over-fall and so deviated from the critical-flow condition, as Hansen et al. (2005) reported.

  • 2.

    Although the model has been primarily developed for the unsteady flow, it is conceivable to validate it using observed data under the steady-state condition. In this regard, the experimental data gathered by Hansen (1992), namely tests GF29d and GF31c (Supplementary Material, Appendix B: Table B.1 and B.2), for the trapezoidal dam, is adopted. In the model, the upstream water level goes up linearly from zero to a known value associated with the steady flow condition (0.143, 0.271 and 0.332 m for test GF29d and 0.157, 0.219 and 0.299 m for test GF31c within 10 s). It is then kept constant until the phreatic surface reaches stable status. At this stage, the water depths and the discharges are set equal to the steady-state values. The results of the comparison for the time step Δt= 0.20 s and space step Δx= 0.0005 m are depicted in Figures 6 and 7.

All equations for the trapezoidal dam, including Equations (16), (18), (19) and (21), were accounted for when considering the DBC. Regarding the critical-flow equation, it was recognized that the correction coefficient K varies from 3 to 5 not only due to the variation of the downstream face angle (as Hansen et al. (2005) previously reported) but also with the upstream water level so that K increases as hup rises. Therefore, it is not proper to consider Equation (18) for the DBC within GVUF. The other downstream boundary equations were investigated and summarized in Table 2 in the form of average MAE for tests GF29d and GF31c. Since the present model is unsteady, wherein the discharge is computed at the end of the simulation as the phreatic surface, in addition to the discharge, remains unchanged.

Figure 5

Comparison of the model with analytical solution (Hansen 1992) using the Neumann-type B.C.

Figure 5

Comparison of the model with analytical solution (Hansen 1992) using the Neumann-type B.C.

Close modal
Figure 6

Comparison of the non-inertia wave model with experimental data GF29d (Hansen 1992) using the Neumann-type B.C.

Figure 6

Comparison of the non-inertia wave model with experimental data GF29d (Hansen 1992) using the Neumann-type B.C.

Close modal
Figure 7

Comparison of the non-inertia wave model with experimental data GF31c (Hansen 1992) using the Neumann-type B.C.

Figure 7

Comparison of the non-inertia wave model with experimental data GF31c (Hansen 1992) using the Neumann-type B.C.

Close modal
Table 1

Comparison of hexit and uexit computed from Neumann B.C. and Hansen et al. (2005) equation with the analytical solution (Hansen 1992)

DBC computed from Equation (15) (Hansen et al. 2005)
DBC computed from Equation (21) (Neumann B. C.)
Analytical (Hansen 1992)
hup (m)hexit (m)MAE hexit (%)uexit (m/s)MAE uexit (%)hexit (m)MAE hexit (%)uexit (m/s)MAE uexit (%)hexit (m)uexit (m/s)
0.0374 1.5 0.4491 1.6 0.0384 1.1 0.4375 1.0 0.0380 0.4421 
0.0803 0.5 0.6588 0.5 0.0813 0.7 0.6506 0.9 0.0807 0.6555 
0.1244 1.9 0.8207 1.9 0.1254 1.1 0.8142 1.1 0.1268 0.8052 
DBC computed from Equation (15) (Hansen et al. 2005)
DBC computed from Equation (21) (Neumann B. C.)
Analytical (Hansen 1992)
hup (m)hexit (m)MAE hexit (%)uexit (m/s)MAE uexit (%)hexit (m)MAE hexit (%)uexit (m/s)MAE uexit (%)hexit (m)uexit (m/s)
0.0374 1.5 0.4491 1.6 0.0384 1.1 0.4375 1.0 0.0380 0.4421 
0.0803 0.5 0.6588 0.5 0.0813 0.7 0.6506 0.9 0.0807 0.6555 
0.1244 1.9 0.8207 1.9 0.1254 1.1 0.8142 1.1 0.1268 0.8052 
Table 2

Evaluation of the various DBC equations

TestEquationMAE of hexit (%)MAE of uexit (%)MAE of Q (%)
GF29d (θdown = 45°) Equation (16): Hansen et al. (2005)  35.9 34.5 1.9 
Equation (19): Sarkhosh et al. (2018)  10.7 8.7 3.3 
Equation (21): Neumann-type B.C. 2.4 1.3 3.4 
GF31c (θdown = 26.56°) Equation (16): Hansen et al. (2005)  9.9 6.3 6.3 
Equation (19): Sarkhosh et al. (2018)  17 24.8 12.25 
Equation (21): Neumann-type B.C. 2.7 4.9 4.9 
TestEquationMAE of hexit (%)MAE of uexit (%)MAE of Q (%)
GF29d (θdown = 45°) Equation (16): Hansen et al. (2005)  35.9 34.5 1.9 
Equation (19): Sarkhosh et al. (2018)  10.7 8.7 3.3 
Equation (21): Neumann-type B.C. 2.4 1.3 3.4 
GF31c (θdown = 26.56°) Equation (16): Hansen et al. (2005)  9.9 6.3 6.3 
Equation (19): Sarkhosh et al. (2018)  17 24.8 12.25 
Equation (21): Neumann-type B.C. 2.7 4.9 4.9 

Regarding the non-Darcy flow, the equation suggested by George & Hansen (1992) for converting the power form to a quadratic polynomial form is adopted, as discussed in Supplementary Material, Appendix C. Moreover, the accuracy of the computed phreatic surfaces by the Neumann-type B.C. compared to the observed data (Hansen 1992) for the tests above is summarized in Figure 8 in the form of root mean square error (RMSE).

  • 3.

    The 2D FEM model created by Larese et al. (2015) was taken into account to compare the accuracy of the present model with a 2D model. They solved two-dimensional NS equations using the Galerkin method to model free surface flow-through porous media by a semi-explicit stabilized scheme. They applied two experimental setups to validate their model for a trapezoidal-shaped rockfill dam with properties summarized in Table B.3 (Supplementary Material, Appendix B). The results of the comparison for the steady-state condition considering the Neumann B.C. for the exit height have been shown in Figure 9(a) and 9(b) and also Table 3. The time step is Δt = 0.10 s, space step Δx = 0.04 m and for computing the non-Darcy coefficients' calculations (a and b in Equation (4)), the relationship proposed by Li et al. (1998) was applied. As can be seen, the results of the present model have good agreement with the 2D FEM model as well as the experimental data.

Figure 8

The error of phreatic surfaces computed from non-inertia wave model considering the Neumann-type B.C. compared to observed data (Hansen 1992) for tests GF29d and GF31c.

Figure 8

The error of phreatic surfaces computed from non-inertia wave model considering the Neumann-type B.C. compared to observed data (Hansen 1992) for tests GF29d and GF31c.

Close modal
Figure 9

(a) Case I: Comparison of the non-inertia wave model using the Neumann-type B.C with experimental data and 2D FEM model (Larese et al. 2015). (b) Case II: Comparison of the non-inertia wave model using the Neumann-type B.C. with experimental data and 2D FEM model (Larese et al. 2015).

Figure 9

(a) Case I: Comparison of the non-inertia wave model using the Neumann-type B.C with experimental data and 2D FEM model (Larese et al. 2015). (b) Case II: Comparison of the non-inertia wave model using the Neumann-type B.C. with experimental data and 2D FEM model (Larese et al. 2015).

Close modal
Table 3

Comparison of the non-inertia wave model and 2D FEM (Larese et al. 2015) with observed data

TestModelMAE of Q (%)Phreatic surface RMSE (m)Velocity RMSE (m/s)
Case I (θdown = 18.44°) 2D FEM (Larese et al. 20150.0090 0.0008 
non-inertia wave model 2.0 0.0095 0.0014 
Case II (θdown = 24.44°) 2D FEM (Larese et al. 20150.0104 0.0010 
non-inertia wave model 5.8 0.0158 0.0013 
TestModelMAE of Q (%)Phreatic surface RMSE (m)Velocity RMSE (m/s)
Case I (θdown = 18.44°) 2D FEM (Larese et al. 20150.0090 0.0008 
non-inertia wave model 2.0 0.0095 0.0014 
Case II (θdown = 24.44°) 2D FEM (Larese et al. 20150.0104 0.0010 
non-inertia wave model 5.8 0.0158 0.0013 

Unsteady-state condition

  1. For unsteady-state conditions, first, the model was compared to the smoothed-particle hydrodynamics (SPH) model, a mesh-free Lagrangian method suggested by Peng et al. (2017). They adopted two-dimensional Navier–Stokes equations that originated from SPH multiphase flow-through coarse porous media so that they considered the solid and fluid phases. They reported that since the porous media is immovable, the equations are applied only to the fluid phase. It is worth mentioning that they neglected the interaction between the solid and fluid phases. The behavior of coarse porous media was taken into account by adding a nonlinear resistance force (Forchheimer equation) as an external force into the momentum equation. They simulated the unsteady flow for a trapezoidal dam with properties mentioned in Supplementary Material, Table B.4.

The comparison between the non-inertia model and SPH model, using Δt = 0.20 s and Δx = 0.0005 m, has been depicted in Figures 10 and 11. Regarding the non-Darcy coefficients' calculations, the relationship suggested by Li et al. (1998) (Supplementary Material, Appendix C) was employed in the non-inertia wave model. Peng et al. (2017) state that at t = 100 s, the pressure stays unchanged, which means the flow becomes steady with discharge equaling to 25.46 L/s. The flow discharge from the FVM model was computed 25.22 L/s at this stationary stage showing that the discharge error is less than 1%. Furthermore, the FVM and SPH velocity profiles are also shown in Figure 10 in the steady-state condition (t = 100 s). The velocity error computed from the present model compared to the SPH model is RMSE = 0.0019 m/s.

  • 2.

    The model was also assessed for the unsteady-state condition using test SU4 conducted by Hosseini & Joy (2007). For this purpose, a rectangular dam with characteristics mentioned in Supplementary Material, Appendix B: Table B.5 was adopted. The UBC is in the form of hup variation with time, as displayed in Figure 12. The overall time simulation is 250 s. Following the collected data from experiments, the results of the comparison for Δx = 0.01 m, and Δt = 0.15 s at three times have been illustrated in Figure 13. Moreover, the results from the dynamic wave model titled: ROCKFLOW, introduced by Hosseini & Joy (2007), were employed in the comparison.

The quantitative errors of the phreatic surfaces computed from both models are shown in Figure 14. As can be seen from Figures 13 and 14, the results show more variability at an earlier time (t= 26 s) when the flow is highly transient so that the non-inertia model overestimates the phreatic surface compared to the observed data, whereas there is an underestimation of the water surface profile in the dynamic wave model. However, the error arising from the two models is close and approximately 1 cm. Furthermore, the present model shows less error than the dynamic wave model (Hosseini & Joy 2007) compared to the experimental data at t = 46 s. It might be within the interpolation process of the discretely observed data for imposing as the inflow hydrograph. There is a slight discrepancy between the observed inflow hydrograph with one adopted as the UBC in the dynamic model proposed by Hosseini & Joy (2007) in the interval between t = 40 and 100 s. In the present model, high-order polynomial interpolation was adopted for obtaining the inflow hydrograph as the UBC. The results of both models show the highest accuracy at t = 106 s with less than 0.004 m error.

Figure 10

Comparison of the FVM non-inertia wave model using the Neumann-type B.C. with the SPH model (Peng et al. 2017).

Figure 10

Comparison of the FVM non-inertia wave model using the Neumann-type B.C. with the SPH model (Peng et al. 2017).

Close modal
Figure 11

The error of phreatic surfaces computed from non-inertia wave model using the Neumann-type B.C. compared to the SPH model (Peng et al. 2017).

Figure 11

The error of phreatic surfaces computed from non-inertia wave model using the Neumann-type B.C. compared to the SPH model (Peng et al. 2017).

Close modal
Figure 12

The UBC for gradually varied unsteady flow (GVUF) obtained from the experimental data (Hosseini & Joy 2007).

Figure 12

The UBC for gradually varied unsteady flow (GVUF) obtained from the experimental data (Hosseini & Joy 2007).

Close modal
Figure 13

Comparison of the non-inertia with dynamic wave model and experimental data (Hosseini & Joy 2007) within GVUF.

Figure 13

Comparison of the non-inertia with dynamic wave model and experimental data (Hosseini & Joy 2007) within GVUF.

Close modal
Figure 14

Error arising from the non-inertia and dynamic wave models compared to the experimental data (Hosseini & Joy 2007) within GVUF.

Figure 14

Error arising from the non-inertia and dynamic wave models compared to the experimental data (Hosseini & Joy 2007) within GVUF.

Close modal

From Figures 13 and 14, it can be inferred that the inertia terms, in general, do not meaningfully affect the prediction of the flow hydraulics for flow-through rockfill dam. It is due to the impact of the high magnitude of the form drag, which is exerted by the coarse porous media, on the local and convective accelerations. The results show more variability between the two models at earlier times when the flow is more transient. Moreover, since the momentum equation is reduced to a non-inertia wave equation, there is no need for computing the velocities at each time step implicitly that is common in the implicit scheme of a dynamic wave model.

  • 3.

    Up to now, all validation stages were conducted in the scope of GVUF. In open-channel flows, the flow is highly convective in rapidly varied unsteady flow (RVUF), which means the gravity force is not dominant over the inertia force. Thus, it is not proper to employ the non-inertia wave model. However, the flow might maintain the non-inertia nature as it interacts with a rockfill structure. In this part, the non-inertia wave model for the flow-through porous structure is assessed within RVUF. Liu et al. (1999) conducted a series of experiments for rapid flow in the presence of a porous structure that has been a benchmark in recent years (Larese et al. 2015; Peng et al. 2017; Xu & Jin 2019; Sun et al. 2019). In this regard, one of their experimental setups mentioned in Table B.6 (Supplementary Material, Appendix B) is presented.

A porous structure is built in the middle of a flume enclosed by walls, and a closed gate is accommodated 2 cm upstream of the porous stricture. At t = 0, the gate is removed immediately, and the wave propagates through the porous structure. Herein, the porous structure's exit region is not a free-flow condition due to the existence of the enclosed wall downstream of the flume. Therefore, the UBC and DBC for both entrance and exit points of the dam should be first obtained from the observed data (Liu et al. 1999) as depicted in Figure 15. It is worth noting that since the observed data at the boundaries are sinusoidal, higher-order polynomial interpolation was applied instead of a power form.

Figure 15

The UBC and DBC for rapidly varied unsteady flow (RVUF) obtained from the experimental data (Liu et al. 1999).

Figure 15

The UBC and DBC for rapidly varied unsteady flow (RVUF) obtained from the experimental data (Liu et al. 1999).

Close modal

Furthermore, unlike the previous comparisons that the equation proposed by Li et al. (1998) was adopted, herein, the non-Darcy flow equation generated by van Gent (1996) (Supplementary Material, Appendix C) is employed. The main reason is that due to the high magnitude of turbulence, the flow regime in RVUF differs from other relationships in the literature, as reported by Liu et al. (1999). The comparison of the present model with the experimental data for Δx = 0.0005 m, and Δt = 0.05 s have been gathered in Figure 16. The total simulation is 2 s, and 10 subfigures are considered. Furthermore, the Froude number (Fr = u/p (g h)0.5) is simulated.

Figure 16

Comparison of the non-inertia with dynamic wave model and experimental (Liu et al. 1999) data within RVUF.

Figure 16

Comparison of the non-inertia with dynamic wave model and experimental (Liu et al. 1999) data within RVUF.

Close modal

Figure 16 shows that the flow behavior differs from the ordinary dam-break problems in open channels due to the high magnitude of the drag force generated by the rockfill materials. The Froude number remains in the vicinity of the subcritical flow condition that guarantees the precondition of the non-inertia wave model. Except for t = 0.2 s, the non-inertia model can predict the rapid wave propagation through the porous structure with a good agreement with the observed data. It is worth mentioning that in both 2D models proposed by Larese et al. (2015) and Peng et al. (2017), under the Eulerian and Lagrangian frameworks, respectively, the discrepancy of their model results with the observed data at t= 0.2 s was reported.

It is of great importance to evaluate the model response to various spatial and temporal steps to ensure the model accuracy, detect the potential restrictions of the model applicability and choose appropriate spatial and temporal intervals. In this part, the model performance is evaluated regarding various space and time steps. First, the Neumann-type boundary condition's sensitivity to various space steps and the mean size of the rockfill particles are analyzed. Then, the numerical stability concerning the simulation time step is assessed.

Sensitivity analysis of the Neumann-type boundary condition

As can be seen from Figure 17, it is evident that the estimation of the exit height using the Neumann-type B.C. requires that Δx would be small. In other words, the magnitude of becomes close to θff provided that Δx is set up tiny. For the tests GF29d and GF31c, various numbers of cells were adopted to investigate the effect of the space step on the exit height obtained from the Neumann-type B.C., as shown in Figure 18(a) and 18(b). The dimensionless exit height versus space step is utilized so that the simulated and observed exit heights are indicated by hN.B.C and hobserved, respectively. Furthermore, the ratio of Δx to the effective length (L2 in Figure 2) is adopted. The sensitivity of the exit height to the space step for the rectangular-type rockfill dam is also assessed. In this regard, a dam with properties as follows is considered: UBC = 0.5, dam's length L = 3 m, dam width B = 1 m, bed slope S0 = 0.0, porosity p = 0.45. Three mean particle sizes D = 0.025, 0.050 and 0.075 are adopted. The analytical solution (Hansen 1992) was applied to compare the present model's results. The discharges and exit heights were obtained 0.0123 m3/s, 0.0925 m, 0.0183 m3/s, 0.1019 m, 0.0227 m3/s and 0.10495 m for the particle sizes described above, respectively. The result of the comparison is shown in Figure 19.

Figure 17

Exit height approximation using Neumann-type BC.

Figure 17

Exit height approximation using Neumann-type BC.

Close modal
Figure 18

(a) The sensitivity of the exit height obtained from the Neumann-type B.C. to the space step for test GF29d (Hansen 1992). (b) The sensitivity of the exit height obtained from the Neumann-type B.C. to the space step for test GF31c (Hansen 1992).

Figure 18

(a) The sensitivity of the exit height obtained from the Neumann-type B.C. to the space step for test GF29d (Hansen 1992). (b) The sensitivity of the exit height obtained from the Neumann-type B.C. to the space step for test GF31c (Hansen 1992).

Close modal
Figure 19

The sensitivity of the exit height obtained from the Neumann-type B.C. to the space step for a rectangular-type dam using the analytical solution (Hansen 1992).

Figure 19

The sensitivity of the exit height obtained from the Neumann-type B.C. to the space step for a rectangular-type dam using the analytical solution (Hansen 1992).

Close modal

As can be seen from Figures 18(a), 18(b) and 19, the simulated exit height converges quickly by decreasing the space step. Furthermore, as the discharge increases (for higher values of hup), the sensitivity of the exit height to the space step decreases since the angle of the emergent seepage field (θff) approaches the downstream dam's sloping face (θdown) as Hansen et al. (2005) reported. Also, the model shows no significant sensitivity to the particle size for a dam with similar properties and UBC. It is worth mentioning that the sensitivity of the model to the space step is high in the steady-state condition as small Δx is required to predict the M2 profile; however, this is not the case in transient condition.

Evaluation of the simulation time step

The numerical solution might be distorted in huge time steps under the implicit framework (Hicks 1997; Fernández-Pato & García-Navarro 2016). The shallow water wave models necessitate a proper numerical scheme without numerical damping and phase shift. The phase shift is a delay in wave propagation, and numerical damping is nonselective numerical diffusion arising from fully implicit or semi-implicit schemes in case large time steps are chosen (Imamura & Goto 1988; Casulli & Cattani 1994; Hicks 1997, Fan et al. 2020). To fulfill this criterion, it is not recommended to select very high values of the Courant in order to avoid the numerical damping and phase shift (Hicks 1997; Litrico & Fromion 2004).

In 1992, Casulli & Cheng developed a semi-implicit scheme for solving three-dimensional shallow water flow in open channels using the Crank–Nicolson method. They stated that the stability of their model is attributed to the horizontal viscous term. In terms of the rockfill structures, the viscous term (first term of the non-Darcy flow) is much larger than the horizontal viscous force (Liu et al. 1999). Therefore, the nonlinear constitutive non-Darcy equation is a limiter on the numerical stability.

Various CFL values were suggested under semi-implicit and fully implicit schemes in the literature according to the models' stabilized performance to avoid numerical damping and phase shift. Casulli (1990) reported the values up to 8 in their semi-implicit solving of 2D SWE. Gross et al. (1998) analyzed the numerical stability of their advection–diffusion equation thoroughly under the semi-implicit method and stated in the range of 1 < CFL < 10; their model performs adequately. Morse (1991) suggests a Courant value between 5 and 15 for sediment transport modeling in multi-channel rivers under the unsteady-state condition. Fernández-Pato & García-Navarro (2016) compared the explicit and implicit schemes of the non-inertia wave models for overland flow. They reported that the time step of the implicit scheme could be 11 times larger than the maximum allowed in the explicit scheme, which means the CFL for the implicit scheme could be around 11.

Regarding the simulation of the present model's time step, the model performance using various CFL was evaluated. For this purpose, the unsteady flow through the rectangular-type rockfill dam was considered to see the effect of the time simulation. The experimental data suggested by Hosseini & Joy (2007) and Liu et al. (1999), as discussed earlier, were adopted. First, the model was compared to the experimental data suggested by Hosseini & Joy (2007) at t = 26 s, in which the flow is highly transient, using two thresholds ε = 0.0001 and 0.001 m for treating the dry bed. The wave propagation results have been illustrated in Figure 20 considering Δx = 0.01 m and CFL = Δt (gh)0.5x. As can be seen from Figure 20, the numerical damping and phase shift are significant at CFL = 73.0 for the threshold value equal to 0.0001 m. The effect of wave damping is still observed at average values of the Courant (24.3 and 35.9); however, with decreasing CFL, the model results show better agreement with observed data. At CFL = 10.4, there is no numerical error. With the threshold = 0.001 m, the model shows less sensitivity so that from CFL = 36.4 down, the results do not have variability. It is in harmony with findings reported by Bellos & Sakkas (1987), Bradford & Sanders (2002) and Duran (2015), stating that the threshold value has a high effect on numerical damping and phase shift in frictional flow with the wet/dry interface.

Figure 20

Comparison of the non-inertia wave model with experimental data (Hosseini & Joy 2007) for various CFL values at t = 26 s using two different thresholds: LHS: 0.0001 m, RHS: 0.001 m.

Figure 20

Comparison of the non-inertia wave model with experimental data (Hosseini & Joy 2007) for various CFL values at t = 26 s using two different thresholds: LHS: 0.0001 m, RHS: 0.001 m.

Close modal

With respect to the second test case, the model was compared to the observed data collected by Liu et al. (1999) at times 0.2 and 2.0, representing the highly transient and quasi-stationary flows, respectively. Figure 21 shows that at t = 0.2 s, the model shows low sensitivity to the , whereas, for the quasi-stationary flow at t = 2.0 s, there is no considerable sensitivity to the Courant number in a way that it could be up to 111. On the other hand, a comparison with the highly transient flow regime between experimental data collected by Hosseini & Joy (2007) and Liu et al. (1999) indicates that lower magnitudes of CFL should be taken into account for a dry bed to avoid numerical damping and phase shift. In other words, the model is more sensitive to a dry bed than a wet one considering similar CFL values.

Figure 21

Comparison of the non-inertia wave model with experimental data (Liu et al. 1999) for various CFL values at two different times: LHS: t = 0.2 s, RHS: 2.0 s.

Figure 21

Comparison of the non-inertia wave model with experimental data (Liu et al. 1999) for various CFL values at two different times: LHS: t = 0.2 s, RHS: 2.0 s.

Close modal

To sum up, the CFL value around 10, as recommended in the literature for semi- or fully implicit schemes, satisfies the numerical stability within different flow regimes for both dry and wet beds in the present model. In the next step, the sensitivity of the preferred CFL is investigated with respect to various thresholds. In this regard, the Hosseini & Joy (2007) test is employed using various thresholds, as depicted in Figure 22. Figure 22 shows the model does not experience any numerical damping and phase shift in case the CFL equals to 10, even considering very small magnitudes of the threshold.

Figure 22

The sensitivity of the non-inertia wave model for various threshold values at wet/dry interface considering CFL = 10.4 and using experimental data (Hosseini & Joy 2007).

Figure 22

The sensitivity of the non-inertia wave model for various threshold values at wet/dry interface considering CFL = 10.4 and using experimental data (Hosseini & Joy 2007).

Close modal

This study's one-dimensional non-inertia wave model fulfills the prediction of unsteady flow-through rockfill dams with various sloping faces. Herein, the flow regime deviates from Darcy's law, and so nonlinear flow equations, so-called Forchheimer's, should be adopted. The flood routing was conducted spatially and temporally through the dam to understand the behavior of rockfill barriers as a detention dam during floods and accurately determine the hydraulic flow characteristics using the proposed model. Moreover, the model can predict the phreatic seepage under the steady-state condition with very good precision and estimate the discharge with a known upstream water level with much lower computational time than two-dimensional models.

Considering the non-inertia wave formulation seems to satisfy the accuracy of the hydraulic routing through rockfill dams. It is unnecessary to account for the whole terms of a dynamic wave that complicates the calculational procedure. This matter is mainly due to the impact of the form drag, created by the coarse material of rockfill embankments on the inertia terms. The numerical scheme presented in this paper reduces the computational complexity arising from iterative methods conventional in fully implicit models. The nonlinear friction term is linearized using a semi-implicit scheme within the time integration. By doing this linearization, the model also does not experience any unphysical oscillation. The Thomas algorithm then solves the system of the equations. As the momentum equation is reduced to a non-inertia wave equation, there is no need for computing the velocities at each step implicitly, which is common in the implicit scheme of dynamic wave models.

A reliable estimation of the exit height is of great importance for dam safety and stability. The previous DBC equations were analyzed. For rectangular-type cases, it was recognized that the exit height is more grounded on the brink-end of a free over-fall, as discussed by Hansen et al. (2005). Among the downstream boundary equations for trapezoidal ones, the critical-flow equation cannot be employed for unsteady flows. The correction factor varies not only with the angle of the downstream dam face but also with the hydraulic head in the exit zone. On the other hand, further studies seem to be required to understand the changes of this correction factor as a function of the downstream face slope and upstream water level to apply the critical-flow condition to steady-state models.

Moreover, unlike the equation suggested by Sarkhosh et al. (2018) that yields better accuracy for rockfill dams with a steep downstream sloping face, the Hansen et al. (2005) equation is more suitable for dams with a mild downstream sloping face that becomes close to the angle of the emergent flow field. Finally, a new downstream boundary condition based on the emergent seepage field was adopted into the model. It was found that the new approach improves the estimation of the exit height significantly for both rectangular and trapezoidal rockfill dams under steady- and unsteady-state conditions.

In the present model, the solution domain is restricted to an effective length that is proportional to the horizontal length of the phreatic surface. That means the mesh is generated exclusively in the vicinity restricted between the upstream and downstream boundaries, which is a benefit for reducing computational time. It is also recommended to research the use of a non-uniform locally refined mesh in the exit area that could increase the resolution in this zone and speed up the computational performance.

All relevant data are included in the paper or its Supplementary Information.

Aureli
F.
Mignosa
P.
Tomirotti
M.
2000
Numerical simulation and experimental verification of dam-break flows with shocks
.
Journal of Hydraulic Research
38
(
3
),
197
206
.
Bari
M. R.
1999
The Hydraulics of Buried Streams
.
Technical University of Nova Scotia
,
Halifax
.
Bari
R.
Hansen
D.
2002
Application of gradually-varied flow algorithms to simulate buried streams
.
Journal of Hydraulic Research
40
(
6
),
673
683
.
Bear
J.
2013
Dynamics of Fluids in Porous Media
.
Courier Corporation
,
New York
.
Bear
J.
Cheng
A. H. D.
2010
Modeling Groundwater Flow and Contaminant Transport
, Vol.
23
.
Springer Science & Business Media
,
New York
.
Bellos
C. V.
Sakkas
J. G.
1987
1-D dam-break flood-wave propagation on dry bed
.
Journal of Hydraulic Engineering
113
(
12
),
1510
1524
.
Bradford
S. F.
Sanders
B. F.
2002
Finite-volume model for shallow-water flooding of arbitrary topography
.
Journal of Hydraulic Engineering
128
(
3
),
289
298
.
Cappelaere
B.
1997
Accurate diffusive wave routing
.
Journal of Hydraulic Engineering
123
(
3
),
174
181
.
Castro-Orgaz
O.
Dey
S.
2014
Second-order shallow-flow theory and Dupuit approximation for phreatic aquifers
.
Journal of Hydraulic Engineering
140
(
9
),
04014040
.
Casulli
V.
Cattani
E.
1994
Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow
.
Computers & Mathematics with Applications
27
(
4
),
99
112
.
Casulli
V.
Cheng
R. T.
1992
Semi-implicit finite difference methods for three-dimensional shallow water flow
.
International Journal for Numerical Methods in Fluids
15
(
6
),
629
648
.
Caviedes-Voullième
D.
Fernández-Pato
J.
Hinz
C.
2020
Performance assessment of 2D zero-inertia and shallow water models for simulating rainfall-runoff processes
.
Journal of Hydrology
584
,
124663
.
Chang
T. J.
Kao
H. M.
Chang
K. H.
Hsu
M. H.
2011
Numerical simulation of shallow-water dam break flows in open channels using smoothed particle hydrodynamics
.
Journal of Hydrology
408
(
1–2
),
78
90
.
Chow
V. T.
1959
Open-Channel Hydraulics
.
McGraw-Hill
,
New York
.
Chow
V. T.
Maidment
D. R.
Mays
L. W.
1988
Applied Hydrology
.
McGraw-Hill
,
New York
.
Defina
A.
2000
Two-dimensional shallow flow equations for partially dry areas
.
Water Resources Research
36
(
11
),
3251
3264
.
El-Amin
M. F.
Salama
A.
Sun
S.
2011
Solute Transport With Chemical Reaction in Single and Multi-Phase Flow in Porous Media
. In:
Mass Transfer in Multiphase Systems and its Applications
.
doi:10.5772/594
.
El-Amin
M. F.
Kou
J.
Sun
S.
Salama
A.
2016
An iterative implicit scheme for nanoparticles transport with two-phase flow in porous media
.
Procedia Computer Science
80
,
1344
1353
.
Ergun
S.
1952
Fluid flow through packed columns
.
Chemical Engineering Progress
48
(
2
),
89
95
.
Eymard
R.
Gallouët
T.
Herbin
R.
2000
Finite volume methods
.
Handbook of Numerical Analysis
7
,
713
1018
.
Fernández-Pato
J.
García-Navarro
P.
2016
2D zero-inertia model for solution of overland flow problems in flexible meshes
.
Journal of Hydrologic Engineering
21
(
11
),
04016038
.
Ferrand
M.
Joly
A.
Kassiotis
C.
Violeau
D.
Leroy
A.
Morel
F. X.
Rogers
B. D.
2017
Unsteady open boundaries for SPH using semi-analytical conditions and Riemann solver in 2D
.
Computer Physics Communications
210
,
29
44
.
Findikakis
A. N.
Tu
S. W.
1985
Flood routing through rock dumps in natural channels
. In:
Hydraulics and Hydrology in the Small Computer Age
.
ASCE
,
San Francisco
, pp.
1213
1218
.
Fuhrmann
J.
Ohlberger
M.
Rohde
C.
2014
Finite Volumes for Complex Applications VII: Methods and Theoretical Aspects
.
Springer
,
Berlin
.
George
G. H.
Hansen
D.
1992
Conversion between quadratic and power law for non-Darcy flow
.
Journal of Hydraulic Engineering
118
(
5
),
792
797
.
Giraldo
F. X.
Restelli
M.
2010
High-order semi-implicit time-integrators for a triangular discontinuous Galerkin oceanic shallow water model
.
International Journal for Numerical Methods in Fluids
63
(
9
),
1077
1102
.
Greenly
B. T.
Joy
D. M.
1996
One-dimensional finite-element model for high flow velocities in porous media
.
Journal of Geotechnical Engineering
122
(
10
),
789
796
.
Gresho
P. M.
1991
Some current CFD issues relevant to the incompressible Navier-Stokes equations
.
Computer Methods in Applied Mechanics and Engineering
87
(
2–3
),
201
252
.
Gross
E. S.
Casulli
V.
Bonaventura
L.
Koseff
J. R.
1998
A semi-implicit method for vertical transport in multidimensional models
.
International Journal for Numerical Methods in Fluids
28
(
1
),
157
186
.
Hansen
D.
1992
The Behavior of Flow Through Rockfill Dams
.
Doctoral dissertation
,
University of Ottawa
,
Canada
.
Hansen
D.
Zhao
W. Z.
Han
S. Y.
2005
Hydraulic performance and stability of coarse rockfill deposits
.
Proceedings of the Institution of Civil Engineers-Water Management
158
(
4
),
163
175
.
Thomas Telford Ltd
.
Hassanizadeh
S. M.
Gray
W. G.
1987
High velocity flow in porous media
.
Transport in Porous Media
2
(
6
),
521
531
.
Henderson
F. M.
1966
Open Channel Flow
.
Macmillan
,
New York
.
Hosseini
S. M.
Joy
D. M.
2007
Development of an unsteady model for flow through coarse heterogeneous porous media applicable to valley fills
.
International Journal of River Basin Management
5
(
4
),
253
265
.
Hromadka
T. V.
II
Yen
C. C.
1986
A diffusion hydrodynamic model (DHM)
.
Advances in Water Resources
9
(
3
),
118
170
.
Hunt
B.
2005
Bank-storage problem and the Dupuit approximation
.
Journal of Hydrologic Engineering
10
(
2
),
118
124
.
Hunter
N. M.
Horritt
M. S.
Bates
P. D.
Wilson
M. D.
Werner
M. G.
2005
An adaptive time step solution for raster-based storage cell modelling of floodplain inundation
.
Advances in Water Resources
28
(
9
),
975
991
.
Imamura
F.
Goto
C.
1988
Truncation error in numerical tsunami simulation by the finite difference method
.
Coastal Engineering in Japan
31
(
2
),
245
263
.
Jandaghian
M.
Shakibaeinia
A.
2020
An enhanced weakly-compressible MPS method for free-surface flows
.
Computer Methods in Applied Mechanics and Engineering
360
,
112771
.
Kazezyılmaz-Alhan
C. M.
Medina
M. A.
Jr
2007
Kinematic and diffusion waves: analytical and numerical solutions to overland and channel flow
.
Journal of Hydraulic Engineering
133
(
2
),
217
228
.
Kells
J. A.
1996
The Analysis of Flow Through and Over A Gabion dam
.
Doctoral dissertation
,
University of Saskatchewan
,
Canada
.
Larese
A.
Rossi
R.
Oñate
E.
2015
Finite element modeling of free surface flow in variable porosity media
.
Archives of Computational Methods in Engineering
22
(
4
),
637
653
.
Li
B.
Garga
V. K.
Davies
M. H.
1998
Relationships for non-Darcy flow in rockfill
.
Journal of Hydraulic Engineering
124
(
2
),
206
212
.
Litrico
X.
Fromion
V.
2004
Frequency modeling of open-channel flow
.
Journal of Hydraulic Engineering
130
(
8
),
806
815
.
Liu
P. L. F.
Lin
P.
Chang
K. A.
Sakakiyama
T.
1999
Numerical modeling of wave interaction with porous structures
.
Journal of Waterway, Port, Coastal, and Ocean Engineering
125
(
6
),
322
330
.
López-Barrera
D.
García-Navarro
P.
Brufau
P.
Burguete
J.
2012
Diffusive-wave based hydrologic-hydraulic model with sediment transport. I: model development
.
Journal of Hydrologic Engineering
17
(
10
),
1093
1104
.
McCorquodale
J. A.
Hannoura
A. A. A.
Sam Nasser
M.
1978
Hydraulic conductivity of rockfill
.
Journal of Hydraulic Research
16
(
2
),
123
137
.
Morse
B.
1991
Mathematical Modelling of Sediment Transport and bed Transients in Multichannel River Networks Under Conditions of Unsteady Flow
.
University of Ottawa
,
Canada
.
Ostad-Ali-Askari
K.
Shayannejad
M.
2015
Usage of rockfill dams in the HEC-RAS software for the purpose of controlling floods
.
American Journal of Fluid Dynamics
5
(
1
),
23
29
.
Peng
C.
Xu
G.
Wu
W.
Yu
H. S.
Wang
C.
2017
Multiphase SPH modeling of free surface flow in porous media with variable porosity
.
Computers and Geotechnics
81
,
239
248
.
Ponce
V. M.
1990
Generalized diffusion wave equation with inertial effects
.
Water Resources Research
26
(
5
),
1099
1101
.
Samani
H. M.
Samani
J. M.
Shaiannejad
M.
2003
Reservoir routing using steady and unsteady flow through rockfill dams
.
Journal of Hydraulic Engineering
129
(
6
),
448
454
.
Samani
J. M. V.
Samani
H. M. V.
Shaiannejad
M.
2004
Reservoir routing with outflow through rockfill dams
.
Journal of Hydraulic Research
42
(
4
),
435
439
.
Samani
J. M.
Moknatian
M.
Heidari
M.
2013
Two-dimensional multi-rockfill detention dam flow model in reservoir routing
.
Proceedings of the Institution of Civil Engineers-Engineering and Computational Mechanics
166
(
2
),
59
67
.
Sarkhosh
P.
Samani
J. M. V.
Mazaheri
M.
2018
A one-dimensional flood routing model for rockfill dams considering exit height
.
Proceedings of the Institution of Civil Engineers-Water Management
171
(
1
),
42
51
.
Thomas Telford Ltd
.
Siddiqua
S.
Blatz
J. A.
Privat
N. C.
2011
Evaluating turbulent flow in large rockfill
.
Journal of Hydraulic Engineering
137
(
11
),
1462
1469
.
Singh
V. P.
1996
Kinematic Wave Modeling in Water Resources
.
Surface-Water Hydrology. Wiley
,
Chichester
, p.
1399
.
Stephenson
D. J.
1979
Rockfill in Hydraulic Engineering
, Vol.
27
.
Elsevier
,
Amsterdam
.
Sun
X.
Sun
M.
Takabatake
K.
Pain
C. C.
Sakai
M.
2019
Numerical simulation of free surface fluid flows through porous media by using the explicit MPS method
.
Transport in Porous Media
127
(
1
),
7
33
.
Talaee
P. H.
Heydari
M.
Fathi
P.
Marofi
S.
Tabari
H.
2011
Numerical model and computational intelligence approaches for estimating flow through rockfill dam
.
Journal of Hydrologic Engineering
17
(
4
),
528
536
.
Townsend
R. D.
Garga
V. K.
Hansen
D.
1991
Finite difference modelling of the variation in piezometric head within a rockfill embankment
.
Canadian Journal of Civil Engineering
18
(
2
),
254
263
.
van Gent
M. R. A.
1996
Wave interaction with permeable coastal structures
.
In International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts
6
(
33
),
277A
.
Weill
S.
di Chiara-Roupert
R.
Ackerer
P.
2014
Accuracy and efficiency of time integration methods for 1D diffusive wave equation
.
Computational Geosciences
18
(
5
),
697
709
.
Wilkins
J. K.
1955
Flow of water through rock fill and its application to the design of dams
.
New Zealand Engineering
10
(
11
),
382
.

Supplementary data