The TiVaSS model developed by An et al. (2016) analyzes slope stability based on the infinite slope-stability model as shown in Figure 1, and the shear and normal stress at the slope surface are calculated as follows: τ f s = T b / cos ⁡ ϕ = W b cos ⁡ ϕ ⁡ sin ⁡ ϕ , (1) σ = P b / cos ⁡ ϕ = W b cos 2 ⁡ ϕ , (2) where W = γ s D b is the soil weight (kg/m) and γ s is the unit weight of the soil (kg/m³), D is the soil depth (m), and b is the slope width (m); T is the shear force (kg/m); P is the normal force (kg/m); and ϕ is the slope angle (radian). According to the Mohr–Coulomb theory, the shear strength at an infinite slope is calculated as S = c + σ ′ ⁡ tan ⁡ φ , (3) where c is the cohesion (kg/m²); σ ′ = σ − δ a + σ s is the effective stress with excessive suction force (kg/m²), which was generalized by Lu and Likos (2006) and δ a is the air pressure and σ s = − δ a − δ w S e is the suction stress, where δ w = γ w ψ is the water pressure, γ w is the unit weight of water (kg/m³), and ψ is the pressure head of the subsurface water (m); and S e = θ − θ r θ s − θ r is the effective saturation and θ is the volumetric moisture content, where θ s is the saturated moisture content, and θ r is the residual moisture content; and φ is the internal friction angle (degrees). Finally, time-varying slope stability proposed by Iverson (2000) F S is defined as follows: F S = S τ f s = ψ S e γ w ⁡ tan ⁡ φ γ s D ⁡ sin ⁡ ϕ ⁡ cos ⁡ ϕ + tan ⁡ φ tan ⁡ ϕ + c γ s D ⁡ sin ⁡ ϕ ⁡ cos ⁡ ϕ . (4)

The only time-variant variable in Eq. 4 is the pressure head because effective saturation is usually given as a function of the pressure head. Moreover, in the TiVaSS model, the subsurface flow in saturated soil due to rainfall is interpreted through the 3D Richard equation, which is as follows: ∂ θ ψ ∂ t − ∇ ⋅ K ψ ∇ ψ + z − q = 0 , (5) where K is the hydraulic conductivity (m/s); t is time (s); z is the vertical dimension, which is assumed to be positive in the upward direction; and q is a general source term (1/s), including rainfall.

In the TiVaSS model, Eq. 5 is discretized within FVM, and the following equation is applied with the Gauss–Green divergence theorem (An and Yu, 2014): ∫ V ∂ θ ∂ t d V − ∫ ∂ V n ⋅ K ∇ ψ + z d ∂ V − ∫ V q d V = 0 (6) where V is the control volume; ∂ V is the control-volume boundary; and θ and ψ are assumed to be the cell-averaged values from the finite-volume approximation. The detailed equation from this model was given by An et al. (2016).

Flow can efficiently be analyzed if it has a small vertical height relative to the area of the horizontal surface, such as a debris flow, using the shallow-water equation based on the Navier–Stokes equation. The Deb2D model, developed by An et al. (2019), analyzes debris flow using a two-dimensional shallow-water equation based on a rectangular grid that utilizes an adaptive mesh-refinement technique. Therefore, this numerical model requires a shorter time than that required by other models to calculate the flow state. The shallow-water equation is as follows: ∂ q ∂ t + ∂ f ∂ x + ∂ g ∂ y = s , (7) where t denotes time; x and y are Cartesian coordinates; and q , f , g , and s are vectors representing conserved variables, fluxes in the x and y directions, and source terms, respectively. The vectors can be written as q = h h u h v , f = h u h u 2 + g h 2 / 2 h u v , g = h v h u v h v 2 + g h 2 / 2 , s = V C S g x − S f x S g y − S f y , (8) where h is the depth of the debris-flow mixture; u and v are the depth-averaged velocity components in the x and y directions, respectively; g is the acceleration of gravity; V c is the volume change of the debris flow mixture; S g x and S g y represent the gravitational acceleration in the x and y directions, respectively; and S f x and S f y represent the driving friction in the x and y directions, respectively. Here, V c and S f are factors for realistically simulating debris flow in numerical analysis.

Debris flows cause various interactions with the ground surface during the flow process, with erosion, entrainment, and deposition being essential mechanisms for simulating the debris flow in a numerical model. Calculating S f is required to implement these processes physically. This study considers the influence of vegetation in the debris flow, and S f , which reflects this, is calculated as follows: S f = S f r + S f v (9) where S f r is commonly estimated as a non-Newtonian fluid in simulated mixtures, such as debris flows, and S f v is the vegetation drag force. First, the Voellmy rheology for calculating S f r in this study is as follows. S f r x = u u μ g h + g u 2 + v 2 ξ , S f r y = v v μ g h + g u 2 + v 2 ξ , (10) where μ and ξ , the Coulomb friction and turbulent friction coefficients, dominate the deceleration behavior when the flow is slow and fast, respectively (Bartelt et al., 2013; Frank et al., 2015). Because these parameters reflect the field indirectly rather than directly, researchers use calibrated values by performing back-analysis. S f v , which reflects the effect of vegetation, is developed as follows (Lee et al., 2004; Hui et al., 2010; Tang et al., 2014; Kang et al., 2022): S f v x = 1 2 n s D s A s C v h v u u 2 + v 2 , S f v y = 1 2 n s D s A s C v h v v u 2 + v 2 , (11) where n s is the number of plants per unit area; D s is the average vegetation diameter; A s is the unit area; C v is the vegetation drag coefficient; and h v is the maximum vegetation height, which is calculated from h v = min h , h v . The mechanism behind these three processes is shown in Figure 2, and the algorithm proposed by Lee et al. (2022b) is as follows: V C x , y , t = d z d t e i f )   τ d f > τ e − d z d t d e l s e   i f )   τ d f < τ d (12) h max x , y , 0 = d z d τ τ d f − τ e = d z d τ ρ g h s − τ e i f )   τ d f > τ e 0 e l s e (13) where d z / d t e is the constant erosion-entrainment rate; d z / d t d is the constant deposition rate; h max x , y , 0 is the maximum potential erosion depth; τ d f is the shear stress of debris flow; τ e and τ d are the critical shear stress of erosion and deposition, respectively; d z / d τ is the average potential erosion depth, ρ is the mass density; and s is the channel slope. An et al. (2019) and Lee et al. (2022b) developed a detailed equation of this model.

In this study, Eq. 11 was adopted to consider the effect of resistance caused by vegetation in the debris flow process. However, according to Julian and Torres (2006) and Shen et al. (2017), vegetation not only affects the resistance in the flow process but also affects the erosion process in the soil. Shen et al. (2017) analyzed that the critical shear stress can increase 1.3–2.4 times in the presence of vegetation. Therefore, in this study, the increase in critical shear stress derived from a previous study was considered by assuming that τ e increase 1.3 and 2.4 times when vegetation is present.

This study simulated the landslide-debris flow event induced by Typhoon Ewiniar in 2006 at the Doam dam basin. However, due to the large scale of the target basin, which covers 144.9 km², an initial analysis of overall slope stability was performed using a 30 m × 30 m resolution DEM, followed by a detailed analysis using a 10 m × 10 m resolution DEM for high-risk areas. The slope stability analyzed for various scenarios of vegetation distribution via the TiVaSS model and low-resolution DEM is shown in Figure 9.

The overall slope stability of the Doam basin analyzed using a low-resolution DEM revealed that the slope would collapse at approximately the 12th hr on the 16th, with many slope failures predicted to occur in the lower part where the Doam dam and lake are located. In scenario 1, which did not consider vegetation, many slopes were analyzed as unstable, particularly near Doam Lake. However, in scenario 2, which considered vegetation, the slopes were analyzed as relatively stable. For precision, we used a high-resolution DEM to perform a detailed analysis of the lower part of the basin, where the areas on the verge of instability were observed.

As shown in Figure 10, the simulation results of the high-resolution DEM indicated that slope failures occurred twice on the lower part of the Doam dam basin at the 10th and 12th hr on the 16th. The collapse occurred near the left bank of the river at the 10th hr and in areas b1 and b2 near the right bank and Doam Lake at the 12 h. In scenario 1, the unstable slopes rapidly increased from the 10th to 12th hr on the 16th. In scenario 2, the influence of vegetation was noticeable, and the slopes previously identified as high-risk areas in scenario 1 were largely stable until the 10th hr. Further, some unstable slopes were observed at the 12th hr. This simulation results reveal that vegetation not only delays the occurrence of slope failures but also reduces their magnitude. Therefore, if vegetation is ignored in the slope stability analysis of areas where vegetation is sufficiently distributed, the risk of slope failures will be overestimated. As illustrated in Figure 6, this phenomenon seems to be due to the increasing cohesion (in Eq. 3) by considering the root reinforcement. As a result, vegetation increased the slope stability (in Eq. 4). The quantitative analysis of these results is summarized in Table 2. For the Doam dam basin, if vegetation were ignored, 297,540 m³ (=2.45 times) of additional slope failures would occur, which is an overestimation. Based on the simulation results shown in Figure 10, this study assumed that debris flow occurred in unstable areas (FS < 1). The a and b2 indicated in Figure 10 demonstrate that TiVaSS well simulated the collapse area in Figures 7A–C, respectively. We analyzed debris flow and the impact on the Doam dam by quantitatively comparing and analyzing the debris flow influx and resulting damage to the Doam dam under different vegetation distribution scenarios.

In this study, we analyzed the influence of vegetation on the stability analysis and the debris flow analysis models. Based on the simulation results for scenarios 1 and 2 obtained using the TiVaSS model, we divided the two scenarios into three by simulating the presence and absence of vegetation in the debris flow using the Deb2D model. Figure 11 shows the flow height simulation results of the debris flow using a high-resolution DEM, and debris flowed into Doam Lake due to the collapse shown in Figure 10. Scenario 1 in Figure 11 represents the analysis of both TiVaSS and Deb2D models without considering vegetation. Scenario 2-1 presents the effects of considering vegetation in the TiVaSS model but not in the Deb2D model. Scenarios 2-2 and 2-3 represent considering vegetation in both TiVaSS and Deb2D models. We set the critical shear stress as 1.3 kPa in scenario 2-2 and 2.4 kPa in scenario 2-3. Scenarios 1 and 2.1 highlighted the importance of vegetation in slope stability analysis by comparing the amounts of sediment entering the Doam Lake. The importance of vegetation in debris flow analysis was demonstrated through scenarios 2.1 and 2.2/2.3. The debris flows under different vegetation distributions were analyzed based on scenarios 2-2 and 2-3. Scenarios 1 and 2-2/2-3 involved comprehensively evaluating the importance of vegetation in both slope stability and debris flow analysis by comparing the simulation results with and without considering vegetation. As shown in Figure 11, the flow depth of the debris flow decreased under scenario 2-1 compared to under scenario 1, and the deposition range slightly decreased. However, scenario 2.2, in which the vegetation blocked debris flow and increased the critical shear stress 1.3 fold, did not visually differ significantly from scenario 2.1. In scenario 2-3, the flow depth of the debris flow decreased significantly compared to that in other scenarios, and the deposition range decreased significantly. Vegetation effectively reduced the damage caused by landslide-debris flow events by suppressing erosion and entrainment by increasing the critical shear stress rather than simply blocking the flow. Thus, vegetation was more effective in mitigating damage by increasing the critical shear stress rather than by blocking the flow. The quantitative analysis of these results is summarized in Table 2.

As listed in Table 2, a slope collapse of 503,190 m³ resulted in 666,110 m³ of sediment inflow into the Doam Lake in scenario 1, leading to a 1.67% decrease in dam performance. In Scenario 2-1, the TiVaSS model considered vegetation in the Doam dam basin but the Deb2D model did not. Moreover, a slope collapse of 205,650 m³ resulted in 575,150 m³ of sediment inflow, causing a deterioration in function of 1.44%. In Scenario 2-2, in which the TiVaSS and Deb2D models considered vegetation, a slope collapse of 205,650 m³ resulted in 493,500 m³ of sediment inflow into the Doam Lake, causing a function deterioration of 1.23%. In Scenario 2-3, a slope collapse of 205,650 m³ resulted in a sediment inflow of 134,710 m³, causing a 0.34% function deterioration. Four scenarios in Figure 11 illustrate that the Deb2D well simulated the debris flow traces shown in Figures 7A–C. Thus, the TiVaSS and Deb2D models successfully simulated the observed damaged areas traced by satellite images.

When vegetation only was considered in the slope stability analysis, as shown in Figure 12, a substantial reduction of 297,540 m³ in slope collapse was observed compared to that in scenarios where vegetation was not considered. However, in the debris flow analysis, ignoring vegetation produced a minor difference of 90,960 m³ of sediment inflow into Doam Lake. Nevertheless, considering vegetation in both simulation processes reduced the sediment inflow into Doam Lake by 172,610–531,400 m³ compared to when vegetation was ignored in both processes, indicating that vegetation is a critical factor in the analysis of sediment inflow due to landslide-debris flow. As demonstrated in the analysis of scenarios 2-2 and 2-3, even when considering vegetation in the debris flow simulation, the mitigation effect can vary greatly depending on the distribution and type of vegetation. As listed in Table 2, the dam function deterioration caused by landslide-debris flow events can seem relatively insignificant. However, in South Korea, heavy rainfall occurs every year with typhoons, and rainfall intensity is increasing (Kim H. et al., 2021; Yeo et al., 2022). Therefore, considering vegetation is essential when analyzing sediment inflow into the dam/reservoir due to landslide-debris flow over a long period, especially in mountainous areas.

Vegetation increases slope stability and reduces the scale of slope collapse, as shown in Figure 10. Furthermore, it delays collapse and impedes flow and erosion-entrainment processes during debris flow. Therefore, vegetation effectively reduces the damage caused by landslide-debris flow events. This study addressed that increasing the critical shear stress of the topsoil was more effective than blocking the debris flow due to vegetation. It is hard to generalize these simulation results. Therefore, additional research on the role of vegetation in debris flow analysis is necessary.

In this study, TiVaSS and Deb2D models were used to analyze the impact of landslide-debris flow events on dam functions. Additionally, we constructed vegetation distribution scenarios and performed simulation analyses to evaluate the influence of vegetation. We effectively identified the deterioration effects of landslide-debris flow events on dam functions and the effects of vegetation in mitigating damage. However, physical and quantitative comparisons were somewhat limited due to insufficient observational data for validating the simulation results. Thus, collecting data to validate simulation results is necessary to analyze and understand these phenomena.

The commonly used DEM for input data is constructed based on topography, which differs from the actual features of water systems, such as rivers, lakes, and seas, as shown in Figure 13. Therefore, the method used in this study is limited in precisely analyzing the volume of the debris flow that enters lakes, and stimulating the waves generated in the lake due to debris flow is challenging. To accurately analyze the deterioration effects on the dam/reservoir functions and analyze incidents such as the 1963 Vajont dam disaster, additional research, such as that involving a two-layer debris-water system is, necessary. This could solve the problem illustrated in Figure 13.