Non-Equilibrium Bedload Transport Model Applied to Erosive Overtopping Dambreach

Abstract

Bedload sediment transport is an ubiquitous process in natural surface water flows (rivers, dams, coast, etc), but it also plays a key role in catastrophic events such as dyke erosion or dam breach collapse. The bedload transport mechanism can be under equilibrium state, where solid rate and flow carry capacity are balanced, or under non-equilibrium (non-capacity) conditions. Extremely transient surface flows, such as dam/dyke erosive collapses, are systems which always change in space and time, hence absolute equilibrium states in the coupled fluid/solid transport rarely exist. Intuitively, assuming non-equilibrium conditions in transient flows should allow to estimate correctly the bedload transport rates and the bed level evolution. To get insight into this topic, a 2D Finite Volume model for bedload transport based on the non-capacity approach is proposed in this work. This non-equilibrium model considers that the actual bedload sediment discharge can be delayed, spatial and temporally, from the instantaneous solid carry capacity of the flow. Furthermore, the actual solid rate and the adaptation length/time is governed by the temporal evolution of the bedload transport layer and the vertical exchange solid flux. The model is tested for the simulation of overtopping dyke erosion and dambreach opening cases. Numerical results seems to support that considering non-equilibrium conditions for the bedload transport improves the general agreement between the computed results and measured data in both benchmarking cases.

1. Introduction

In surface water flows, such as river or coasts, the bed material particles are transported by two main mechanisms: suspended load transport accounts for the sediment particles moving within the flow column, whereas bedload transport involves the particles moving in a relatively thin layer over the static bed stratum. In this bedload transport layer, the solid particles are in contact with the underlying static stratum and hence the particle velocity is much lower than that of the the surrounding fluid [1]. The bedload transport plays a key role in the riverbed morphology since it can lead to marked and rapid variations in the bed level. Bedload transport can occur under a steady equilibrium state or, contrarily, under transient non-equilibrium conditions. The classical equilibrium approach considers that the actual bedload rate is the capacity of the flow to carry solid weight and hence only depending on instantaneous local flow features. The capacity bedload rate can be estimated using different empirical closure relations found in literature [2] and is the most widespread approach in current bedload transport models [3,4,5,6,7,8,9,10].

Contrarily, in bedload models based on the non-equilibrium assumption, the actual transport rates must be computed through the advection of the solid particles in the bedload layer and the mass exchange with the underlying static riverbed [11,12,13,14,15], and usually it does not agree with the local load carry capacity of the flow. Intuitively, assuming non-equilibrium conditions in transient flows should allow to estimate correctly the solid transport rates and the bed evolution, but this topic remains uncertain nowadays. Van-Rijn [16] suggested that the bedload rate can be assumed equal to the solid carry capacity event in moderate unsteady flows, since the adjustment of the transported solid load to new local flow conditions proceeds enough rapid close to the riverbed surface. However, this assumption was only conceptual and not supported by theoretical or numerical results. Later, analysing numerically the time scale of the processes involved in bedload transport, Cao et al. [13,17] found that the bedload rate was able adapt rapidly to the local carry capacity in flood cases. These works seem to justify the widespread application of the equilibrium/capacity bedload models [17]. However, this fact is not supported by any results in case of highly transient erosive dambreach flows.

Assuming non-equilibrium conditions for the bedload transport requires to compute locally both the thickness of the bedload transport layer and the bulk velocity of the solid particles. Furthermore, the net exchange flux of solid particles between static riverbed and moving bedload layer also has to be estimated, since it works as a source/sink of solid load. Therefore, non-equilibrium bedload models for transient flows involve a high uncertainty in the estimation of the parameters associated to entrainment/deposition vertical flux and the transport layer thickness and velocity. In the specialized literature, different strategies have been reported to overcome this uncertainty. El Kadi Abderrezzak et al. [18] assumed non-storage conditions for the transport layer and approximated the net solid flux through the interface between the static bed and the moving transport layer by the divergence between the actual and capacity transport rates. This theoretical model also introduced the adaptation length concept, fixed in time and space, leading to the quasi-steady relaxation flux approach for bedload models. Based on the same adaptation length concept, refs. [11,19,20] proposed more complex approaches to determine the net exchange flux through the interface but accounting now for the mass storage in the bedload transport layer. These models assume the actual transport rate as one of the unknown conservative variables, but they estimate the bulk bedload layer velocity using the closure relations found in literature [21,22,23] developed in laboratory under solid carry capacity conditions.

Recently, novel complex developments applying the non-capacity approach in bedload transport models have been derived, most of them based on the pioneering work Charru [24]. This work studies, at a grain scale, the mechanical interaction between moving fluid (clear water) and sediment particles at the interface between the static riverbed and the moving bedload layer. Based on this particle-scale momentum balance, Zech et al. [25] modelled the net exchange flux through the interface as a function of the shear stress exerted by the flow in the bedload transport layer and the shear strength for the incipient motion at the upper boundary of the static underlying stratum. Based also on [24], Fernández-Nieto et al. [14,26] derived closure relationships for the erosion and deposition rates at the static-moving interface, leading to a novel expression for the actual solid transport rate under non-equilibrium conditions. Under equilibrium conditions, this new expression for the solid transport reduces to the classical formula for bedload rate. Some interesting works have been recently published modelling bedload transport rates in cohesive clay–silt–gravel mixtures [27,28]. Finally, Bohorquez and Ancey [15] performed a grain-scale analysis of the particle activity at the bedload transport layer and proposed expressions for the entrainment/deposition rates including an additional diffusion term in the non-capacity bedload transport equation.

In this work, the novel generalized non-equilibrium bedload transport formulation developed by Martínez-Aranda et al. [29,30] for 2D bedload transport models is used to simulate dyke overtopping erosion and dam breach collapse processes. The equilibrium and non-equilibrium approaches are compared for different synthetic tests and experimental benchmarking cases. The main goal of this work is to test the suitability of the non-equilibrium approach for these highly unsteady erosive flows. This manuscript organized as follows: in Section 2 the system of conservation equations for 2D bedload transport is introduced; the numerical method is detailed in Section 3; the model results for one synthetic test and two experimental benchmarking cases are analysed in Section 4; and finally the main conclusions are drawn in Section 5.

