Predictive Simulation of Concurrent Debris Flow: How Slope Failure Locations Affect Predicted Damage

Predictive simulation of concurrent debris flow using only pre-disaster information has proven to be difficult as a result of problems in predicting the location of debris-flow initiation (i.e., slope failure). However, because catchment topography has concave characteristics, with all channels in a catchment joining each other as they flow downstream, it is possible to predict damage to downstream area using relatively inaccurate initiation points. Based on this, this paper presents methodologies employing debris-flow initiation points generated randomly using statistical slope failure prediction. A many-case simulation across numerous initiation points was performed to quantify the effect of slope-failure location in terms of deviations in the predicted water level and terrain deformation. It was found that the relative standard deviation diminished as the points approached the downstream area, indicating a location-based predictability effect. (131 words)


Introduction
Debris flow is a primary hazard to human life. When significant rainfall events impact watershed topography, multiple debris flows take place concurrently within a single catchment. A case occurred on Northern Kyushu Island, Japan, in July 2017, when heavy rains induced multiple landslides, producing debris flows in the Akatani River catchment with areas of approximately 20km2. Subsequently, massive amounts of sediment covered a residential area located at the bottom of the valley. Sediment was also transported several kilometers downstream and deposited at widths of approximately 150m through a channel with an original width of only about 10m. The sediment-laden flood damaged both channel sides and areas far from the main channels. The disaster killed 33 people, and it damaged 1,046 buildings throughout Asakura city according to the report by Fukuoka pref., Japan.
The event was remarkable in that it produced sediment from multiple landslides, which strongly affected the inundation processes. Large terrain deformation around channels caused by the flood gave severe damage to the areas outside of the hazard map.
Numerical simulation is an important method for assessing risk and designing counter-measures against debris-flow-and sediment-laden floods. Such simulation models (e.g., Naef et al. (2006); O'Brien et al. (1993); Hungr and McDougall (2009)) normally depend on theory derived from flume experiments or field investigations. By employing these techniques, reproductive debris-flow simulation is possible using post-disaster information, including field survey data such as peak discharge, area of deposition, material size, and depth of erodible soil. However, the limited availability of prior field data means that the precise risk of disaster cannot be obtained via simulation during the pre-disaster phase and, to date, no effective predictive prior-disaster simulation schemes have been developed.
Given this limitation, model simulations are typically carried out using a hydrograph that flows from the outlet of the debris-flow-prone catchment (e.g., Rickenmann et al. (2006); Nakatani et al. (2016)). In general, such hydrographs are modeled using rainfall-runoff simulation or rainfall data processing. However, the prediction of the potential scales of heavy rainfall events requires significant effort. Furthermore, because debris flow can grow by eroding and consuming bed material, the debris-flow peak discharge tends to be higher than that of floods, which can be predicted using a simple rainfall runoff simulation (Hungr et al., 2001). Thus, hydrograph-based condition-setting for debris flow simulation usually requires empirically based modification.
Focusing on landslide-induced debris flow alone enables the use of an alternative approach in which the debris flow is simulated from its initiation based on its development process (Revellino et al., 2004). This requires no assumptions regarding the rainfall-runoff process, although some assumptions regarding how the debris flow initiates at the landslide are needed. Furthermore, precise landslide data are needed to carry out the simulation. Unfortunately, precise deterministic landslide prediction is difficult using current techniques because of inequalities in underground conditions (cohesion, inside friction angle, water retention, and hydraulic conductivity) and observational difficulties in obtaining such data.
However, because watershed topography has a concave characteristic, in which tributaries aggregate as they flow downward, similar scales of debris flow can result even if the landslide initiation point changes within the catchment. This suggests that a relatively low prediction accuracy in terms of debris flow initiation location might not be unacceptable in determining the risk to downstream areas. Using this hypothesis, we could quantify the uncertainty in debris-flow initiation location by conducting a multi-case simulation in which the initiation point was varied using a statistical method. This paper presents following findings: A parallelized debris-flow simulation code that can be applied using a supercomputer to conduct simultaneous multi-case simulations was developed. The code was validated by reproductively simulating the 2019 Northern-Kyushu heavy rainfall event. To obtain the slope failure likelihood, we applied a simultaneous multi-case simulation using 60 patterns of artificial landslide initiation point data produced using logistic regression. The 60 results were assessed by calculating the average value, standard deviation, and relative standard deviation to quantify the effects of slope failure location on predicted damage.

Governing equations
To treat debris flow, we selected a two-dimensional simulation model based on the dilatant fluid model developed by Takahashi (2007), which has been frequently applied in practice (Wu et al., 2013;Nakatani et al., 2016). The governing equations can be written as , ( where is the conservative variable vector, and are the flux vectors in the xand y-directions, respectively, and is the source vector. This equation is the integrated form of shallow-water equations and equations for sediment transport comprising the volumetric conservation law of the fluid (debris flow), momentum equations for the x-and y-directions, the volumetric conservation law of concentration of sediment in fluid, and the land-surface material conservation law.
The vectors are defined as where ℎ is the water depth, and are the depth-averaged velocities in the anddirections, respectively, is the concentration of sediment in the fluid, is the landsurface elevation, is the acceleration of gravity, and * is the sediment concentration in the soil at the ground surface. The topographical gradients in the x-and y-directions, and , respectively, are calculated as / , / . is the eddy momentum diffusivity, which is calculated as where is von Karmans' constant, * is the frictional velocity, and and are the frictional gradients between the fluid and bed surface in the x-and y-directions, respectively. To consider sediment-laden flooding in addition to debris-flow, the threephase (stony debris, hyper-concentrated, and water flow) model proposed by Takahashi and Nakagawa (1991) is used: where is the representative sediment diameter, is Manning's roughness coefficient, and is the erosion/deposition velocity in the vertical direction, which is calculated as where and are the coefficients for erosion and deposition, respectively, and and are the specific weights of water and sediment, respectively. The positively valued expresses the erosion speed [m 2 /s], where is the equilibrium sediment concentration obtained by 0.9 * , t a n tan tan tan tan , t a n tan 0.138 6.7 tan tan tan , 0.138 tan 0.03 where is the water surface gradient, is the internal friction angle, and * is the critical non-dimensional tractive force calculated by * 0.04 10 . , where * is the non-dimensional tractive force obtained by * ℎtan , where is calculated by 2 0.425 tan 1 tan (9)

Numerical modeling
To prevent the occurrence of numerical instability at the boundary of the supercritical and subcritical flows, the MacCormack scheme with artificial viscosity is used to implement the discretization equations. In this method, calculation is carried out in two steps-prediction and collection. If the backward and forward differences are used to carry out the predictor and collector steps, respectively, the MacCormack scheme can be expressed as follows. Predictor: Collector: where, where is the coefficient vector for artificial viscosity.

Code Parallelization
A numerical simulation code for carrying out the calculations above was implemented in Fortran. As changing conditions to apply the program at the watershed scale could lead to considerable computational cost, we parallelized the code to enable the use of high-performance computing techniques (supercomputing).
To do so, we introduced thread-parallelization to the program using OpenMP, a technique for parallelizing program loops. The use of OpenMP enabled the parallelization of most of the loops that did not cause problems related to dependencies. To enable execution of the program on supercomputers comprising multiple nodes, we implemented MPI parallelization, in which a target area is decomposed into multiple areas overlapping two arrays of numerical grids. Each area is allocated to corresponding nodes, and the simulation is executed independently on each node. Using MPI, the variables on overlapping grids were synced just after the predictor and collector steps to retain collective variables on the boundary area. As a result of this hybrid parallelization process, a one-case simulation of a 8km × 9km domain could successfully be conducted on 2,304 CPU cores (288 nodes, eight cores per node); for the full study, a 60-case parallel simulation was conducted using 138,240 cores (17,280 nodes, eight cores per node).

Treatment of debris flow initiation
The simulation framework used in the study was based on fluid dynamics; to assess debris flow using the framework, the flow had to be initiated as a state of the fluid. In general, there are three methods for debris flow initiation: 1) fluidization of a shallow landslide; 2) collapse of a natural dam composed of landslide mass; and 3) erosion of river bed material by overland flow (Takahashi, 2007). As the initiation process includes many unknowns resulting from a lack of in-situ observations, a number of assumptions are made in conducting debris-flow simulations. The first is to assume the hydrograph and sediment concentration values at the inlet of the target domain Chen et al., 2017;Han et al., 2018;Bao et al., 2019;Nakatani et al., 2016;Frank et al., 2015;Gao et al., 2016). In this method, the inlet point is usually set as the initiation point or as the location at which the debris flow is sufficiently developed. Using this method, the simulation can be simplified even if all elementary processes (e.g., rainfall infiltration, landslides, and development of debris flow) are neglected. However, the hydrograph setting requires an insitu trace of debris flow, which makes predictive simulation difficult.
Another approach involves connecting a rainfall-infiltration-runoff simulation to the surface erosion model (van Asch et al., 2014;Melo et al., 2018). Although this can precisely treat debris flow initiation types 1) and 2), it requires multiple parameters related to water transportation (e.g., coefficient of saturated/unsaturated permeability) and soil strength (e.g., cohesion coefficient, internal friction angle). In general, these widely-distributed physical underground parameters cannot be uniquely determined at the watershed scale, and there are no well-established remote observation techniques for obtaining them. Thus, this method is not suitable for predictive simulations. A third approach is to set the free-flowing landslide mass at the location of slope failure (Revellino et al., 2004;Schraml et al., 2015;Rodríguez-Morata et al., 2019). Although this can treat only initiation type 1), it requires fewer parameters than the second approach. As the effect of landslide volume on debris-flow discharge is limited, the distance between the slope failure and the damaged area will be sufficiently large if there is erodible material in the slope surface; therefore, accurate simulation should be obtainable using this methodology if sufficient slope failure location data are available. We adopted this method for the modeling of landslide initiation in our study. In the model, the depth of the freeflowing landslide mass ℎ was assumed to be equal to the landslide depth , while the sediment concentration of the mass was assumed to be equal to a specific value, .

