## 1 Introduction  
Landslides are common geological disasters in mountainous areas that pose serious threats to human life, property, and the natural environment [1; 2; 3]. Landslides can involve complex evolutions of dynamic characteristics of the mobile mass, such as the traveled distance, velocity, and extent of affected areas. Evaluation of the runout dynamic characteristics of landslides has significant geotechnical significance, especially in landslide hazard assessment or identifying potential mitigation measures [4; 5; 6; 7]. Various methods have been proposed for landslide runout analysis, which can be categorized into empirical-statistical and analytical methods [3]. The empirical-statistical methods are typically based on observations to establish correlations between the characteristics of the source area and the runout behaviors (e.g., [8; 9; 10]). Analytical methods are mainly based on the analysis of landslide physical processes and can be categorized into numerical (dynamic) methods [1; 11; 12; 13] and closed analytical equations. The application of the statistical method is limited by the representativeness of source area features, and the accuracy of the method is affected by the amount of statistical data. In contrast, the dynamic methods, which are based on physical process analysis, can better reproduce the real process of landslides via numerical simulation and provide more dynamic information about landslides [3; 6]. Dynamic models with simplified rigid body assumptions have been proposed early, such as the mass point method (or lumped mass model) [14; 15] or the energy-based approach of [PERSON] [16]. All these methods assume that the entire mass is concentrated at the centroid of the landslide body. During the analysis, the internal deformation of the sliding mass will be ignored, so these models are more suitable for analyzing landslides with high integrity. After that, methods using the deformable hypothesis were proposed. In these methods, a sliding body will be discretized into deformable blocks adjacent to each other [17] or rigid bodies connected by springs, through which the deformation of the block can be approximated equivalently [18; 19]. These methods are suitable for analyzing slopes with high initial integrity and whose integrity has not experienced a significant decline during sliding. In the field of geophysical mass flow, [PERSON] and [PERSON] [20] proposed the first depth-averaged model for analyzing the avalanches of dry granular mass. Its extended model takes into account the influence of complex terrain and other factors [21; 22]. The depth-averaged dynamic models have been extended to solid-fluid two-phase flow [1; 12; 13], two-layer flow [13; 23], two-fluid flow [24], and multiphase flow [25]. For dry landslides with low integrity or high mobility landslides with high water content, the models based on hydrodynamics are more suitable. Most depth-averaged models are based on hydrodynamics [26]. They ignore the stiffness of sliding materials and assume that landslides become fluidized instantly when failure occurs [27]. Some landslides move like a rigid block at the initial stage and then disintegrate into a flow-like debris avalanche progressively due to friction and vibration [28; 29]. In these cases, the application of the depth-averaged model often results in an excessive impact area at the initial stage [3; 30; 31; 32]. In order to obtain results closer to the field within the same dynamic model, some researchers (e.g., [33]) used zone-wised parameter values at the initiation, avalanche, and accumulation zones. However, due to the limitations of the model's own characteristics, the differences in dynamic characteristics that can be reflected by adjusting parameter values are limited. To overcome this issue, [PERSON] and [PERSON] [32] developed a modified version of DAN3D that allows for delaying fluidization based on a user-specified index.