Optimizing the Numerical Simulation of Debris Flows: A New Exploration of the Hexagonal Cellular Automaton Method

: Debris flow, driven by natural events like heavy rainfall and snowmelt, involves sediment, rocks, and water, posing destructive threats to life and infrastructure. The accurate prediction of its activity range is crucial for prevention and mitigation efforts. Cellular automata circumvent is the cumbersome process of solving partial differential equations, thereby efficiently simulating complex dynamic systems. Given the anisotropic characteristics of square cells in the simulation of dynamic systems, this paper proposes a novel approach, utilizing a hexagonal cellular automaton for the numerical simulation of debris flows, where the direction judgment efficiency increased by 25%. Employing cubic interpolation, the model thereby determines the central elevation of each hexagonal cell. By modifying the flow direction function and stopping conditions, it achieves more accurate predictions of the debris flow run-out extent. This method was applied to the 2010 Yohutagawa debris flow event and the flume test. To evaluate the simulation’s accuracy, the Ω value and F β score were used. The Ω value is a comprehensive evaluation factor that takes into account missed or misjudgment areas. On this basis, the F β score emphasizes that the missed identification of debris flow areas will bring greater harm. Research indicates that the Ω value showed improvements of 6.47% and 3.96%, respectively, while the F β score improved by 3.10% and 4.61%.


Introduction
The onset of a debris flow requires large amounts of loose solid material and steep terrain [1,2].Rain or snowfall acts as a trigger, pushing water, sand, and stones forward [3].Mountainous areas offer ideal conditions for debris flow formation [4].Around 66% of China's landmass consists of mountainous landscapes, especially in the transition zone between the second and third topographic levels [5,6].Active geological processes significantly contribute to the accumulation of loose materials, facilitating the formation of debris flows [7].Moreover, China's predominantly monsoon climate, with heavy rainfall in the wet season, fosters favorable conditions for debris flow initiation.Debris flows cause severe damage to residences, agricultural land, and transportation infrastructure, leading to casualties and extensive property damage.The urgency of enhancing debris flow mitigation and reducing their impacts is increasingly prompting scholars to focus on predicting their impact range [8,9].Precisely predicting the path and scale of debris flows allows for the rational design and planning of transportation infrastructure like roads and railways, preventing the development of vital transport hubs in debris flow-prone areas [10][11][12].
In the field of debris flow disaster mitigation, numerical technologies play a crucial role in disaster prevention [13,14].Many studies focus on numerically solving the Navier-Stokes (N-S) equation.These studies encompass approaches such as the shallow water equations (SWEs) and the smoothed particle hydrodynamics (SPH) method.
Water 2024, 16, 1536 2 of 18 The shallow water equations are implemented within a grid-based numerical framework specifically designed to simulate hydrodynamic flows under the assumption of hydrostatic pressure distribution.This methodology simplifies the Navier-Stokes equations by averaging vertical field variables, thereby efficiently modeling the flow over complex terrain with reduced computational demand [15][16][17][18].On the other hand, the smoothed particle hydrodynamics method employs a particle-based Lagrangian framework, eschewing traditional grid constraints for a more flexible representation of fluid dynamics [19][20][21][22].Through the use of kernel interpolation and particle discretization, the smoothed particle hydrodynamics method possesses the capability to handle physical processes such as large deformations and free-debris flow characteristics such as velocity, depth, and spatial positioning.Nonetheless, they consistently require significant computational resources [23,24].Additionally, acquiring parameters related to rheological properties, such as viscosity, presents notable costs [25].The above issues make physical-based methods still problematic for field-scale applications [26,27].
The cellular automaton (CA) model represents a numerical framework that partitions the entire area into discrete units, known as cells, each characterized by specific variables [28][29][30].Since the 1980s, CA models have been widely used in the simulation of a variety of natural disasters, including debris flows, landslides, avalanches, and rockslides [31,32].These applications underscore the efficacy of the CA model in simulating the dynamic processes of natural disasters and highlight its potential in disaster risk assessment and management strategies [33][34][35].Within the CA framework, the transition function incorporated with the probabilistic iterative algorithm is the key to controlling the dynamic evolution of the debris flow.Gamma [36] initially designed this framework, which is based on the Monte Carlo iterations, to advance the integration of the flow path through a topography-based model.Subsequently, other refinement mechanisms, including friction models [37], bedload transfer [38,39], and stopping conditions [40,41], have been continuously integrated into the paradigm to better constrain the flow dynamics of the debris flow.Han [42] revised the persistence weight parameters based on the flume experiments and quantified the inertial effect based on the persistence functions.Ma [43] further introduced a roughness function to depict the impact of bed roughness on the flow direction.
Historically, research applications of the CA models have predominantly utilized square cells.However, this architecture encounters challenges related to inconsistent neighbor relationships and directional biases, ultimately leading to anisotropic issues [44].Therefore, introducing a lattice structure that aligns more consistently and compatibly with the motion of debris flows is imperative to accurately characterize the evolution of the flow direction.Sousa [45] found that hexagonal digital elevation model data provide a higher spatial resolution compared to a square-cell structure.To simulate the movement of debris flows, Frisch [46] discrete the Navier-Stokes equation in the context of a hexagonal lattice.Ambrosio [47] and Avolio [48] developed a hexagonal-based CA model named SCIDDICA and applied it to real-world cases.Additionally, the hexagonal CA has fewer neighbors than the square-cell structure, which includes only six neighbors and leads to a 25% efficiency improvement in efficiency when determining the flow direction [49].
This study integrates the persistence function, previously validated in square cellular automaton (SCA) models, into the hexagonal cellular automaton (HCA) scheme.This modification aims to refine the model's ability to accurately simulate the flow state of debris flows.The integration of the persistence function into this novel context is expected to enhance both the robustness and simulation accuracy.Furthermore, the stopping conditions, a cornerstone of this framework, primarily address the conditions under which debris flows halt based on the maximal reach of debris accumulation.By accounting for both the "sink-filling" [36] and "maximum path length" considerations, the proposed method controls the dynamic evolution of the debris flow.
In general, this paper presents the HCA model for the numerical simulation of debris flow, which involves dividing the terrain into a network resembling a hexagonal honey-comb.This approach focuses on investigating the dynamics of the debris flow direction and the stopping conditions applicable to this cellular structure.Finally, the run-out extent of the debris flow is determined through Monte Carlo iterations, integrating all flow paths to effectively address the anisotropic challenges that SCA encounters in simulating debris flow dynamics.