Calculation condition
To validate the proposed simulation scheme, we applied it to an actual disaster event that happened on Northern Kyushu Island, Japan in 2017. An area spanning 9 and 8km in the east-west and north-south directions, respectively, around the Akatani River catchment area was selected as the target region. We used the following method to convert the polygonal Geospatial Information Authority of Japan debris-flow traces into a point data set of debrisflow initiation. We defined the debris-flow initiation point in each debris flow trace as the highest point. As most of the debris-flow traces joined each other, we assumed that the highest-elevation cells were the high points within a specific diameter and for the debris-flow trace. The resulting initiation points and the elevation data for the simulation are shown in Fig. 1 (b), which also shows the trace of the flood. The simulation was conducted under the assumptions that all areas were completely saturated, that the erodible soil depth was equally distributed at 1 m, and that all debris flows initiate simultaneously upon commencement of the simulation (T=0 [s]). Other parameters used in the simulation were 1.0, * 0.6, 0.5, 2.65, 5 cm , 35 deg. , 0.05,, and 0.0007.

Calculation results
The maximum water level and the terrain deformation (i.e., erosion/deposition) are the critical variables for representing damage caused by the inundation of sediment-laden flood. In this study, these factors were determined by the relative elevations of the water level and ground surface from the initial ground surface, respectively. The maximum water level during the simulation and the final terrain deformation are shown in Fig. 2 (a). A rough comparison reveals that the distribution of each variable agrees with the trace of the disaster shown in Fig. 1 (a).
For the quantitative comparison, we selected three major internal catchments to neglect any areas outside of the target area that might be affected by the phenomenon. The available trace data (shown in Fig. 1 (a)) were categorized into flooded and landslide areas. The trace was obtained by manual categorization based on interpreting an aerial photograph, suggesting that the data included areas with decreases and increases in elevation (erosion and deposition, respectively) and water-covered (inundated) areas. To account for this, three binary maps (flooded, landslide, and both) were compared to obtain the "true" map. To obtain a predictive map, three thresholdings using the maximum water level and terrain surface deformation were used to produce flooded, landslide, and damaged areas. The landslide area was first compared to the area in which the simulated absolute value of terrain deformation was greater than the threshold value of 0.01 m, as shown in Fig. 2 (c). The comparison produced a high recall value (= 0.744), indicating that the simulated and actual landslide areas mostly coincided. The flooded area was then compared with the area in which the simulated maximum water level was greater than the threshold value (0.01 m), as shown in Fig. 2 (d). This result also produced a high recall value (= 0.747), indicating that the simulated high water-level and inundated areas mostly coincided. Finally, the damaged area itself was determined by comparing the combined flood and landslide area with the area in which the absolute value of was greater than 0 , as this should correspond to locations in which water was present at least once (Fig. 2 (e)). The results had high recall, accuracy, and precision values (0.794, 0.881, and 0.553, respectively), indicating a high degree of coincidence between actual results and prediction. These comparison results allowed us to conclude that the proposed simulation method has a reproductivity suitable for disaster prediction work.