2. Mathematical Modelling of 2D Generalized Bedload Transport

The 2D governing equations for bedload transport assuming non-equilibrium approach can be written as a compact conservation laws system [30], which includes the conservation laws for the water mass and the depth-averaged momentum of the flow, together with the transport equation for the erodible bed material:

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}\left(\mathbf{U}\right)}{\partial x}+\frac{\partial \mathbf{G}\left(\mathbf{U}\right)}{\partial y}={\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)+{\mathbf{S}}_{\mathbf{\tau }}\left(\mathbf{U}\right)$$

(1)

where $\mathbf{U}={(h,hu,hv,z_b)}^T$ are the conservative variables, $\mathbf{F}\left(\mathbf{U}\right)$ and $\mathbf{G}\left(\mathbf{U}\right)$ are the conservative fluxes along the x- and y-coordinates respectively:

$$\begin{array}{c}\hfill \mathbf{U}=\left(\begin{array}{c}h\\ hu\\ hv\\ z_b\end{array}\right)\mathbf{F}\left(\mathbf{U}\right)=\left(\begin{array}{c}hu\\ h{u}^2+\frac{1}{2}g{h}^2\\ huv\\ q_{bx}\end{array}\right)\mathbf{G}\left(\mathbf{U}\right)=\left(\begin{array}{c}hv\\ huv\\ h{v}^2+\frac{1}{2}g{h}^2\\ q_{by}\end{array}\right)\end{array}$$

(2)

being h the flow depth, ($u,v)$ the components of the depth-averaged velocity $\mathbf{u}$ along the x- and y-coordinates respectively, $z_b$ is the erodible bed elevation, $(q_{bx},q_{by})$ are the components of the non-equilibrium bedload discharge ${\mathbf{q}}_{\mathbf{b}}$ along the x- and y-coordinates respectively and g is the gravity acceleration. On the right hand side, the terms ${\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)$ includes the bed-pressure momentum source and ${\mathbf{S}}_{\mathbf{\tau }}\left(\mathbf{U}\right)$ accounts for the frictional momentum dissipation along the x- and y-coordinates:

$${\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)=\left(\begin{array}{c}0\\ -gh\frac{\partial z_b}{\partial x}\\ -gh\frac{\partial z_b}{\partial y}\\ 0\end{array}\right){\mathbf{S}}_{\mathbf{\tau }}\left(\mathbf{U}\right)=\left(\begin{array}{c}0\\ -{\tau }_{bx}/{\rho }_w\\ -{\tau }_{by}/{\rho }_w\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ -gh{S}_{fx}\\ -gh{S}_{fy}\\ 0\end{array}\right)$$

(3)

being $({\tau }_{bx},{\tau }_{by})$ the components of the boundary shear stress vector ${\mathbf{\tau }}_{\mathbf{b}}$ between the flow and the movable bed along the x- and y-coordinates, respectively, and ($S_{fx},S_{fy}$) the corresponding frictional energy loss. For the sake of simplicity, in this work the turbulent Manning friction formulation is adopted:

$$S_{fx}=\frac{n^2\left|\mathbf{u}\right|u}{h^{4/3}}{S}_{fy}=\frac{n^2\left|\mathbf{u}\right|v}{h^{4/3}}$$

(4)

where $\left|\mathbf{u}\right|=\sqrt{u^2+v^2}$ is the depth-averaged flow velocity modulus and n is the bulk Manning roughness parameter depending on the bed material.

The bedload solid discharge in classical morpho-dynamical models is considered under equilibrium (capacity) conditions and usually modelled using the Grass law [31], which scales exponentially the solid transport rate with the flow velocity as ${\mathbf{q}}_{\mathbf{b}}=A_g{\left|\mathbf{u}\right|}^{m-1}\mathbf{u}$, being $A_g$ a flow-bed interaction constant factor and m an behaviour exponent (usually $m=3$). More complex equilibrium models applied to realistic bedload transport tests estimate the instantaneous local flow-bed interaction factor from the hydro-dynamical flow features, hence $A_g\equiv G(h,\theta )$ [2,10,32], being $\theta$ the dimensionless Shields stress at the bed surface. According to Martínez-Aranda et al. [29] and Martínez-Aranda et al. [30], the bedload transport rate assuming a non-equilibrium approach can be expressed following a modified Grass-type law as:

$$q_{bx}=\frac{1}{1-\xi }{G(h,\theta ,\eta )\left|\mathbf{u}\right|}^{2}u{q}_{by}=\frac{1}{1-\xi }G(h,\theta ,\eta ){\left|\mathbf{u}\right|}^{2}v$$

(5)

being $\xi$ the bed layer porosity and G a interaction factor between the flow and erodible riverbed. This interaction factor depends not only on the flow depth h and the dimensionless Shields stress $\theta$, but also on the bedload transport layer thickness $\eta$. Figure 1 shows a schematic representation of the bedload generalized bedload mechanism at grain-scale.

It is useful to split G in a product of three factors ${\Gamma }_1\left(h\right)$, ${\Gamma }_2\left(\theta \right)$ and ${\Gamma }_3\left(\eta \right)$ as:

$$G={\Gamma }_1\left(h\right){\Gamma }_2\left(\theta \right){\Gamma }_3\left(\eta \right)$$

(6)

Table 1. Grass-type interaction factor G for transport rate formulations.

