Introduction

Debris flows are defined as massive flows of debris mixed with water, which occur when debris, such as soils, gravels, boulders, and drift-woods, is entrained by surface water flow due to heavy rainfall. Debris flows are coupled geo-hydro-mechanical phenomena that initiate on steep slopes, and involve fast flow along a confined channel, and deposition in flat areas (Jakob and Hungr, 2005). Debris flows are classified as one of the perilous landslide-related hazards due to their fast flow velocity, large impact force, and long runout, in association with poor predictability (Jakob and Hungr, 2005). In many countries, debris flows have been recorded to have caused massive damages in urban areas, e.g., in the Hongchun gully, China, on August 12–14, 2010 (Tang et al., 2012; Shen et al., 2019); in Southern Leyte Province, Philippines, on February 17, 2006 (Orense and Sapuay, 2006); and in Mt. Woomyeon, Seoul, South Korea, on July 26–27, 2011 (Yune et al., 2013). Accordingly, debris-flow barriers that can dissipate the energy of such flows have been widely implemented to mitigate potential damages.

Various physical and numerical modeling approaches have been taken for assessment of barrier performance against debris flows and optimization of barrier design. Scaled physical modeling experiments have the advantage of being repeatable; however, they also face challenges in the reproduction of debris flows that exactly simulate real-world conditions due to the difficulties in replicating rheological characteristics and entrainment processes and because of scaling effects (e.g., Ikeya and Uehara, 1980; Takahara and Matsumura, 2008; Ng et al., 2015; Choi et al., 2018). Debris flows have also been numerically modeled to identify flow characteristics using various techniques (e.g., Beguería et al., 2009; Christen et al., 2010; McDougall and Hungr, 2004, O’Brien et al., 1993; Pastor et al., 2009; Pirulli and Mangeney, 2008; Quan, 2012). For instance, Remaître et al. (2008) used 1D flow modeling based on Bingham plastic rheology to assess where the greatest reduction in debris-flow intensity would occur upon the installation of check dams. Jeong et al. (2015) used the Eulerian-Lagrangian finite element method (FEM) to examine the impact pressure of debris flows on a closed-type barrier. Table 1 summarizes the salient observations made on the barrier-debris flow interactions via numerical modeling approaches. However, previous studies have barely considered the entrainment process (i.e., surface erosion by debris flows) and have been limited to one- or two-dimensional channel flow analysis.

Table 1 Previous numerical studies on debris-flow barrier effect
Full size table

Numerical modeling of debris flows requires (i) the selection of rheological model and determination of associated parameters, (ii) consideration of entrainment process, and (iii) incorporation of 3D terrain morphology. The rheological model is used to compute the basal shear stress which operates opposite to the flow direction. Thereby, the rheological model and the associated parameters fundamentally determine how fast and how far the flow travels. Entrainment occurs between the moving debris and the original channel bed, resulting in an increase in debris volume and a change in its composition (Revellino et al., 2004; Hungr and Evans, 2004; Mcdougall, 2006; Iverson 2012; Mcdougall, 2016; Shen et al., 2019). The flowing debris generates the basal shear stress to the bed materials and causes failure of these bed materials. Therefore, the flowing debris erode and entrain the failed and mobilized bed materials (McDougall, 2006). Because the entrainment plays a critical role in determining the volume of moving debris, it also influences the overall flow characteristics (Shen et al., 2019). In addition, it is essential to model debris flows over realistically complex three-dimensional (3D) terrain, in which the natural slope angle and channel width change in an unpredictable manner. For this reason, the interpretation of depositional characteristics of debris flows is fairly limited in two-dimensional (2D) channel models. By contrast, there are only limited studies available which fully consider the rheological characteristics, entrainment, and 3D morphology of the terrain over which debris flows occur.

This study investigated the effect of barrier locations on the characteristics of debris flows by using numerical modeling based on the smoothed particle hydrodynamics (SPH) method. We performed back-analysis of debris-flow events recorded in Korea using the DAN3D code (Dynamic Analysis of Landslides in Three Dimensions; McDougall and Hungr, 2004; McDougall and Hungr 2005; McDougall, 2006). For this study, three events were chosen, where (i) a debris-flow event occurred along a short path, (ii) a debris-flow event occurred along a long path, and (iii) an event with multiple source locations. The input parameters for modeling were obtained from field surveys. On 3D terrain models of the field sites, closed-type barriers were created along the debris-flow channel, while varying the source-to-barrier distance. Changes in velocity and volume were monitored as well as the volume of debris trapped and deposited by the barrier. In addition, overflow behavior was analyzed in relation to the barrier design and installation method. Our results contribute better understanding of debris-flow behaviors associated with barriers installed as a mitigation measure, and can be used to determine optimum designs and installation locations for debris-flow barriers.

Study methods

Numerical code and governing equations

This study used a smoothed particle hydrodynamics (SPH) base code called DAN3D (Dynamic Analysis of Landslides in Three Dimensions; McDougall and Hungr, 2004; McDougall and Hungr 2005; McDougall, 2006) to simulate the flow dynamics of viscous, liquid-like debris. This code adopts the concept of equivalent fluid, which homogenizes a multi-component moving mass to a single-phase and single-component fluid (Hungr, 1995). The characteristics of the debris flow are expressed as particles which can be split apart by topographic changes or obstacles, avoiding the problems of mesh distortion (McDougall, 2016). The main governing equations are expressed, as follows:

$$ \frac{\partial h}{\partial t}+h\left(\frac{\partial {v}_x}{\partial x}+\frac{\partial {v}_y}{\partial y}\right)=\frac{\partial b}{\partial t} $$
(1)
$$ \rho h\frac{\partial {v}_x}{\partial t}=\rho {hg}_x+{k}_x{\sigma}_{z_{\left(z=b\right)}}\left(-\frac{\partial h}{\partial x}\right)+{k}_{yx}{\sigma}_z\left(-\frac{\partial h}{\partial y}\right)+{\tau}_{zx_{\left(z=b\right)}}-\rho {v}_x\frac{\partial b}{\partial t} $$
(2)
$$ {\sigma}_{z_{\left(z=b\right)}}=\rho h\left(g\cos \alpha +\frac{v_x^2}{R}\right) $$
(3)
$$ {\tau}_{zx_{\left(z=b\right)}}=-\left({\sigma}_{z_{\left(z=b\right)}}f+\frac{\rho {gv}_x^2}{\xi}\right) $$
(4)
$$ \frac{\partial b}{\partial t}={E}_t={E}_s{hv}_x $$
(5)

where σz is the bed-normal stress, ρ is the material bulk density, h is the bed-normal flow depth, α is the inclination of the bed from the horizontal, R is the bed-normal radius of curvature of the path in the direction of motion, vx is the flow velocity in the direction of motion, b is the bed-normal erosion-entrainment depth, g is the acceleration due to gravity, k is the stress coefficients, τzx is the basal shear stress, f is the friction coefficient, ξ is the turbulence parameter, Et is the time-dependent erosion rate, and Es is the displacement-dependent erosion rate.

DAN3D was chosen for this study because of the following advantages. (a) The code can simulate debris flows occurring over complex 3D digital elevation model (DEM). However, it is worth noting that DAN3D employs a depth-averaged method for z-direction; therefore, it is not strictly 3D modeling but 2D modeling. Instead, the pillar-like particles flow over a 3D DEM as the particles cannot overlap each other, as shown in Eqs. (1–3). For a given number of particles for each simulation, the initial particle volume is determined by dividing the debris source volume by the number of particles. Therefore, the debris particles in this SPH model start with the uniform volume in the beginning. But soon, when these particles begin to flow, the particle volume increases with the entrainment according to its own path based on Eq. (5) while keeping the number of particles constant (McDougall and Hungr, 2004; McDougall, 2006). Difference in flow path of each particle results in the different entrained volume.

(b) The code employs an erosion rate to take the entrainment process into account. Herein, an exponential growth in the debris volume with respect to flow displacement is assumed, based on the effects of momentum transfer between the main flow and the bed materials, as shown in Eq. (5). The erosion rate Es is calculated from the initial debris-flow volume Vi, the final debris-flow volume Vf, and the path length S, as follows (McDougall and Hungr, 2004):

$$ {E}_s=\ln \left(\frac{V_f/{V}_i}{S}\right) $$
(6)

In this study, we measured these values (Vi, Vf, and S) from the field surveys to estimate the erosion rate Es. Therefore, the debris-flow volume growth due to the entrainment is determined by the the final volume Vf, the initial volume Vi, and the travel distance S. Meanwhile, the time-dependent erosion rate Et, which is equivalent to the rate of the bed-normal erosion depth increment ∂b/∂t, is proportional to the debris depth h and the flow velocity vx (see Eq. (5)). This implies that the reductions in the debris height and velocity caused by barrier installation would retard the entrainment.

(c) Various rheological models, including Newtonian, plastic, Bingham, frictional, or Voellmy models, are available. This allows a wide range of the flow characteristics associated with landslides and debris flows to be modeled. Herein, the Voellmy model was chosen to model the rheology of debris materials, as shown in Eq. (4). (d) DAN3D requires fairly affordable computational resources and time by using simplified governing equations and simplified boundary conditions.

In this numerical simulation, a barrier is created and considered a terrain of which the elevation rapidly changes, and no specific model for debris-flow-barrier interactions is implemented. The radius of curvature R in Eq. (3) becomes negative at the upstream face of the barrier. In addition, the term cos(α) in Eq. (3) approaches to zero due to the vertically installed barrier. Thereby, the bed-normal stress σz(z = b) in Eq. (3) can be even negative when gcosα + vx2/R < 0. Note that the bed-normal stress σz(z = b) acts against the gravity-driven term ρhgx in Eq. (2). Accordingly, the particle acceleration significantly increases with relatively small or negative bed-normal stress. As a result, the particles may show unrealistically rapid velocity and sometimes burst, becoming airborne. With such bursting or airborne particles, the simulation is no longer interpretable. Therefore, we flattened the terrain around the barrier as a deposition zone to reduce the flow velocity and minimize particle bursting even though the curvature R is still negative. Detailed description on the deposition zone follows (see the “Barrier installation in 3D morphological DEM” section). More details on the DAN3D code can be found in McDougall and Hungr (2004, 2005) and McDougall (2006).

Barrier installation in 3D morphological DEM

