TiVASS Software Article

Environmental Modelling & Software 85 (2016) 172e 172 e183

Contents lists available at  ScienceDirect

Environmental Modelling & Software j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m/ m/ l o c a t e / e n v s o f t

Development of time-variant landslide-prediction software considering three-dimensional subsurface unsaturated �ow Hyunuk An  a , Tran ran The Viet  b , Giha Lee b, Yeonsu Kim  c , Minseok Kim  c , Seongjin Noh d,  Jaekyoung Noh a  * ,

a

Dept. of Agricultural  &  Rural Eng., Chungnam National University, Daejeon, South Korea Dept. of Construction  &  Disaster Prevention Eng., Kyungpook National University, Sangju, South Korea c International Water Resources Research Institute, Chungnam National University, Daejeon, South Korea d Department of Civil Engineering, University of Texas at Arlington, Arlington, USA

b

a r t i c l e

i n f o

 Article history: Received 12 May 2016 Received in revised form 9 August 2016 Accepted 15 August 2016

Keywords: Landslide Slope stability 3D subsurface unsaturated  � ow Numerical model Mt. Umyeon Umyeon

a b s t r a c t

An accurate landslide-susceptibility landslide-susceptibility assessment is fundamental for preventing landslides and minimizing damage. damage. In this study, a new time-variant time-variant slope-stabili slope-stability ty (TiVaSS) (TiVaSS) model for landslide prediction prediction is developed. A three-dimensional (3D) subsurface  � ow model is coupled with the in �nite slope-stability model to consider the effect of horizontal water movement in the subsurface. A 3D Richards' equation is solved numerically for the subsurface  � ow. To overcome the massive computational requirements of  the 3D subsurface  � ow module, partially implicit temporal discretization and the simpli �cation of  �  � rstorder order spatial spatial discretizati discretization on are proposed and applied applied in TiVaSS. TiVaSS. A graphical graphical user interface interface and twodimensional data visualization are supported in TiVaSS. The model is applied to a 2011 Mt. Umyeon landslide in the Republic of Korea, and its overall performance is satisfactory. ©  2016 Elsevier Ltd. All rights reserved.

Software/data availability  Name of software software TiVaSS TiVaSS Developer's email address  [email protected] þ82-42-821-5797 Tel þ82-42-821-5791 Fax Year  � rst available available It will be open in public within within this year Software required required windows OS, windows 7 or later Program language Cþþ Program size 16.3 MB zip  � le Download link https://sites.go https://sites.google.com/site/c ogle.com/site/cnuaehelab/ nuaehelab/ software/tivass/TiVaSS_release.zip? ¼0&d¼1  User manual  &  tutorial attredirects¼ attredirects is now in preparation.

*

Corresponding author. author. E-mail address:  address:   [email protected]  (J. Noh).

http://dx.doi.org/10.1016/j.envsoft.2016.08.009 1364-8152/©  2016 Elsevier Ltd. All rights reserved.

1. Introduction Landslides are common natural hazards in mountainous areas, causing signi�cant damage and casualties worldwide every year. A high rainfall intensity or frequency can trigger landslides. An accurate landslide prediction is fundamental for preventing such disasters and minimizing casualties and property damage (Liao ( Liao et al., 2011; Borga et al., 2002; PapathomaeKohle et al., 2015; Pradhan and Lee, 2010; Ciabatta et al., 2016). 2016). Shallow landslide modeling falls into three categories: empirical (Sirangelo ( Sirangelo et al., 2003; Aleotti, 2004; Guzzetti et al., 2007, 2008; Bezak et al., 2016), 2016), statistical statistical (Dickson and Perry, 2016; Guzzetti et al., 1999, 2005, 2006; 2006;  Pham et al., 2016; Soeters and van Westen, 1996), 1996), and physically based models mod els.. Among Among these, these, physi physical cally ly based based mod models els are prefer preferred red because of their accuracy and ability to forecast the spatial and temporal occurrence of landslides (Raia ( Raia et al., 2014). 2014). They describe the triggering processes of shallow landslides and provide spatially variable slope-stability information, often as a safety factor (FS  (FS ), ), which is an index expressing the ratio between the local resisting and driving force. The in�nite slope-st slope-stability ability model and simple simple in�ltration-runoff process are commonly combined for computing the FS  the FS  in   in physically based slope-stability models. Those models are €

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

