Numerical simulation method for predicting a flood hydrograph due to progressive failure of a landslide dam

Reducing the damage due to landslide dam failures requires the prediction of flood hydrographs. Although progressive failure is one of the main failure modes of landslide dams, no prediction method is available. This study develops a method for predicting progressive failure. The proposed method consists of the progressive failure model and overtopping erosion model. The progressive failure model can reproduce the collapse progression from a dam toe to predict the longitudinal dam shape and reservoir water level when the reservoir water overflows. The overtopping erosion model uses these predicted values as the new initial conditions and reproduces the dam erosion processes due to an overtopping flow in order to predict a flood hydrograph after the reservoir water overflows. The progressive failure model includes physical models representing the intermittent collapse of a dam slope, seepage flow in a dam, and surface flow on a dam slope. The intermittent collapse model characterizes the progressive failure model. It considers a stabilization effect whereby collapse deposits support a steep slope. This effect decreases as the collapse deposits are transported downstream. Such a consideration allows the model to express intermittent, not continuous, occurrences of collapses. Field experiments on the progressive failure of a landslide dam were conducted to validate the proposed method. The progressive failure model successfully reproduced the experimental results of the collapse progression from the dam toe. Using the value predicted by the progressive failure model, the overtopping erosion model successfully reproduced the flood hydrograph after the reservoir water started to overflow.