Model Introduction
A standard CA model comprises several core components as follows: cells, lattice, neighborhoods, and transition functions.

•
Cell and Lattice The cell is the most fundamental unit of a CA model, serving as the system's basic building block.The lattice describes how cells are spatially arranged, including the type, size, structure, and boundary conditions of the cellular network.
In past research, square cells were more popular due to their inherent advantages in data representation.Initially, it is crucial to acknowledge that original terrain data frequently take the form of digital elevation model (DEM) raster data.In the DEM, the map is partitioned into a matrix of pixels, each pixel delineating a square shape and encoding specific values to represent the features or attributes of the corresponding location on the map.The alignment of the raster data's square pixels with the square-cell data highlights the inherent compatibility between DEM formats and the structure framework of CA models, thereby facilitating straightforward data utilization without requiring extensive processing.Horton [37] and Ma [43]'s research shows that the run-out extent of the debris flow can be accurately simulated using a cell size of 2.5 m.We adopted hexagonal cells as the lattice format, which requires additional processing of the cells and lattices before discussing the transition functions.Firstly, there is the selection of cell orientation, which can be either flat-topped or pointy-topped, as seen in Figure 1a.Next, is the choice of the lattice offsetting method.For the flat-topped cells in Figure 1c, the options include odd-column-down offsetting or odd-column-up offsetting.(If using pointy-topped cells, one must choose between odd-row-left offsetting or odd-rowright offsetting).The above processing flow will affect the matrix-coordinate mapping relationship.It is important to note that the coordinate mapping relationship varies for different cell and lattice forms.Equations ( 1) and (2) depict the scenario of "flat-topped cells" and "odd-column-up offsetting".The coordinates (x, y) of the lattice with index (m, n) in a two-dimensional matrix are given by: When n is an odd number, While n is an even number, In the formula, S x and S y refer to the distance in the x and y directions between the centers of the hexagonal grid cells, respectively.The relationship between S x , S y and the cell size of the neighboring unit center (D Hex ) is shown in Figure 1b and Equations (3) and (4): After completing the coordinate mapping, obtaining the elevation of the cells became an urgent issue.Considering both interpolation accuracy and efficiency, we opted for cubic interpolation.This aspect will be further elaborated in Section 4.

•
Neighborhood: The concept of neighborhood refers to the range of cells in which the center of a cell is affected by surrounding cells, usually defined by its radius.Figure 2 illustrates the cell arrangement and neighborhood relationship of the SCA models and HCA models in (a) and (b), respectively.In Figure 2a, the distance between the adjacent edge grids is D Square , while the distance between the adjacent corner grids is √ 2 × D Square .
Water 2024, 16, x FOR PEER REVIEW 4 of 19 After completing the coordinate mapping, obtaining the elevation of the cells became an urgent issue.Considering both interpolation accuracy and efficiency, we opted for cubic interpolation.This aspect will be further elaborated in Section 4.

•
Neighborhood: The concept of neighborhood refers to the range of cells in which the center of a cell is affected by surrounding cells, usually defined by its radius.Figure 2 illustrates the cell arrangement and neighborhood relationship of the SCA models and HCA models in (a) and (b), respectively.In Figure 2a, the distance between the adjacent edge grids is   , while the distance between the adjacent corner grids is √2 ×  .
After completing the coordinate mapping, obtaining the elevation of the cells became an urgent issue.Considering both interpolation accuracy and efficiency, we opted for cubic interpolation.This aspect will be further elaborated in Section 4.

•
Neighborhood: The concept of neighborhood refers to the range of cells in which the center of a cell is affected by surrounding cells, usually defined by its radius.Figure 2 illustrates the cell arrangement and neighborhood relationship of the SCA models and HCA models in (a) and (b), respectively.In Figure 2a, the distance between the adjacent edge grids is   , while the distance between the adjacent corner grids is √2 ×  .• Transition Functions: The CA model cell states transition to another state through transition functions, which serve to employ rule-based simulation processes rather than traditional mathematical functions for simulation.The CA method derives the development of complex phenomena with simple transition functions and a massive iteration.Instead of solving the complex PDEs, this approach employs interaction relations among neighboring cells [36,42,43].These transition functions are not only simpler and more straightforward to express and compute than the physical solutions required in dynamic models, but they also facilitate Water 2024, 16, 1536 5 of 18 a more accessible understanding of the system's processes.Refer to Sections 2.2-2.4 for detailed information.

