## Abstract

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.

## HIGHLIGHTS

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

## NOTATION

*a*nonlinear flow constant [L

^{−1}T]*b*nonlinear flow constant [L

^{−2}T^{2}]*B*dam width [L]

*a*′nonlinear flow constant [L

^{−1}T]*b*′nonlinear flow constant [–]

*c*dimensional coefficient [L

^{−1}T^{2}]*d*particle diameter [L]

*e*void ratio [–]

- Fr
Froude number [–]

*g*gravity acceleration [LT

^{−}²]*h*water depth [L]

*h*_{0}water depth at the start point for the analytical solution [L]

*h*_{down}downstream water depth [L]

*h*_{exit}exit height [L]

*h*_{up}upstream water depth [L]

*h*_{x}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 [L

^{2}T^{−1}]*r*rating curve constant [–]

*s*rating curve constant [–]

*Q*discharge [L

^{3}T^{−1}]*r*_{e}constant representing the particle surface area efficiency [–]

*S*_{f}friction slope [–]

*S*_{0}bed slope [–]

*t*time [T]

*u*bulk velocity [LT

^{−1}]*u*_{exit}exit velocity [LT

^{−1}]*W*Wilkins' empirical constant [L

^{0.5}T^{−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 [L

^{2}T^{−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 [–]

## INTRODUCTION

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.

## METHODOLOGY

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

*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),

*S*and

_{f}*S*

_{0}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:

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).

*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):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):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:

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 *L*_{2}.

### 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

*u*

^{2}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:

#### The continuity equation

*et al.*2005; López-Barrera

*et al.*2012). Replacing the velocity obtained from Equation (7) into the continuity (Equation (1)) yields:

*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:

*A*,

*B*and

*C*are defined in the following forms: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:

### 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 *S*_{0}/*a*.

*Upstream boundary condition:* the upstream boundary condition (UBC) based on the discharge *Q* (*x*_{in}, *t*) is automatically obtained considering an inflow hydrograph, and then using the upstream rating curve (i.e., *h*_{up} = *rQ ^{s}* where

*r*and

*s*are the rating curve constants), the UBC with respect to the hydraulic head is determined (i.e.,

*h*(

*x*

_{in},

*t*) =

*h*

_{up}).

*Downstream boundary condition:**h* (*x*_{exit}*, t*), defined here as *h*_{exit}, 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.

- 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:where
*q*is the discharge per unit width (L^{2}T^{−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: - 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: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 (L^{0.5}T^{−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:where*e*is the void ratio (dimensionless),*r*is a constant (dimensionless) representing the particle surface area efficiency that is around 1.3 for crushed stones._{e}

*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

*θ*in Figure 3(a), is approximately equal to

_{ff}*θ*

_{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: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:where

*h*

_{down}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:

#### 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.

### 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*L*_{1},*L*_{2}and*L*_{3}as illustrated in Figure 2 - 2.
At first attempt, the length of the solution domain is considered to be

*L*_{2}as an effective length computed: where is the horizontal distance from the point of seepage entrance to the dam heel and*L*_{3}is the length between the dam toe and the final cell face at*n +*1. At the initial time,*L*_{3}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.
*L*_{2}and*L*_{3}are updated: . - 6.
The steps 3–5 are repeated until gets less than a desirable value (equal to ).

*L*

_{2}. 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.

## RESULTS

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):where
*h*_{0}is the water depth at the start point (the exit height),*h*is the water depth at the distance_{x}*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*S*_{0}= 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.

*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

*h*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.

_{up}. | DBC computed from Equation (15) (Hansen et al. 2005). | DBC computed from Equation (21) (Neumann B. C.) . | Analytical (Hansen 1992) . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

h_{up} (m)
. | h_{exit} (m)
. | MAE h_{exit} (%)
. | u_{exit} (m/s)
. | MAE u_{exit} (%)
. | h_{exit} (m)
. | MAE h_{exit} (%)
. | u_{exit} (m/s)
. | MAE u_{exit} (%)
. | h_{exit} (m)
. | u_{exit} (m/s)
. |

1 | 0.0374 | 1.5 | 0.4491 | 1.6 | 0.0384 | 1.1 | 0.4375 | 1.0 | 0.0380 | 0.4421 |

2 | 0.0803 | 0.5 | 0.6588 | 0.5 | 0.0813 | 0.7 | 0.6506 | 0.9 | 0.0807 | 0.6555 |

3 | 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) . | |||||||
---|---|---|---|---|---|---|---|---|---|---|

h_{up} (m)
. | h_{exit} (m)
. | MAE h_{exit} (%)
. | u_{exit} (m/s)
. | MAE u_{exit} (%)
. | h_{exit} (m)
. | MAE h_{exit} (%)
. | u_{exit} (m/s)
. | MAE u_{exit} (%)
. | h_{exit} (m)
. | u_{exit} (m/s)
. |

1 | 0.0374 | 1.5 | 0.4491 | 1.6 | 0.0384 | 1.1 | 0.4375 | 1.0 | 0.0380 | 0.4421 |

2 | 0.0803 | 0.5 | 0.6588 | 0.5 | 0.0813 | 0.7 | 0.6506 | 0.9 | 0.0807 | 0.6555 |

3 | 0.1244 | 1.9 | 0.8207 | 1.9 | 0.1254 | 1.1 | 0.8142 | 1.1 | 0.1268 | 0.8052 |

Test . | Equation . | MAE of h_{exit} (%)
. | MAE of u_{exit} (%)
. | 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 |

Test . | Equation . | MAE of h_{exit} (%)
. | MAE of u_{exit} (%)
. | 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.

Test . | Model . | MAE of Q (%)
. | Phreatic surface RMSE (m) . | Velocity RMSE (m/s) . |
---|---|---|---|---|

Case I (θ_{down} = 18.44°) | 2D FEM (Larese et al. 2015) | 0 | 0.0090 | 0.0008 |

non-inertia wave model | 2.0 | 0.0095 | 0.0014 | |

Case II (θ_{down} = 24.44°) | 2D FEM (Larese et al. 2015) | 0 | 0.0104 | 0.0010 |

non-inertia wave model | 5.8 | 0.0158 | 0.0013 |

Test . | Model . | MAE of Q (%)
. | Phreatic surface RMSE (m) . | Velocity RMSE (m/s) . |
---|---|---|---|---|

Case I (θ_{down} = 18.44°) | 2D FEM (Larese et al. 2015) | 0 | 0.0090 | 0.0008 |

non-inertia wave model | 2.0 | 0.0095 | 0.0014 | |

Case II (θ_{down} = 24.44°) | 2D FEM (Larese et al. 2015) | 0 | 0.0104 | 0.0010 |

non-inertia wave model | 5.8 | 0.0158 | 0.0013 |

### Unsteady-state condition

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

*h*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 Δ_{up}*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.

*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.

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.

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 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.

## MODEL SENSITIVITY AND STABILITY ANALYSIS

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

*h*

_{N.B.C}and

*h*

_{observed}, respectively. Furthermore, the ratio of Δ

*x*to the effective length (

*L*

_{2}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

*S*

_{0}= 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 m

^{3}/s, 0.0925 m, 0.0183 m

^{3}/s, 0.1019 m, 0.0227 m

^{3}/s and 0.10495 m for the particle sizes described above, respectively. The result of the comparison is shown in Figure 19.

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 *h*_{up}), 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.5}/Δ*x*. 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.

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.

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.

## CONCLUSION

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.

## DATA AVAILABILITY STATEMENT

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

## REFERENCES

*Doctoral dissertation*

*Doctoral dissertation*