PREDICTIVE SIMULATION EMPLOYING STATISTICAL LANDSLIDE PREDICTION
Deterministic prediction of slope failure location is difficult using current techniques because underground physical parameter conditions are unobservable and unpredictable and there are no robust techniques for the simulation of complex actual topographies. However, watershed topography in mountainous regions has a "concave" characteristic in which all streams converge downstream. In other words, unless the debris flow occurs within different channels in a single Strahler stream order, a given scale of terrain deformation or water level might be reproducible at two or larger orders. If so, the error in predicted location might not strongly affect the downstream area damage. This suggests the possibility of a deterministic, statistically based approach to slope failure prediction. In this section, we describe a statistical landslide prediction model based on logistic regression and use the model to generate artificial landslide datasets. The generated landslide data were applied as input data to a predictive simulation to quantify landslide location uncertainty.

Statistical landslide prediction model for generation of artificial landslide dataset
Debris flow initiation points, i.e., slope failures or landslides, are essential inputs to the proposed simulation approach. Note that a fully predictive simulation requires data sourced solely from pre-disaster information. Approaches to generating landslide susceptibility maps include geomorphological mapping, analysis of landslide inventories, heuristic or index-based approaches, process-based (physically-based) methods, and statistically based modeling (Reichenbach et al., 2018). Of these, physically and statistically based modeling are preferred for quantitative evaluation (Reichenbach et al., 2018). Among the physically based models, which have advantages in terms of physical validity, the infinite slope stability model is a simplified, widely used approach Raia et al., 2014). Although obtaining reliable results requires accurately modeling complicated, heterogeneous underground structures based on field observation or parameter fitting using landslide records, it is not always necessary to use difficult-to-obtain physical parameters in statistically based approaches, and as long as reliable landslide records are used to obtain the statistical relationships, no abnormal results should be produced.
We therefore used a logistic regression model-a preferred approach in statistical modeling-to derive the initiation point likelihood at each top point in the debris-flow trace. The input variables for this model are shown in Figs. 3 (a)-(e). To improve prediction quality, several variables are generally selected as explanatory variables, including curvature, slope, aspect, accumulation, elevation, lithology, precipitation, etc. However, the relationship between landslide occurrence and aspect or elevation is location-dependent and, to ensure generality, we used only four variables: slope, accumulation, and profile/tangential curvature. The predicted probabilities at a 10-m resolution are shown in Fig. 3 (f).
To enable simulation inputs, the probability map had to be converted into landslide maps. To do this, we generated pseudorandom numbers between zero and one for all meshes in the target areas. For results smaller than the corresponding cell probabilities, the mesh was selected as a debris-flow initiation point. By changing the random-seed generator, 60 sets of debris-flow initiation points were created; these are partially shown in Figs. 4 (a) and (d).  (j), (k), and (l) show, respectively, average value, standard deviation, and relative standard deviation for terrain deformation.