Formulation	${\mathbf{\Gamma }}_1\left(\mathit{h}\right)$	${\mathbf{\Gamma }}_2\left(\mathit{\theta }\right)$	${\mathbf{\Gamma }}_3\left(\mathit{\eta }\right)$	${\mathit{\theta }}_{\mathit{c}}$
MPM	$\frac{n^3\sqrt{g}}{(r_s-1)\sqrt{h}}$	$8\sqrt{\Delta \theta }{\theta }^{-3/2}$	$\frac{r_s{k}_d}{k_e}\frac{\eta }{d_s}$	0.047
Nielsen	$\frac{n^3\sqrt{g}}{(r_s-1)\sqrt{h}}$	$12{\theta }^{-1}$	$\frac{r_s{k}_d}{k_e}\frac{\eta }{d_s}$	0.047
Fernandez-Luque	$\frac{n^3\sqrt{g}}{(r_s-1)\sqrt{h}}$	$5.7\sqrt{\Delta \theta }{\theta }^{-3/2}$	$\frac{r_s{k}_d}{k_e}\frac{\eta }{d_s}$	0.037
Wong	$\frac{n^3\sqrt{g}}{(r_s-1)\sqrt{h}}$	$3.97\sqrt{\Delta \theta }{\theta }^{-3/2}$	$\frac{{\rho }_s/{\rho }_f{k}_d}{k_e}\frac{\eta }{d_s}$	0.0495
Smart	$\frac{n^2}{(r_s-1)h^{1/3}}$	$4.2S_b^{0.6}{\theta }^{-1}$	$\frac{r_s{k}_d}{k_e}\frac{\eta }{d_s}$	0.047
Wu	$\frac{n^3\sqrt{g}}{(r_s-1)\sqrt{h}}$	$0.0053{\theta }_c^{-2.2}\sqrt{\Delta \theta }{\theta }^{-3/2}$	$\frac{r_s{k}_d}{k_e}\frac{\eta }{d_s}$	0.030

In the Smart formula, the parameter $S_b$ accounts for the bed slope along the bedload transport direction.

provides expressions for functions ${\Gamma }_1$, ${\Gamma }_2$ and ${\Gamma }_3$ using empirical relations found in literature, where $r_s$ is the relative density of the sediment particles $r_s={\rho }_s/{\rho }_w$, $d_s$ is the characteristic diameter of the bed sediment, $k_d$ and $k_e$ are the two constants associated to the deposition and erosion exchange rates, respectively, between the underlying static bed stratum and the bedload transport layer (moving layer). The term $\Delta \theta$ accounts for the excess of Shields stress over the critical threshold ${\theta }_c$ for the incipient motion threshold of the solid particles, calculated as:

$$\begin{array}{cc}\hfill & \Delta \theta =\left\{\begin{array}{c}\hfill \theta -{\theta }_c\mathrm{if}\theta >{\theta }_c\hfill \\ \hfill 0\mathrm{Otherwise}\hfill \end{array}\right.\hfill \\ \hfill & \mathrm{with}:\theta =\frac{|{\mathbf{\tau }}_{\mathbf{b}}|}{({\rho }_s-{\rho }_w)g{d}_s}=\frac{n^2{\left|\mathbf{u}\right|}^2}{(r_s-1)d_s{h}^{1/3}}\hfill \end{array}$$

(7)

This generalized Grass-type model (6) allows to consider capacity (equilibrium) or non-equilibrium approaches for the bedload transport only by selecting the proper estimation procedure for the bedload layer thickness $\eta$. On one hand, assuming the equilibrium approach leads to the instantaneous adaptation of the bedload transport rate ${\mathbf{q}}_{\mathbf{b}}$ to its capacity value ${\mathbf{q}}_{\mathbf{b}}^{*}$ and the equilibrium transport layer thickness ${\eta }^{*}$ can be estimated as:

$$\eta \equiv {\eta }^{*}=\frac{k_e}{r_s{k}_d}\Delta \theta d_s$$

(8)

On the other hand, the non-equilibrium assumption forces the necessity of estimating the temporal evolution of the bedload transport layer thickness $\eta$ using the conservation equation:

$$\frac{\partial \eta }{\partial t}+\frac{\partial q_{bx}}{\partial x}+\frac{\partial q_{by}}{\partial y}=({\dot{\eta }}_e-{\dot{\eta }}_d)$$

(9)

where ${\dot{\eta }}_e$ and ${\dot{\eta }}_d$ are the entrainment and deposition vertical rates, respectively, between the the underlying static stratum and the bedload transport layer. Following [24,29], the volumetric entrainment and deposition rates can be expressed as:

$${\dot{\eta }}_e=k_e\frac{\Delta \theta }{r_s}\sqrt{(r_s-1)g{d}_s}{\dot{\eta }}_d=k_d\frac{\eta }{d_s}\sqrt{(r_s-1)g{d}_s}$$

(10)

3. Numerical Scheme for the 2D Bedload Transport Model

In order to obtain the numerical solution applying Finite Volume (FV) methods, the spatial domain is discretized into a fixed-in-time mesh of triangular cells and system (1) is integrated in each computational cell using the Gauss theorem:

$$\frac{d}{dt}{\int }_{{\Omega }_i}\mathbf{U}d\Omega +\sum _{k=1}^{NE}{\left[\mathbf{F}\left(\mathbf{U}\right)n_x+\mathbf{G}\left(\mathbf{U}\right)n_y\right]}_k{l}_k={\int }_{{\Omega }_i}{\mathbf{S}}_{\mathbf{b}}\left(\mathbf{U}\right)d\Omega +{\int }_{{\Omega }_i}{\mathbf{S}}_{\mathbf{\tau }}\left(\mathbf{U}\right)d\Omega$$

(11)

being ${\Omega }_i$ the inner domain and $NE$ the number of edges for the i cell, $\mathbf{n}=(n_x,n_y)$ the outward unit normal vector and $l_k$ the length of the kth edge. Assuming a piecewise representation of the conservative variables $\mathbf{U}$ at each i cell, ${\mathbf{U}}_i$, the updating formulation between from time $t=t^n$ to the next step $t=t^{n+1}$ is expressed as:

$$\begin{array}{c}\hfill {\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n-\frac{\Delta t}{A_i}\sum _{k=1}^{\mathrm{NE}}{\mathbf{R}}_k^{-\mathbf{1}}{\mathcal{F}}_k^{\downarrow -}{l}_k\end{array}$$

(12)

where $A_i$ is the i cell area, $\Delta t=t^{n+1}-t^n$ is the time step, ${\mathbf{R}}_k^{-\mathbf{1}}$ is the inverse of the rotation matrix ${\mathbf{R}}_k$ at kth edge [33] and ${\mathcal{F}}_k^{\downarrow -}$ is the numerical normal flux throughout the kth edge. This numerical normal flux includes the convective flux augmented with the non-conservative bed-pressure contribution and the frictional momentum dissipation integrated thought the edge.

In this work, the upwind computation of ${\mathcal{F}}_k^{\downarrow -}$ is achieved using a fully-coupled Roe-type Riemann solver for bedload transport (FCM solver) [34]. This augmented 5-waves solver is based on computing the approximated solution of the linear Riemann problem:

$$\begin{array}{c}{\displaystyle \frac{\partial {\mathbf{R}}_k\mathbf{U}}{\partial t}}+{\tilde{\mathbf{M}}}_k{\displaystyle \frac{\partial {\mathbf{R}}_k\mathbf{U}}{\partial \widehat{n}}}={\mathbf{R}}_k{\tilde{\mathbf{T}}}_k\end{array}$$

(13)

where $\widehat{n}$ denotes the normal axis to the cell edge and superscripts ^˜ indicates edge-averaged quantities. The pseudo-Jacobian matrix ${\tilde{\mathbf{M}}}_k={(\tilde{\mathbf{J}}-\tilde{\mathbf{H}})}_k$ of the coupled system includes the approximated Jacobian matrix of the conservative convective flux ${\tilde{\mathbf{J}}}_k$ and the non-conservative matrix ${\tilde{\mathbf{H}}}_k$ of the bed-pressure momentum contribution at the edge. On the right hand side, the term ${\tilde{\mathbf{T}}}_k$ denotes the momentum dissipation due to flow-bed frictional stress, integrated through the kth edge. Details on the flux computation have been extensively reported in Martínez-Aranda et al. [34] but, for the sake of brevity, ${\mathcal{F}}_k^{\downarrow -}$ at the kth edge is computed as:

$${\mathcal{F}}_k^{\downarrow -}=\mathbf{F}\left({\mathbf{R}}_k{\mathbf{U}}_i^n\right)+\sum _{m-}{\left[({\tilde{\lambda }}_m{\tilde{\alpha }}_m-{\tilde{\sigma }}_m){\tilde{\mathbf{e}}}_m\right]}_k^n$$

(14)

where ${\tilde{\lambda }}_{m,k}$ are the wave celerities at the edge, i.e., the eigenvalues of the pseudo-Jacobian matrix ${\tilde{\mathbf{M}}}_k$ whereas ${\left({\tilde{\mathbf{e}}}_m\right)}_k$ are the corresponding eigenvectors. The coefficients ${\tilde{\alpha }}_{m,k}$ denotes the wave strengths which accounts for the discontinuity on the conservative variables, whereas the coefficients ${\tilde{\sigma }}_{m,k}$ are the source strengths involving the frictional momentum contribution at the edge. The subscript $m-$ under the sums indicate that only the waves travelling inward the i cell are considered, leading to the upwind computation of the flux at the edge.

Furthermore, the stability of the explicitly computed numerical solution is addressed by the dynamical limitation of the global time step using a CFL condition. The local time step allowed for the kth edge is estimated assuming that the absolute maximum of the eigenvalues corresponds to the fastest wave celerity. Hence the global time step $\Delta t$ is limited by minimum of the local time steps, i.e., the more restrictive instantaneous local flow features (see Martínez-Aranda et al. [34]).

3.1. Transport Layer Updating with Capacity and Non-Capacity Approaches

Once the conservative variables have been updated using (12), the bedload transport rate ${\mathbf{q}}_{\mathbf{b}}$ is computed at the next time level $t^{n+1}$ as:

$$\begin{array}{cc}\hfill & {\left({\mathbf{q}}_{\mathbf{b}}\right)}_i^{n+1}=\frac{1}{1-\xi }{G}_i^{n+1}\left[{\left(u_i^{n+1}\right)}^2+{\left(v_i^{n+1}\right)}^2\right]{\mathbf{u}}_i^{n+1}\hfill \\ \hfill & \mathrm{with}:G_i^{n+1}={\Gamma }_1\left(h_i^{n+1}\right){\Gamma }_2\left({\theta }_i^{n+1}\right){\Gamma }_3\left({\eta }_i^{n+1}\right)\hfill \end{array}$$

(15)

where the hydrodynamic quantities $h_i^{n+1}$, $u_i^{n+1}$, $v_i^{n+1}$ and ${\theta }_i^{n+1}$ are directly extracted from ${\mathbf{U}}_i^{n+1}$. Nevertheless, Equation (15) also requires to compute the updated bedload transport layer thickness ${\eta }_i^{n+1}$ at the next time level. The updating of the transport layer thickness depends on the assumption made for the bedload transport:

• Equilibrium hypothesis: The new transport layer thickness is directly computed as:

  $${\eta }_i^{n+1}=\frac{k_e}{r_s{k}_d}{(\Delta \theta )}_i^{n+1}{d}_s$$

  (16)

  where ${(\Delta \theta )}_i^{n+1}$ is the Shields excess (7) at the i cell computed from the updated conservative variables ${\mathbf{U}}_i^{n+1}$.
• Non-capacity approach: This leads to the necessity of solving Equation (9) each time step. The updating formula for the transport layer thickness $\eta$ is expressed as:

  $${\eta }_i^{n+1}={\eta }_i^n-\frac{\Delta t}{A_i}\sum _{k=1}^{NE}{F}_k^{\eta \downarrow }{l}_k+\Delta t{({\dot{\eta }}_e-{\dot{\eta }}_d)}_i^n$$

  (17)

  being $F_k^{\eta \downarrow }$ the numerical normal flux at the kth intercell edge for the transport layer equation and ${\dot{\eta }}_{e,i}^n$ and ${\dot{\eta }}_{d,i}^n$ the values of the entrainment and deposition rates for the i cell at time n, respectively. To compute the numerical flux $F_k^{\eta \downarrow }$, the bedload transport layer Equation (9) is projected along the normal direction to the edge $\widehat{n}$ and approximated by the local scalar Riemann problem:

  $$\begin{array}{cc}\hfill & \frac{\partial \eta }{\partial t}+{\tilde{\lambda }}_{\eta ,k}\frac{\partial \eta }{\partial \widehat{n}}=0\hfill \end{array}$$

  (18)

  $$\begin{array}{cc}\hfill & \mathrm{with}:\eta (\widehat{n},0)=\left\{\begin{array}{ccc}{\eta }_i^n& \mathrm{if}& \widehat{n}<0\\ {\eta }_j^n& \mathrm{if}& \widehat{n}>0\end{array}\right.\hfill \end{array}$$

  (19)

  where the virtual bedload wave celerity ${\tilde{\lambda }}_{\eta ,k}$ is defined as:

  $${\tilde{\lambda }}_{\eta ,k}={\left(\frac{\delta q_{bn}}{\delta \eta }\right)}_k$$

  (20)

  being $\delta q_{bn}={(q_{bx}{n}_x+q_{by}{n}_y)}_j^n-{(q_{bx}{n}_x+q_{by}{n}_y)}_i^n$ and $\delta \eta ={\eta }_j^n-{\eta }_i^n$, and the bedload numerical flux $F_k^{\eta \downarrow }$ is computed as:

  $$\begin{array}{cc}\hfill & F_k^{\eta \downarrow }=\left\{\begin{array}{ccc}{(q_{bx}{n}_x+q_{by}{n}_y)}_i^n& \mathrm{if}& {\tilde{\lambda }}_{\eta ,k}>0\\ {(q_{bx}{n}_x+q_{by}{n}_y)}_j^n& \mathrm{if}& {\tilde{\lambda }}_{\eta ,k}<0\end{array}\right.\hfill \end{array}$$

  (21)
  The cell-centered exchange rates ${\dot{\eta }}_{e,i}^n$ and ${\dot{\eta }}_{d,i}^n$ between the underlying static stratum and the transport layer are computed as:

  $${\left({\dot{\eta }}_e\right)}_i^n=k_e\frac{{(\Delta \theta )}_i^n}{r_s}\sqrt{(r_s-1)g{d}_s}{\left({\dot{\eta }}_d\right)}_i^n=k_d\frac{{\eta }_i^n}{d_s}\sqrt{(r_s-1)g{d}_s}$$

  (22)

3.2. Morphological Collapse Mechanism

In realistic applications, the maximum bed level slope is limited due to morphological collapse processes. In order to implement an efficient bed slope control mechanism, we propose the non-iterative limitation of the numerical beload flux provided for the augmented Riemann solver at each cell edge imposing geometrical restrictions. The bed level updating expressions for the left L and right R cells sharing the edge $k_0$ read:

$$\begin{array}{cc}\hfill & z_{bL}^{n+1}=z_{bL}^n-\frac{\Delta t}{A_L}\left[{\left(q_b^{\downarrow }l\right)}_{k_0}+\cdots +{\left(q_b^{\downarrow }l\right)}_{k_{NL}}\right]\hfill \\ \hfill & z_{bR}^{n+1}=z_{bR}^n-\frac{\Delta t}{A_R}\left[{\left(q_b^{\downarrow }l\right)}_{k_0}+\cdots +{\left(q_b^{\downarrow }l\right)}_{k_{NR}}\right]\hfill \end{array}$$

(23)

where $k_{NL}$ and $k_{NR}$ refer to the non-shared edges of the left and right cells respectively and ${\left(q_b^{\downarrow }\right)}_{k_0}$ denotes the forth component of the numerical flux vector ${\mathcal{F}}^{\downarrow -}$ at the common edge $k_0$ between cells L and R. This value is provided by the Riemann solver, i.e., the numerical upwind bedload normal flux at the edge ($q_b^{\downarrow }\equiv {\mathcal{F}}_4^{\downarrow -}$). Therefore, the evolution of the local bed slope at edge $k_0$ can be expressed as:

$${(z_{bR}-z_{bL})}^{n+1}={(z_{bR}-z_{bL})}^n+\frac{A_L+A_R}{A_L{A}_R}{\left(q_b^{\downarrow }l\right)}_{k_0}\Delta t$$

(24)

Assuming a collapse bed slope angle ${\delta }_b$ along the edge-normal distance between cell centres $d_n$, two conditions arise for the morphological limiting of the numerical bedload flux depending on the sign of the local bed slope:

• Positive bed slope if $0<{(z_{bR}-z_{bL})}^{n+1}<tan{\delta }_b{d}_n$

  $$q_b^{\downarrow }\le \frac{tan{\delta }_b{d}_n-{(z_{bR}-z_{bL})}^n}{l\Delta t}\frac{A_L{A}_R}{A_L+A_R}$$

  (25)
• Negative bed slope if $-(tan{\delta }_b{d}_n)<{(z_{bR}-z_{bL})}^{n+1}<0$

  $$q_b^{\downarrow }\ge \frac{-\left[tan{\delta }_b{d}_n+{(z_{bR}-z_{bL})}^n\right]}{l\Delta t}\frac{A_L{A}_R}{A_L+A_R}$$

  (26)

Actually, this limitation on the numerical bedload flux $q_b^{\downarrow }$ acts as a bed form control mechanism, including artificial diffusion in the numerical bed level solution and avoiding the appearance of shock waves in the bed level with local slopes larger than the morphological collapse limit.

4. Numerical Results and Discussion

4.1. Adaptation of the Non-Equilibrium Bedload Rate to Equilibrium States

This new generalised model for the non-equilibrium bedload transport (5)–(6) requires to compute the thickness $\eta$ of the transport layer using (9), but it is able to account for the spatial and temporal delay of the actual bedload rate with respect to the solid transport capacity in highly unsteady flows. Moreover, if the local flow features allow to develop steady states in the bedload transport process, the entrainment and deposition rates tend to balance (${\dot{\eta }}_e={\dot{\eta }}_d$) and the solid rate recovers the equilibrium condition, characterized by the transport layer thickness ${\eta }^{*}$ (8). Therefore, the equilibrium state is a particular case of the generalised non-equilibrium model where $\eta ={\eta }^{*}$.

In order to illustrate the temporal-spatial behaviour of the generalised non-equilibrium bedload model, a simple idealised test is included here forcing the following flow conditions

$$\begin{array}{cc}\hfill & h\left|\mathbf{u}\right|(t,x,y)=\begin{array}{cccc}2{\mathrm{m}}^2{\mathrm{s}}^{-1}\hfill & & \forall t,& \forall x\end{array}\hfill \\ \hfill & \left|\mathbf{u}\right|(t,x,y)=\left\{\begin{array}{cccc}2\mathrm{m}/\mathrm{s}\hfill & \mathrm{if}& t=0\mathrm{s},& \forall x\\ 4\mathrm{m}/\mathrm{s}\hfill & \mathrm{if}& t>0\mathrm{s},& 20\mathrm{m}\le x\le 70\mathrm{m}\\ 2\mathrm{m}/\mathrm{s}\hfill & \mathrm{if}& t>0\mathrm{s},& \mathrm{otherwise}\end{array}\right.\hfill \end{array}$$

with ${\rho }_s=2650$ kg/m${}^3$, $\xi =0.35$, $d_s=1$ mm, $n=0.035$, $k_e/k_d=20$ and the MPM relation (see Table 1). Using (15)–(16), these local flow features lead to a uniform bedload discharge at the initial time $t=0$ s and a stepped capacity transport rate $|{\mathbf{q}}_{\mathbf{b}}^{*}|$ for any $t>0$ s, as depicted in Figure 2.

The generalised non-equilibrium model introduces a temporal and spatial delay of the actual transport rate with respect to the capacity value. Figure 3 shows the temporal evolution of the bedload rate $|{\mathbf{q}}_{\mathbf{b}}|$ in the spatial domain computed using the non-equilibrium model (15)–(17), with (top) $k_e=0.06$ and (bottom) $k_e=0.015$. As time increases, the actual solid rate adapts progressively to the capacity rate. Nevertheless, sudden spatial changes of the local flow features, as occurs at $x=20$ m and $x=70$  m, need a length for the actual transport rate to adapt to the transport capacity value even with $t\to \infty$. Comparing the adaptation of the beload rate to the flow carry capacity with (top) $k_e=0.06$ and (bottom) $k_e=0.015$, it is worth noting that the temporal and spatial delay of the non-equilibrium solid discharge respect to the corresponding equilibrium state increases as the entrainment constant $k_e$ decreases.

Furthermore, the net exchange flux through the static-moving bed layers interface can be expressed as $N_{exch}={\dot{\eta }}_d-{\dot{\eta }}_e$ (22). Previous non-equilibrium (non-capacity) bedload models [12,20,35,36] assumed a spatial length $L_b$ for the adaptation of the actual bedload solid rate to its equilibrium state and estimate the net exchange vertical flux as

$$N_{exch}=\frac{1}{1-\xi }\frac{|{\mathbf{q}}_{\mathbf{b}}|-|{\mathbf{q}}_{\mathbf{b}}^{*}|}{L_b}$$

(27)

where $L_b$ is a constant parameter which needs to be calibrated for each case. Comparison of method (27) to determine the entrainment/deposition flux $N_{exch}$ with the method proposed in this work ($N_{exch}={\dot{\eta }}_d-{\dot{\eta }}_e$) demonstrates that the current model leads to a dynamic adaptation length $L_b$ which depends on the local flow features and scales following

$$L_b\propto \frac{\sqrt{\theta }}{k_e}{d}_s$$

(28)

According to (28), higher boundary Shields stresses $\theta$ are associated longer distances for the adaptation of the actual bedload discharge to its equilibrium state (see Figure 3). Contrarily, the lower the entrainment parameter $k_e$, the larger the dynamic adaptation length. Therefore, high Shields stresses $\theta$ and low entrainment factors $k_e$ enhance non-equilibrium bedload transport states, i.e., highly erosive flow features lead to larger adaptation lengths. Instead, most of the previous models assume a global constant value based on the dominant bed form [36].

To analyse the influence of the Shields stress on the dynamic adaptation of actual non-equilibrium bedload rate to the equilibrium state, we use the above case but setting a entrainment constant $k_e=0.025$ and increasing flow discharges $h\left|\mathbf{u}\right|(t,x,y)=[1,2,3,4]$ m${}^2$s${}^{-1}$. For $t=0$ s the flow is uniform, with $h=1$ m along the whole spatial domain regardless of the discharge, whereas for $t>0$ s a steady step in the flow features is set between 20 m $\le x\le$ 70 m so that $|{\mathbf{q}}_{\mathbf{b}}^{*}|=2|{\mathbf{q}}_{\mathbf{b}}^{*}{|}_{t=0}$ regardless of the flow discharge. The shear stress at the bed surface, hence the Shields stress $\theta$, increases progressively with the flow discharge as $\theta =[1.18,4.71,10.6,18.9]$ respectively. Figure 4 shows the non-equilibrium bedload rate $|{\mathbf{q}}_{\mathbf{b}}|$ for $t=100$ s, normalised by the initial uniform equilibrium value $|{\mathbf{q}}_{\mathbf{b}}^{*}{|}_{t=0}$, along the whole spatial domain. As the Shields stress at the bed surface increases, the adaptation length increases and enhances the non-capacity state in the bedload transport.

4.2. Dike Breaking by Overtopping Flow Erosion

This experimental test case was carried out by Tingsanchali and Chinnarasri [3] in a 35 m long 1 m wide straight rectangular cross-section flume. In the middle of the flume, a trapezoidal dyke was constructed with a non-cohesive sand of medium diameter $d_s=0.86$ mm, solid particle density ${\rho }_s=2650$ kg/m${}^3$ and associated bulk porosity $\xi =0.34$. The dyke height over the flume floor was $0.80$ m and the crest width was $0.30$ m, whereas the upward slope $S_u$ and downward slope $S_d$ depends on the laboratory case. The initial water surface elevation was $0.83$ m over the flume floor upstream the dyke. During each experimental case, a constant inlet discharge $q_{inlet}$ was set at the flume upstream section. In the original work [3], the temporal evolution of the flow discharge at the dyke crest section and the upstream reservoir level were reported during the whole experiment. Temporal bed elevation data were also provided at probes $P1=15$ cm, $P2=65$ cm and $P3=115$ cm downstream the upward vertex of the dyke crest.

For this case, simulations have been performed using a single-row mesh of 3500 rectangular cells with $\Delta x=0.01$ m and $\Delta y=1$ m. The Manning’s roughness coefficient used has been $n=0.010$ for the flume floor and $n=0.015$ for the sediment material. The CFL = $0.50$ was set for a full simulated period of 140 s. The Smart formula (Table 1) has been used to compute the bedload transport rate whereas the collapse mechanism has been neglected in this case. Two different benchmarking experiments, C1 and C2, have been used to study the capabilities of both the equilibrium and the non-equilibrium approaches (see

Table 2. Setup for benchmarking cases C1 and C2.

Case	${\mathit{S}}_{\mathit{u}}$	${\mathit{S}}_{\mathit{d}}$	${\mathit{q}}_{\mathbf{inlet}}$ (L/s)	${\mathit{k}}_{\mathit{e}}$	${\mathit{k}}_{\mathit{d}}$
C1	1V:3H	1V:5H	1.05	0.24	0.012
C2	1V:3H	1V:3H	1.23	0.34	0.017

).

The left column in Figure 5 shows the computed dam profile evolution for both benchmarking tests C1 and C2, considering the equilibrium and non-equilibrium bedload transport assumption, whereas the right column depicts the temporal evolution of the dam level at probes P1, P2 and P3. The black lines show the available experimental data for each test. Analysing the results obtained in both cases C1 and C2, the assumption of non-equilibrium bedload transport leads to a delay on the erosive process in the dam region and hence the change on the bed level is reduced for early stages after the flow starts. However, as the dam erosion progresses, the non-equilibrium bedload transport tends to approach to the equilibrium state and the dam profile level predicted by the non-equilibrium model converges to that obtained by the equilibrium model. It is worth noting that, in case C2, a marked scouring is observed in the experiments just downstream of the dam crest at time $t=60$ s. This scouring was caused likely by 3D behaviour of the flow at the upstream region of the eroded dam. The proposed depth-averaged bedload model is not able to capture the appearance and development of this feature in the bed profile, regardless of the equilibrium or the non-equilibrium approach is selected, due to its 3D nature.

Table 3. RMSE for dam level evolution data in benchmarking cases C1 and C2.

Data Series	Equil.	No-Equil.
C1—Bed level in P1	0.057 m	0.023 m
C1—Bed level in P2	0.053 m	0.022 m
C1—Bed level in P3	0.033 m	0.020 m
C2—Bed level profile at $t=30$ s	0.112 m	0.061 m
C2—Bed level profile at $t=60$ s	0.029 m	0.039 m

shows the Root Mean Square Error (RMSE) of the computed results with both the equilibrium and non-equilibrium models respect to the measured bed level evolution at P1, P2 and P3 in case C1 (available experimental data). Also, for case C2, the RMSE of the computed bed profiles at times $t=30$ s and $t=60$ s respect to the dam profiles measured in the laboratory for the same times are shown in Table 3. Generally, the computed dam level changes with the non-equilibrium approach show a much better agreement with the measured data than those obtained assuming the equilibrium restriction.

In order to get insight into the non-equilibrium bedload behaviour, the temporal evolution of the thickness of the transport layer $\eta$ at probes P1, P2 and P3, normalised by the corresponding equilibrium value ${\eta }^{*}$ computed from the instantaneous local hydrodynamic features using (8), is depicted in Figure 6 for cases C1 and C2. During the early stages after the overtopping flow starts, non-equilibrium states in the bedload transport layer occur all along the downward dam slope. The closer the dyke crest, the more enhanced is the non-equilibrium state of bedload transport and hence the normalised value $\eta /{\eta }^{*}$ generally shows lower values at P1 location than that at probes P2 and P3 during the first flow stages. Then, as the dam erosion progresses, the thickness of the transport layer tends progressively to the equilibrium value ($\eta /{\eta }^{*}\approx 1$) all along the dyke downward side.

The temporal evolution of the computed and measured discharge at the dyke crest and the water surface level at the the upstream reservoir have been plotted in Figure 7 for cases C1 and case C2. Also the Root Mean Square Error (RMSE) of the computed data with respect to the available measured data have been reported in

Table 4. RMSE for the crest discharge and reservoir level data in benchmarking cases C1 and C2.

Data Series	Equil.	No-Equil.
C1—Dam discharge	10.94 L/s	3.10 L/s
C1—Reservoir level	0.055 m	0.019 m
C2—Dam discharge	45.02 L/s	11.76 L/s
C2—Reservoir level	0.116 m	0.034 m

considering both equilibrium and non-equilibrium approaches. Considering non-equilibrium conditions for the bedload transport improves clearly the general agreement between the computed and measured evolution of the crest discharge and reservoir level in both benchmarking cases.

4.3. Breach Opening in Homogeneous Dam

This benchmark case aims to explore qualitatively the differences on the numerical modelling of the dambreach opening process using both bedload equilibrium and non-equilibrium approaches. Here we consider one of the experimental dambreach tests carried out at the Laboratório Nacional de Engenharia Civil (LNEC), Lisbon, Portugal, as part of the 6th Workshop on River and Sedimentation Hydrodynamics and Morphodynamics (Online, 24–26 November 2021), supported by Experimental Methods and Instrumentation Committee of the IAHR. The main results obtained during this laboratory campaign are unpublished yet. The experimental data available for this work were provided to the numerical modellers during the first blind stage and hence they are limited. Figure 8 shows a 3D sketch of the initial experimental setup. A homogeneous trapezoidal dam was constructed at the end of a 6 m length $1.2$ m width rectangular cross-section channel. The dam was $0.45$ m height, 1V:2H upward slope, 1V:2.5H downward slope, $2.20$ m base length and $0.17$ m crest width. Downstream the dam, there exists a $3.6 \mathrm{m}\times 2.0 \mathrm{m}\times 0.55 \mathrm{m}$ settling basin, where the elevation of the free surface is zero relatively to the bottom of the upstream reservoir (dam base).

The dam was composed of a non-cohesive silty sand with 25% of fines, solid density 2650 kg/m${}^3$, median diameter $d_s=0.320$ mm, relative compaction of 90% of the Standard Proctor reference value, 10% optimal water content for compaction. At the centerline of the dam crest, a triangular cross-section notch with 1V:1H side slopes and $1.5$ cm deep (see Figure 8) was carved to guide the dambreach process during the early flow stages. Initially, the water level at the upstream reservoir was equal to the height of the model dam ($0.45$ m) and the free surface at the downstream basin is at the dam toe elevation. Note that the water at the downstream basin is spilt if the elevation is higher than the dam toe. Once the flow started at the notch region and the breach opening process began, the water level at the upstream reservoir was maintained approximately constant by changing the upstream inlet discharge. The maximum inlet discharge allowed in the upstream reservoir is $0.030$ m${}^3$/s. The test was conducted until the maximum discharge at the upstream reservoir was attained. During the test, the temporal evolution of the inlet discharge at the upstream reservoir was recorded using an electromagnetic flowmeter whereas the water level was obtained from an acoustic probe following the free-surface.

This test has been simulated using an unstructured mesh of almost 16,500 triangular cells, with 50 cm${}^2$ of base cell area and a progressively refined region at the dam and notch zones of $0.5$ cm${}^2$ minimum cell area. At the inlet section of the upstream reservoir, the time series of inlet discharge recorded in the experiment has been imposed during the simulation, whereas free outflow has been allowed at the outlet section at the downstream basin. The Manning’s roughness coefficient $n=0.0125$ has been used for the sediment material in the Manning-Strickler relation $n=1/21d_s^{1/6}$. The porosity of the dam layer has been estimated as $\xi =0.335$ from the bulk weight of the compacted dam material (90% of the Standard Proctor) and the entrainment and deposition constants have been set to $k_e=0.06$ and $k_d=0.003$ respectively. The Smart formula (Table 1) has been used to compute the bedload transport rate with the critical Shields stress for the incipient motion ${\theta }_c=0.047$. The morphological mechanism plays a key role in this case to model the breach opening due to the collapse of the lateral sides. Based on experimental qualitative observations, a collapse bed slope angle ${\delta }_b=75{}^{\circ }$ has been set (see Section 3.2). CFL = $0.95$ for a full simulated period of 4500 s.

Figure 9, Figure 10, Figure 11 and Figure 12 show a 3D view of the dam topography and free water surface computed with the non-equilibrium model at $t=1500$ s, $t=2650$ s, $t=3566$ s and $t=4317$ s respectively. Also the thickness of the bedload transport layer and the flow velocity field at the breach region are shown. The breach opens progressively, increasing its cross-wise width due to the side wall collapse thanks to the collapse mechanism. A couple of coherent vortices appear downstream the dam, as the flow through the breach goes into the downward basin. The thickness of the transport layer increases progressively along the breach, allowing that the mobilized material incorporated to the flow from the breach sides can be transported downstream and reaches the settling basin.

Differences appears in the breach opening process depending on the equilibrium or non-equilibrium approaches selected. Figure 13 depicts (left) the total breach discharge and (right) the upstream reservoir level during the test computed with both equilibrium and non-equilibrium models. The total water discharge at the breach crest cross-section is almost similar using both approaches. Nevertheless, as the breach opens, the reservoir level decreases more markedly with the equilibrium assumption than that using the non-equilibrium model. Note that the experimental test was carried out maintaining an almost-constant water level at the upstream reservoir (solid black line in Figure 13 right). Therefore, the non-equilibrium model seems to approximate better the breach opening process than the equilibrium approach.

Finally, the left column of Figure 14 shows cross-section wetted-area at the dam crest section whereas the right column depicts the averaged flow velocity at the same section. Assuming non-equilibrium bedload transport leads to a slower flow velocity during the early stages of the breach opening and a higher wetted area along the breach region. This hydrodynamic combination reduces the erosive capability of the flow at this region and delays the breach opening process.

5. Conclusions

In this work, a recently developed 2D generalized bedload transport model [30], able to account for equilibrium and non-equilibrium transport conditions, is applied to study the dam breach and dyke erosion processes. This generalized model is based on the grain-scale analysis of the net exchange flux through the static-moving bed layers interface and provide new expressions for the erosion and deposition rates ${\dot{\eta }}_e$ and ${\dot{\eta }}_d$, respectively. The proposed bedload transport model is supported by a generalized Grass-type formula for the non-equilibrium 2D bedload rate ${\mathbf{q}}_{\mathbf{b}}$ as a function of the transport layer thickness $\eta$, which must be solved depending on equilibrium or non-equilibrium approach is considered. The 2D morpho-dynamical system of flow-bed conservation equations is solved using and fully-coupled augmented Roe solver [34] whereas the bedload transport layer equation is solved using an upwind scheme. Furthermore, the maximum bed slope is limited using an efficient non-iterative morphological collapse mechanism, which imposes geometrical restrictions on each cell edge.

The proposed model has been benchmarked against dyke erosion and dambreach opening laboratory experiments, involving highly unsteady flows. The non-equilibrium model is able to predict the measured data in all the cases tested, improving the results obtained with the classical equilibrium approach. Furthermore, the effects of the entrainment and deposition constants, $k_e$ and $k_d$ respectively, have also been analysed. When non-equilibrium model is applied, the adaptation of the actual bedload solid rate to the flow carry capacity becomes non-instantaneous and the appearance of bed changes is delayed in time and space. The temporal and spatial delay of the non-equilibrium bedload rate respect to the corresponding flow carry capacity increases as the entrainment constant $k_e$ decreases. Moreover, as the local Shield stress $\theta$ exerted by the flow in the bedload transport layer grows, the equivalent dynamic adaptation length $L_b$ is increased and hence the non-equilibrium states are promoted (see Section 4.1).

In general, considering the unsteady dambreach flow configuration, better predictions of the breach opening mechanism can be obtained using the non-equilibrium approach for the bedload transport rate with moderate values of the exchange entrainment and deposition parameters respect those obtained using the equilibrium model. Numerical results seems to support that considering non-equilibrium conditions for the bedload transport improves the general agreement between the computed results and measured data for overtopping dyke erosion and dambreach opening cases, delaying the breach opening process. However, the non-equilibrium approach requires a careful calibration of the exchange parameters in regions with complex transient processes.