Estimation of the area of sediment deposition by debris flow using a physical-based modeling approach

Quaternary International 503 (2019) 59–69

Contents lists available at ScienceDirect

Quaternary International journal homepage: www.elsevier.com/locate/quaint

Estimation of the area of sediment deposition by debris flow using a physical-based modeling approach

T

Hyunuk Ana, Minseok Kimb,∗, Giha Leec, Yeonsu Kimd, Hyuntaek Lime a

Department of Agricultural and Rural Engineering, Chungnam National University, Daejeon 305-764, Republic of Korea Geologic Environment Division, Korea Institute of Geoscience and Mineral Resources, 124 Gwahak-ro, Yuseong-gu, Daejeon, 34132, Republic of Korea c Department of Constructional Disaster Prevention Engineering, Kyungpook National University, Sangju 742-711, Republic of Korea d K-water Convergence Institute, Korea Water Resources Corporation, 125 1689beon-gil, Yuseong-daero, Yuseong-gu, Daejeon 34045, Republic of Korea e Department of Regional Infrastructure Engineering, Kangwon National University, Chuncheon 24341, Republic of Korea b

ARTICLE INFO

ABSTRACT

Keywords: Adaptive mesh refinement Debris flow Finite volume model Numerical model Voellmy rheological model

Debris flow triggered by shallow landslides in hillslope catchments is a main geological phenomenon driving landscape changes, and represents an important natural hazard. Numerous studies have assessed sediment transport and deposition by debris flows in hillslope catchments. Thus, the objective of this study is a development of two-dimensional debris flow model to estimate sediment transport and deposition in hillslope catchments. To simulate debris flow, we implemented a vertically integrated shallow-water governing equation based on the Voellmy rheological model and a simple entrainment model. In addition, we applied a quadtree grid structure to support adaptive mesh refinement, where the mesh for the simulation was automatically generated as the debris flow proceeded. Finally, a well-balanced numerical scheme for wet–dry transition treatment was implemented and implicit discretization of the general source terms, including the rheological term, was included for numerical stability. The developed model was verified based on a debris flow triggered by the 2011 Mt. Umyeon landslides in the Republic of Korea. Sediment transport was successfully generated and the sediment deposition area generally matched the field survey data well. Overall, the simulated sediment volume was in good agreement with the survey results, with an error below 1%.