Multi-case predictive simulation for predictive simulation using artificial landslide datasets
Sixty simulation cases were conducted using the artificial debris-flow initiation points as input data. Other than the initiation points, all conditions and parameters were held constant across the simulations. Some of the obtained maximum water levels and terrain deformations are shown in Figs. 4 (b), (c), (e), and (f). Roughly similar values of maximum water level and terrain deformation occur in the downstream area, suggesting that the downstream area damage is not sensitive to the landslide initiation point distribution. By contrast, there is high variability in the upstream area closer to the debris-flow initiation points. The variability among the 60 results was quantified using the averaged value (AVE) and standard deviation (SD); relative standard deviation (RSD; =SD/|AVE|)) values were calculated for maximum water level and terrain deformation (Figs. 4 (g)-(l)).
Figs. 4 (g), (h), and (i) show the AVE, SD, and RSD for the maximum water level, respectively. At the outlet of the debris-flow prone basins, the values of AVE range from 0.5-1.0 m; at the bottom of the valley in the downstream area, the value is 1.0-3.0 m, indicating a general increase as the streams flow downstream. By contrast, there are only modest variations in SD of between 0.0-0.2 m. Thus, the RSD values are larger in the upstream space and smaller downstream. As the ratio of the standard deviation and averaged value to the predicted value of the water level, these RSD values indicate the prediction uncertainty owing to landslide location and the predictability of the maximum water level when parameters/conditions other than landslide location are appropriately provided. These results suggest a water-level "predictability" areas located in the valley bottom or downstream alluvial fan that decreases upstream. By contrast, the areas around the outlets of small debris-flow-prone streams have RSDs of approximately 0.5-1.0, indicating a higher degree of predicted-damage uncertainty. Overall, the required precision of landslide location for debris-flow damage prediction differs by location.
The SD, AVE, and RSD values for terrain deformation are shown in Figs. 4 (j), (k), and (l), respectively. Each index follows a trend similar to the corresponding maximum water level index, although the RSD value in the residential area at the bottom of the valley is slightly larger than that for maximum water level.