Flow Direction Function 2.2.1. Terrain Probability
In the dynamics of debris flow, there is a strong tendency for the flow to follow the path of maximum slope descent.In SCA models, the D8 algorithm [50] is a common method for determining the flow direction.It primarily identifies and utilizes the steepest slope direction of each unit flow path as the debris flow direction.However, the D8 algorithm, with its single flow direction approach, does not adequately handle the divergent nature of the debris flows observed in real-world scenarios, especially in areas with gentle slopes where divergence is pronounced, resulting in significant inaccuracies.
To overcome these limitations, our terrain function calculates slopes in all directions within a hexagonal honeycomb network, considering the divergence phenomenon inherent to debris flows.We propose that cells with a positive slope gradient have a probability of directing the flow.Since a debris flow is unlikely in directions with a negative longitudinal slope, these cells are assigned a flow direction probability of zero.
The probability of the flow direction, denoted as P t f i , is calculated using Equation (5).
As shown in Figure 3c, β i and β j represent the longitudinal slope descent in directions i and j, respectively.They are set to zero when the longitudinal slope in those directions is negative.At a positive slope, the Equation ( 6) calculation formula is as follows: ∆H i represents the vertical slope in the i direction relative to the central cell.

Persistence Probability
When debris flows on flat terrain, the movement is mainly forward, aligned with the velocity direction.However, as the average flow velocity decreases, the debris flow not only continues to advance but also begins to spread.To describe this behavior, we introduce the concept of "persistence probability," which captures this phenomenon.
In Section 2.2.1, our analysis of the flow direction function predominantly focuses on the static properties of the terrain.Moving forward, this section delves into the dynamic attributes of debris flow.By integrating these dynamic aspects, we aim to enhance the

Persistence Probability
When debris flows on flat terrain, the movement is mainly forward, aligned with the velocity direction.However, as the average flow velocity decreases, the debris flow not only continues to advance but also begins to spread.To describe this behavior, we introduce the concept of "persistence probability," which captures this phenomenon.
In Section 2.2.1, our analysis of the flow direction function predominantly focuses on the static properties of the terrain.Moving forward, this section delves into the dynamic attributes of debris flow.By integrating these dynamic aspects, we aim to enhance the precision of our debris flow modeling.In Figure 4a, the debris flow enters the central cell along a trajectory shown by a red arrow.The persistence probability suggests that in the next time step, the debris flow could transition into the cell pointed by a black arrow.We assign persistence weights to each neighboring cell, reflecting the dynamic characteristics and momentum of the debris flow.

Persistence Probability
When debris flows on flat terrain, the movement is mainly forward, aligned with the velocity direction.However, as the average flow velocity decreases, the debris flow not only continues to advance but also begins to spread.To describe this behavior, we introduce the concept of "persistence probability," which captures this phenomenon.
In Section 2.2.1, our analysis of the flow direction function predominantly focuses on the static properties of the terrain.Moving forward, this section delves into the dynamic attributes of debris flow.By integrating these dynamic aspects, we aim to enhance the precision of our debris flow modeling.In Figure 4a, the debris flow enters the central cell along a trajectory shown by a red arrow.The persistence probability suggests that in the next time step, the debris flow could transition into the cell pointed by a black arrow.We assign persistence weights to each neighboring cell, reflecting the dynamic characteristics and momentum of the debris flow.The concept of persistence weight, denoted as  , plays a crucial role in the debris flow simulations by considering the position of the inflow cell.For a given incoming cell, we derive the forward and oblique direction numbers from Table 1 and then obtain their persistence probability from Equation (7).Specifically, persistence probabilities such as The concept of persistence weight, denoted as P d f i , plays a crucial role in the debris flow simulations by considering the position of the inflow cell.For a given incoming cell, we derive the forward and oblique direction numbers from Table 1 and then obtain their persistence probability from Equation (7).Specifically, persistence probabilities such as w 0 and w 60 represent the likelihood of the debris flow continuing in the forward direction (0 • ) and in an oblique direction (60 • ) relative to the inflow trajectory, respectively.This framework provides a nuanced understanding of how the debris flow momentum affects its trajectory.Furthermore, persistence probabilities at angles of 120 • and 180 • from the inflow direction (w 120 and w 180 ) are practically assigned a persistence probability of zero.This decision is based on the observation that debris flows are unlikely to disperse in the reverse speed direction due to the momentum and physical characteristics of the flow.Figure 4b,c depict scenarios where the inflow cells of a SCA model are edge neighbors and corner neighbors, respectively.In our previous research on SCA models [42], we suggested the values ω 0 = 0.42, ω 45 = 0.29, and ω 90 = 0, which yielded optimal results in the inversion scenario.Building on this successful value scheme, we computed the unit angle weight for each region and implemented it in the HCA models.This resulted in a persistence weight scheme applicable to HCA models, which is ω 0 = 0.52, ω 60 = 0.24.