widely used and packaged as computer programs such as SINMAP (Stability Index MAPping), SHALSTAB (Shallow Landsliding Stability Model), TRIGRS (Transient Rainfall In�ltration and Grid-Based Regional Slope stability), and GEOtop-FS. SINMAP (Pack et al., 1999) and SHALSTAB (Dietrich and Montgomery, 1998; Montgomery and Dietrich, 1994) consider a simple steady-state hydrological process under constant rainfall. This means that both models are useful for producing a spatially distributed slope stability but are limited to the temporal prediction of the slope stability because of the steady-state description of  hydrological  � uxes. On the other hand, TRIGRS (Baum et al., 2008; Alvioli and Baum, 2016) computes the transient pore-water pressure using the solution of the one-dimensional (1D) Richards equation with the assumption of a simple exponential soil-water retention relationship (Gardner, 1958) and can produce a timevariant distribution under time-variant rainfall conditions. Thus, this model can forecast the timing and distribution of shallow landslides (Baum et al., 2010) and has been widely applied in the past decade (Baum et al., 2010; Liao et al., 2011; Liu and Wu, 2008; Luan et al., 2010; Montrasio et al., 2011; Saadatkhah et al., 2014; Peres and Cancelliere, 2016). In recent years, the predictive accuracy and computational ef �ciency of TRIGRS have been improved via the addition of a probabilistic approach (Raia et al., 2014) and the parallelization of the code (Alvioli and Baum, 2016). GEOtop-FS (Simoni et al., 2008) solves the three-dimensional (3D) subsurface �ow numerically on the basis of the 3D Richards equation and considers the physical hydrological process for estimating the slope stability. GEOtop is also used for predicting soil erosion and deposition (Zi et al., 2016). In the context of landslide early-warning systems, practical physically based slope-stability models are essential for determining the in�uence of rainfall recharge on the subsurface hydrological behavior and soil mechanics in triggering landslides. That is, advanced slope-stability tools should account for the time-varying subsurface water  � ow during a landslide event and consequently provide accurate spatio-temporal information about the slope failure (Apip et al., 2010). This study aims to develop a new time-variant slope stability (TiVaSS) model, which combines a 3D subsurface  � ow model and an in�nite slope stability model. An analytical solution is not available for the 3D subsurface  � ow, and a numerical technique is used to solve the problem. A numerical solver establishes the general soil-water retention relationships between the water content and the pore-water pressure such as van Genuchten model (1980) or Brooks-Corey model (1964). 3D subsurface  � ow usually requires a large amount of computational resources because fully implicit temporal discretization is essential for the numerical stability owing to the high nonlinearity of the Richards equation. The large computational requirement can hinder the practical application of landslide-occurrence assessment with the 3D subsurface �ow. To reduce this requirement, the 3D subsurface  � ow is solved in a partially implicit manner in TiVaSS, with the assumption that the horizontal  � ow is far slower than the vertical  � ow. The vertical �ow is solved implicitly, and the horizontal �ow is solved explicitly with a low-order discretization. Further, TiVaSS supports a graphical user interface (GUI). The developed model is applied to the 2011 Mt. Umyeon landslides in Seoul, Korea for the evaluation of its performance and functions. This paper is organized as follows. The in�nite slope-stability model and 3D subsurface �ow model, including the numerical scheme used in TiVaSS, are described in the next section. Details about the Mt. Umyeon landslide, including the data set, are provided in Section  3. The application results and discussion are also presented in Section 3. We draw our conclusions in Section 4.

173

2. Model description  2.1. In �nite slope-stability model Fig. 1 shows a schematic of an in �nite slope-stability model. Assuming an in�nite slope and failure parallel to the slope surface yields the following shear and normal stress:

t ¼

s

 ¼

T   W  ¼ cos f  sin  f ; b=cos  f b

(1)

P   W  ¼  cos2 f; b=cos f b

(2)

where W  ¼  g sDb is the soil weight (kg/m), T   is the shear force (kg/ m), P  is the normal force (kg/m), D is the soil depth (m), b  is the slope width (m), 4 is the slope angle (rad), and  gs is the unit weight of the soil (kg/m3). According to Mohr eCoulomb theory, the shear strength of an in�nite slope is

S  ¼ c  þ s0 tan 4;

(3) 0

where  c  is the cohesion (kg/m2),  s is the effective stress (kg/m2), 0 which was expressed as  s ¼  su by Terzaghi (1943), u is the pore water pressure (kg/m2), and  4  is the angle of internal friction (rad). Then, the FS   is de�ned as

 S  FS  ¼ : t

(4)

Substituting Eqs. (1)e(3)   into Eq. (4)   yields the following equation:

 c  þ gs D cos2 f  u tan 4 : FS  ¼ gs D sin  f  cos  f





(5)

Substituting u  ¼  g wj  into Eq. (5) gives

 c  þ Dgs  cos 2 f  gw j tan  4 ; FS  ¼ gs D sin  f  cos  f





(6)

where  gw is the unit weight of water (kg/m3), and j  is the pressure head of the subsurface water (m). Eq. (6) is rewritten by separating the time-variant term and steady terms as (Iverson, 2000)

jgw  tan  4 tan  4 c  ; þ þ FS  ¼  gs D sin  f  cos  f tan f gs D sin  f  cos  f

(7)

Eq. (7)  is valid for saturated soil. However, Terzaghi's effective stress is unsatisfactory for unsaturated soil ( Lu and Likos, 2006; Lu

Fig. 1.   Schematic of the in �nite slope-stability model.

174

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

et al., 2013). If the soil is highly dried, the large negative pressure head causes an excessive suction force with a Terzaghi's effective 0 stress of  s ¼ su.   Lu and Likos (2006)   generalized the effective stress as follows: 0

s

¼  s   ua  þ  s s ;

(8)

where ua  is the air pressure, and  s s is the suction stress, which is a characteristic function of the saturation or matric suction and is expressed in a closed form as follows: s

(9)

¼ ðua   uw ÞS e :

s

Here, uw is the water pressure, and  S e is the effective saturation, which is de�ned as

S e  ¼

q  qr  ; qs   qr 

 

(10)

where q   is the volumetric moisture content, qs   is the saturated moisture content, and qr  is the residual moisture content. Assuming 0 ua ¼ 0 and u w  ¼  jg w  yields an effective stress of  s ¼  s    j S egw. Therefore, Eq. (7) is generalized for variably saturated soil as

jS e gw  tan 4 tan  4 c  : FS  ¼  þ þ gs D sin  f  cos  f tan  f gs D sin  f  cos  f

 

(11)

If the soil is saturated, the effective saturation becomes zero, and Eq. (11) is equivalent to Eq. (7). The only time-variant value of Eq. (11) is the pressure head because the effective saturation is usually given as the function of the pressure head.  2.2. 3D subsurface �ow model The 3D subsurface �ow system of partially saturated soil is formulated by the Richards equation as follows: vqðjÞ vt 

 

  V $ðK ðjÞVðj þ z ÞÞ  q ¼ 0 ;

vq

dV   vt 

V 

Z 

n$K Vðj þ z ÞdvV  

vV 

Z 

qdV  ¼ 0 ;

 

vqi; j;k vt 



X  p

V 

n K  p Vðj þ z Þ p S  p   V i; j;k qi; j;k  ¼  0; ,

Vðj

þ z Þr  ¼

Vðj

þ z Þ f  ¼

Vðj

þ z Þt  ¼

ðj þ z Þiþ1; j;k   ðj þ z Þi;k;k D x

ðj þ z Þi; jþ1;k   ðj þ z Þi; j;k D y

; ;

 

(14)