1. Introduction Debris flows are common phenomena in mountainous areas worldwide. When caused by shallow landslides, they often result in extensive property damage and the loss of human life due to their high mobility and impact energy (Iverson, 2000). Therefore, there have been many efforts to reduce the risk of debris flows, ranging from the implementation of countermeasures to the avoidance of hazardous areas based on hazard mapping. Among such efforts, numerical modeling of flows has an important role in risk mitigation for its ability to predict and simulate the movement of flows over irregular topographic terrains. Many numerical models of debris flow have been developed in the past two decades (O'Brien et al., 1993; Beguería et al., 2009; Blanc and Pastor, 2009; Mergili et al., 2012). Such models can be either one-dimensional, in which flow movement in one spatial dimension is a crosssection of a single predefined width (Imran et al., 2001; Shrestha et al.,

2008; D'Aniello et al., 2015; Paik, 2015), or two-dimensional, in which flow movement is in two spatial dimensions and can take into account irregular topography (Hsu et al., 2010; Hussin et al., 2012; Kim et al., 2013). Most models use uniform rectangular grids or triangular grids. However, uniform grid-based models are often inefficient, particularly when paired with high-resolution elevation data, which is becoming increasingly popular due to the rapid development of LiDAR [Laser Radar, Light Detection And Ranging] and unmanned aerial vehicle technology (Scheidl et al., 2008; Bull et al., 2010; Adams et al., 2016; Kim et al., 2016). In debris flow simulations, the flow path is determined by the elevation model; therefore, resolution of the elevation model can critically affect the results of simulations. Compared to rectangular grid-based models, triangular ones require greater effort for grid generation, because they can generate efficient grids for high-resolution elevation data. Elevation data should be properly interpolated when the resolution of the original input data differs from the grid resolution in the numerical model. However, interpolation is more

Corresponding author. Tel:+82 42 868 3253; fax: +82 42 868 3414. E-mail addresses: [email protected] (H. An), [email protected] (M. Kim), [email protected] (G. Lee), [email protected] (Y. Kim), [email protected] (H. Lim). ∗

https://doi.org/10.1016/j.quaint.2018.09.049 Received 30 May 2018; Received in revised form 14 September 2018; Accepted 30 September 2018 Available online 14 October 2018 1040-6182/ © 2018 Published by Elsevier Ltd.

Quaternary International 503 (2019) 59–69

H. An et al.

difficult in triangular grids than rectangular grids because spatial data are usually provided in a rectangular or uniform grid format. The application of rectangular grids using the adaptive mesh refinement (AMR) technique is an attractive alternative. Using this method, the grid resolution can be adjusted based on the digital elevation model (DEM) or flow features. The AMR technique can improve the efficiency of numerical simulations (Berger and Oliger, 1984; Berger and Colella, 1989; Borthwick et al., 2001) and has been widely adopted in numerical simulations in various engineering and scientific fields including in the modeling of tsunamis (Popinet, 2011; Chida et al., 2014; de la Asunción and Castro, 2017), floods (An and Yu, 2011; An et al., 2015; Huang et al., 2015), magnetohydrodynamics (Balsara, 2001; Jiang et al., 2012), subsurface flow (Molins et al., 2015), storm surge (Mandli and Dawson, 2014), and ocean dynamics (Santilli and Scotti, 2015). In this study, we developed a two-dimensional debris flow model (Deb2D) based on a rectangular grid with AMR to estimate the sediment transport area in a hillslope catchment. We used a vertically integrated shallow-water governing equation and implemented a finitevolume method for discretization. The choice of the rheological model is essential in debris flow modeling (Rickenmann et al., 2006). The most common rheologies used in dynamic models include (Hussin et al., 2012) frictional (or coulomb) resistance (Hungr and McDougall, 2009), frictional-turbulent Voellmy resistance (Voellmy, 1964), visco-plastic Bingham resistance (Coussot, 1997; Malet et al., 2004), and quadratic resistance (O'Brien et al., 1993). We implemented the Voellmy resistance model for its simplicity, as it only requires two parameters, which makes the model easy to use and calibrate. It is also one of the most widely used rheologies in debris flow simulations (Hürlimann and Graf, 2003; Rickenmann et al., 2006; Pirulli and Sorbino, 2008; Hungr and McDougall, 2009; Klaus et al., 2015), and has been found to be in good agreement with historical debris flow events. The entrainment process is an important factor determining the magnitude and intensity of debris flows (Hussin et al., 2012). For example, entrained materials may accumulate 10–50 times their initial soil water volume (Remaître et al., 2005). Several entrainment models for debris flow simulations have been proposed in the past two decades, most of which are process-based entrainment rate approaches (Crosta et al., 2003, 2009; D'Ambrosio et al., 2003; Medina et al., 2008) or defined entrainment rate approaches (Beguería et al., 2009; Hungr and McDougall, 2009; Pastor et al., 2009). In this study, we used the defined entrainment rate approach of Sovilla et al. (2006), which is used in the RAMMS commercial software package. We also equipped a graphical user interface (GUI) in the model to allow users to easily run debris flow simulations and interpret the results. The remainder of this paper is organized as follows. The governing equations of the debris flow system and finite volume discretization of the governing equation are presented in Section 2 and the adaptive mesh generation and GUI features of the Deb2D are also described in Section 2. A verification of the developed model based on a debris flow event triggered by the 2011 Mt. Umyeon landslides in Korea is presented in Section 3, and the conclusions are presented in Section 4.

Fig. 1. Conceptual diagram of (a) the quadtree discretization and (b) its logic structure.

directions, and source terms, respectively. The vectors can be written as:

g=

hv huv , hv 2 + gh2/2

s

Qe = Sgx Sgy

Sfx , Sfy

(2)

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, Qe is the entrainment (Qe > 0) or deposition rate (Qe < 0), Sgx and Sgy represent gravitational acceleration in the x and y directions, respectively, and Sfx and Sfy represent driving friction in the x and y directions, respectively. The driving friction is given by the Voellmy model as follows:

Sfx =

g (u2 + v 2) u µgh + , u

Sfy =

g (u2 + v 2) v µgh + , v

(3)

where u is the Coulomb friction coefficient and ξ is the turbulent friction coefficient. Qe with one layer is given as described by Sovilla et al. (2006) as follows:

Qe (x , y, t ) =

K (u2 + v 2 ) 0

if hs (x , y, 0) else

t 0

Qe (x , y, ) d > 0

, (4)

where hs (x, y, 0) is the initial height of the entrainment layer at position (x,y) and K is a dimensionless entrainment coefficient. Note that hs (x, y, 0) should be given in spatial form as a simulation input. 2.1.2. Spatial and temporal discretization The governing equation is discretized using the finite-volume method. Numerical models with shallow-water governing equations often suffer from an imbalance between the gradient of the water depth and bed slope in irregular topography. This imbalance may cause unphysical perturbations and instability of the simulation near shocks or wet–dry transitions. A number of well-balanced numerical schemes have been proposed (e.g., Audusse et al., 2004; Audusse and Bristeau, 2005; Caleffi et al., 2006; Kim et al., 2008; Caleffi and Valiani, 2009). This model implements the hydrostatic reconstruction technique proposed by Audusse et al. (2004), which has been successfully applied to quadtree adaptive grid-based shallow-water models. The discretized governing equation is written as:

2. Materials and methods 2.1. Numerical modeling 2.1.1. Governing equation The model employs two-dimensional depth-integrated shallowwater equations to simulate debris flow. The hyperbolic conservation form of the mass and momentum balance equation is written as follows:

g q f + + = s, t x y

hu f = hu2 + gh2/2 , huv

h q = hu , hv

Ai, j t

qni, j+ 1

qni, j + Fi, j = Si, j + Sci, j + Sfi, j

(5)

Where Ai, j is the area of the cell (i j ), Sci, j is an additional source term to balance the presence of a source, Sfi, j is the entrainment source term and driving friction term, and Fi, j is the numerical flux term, which is defined as:

(1)

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 60

Quaternary International 503 (2019) 59–69

H. An et al.

Fig. 2. Main window in the Deb2D interface.

Fi, j = li, j, r Fi + 1/2, j

li, j, l Fi

1/2, j

+ li, j, t Fi, j + 1/2

li, j, b Fi, j

(6)

1/2

where li, j, p is the length of cell face p in cell (i j ), and the subscripts p = {r , l , t , b} indicate that the respective value is constructed on the right, left, top, or bottom side of the cell, respectively. Please note that (i j ) indices are used on a uniform template and i ± 1/2 and j ± 1/2 indicate the cell face on the uniform grid to clearly describe the numerical scheme in a two-dimensional space. Each flux term is calculated as follows:

Fi + 1/2, j =

(qi + 1/2 , j , q i + 1/2+, j),

Fi

Fi, j + 1/2 =

(qi, j + 1/2 , qi, j + 1/2 +),

Fi, j

1/2, j

=

(q i

1/2

=

( q i, j

1/2 , j, 1/2

qi

1/2+, j),

, q i, j

1/2 +),

(7) where is the numerical flux solver. The interface conserved values are defined as:

qi + 1/2

,j

hi + 1/2 , j h = i + 1/2 , j ui, j, r , hi + 1/2 , j vi, j, r

q i + 1/2+, j

hi, j + 1/2 qi, j + 1/2 = hi, j + 1/2 ui, j, t , hi, j + 1/2 vi, j, t

Fig. 3. Flow chart of the Deb2D simulation process.

hi + 1/2+, j h = i + 1/2+, j ui + 1, j, l , hi + 1/2+, j vi + 1, j, l hi, j + 1/2+

qi, j + 1/2 + = hi, j + 1/2+ui, j + 1, b , hi, j + 1/2 +vi, j + 1, b (8)

with the following hydrostatic reconstruction: 61

Quaternary International 503 (2019) 59–69

H. An et al.

Fig. 4. Target catchment of the Raemian APT catchment debris flow triggered by the 2011 Mt. Umyeon landslides.

Fig. 5. (a) Digital elevation model of the target area, where the yellow point represents the debris flow initiation area and brown polygon represents the entrainment area, and (b) initial debris flow volume.

62

Quaternary International 503 (2019) 59–69

H. An et al.

Table 1 Observed and simulated debris flow results with different parameters. (m/s2)

µ

Entrainment volume (m3)

Inundated deptha (m)

Max. velocityb (m/s)

0.01

100 200 500 1000

40519 42051 44154 45275

5.1 7.5 10.3 13.6

7.7 10.3 15.5 19.6

0.05

100 200 500 1000

38554 40041 42012 43145

3.6 4.5 8.4 10.9

6.3 8.7 12.9 17.9

0.06

1000

42581

10.3

17

0.07

1000

41152

9.7

16.2

0.08

1000

41589

8.7

15.7

0.09

1000

40133

7.2

14.6

0.1

200 500 1000

36889 39217 40394

2.4 3.2 5.8

6.1 10 14.1

Field survey

42,500

10

28

a b

,j

= max(0, hi, j, r + z i, j, r

hi, j + 1/2 = max(0, hi, j, t + z i, j, t z i + 1/2, j = max(z i, j, r , z i + 1, j, l ),

Si + 1/2, j = ghi2+ 1/2

ghi2, j, r /2 ,

0 = ghi2 1/2+, j /2 0 0 0 = ghi2, j + 1/2 /2

ghi2, j

0 0 1/2 +/2

+ li, j, t Si, j + 1/2

li, j, b Si, j

Si

1/2, j

Si, j

1/2

1/2,

ghi2, j, l /2 ,

ghi2, j, t /2

ghi2, j, b/2

,

. (10)

The additional source terms to ensure the balance for the secondorder scheme are defined as:

Sci, j =

li, j, r Sci, j, r + li, j, l Sci, j, l

Sci, j, r = g (hi, j, r

0 + hi, j )(z i, j, r 0

li, j, t Sci, j, t + li, j, b Sci, j, b, z i, j )/2 ,

0 Sci, j, l = g (hi, j, l + hi, j )(z i, j, l z i, j )/2 , 0 0 0 Sci, j, t = , g (hi, j, t + hi, j )(z i, j, t z i, j )/2

z i + 1/2, j ),

Sci, j, b =

z i, j + 1/2),

hi, j + 1/2 + = max(0, hi, j + 1, b + z i, j + 1, b

1/2, j

0 , j /2 0

=

z i + 1/2, j ),

hi + 1/2+, j = max(0, hi + 1, j, l + z i + 1, j . l

li, j, l Si

Si, j + 1/2

Estimated at Reamian APT. Estimated where the debris surged to the road.

hi + 1/2

Si, j = li, j, r Si + 1/2, j

g (hi, j, b

0 0 + hi, j )(z i, j, b

z i, j )/2

. (11)

(q i + 1/2 , j , qi + 1/2+, j) are calculated with the The fluxes Harten–Lax–van Leer contact approximate Reimann solver in the developed model. MUSCL-type unsplit discretization is combined with a predictor-corrector time-stepping scheme to construct an accurate second-order Godunov-type scheme. The predicator step computes a temporary cell-center value over the half-time step t/2 :

z i, j + 1/2), z i, j + 1/2 = max(z i, j, t , z i, j + 1, b). (9)

where ui, j, p , vi, j, p , hi, j, p and z i, j, p are the reconstructed values for the second-order accuracy at the cell interface p of cell (i , j ). The source terms are defined as:

Ai, j qni, j+ 1/2 = Ai, j qni, j

t n (Fi, j 2

Sin, j

Scin, j ),

(12)

and the corrector step updates the conservative values over the time

Fig. 6. Simulated debris flow area with (a) different turbulent friction coefficients and (b) different Coulomb friction coefficients. 63

Quaternary International 503 (2019) 59–69

H. An et al.

Fig. 7. Simulated flow depth over time, where u = 0.06 and ξ = 1000 m/s2.

step

t as:

series expansion as:

Ai, j qni, j+ 1 = Ai, j qni, j

t (Fin, j+ 1/2

Sin, j+ 1/2

Scin, j+ 1/2).

(13)

The treatments of additional source/sink terms from the entrainment and driving friction force are described in the following subsection.

qn t

where,

D( t) = I

t

n

Sf q

(19)

where I is a unit matrix and

(14)

Sf = q

Qe h

0

0

Sfx

0

(hu)

0

0 0

.

Sf y

(20)

(hv )

Equation (18) ignores the effects of off-diagonal entries. This method can be successfully applied to shallow water systems with Manning's friction (An et al., 2014). Denoting Sf n ( t ) = D 1 ( t ) Sf n yields the implicit scheme written as:

(15)

= Sf n + 1 .

(18)

qn+ 1 = qn + t D 1 ( t ) Sf n ,

Ai, j qni, j+ 1/2 = Ai, j qni, j

The Euler implicit discretization of the above equation yields:

qn + 1

(17)

where q = Ignoring high-order terms and substituting Equation (12) into Equation (11) yields:

The inclusion of a friction term can cause numerical instability within a discretized domain, particularly when a dry–wet transition exists. In other words, if the flow depth approaches zero, driving friction can result in an exaggerated force, which reverses the flow in a non-physical manner. Therefore, the friction term should be properly treated. The developed model implements an implicit discretization based on the splitting method (Liang and Marche, 2009), which is equivalent to solving the following equation:

dq = Sf . dt

( q) + o ( q2),

qn .

qn+ 1

2.1.3. Treatment of entrainment source and driving friction terms The additional source/sink terms from the entrainment and driving friction force are written as:

Qe Sf = Sfx , Sf y

n

Sf q

Sf n + 1 = Sf n +

t (Fin, j 2

qni, j+ 1/2 = qni, j+ 1/2 + Ai, j qni, j+ 1 = Ai, j qni, j

(16)

qin, j+ 1 = qni, j+ 1 + tSf

The right side of Equation (11) can be represented by the Taylor 64

t Sf n 2

Sin, j

( ),

t (Fin, j+ 1/2 n + 1/2

Scin, j ),

( t ).

t 2

Sin, j+ 1/2

Scin, j+ 1/2), (21)

Quaternary International 503 (2019) 59–69

H. An et al.

Fig. 8. Simulated debris flow velocity vectors over time, where u = 0.06 and 1000 m/s2.

of the flow domain. An ESRI *.asc raster format file is used for the DEM input. Buildings are optionally represented by ESRI polygon shape files by raising the bathymetry within the polygon. The simulation period and time step for saving are set at the bottom of the main windows. Note that the actual time step for computation is automatically determined by the conditions as follows:

Although the friction terms are discretized implicitly, the friction force in Equation (21) can cause exaggerated reverse flow, which makes no physical sense. Therefore, if the driving friction results in reverse flow, the debris flow stops. 2.1.4. Mesh generation and graphical user interface The adaptive mesh in the developed model is generated based on a quadtree grid structure. Fig. 1 shows a conceptual diagram of the quadtree structure and its logical tree structure. In this structure, a “parent cell” can have four “child cells,” where each child cell can be a parent cell of four children cells (Fig. 1a). The base of the tree is the “root cell,” and any cell without child cells is a “leaf cell”.” The “level” of a cell is defined by starting from zero for the root cell and adding one when children are added (Fig. 1b). The entire quadtree grid can be conveniently described based on the level and parent–child relationship (An and Yu, 2014). In Fig. 1, the refinement level of the smallest square is 4. The general rectangular domain is allowed in the developed model with a continuous connection of squares, which is the shape of the cell in the quadtree grid. Mesh refinement is automatically conducted until the mesh level reaches a predefined maximum level, where the difference in flow depth exceeds the predefined tolerance. The maximum and minimum grid levels can be adjusted by users. Cell merging is not conducted in the current version of the developed model. The current version of Deb2D runs on a Windows operating system. Fig. 2 shows the main windows of the program. Users can create, open, save, or close a project using the icons in the left toolbar. The flow domain can be set using five parameters, Nx, Ny, L, x0, and y0, which represent the number of root cells along the x- and y-directions, length of root cell, and location of x and y at the left-bottom edge

t < C min min i, j

x ij uij +

ghij

, min i, j

yij vij +

ghij

,

(22)

where the Courant number C is set as 0.5 in Deb2D. The simulation results (i.e., flow depth, maximum flow depth, flow velocity, and entrainment depth) are plotted at the center of the main windows as the simulation progresses. The results at different times can be plotted in the result tab at the top of the main window. The mesh level and refinement tolerance can be controlled at the bottom of the plotted graph tab. Fig. 3 shows a flow chart of the simulation process in Deb2D. 2.2. Study area Mt. Umyeon, with a peak elevation of 293 m, is located in the southern part of Seoul, the capital city of Korea (Fig. 4). The mountain is composed of highly weathered banded gneiss with subordinate granitic gneiss and has hill slopes with an average angle of 15°. There are many manmade, ecofriendly facilities on the mountain, such as a natural ecological park and mineral springs, which local residents can easily access. Before the 2011 Mt. Umyeon landslides, a few landslides were caused by torrential rainfalls produced by Typhoon Kompasu in 2010. Some of these induced debris flows in the northern valleys of the mountain, but the damage was less serious than that caused by the 2011 65

Quaternary International 503 (2019) 59–69

H. An et al.

Fig. 9. Adaptive mesh refinement over time with a maximum refinement level of 10.

raising the elevation at ground level by 30 m. The initial debris flow volume was estimated by comparing the LiDAR data before and after the event (Fig. 5b). The total entrainment volume was estimated to be 42,500 m3 based on a field survey, and the velocity of the flow when the debris surged to the road was estimated to be 28 m/s based on video analyses (Seoul City, 2014). The Raemian APT was reportedly inundated up to the third floor of buildings, which was assumed to be 10 m above ground level in this study.

Table 2 Simulation results with different mesh refinement levels. Level of refinement

10 9 8 7 a b

x (m)

0.88 1.76 3.52 7.03

Entrainment Volume (m3)

Inundated deptha (m)

Max. velocityb (m/s)

CPU Time (sec)

CPU Time in uniform grid (sec)

42581 42271 41467 36249

10.3 8.4 7.6 3.5

17 15.9 13.2 9.4

1258 188 31 10

– 1935 294 60

3. Results and discussion

Estimated at Reamian APT. Estimated where the debris surged to the road.

3.1. Application of the model

landslides. Despite the reinforcement of structural countermeasures against both landslides and debris flows, the same areas damaged by the landslides in 2010 were damaged again by the severe landslides in 2011. The Mt. Umyeon Landslides Research Report identified a number of damaged sites, with 140 slope failures and 33 debris flows (Seoul City, 2014).

Setting rheological parameter values is challenging because they are difficult to directly estimate from field surveys; however, they have significant effects on the results. Therefore, rheological parameters are usually calibrated by running back-analyses to fit the results with the observations. This method is often used to predict potential events (Pirulli and Sorbino, 2008). We tested 15 combinations of rheological parameters within a range of u = 0.01–0.1 and ξ = 100–1000 m/s2. The calibrated results are listed in Table 1. The height of the entrainment layer, hs (x , y ), was set uniformly as 2 m, considering the relationship between the entrainment volume and area. Fig. 6 shows the flow area based on different parameter combinations. Simulated flow areas were similar in the mountain area but differed in the residential area. Unfortunately, the observations of the inundated area in the residential area were less precise than those in the mountain. Based on a comparison of the values listed in Table 1, the optimal parameter set was selected as u = 0.06 and ξ = 1000 m/s2. The simulated entrainment volume and inundated depth at Raemian APT agreed well with

2.3. Parameter construction A target debris flow event occurred in the Raemian APT catchment (Fig. 4). It started at two locations near the top of the mountain, and split into two downstream (Fig. 5a). The flow damaged the residential area of Raemian APT, causing three casualties and damaging commuter vehicles on one of the main roads in the city. A number of field surveys and observations were conducted in the target catchment as a result of the damage. A DEM with 1 × 1 m resolution captured by LiDAR was used in this study. Buildings in the Raemian APT are represented by 66

Quaternary International 503 (2019) 59–69

H. An et al.

Fig. 10. Maximum depth of the debris flow with different mesh refinement levels.

the field survey; however, the maximum velocity where the debris surged over the road significantly differed from the reported value. It should be noted that the maximum velocity was not in good agreement with the reported value, regardless of the rheological parameter set. This was presumed to be due to limitations stemming from the assumption of a uniform entrainment depth, limited range of parameters, or simple entrainment process employed (Eq. (4)). These factors are important in debris flow modeling; however, there are still many uncertainties associated with these parameters and processes. Uncertainty analyses of the modeling process were beyond the scope of this study and will be addressed in future work. Fig. 7 describes the simulated flow over time using the calibrated parameters. The debris flow reached the residential area 50 s after starting the simulation and deposition continued until 90 s. The volume of the flow increased as it proceeded due to entrainment at the bottom. Fig. 8 shows the simulated flow velocity. The velocity increased as the flow proceeded, until reaching the residential area, at which point it was blocked by buildings and was deposited. Fig. 9 describes the adaptively refined mesh over time. The maximum and minimum refinement levels were predefined as 10 (Δx; 0.88 m) and (2Δx; 225 m), respectively, and the tolerance of the refinement was set as 0.01 m. The mesh was successfully refined along the flow path. The effects of the mesh refinement level were tested, and the results are presented in Table 2. The simulated results significantly differed depending on the refinement level. In particular, a refinement level of 7 (Δx; 7 m) yielded much less accurate simulation results compared to higher refinement levels. Fig. 10 shows the maximum depth of the debris flow based on different refinement levels. Overall, the simulation results were affected by refinement. Moreover, the computation times of the simulations

dramatically differed depending on refinement. For comparison, the computation times of the simulations using uniform grids are listed in Table 2. Overall, the use of AMR required 6–10 times less CPU time than uniform grids. Although the flow domain of the simulation was not optimal in this study, the results confirm the efficiency of the AMR technique. 3.2. Further study Debris flow initiation is usually triggered by landslides with subsurface water flow. Several spatial and temporal landslide prediction methods have been created including empirical, statistical, and physical approaches (An et al., 2016a). Among them, physical models can predict the location of landslides, which can be converted into debris flow initiations. Combining the developed model with slope stability prediction models, such as that of An et al. (2016b) or Baum et al. (2008), will be an integral part of our future work. 4. Conclusions We developed a finite-volume two-dimensional debris flow model, Deb2D, to analyze the sediment transport area on hillslope catchments. This model used a quadtree grid structure to support AMR, where the mesh for the simulation was automatically generated as the debris flow proceeds, and mesh refinement could be controlled by selecting the mesh refinement level in the developed model. A vertically integrated shallow-water governing equation based on the Voellmy rheological model and a simple entrainment model, with reference to previous research, were used to simulate the debris flow. We validated the 67

Quaternary International 503 (2019) 59–69

H. An et al.

developed model via analyses of the debris flow triggered by the 2011 Mt. Umyeon landslides in the Republic of Korea. The mesh was successfully adaptively generated along the area inundated by the debris flow over time. Furthermore, the simulation results were generally in good agreement with the field survey data. Based on a test of the sensitivity of the mesh refinement level, the size of the mesh significantly affected the accuracy of the simulation. In addition, the computation times of models based on adaptively refined mesh and uniform mesh were compared, which revealed that the adaptively refined mesh required 6–10 times less CPU time than the uniform mesh. These results support the use of the developed model in the analyses of physical surface processes in many fields, particularly those related to geomorphic changes and natural hazards.

Caleffi, V., Valiani, A., Bernini, A., 2006. Fourth-order balanced source term treatment in central WENO schemes for shallow-water equations. J. Comput. Phys. 218 (1), 228–234. Chida, Y., Mori, N., Yasuda, T., Mase, H., 2014. Optimum scheme of adaptive mesh refinement for tsunami simulation. Coastal Eng. Proc. 1, 11. https://doi.org/10.9753/ icce.v34.currents.11. Coussot, P., 1997. Mudflow Rheology and Dynamics, 1 edition. CRC Press, Rotterdam. Crosta, G.B., Imposimato, S., Roddeman, D.G., 2003. Numerical modelling of large landslides stability and runout. Nat. Hazards Earth Syst. Sci. 3, 523–538. https://doi. org/10.5194/nhess-3–523–2003. Crosta, G.B., Imposimato, S., Roddeman, D., 2009. Numerical modelling of entrainment/ deposition in rock and debris-avalanches. Eng. Geol. Mech. Velocity of Large Landslides 109, 135–145. https://doi.org/10.1016/j.enggeo.2008.10.004. de la Asunción, M., Castro, M.J., 2017. Simulation of tsunamis generated by landslides using adaptive mesh refinement on GPU. J. Comput. Phys. 345, 91–110. https://doi. org/10.1016/j.jcp.2017.05.016. D'Ambrosio, D., Di Gregorio, S., Iovine, G., 2003. Simulating debris flows through a hexagonal cellular automata model: SCIDDICA S3–hex. Nat. Hazards Earth Syst. Sci. 3, 545–559. https://doi.org/10.5194/nhess-3–545–2003. D'Aniello, A., Cozzolino, L., Cimorelli, L., Morte, R.D., Pianese, D., 2015. A numerical model for the simulation of debris flow triggering, propagation and arrest. Nat. Hazards 75, 1403–1433. https://doi.org/10.1007/s11069–014–1389–8. Hsu, S.M., Chiou, L.B., Lin, G.F., Chao, C.H., Wen, H.Y., Ku, C.Y., 2010. Applications of simulation technique on debris-flow hazard zone delineation: a case study in Hualien County, Taiwan. Nat. Hazards Earth Syst. Sci. 10, 535–545. https://doi.org/10.5194/ nhess-10–535–2010. Huang, W., Cao, Z., Pender, G., Liu, Q., Carling, P., 2015. Coupled flood and sediment transport modelling with adaptive mesh refinement. Sci. China Technol. Sci. 58, 1425–1438. https://doi.org/10.1007/s11431–015–5880–6. Hungr, O., McDougall, S., 2009. Two numerical models for landslide dynamic analysis. Comput. Geosci. 35, 978–992. https://doi.org/10.1016/j.cageo.2007.12.003. Hürlimann, M.D.R., Graf, C., 2003. Field and monitoring data of debris-flow events in the Swiss Alps. Can. Geotech. J. 40, 161–175. https://doi.org/10.1139/t02–087. Hussin, H.Y., Quan Luna, B., van Westen, C.J., Christen, M., Malet, J.-P., van Asch, T.W.J., 2012. Parameterization of a numerical 2-D debris flow model with entrainment: a case study of the Faucon catchment, Southern French Alps. Nat. Hazards Earth Syst. Sci. 12, 3075–3090. https://doi.org/10.5194/nhess-12–3075–2012. Imran, J., Parker, G., Locat, J., Lee, H., 2001. 1D numerical model of muddy subaqueous and subaerial debris flows. J. Hydraul. Eng. 127, 10. 2001)127:11959. https://doi. org/10.1061/(ASCE)0733–9429. Iverson, R.M., 2000. Landslide triggering by rain infiltration. Water Resour. Res. 36 (7), 1897–1910. https://doi.org/10.1029/2000WR900090. Jiang, R.L., Fang, C., Chen, P.F., 2012. A new MHD code with adaptive mesh refinement and parallelization for astrophysics. Comput. Phys. Commun. 183, 1617–1633. https://doi.org/10.1016/j.cpc.2012.02.030. Kim, D.H., Cho, Y.S., Kim, H.J., 2008. Well-balanced scheme between flux and source terms for computation of shallow-water equations over irregular bathymetry. J. Eng. Mech. 134 (4), 277–290. Kim, S., Paik, J., Kim, K., 2013. Run-out modeling of debris flows in Mt. Umyeon using FLO-2D. J. Korean Soc. Civ. Eng. 33. https://doi.org/10.12652/Ksce.2013.33.3.965. Kim, Y., Jun, K., Jun, B., 2016. Analysis of debris flow disaster area using terrestrial LiDAR and Kanako-2D. J. Korean Soc. Hazard Mitig. 16, 113–118. https://doi.org/ 10.9798/KOSHAM.2016.16.5.113. Klaus, S., Barbara, T., Brian, M., Christoph, G., Oldrich, H., Roland, K., 2015. Modeling debris-flow runout pattern on a forested alpine fan with different dynamic simulation models. In: Engineering Geology for Society and Territory–volume 2. Springer, Cham, pp. 1673–1676. https://doi.org/10.1007/978–3-319–09057–3_297. Liang, Q., Marche, F., 2009. Numerical resolution of well-balanced shallow-water equations with complex source terms. Adv. Water Resour. 32, 873–884. Malet, J.-P., Remaître, A., Olivier, M., 2004. Runout modelling and extension of the threatened area associated with muddy debris flows/Modélisation de l’écoulement et extension de la zone exposée à des laves torrentielles boueuses. Geomorphologierelief Processus Environnement–GEOMORPHOLOGIE 10, 195–209. https://doi.org/ 10.3406/morfo.2004.1218. Mandli, K.T., Dawson, C.N., 2014. Adaptive mesh refinement for storm surge. Ocean Model. 75, 36–50. https://doi.org/10.1016/j.ocemod.2014.01.002. Medina, V., Hürlimann, M., Bateman, A., 2008. Application of FLAT Model, a 2D finite volume code, to debris flows in the northeastern part of the Iberian Peninsula. Landslides 5, 127–142. https://doi.org/10.1007/s10346–007–0102–3. Mergili, M., Fellin, W., Moreiras, S.M., Stötter, J., 2012. Simulation of debris flows in the Central Andes based on open source GIS: possibilities, limitations, and parameter sensitivity. Nat. Hazards 61, 1051–1081. https://doi.org/10.1007/ s11069–011–9965–7. Molins, S., Day, M., Trebotich, D., Graves, D.T., 2015. Adaptive mesh refinement in reactive transport modeling of subsurface environments. AGU Fall Meeting Abstracts 21, H21A–H1327. O'Brien, S.J., Julien, Y.P., Fullerton, T.W., 1993. Two-dimensional water flood and mudflow simulation. J. Hydraul. Eng-ASCE 119. https://doi.org/10.1061/(ASCE) 0733–9429(1993)119:2244. Paik, J., 2015. A high resolution finite volume model for 1D debris flow. J. Hydroenviron. Res. 9, 145–155. https://doi.org/10.1016/j.jher.2014.03.001. Pastor, M., Haddad, B., Sorbino, G., Cuomo, S., Drempetic, V., 2009. A depth-integrated, coupled SPH model for flow-like landslides and related phenomena. Int. J. Numer. Anal. Methods GeoMech. 33, 143–172. https://doi.org/10.1002/nag.705. Pirulli, M., Sorbino, G., 2008. Assessing potential debris flow runout: a comparison of two simulation models. Nat. Hazards Earth Syst. Sci. 8, 961–971. https://doi.org/10.

Acknowledgments The developed software supports a simple GUI on Windows OS and is available to the public at the following website: https://sites.google. com/site/cnuaehelab/software/deb2d. This work was supported by a National Research Foundation of Korea (NRF) grant funded by the Korean government (MSIT) (No. NRF-2017R1C1B2005750) and a Basic Research Project of the Korea Institute of Geoscience and Mineral Resources funded by the Ministry of Science and ICT of Korea. References Adams, M., Fromm, R., Lechner, V., 2016. High-resolution debris flow volume mapping with unmanned aerial systems (UAS) and photogrammetric techniques. ISPRS J. Photogrammetry Remote Sens. 749–755. http://dx.doi.org/10.5194/isprsarchivesXLI-B1-749-2016. An, H., Yu, S., 2011. Numerical simulation of urban flash flood experiments using adaptive mesh refinement and cut cell method. J. Korea Water Resour. Assoc. 44. https://doi.org/10.3741/JKWRA.2011.44.7.511. An, H., Yu, S., 2014. Finite volume integrated surface-subsurface flow modeling on nonorthogonal grids. Water Resour. Res. 50 (3), 2312–2328. https://doi.org/10. 1002/2013WR013828. An, H., Yu, S., Lee, G., Kim, Y., 2015. Analysis of an open source quadtree grid shallow water flow solver for flood simulation. Quat. Int. 384, 118–128. https://doi.org/10. 1016/j.quaint.2015.01.032. An, H., Kim, M., Lee, G., The Viet, T., 2016a. Survey of spatial and temporal landslide prediction methods and techniques. Korean J. Agri. Sci. 43, 507–521. https://doi. org/10.7744/kjoas.20160053. An, H., Viet, T.T., Lee, G., Kim, Y., Kim, M., Noh, S., Noh, J., 2016b. Development of timevariant landslide-prediction software considering three-dimensional subsurface unsaturated flow. Environ. Model. Software 85, 172–183. https://doi.org/10.1016/j. envsoft.2016.08.009. Audusse, E., Bristeau, M.-O., 2005. A well-balanced positivity preserving “second-order” scheme for shallow water flows on unstructured meshes. J. Comput. Phys. 206 (1), 311–333. Audusse, E., Bouchut, F., Bristeau, M.-O., Klein, R., Perthame, B., 2004. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput. 25 (6), 2050–2065. Balsara, D.S., 2001. Divergence-free adaptive mesh refinement for magneto hydrodynamics. J. Comput. Phys. 174, 614–648. https://doi.org/10.1006/jcph.2001.6917. Baum, R.L., Savage, W.Z., Godt, J., 2008. TRIGRS: a Fortran program for transient rainfall infiltration and grid-based regional slope-stability analysis. US Geol. Survey Open File Rep. 2008–11592. Beguería, S., Van Asch, T.W.J., Malet, J.-P., Gröndahl, S., 2009. A GIS-based numerical model for simulating the kinematics of mud and debris flows over complex terrain. Nat. Hazards Earth Syst. Sci. 9, 1897–1909. https://doi.org/10.5194/nhess9–1897–2009. Berger, M.J., Colella, P., 1989. Local adaptive mesh refinement for shock hydrodynamics. J. Comput. Phys. 82, 64–84. https://doi.org/10.1016/0021–999189)90035–1. Berger, M.J., Oliger, J., 1984. Adaptive mesh refinement for hyperbolic partial differential equations. J. Comput. Phys. 53, 484–512. https://doi.org/10.1016/ 0021–999184)90073–1. Blanc, T., Pastor, M., 2009. Numerical simulation of debris flows with the 2D–SPH depth integrated model. In: Presented at the EGU General Assembly Conference Abstracts, vol. 1978. Borthwick, A., Cruz Leon, S., Józsa, J., 2001. Adaptive quadtree model of shallow-flow hydrodynamics. J. Hydraul. Res. 39, 413–424. https://doi.org/10.1080/ 00221680109499845. Bull, J.M., Miller, H., Gravley, D.M., Costello, D., Hikuroa, D.C.H., Dix, J.K., 2010. Assessing debris flows using LIDAR differencing: 18 May 2005 Matata event, New Zealand. Geomorphology 124, 75–84. https://doi.org/10.1016/j.geomorph.2010.08. 011. Caleffi, V., Valiani, A., 2009. Well-balanced bottom discontinuities treatment for highorder shallow-water equations WENO scheme. J. Eng. Mech. 135 (7), 684–696.

68

Quaternary International 503 (2019) 59–69

H. An et al. 5194/nhess-8–961–2008. Popinet, S., 2011. Quadtree-adaptive tsunami modelling. Ocean Dynam. 61, 1261–1285. https://doi.org/10.1007/s10236–011–0438-z. Remaître, A., Malet, J.-P., Olivier, M., 2005. Morphology and sedimentology of a complex debris flow in a clay–shale basin. Earth Surf. Process. Landforms 30, 339–348. https://doi.org/10.1002/esp.1161. Rickenmann, D., Laigle, D., McArdell, B.W., Hübl, J., 2006. Comparison of 2D debris-flow simulation models with field events. Comput. Geosci. 10, 241–264. https://doi.org/ 10.1007/s10596–005–9021–3. Santilli, E., Scotti, A., 2015. The stratified ocean model with adaptive refinement (SOMAR). J. Comput. Phys. 291, 60–81. https://doi.org/10.1016/j.jcp.2015.03.008. Scheidl, C., Rickenmann, D., Chiari, M., 2008. The use of airborne LiDAR data for the

analysis of debris flow events in Switzerland. Nat. Hazards Earth Syst. Sci. 8, 1113–1127. https://doi.org/10.5194/nhess-8–1113–2008. Seoul City, 2014. Research Contract Report: Addition and Complement Causes Survey of Mt, vol. 2011 Woomyeon landslide (In Korean). Shrestha, B.B., Nakagawa, H., Kawaike, K., Baba, Y., 2008. Numerical Simulation on Debris-flow Deposition and Erosion Processes Upstream of a Check Dam with Experimental Verification. Annuals of the Disaster Prevention Research Institute, Kyoto University, No. 51 B, pp. 613–624. Sovilla, B., Burlando, P., Bartelt, P., 2006. Field experiments and numerical modeling of mass entrainment in snow avalanches. J. Geophys. Res. 111, F03007. https://doi. org/10.1029/2005JF000391. Voellmy, A., 1964. On the Destructive Force of Avalanches. Alta Avalanche Study Center.

69