Development of an open source Debris Flow Simulator for “Sabo” (DFSS)

. We report the development of an open source “Debris Flow Simulator for Sabo” (DFSS). DFSS consists of three models corresponding to three geomorphological regions (mountain hillside slopes, mountain-to-plain river channels, and lower plains). A 2D rainfall-runoff (RR) model is adopted mainly for hillside slopes; a 1D debris flow runoff (DR) model is used for the river channel; and a 2D debris flow flooding or inundation (DF) model is adopted for the lower plains. These models are weakly coupled such that one can modify the overall model to fit diverse problems. In this article, we provide a detailed introduction to the models comprising DFSS, and demonstrate its utility.


Introduction
The purpose of this article is to introduce our recently developed open-source "Debris Flow Simulator for Sabo" (hereafter, DFSS; [1][2]). The Japanese word "Sabo" is interpreted to mean not only Sabo works such as check dams, but also an entire Sabo Project, including structural and non-structural measures against sediment and debris flow disasters [3]. Among the important nonstructural measures adopted is setting a debris flow evacuation zone. A range of physical models are used to predict areas at the risk of hazard and runout area associated with debris flows. It is crucial to select appropriate governing equations, geometric models, and spatial dimensions to adequately approximate specific problems. Furthermore, the recent expanding use of satellite-based spatiotemporal high-resolution data has enable the occurrence and magnitude of debris flows to be determined more accurately. Given this background, the development of DFSS was conceptualized to use a variety of distributed data, such as rainfall and land use, and to adopt a different spatial dimensional models to achieve a reasonable calculation time. DFSS reflects the heavy rainfall conditions and topography typical of Japan, which has many mountains and few plains. In this article, we provide a detailed introduction to DFSS and concisely present results using this approach, since [1] is currently only available in Japanese. It is hoped that this will encourage other researchers to use the open source code in [2]. Note that DFSS is currently in the early development stage.

Models and Data
DFSS consists of three models corresponding respectively to three regions (mountain hillside slopes, mountain-to-plain river channels, and lower plains) as shown in Fig. 1. A 2D rainfall-runoff (RR) model is adopted for hillside slopes; a 1D debris flow runoff (DR) model is used for the river channel; and a 2D debris flow flooding or inundation (DF) model is used for the lower plains. Standard calculations are performed independently in this order, and information regarding the flow is weakly coupled, such that the time-series output files of the previous model are used as input data for the next model. These output data to the next model could be regarded as the data such that real-world data which are observed in discrete time intervals. Although less accurate in its calculations than a strongly coupled model, where data connects models at each calculation time step (dt), the approach presented herein has the advantage of a configuration and weak coupling mode that are able to be controlled to fit a diverse range of problems. Furthermore, the calculation times associated with DFSS are particularly short when a 1D model is used. For example, if the target area features only a steep mountainous topography, the 1D DR model would provide a more appropriate model than the 2D DF model in terms of geometry and calculation time. When developing DFSS, we adopted the concept of easy modification of the model to cover a wide range of problems. The governing equations and constitutive equations of each model are based on previous research [4][5][6][7][8].

Configuration of topographic models
In the 2D region (mountain-side slopes and lower plains), the region of interest is divided into square grids according to the north-south and east-west directions as shown in Fig. 1. When setting up a non-calculation grid, a binary mask of any shape may be used. The 1D region (river channel) is represented by lines connecting the center of the grid (node) with the center of the adjacent grids in eight directions with the steepest downward gradients. In other words, the river channel is portrayed as a graph consisting of nodes and lines. The nodes store information including time and discharge. The upstream end of a river channel is determined by a grid with a catchment area that exceeds a user-defined threshold. The width and depth of the river channel are determined according to regime theory, and do not depend on the size of the square grid.

2D rainfall-runoff (RR) model
In the RR model, the amount of water supplied by rainfall is divided into surface water and subsurface (ground) water. As shown in Fig. 2, the continuity equation of water in topsoil layer holds as follows: where is the volumetric water content of the topsoil, is time, is the depth of topsoil layer, is rainfall intensity, and is evapotranspiration rate.
In the case of seepage flow, because the topsoil layer is a solid-water-air phase mixture, the porosity of topsoil layer can be integrated into the continuity equation of seepage flow as follows: where , are the axis parallel to slope and the vertical axis, ℎ is water depth in the topsoil, , are the and components of seepage flux , respectively, is the gradient of the slope, and 1 , 2 are the equilibrium infiltration rate of the surface layer and the bottom part of the topsoil layer, respectively. Assuming that the seepage flow is steady and uniform, the equation of motion is the same as Darcy's law as follows : where is the steepest gradient of the water surface and is saturated hydraulic conductivity. In a surface flow, the following continuity equation holds : where , are the and components of discharge flux , respectively, is the unit width discharge from the slopes to the river channel, and is river channel length : which is given the same value as grid length. Assuming a steady and uniform flow, the equation of motion using Manning's law holds as follows : where is the equivalent roughness coefficient of the slope. Discharge from the slopes to the river channel is calculated using overflow formula as follows : where is constant value of 0.35 which value is based on the design manual of river works of Japan ("Technical Criteria for River Works") and is acceleration due to gravity. For weak coupling, the time series data of the discharge from the slope to the river channel are stored in the nodesl, which output text files passing the information to the next model. When the RR model is used alone, runoff analysis in the river channel without sediment is performed ; this concept is not explained in detail herein.

