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−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 [–]
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
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).
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.
Rectangular rockfill dam: the head values and the properties are defined in the cells, whereas the velocity and harmonic head defined at the faces.
Rectangular rockfill dam: the head values and the properties are defined in the cells, whereas the velocity and harmonic head defined at the faces.
Trapezoidal rockfill dam: the solution domain (L2) varies with time.
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
The continuity equation


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.
- 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 (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:
- 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 (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: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.
(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)).
(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)).
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 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
).
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 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.
Comparison of the model with analytical solution (Hansen 1992) using the Neumann-type B.C.
Comparison of the model with analytical solution (Hansen 1992) using the Neumann-type B.C.
Comparison of the non-inertia wave model with experimental data GF29d (Hansen 1992) using the Neumann-type B.C.
Comparison of the non-inertia wave model with experimental data GF29d (Hansen 1992) using the Neumann-type B.C.
Comparison of the non-inertia wave model with experimental data GF31c (Hansen 1992) using the Neumann-type B.C.
Comparison of the non-inertia wave model with experimental data GF31c (Hansen 1992) using the Neumann-type B.C.
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) . |
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) . | |||||||
---|---|---|---|---|---|---|---|---|---|---|
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) . |
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 |
Evaluation of the various DBC equations
Test . | Equation . | MAE 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 |
Test . | Equation . | MAE 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.
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.
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.
(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).
(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).
Comparison of the non-inertia wave model and 2D FEM (Larese et al. 2015) with observed 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 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.
Comparison of the FVM non-inertia wave model using the Neumann-type B.C. with the SPH model (Peng et al. 2017).
Comparison of the FVM non-inertia wave model using the Neumann-type B.C. with the SPH model (Peng et al. 2017).
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).
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).
The UBC for gradually varied unsteady flow (GVUF) obtained from the experimental data (Hosseini & Joy 2007).
The UBC for gradually varied unsteady flow (GVUF) obtained from the experimental data (Hosseini & Joy 2007).
Comparison of the non-inertia with dynamic wave model and experimental data (Hosseini & Joy 2007) within GVUF.
Comparison of the non-inertia with dynamic wave model and experimental data (Hosseini & Joy 2007) within GVUF.
Error arising from the non-inertia and dynamic wave models compared to the experimental data (Hosseini & Joy 2007) within GVUF.
Error arising from the non-inertia and dynamic wave models compared to the experimental data (Hosseini & Joy 2007) within GVUF.
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.
The UBC and DBC for rapidly varied unsteady flow (RVUF) obtained from the experimental data (Liu et al. 1999).
The UBC and DBC for rapidly varied unsteady flow (RVUF) obtained from the experimental data (Liu et al. 1999).
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.
Comparison of the non-inertia with dynamic wave model and experimental (Liu et al. 1999) data within RVUF.
Comparison of the non-inertia with dynamic wave model and experimental (Liu et al. 1999) data within RVUF.
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 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.
(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).
(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).
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).
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).
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.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.
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.
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.
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.
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.
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.
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.
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).
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).
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.