Introduction
Heavy rainfall and strong earthquakes can cause large-scale mass movements such as landslides and debris flows in mountain basins. Such mass movements can block river valleys, forming landslide dams (Costa and Schuster 1988;Fan et al. 2020). Most landslide dams fail within a year (Peng and Zhang 2012). The landslide dam failure triggers an enormous flood flow, causing catastrophic damage downstream. Two successive landslides occurred in Baige village, the border between Sichuan Province and Tibet Autonomous Region, in China, on October 10 and November 3, 2018 (Zhong et al. 2020). These landslides completely blocked the Jinsha River on both occasions, resulting in landslide dam failure. The failures of the landslide dams formed on October 10 and November 3 produced flood flows of approximately 10,000 and 31,000 m 3 /s peak discharge, respectively.
Reducing the damage due to a landslide dam failure requires the prediction of the flood flow scale. To this end, prediction methods, including statistical and numerical analysis methods, have been developed. Studies on statistical analysis methods have investigated correlations between parameters, such as dam height, reservoir water volume, and peak flood discharge, and relevant regression equations have been proposed (Costa 1985;Walder and O'Connor 1997;Peng and Zhang 2012). Studies on numerical analysis methods have clarified dam failure processes by conducting experiments or observations, and numerical simulation methods that can reproduce the relevant processes have been developed (e.g., Takahashi and Nakagawa 1994;Chang and Zhang 2010;Cao et al. 2011). Numerical simulation methods can predict a flood hydrograph caused by a landslide dam failure, unlike regression equations, which can only predict the peak flood discharge. A flood hydrograph is necessary for evaluating inundation areas with high accuracy.
The development of numerical simulation methods requires sufficient understanding of landslide dam failure processes. The modes of landslide dam failure include overtopping failure, slope failure, and progressive failure (Takahashi 2014). In overtopping failure, an overtopping flow erodes a dam after the reservoir water overflows. In slope failure, the sliding of the upper part of a dam induces the release of the reservoir water. In progressive failure, a partial collapse at the dam toe proceeds upward, leading to the collapse of the entire dam. Global landslide dam databases show that overtopping, progressive, and slope failures account for 91%, 8%, and 1%, respectively, of all dam failures (Peng and Zhang 2012).
Previous studies have proposed numerical simulation methods for overtopping and slope failures of a landslide dam. Numerical simulation methods for overtopping failure involve a surface flow model and side bank collapse model, which successfully reproduce the experimental or observed flood hydrographs (Takahashi and Nakagawa 1994;Chang and Zhang 2010;Satofuka et al. 2010;Cao et al. 2011;Akazawa et al. 2014;Zhong et al. 2018). Numerical simulation methods for slope failure involve a seepage flow model and dam slope stability model, which predict slope failure occurrence (Awal et al. 2007). Combining slope failure methods with overtopping failure methods enables the prediction of a flood hydrograph due to slope failure (Awal et al. 2008).
Experimental studies have proposed equations that can predict whether progressive failure will occur or not (Gregoretti et al. 2010;Jiang et al. 2020). Gregoretti et al. (2010) showed that a downstream dam slope gradient, channel gradients, and the median particle diameter of a landslide dam are the main parameters for predicting the minimum reservoir water level that triggers the progressive failure, and they proposed a regression equation using these dimensionless parameters. Jiang et al. (2020) indicated that both seepage flow stress and surface flow stress affect progressive failure occurrence and proposed a semi-theoretical equation by considering both effects.
Previously, we developed a numerical simulation method for progressive failure, which can predict a flood hydrograph (Takayama et al. 2020). This method includes a seepage flow model, dam slope stability model, and surface flow model, which successfully reproduced the experimental flood hydrographs obtained from small-scale flume experiments. These experiments created 0.15-m-high landslide dams, which are much smaller than real landslide dams. Validating the method using field-scale data is essential for predicting flood flows due to real landslide dam failures.
This study aims to develop a numerical simulation method for the progressive failure of a landslide dam, which can reproduce flood hydrographs obtained from field-scale experiments. Thus, experiments are conducted by considering the progressive failure of a landslide dam in a mountain stream. The experimental results are subsequently used to validate the proposed numerical simulation method.
Field experiments on progressive failure of a landslide dam Experimental setup and procedure The field experiments were conducted under two conditions (cases 1 and 2) with differences in the longitudinal dam shapes in the Hirudani stream of the Ashiaraidani Basin, Takayama City, Gifu Prefecture, Japan (Supplementary Fig. 1 and 2). Figure 1 shows schematic diagrams of the landslide dam shapes. The longitudinal dam shapes of both cases 1 and 2 were trapezoids (dam height, 1.0 m; upstream dam slope, 29°, downstream dam slope, 27°). The dam crest length (L) of case 1 and case 2 was 0.8 m and 1.4 m, respectively. The longitudinal dam shape of case 2 was thus longer in the flow direction than that of in case 1. The landslide dams were constructed on highly permeable layers (Fig. 1). The layer thickness of case 1 and case 2 was 0.2 m and 0.3 m, respectively. The inflow discharge into the dam in case 1 and case 2 was 0.052 m 3 /s and 0.060 m 3 /s, respectively. Peng and Zhang (2012) suggested that the landslide dam height (H), longitudinal dam bottom length (L b ), dam volume (V d ), and reservoir volume (V w ) are the key geometric parameters of landslide dams, and they proposed the following set of dimensionless w =H. These parameters allow verification of whether an experimental model concurs with real topography including a landslide dam (Zhou et al. 2019). In case 1, H/L b was 0.20, V 1=3 d =H was 2.0, and V 1=3 w =H was 2.2. In case 2, H/L b was 0.18, V 1=3 d =H was 2.2, and V 1=3 w =H was 2.3. These values are within the acceptable ranges of the dimensionless parameters of real landslide dams, as shown by Zhou et al. (2019). Thus, the current experimental models represent real topography including a landslide dam.
Laboratory soil tests of the dam and highly permeable layer were conducted. Figure 2a shows the results of the particle size distribution tests. The dam and highly permeable layer were composed of homogeneous cohesionless soils with mean particle diameters of 5 and 14 mm, respectively. The constant head permeability test of the dam material showed that its saturated hydraulic conductivity was 0.001 m/s. Figure 2b shows the results of the water retention test of the dam material and the fitting result of the water retention curve function proposed by Tani (1982). An amount of water corresponding to 10% of the water content of the dam was added into its material while constructing the dam. Casagli et al. (2003) reported the particle size distributions of 42 landslide dams in the Northern Apennines in Italy, the range of which is shown in Fig. 2a. Most of the particle size distributions were bimodal, which suggests that real landslide dams are composed of heterogeneous soils. In our experiments, the dam and highly permeable layer were composed of homogeneous soils without fine particles, in contrast to real landslide dams (Fig.  2a). Removing fine particles from a landslide dam increases its permeability and decreases its apparent cohesion due to matric suction. Such a discrepancy between our experimental models and real landslide dams was necessary to guarantee the occurrence of progressive failure.
Our experiments measured the inflow discharge into a dam, dam failure process, reservoir water level, groundwater level in a dam, and outflow discharge from a dam. An overtopping weir approximately 50 m upstream of the dam construction site measured the inflow discharge. Video cameras around the landslide dam recorded the dam failure process. These videos were used for structure from motion (SfM) analysis to obtain three-dimensional coordinates of the dam shapes. Pressure water level gauges (S&DL mini; OYO Corporation, Japan) measured the reservoir and groundwater levels. These gauges were installed at the bottom of the reservoir and dam (Fig. 1). The outflow discharge from a dam was calculated by subtracting the temporal change in the reservoir water volume from the inflow discharge. The reservoir water volume was given on the basis of the measured reservoir water level and river topography.
Experimental results of the dam failure processes Figure 3 shows the experimental results of the dam failure processes of case 1. Before 78 s, although seepage from the dam toe continued to occur, no collapse was observed visually. At 78 s, a small-scale collapse occurred at the center of the dam toe. From 78 to 366 s, the downstream dam slope collapsed intermittently, resulting in collapse progression from the dam toe. During the collapse progression, the surface flow continuously transported the collapse deposits downstream, inducing intermittent collapse. The reservoir water overflowed at 366 s, and the overtopping flow subsequently eroded the dam rapidly, triggering a large-scale flood flow.
The dam failure process and mechanism of case 2 were similar to those of case 1, except for the starting times of collapse progression and reservoir overflow. The starting times of case 2 were later than those of case 1 because the longitudinal dam shape of case 2 was longer in the flow direction than that of case 1.
The partial collapse at the center of the dam toe extended not only in the flow direction but also in the cross-flow direction. This extension is different from the results of experiments using a flume with a rectangular cross section (Takayama et al. 2020). In the flume experiments, the collapse at the dam toe extended only in the flow direction. Such a difference may arise because of the difference between the channel cross-sectional geometries of the field and flume experiments.

Discussion on the dam failure processes
In the present experiments, the collapse progression from the dam toe induced reservoir overflow. Subsequently, the overtopping flow eroded the dam rapidly. Such processes can be divided into progressive and overtopping erosion processes (Fig. 4). The progressive failure process is defined as the process in which the collapse at a dam toe proceeds upward until the reservoir water overflows.

Technical Note
The overtopping erosion process is defined as the process in which an overtopping flow erodes a dam after the reservoir water overflows.
The progressive failure process includes three physical processes: intermittent collapse of a dam slope, seepage flow in a dam, and surface flow on a dam slope (Fig. 4a). Reservoir water level fluctuations cause unsteady seepage flow in a dam, generating surface flow on a dam slope. Both the seepage and the surface flow induce the intermittent collapse of a dam slope.
The overtopping erosion process is dominated by surface flow on a dam slope (Fig. 4b). The surface flow is larger and more unsteady than that of the progressive failure process. The seepage flow in a dam does not affect the overtopping erosion process significantly because the downstream dam slope surface is generally saturated when the reservoir water overflows (case 1 at 366 s, as shown in Fig. 3).

Flowchart of the numerical simulation method
The numerical simulation method proposed in this study involves a progressive failure model and an overtopping erosion model. The progressive failure model considers the intermittent collapse of a dam slope, seepage flow in a dam, and surface flow on a dam slope. This model expresses the progressive failure process (Fig. 4a). The overtopping erosion model considers only the surface flow on a dam slope. It expresses the overtopping erosion process (Fig. 4b).  (1) The initial and boundary conditions, including the initial river and dam geometries and inflow discharge into a dam, are set. (2) The progressive failure model reproduces the progressive failure process to predict the longitudinal dam shape and reservoir water level when the reservoir water overflows. (3) The overtopping erosion model uses these predicted values as the new initial conditions and reproduces the overtopping erosion process to predict a flood hydrograph after the reservoir water overflows.
Progressive failure model The progressive failure model considers the intermittent collapse of a dam slope, seepage flow in a dam, and surface flow on a dam slope. In the intermittent collapse, the two-dimensional simplified Janbu's method predicts the collapse occurrences, and the equilibrium sediment transport equation represents the transport of collapse deposits downstream. The seepage flow is described by the two-dimensional Richards' equation. The surface flow is represented by a one-dimensional conventional simulation model of a debris flow. The intermittent collapse model characterizes the progressive failure model. This collapse model considers a stabilization effect whereby collapse deposits support a steep slope. This effect decreases as the collapse deposits are transported downstream. Such a consideration allows the model to express intermittent, not continuous, occurrences of collapses.

Intermittent collapse of a dam slope
In the progressive failure process, both the seepage and the surface flow induce the intermittent collapse of a dam slope. This collapse mechanism is similar to riverbank erosion mechanisms. Riverbank erosion is due to both fluvial erosion and mass failure. Fluvial erosion changes the bank geometry and destabilizes it. This destabilization induces mass failure and cooperated with increasing the pore water pressure, decreasing the matric suction, and losing the hydrostatic confining pressure (Rinaldi et al. 2004(Rinaldi et al. , 2008Rinaldi and Casagli 1999;Simon et al. 2000).
Some numerical simulation methods considering fluvial erosion and mass failure of a riverbank have been proposed . The fluvial erosion rate is typically evaluated by assuming that the excess shear stress acting on a riverbank is proportional to the erosion rate (Partheniades 1965;Arulanandan et al. 1980). The mass failure occurrence is typically predicted by limit equilibrium methods (LEMs), which calculate a safety factor (F s ) defined as the ratio of stabilizing forces to destabilizing forces. The LEMs can consider the effect of the differences in bank geometry on stability. This consideration allows the expression of destabilization due to changes in the bank geometry caused by fluvial erosion (Osman and Thorne 1988). The LEMs can also consider the effects of pore water pressure, matric suction, and hydrostatic confining pressure on bank stability (e.g., Amiri-Tokaldany et al. 2003;Thorne 1996, Darby et al. 2007;Rinaldi and Casagli 1999;Rinaldi et al. 2004Rinaldi et al. , 2008Simon et al. 2000).  Collapse deposits support a steep riverbank slope; however, no method considers this stabilization effect. Bank erosion analysis may ignore this effect because hydrostatic confining pressure strongly affects riverbank stability (Rinaldi et al. 2004(Rinaldi et al. , 2008. In the progressive failure process of a landslide dam, the earth pressure from collapse deposits can primarily affect the slope stability. This stabilization effect decreases as the collapse deposits are transported downstream, triggering the collapse of a dam slope. This collapse produces deposits to recover the stabilization effect. Such a cycle of the stabilization effect can have a major impact on the progressive failure process. This study proposes an intermittent collapse model by considering the stabilization effect whereby collapse deposits support a steep slope (Fig. 6). The LEM used here considers the earth pressure from collapse deposits (P p ), representing the stabilization effect. The collapse deposits are assumed to be transported by the equilibrium sediment discharge rates of the surface flow. This assumption enables the expression of the gradual loss of the stabilization effect.
Equations of slope stability Slice methods, which are representative LEMs, divide a slip mass into multiple slices (Fig. 6a) and then calculate a safety factor (F s ) that can determine whether a collapse occurs (F s < 1.0) or not (F s ≥ 1.0) using force and moment equilibrium equations and safety factor equations for each slice. Slice methods can physically consider the spatial differences in soil strength and external forces. However, these methods involve assumptions because the number of equations is smaller than the number of unknown variables. Various assumptions generate various slice methods (Duncan and Wright 1980).
The simplified Janbu's method is a typical slice method applicable to any slip surface. This method discards the moment equilibrium equations and is thus less rigorous than the methods in which all equilibrium equations are satisfied (e.g., rigorous Janbu's method, Morgenstern and Price's method, and Spencer's method). However, the simplified Janbu's method has high numerical stability. This is necessary for calculating the safety factors (F s ) multiple times under various conditions. The simplified Janbu's method has been applied to many slope stability problems, and it has successfully reproduced the collapse occurrence time and slip surface shape (e.g., Awal et al. 2007;Tsutsumi and Fujita 2008).
This study used the simplified Janbu's method. The earth pressure from collapse deposits (P p ) was introduced into the commonly used slope stability equations to consider the stabilization effect provided by collapse deposits supporting a steep slope. The earth pressure from the collapse deposits (P p ) is assumed to act on the slice side in the horizontal direction (Fig. 6a). The slope stability equations are thus given by where F s is the safety factor, R i is the resistance component, T i is the sliding component, P p is the earth pressure from collapse deposits, c is the cohesive strength, τ suc i is the apparent cohesive strength due to matric suction, l i is the length of the slice base, α s i is the slice base gradient, W i is the slice weight, u i is the pore water pressure acting on the slice base, and ϕ is the internal friction angle. Equation (1) shows that the larger the value of P p , the larger is the value of F s .
The earth pressure from collapse deposits (P p ) is calculated using Coulomb's widely accepted earth pressure theory. A slice is acted by passive, but not active, earth pressure from collapse deposits because the slip failure occurrence results in lifting of the collapse deposits. The earth pressure from collapse deposits (P p ) is thus given by Coulomb's earth passive pressure formula as where σ is the soil density, ρ is the water density, C * is the sediment concentration of deposits, θ c is the volumetric water content of deposits, g is the gravitational acceleration, H c is the collapse deposit thickness (Fig. 6a), K p is Coulomb's passive earth pressure coefficient, and ζ is the deposition gradient of the collapse deposits. In this study, ζ was assumed to be equal to the internal friction angle (ϕ) because the collapse deposits in the present experiment are cohesionless. Further, H c is calculated geometrically using the volume of collapse deposits (V c ) obtained from Eq. (4).
The apparent cohesive strength due to matric suction (τ suc i ) and the pore water pressure acting on the slice base (u i ) are calculated using the profiles of the pressure head (ψ) and volumetric water content (θ), which are calculated using seepage flow models described later. Some methods for evaluating τ suc i have been proposed (e.g., Fredlund et al. 1978;Vanapalli et al. 1996), and the empirical analytical model proposed by Vanapalli et al. (1996) has been introduced into riverbank stability analyses (Arai et al. 2018;Mizutani et al. 2013). In the present study, the following model was used: where ψ sb i is the average pressure head along the slice base, θ sb i is the average volumetric water content along the slice base, θ s is the saturated volumetric water content, and θ r is the residual water content.

Searching for a slip surface
To determine the slip surface shape, it is necessary to search for the slip surface having the minimum safety factor. To this end, several algorithms have been proposed (e.g., Baker 1980;Yamagami and Ueta 1986). Slopes at water edges, such as riverbank slopes, tend to fail along a planar slip surface (Osman and Thorne 1988); therefore, their stabilities are typically analyzed by limiting the analysis to only planar slip failures. Despite this limitation, collapse occurrences at water edges are successfully reproduced (e.g., Rinaldi and Casagli 1999;Simon et al. 2000). The present study thus limits the search for slip surfaces to only planar slip surfaces having a scarp (Fig. 6a).
The procedure for the slope stability analysis is as follows. (I) The safety factors of every conceivable planar slip surface are calculated using Eq. (1), and the slip surface having the minimum safety factor is identified. (II) If the minimum safety factor is less than 1.0 (F s < 1.0), it is determined that a slip failure occurs along the corresponding slip surface. If the minimum safety factor is greater than 1.0 (F s ≥ 1.0), it is determined that a slip failure does not occur; thus, the slope stability analysis procedure is completed. (III) If F s < 1.0, the dam surface height decreases to the slip surface height. The volume (V c ) and volumetric water content (θ c ) of the collapse deposits are also updated using Eqs. (4) and (5), respectively.

Equations for fluidization of collapse deposits
The volume of the collapse deposits continuously decreases as they are transported downstream by a surface flow and intermittently increases with the occurrence of slip failures. The mass conservation equation is thus expressed as where V c is the volume of the collapse deposits, t is the time, q dps s is the sediment discharge from the collapse deposits, and V slip is the volume of a new sliding mass.
The volumetric water content of the collapse deposits (θ c ) varies spatially; however, it is assumed to be spatially uniform for simplicity. Owing to this assumption, the mass conservation of water in the collapse deposits can be expressed as where θ slip is the volumetric water content in a new sliding mass.
The surface flow in the area just below the scarp (called the fluidization area, as shown in Fig. 6b) continuously transports the collapse deposits downstream. The sediment discharge from the collapse deposits (q dps s ) is given by the equilibrium sediment discharge of the surface flow in the fluidization area as follows: where q dps s is the sediment discharge from the collapse deposits, q fld is the surface flow discharge in the fluidization area, C ∞fld is the equilibrium sediment concentration of the surface flow in the fluidization area, θ w fld is the average water surface gradient in the fluidization area, h fld is the surface flow depth in the fluidization area, and C fld is the sediment concentration of the surface flow in the fluidization area. Further, q fld and C ∞fld are  12) and (15), respectively, using θ w fld , h fld , and C fld . In addition, h fld and C fld are assumed to be spatially uniform and are obtained from Eq. (10).
The length of the fluidization area (L fld ) may be correlated with the scarp height (h scarp ) shown in Fig. 6b. Thus, L fld is given by L fld = k fld h scarp , where k fld is a coefficient. In this study, the value of k fld was set to 3.0 for convenience.

Seepage flow in a dam
Richards' equation, which describes the flow in both saturated and unsaturated porous media, is used in hydrosciences such as forest hydrology and hydraulic engineering. It has successfully reproduced the seepage flow in riverbanks and landslide dams (e.g., Awal et al. 2007;Rinaldi et al. 2008). However, it is one of the most challenging equations to solve reliably and accurately among all the hydrosciences (Farthing and Ogden 2017).
Richards' equation solvers commonly use a finite difference method (FDM), finite volume method (FVM), or finite element method (FEM) for spatial discretization (Farthing and Ogden 2017). The FVM and FEM, which use polygons for calculation grids (i.e., unstructured grids), can accurately express complex boundary shapes. The FDM, which is limited to rectangular calculation grids (i.e., structured grids), requires a grid size that is sufficiently small to express boundary shapes, but its grid partitioning is relatively easy. This simplicity facilitates the solution of problems that need to adjust the grid partitioning to fit changing boundary shapes over time, as in the present problem.
Richards' equation solvers typically use implicit methods, represented by modified Picard iterations, for temporal discretization (Celia et al. 1990;Farthing and Ogden 2017). These implicit solutions, which require iterative computations, often converge rapidly and easily; however, some solutions do not converge for a variety of reasons (Niswonger and Prudic 2009). Satofuka and Mizuyama (2009) proposed a numerical simulation model representing the development of debris flows on unsaturated deposits. This model coupled Richards' equation, which describes the seepage flow in deposits, with the shallow water equation, which describes surface flows on deposits. Both equation solvers used the same spatial and temporal discretization methods to facilitate the coupling of the two equations. The spatial and temporal discretizations used an FDM and an explicit method, respectively. The explicit method does not require iterative computations and thus provides a stable solution as long as the time step is sufficiently small to satisfy the Courant condition.
Based on the work of Satofuka and Mizuyama (2009), Richards' equation solver in this study adopted an FDM using a staggered grid for spatial discretization and an explicit method for temporal discretization. In the staggered grid, the definition positions of vector quantities are shifted by half the grid size from those of scalar quantities (Fig. 7a), increasing the calculation stability. The same methods for discretization are used in the surface flow model, as described later.

Equation governing the seepage flow
For the vertical two-dimensional field shown in Fig. 7a, Richards' equation describing a seepage flow in a dam is expressed as where θ is the volumetric water content, ψ is the pressure head, β is a coefficient having a value of 1.0 or 0.0 in a saturated or unsaturated zone, respectively, S s is the specific storage, M is the seepage flow flux in the x-direction, N is the seepage flux in the zdirection, and K is the hydraulic conductivity.
The solver for Eq. (7) adopted an FDM using the staggered scheme for spatial discretization and an explicit method for temporal discretization. The solution procedure is as follows. (I) The flow field (M and N) of the next time step is calculated using the pressure field (ψ) of the previous time step. (II) The pressure field of the next time step is calculated using the flow field calculated above.
Function of water retention curve and unsaturated hydraulic conductivity Solving Eq. (7) requires the functions for the relation between θ-ψ and K-θ. Many such functions have been proposed, and their validity has been verified by comparing them with soil test data (Leij et al. 1997). The function for the θ-ψ relation proposed by Tani (1982) is a simple function; it has only a few parameters but accurately represents the θ-ψ curves of three typical soil types, as shown by Hillel and van Bavel (1976). This function also sufficiently represents the results of the soil water retention test of the landslide dam material used in our experiments, as shown in Fig.  2b. The function for the K-θ relation proposed by Brooks and Corey (1964), which similarly has only a few parameters, also represents the soil test data relatively well (Leij et al. 1997).
Thus, this study used the function for the θ-ψ relation proposed by Tani (1982) and that for the K-θ relation proposed by Brooks and Corey (1964) as follows: where θ s is the saturated volumetric water content, θ r is the residual volumetric water content, ψ 0 is the pressure head at the inflection point of the soil water retention curve, K s is the saturated hydraulic conductivity, and m is a coefficient representing the K-θ relation.

Boundary condition of a dam surface
Two types of methods for setting the boundary condition between a surface and a seepage flow have been proposed: one uses the water exchange flux equation, and the other assumes pressure continuity (Kollet and Maxwell 2006). In the former method, Darcy's law gives the flux across the boundary using the hydraulic gradient calculated from the surface flow depth and seepage flow pressure head (Satofuka and Mizuyama 2009;Ogasawara and Sekine 2010). In the latter method, the surface and seepage flow pressures acting on the boundary are assumed to be equal constantly (Kollet and Maxwell 2006). The pressure continuity method cannot express unsteady phenomena, such as surface flow on unsaturated deposits. By contrast, the water exchange equation method is applicable to such phenomena.
This study used the water exchange flux equation proposed by Satofuka and Mizuyama (2009). Darcy's law provides the water exchange flux across the dam surface (w i in Fig. 7b) using the hydraulic gradient calculated from the depth of a surface flow on a dam slope (h i ) and the pressure head of a seepage flow near a dam surface (ψ i, jb ) as follows: where w i is the water exchange flux, h i is the depth of a surface flow on a dam slope, ψ i, jb is the pressure head of a seepage flow near a dam surface, and Δz ′ is the distance between the calculation point of ψ i, jb and the dam surface (Fig. 7b).

Surface flow on a dam slope
In the progressive failure process, a surface flow transports the collapse deposits, produced by the intermittent collapse of a dam slope, to a place where the riverbed gradient is sufficiently gentle. In such sediment transportation, riverbed gradients vary in time and space, triggering various sediment transport modes, including debris flow and bedload.
Many numerical simulation models for debris flow have been proposed (e.g., Takahashi et al. 1992;Iverson 1997;Egashira et al. 2001;Liu and Huang 2006;Rickenmann et al. 2006;Wu et al. 2013;Han et al. 2015). These models commonly include flow continuity equations for the total volume, momentum equations, and constitutive equations. Considering debris flow erosion/deposition also requires flow continuity equations for the sediments, riverbed continuity equations, and erosion/deposition velocity equations (Takahashi et al. 1992;Egashira et al. 2001). Various constitutive equations and erosion/deposition velocity equations have been proposed; however, they have not been interpreted uniformly among researchers, whereas the other equations are generally common to all the models.
The numerical simulation model proposed by Takahashi et al. (1992) (Takahashi model) considers various sediment transportation modes, including stony debris flow, immature debris flow, and bedload. Stony debris flow occurs with gravel dispersed throughout the flow on a steep riverbed, and this flow is assumed to be a dilatant fluid that depends on the particle collision stress (Takahashi 1978). Immature debris flow occurs on a gentler riverbed because the gravel cannot be dispersed throughout the flow; there is a two-layer flow, with the lower layer containing the gravel and water and the upper layer containing only water. Bedload occurs on a much gentler riverbed because the flow cannot transport the mass of gravel. This results in an ordinary turbulent flow.
The Takahashi model is widely accepted in Japan. On the basis of this model, Nakatani et al. (2008) developed a two-dimensional debris flow simulator with a graphical user interface (GUI). The simulator has been improved for practical purposes and has successfully reproduced the observed data of debris flow disasters (Nakatani et al. 2012;Liu et al. 2013;Uchida et al. 2013;Nakatani et al. 2016aNakatani et al. , 2016b. Thus, the Takahashi model has high reliability. This study used the Takahashi model to represent a surface flow on a dam slope. The solver of this model adopted an FDM using the staggered scheme for spatial discretization and an explicit method for temporal discretization, as with the debris flow simulator developed by Nakatani et al. (2008). This solver is the same as the solver of Richards' equation (Eq. (7)), facilitating the coupling of the surface and seepage flow models.

Flow continuity equations
The flow continuity equations for the total volume and sediments are given by adding the water exchange flux through a dam surface (w i ) calculated from Eq. (9) to the equations generally used in the Takahashi model as where h is the surface flow depth, q is the surface flow discharge, w i is the water exchange flux through a dam surface, i b is the erosion/deposition velocity (erosion is positive), C * is the sediment concentration of deposits, θ bed is the volumetric water content of a riverbed, θ s is the saturated water content, and C is the sediment concentration of a surface flow. Further, w i enables coupling of the seepage and surface flow models.
The riverbed continuity equation is generally expressed as (Takahashi et al. 1992) where z b is the riverbed height.

Momentum and constitutive equations
The uniform flow approximation of the momentum conservation laws derives the uniform flow equation as where R is the flow resistance coefficient. The Takahashi model gives the equations representing the flow resistance coefficient (R) of stony debris flow, immature debris flow, and ordinary turbulent flow. These sediment transport modes change with the differences in the sediment concentration of a surface flow (C) and the relative flow depth (h/d, where d is the sediment particle diameter) (Takahashi et al. 1992;Takahashi and Nakagawa 1994;Uchida et al. 2013). These critical values were empirically determined on the basis of the experimental results. If C ≥ 0.4C * , stony debris flow occurs, where R is given by the theoretical equation derived from assuming the dilatant fluid. If 0.01 < C < 0.4C * , immature debris flow occurs, where R is given by the empirical equation obtained from the experimental results. If 0.01 ≤ C or h/d ≥ 30, ordinary turbulent flow occurs, where R is given by Manning's equation. The flow resistance coefficient equations are thus expressed as where d is the sediment particle diameter and n m is Manning's coefficient.

Erosion/deposition velocity equations
The Takahashi model expresses the erosion/deposition velocity of a debris flow in a form proportional to the difference between the equilibrium sediment concentration (C ∞ ) and the sediment concentration (C) of a surface flow (Takahashi et al. 1992). The Takahashi model has two types of erosion velocity equations for saturated and unsaturated beds. This study used only the erosion velocity equation for a saturated bed because a saturated bed is common in the progressive failure process. The erosion/ deposition velocity equation is given by where C ∞ is the equilibrium sediment concentration of a surface flow, δ e is the erosion velocity coefficient, and δ d is the deposition velocity coefficient.
The Takahashi model gives the equilibrium sediment concentration (C ∞ ) equations for stony debris flow, immature debris flow, and bedload. These sediment transportation modes are distinguished by the difference in the water surface gradients (θ w ) (Takahashi et al. 1992;Takahashi and Nakagawa 1994;Uchida et al. 2013). If tanϕ > tan θ w > 0.138, stony debris flow occurs, where C ∞ is given by a theoretical equation derived from the equilibrium of forces acting on the riverbed. If 0.138 ≥ tan θ w > 0.03, immature debris flow occurs, where C ∞ is given by an empirical equation obtained from the experimental results. If tanθ w ≤ 0.03, bedload occurs, where C ∞ is given by a theoretical equation derived from the bedload discharge formula. The critical value of 0.138 is decided to become the continuous transition of the C ∞ between the stony and immature debris flow. The critical value of 0.03 was empirically determined. The equilibrium sediment concentration equations are thus given by where τ * is the dimensionless bed shear stress and τ * c is the critical dimensionless bed shear stress. The maximum value of C ∞ is known to be 0.9C * based on t experimental results (Takahashi 2014). If the calculated value of C ∞ is greater than 0.9C * , the value is modified to 0.9C * .

Critical condition of sediment movement
Assuming a deposit layer composed of particles having a uniform grain size, a surface layer that is shallower compared to the grain size cannot start to move unless the particles are fractured, which is not usually considered. Starting the sediment movement requires the shear stress acting on the face located one grain size below the deposit surface (τ d ) to be greater than its shear resistance (τ dc ). The critical condition of sediment movement is given by (Takahashi 2014) where τ d and τ dc are the shear stress and shear resistance, respectively, acting on the face located one grain size below the deposit surface. If this critical condition is satisfied, no debris flow erosion occurs (i b ≤ 0).

Overtopping erosion model
In the overtopping erosion process (Fig. 4b), the reservoir water overflows from the lowest point of the dam crest, and the overtopping flow then erodes the dam, forming a channel on its surface. The channel is vertically eroded owing to overtopping flow (i.e., vertical erosion) and is horizontally expanded owing to side bank erosion or collapse (i.e., lateral erosion). The surface flow in the channel is larger and more unsteady than that in the progressive failure process (Fig. 4a).
Many numerical simulation models representing the overtopping erosion processes have been proposed, and they have successfully reproduced experimental or observed flood hydrographs (e.g., Takahashi and Nakagawa 1994;Chang and Zhang 2010;Satofuka et al. 2010;Cao et al. 2011). Most of the studies have introduced models representing lateral erosion in conventional surface flow models representing vertical erosion.
In this study, the equations representing surface flow in the overtopping erosion process are the same as those in the progressive failure process (Eqs. (10), (11), (12), (13), (14), (15), and (16)), except for the momentum equation. The overtopping erosion model uses the momentum equation including the acceleration term, instead of the uniform flow equation (Eq. (12)) to consider a highly unsteady surface flow. The momentum equation including the acceleration term is expressed as where u is the surface flow velocity, z b is the riverbed height, h is the surface flow depth, τ b is the shear stress acting on the riverbed, and ρ m is the surface flow density. Most numerical simulation models for the overtopping erosion process have also used the momentum equation including the acceleration term (e.g., Takahashi and Nakagawa 1994;Satofuka et al. 2010;Cao et al. 2011).
The Takahashi model gives the shear stress acting on a riverbed (τ b ) of stony debris flow, immature debris flow, and bedload, similar to the resistance coefficient equations of the progressive failure model (Eq. (13)). These sediment transportation modes were distinguished in the same way as Eq. (13). The equations of shear stress acting on a riverbed are thus given by where d is the sediment particle diameter, C is the sediment concentration of the surface flow, ρ is the water density, σ is the soil density, C * is the sediment concentration of deposits, and n m is Manning's coefficient.
In this study, no horizontal erosion model was used for simplicity. The overtopping flow channel width was set to a constant value of 1.0 m by referring to the present experimental videos.

Numerical simulation conditions and procedures
Numerical simulation conditions The numerical simulation method proposed herein was validated by reproducing the experimental results of cases 1 and 2. Figure 8 shows the initial longitudinal topography. The initial topography of case 1 used the topographic measurements obtained via the total station and SfM, whereas that of case 2 used the topographic measurement obtained via the total station only because the SfM for case 2 has low accuracy.
The actual dam cross section has a bottom width of around 2.0 m and a top width of 4.0 m ( Fig. 1c and Supplementary Fig. 3a). In the progressive failure model, the dam cross section has a uniform width of 3.0 m ( Supplementary Fig. 3b) because this model is onedimensional. In the actual overtopping erosion process, the overtopping flow width (B f ) gradually increased ( Supplementary  Fig. 3c). For simplicity, the model assumes that B f is constantly 1.0 m by referring to the present experimental videos ( Supplementary  Fig. 3d).
The relationship between the reservoir water volume and level was obtained from the topographic measurement via the total station ( Supplementary Fig. 4). This relationship was used to transform the calculated reservoir water volume into the level. Table 1 lists the parameters for the numerical simulation. The representative sediment diameter (d) used the mean diameter of the particle size distribution test results shown in Fig. 2a. The parameters of the soil water retention curve function (Eq. (8)) were determined by fitting the water retention test results of the dam material, as shown in Fig. 2b. The saturated hydraulic conductivity (K s ) used the value that reproduced the experimental results most accurately. Further, K s for the dam material and highly permeable layer material was set to 0.030 m/s and 0.150 m/s, respectively. However, the constant head permeability test of the dam material showed that its K s was 0.001 m/s, which is much smaller than our setting. This may be because a partial fast flow through a dam actually exists. The erosion velocity coefficient (δ e ), including Eq. (14), used the value commonly used in the progressive failure model and the value that reproduced the experimental results most accurately in the overtopping erosion model. The other parameters used general values.

Numerical simulation procedures
The procedure of the numerical simulation was as follows: 1. The steady-state pressure head (ψ) profile of the dam was precomputed. First, the pressure head in the dam was uniformly given by ψ = − 0.051 m, corresponding to the experimental value. Subsequently, the progressive failure model was analyzed for 10,800 s under the constant reservoir water level observed before starting the experiments (0.00 and 0.10 m in cases 1 and 2, respectively). 2. The precomputed results were used as the initial condition.
The progressive failure model reproduces the progressive failure process (Fig. 4a) under the inflow condition observed in the experiments (0.052 and 0.060 m 3 /s in cases 1 and 2, respectively), predicting the riverbed height (z b ) and reservoir water level when the reservoir water overflows. 3. The new initial conditions used the values predicted by the progressive failure model. The overtopping erosion model reproduced the overtopping erosion process (Fig. 4b) under the inflow condition observed in the experiments, predicting the flood hydrograph after the reservoir water overflows.>

Verification of the validity of the numerical simulation model
Progressive failure processes of a landslide dam Figure 9 shows the simulation results for the progressive failure process of case 1. From 0 to 70 s, the groundwater level in the dam rose with the reservoir water level. At 70 s, the dam toe collapsed. From 70 to 384 s, the downstream dam slope collapsed intermittently, resulting in collapse progression from the dam toe. At 384 s, the reservoir water started to overflow the dam. This result indicates that the progressive failure model can express a series of progressive failure processes.
Temporal changes in the water levels Figure 10 compares the simulation and experimental results regarding the temporal changes in the water levels until the reservoir water overflows. The results of x = −1.78 m for case 1 and x = −1.30 m for case 2 indicate temporal changes in the reservoir water levels, and the other results indicate the temporal change in the groundwater level in the dam. In case 1, the simulation results after approximately 100 s are in agreement with the experimental results, except for x = −0.46 m (Fig. 10a). In case 2, the simulation results after approximately 150 s are in agreement with the experimental results, except for x = 2.83 m (Fig. 10b).
In the results of case 1 before 100 s and case 2 before 150 s, the groundwater level of the simulation results increased more gradually than that of the experimental results. This can be attributed to the differences between the actual and model dam cross sections ( Supplementary Fig. 3a and 3b). The actual dam cross section is a trapezoid whose bottom width is smaller than the top width ( Supplementary Fig. 3a), whereas the model dam cross section is rectangular (Supplementary Fig. 3b). These differences will change the rate of increase in the groundwater level.
The simulation results of case 1 for x = −0.46 m and case 2 for x = 2.83 m differ significantly from the experimental results. The former should have approached the reservoir water level (x = −1.78 m); however, the experimental result did not. These discrepancies may arise from measurement errors.
Collapse progression from dam toe Figure 11 compares the simulation and experimental results regarding the collapse progression from the dam toe. The vertical axis represents the horizontal distance between the toe of the initial dam shape and the top of the slip surface (L c ), as shown in Fig. 11. Although the experimental value of L c changed in the cross-flow (y-axis) direction, its maximum values are shown in Fig. 11. The simulation results for cases 1 and 2 are in agreement with the experimental results, except for case 1 before approximately 100 s. In the case 1 simulation before approximately 100 s, although the dam toe collapsed, the longitudinal dam shape was nearly the same as that of the experimental results. Furthermore, the simulation results and experimental results of the longitudinal dam shapes for case 1 are in agreement with the experimental results ( Supplementary Fig. 5).
The above-mentioned results show that the progressive failure model successfully reproduced the experimental results of the one-dimensional collapse progression from the dam toe. The experimental results of the progressive failure processes (Fig. 3) demonstrate that the collapse at the dam toe expanded not only in the flow direction but also in the cross-flow direction. Reproducing such a collapse expansion will require extension of (i) the two-dimensional intermittent collapse model to three dimensions, (ii) the two-dimensional seepage flow model to three dimensions, and (iii) the one-dimensional surface flow model to two dimensions.
Outflow hydrographs after reservoir water overflows Figure 12 compares the simulation and experimental results of the outflow hydrograph after the reservoir water overflows. The experimental results were given by the outflow discharge that included only water. The simulation results were thus given by the outflow discharge that includes only the water at x = 0.00 m. Figure 12 shows the simulation results for the two failure modes of progressive failure and overtopping failure. The simulation for progressive failure used both the progressive failure model and the overtopping erosion model. The simulation for the overtopping failure used only the overtopping erosion model, where the initial condition was given by the initial topography shown in Fig. 8.
The simulation results for the progressive failure were in agreement with the experimental results. The simulation results for the overtopping failure showed a larger peak discharge and later peak occurrence time than the experimental results. These results suggest that the progressive failure model is necessary for predicting a flood hydrograph due to the progressive failure of a landslide dam.  1. The landslide dam failure process observed in the present experiments was divided into two processes: progressive failure and overtopping erosion (Fig. 4). The progressive failure process is defined as the process in which the collapse at a dam toe proceeds upward until the reservoir water overflows. It includes the physical processes of the intermittent collapse of a dam slope, seepage flow in a dam, and surface flow on a dam slope (Fig. 4a). The overtopping erosion process is defined as the process in which an overtopping flow erodes a dam after the reservoir water overflows. It depends on the surface flow on the dam slope (Fig. 4b). 2. The proposed numerical simulation method involves the progressive failure model and overtopping erosion model (Fig. 5). The progressive failure model reproduces the progressive failure process to predict the longitudinal dam shape and reservoir water level when the reservoir water overflows. The overtopping erosion model uses these predicted values as the new initial conditions and reproduces the overtopping erosion process, predicting a flood hydrograph after the reservoir water overflows. 3. The progressive failure model includes physical models representing the intermittent collapse of a dam slope, seepage flow in a dam, and surface flow on a dam slope. In the intermittent collapse, the two-dimensional simplified Janbu's method predicts the collapse occurrences, and the equilibrium sediment transport equation represents the transportation of the collapse deposits downstream. The seepage flow is described by the two-dimensional Richards' equation. The surface flow is represented by the one-dimensional conventional simulation model of debris flow. The intermittent collapse Fig. 9 Numerical simulation results of the progressive failure processes of case 1 Fig. 10 Comparison between the simulation and experimental results regarding the temporal changes in the water levels until the reservoir water starts to overflow. a Case 1. b Case 2. The horizontal axis indicates the time elapsed since the reservoir water level started to rise. The simulated groundwater levels in the dam indicate the height at which the pressure head (ψ) has a value of 0.0 m model characterizes the progressive failure model. This model considers a stabilization effect in which the collapse deposits support a steep slope. This effect decreases as the collapse deposits are transported downstream. Such a consideration allows the model to express intermittent, not continuous, occurrences of collapses. 4. The progressive failure model successfully reproduced the experimental results of the one-dimensional collapse progression from the dam toe ( Fig. 11 and Supplementary Fig. 5). Using the value predicted by the progressive failure model, the overtopping erosion model successfully reproduced the flood hydrograph after the reservoir water started to overflow (Fig. 12). 5. The numerical simulation models proposed in the present study do not consider two-or three-dimensional phenomena. In the future, we will focus on the following three issues.
(1) The experimental results demonstrated that the collapse at the dam toe extended not only in the flow direction but also in the cross-flow direction (Fig. 3); however, the progressive failure model ignores the extension in the cross-flow direction. (2) The groundwater level of the simulation results increased more gradually than that of the experimental results (Fig. 10), which may be due to the difference between the actual and model dam cross sections (Supplementary Fig. 3a and 3b).  Availability of data and material The datasets during and/or analyzed during the current study are available from the corresponding author on reasonable request.
Code availability The codes during and/or analyzed during the current study are available from the corresponding author on reasonable request.