where the subscripts  p  ¼  {r ,l, f ,b,t ,bot } denote the right, left, front, back, top, and bottom surfaces of the control volume (i, j,k), respectively,  S  p is the surface area of the control volume, and n is a normal vector along the control-volume surface. The subscripts  p ¼  {r ,l, f ,b,t ,bot } can also be expressed as ( i þ  1/2, j,k), (i   1/2, j,k), (i, j þ 1/2,k), (i, j  1/2,k), (i, j,k þ 1/2), and (i, j,k  1/2), respectively. A

(15)

ðj þ z Þi; j;kþ1  ðj þ z Þi; j;k : D z 

The gradients for the left, back, and bottom surfaces are estimated in the same manner. Note that D x(¼D y) is far larger than  D z  and that the vertical gradient is usually greater than the horizontal gradient. Eq. (14) is usually discretized in a fully implicit manner as

V i; j;k

n qnþ1 i; j;k  qi; j;k

Dt 



X  p

¼ 0 ; n$K  pnþ1 Vðj þ z Þ pnþ1 S  p   V i; j;k qnþ1 i; j;k (16)

(13)

where  v V   is the control-volume boundary.  q  and  j  are assumed to be the cell-averaged values according to the �nite-volume approximation. Thus, Eq. (13) can be rewritten as

V i; j;k

second-order estimation of the gradient at the surface can be obtained by the coordinate transformation method ( Jie et al., 2004; Ruhaak et al., 2008; An and Yu, 2012). Nevertheless, to reduce the computational cost, a simple �rst-order approximation is used in TiVaSS, as follows:

(12)

where j   is the pressure head (m); q  is the volumetric moisture content; 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 the general source term (1/s), including rainfall. Eq.   (12)  is discretized within a �nite-volume framework. The mesh used in TiVaSS is orthogonal horizontally and non-orthogonal vertically, as shown in Fig. 2. Integrating Eq. (12) over the control volume V and applying the GausseGreen divergence theorem gives (An and Yu, 2014)

Z 

Fig. 2.  Mesh used in TiVaSS: orthogonal horizontally and non-orthogonal vertically.

where the superscript n  is the time index. This scheme requires considerable computational resources because the large linear matrix should be solved iteratively, and matrix solvers have a computational cost of the order of  N m, where N  is the number of  cells, and m is a value between 1 and 2 that depends on the linear solver. To reduce the computational requirement, partially implicit temporal discretization is implemented in TiVaSS, as follows:

V i; j;k 

n qnþ1 i; j;k  qi; j;k

Dt 

X



X

n$K  pn Vðj þ z Þ pn S  p

 p¼r ;l; f ;b

n$K  pnþ1 Vðj þ z Þ pnþ1 S  p   V i; j;k qnþ1 i; j;k

 p¼t ;bot 

¼ 0 :

 

(17)

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

Eq. (17) results in multiple tridiagonal matrices along the x- and  y-directions, and each matrix system can be solved in separate manner. Therefore, the resulting system consumes computational resources on the order of  N . In this model, van Genuchten's [1980] soil-water retention curve and Mualem's [1976] unsaturated hydraulic conductivity function are used, as follows:

q  qr  ¼ S e  ¼ qs   qr  1=2 K  ¼  K s S e

(

11=n

)  

1 1 þ ajjjn

 

;

v

  1

v

n =ðn 1Þ 11=n 1  S e v

2

  v

v

(18)

 

;

(19)

where a  and nv are van Genuchten parameters that depend on the soil properties, and K s  is the saturated hydraulic conductivity (m/s). The maximum rainfall intensity, r max, is estimated using Darcy's law as

r max  ¼  K s



ðj þ z Þsur   ðj þ z Þtop jtop :  ¼ K s 1 þ D z =2 D z =2

 



 

(20)

Here, the subscript sur   represents the surface of the ground, where it is assumed that the pore-water pressure is zero, and the subscript top represents the cell immediately below the surface. If  the rainfall intensity is larger than r max, the over�ow is eliminated from the  � ow domain. The discretized system of Eq.  (17) is highly nonlinear, necessitating iterative calculations. For this, we employ the Newton iteration, which is expressed as



 J jnþ1;m



nþ1;m Dj





¼ F jnþ1;m ;

 

(21)

where m  denotes the iteration level,  j  represents the vector of  j ,  Jð¼  v F=vjÞ   is the Jacobian matrix, and nþ1;m Dj

   

for all cells (

nþ1;m Dj

nþ1;m Dj

  is less than the user-speci �ed tolerance

   d j , where dj  denotes the convergence

tolerance). In this study, a tolerance of  dj ¼ 104 m is used. Because the system of Eq. (17) is implicit only for one dimension, the Jacobian matrix is a tridiagonal matrix and can be solved by the Thomas algorithm.  2.3. Other features of TiVaSS  The combination of Eqs.   (11) and (17)  gives a spatiotemporal slope-stability distribution for time-variant rainfall. TiVaSS is coded in Cþþ and supports a GUI using the Qt framework. An ESRI arc format �le or a constant value is recognized as spatially distributed input data such as a digital elevation model (DEM), the soil depth, friction angle, soil cohesion, root cohesion, saturated soil weight, and saturated hydraulic conductivity without a zoning process. Two-dimensional (2D) visualization for simulation results, including the FS , the pore-water pressure, and the input data, is supported. The 3D subsurface �ow module and the spatially uniform 1D in�ltration module can be optionally selected in TiVaSS. The 1D in�ltration module is based on the 1D version of Eq. (12), as follows

vt 



The

discretization are applied for Eq.  (22). The 1D in�ltration model is only valid for uniform soil properties and the initial pore-water pressure without a horizontal �ow. If soil properties and the initial pore-water pressure are uniform and the horizontal �owcan be ignored, the pore-water pressure becomes uniform throughout the domain. This means that we are not required to compute the subsurface �ow for all domain cells; the pore-water pressure computed in one cell can be applied to all the study areas. TiVaSS does not account for surface �ow routing and excess surface water, which is removed from the �ow domain because the downstream is hardly ever unsaturated while the upstream is saturated. Only in this situation can the surface water affect the subsurface  � ow and the FS . The current version of TiVaSS runs on a Windows operating system. Fig. 3 shows the main windows of the program. Users can create, open, save, or close a project using icons in the left toolbar. Data are input on the left side of the main windows. An ESRI *.asc raster format  � le is used for the DEM input. A constant value or an ESRI *.asc raster format  � le can be used for the soil-depth setting. After the input of the DEM and soil depth, a simulation domain is created by setting the number of vertical soil layers. The numbers of  cells along the x- and y-directions are automatically set according to the raster �le. Unless all the input raster �les have the same domain, TiVaSS gives an error message. The other slope-stability parameters, friction angle, soil cohesion, saturated soil weight, and hydraulic conductivity can be set using two options: a constant value or an ESRI *.asc raster format �le. The soil-water retention curve is uniform throughout the domain. The van Genuchten soilwater retention model is the only soil-water retention model available in the current version of TiVaSS. Input data dincluding distributed 2D data, rainfall, and soil-water retention curvesdcan be plotted at the top of the main windows. The initial and end times are set at the bottom of the main windows. The simulation results, i.e., the FS , 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 windows.

¼  j nþ1;mþ1  jnþ1;m . The iteration process of Eq.   (21)

continues until

vqðjÞ

175

v



v z 

K ðjÞ

vðj



þ z Þ  ¼ 0 : v z 

conventional

�nite

volume

  and

(22) backward

Euler

3. Model application  3.1. Study area and landslide event  As shown in   Fig. 4, Mt. Umyeon, with a height of 293 m, is located in the southern part of Seoul, which is the capital city of  Korea. The mountain consists of highly weathered banded gneiss with subordinate granitic gneiss and has hill slopeswith an average angle of 15 . There are many manmade, eco-friendly facilities within the mountain, such as a natural ecological park and mineral springs that local residents can easily access. Before the Mt. Umyeon landslides i n 2011, torrential rainfal l due to Typhoon Kompasu in 2010 caused a few landslides, and some of  these induced debris �ows in the northern valleys of the mountain, but the damage was less serious than that caused by the 2011 landslides. Despite the reinforcement of structural countermeasures against both landslides and debris �ows, the same areas damaged by the landslides in 2010 were damaged again by the severe landslides in 2011. The number of damaged sites was reported as 140 slope failures and 33 debris  � ows by the Mt. Umyeon landslides Research Report of Seoul City (2014). Two automatic weather stations (Namhyeon and Seocho) are operated by the Korea Meteorological Administration near Mt. Umyeon. Fig. 5 shows the rainfall intensity from July 26th at 00:00 to the 28th at 24:00 at the Nahyeon and Seocho stations. According to Yune et al. (2013), the cumulative rainfall for two months before the landslide event was 1498.5 mm at Namhyeon (1105.0 mm at Seocho), which is greater than the annual average rainfall in Seoul

176

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

Fig. 3.  Main windows of TiVaSS.

(1450.5 mm). Antecedent rainfall led to soil saturation on Mt. Umyeon, and several sites had already become vulnerable to slope failure. Field observation by  Son et al. (2012)  evidenced that the groundwater level at locations where debris �ows were initiated almost reached the surface. Additional torrential rainfall of 306 mm for the 24-h period from 9 a.m. on July 26th to 9 a.m. on the 27th produced surface runoff continuously and weakened the ground, which is made up of gneiss, eventually resulting in the tremendous landslides. As shown in  Fig. 6, most areas of the mountain collapsed. In particular, the residential areas near the mountain were damaged by the debris �ow: 16 people were killed, and basements were inundated by mud  � ows. The points (except for two green points, which represent the air force base at Mt. Umyeon) and yellow lines indicate 140 landslide-initiation coordinates and 33 debris-�ow paths, respectively. In this study, we generate time-variant slope stability maps for the three days from July 26th to 28th using TiVaSS and then validate the model accuracy with 140 landslide observations.

 3.2. Input data A 10 m  10 m resolution DEM of Mt. Umyeon is generated from a 1:5,000 digital topography map of the National Geographic Information Institute in Korea.   Tarolli and Tarboton (2006)  considered a 10-m DEM optimal for shallow landslide modeling, as it overcomes the local slope noises detected with a higher-resolution topography. Tarolli (2014) also suggested that a coarse DEM (10 m) has a resolution more suitable for analyzing major hydrogeomorphic processes than a �ne DEM (less than 10 m). The slope angle of each individual cell within the study domain was calculated using the same DEM data. The geologic parameters of  Eqs.   (7) and (11), which are related to the soil properties, were determined according to Table 1 following Park et al. (2013), who investigated the same landslide event using TRIGRS. All available data were obtained from the �eld investigation for landslide hazard restoration work conducted by the National Forestry Cooperative Federation, the Korean Society of Civil Engineers, and the Korean Geotechnical Society. After the landslides occurred on 27 July, a

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

177

Fig. 4. Topography of Mt. Umyeon and slope failure locations (blue dots). (For interpretation of the references to colour in this �gure legend, the reader is referred to the web version of this article.)

Genuchten parameters a and nv   were set as 20 (m 1) and 1.4, respectively, as shown in Fig. 7 following the retention curve in Park et al. (2013). Rainfall data is generated by Thiessen polygon method using data from Nahyeon and Seocho gauge stations. A uniform soil depth of 2 m was considered according to site investigation reports (Korean Society of Civil Engineering, 2011). The landslides mostly occurred between depths of 1 and 3 m. The setting of the initial condition is dif �cult because there is usually no available soil-water observation to refer to, while the initial pore-water condition signi�cantly affects the results of the slope-stability model. Raia et al. (2014) found that the modeling of  the stability condition in low-gradient terrain is very sensitive to the initial conditions, which are uncertain and dif �cult to determine spatially. Therefore, the initial conditions should be set carefully. TiVaSS provides two options for the initial conditions: spatially constant conditions and distributed pore-water conditions. The latter can be input by an ESRI raster format �le (*.asc). Here, constant initial conditions (jinit  ¼ 0.3 m) are calibrated manually. Fig. 5.  Rainfall intensity and cumulative rainfall from 26th 00:00 to 28th 24:00 at the Namhyeon and Seocho weather stations.

total of 58 geotechnical investigation boreholes were drilled for collecting soil and hydrologic and geological information. Data from 13 boreholes and 19 soil samples were used in our analysis. With the large database, the average values were used. The van

 3.3. Results The results of TiVaSS are shown in Fig. 8. The distribution of the pore-water pressure presented in  Fig. 8 is obtained for the bottom subsurface layer, and the pore-water pressure used in Eq.  (11) is that at the soil bottom. Most of areas were in a stable condition

178

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

Fig. 6.   Landslide and debris-�ow inventories.

 Table 1 Parameter setting. K s  (m/s)

gs  (kg/m3)

4  (

1.3E-05

1874

29.63

)

c s  (kg/m2)

c r  (kg/m2)

qs  (m 3/m3)

qr  (m 3/m3)

1037

0

0.5

0.18

Fig. 7.  Soil-water characteristic curve for the study area (Park et al., 2013).

until the heavy rainfall on the 27th from 7 to 10 a.m., and the porewater from previous rainfall barely reached the bottom soil layer at 7 a.m. on the 27th. After 7 a.m. on the 27th, the pore-water

pressures in the bottom soil layer changed drastically, causing a rapid change in the FS . At 9 a.m. on the 27th, the bottom soil layer became saturated, and the unstable areas began to increase. At 10 a.m. on the 27th, almost all the areas became saturated, and the minimum FS   values were estimated. For comparison, TiVaSS with the 1D-in�ltration option was performed using the same parameters setting and data. The results are shown in Fig. 9. The pore-water pressure is not presented here, because the 1D-in�ltration option assumes the spatially uniform distribution of the pore water. The results of the 1D in �ltration are very similar to those of the 3D subsurface �ow. This is attributed to the heavy rainfall during a relatively short time, which caused saturation throughout the study area. To precisely validate the model performance, we employed slope-stability classes (unstable, FS  <   1: stabilizing factors are needed for stability; quasi-stable, 1 < FS < 1.25: minor destabilizing factors can lead to instability; moderately stable, 1.25 < FS  <   1.5: moderate destabilizing factors lead to instability; stable, 1.5 <  FS : only major destabilizing factors lead to instability) at the 140 observed landslide initiations shown in Fig. 4, as listed in Table 2. Before the heavy rainfall, 84 locations(points) were estimated as stable, and 48 locations were estimated as moderately stable. After the heavy rainfall, 95 locations were estimated as unstable and 18 locations were estimated as quasi-stable among the 140 locations. We found that 15 locations were unconditionally stable regardless of the pore-water pressure. Most of the areas were saturated and the maximum pore-water pressure, i.e., the minimum FS value, was obtained at 10 a.m. on the 27th, as shown in Fig. 8. However, the landslides were observed by the �eld survey in unconditionally

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

Fig. 8.  Safety-factor (FS ) maps and pore-water pressure for the bottom subsurface layer obtained by TiVaSS on the 27th from 7 to 10 a.m.

179

180

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

Fig. 9.  Safety-factor (FS ) maps obtained by TiVaSS with the 1D-in �ltration option on the 27th from 7 to 10 a.m.

 Table 2 Slope-stability classes computed by TiVaSS on the 27th from 7 to 10 a.m. at 140 observed landslide locations. FS 

TiVaSS 3D subsurface � ow

TiVaSS 1D in�ltration

2 7th 7 a .m.

27 th 8 a. m.

27 th 9 a .m.

2 7th 10 a .m.

2 7th 7 a. m.

27 th 8 a. m.

2 7t h 9 a .m.

2 7th 10 a .m.

FS  < 1 (Unstable) 1  <  FS  < 1.25 (Quasi stable) 1.25  <  FS  < 1.5 (Moderately Stable) FS  > 1.5 (Stable)

0 8 48 84

0 17 52 71

65 43 14 18

95 18 12 15

0 8 48 84

0 15 54 71

64 46 13 18

93 20 13 14

Sum

140

140

140

140

140

140

140

140

stable areas estimated by TiVaSS. The input data or parameters such as a constant soil depth, soil cohesion, soil friction angle, and initial condition are considered as reasons for this disagreement. The slope-stability classes obtained by TiVaSS with the 1D-in �ltration option are listed in Table 2. The results are similar, but TiVaSS with the 3D subsurface  � ow outperforms that with 1D in�ltration. The accuracy of landslide-prediction models is typically evaluated by comparing observed landslides with the predicted results obtained using the model. According to previous studies ( Crosta and Dal Negro, 2003;  Salciarini et al., 2006; Vieira et al., 2010; Raia et al., 2014), the agreement of cells between predicted and actual landslides can be used to evaluate the performance of  landslide-forecasting models. However, there are no precise landslide observation data for the 2011 Mt. Umyeon landslide event. The only available observations are of the landslide-initiation locations shown in Fig. 4. The aforementioned method is dif �cult to apply with landslide observation data points. It is expected that the denser grid causes a lower accuracy with the cell-agreement method because the observed landslide areas become narrow

with point data contained within only one pixel for any grid resolution. Therefore, in this study, the landslide ratio of each predicted FS   class (LR class) proposed by  Park et al. (2013) is implemented to estimate the model accuracy as follows:

LR class ¼

% of observedslopefailurelocations in each class of FS : % of predictedslope failure areas in each class of FS (23)

Table 3 shows the results of the LR class  analysis. 19.257% of the area is classi�ed as unstable, and 67.86% of the actual slope failures are correctly predicted. The percentage of LR FS<1  is 72.54%, which means that if a landslide occurs, the predicted unstable area (FS < 1) has a 72.54% chance of including the actual landslides. Table 4 shows the results of the LR class  analysis with 1D in�ltration. The performance of the 3D subsurface  � ow outperformed that with 1D in�ltration in this landslide simulation, although the difference is not critical. The Mt. Umyeon landslide was an extreme event, and all of the mountainous area became saturated, limiting

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

181

 Table 3 LR class  of results obtained with the 3D subsurface  � ow option at 10 a.m. on the 27th. FS classes

Slope failure location (a)

% of location (c )  ¼  a/b

% of predicted area (d)

LR class (e) ¼  c /d

% of LR class  ¼  e/ f 

FS  < 1 1  <  FS  < 1.25 1.25  <  FS  < 1.5 1.5  <  FS 

95 18 12 15

67.86 12.86 8.57 10.71

19.25 20.81 18.79 41.15

3.52 0.61 0.46 0.26

72.54 12.72 9.39 5.36

sum

140 (b)

100

100

4.86 ( f )

100

 Table 4 LR class  of results obtained with the 1D in �ltration  � ow option at 10 a.m. on the 27th. FS classes

Slope failure location (a)

% of location (c )  ¼  a/b

% of predicted area (d)

LR class (e) ¼  c /d

% of LR class  ¼  e/ f 

FS  < 1 1  <  FS  < 1.25 1.25  <  FS  < 1.5 1.5  <  FS 

93 20 13 14

66.43 14.29 9.29 10.0

19.26 20.88 18.81 40.06

3.45 0.68 0.49 0.24

70.81 14.05 10.14 5.00

sum

140 (b)

100

100

4.87 ( f )

100

Fig. 10.   TiVaSS plots of 1D and 2D data.

the effect of horizontal �ow. It is expected that a partially saturated area or a longer simulation period yields a more signi �cant difference between the 3D subsurface �ow and the 1D in�ltration

method. Nevertheless, if the local slope is mild and the simulation period is relatively short, the 1D-in�ltration method has an advantage with regard to the computational cost.

182

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183

The computation times of 3D subsurface �ow and 1Din�ltration approaches are 931 s and 97sec, respectively, on a laptop computer with intel core i5 2.20 GHz CPU. As stated in Section 3.2, Park et al. (2013) investigated the same landslide event using TRIGRS, which estimated 33.33% of landslide locations and 70.3% of LR class  with similar simulation condition. Although a direct comparison of two models is dif �cult because of  differences of model structures, initial condition, soil water retention model, and etc. it can be found that the performance of TiVaSS is comparable to that of TRIGRS.  3.4. Visualization in TiVaSS  Visualization in TiVaSS is performed using the QCustomPlot library. The simulation results, i.e., the FS , are plotted at the center of  the main windows during the simulation, as shown in  Fig. 3. The results at a speci�c timecan beplottedin the result tab atthe top of  the main windows. Fig. 10 shows an example of plotting 1D and 2D data using TiVaSS. Fig.10(a) and (b) show the rainfall and soil-water retention curves, respectively.   Fig. 10(c) and (d) have the same interface, which supports a 2D map with several color gradients. Simple x and y point data and a polygon or polyline �le (ESRI shape format) can be plotted in a 2D map, as show in Fig. 10(d). In addition, if the point data are loaded on an FS map, the FS class and LR class, as shown in Tables 2 and 3, can be calculated in TiVaSS. The 2D map supports the asc, png, and jpg formats for saving  � les.

research scope.  Acknowledgment This study was supported by a grant (2014-2061-01) from Chungnam National University in 2015. References An, H, Yu, S., 2012. Finite volume integrated surface-subsurface  � ow modeling on non-orthogonal grids. Water Resour. Res. 50, 2312 e2328.  http://dx.doi.org/ 10.1002/2013WR013828. An, H., Yu, S., 2014. Finite volume integrated surface-subsurface  � ow modeling on nonorthogonal grids. Water Resour. Res. 50 (3), 2312 e2328. http://dx.doi.org/ 10.1002/2013WR013828. Aleotti, P., 2004. A warning system for rainfall-induced shallow failures. Eng. Geol. 73, 247e265. Alvioli, M., Baum, R.L., 2016. Parallelization of the TRIGRS model for rainfall-induced landslides using the message passing interface. Environ. Model. Softw. 81, 122e135.  http://dx.doi.org/10.1016/j.envsoft.2016.04.002. Apip, Takara K., Yamashiki, Y., Sassa, K., Agung, A.B., Fukuoka, H., 2010. A distributed hydrological egeotechnical model using satellite-derived rainfall estimates for shallow landslide prediction system at a catchment scale. Landslides 7 (3), 237e258.  http://dx.doi.org/10.1007/s10346-010-0214-z. Baum, R.L., Godt, J.W., Savage, W.Z., 2010. Estimating the timing and location of  shallow rainfall-induced landslides using a model for transient, unsaturated in�ltration. J. Geophys. Res. Earth Surf. 115 (F3), F03013.  http://dx.doi.org/ 10.1029/2009JF001321. Baum, R.L., Savage, W.Z., Godt, J.W., 2008. TRIGRS a Fortran Program for Transient Rainfall In�ltration and Grid-based Regional Slope-stability Analysis, Version 2.0. U.S. Geological Survey, Reston, Va. Retrieved from. http://pubs.usgs.gov/of/ 2008/1159. Bezak, N., Sraj, M., Miko s, M., 2016. Copula-based IDF curves and empirical rainfall thresholds for �ash �oods and rainfall-induced landslides. J. Hydrol.   http:// dx.doi.org/10.1016/j.jhydrol.2016.02.058  (in press). Borga, M., Fontanaa, G.D., Cazorzi, F., 2002. Analysis of topographic and climatic control on rainfall-triggered shallow landsliding using a quasi-dynamic wetness index. J. Hydrol. 268, 56 e71. Ciabatta, L., Camici, S., Brocca, L., Ponziani, F., Stelluti, M., Berni, N., Moramarco, T., 2016. Assessing the impact of climate-change scenarios on landslide occurrence i n Umb ri a R egi on , I ta ly. J . Hyd rol .   http://dx.doi.org/10.1016/j.jhydrol.2016.02.007  (in press). Crosta, G.B., Dal Negro, P., 2003. Observations and modelling of soil slip-debris  � ow initiation processes in pyroclastic deposits: the Sarno 1998 event. Nat. Hazards Environ. Syst. Sci. 3, 53 e69. Dickson, M.E., Perry, G.L.W., 2016. Identifying the controls on coastal cliff landslides using machine-learning approaches. Environ. Model. Softw. 76, 117 e127. Dietrich, W.E., Montgomery, D.R., 1998. SHALSTAB: a Digital Terrain Model for Mapping Shallow Landslide Potential. NCASI (National Council of the Paper Industry for Air and Stream Improvement) Technical Report, 1998 . Gardner, W., 1958. Some steady-state solutions of the unsaturated moisture �ow equation with application to evaporation from a water table. Soil Sci. 85, 228e232. Guzzetti, F., Carrara, A., Cardinali, M., Reichenbach, P., 1999. Landslide hazard evaluation: a review of current techniques and their application in a multi-scale study, Central Italy. Geomorphology 31 (1 e4), 181e216. Guzzetti, F., Reichenbach, P., Cardinali, M., Galli, M., Ardizzone, F., 2005. Probabilistic landslide hazard assessment at the basin scale. Geomorphology 72 (1 e4), 272e299. Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., Galli, M., 2006. Estimating the quality of landslide susceptibility models. Geomorphology 81 (1 e2), 166e184. Guzzetti, F., Peruccacci, S., Rossi, M., Stark, C.P., 2007. Rainfall thresholds for the initiation of landslides in central and southern Europe. Meteorol. Atmos. Phys. 98, 239e267. Guzzetti, F., Peruccacci, S., Rossi, M., Stark, C.P., 2008. The rainfall intensity eduration control of shallow landslides and debris  � ows: an update. Landslides 5, 3 e17. Iverson, R.M., 2000. Landslide triggering by rain in �ltration. Water Resour. Res. 36 (7), 1897e1910. http://dx.doi.org/10.1029/2000WR900090.  Jie, Y., Jie, G., Mao, G., Li, G., 2004. Seepage analysis based on boundary- �tted coordinate transformation method. Comput. Geotech. 31, 279 e283. Korean Society of Civil Engineering, 2011. Research Contract Report: Causes Survey and Restoration Work of Mt. Woomyeon Landslide. Seoul . Liao, Z., Hong, Y., Kirschbaum, D., Adler, Robert F., Gourley, Jonathan J., Wooten, R., 2011. Evaluation of TRIGRS (transient rainfall in �ltration and grid-based regional slope-stability analysis) s predictive skill for hurricane-triggered landslides: a case study in Macon County, North Carolina. Nat. Hazards 58 (1), 325e339.  http://dx.doi.org/10.1007/s11069-010-9670-y. Liu, C.N., Wu, C.C., 2008. Mapping susceptibility of rainfall-triggered shallow landslides using a probabilistic approach. Environ. Geol. 55 (4), 907 e915. http:// dx.doi.org/10.1007/s00254-007-1042-x . Luan, T., Long, N., Shibayama, M., Cannata, M., 2010. Area Informacitc and TRIGRS 

4. Conclusions A new time-variant slope stability for landslide prediction is developed. A 3D subsurface  � ow model is coupled with an in �nite slope-stability model to consider the effect of horizontal water movement in the subsurface. To overcome the large computational requirements of the 3D subsurface  � ow module, partially implicit temporal discretization and the simpli�cation of  � rst-order spatial discretization are proposed and applied in TiVaSS. A GUI and 2D data visualization aresupported with the Qt framework. Thewidely used ESRI arc �le format is supported for the input and output �les. The performance of TiVaSS is validated against data from the 2011 Mt. Umyeon landslide. Despite many uncertainties, such as the model structure, valid soil depth, and soil cohesion, the performance of TiVaSS was satisfactory overall. Among 140 landslide initiation locations, 78 locations were judged as unstable, and 34 locations were judged as quasi-stable, which means that minor destabilizing factors can lead to instability. The time of occurrence of the landslide was successfully reproduced by TiVaSS. The performance of the 1D-in�ltration option was also investigated. The 3D subsurface �ow module and the spatially uniform 1D in�ltration module can be optionally selected in TiVaSS. The 3D subsurface �ow showed a slightly better performance than the 1D in �ltration, whereas the overall FS  distributions were very similar. If the local slope is mild or the simulation period is relatively short, 1D in �ltration is an attractive option with regard to the computational cost. The main advantage of TiVaSS is the numerical modeling of the 3D subsurface �ow, which allows more accurate and �exible simulation. Further applications for several landslide events are ongoing, and additional functions such as the estimation of the landslide-prediction accuracy and overland �ow will be added, which is not currently considered in TiVaSS. The overland  � ow or debris �ow could be included in the next version of TiVaSS. The Monte Carlo approach for TiVaSS with a 1D in �ltration module (Raia et al., 2014; Peres and Cancelliere, 2016) and the parallelization of the model (Alvioli and Baum, 2016) are also in our future



’

H. An et al. / Environmental Modelling & Software 85 (2016) 172e183 model for study shallow landslide vulnerability assessment. A case study in Bavi area, Hanoi region, Vietnam. In: GeoInformatics for Spatial-infrastructure Development in Earth and Allied Sciences (GIS-IDEAS) 2010. Hanoi University of Technology. Retrieved from. http://gisws.media.osaka-cu.ac.jp/gisideas10/ viewabstract.php?id¼369. Lu, N., Likos, W.J., 2006. Suction stress characteristic curve for unsaturated soil.  J. Geotech. Geoenvironmental Eng. 132 (2), 131 e142.  http://dx.doi.org/10.1061/ (ASCE)1090-0241(2006)132:2(131) . Lu, N., Wayllace, A., Oh, S., 2013. In �ltration-induced seasonally reactivated instability of a highway embankment near the Eisenhower Tunnel, Colorado, USA. Eng. Geol. 162, 22 e32.  http://dx.doi.org/10.1016/j.enggeo.2013.05.002. Montgomery, D.R., Dietrich, W.E., 1994. A physically based model for the topographic control on shallow landsliding. Water Resour. Res. 30 (4), 1153 e1171. http://dx.doi.org/10.1029/93WR02979 . Montrasio, L., Valentino, R., Losi, G.L., 2011. Towards a real-time susceptibility assessment of rainfall-induced shallow landslides on a regional scale. Nat. Hazards Earth Syst. Sci. 11 (7), 1927 e1947.   http://dx.doi.org/10.5194/nhess-111927-2011. Pack, R.T., Tarboton, D.G., Goodwin, C.N., 1999. SINMAP 2.0-A Stability Index Approach to Terrain Stability Hazard Mapping, User's Manual. Retrieved from. http://works.bepress.com/robert_pack/8. Park, D.W., Nikhil, N.V., Lee, S.R., 2013. Landslide and debris �ow susceptibility zonation using TRIGRS for the 2011 Seoul landslide event. Nat. Hazards Earth Syst. Sci. Discuss. 1 (3), 2547 e2587.  http://dx.doi.org/10.5194/nhessd-1-25472013. Papathoma-K ohle, M., Zischg, A., Fuchs, S., Glade, T., Keiler, M., 2015. Loss estimation for landslides in mountain areas e  an integrated toolbox for vulnerability assessment and damage documentation. Environ. Model. Softw. 63, 156 e169. http://dx.doi.org/10.1016/j.envsoft.2014.10.003 . Peres, D.J., Cancelliere, A., 2016. Estimating return period of landslide triggering by Monte Carlo simulation. J. Hydrol.   http://dx.doi.org/10.1016/j.jhydrol.2016.03.036  (in press). Pham, B.T., Pradhan, B., Tien Bui, D., Prakash, I., Dholakia, M.B., 2016. A comparative study of different machine learning methods for landslide susceptibility assessment: a case study of Uttarakhand area (India). Environ. Model. Softw. 84, 240e250. Pradhan, B., Lee, S., 2010. Landslide susceptibility assessment and factor effect analysis: backpropagation arti �cial neural networks and their comparison with frequency ratio and bivariate logistic regression modelling. Environ. Model. Softw. 25, 747 e759.  http://dx.doi.org/10.1016/j.envsoft.2009.10.016. Raia, S., Alvioli, M., Rossi, M., Baum, R.L., Godt, J.W., Guzzetti, F., 2014. Improving €

183

predictive power of physically based rainfall-induced shallow landslide models: a probabilistic approach. Geosci. Model Dev. 7, 495 e514. Ruhaak, W., Rath, V., Wolf, A., Clauser, C., 2008. 3D  � nite volume groundwater and heat transport modeling with non-orthogonal grids, using a coordinate transformation method. Adv. Water Resour. 31, 513 e524. Saadatkhah, N., Kassim, A., Lee, M., 2014. Hulu Kelang, Malaysia regional mapping of  rainfall-induced landslides using TRIGRS model. Arabian J. Geosciences 1 e12. http://dx.doi.org/10.1007/s12517-014-1410-2. Salciarini, D., Godt, J.W., Savage, W.Z., Conversini, P., Baum, R.L., Michael, J.A., 2006. Modeling regional initiation of rainfallinduced shallow landslides in the eastern Umbria Region of central Italy. Landslides 3, 181 e194. Simoni, S., Zanotti, F., Bertoldi, G., Rigon, R., 2008. Modelling the probability of  occurrence of shallow landslides and channelized debris  � ows using GEOtopFS. Hydrol. Process. 22 (4), 532 e545.  http://dx.doi.org/10.1002/hyp.6886. Sirangelo, B., Versace, P., Capparelli, G., 2003. Forewarning model for landslides triggered by rainfall based on the analysis of historical data �le. In: Hydrology of  the Mediterranean and Semiarid Regions, Proceedings International Symposium Held at Montpellier, 1Al IS Publ., 278, pp. 298 e304. Soeters, R., van Westen, C.J., 1996. Slope instability rec-ognition, analysis, and zonation. In: Turner, A.K., Schuster, R.L. (Eds.), Landslides Investigation and Mitigation. National Academy Press. National Research Council, Washington, D.C., pp. 129e177. Transportation Research Board Special Report 247. Son, S.H., Choi, B., Paik, J., 2012. Characteristics of rainfall and groundwater table in a small forested watershed of Mt. Umyeon. In: Proceedings of 38th Korea Society of Civil Engineering Conference and Civil Expo 2012, pp. 769 e772 (In Korean). Tarolli, P., 2014. High-resolution topography for understanding Earth surface processes: opportunities and challenges. Geomorphology 216, 295 e312. Tarolli, P., Tarboton, D.G., 2006. A new method for determination of most likely landslide initiation points and the evaluation of digital terrain model scale in terrain stability mapping. Hydrol. Earth Syst. Sci. 10, 663 e677. Terzaghi, K., 1943. Theoretical Soil Mechanics. J. Wiley . Vieira, B.C., Fernandes, N.F., Filho, O.A., 2010. Shallow landslide prediction in the Serra do Mar, S ao Paulo, Brazil. Nat. Hazards Earth Syst. Sci. 10, 1829 e1837. http://dx.doi.org/10.5194/nhess-10-1829-2010 . Yune, C.Y., Chae, Y.K., Paik, J., Kim, G., Lee, S.W., Seo, H.S., 2013. Debris �ow in metropolitan area  d  2011 Seoul debris  � ow. J. Mt. Sci. 10 (2), 199 e206. http:// dx.doi.org/10.1007/s11629-013-2518-7. Zi, T., Kumar, M., Kiely, G., Lewis, C., Albertson, J., 2016. Simulating the spatiotemporal dynamics of soil erosion, deposition, and yield using a coupled sediment dynamics and 3D distributed hydrologic model. Environ. Model. Softw. 83, 310e325. ~