Probability Combination
In Sections 2.2.1 and 2.2.2, this study introduced terrain probability and persistence probability into the HCA model to further determine the flow direction function of the debris flow.There are two mainstream coupling methods, which are the multiplication method [37] in Equation ( 8) and the addition method [42] in Equation (9).
Using the multiplication method of Equation ( 8), when the terrain and flow inertia are combined probabilistically-when any factor produces a probability of 0-the comprehensive flow probability of the cell unit will also be zero.Under the complex interaction of speed and topography, reverse slope flow and reverse velocity flow will occur in debris flows under realistic physical systems.However, this probabilistic combination method considers it impossible to flow in these directions, resulting in the excessive concentration of flow directions and deviating from the behavioral characteristics of a real-world debris flow.
The addition method of Equation ( 9) overcomes the above problems.Han [42] discussed the impact of different combination methods on model accuracy, and the experimental results showed that the additive method is better than the multiplicative method in model accuracy.

The Sink-Filling Approach
The sink-filling approach, utilized to simulate the traversal of troughs or dams, was first introduced by Gamma [36] and later refined by Wichmann [51].It aims to tackle the challenge posed by negative slopes that inherently impede debris flow movement.Primarily employed in geographic information systems (GIS) and terrain analysis, the sink-filling approach improves the simulation accuracy of debris flow characteristics.In natural settings, the debris flow typically flows towards the lowest points in the terrain.However, encountering a groove characterized by surrounding cells at elevations higher than that of the current cell creates a negative slope scenario that impedes debris flow movement.The sink-filling models are designed to overcome these barriers by "filling" grooves.
When the debris flow reaches a particular cell, and all neighboring cells are found to have negative slopes, the sink-filling model is activated.Debris flows cease movement and deposit sediment here, altering the elevation of the cell.The sedimentation-filling effect within the "groove" enables it to overcome the obstacle after numerous Monte Carlo iterations.This approach is primarily viewed as a computational estimate, employed in scene modeling and aiding further model iterations, thereby facilitating the debris flow to enact the over-dam accumulation process ahead of the sediment dam and enabling the exploration of debris flow pathways.

Maximum Length Function
As the debris flow transitions from steep slopes to gentle alluvial flats, its velocity progressively diminishes.Upon the velocity reaching zero, the movement of the debris flow ceases.The maximum length of debris flow accumulation is an important parameter for predicting the range of debris flow.Wang [52]'s debris flow accumulation experiment showed that when the debris flow concentration (C Vt ) and the downstream flow path slope (β) remain unchanged, the maximum length of the debris flow (L) is related to the impact of the debris flow.It is proportional to the total runoff (V), and the empirical formula form of the maximum length of debris flow is given, as shown in Equation (10).
where λ and µ are the parameters related to the concentration of the debris flow.The recommended values for the regression results of the debris flow experiments are λ = 0.42 and µ = 0.935.
When the debris flow concentration changes, this formula does not describe the maximum length accurately enough.For this purpose, he conducted a series of flume tests.For different debris flow concentrations C Vt and the debris flow outflow volume V, the maximum length formula considering the concentration factor is obtained through Equation (11).When the debris flow concentration is known, it is recommended to use Equation (11), otherwise, Equation ( 10) can be used.
In the realm of HCA models, the stopping of debris flows emerges from a nuanced interplay between static terrain factors and dynamic momentum indicators.The terrain aspect scrutinizes each potential flow direction for the presence of a 'sink-filling' scenario, which triggers the stopping conditions and subsequent deposition of debris flows.Dynamically, the focus shifts to the length of the debris flow's path, particularly whether it has reached its predefined maximum.Achieving this milestone prompts the halt of the debris flow.The number of iterations of the debris flow direction function, which is denoted as step max , is determined by the maximum length.Based on the maximal length L of the debris flow, the maximum number of iterations of the debris flow direction function for the CA models is determined according to Equation (12).
The coefficient ξ serves as a crucial correction factor for the maximum length L, as delineated in Equations ( 10) and (11), compensating for the empirical derivation of L. Our findings suggest that the integration of ξ, particularly set at 1.2, significantly enhances the simulation's accuracy.Additionally, Cellsize denotes the dimension of each hexagonal cell in the model.The simulation of the debris flow path is designed to cease once the iteration count reaches or exceeds the predefined maximum number of steps, step max , thereby delineating a comprehensive and plausible debris flow trajectory.