Debris flow runoff (DR) model
The DR model inputs only the output information reflecing the RR flowing from the slope into the river channel to the nodes at the appropriate time during the DR calculation. In the DR model, the amount of sediement yield is calculated by analysis of slope stability in the river channel. The driving force is a body force from gravity and the resistance force represents Coulomb friction and adhesive force as follows: where is soil particle density, is water density, * is volume concentration of the topsoil layer, is the internal friction angle, and is the adhesive force at the bottom of the topsoil layer. If under state of limit equilibrium, = and solving for yields the limit equilibrium angle as follows: where is the critical gradient for the occurrence of collapse, and ′ is the dimensionless adhesive force : = ℎ�2 ℎ (6) https://doi.org/10.1051/e3sconf/202341502020 , 02020 (2023) E3S Web of Conferences 415 DFHM8 regarded as river bed sediment which is contributed to the debris flow through the source term.
In the DR model, the debris flow is a mixed flow comprising coarse and fine sediment, as well as water. Fine particles are treated as "particles dissolved in interstitial fluid" and, as a result the bulk density of the interstitial fluid is greater than that of water. Continuity equations of mixed flow as a whole, the flow of coarse sediment, and the flow of fine sediment hold as follows : where ℎ is the depth of the mixed flow, is the river width, is the depth averaged velocity of flow, is erosion/deposition velocity (described in detail below), * is the volume concentration of the static riverbed, is the discharge from the slopes to the river channel as RR output, is the overflow discharge from the river channel to the plain (which is passed to the DF model), and f are the depth-averaged volume concentrations of coarse particles in whole space and fine particles in interstitial space, respectively, is the sediment flux correction factor, and are the volume fractions of coarse and fine particles on the riverbed, respectively, and * is the volume concentration of coarse particles on the static riverbed after deposition. [5] is given by : where is the equilibrium gradient of the debris flow expressed as Takahashi's equilibrium condition for a debris flow of non-adhesive material [6]: where is the bulk (quasi-homogenous) density of an interstitial fluid with fine particle and water as follows: To describe the change of a river bed due to , the continuity equation of riverbed holds as follows: The equation of motion of the DR model is a depthaveraged one-fluid model as follows: where is the momentum correction factor, is the water level of the debris flow expressed as = + ℎ cos , and is basal stress acting on the riverbed, which is expressed as follows [7]: where is yield stress and is the flow resistance coefficient. These are expressed as follows: where is the representative diameter of a coarse particle, and are constant parameters with values of 0.16 and 0.0828, respectively [5], is the restitution coefficient between coarse particles, and b is as follows [8]: where is constant (8.5), is the Karman constant (0.4), is the diameter of coarse sediment, . These constants are values usually used in the log law of velocity profiles on a movable bed. Note that the first term of right-hand side of eq. (17) could be expressed as follows: where / = sin . In this expression, the term of static pressure and gravity acceleration is clarified.

Debris flow flooding (DF) model
The DF model calculates the inundation overflowing from river channels to the plain. Places prone to flooding along the river channel are determined by DR calculations. The governing equations of the DF model are similar to those of the DR model described above, with the exception of dimensions (1D and 2D). From here onward, only continuity equations and equations of motion in 2D form are described. Continuity equations of the mixed flow as a whole, a flow of coarse and a flow of fine sediment are, respectively, as follows: where and are the depth averaged velocity of mixed flow along the x and y directions, respectively, , and are the total, coarse, and fine sediment fluxes from the river channel to the plain, respectively, which are given as outputs of the DR model calculation.
Equations of motion of x and y direction are, respectively, as follows: where and are basal stresses (ref. eq. 18) along the x and y directions, respectively. Note that the relationship of eq. (24) can be applied to first term of the right-hand side of these equations.

Data
On July 2 nd , 1992, a large scale debris flow disaster occurred at Aso-cho, Kyusyu, Japan (Fig. 3). Aso is famous for its large caldera and the surrounding hillsides feature large amount of unstable sediment.
Data such as soil parameters and deposition . area from [9] were used for confirmation of the applicability of DSFF to debris flow disasters in volcanic area. Rainfall data from [9] was also used by the Automated Meteorological Data Acquisition System of Japanese Meteorological Agency (see [1] for further detail).

Results
The results of our calculations are shown in Fig.  4. The calculated area was got by the DF model using the output of the RR and DR model as input files. As shown in Fig. 4, comparison between the calculated area of maximum depth and the area of deposition determined from aerial photography [9], shown in the dotted area, indicates that the calculated results match the observed depositional area near the mountain, but do not match the broader observed depositional area, especially in the lower region. The reason for this inaccuracy is thought to be the topography of the river channel in the mountain, which is confined by the surrounding hillsides. As such, the results of DFSS, especially in 2D calculation, were found to be sensitive to digital elevation models.
In order to accurately estimate the inundation area across the plains, it will be necessary to improve the DFSS by considering the presence of various detailed features such as topography and Sabo structures, these as elevation data or boundary conditions.