A closed-type barrier was designed and numerically modeled in the code. We varied the source-to-barrier distance from 0.3 to 0.9L to examine the effect of barrier location along a channel. The barrier was placed either at the uppermost upstream location close to a debris-flow initiation point (0.3L; case 0.3L), at the mid-upper upstream location (0.5L; case 0.5L), at the mid-lower downstream location (0.7L; case 0.7L), or at the lowermost downstream location close to a village or depositional area (0.9L; case 0.9L). The distance L is defined as the length of the actual flow path between the initiation point to the deposition region (e.g., flat area or village in Fig. 1a).

Fig. 1
figure 1

Schematics of a debris-flow barrier installed: a barrier location, b the top view, and c the side view

Full size image

A barrier was created by numerically adjusting the elevation values in the topographic map to achieve the desired shape and size, using ArcGIS software (ver. 10.5). The height and thickness of the barrier were set to be 5 m and 3 m, respectively. The barrier width was determined by considering the width of the debris flow simulated without barriers, and was set to be twice as wide as the width of the frontal part of the debris flow. The barrier was placed to be perpendicular to the flow direction.

Prior to the barrier placement, the ground surface was flattened to the lowest elevation of the channel, following the typical flattening earthwork practices. Without such flattening, debris flows readily overflow when a barrier is installed, because of the concave topography of the flow channels. The presence of this flattened area facilitates the entrapment and deposition of debris. The flattened area was 10 m long and its width was equal to the width of the barrier. Barrier shoulders (or wings) were often added to both the left and right sides of the barrier when the basin was too wide to block with a single barrier, as shown in Fig. 1c. The cell size of 1 m was used throughout this study. In practices, the ground around a barrier is typically treated to prevent erosion by debris flow. For this reason, surface erosion in the channel near the barrier was restricted in the model, over an area 63 m long and 100–120 m wide, with the barrier located at the center. Flow simulations were conducted with barriers installed at the different locations (i.e., cases 0.3L, 0.5L, 0.7L, and 0.9L), as well as one without a barrier (i.e., case REF); i.e., each debris-flow event was simulated with a total of five cases.

Study areas

Three sites with different characteristics were selected in this study: a short-path basin (referred to as Site 1), a long-path basin (Site 2), and a short-path basin with multiple sources (Site 3). These sites are the areas in which debris-flow events have occurred in the past. Numerical simulation of these debris-flow events using DAN3D requires the determination of input parameters, including the friction coefficient, the turbulence parameter, and the channel erosion rate. These parameters were back-calculated from the source volume, final debris volume, and depositional area. Table 2 summarizes the details of the debris-flow events and the input properties used for Sites 1, 2, and 3. In addition, the geotechnical engineering parameters were obtained from laboratory tests of soil samples from the study sites, as listed in Table 2. The internal friction angle and the unit weight were used to calculate the Rankine pressure coefficient and the mass of the debris flow, respectively. The number of particles was kept constant as 1,000. The initial particle size in our SPH simulations was determined by dividing the source volume by the number of particles; thus, the initial volume of particles was kept uniform. As a result, each particle had the initial volume of 0.164 m3 in Site 1, 0.057 m3 in Site 2, and 0.404 m3 in Site 3. Note that the entrainment increases the particle volume based on Eq. (5) while keeping the number of particles constant.

Table 2 Input properties of debris-flow events and barrier conditions in Site 1,2, and 3
Full size table

Table 2 also lists the geometry and location of the barriers installed in all cases at the studied sites. Closer to the downstream area, the barrier size became greater due to increases in debris volume and channel width. At each site, the baseline case, in which the debris flow occurred without a barrier (i.e., case REF), was conducted first, followed by the four cases with different barrier locations (cases 0.3L, 0.5L, 0.7L, and 0.9L).

Site 1: short-path basin

A short-path basin located in Bungnae-myeon, Yeoju, Korea (37° 24′ 14.80″ N and 127° 39′ 25.90″ E; Fig. 2) was chosen as Site 1. The elevations of flow initiation and deposition location were 197 m and 155 m above sea level, respectively, and the total length of the debris-flow path was estimated to be 262 m. Figure 2b shows clear evidence of the landslide occurrence, which is presumed to be the initiation location of the debris-flow event associated with the landslide (i.e., referred to as source location or initiation location). Figure 2c shows the basal channel eroded by the debris flow. The initial debris volume (or source volume) was estimated to be 164 m3 and the final debris volume to be 2,304 m3, based on the field survey. All debris was observed to have been deposited downstream close to a village (i.e., end point in Fig. 2a). The Site 1 debris-flow event was reproduced via numerical simulation using the following input parameters: a friction coefficient of 0.03, a turbulence parameter of 500 m/s2, and an erosion rate of 1.200×10-2%/m, as shown in Table 2. At Site 1, the barriers were placed downstream of the initiation zone in all cases (Fig. 2a).

Fig. 2
figure 2

Field condition of Site 1. a Locations of initiation, end point, and barriers; b a field survey evidence of a landslide event; c the erosion occurred along the basal channel

Full size image

Site 2: long-path basin