Path Simulation Function Based on Monte Carlo Iterative
The Monte Carlo iteration method stands as a pivotal computational technique that leverages random sampling to approximate solutions for intricate problems.At its core, this method employs a random number generator to produce an extensive series of random samples.These samples serve as a basis for inferring the potential outcomes of complex phenomena.Following the generation of a substantial dataset of random numbers, the derived potential outcomes are subjected to rigorous statistical analysis.
Using the flow direction function, we determined the probability of the debris flow to each neighboring cell through Equation ( 13), along with the random range P(j) of the neighbor cell through Equation ( 14): P(j) ϵ ( S(j − 1) , S(j) ) where S(j) denotes the maximum probability of the flow direction for the neighboring cell j, P(j) represents the interval of the flow direction probabilities for the neighboring cell j.When P i + = 0, it follows that S(j) = S(j − 1); under these conditions, P(j) for the neighboring cell j becomes an empty set, indicating that the debris flow will not enter this cell.Utilizing the Monte Carlo iterative algorithm principle, a random number n ϵ (0 , 1] is generated.Should n lie within P(i) 's interval, the subsequent flow direction of the debris is designated to neighbor cell i.
Figure 3b illustrates the debris flow moving from cell No. 1 into the central cell.Known cell elevations enable the determination of the slope and terrain probability P t f i for each neighboring cell via the terrain probability function.The inflow is identified as originating from cell No. 1.The persistence probability function reveals that P d f 4 = w 0 , P d f 3 = P d f 5 = w 60 .Subsequently, the aggregate flow direction probability for each cell is calculated through the additive method.
During the simulation, Figure 3b demonstrates the calculation of the comprehensive flow probabilities for each neighboring cell of the central cell as follows: P 1 + = P 2 + = P 6 + = 0, P 3 + = 0.43, P 4 + = 0.45, P 5 + = 0.12.Applying the random range method, the probability intervals are determined as P(1) = P(2) = P(6) = ∅, P(3) ϵ (0 , 0.43], P(4) ϵ (0.43 , 0.88], P(5) ϵ (0.88 , 1].The magnitude of the random number n dictates the subsequent direction of the debris flow.Integrating the debris flow stopping conditions in Sections 2.3 and 2.4 to ascertain the termination point of the debris flow, a complete path of the debris flow was delineated.Figure 3a describes the start, generation, and end process of such a path.Based on the Monte Carlo iterative algorithm, all possible paths of the debris flow are obtained, and the spatial set of all paths is the flow range of the debris flow.

Yohutagawa Debris Flow
Yohutagawa is located on Amami Oshima in the southern part of Kyushu Island, Japan, at the coordinates 28 • 24 ′ N, 129 • 32 ′ E. It covers an area of approximately 0.24 km 2 , with elevations ranging from 20 to 250 m.The region is predominantly forested.On 20 October 2010, the area was struck by a torrential rainstorm, reaching a peak rainfall intensity of 131 mm/h.The heavy rainfall triggered a rockslide at the highest point of the area, with a volume of 5843 m 3 .A notable characteristic of this debris flow was the relatively narrow flow channel compared to the initial collapse area, with a clear demarcation between the erosion and deposition zones.The average slope of the incised river channel in Yohutagawa is about 16 • .The total length of the debris flow was approximately 750 m. Figure 5 illustrates the different zones of the debris flow in Yohutagawa.
The CA model is based on the DEM data of the debris flow watershed, the location of the material source points, and the total amount of debris flow collapse.It rapidly predicts the potential covered area of the debris flow.The size of the cells in a SCA model is determined by the side length D Square in the x and y directions, whereas for a hexagonal cell automaton, the cell size is defined by the center-to-center distance D Hex between adjacent units.In this study, D Square is set at 2.5 m based on the recommended values [37,43].Additionally, this paper uses D Hex = 2.5 m to conduct a comparative analysis between the HCA and SCA models.Based on the DEM data of the Yohutagawa event (1425 m × 1138 m), this paper performed interpolation processing on the terrain data of 455 rows and 570 columns under the square grid to obtain the hexagonal grid terrain data of 455 rows and 658 columns.
The simulation validation of the Yohutagawa debris flow extent was conducted using both the HCA and SCA models.The parameters for these models are shown in Table 2.The predicted extents of the debris flow from both models are illustrated in Figures 6 and 7.Both models performed well in simulating the upper and middle reaches of the channel, with the primary differences observed in the downstream deposition area.
area, with a volume of 5843 m 3 .A notable characteristic of this debris flow was the relatively narrow flow channel compared to the initial collapse area, with a clear demarcation between the erosion and deposition zones.The average slope of the incised river channel in Yohutagawa is about 16°.The total length of the debris flow was approximately 750 m. Figure 5 illustrates the different zones of the debris flow in Yohutagawa.The CA model is based on the DEM data of the debris flow watershed, the location of the material source points, and the total amount of debris flow collapse.It rapidly predicts the potential covered area of the debris flow.The size of the cells in a SCA model is determined by the side length  in the x and y directions, whereas for a hexagonal cell automaton, the cell size is defined by the center-to-center distance  between adjacent units.In this study,  is set at 2.5 m based on the recommended values [37,43].Additionally, this paper uses  = 2.5 m to conduct a comparative analysis between the HCA and SCA models.Based on the DEM data of the Yohutagawa event (1425 m × 1138 m), this paper performed interpolation processing on the terrain data of 455 rows and 570 columns under the square grid to obtain the hexagonal grid terrain data of 455 rows and 658 columns.
The simulation validation of the Yohutagawa debris flow extent was conducted using both the HCA and SCA models.The parameters for these models are shown in Table 2.The predicted extents of the debris flow from both models are illustrated in Figures 6 and  7.Both models performed well in simulating the upper and middle reaches of the channel, with the primary differences observed in the downstream deposition area.

Flume Test
The objective of the debris flow flume test is to explore the initiation, flow dynamics, and erosion-deposition processes related to debris flows.By simulating variables such as

Flume Test
The objective of the debris flow flume test is to explore the initiation, flow dynamics, and erosion-deposition processes related to debris flows.By simulating variables such as diverse terrain configurations, soil-water mixture ratios, and flow velocities, we can gain insights into dynamic behavior and impacts of debris flows, providing a scientific basis for predicting, mitigating, and managing them.Figure 8a represent the test equipment used by Liu et al. [53].And Figure 8b depicts the simulation result of debris flow.A detailed analysis will be provided in the following section.