Conclusions
This paper discussed the slope failure condition settings needed to predict debrisflow and sediment-laden flooding. As accurate slope failure prediction is difficult using current techniques, we proposed a predictive simulation method employing conventional statistics-based landslide prediction. We generated artificial initiation points using the probability of slope failure (i.e., the probability that a given point is a debris-flow initiation point) and pseudorandom numbers to change the slope failure location based on the random seed. In concave catchment topographies, individual debris flows merge as they flow downstream. Under such conditions, different debris-flow initiation points result in similar damage locations. To quantify the similarities, we calculated the RSD among 60 cases with different random seeds. The RSD diminished in downstream regions with large drainage areas.
Our results suggest the following: The damage produced by debris-flow-or sedimentladen floods in the downstream regions of watershed-like topographies is predictable under certain conditions regardless of the location of the sediment sources (slope failures). The dominant factors in terms of predictability are the presence of a catchment-like concave topography, sufficient water to develop a debris flow, and sufficient density of slope failure to join multiple debris flows. Thus, for concave topographies, it is possible to predictively simulate damage resulting when high rainfall yields multiple slope failures and debris flows.
The proposed method relies on several assumptions: • All debris flows initiate simultaneously.
• All soils are completely saturated.
• Erodible soil depth is limited to 1m in all regions.
• The debris flow material has homogeneous grains.
• Only water in the soil layer, and not additional rainfall, is considered. Further studies will be required to estimate how each assumption affects the simulation results. However, this study demonstrated the effectiveness of using HPC-based multi-case simulation to evaluate the effects of debris-flow initiation location. A similar approach involving changing each parameter or condition will be used in a future study to quantitatively evaluate each assumption.

Acknowledgments
This work was supported by JSPS KAKENHI Grant Number 19K15105 and FOCUS Establishing Supercomputing Center of Excellence (COE). The data that support the findings of this study are available from Geospatial Information Authority of Japan. Restrictions apply to the availability of these data, which were used under license for this study. Data are available https://www.gsi.go.jp/kiban/ and https://www.gsi.go.jp/BOUSAI/H29hukuoka_ooita-heavyrain.html with the permission of Geospatial Information Authority of Japan.