A long-path basin located in Toechon-myeon, Gwangju, Korea (37° 25′ 36.30″ N and 127° 21′ 48.20″ E; Fig. 3) was chosen as Site 2. The elevations of flow initiation and deposition locations were 318 m and 145 m above sea level, respectively, and the debris-flow path was 589 m long. The scar where the landslide occurred (Fig. 3b), the eroded basal channel in the upstream region (Fig. 3c), and the debris deposited (e.g., soils, boulders, and driftwood) in the runout zone (Figs. 3c and 3d)  were identified in the field survey. The field survey of Site 2 estimated an initial debris volume of 57 m3 and a final debris volume of 3,914 m3. Most of the debris was observed to have been deposited downstream, close to a roadway. The input parameters were determined based on the source volume and final debris volume, as follows: a friction coefficient of 0.03, a turbulence parameter of 750 m/s2, and an erosion rate of 7.800×10-3%/m (Table 2). At Site 2, the barriers were placed downstream of the initiation zone in all cases.

Fig. 3
figure 3

Field condition of Site 2. a Locations of initiation, end point, and barriers; b an evidence of a landslide at the initiation zone; c the erosion occurred along the basal channel; d the deposition of debris flow; and e the deposition near the village

Full size image

Site 3: short-path basin with multiple sources

Another short-path basin with multiple sources, located in Heungcheon-myeon, Yeoju, Korea (37° 20′ 59.80″ N and 127° 31′ 43.90″ E; Fig. 4a), was chosen as Site 3. The elevations of the highest initiation and the end point locations were 187 m and 100 m above sea level, respectively, and the debris-flow path was 197 m long. The debris-flow event involved five initiation zones (Fig. 4a), where small five landslides occurred and the resultant debris merged downstream as it flowed down the basin channel. The scar of a landslide in one of the initiation zones was observed during the post-event field survey, as shown in Figure 4b. The eroded basin channel and basin floor (Fig. 4c), a large volume of the debris (e.g., soil, boulders, and driftwood) deposited in the runout zone (Fig. 4d), and the final debris deposited in the flat area close to a village (Fig. 4e) were also identified. The debris flow was deposited over a particularly wide area (Fig. 4e); and optical satellite imagery from Google Earth allowed the quantitative examination of the runout area, as shown in Fig. 5.

Fig. 4
figure 4

Field condition of Site 3. a Locations of multiple initiations, end point, and barriers; b an evidence of landslide at the initiation zone; c an eroded channel; d the deposition of debris including soils, boulders, and driftwoods; and e the deposited fan of debris in the downstream

Full size image
Fig. 5
figure 5

Effect of Voellmy rheology parameters in Site 3. a The satellite image after the event. The yellow line indicates the deposited area of the debris flow, and b the simulated deposited area with respect to the rheology parameters. Note that the values in the matrix indicate the final volume

Full size image

The initial debris volume and the final volume were estimated to be 405 m3 and 8,064 m3, respectively, based on the field survey and satellite imagery of Site 3. The model input parameters were determined based on the source volume, final debris volume, and depositional area, as follows: a friction coefficient of 0.1, a turbulence parameter of 350 m/s2, and an erosion rate of 1.185×10-2%/m (Table 2). Figure 5b shows a visual matrix expressing the simulated deposited areas. The final volume and the deposited areas were matched with the satellite images (see Fig. 5). The width of the barrier was fixed at 40 m because it was difficult to select one reference debris flow, due to multiple debris-flow sources. Note that in cases 0.3L, 0.5L, and 0.7L, the barriers were placed in the upper part of the lowermost source, again owing to the multiple source locations.

Determination of rheological and entrainment parameters

The choice of rheological model has a pronounced effect on debris-flow characteristics (i.e., velocity and deposition). Among various rheological models available, Voellmy model was chosen in this study (Voellmy, 1955). The Voellmy model was originally proposed to model snow avalanches. Because of similarities in velocity and thickness between snow avalanches and debris flows, the model is widely used for debris-flow analysis (Körner, 1976; McDougall, 2016). The Voellmy model requires the friction coefficient f and turbulence parameter ξ as input parameters (see Eq. (4)). The friction coefficient governs the deposition characteristics and thus affects the extent of the deposition zone and runout distance. On the other hand, the turbulence parameter heavily affects the mobility of a debris flow and determines the range of the flow velocity (McDougall, 2016). Figure 6 and Table 3 compile the friction coefficients and turbulence parameters used in modeling of debris flows and rock avalanches in previous studies. It can be seen that debris-flow analysis typically uses a friction coefficient of 0–0.2, and a turbulence parameter between 0 and 1,000 m/s2. Accordingly, a friction coefficient less than 0.2 and a turbulence parameter less than 1,000 m/s2 were used as base guidelines, and the specific values for each site were back-calculated based on the recorded velocities, runout distances, and depositional areas of historical debris flows, as suggested by previous studies (e.g., Cepeda et al., 2010; Aaron et al., 2016; McDougall, 2016; Aaron et al., 2019).

Fig. 6
figure 6

Friction coefficient and turbulence parameter of the Voellmy model for debris flows and rock avalanches. The data were gathered from the literature in Table 3

Full size image
Table 3 Voellmy model parameters used for landslide modeling
Full size table

Aforementioned, we determined the friction coefficients of the studied sites by comparing the deposition areas between the simulation results and the satellite images while satisfying the final volume obtained from field surveys. The friction coefficients were 0.03 for Sites 1 and 2, and 0.1 for Site 3. The turbulence parameters were 500 m/s2 for Site 1, 750 m/s2 for Site 2, and 350 m/s2 for Site 3. These values are within the range of the previously reported values for debris flows (see Fig. 6). The turbulence parameter was arbitrarily determined from the range suggested in previous literature (Fig. 6) because no records on the debris-flow velocities in the studied sites are available. The debris-flow event at Site 2 also has a long runout distance of 589 m, and the Voellmy parameters are consistent with the previous work by Choi et al. (2019). In their study, they used a turbulence parameter of 800 m/s2 and a friction coefficient of 0.03 for a particular class of the debris flows in Korea, which has a markedly rapid flow velocity up to ~ 18 m/s and a long runout distance up to 627 m.

The displacement-dependent erosion rate Es in our study was mathematically determined by the path length, and initial and final debris volume of each event. The erosion rates Es in Sites 1 and 3 are relatively greater than those in Site 2, presumably owing to the relatively short path length and large volume growth. By contrast, Site 2 has the smaller Es value mainly due to the relatively long path length. For the debris-flow events that occurred in Korea, Vasu et al. (2018) used the erosion rate Es of 4.400×10-3%/m and Choi et al. (2019) used the erosion rate Es of 7.800×10-3%/m. Both values are close to or at least in the same order of magnitude with our values. Meanwhile, various physical factors affect the entrainment, including physical characteristics of bed soils, soil depth, and slope, as listed in Table 4. Proper characterizations of debris materials and bed soils are expected to provide physical meaning of the back-calibrated erosion rate.

Table 4 Review on factors affecting the entrainment
Full size table

Simulation results

Figures 7, 8, and 9 show the snapshots of the simulated debris flows with different barrier locations for Sites 1, 2, and 3, respectively. Figure 10 shows the changes in velocity and volume of moving debris flows with respect to barrier location. In the SPH method, the viscous fluid is modeled as particles, and thus the location, velocity, and volume of each particle can be tracked over time. In this work, the velocity values of all moving particles were averaged to determine the representative velocity of the flowing debris. The volume of flowing debris was also determined by summing the volume of all moving particles. In the cases with a barrier, the particles which flowed over the barrier were only considered when computing the velocity and volume. For Sites 1 and 2, these values were plotted with respect to the debris-flow distance (or runout distance) from the initiation location, and the runout distance was computed with respect to the front edge of the flowing debris.

Fig. 7
figure 7

Final results of the simulated debris flow for Site 1. a Before debris flow, b case REF, c case 0.3L, d case 0.5L, e case 0.7L, and f case 0.9L

Full size image
Fig. 8
figure 8

Final results of the simulated debris flow for Site 2. a Before debris flow, b case REF, c case 0.3L, d case 0.5L, e case 0.7L, and f case 0.9L

Full size image
Fig. 9
figure 9

Final results of the simulated debris flow for Site 3. a Before debris flow, b case REF, c case 0.3L, d case 0.5L, e case 0.7L, and f case 0.9L

Full size image
Fig. 10
figure 10

Effect of barrier locations on the flow velocity and volume of moving debris: a Site 1, b Site 2, and c Site 3. Note that the filled circles indicate the location of each barrier. In Site 3, the filled circles and empty circles indicate the time when debris flows approach each barrier and the village, respectively

Full size image

In addition, several volume-related characteristics of debris flows were obtained to examine the barrier effectiveness, which includes the entrained volume Ventrained and deposited volume Vdeposited as well as the barrier capacity Vbarrier. The entrained volume Ventrained is the debris volume that added to the debris flow, and it was determined at the moment when the debris were all trapped. If overflowed, the value was determined at the moment when the group of overflowed debris entered the village. The total debris-flow volume Vdebris is defined as the debris-flow volume passing through the corresponding barrier location, and it is identical to the summation of the entrained volume and initial source volume, i.e., Vdebris = Ventrained + Vinitial. The volume of deposited debris Vdeposited was obtained by summing the volume of particles that had zero velocity. Lastly, the barrier capacity Vbarrier (or the trap capacity of a barrier) was estimated by assuming the slope ratio (H:V) of 2:1 for the debris deposit close to the wall, i.e., Vbarrier = Wb × Hb2, where Wb is the barrier width and Hb is the barrier height. This 2:1 slope ratio is consistent with the experimental observations in previous studies—the debris are usually piled up near the barrier, forming a slope; thereby, the follow-up debris can overflow over the deposited slope (Ng et al., 2014 and 2015; Zhou et al., 2018; Song et al., 2019; Siyou et al., 2020; Tan et al., 2020). Accordingly, the retained volume ratio can be estimated, i.e., Rret = Vdeposited/Vdebris. And the debris-to-barrier volume ratio can be defined, i.e., Rd/b = Vdebris/Vbarrier. Figure 11 shows the debris volumetric characteristics in association with the barrier location.

Fig. 11
figure 11

a The volumetric characteristics of debris flows with respect to barrier location and capacity. b The retained volume ratio (Rret) and the debris-to-barrier volume ratio (Rd/b). The initial, entrained, and deposited volumes indicate the debris source volume, the final entrained volume which adds to the debris flow by erosion, and the final volume of debris deposited by each barrier, respectively. The deposited volume and entrained volume were estimated at the moment either when the debris were all trapped or when the debris entered the village if overflowed. The barrier capacity Vbarrier was defined by assuming the 2:1 slope ratio of debris deposition. The retained volume ratio Rret is defined as Vdeposited/Vdebris, and the debris-to-barrier volume ratio Rd/b is defined as Vdebris/Vbarrier

Full size image

Site 1: short-path basin

In case REF, the velocity increased to its maximum of ~ 7 m/s due to the steep slope angle in the uppermost part of the flow path, and, thereafter, gradually decreased as the debris flowed downslope (Fig. 10a). The debris volume gradually grew, reaching ~ 900 m3 as it approached the village. Upon barrier installation, sudden drops in velocity and volume were observed upon the debris collision with the barrier. Immediately after the collision with the barrier, the debris velocity and volume dropped to zero in all cases (0.3L, 0.5L, 0.7L, and 0.9L cases).

Owing to the erosion-induced entrainment process, the total debris volume gradually increased as the source-to-barrier distance increased (e.g., 270 m3 for case 0.3L, 490 m3 for case 0.5L, 700 m3 for case 0.7L, and 970 m3 for case 0.9L; Fig. 11a). Meanwhile, the retained volume ratios Rret were unity; the deposited volumes were the same with the total debris volumes in all cases. This indicates that most of the debris was trapped and deposited by the barrier, and no overflow occurred regardless of barrier location. It is because the barrier capacities were significantly greater than the total debris volume (Fig. 11a), in which the debris-to-barrier volume ratios Rd/b were less than 0.7 (Fig. 11b). This is also consistent with the debris depth against the barrier height. The debris depth in the vicinity of a barrier was kept less than 2 m, which was significantly lower than the wall height (5 m), as shown in Fig. 12. This depicts the debris deposition behavior close to a barrier when no overflow takes place.

Fig. 12
figure 12

Debris depths in the vicinity of a barrier: case 0.9L in Site 1 with no overflow (a–c) and case 0.9L in Site 2 with an overflow (d–f). The cross-sectional views in a and d are obtained along the central section of the barrier

Full size image

Site 2: long-path basin

In case REF, the debris-flow velocity increased to its maximum of ~ 13 m/s in the uppermost part of the flow path, and thereafter gradually decreased while flowing downslope (Fig. 10b). The debris volume grew exponentially due to the entrainment of material, and reached ~ 2,100 m3 when approaching the village. In cases 0.3L and 0.5L, where the barrier was placed at the uppermost and the mid-upper upstream locations, respectively, the velocity and volume were observed to drop and remain at zero after the barrier, which implies that all the debris was trapped by the barrier. By contrast, in cases 0.7L and 0.9L, where the barrier was installed at the mid-lower and the lowermost downstream locations, respectively, both the velocity and volume were seen to recover from zero after colliding with the barrier, though the recovery was minor (Fig. 10b). Approximately 220 m3 and 55 m3 of debris volume overflowed to the downstream in cases 0.7L and 0.9L, respectively, because the total debris volume exceeded the deposition capacity of the barrier. Thus, while the initial section of the debris flow was trapped, its tail overflowed and passed over the barrier.

This type of overflow is widely observed when the barrier is installed in downstream regions because the debris volume can exponentially grow by entrainment. The total debris volume exponentially increased due to the entrainment as the source-to-barrier distance increased (e.g., 200 m3 for case 0.3L, 480 m3 for case 0.5L, 1,080 m3 for case 0.7L, and 1,990 m3 for case 0.9L; Fig. 11a). Although the barrier capacity also increased toward the downstream (Fig. 11a and Table 2), it appeared to be insufficient. For cases 0.7L and 0.9L, the barrier capacities were comparable to or less than the total debris volume (Fig. 11a); therefore, the debris-to-barrier volume ratios Rd/b were greater than 0.8 (Fig. 11b). Therefore, the overflow occurred in those cases and the retained volume ratios Rret were less than 1. The debris depth close to the barrier wings was often higher than the wall height (5 m), which also corroborates the overflow occurrence, as shown in Fig. 12.

This result demonstrates that the installation of barriers nearer the initiation (or source) location is more effective as it interrupts flow before the debris grows substantially by entrainment, when the location of initiation can be predicted. Similar results have been reported in other studies (Remaître et al., 2008; Choi et al., 2019; Shen et al., 2019).

Site 3: short-path basin with multiple sources

In case REF, the debris-flow velocity increased to a maximum of ~ 5 m/s and decreased to ~ 3 m/s while flowing downslope (Fig. 10c). The debris volume grew exponentially due to the entrainment of material, and reached ~ 4,100 m3 when approaching the village. In this site, debris flows were initiated simultaneously at five different locations, as shown in Fig. 4; two flows were initiated upstream of 0.3L, one between 0.3L and 0.5L, one between 0.5L and 0.7L, and the last one between 0.7L and 0.9L. Therefore, part of the flow began and moved downslope without any collision with the barrier in cases 0.3L, 0.5L, and 0.7L. As a result, only minor reductions in velocity were caused by the interactions with the barrier. For the same reason, the debris volume remained consistently greater than 450 m3 despite partial entrapment and deposition by the barriers (Fig. 10c). In case 0.9L, however, the debris velocity and volume were significantly reduced (to zero) by the barrier. Most of debris was effectively trapped, although a minor portion of the debris overflowed. Unlike the other sites above, the entrained volume decreased as the source-to-barrier distance increased and it was the most in case 0.3L, as shown in Fig. 11a. This is a unique feature caused by the multiple source locations, which is not comparable to the single source cases. This highlights that identifying debris source locations is a prerequisite for barrier location determination.

Discussion

Effect of barrier location and capacity

Growth of the debris volume was generally observed due to entrainment of channel material. A longer distance from initiation to deposition led to a greater volume of debris flow. Closer to the initiation zone, the basin width is likely to be small, and the debris volume growth by entrainment is thus likely to be minimal. The size of the required barrier becomes larger toward the downstream with an increase in the source-to-barrier distance owing to the increased debris volume by entrainment.

Upon the installation of a properly sized barrier, our simulation results show that the scale of the debris flows is significantly reduced. For a debris-flow event with a short flow path (e.g., Site 1), the debris flows can be readily terminated by the barrier regardless of barrier location, owing to the sufficiently large barrier compared to the small debris volume (see Fig. 11). In Site 1, the debris-to-barrier volume ratios Rd/b were consistently less than 0.7; thus, no overflow occurred.

For a debris-flow event with a long flow path (e.g., Site 2), however, part of the debris cannot be captured by the barrier in cases 0.7L and 0.9L when the debris-to-barrier volume ratios Rd/b are greater than 0.8 (Fig. 11b). With overflows, the retained volume ratios Rret are less than 1. It is because the debris volume exponentially grows and also the debris depth increases with entrainment, of which the rate is faster than the rate of barrier size increase toward downstream in our study. In an additional simulation with the greater erosion rate (Site 2 with Es = 1.000×10-2%/m), the overflows occur more significantly as the debris-to-barrier volume ratios Rd/b increase further. Such overflowing debris can grow by entrainment while flowing a remaining path and can cause further damage in areas downstream of a barrier (Liu et al., 2013; Cuomo et al., 2019; Shen et al., 2019; Choi et al., 2019). Only in the case with the barrier installed at the uppermost location (e.g., case 0.3L), the debris-to-barrier volume ratio Rd/b was 0.3, and hence no overflow occurs with the retained volume ratio of unity.

These observations strongly suggest that it is most effective to install the barrier close to the location of landslide initiation before entrainment-induced debris volume growth. This result is consistent with previous studies by Remaître et al. (2008), Shen et al. (2019), and Choi et al. (2019). However, it is worth noting that barrier installation at upstream locations is often not financially viable because of problems associated with access road construction and maintenance difficulties. In addition, our results consistently imply that the threshold debris-to-barrier volume ratio Rd/b for effective debris trapping with no overflow ranges approximately 0.6–0.7.

Multiple landslides can occur simultaneously in a single watershed when heavy rainfall or earthquakes occur. In such a case with multiple sources, part of the debris flow may reach the downstream area without any energy dissipation when the barrier is located upstream of one of the sources, as shown in cases 0.3L, 0.5L, and 0.7L of Site 3. The retained volume ratio Rret increases as the installed downstream through Rret is still less than 0.7 in case 0.9L (Fig. 11b). Therefore, barrier location needs to be carefully determined. If the locations of multiple sources are predictable, either a series of barriers can be installed close to each source, or one large barrier can be placed at the lower part of the main channel where the multiple debris flows converge (e.g., case 0.9L of Site 3). In the latter, the size of the barrier should be sufficiently large enough to trap the cumulative debris volume.

The barrier installation close to the source is expected as the most effective when the accurate prediction of such landslide initiation zones is possible. Therefore, accurate prediction of locations of landslide initiation as well as source volumes and entrainment rates is a pre-requisite to determine the optimum location and size of barriers. However, we currently lack the means to accurately predict initiation locations of debris flows. When one or multiple initiation locations are difficult to predict, either multiple barriers along potential flow channels or an individual large downstream barrier can effectively mitigate against debris flows.

Effect of deposition zone size

Closed-type barriers are designed to trap the majority of the debris, and thus the volumetric capacity of the barrier (or deposition capacity) is an important factor to optimal barrier design. However, an overflow inevitably occurs when the debris-flow volume exceeds the deposition capacity of the barrier (Remaître et al., 2008; Liu et al., 2013; Jeong et al., 2015, Kwan et al., 2015; Dai et al., 2017; Chen et al., 2019; Cuomo et al., 2019; Lee et al., 2019; Shen et al., 2019). The deposition capacity is determined primarily by barrier size (i.e., width and height), and then by the size of the flattened deposition zone, particularly the length perpendicular to the barrier (Figs. 1b and 1c).

For case 0.9L in Site 3, in which overflow occurs, the effect of deposition zone length is further examined by varying the length of the deposition zone from 10 to 20 m, as shown in Fig. 13. It appears that enlarging the deposition zone size can effectively reduce the overflow of debris. Overflow still occurs with a deposition zone length of 15 m, although the overflow volume is less than the case with a deposition zone length of 10 m (Figs. 13b and 13d). A deposition zone length of 20 m, however, results in the complete deposition of all debris with no overflow (Figs. 13c and 13d). It is also observed that when the deposition zone length is 20 m, the maximum deposition depth stays less than 4 m. By contrast, with the smaller deposition zones, it often becomes higher than 4 m, close to the barrier height of 5 m (Fig. 13e). As a result, it is concluded that the size of the deposition zone as well as the barrier size (height, width) determines the deposition capacity of the barrier and affects the occurrence of overflow (Fig. 13f; Liu et al., 2013; Choi et al., 2019; Lee et al., 2019).

Fig. 13
figure 13

Effect of deposition zone area on overflow behavior for case 0.9L of Site 3. a The deposition zone length of 10 m, b 15 m, c 20 m, d the volume of moving debris, e the maximum depth of the deposited debris by the barrier, and f the deposited volume of debris by the barrier. The filled and empty circles indicate the time when debris flows approached each barrier and the village, respectively

Full size image

Implication and limitation of the presented approach

The numerical approach presented here has several noteworthy implications and limitations, as follows. (a) The depositional characteristics of the debris flow can be effectively captured by incorporating a complex 3D terrain morphology. In particular, 3D terrain morphology is required to be incorporated for the assessment of barrier effectiveness because of the complex interactions between debris, terrain, and barrier as the debris collides with, bounces back from, or slides to the side of the barrier (e.g., Iverson et al., 2016). (b) It is possible to create barriers with various geometry, by modifying the digital elevation map (DEM) of the 3D terrain as desired. The cell size of the 3D DEM needs to be carefully manipulated to accurately reflect the barrier geometry and the surrounding terrain morphology. We believe that the SPH method offers possibility for modeling of particles which constitute debris materials; on the other hand, care should be taken for determination of the cell size as it would heavily affect the interactions among the SPH particles, barriers, and terrain. Further research on SPH-base modeling of debris particles is warranted as the entrainment or rheology models become challenging. (c) Emphasis should be given to the growth of debris volume by entrainment processes, because entrainment has a significant influence on the characteristics of debris flows (Revellino et al., 2004; Hungr and Evans, 2004; Mcdougall, 2006; Iverson 2012; Mcdougall, 2016; Choi et al., 2019; Shen et al., 2019). If the entrainment is not considered, the performance of a debris-flow barrier can be easily overestimated. (d) The rheology of the debris flow governs its mobility and depositional characteristics. In this study, the Voellmy model, which uses the friction coefficient and turbulence parameter, was chosen and the velocity and depositional characteristics of historical debris-flow events were re-created by using the back-calibrated coefficients (Figs. 5 and 6).

Some limitations are as follows. The determination of input parameters, such as erosion rate and rheological coefficients, based on field surveys and theoretical/experimental works is a critically important but daunting task which must precede numerical modeling; otherwise, the models may produce significantly different results. Therefore, caution should be given because some of the results presented here can be site-dependent due to the input parameters and terrain morphology. While the methodology presented here can be readily extended to other sites of interest, further extensive simulation study using simplified topology and various rheological parameters is warranted to draw more generalized conclusions. In addition, comparison among the modeling methods and codes, such as SPH-base models (McDougall and Hungr, 2004; Dai et al., 2017; Coumo et al., 2019), FVM-base models (Christen et al., 2010), FEM-base models (Jung et al., 2015; Kwan et al., 2015; Lee et al., 2019), and FDM-base models (Liu et al., 2013; Shen et al., 2019), is presumed to be beneficial to identify the limitation and applicability of a variety of modeling methods.

In the approach presented here, the barrier is installed by modifying the terrain morphology (or 3D DEM), and the barrier is thus considered not as a structural element but as geo-terrain. Therefore, direct assessment of the structural stability of the barriers is not possible. Cuomo et al. (2019) previously performed the analysis of debris flows considering the barrier using SPH method, in which the authors stated that the debris-flow velocity and height could be overestimated as the flow energy loss caused by the impact with the barrier is underestimated. The same conditions apply to our results, and this gives a conservative design option to engineers. Instead, for structural stability analysis, other analytical or numerical methods can be coupled with SPH-based approach by extracting debris-flow characteristics (e.g., depth, velocity, and volume) and using them as inputs (Dai et al., 2017; Chen et al., 2019).

Conclusion

This study investigates the influence of the source-to-barrier distance on the characteristics of debris flows via a numerical approach using DAN3D. As the debris-flow events are simulated, the velocity and volume of the debris flow are drastically reduced at the moment when the debris collide with the closed-type barrier. Meanwhile, that erosion-induced entrainment has profound effects on debris-flow behaviors, barrier performance, and likelihood of overflows. The debris flow volume in a longer flow path can grow further with the entrainment, and hence, the likelihood of overflow also increases for a given size of barrier. Therefore, installation of the barrier closer to the initiation location is found to be more effective by preventing the debris volume from growing by entrainment. Furthermore, the barrier deposition capacity needs to be sufficiently greater than the expected debris volume so as to maintain the debris-to-barrier volume ratio less than ~ 0.6–0.7. The overflow can be minimized by increasing the deposition capacity of the barrier, possibly through increasing barrier height or depositional area size. By contrast, for a debris-flow event with a small initial volume or with a short path, a closed-type barrier effectively blocks the debris flow regardless of the barrier location. In a debris-flow event with multiple sources, identification of such multiple initiation points is a prerequisite prior to the barrier installation in order not to miss one of the debris sources. Therefore, in such a case, it is necessary to install either small barriers close to each initiation point or a huge barrier at a location where the multiple debris sources converge. Our results provide the baseline data on the effect of source-to-barrier distance on debris flow as well as insight into determination of closed-type barrier location in a channel along which the debris flow is predicted. However, it should be noted that some of the results may be site-dependent due to the input parameters and terrain morphology. The presented methodology can be readily applicable in sites of interest and extended to determine appropriate location and size of debris-flow barriers and sizes of deposition zones, and to evaluate the effectiveness of different types of barrier.