Result Analysis
Precision (Pr) is defined as the ratio of the true-positive samples to the total number of samples identified as positive by the model by the model as Equation (15).Recall (Re) is quantified as the ratio of the true-positive samples to the entire set of samples that are genuinely positive as Equation (16).
TP (true positives): The regions that are accurately identified as being subject to debris flows are termed 'overlap areas'.FP (false positives): Areas erroneously predicted as debris flow zones are designated as 'misjudgment areas'.FN (false negatives): Areas inaccurately classified as non-debris flow zones are referred to as 'missed judgment areas.Table 3 details the calculation results of the TP, FP, and FN ranges for the two cases.diverse terrain configurations, soil-water mixture ratios, and flow velocities, we can gain insights into the dynamic behavior and impacts of debris flows, providing a scientific basis for predicting, mitigating, and managing them.Figure 8a represent the test equipment used by Liu et al. [53].And Figure 8b depicts the simulation result of debris flow.A detailed analysis will be provided in the following section.

Result Analysis
Precision () is defined as the ratio of the true-positive samples to the total number of samples identified as positive by the model by the model as Equation (15).Recall () is quantified as the ratio of the true-positive samples to the entire set of samples that are genuinely positive as Equation (16).An evaluation factor, represented by Ω in Equation ( 17), was employed to comprehensively evaluate the model's performance, taking into account both misjudgment and missed judgment [43].
The evaluation metric Ω systematically integrates the overlapping area, misjudgment area, and missed judgment area within the debris flow range prediction, providing a wellstructured quantification approach.However, this metric assumes a linear influence of the three parameters α, β, and γ on the assessment of predictive accuracy for debris flow ranges, thereby neglecting the differing impacts of misjudgments and missed judgments on the detrimental outcomes associated with debris flow range predictions.
Zhang [54] employs the F1 score, as depicted in Equation ( 18), to evaluate model simulation accuracy, defined as the harmonic mean of precision and recall.To account for the heightened adverse impact of the missed judgments over misjudgments, this study adopts the F-beta score.The F-beta score, presented in Equation ( 19), serves as a performance evaluation metric that reflects the model's precision and recall rates.It adjusts the recall rate via a parameter β, thus modulating the significance of recall relative to precision.For debris flows, β = 2 is selected to emphasize the importance of recall.Table 4 shows the run-out extent evaluation index in the Yohutagawa debris flow and flume test.Compared to the SCA model, the evaluation factors Ω and Fβ of the HCA model have improved.In the Yohutagawa debris flow, Ω and Fβ increased from 0.5424 and 0.8304 to 0.5642 and 0.8687, respectively, resulting in an increase in the simulation accuracy by 4.02% and 4.61%.In the flume test, Ω and Fβ increased from 0.5054 and 0.8191 to 0.5381 and 0.8488, respectively, leading to a simulation accuracy improvement of 6.47% and 3.63%.However, it is noteworthy that the evaluation factor F1 has not shown significant improvement, and in some cases, the effect in the Yohutagawa debris flow was not as favorable as the original SCA model.The reason for this occurrence lies in the inherent conservative nature of the HCA model.To minimize missed judgments, it tends to prioritize overestimating the debris flow extent, which inadvertently results in an increase in misjudgments.

Interpolation Accuracy
In executing the slope computations within the flow direction function delineated in Section 2.2, acquiring elevation data for each cell is imperative.This study employs DEM raster data to discretely depict the Earth's topography, furnishing elevation data for a finite set of known points.Within the HCA framework of the network, the position information of all cells is ascertainable via the lattice division method; however, the elevation for these cells remains undetermined.To deduce the elevation values of these indeterminate points, selecting a suitable interpolation method becomes essential.
Based on the first law of geography, widely employed terrain interpolation techniques encompass inverse distance weighting (IDW), bilinear interpolation, and cubic interpolation.This study evaluates the interpolation accuracy of the aforementioned terrain interpolation methods, with the objective of identifying the most appropriate technique for mountainous terrain for future research and application.
The "peaks" function in MATLAB exemplifies a tool for mountain terrain simulations, generating landscapes with notable peaks and valleys.By adjusting the parameters, researchers can customize the terrain features to align with their study's topographic criteria, obtaining DEM raster data, which includes elevation values at the grid points.We compared terrain generated by the bilinear, cubic, and IDW interpolation methods.The bilinear and cubic interpolations produced smoother terrains, whereas IDW appeared to have jagged features.DEM accuracy is assessed using the root mean squared error (RMSE).The calculation method of RMSE is shown in Equation (20).
where Z is the true value of the terrain, z is the observed value or calculated value, ε = z − Z is the error, and n is the number of errors.
Better smoothness can be obtained by interpolating the terrain using the bilinear interpolation method or the cubic interpolation method.When mainly considering the interpolation accuracy, it is recommended to use the cubic interpolation method.This conclusion is based on a MATLAB terrain interpolation experiment, which compared the computational efficiency and accuracy of different interpolation methods.The cubic interpolation method proved superior to the bilinear interpolation method in terms of interpolation accuracy, as reflected in Table 5.The RMSE in the cubic interpolation was an order of magnitude lower, but the results of both methods were within acceptable limits.The inverse distance weighting interpolation method was excluded first due to its low accuracy.The bilinear interpolation method has an efficiency advantage over the cubic interpolation method, especially when high accuracy in terrain interpolation is required.The efficiency analysis of the interpolation methods in this paper is based on the terrain conditions of the Yohutagawa debris flow event, which covered a small area.Even with high-precision interpolation, the increase in computation time was only about two minutes.This is acceptable for an order of magnitude improvement in RMSE precision, so the cubic interpolation method was chosen.

Interpolation Efficiency
Interpolation calculations were conducted on the Yohutagawa DEM data, with the calculation durations for each of the two interpolation methods recorded.A comparative analysis revealed that across varying lattice sizes, the bilinear interpolation method consistently necessitates less time than the cubic interpolation method.Specifically, with a cell size of 3 m, the efficiency of bilinear interpolation significantly surpasses that of cubic interpolation.The interpolation duration decreased from 345.85 s to 245.69 s, marking a reduction of 100.16 s and an enhancement in computational efficiency by 28.96%.However, at a cell size of 5 m, the computational efficiency advantage of the bilinear interpolation method diminishes, with the reduction in interpolation time amounting to merely about 10 s.For cell sizes of 7 m and 10 m, the time savings achieved through interpolation are even more marginal, amounting to only 1-2 s.In scenarios with fine interpolation grids and extensive computational demands, the bilinear interpolation method markedly enhances the calculation speed.Consequently, the bilinear interpolation method is advisable.Nevertheless, considering the relatively short overall interpolation duration and the minimal time differences between commonly utilized lattice sizes (5 m, 10 m), cubic interpolation, which offers superior accuracy, is recommended.This approach is also employed in the method presented in this study for interpolating cell elevations.

Dam-Crossing Test
In Sections 2.3 and 2.4, we introduced the sink-filling model and the maximum length function, respectively.This section conducts a series of model tests to verify the performance of these two models in the HCA, including their advantages and necessity.The tests simulate debris flows of various sizes in the Yohutagawa terrain, originating from the collapse point and flowing downstream to the dam. Figure 9a displays the results when the simulated flushing volume is 1000 m 3 .At this point, the debris flow fails to cross the dam.In Figure 9b, with an overflow volume of 2000 m 3 , the debris flow begins to cross the dam body, posing risks to human habitation areas.Figure 9c depicts the scenario where the flushing volume is set at 3000 m 3 .The simulation results show that the debris flow successfully traverses the dam body.As the outflow reaches 5000 m 3 , shown in Figure 9d, more debris accumulates on the alluvial fan.threshold capacity, in the HCA, this is depicted as a breach volume adequate for complete "dam-front filling," with the surplus flow overflowing the "leveled" dam.Subsequently, the maximum length function determines the debris flow's deposition pattern and extent of the flow.At this juncture, as the debris flow scale increases, so does its maximum length (Equation ( 10)), signifying its impact on more distant areas.However, a notable aspect is the triggered lateral expansion.With the escalation in breach volume, compared to Figure 9c, Figure 9d exhibits an elongation in breach length, and the debris flow's "width" becomes more pronounced.The experimental findings illustrate that the HCA model offers an appropriate depiction of this phenomenon and even conservatively predicts it to some degree.

Conclusions
This study presents a HCA numerical model specifically designed for simulating debris flows.In this framework, digital elevation model data are used to divide the terrain into hexagonal cell spaces and perform elevation interpolation.After debris flow simulations, the case inversion analysis confirms the model's suitability.
1.The hexagonal honeycomb network configuration closely resembles a circular shape compared to quadrilateral grid structures and demonstrates isotropy, guaranteeing uniform properties in all directions.This feature is beneficial for maintaining the This experiment corresponds to the stopping conditions following debris flows of varying scales.When the debris flow scale falls within the dam's capacity, it is entirely obstructed, posing no threat to the human habitation facilities downstream.In our HCA model, such scenarios trigger a sink-filling model, where all debris flow sediments accumulate at the site, raising the terrain elevation to facilitate the potential overtopping of the dam by subsequent debris flows.Once the debris flow scale surpasses the dam's threshold capacity, in the HCA, this is depicted as a breach volume adequate for complete "dam-front filling", with the surplus flow overflowing the "leveled" dam.Subsequently, the maximum length function determines the debris flow's deposition pattern and extent of the flow.At this juncture, as the debris flow scale increases, so does its maximum length (Equation ( 10)), signifying its impact on more distant areas.However, a notable aspect is the triggered lateral expansion.With the escalation in breach volume, compared to Figure 9c, Figure 9d exhibits an elongation in breach length, and the debris flow's "width" becomes more pronounced.
The experimental findings illustrate that the HCA model offers an appropriate depiction of this phenomenon and even conservatively predicts it to some degree.

Conclusions
This study presents a HCA numerical model specifically designed for simulating debris flows.In this framework, digital elevation model data are used to divide the terrain into hexagonal cell spaces and perform elevation interpolation.After debris flow simulations, the case inversion analysis confirms the model's suitability.

1.
The honeycomb network configuration closely resembles a circular shape compared to quadrilateral grid structures and demonstrates isotropy, guaranteeing uniform properties in all directions.This feature is beneficial for maintaining the model's geometric and physical coherence during the simulation of complex terrain and flow dynamics.2.
We compared the accuracy and efficiency of the IDW, bilinear, and cubic interpolation methods.The results show that cubic interpolation has the highest interpolation accuracy, and IDW interpolation accuracy is poor.When the lattice is finely divided, bilinear interpolation has obvious efficiency advantages over cubic interpolation.We recommend using bilinear or cubic interpolation methods as appropriate in the specific case.

3.
Building upon cellular automaton theory, transition rules define the process of debris flow.These transition rules focus on interactions among neighboring cells, avoiding the necessity of solving complex partial differential equations.The model utilizes the flow direction function, the sink-filling approach, the maximum length function, and the Monte Carlo iterative method to simulate the debris flow run-out extent, highlighting the model's usefulness.4.
The model has been applied to both simulate a water tank and replicate the Yohutagawa debris flow incident.

Figure 1 .
Figure 1.Hexagonal cells, lattice form: (a) cell neighbor relations and distance, (b) cell orientation and size, and (c) lattice form in offset coordinate system.

Figure 2 .
Figure 2. Neighbor distance in different situations: (a) SCA model using the Moore neighborhood, (b) HCA model using the hexagonal neighborhood.

Figure 1 .
Figure 1.Hexagonal cells, lattice form: (a) cell neighbor relations and distance, (b) cell orientation and size, and (c) lattice form in offset coordinate system.

Figure 1 .
Figure 1.Hexagonal cells, lattice form: (a) cell neighbor relations and distance, (b) cell orientation and size, and (c) lattice form in offset coordinate system.

Figure 2 .
Figure 2. Neighbor distance in different situations: (a) SCA model using the Moore neighborhood, (b) HCA model using the hexagonal neighborhood.

Figure 2 .
Figure 2. Neighbor distance in different situations: (a) SCA model using the Moore neighborhood, (b) HCA model using the hexagonal neighborhood.

Figure 4 .
Figure 4. Relationship between persistence probability and inflow direction: (a) the HCA model, (b) edge neighbor inflow in the SCA model, and (c) corner neighbor inflow in the SCA model.

Figure 4 .
Figure 4. Relationship between persistence probability and inflow direction: (a) the HCA model, (b) edge neighbor inflow in the SCA model, and (c) corner neighbor inflow in the SCA model.

Figure 5 .
Figure 5. Yohutagawa debris flow: (a) Amami Oshima location, (b) debris flow source area, (c) channel of the Yohutagawa debris flow, (d) dam downstream of the channel, and (e) observation area.

Figure 5 .
Figure 5. Yohutagawa debris flow: (a) Amami Oshima location, (b) debris flow source area, (c) channel of the Yohutagawa debris flow, (d) dam downstream of the channel, and (e) observation area.

Figure 6 .
Figure 6.SCA model run-out extent result.Figure 6. SCA model run-out extent result.

Figure 6 .
Figure 6.SCA model run-out extent result.Figure 6. SCA model run-out extent result.

Figure 8 .
Figure 8.(a) the test equipment used by Liu et al., (b) Comparison of the simulated area with the real area in flume test.

Figure 8 .
Figure 8.(a) the test equipment used by Liu et al., (b) Comparison of the simulated area with the real area in flume test.

Figure 9 .
Figure 9. Model results for debris flows of different scales: (a) Sedimentation in front of the dam, (b) debris flow over the dam, (c) a large number of debris flows rushed out and reached their maximum length, and (d) more outbursts strengthen the lateral expansion of debris flows.

Figure 9 .
Figure 9. Model results for debris flows of different scales: (a) Sedimentation in front of the dam, (b) debris flow over the dam, (c) a large number of debris flows rushed out and reached their maximum length, and (d) more outbursts strengthen the lateral expansion of debris flows.

Table 1 .
Forward and oblique directions under different incoming cells.

Table 2 .
The parameters for CA models.
Water 2024, 16, x FOR PEER REVIEW 11 of 19

Table 3 .
Areas of each area in the Yohutagawa debris flow and flume test.

Table 4 .
Evaluation index in the Yohutagawa debris flow and flume test.

Table 5 .
RMSE of different interpolation methods.
In the Ω-based evaluation framework, the Ω values for the flume test simulation and the Yohutagawa debris flow incident are 0.5381 and 0.5642, respectively.These values represent improvements over the previous CA model's values of 0.5054 and 0.5424.The Fβ score evaluation framework considers the varying impacts of missed judgments and misclassifications.The Fβ scores for the HCA model, 0.8488 and 0.8687, respectively, surpass those of the SCA model (0.8191 and 0. 8304) by 3.63% and 4.41%, respectively.Z.H. led the research program.Y.L. and N.J. designed the simulation.Q.F. and X.Z.wrote the manuscript.Z.H. and Y.M. reviewed and edited the manuscript.All authors have read and agreed to the published version of the manuscript.This study was financially supported by the National Natural Science Foundation of China (Grant No. 52078493), the Natural Science Foundation of Hunan Province (Grant No. 2022JJ30700), the Natural Science Foundation for Excellent Young Scholars of Hunan (Grant No. 2021JJ20057), the Science and Technology Plan Project of Changsha (Grant No. kq2305006), and the Innovation Driven Program of Central South University (Grant No. 2023CXQD033).These financial supports are gratefully acknowledged.
Author Contributions:Funding: