Environmental Modelling & Software 105 (2018) 1e16

Contents lists available at ScienceDirect

Environmental Modelling & Software

journal homepage: www.elsevier.com/locate/envsoft

Integration of two-phase solid fluid equations in a catchment model

for flashfloods, debris flows and shallow slope failures
B. Bout a, *, L. Lombardo b, c, d, C.J. van Westen a, V.G. Jetten a
a
University Twente, Faculty of Geo-Information Science and Earth Observation (ITC), The Netherlands
b
CEMSE Division, King Abdullah University of Science and Technology, Thuwal, Jeddah, Saudi Arabia
c
PSE Division, King Abdullah University of Science and Technology, Thuwal, Jeddah, Saudi Arabia
d
Department of Geosciences, Tübingen University, Tübingen, Germany

a r t i c l e i n f o a b s t r a c t

Article history: An integrated, modeling method for shallow landslides, debris flows and catchment hydrology is
Received 15 May 2017 developed and presented in this paper. Existing two-phase debris flow equations and an adaptation on
Received in revised form the infinite slope method are coupled with a full hydrological catchment model. We test the approach on
20 March 2018
the 4 km2 Scaletta catchment, North-Eastern Sicily, where the 1-10-2009 convective storm caused debris
Accepted 23 March 2018
flooding after 395 shallow landslides. Validation is done based on the landslide inventory and photo-
graphic evidence from the days after the event. Results show that the model can recreate the impact of
both shallow landslides, debris flow runout, and debris floods with acceptable accuracy (91 percent
Keywords:
Spatial numerical modeling
inventory overlap with a 0.22 Cohens Kappa). General patterns in slope failure and runout are well-
Physically-based modeling predicted, leading to a fully physically based prediction of rainfall induced debris flood behavior in the
Flash floods downstream areas, such as the creation of a debris fan at the coastal outlet.
Shallow landslides © 2018 Elsevier Ltd. All rights reserved.
Debris flow
OpenLISEM

1. Introduction and surface properties. A large number of runout models exist,

varying both in modeling approach and used equations. Both one-
Shallow landslides are frequently occurring natural hazards, dimensional channel-simulations, full depth-varied grid methods
which may be triggered by extreme rainfall events, snow melt and and Smooth Particle Hydrodynamics have been used to estimate
earthquakes, and are characterized by a relatively small and debris flow behavior (Pudasaini, 2012; Huang et al., 2012). Other
shallow failure plane. Because the triggering of these landslides implementations, such as SCIDDICA S4c (D'Ambrosio et al., 2007),
takes place predominantly during intense precipitation, the sliding approximate debris flow behavior by means of cellular automata.
soil, mixed with water, may evolve into debris flows, that have a Depth averaged equations are used in a large variety of two-
devastating impact on villages, roads and other elements-at-risk. In dimensional models (Scheidl et al., 2013). A variety of depth aver-
order to understand and predict the behavior of debris flows, nu- aged finite-element models, such as Ramms (Bartelt et al., 2013),
merical models have been frequently used as both predictive and Flo-2D (O'Brien, 2006) and use a fixed volume as input for the
analytical tools. However, in current modeling approaches, pro- debris flow. Others, such as MassMove2D (Beguería et al., 2009),
cesses that relate to debris flows, such as hydrology, shallow Debris Mobility Model (Kwan and Sun, 2006) and AschFlow (Luna
landslides and runout, are mostly simulated separately. et al., 2015), include entrainment and the addition of water flow.
The simulation of debris flow dynamics is performed by debris The processes that cause shallow slope failures and their tran-
flow runout models. These models use (semi-) physically-based sition into debris flows are often also numerically simulated,
estimations of the internal forces in debris flows to numerically although empirical methods are also often used. Hydrological
calculate flow depths, velocities and routing based on topography models are frequently used to predict behavior of both surface and
sub-surface hydrological processes. Through flow simulations,
overland flow and the resulting infiltration patterns can be esti-
mated (van Beek, 2003). Similar to debris flows, hydrology is
* Corresponding author. University Twente, Faculty of Geo-Information Science
and Earth Observation (ITC), The Netherlands. simulated in spatially distributed numerical models such as GEOtop
E-mail address: [email protected] (B. Bout). (Rigon et al., 2006). From the available catchment-scale

https://doi.org/10.1016/j.envsoft.2018.03.017
1364-8152/© 2018 Elsevier Ltd. All rights reserved.
2 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

hydrological models, some are open-source, such as JGrass- A common obstacle in the modeling and prediction of shallow
NewAge (Formetta et al., 2014). Furthermore, hydrological models landslides and debris flows is the accuracy of input data and pa-
can predict soil water flow, and the resulting soil water has been rameters (van Westen et al., 2006; 2008; Mergili et al., 2011;
used as a direct proxy for slope stability (van Asch et al., 1999). Nikolopoulos et al., 2015). Physically-based models are often limited
Slope stability models are used to predict both volume and timing in accuracy by the spatial resolution of elevation and soil data. Slope
of shallow landslides. A variety of regional slope stability models instability caused by structural weakness provides a particularly
have been developed based on the limit equilibrium, which uses large challenge for current modeling methods since the sub-surface
finite elements to estimate the forces acting on a failure plane. The structure is generally unknown over larger areas (van Westen et al.,
infinite slope models furthermore assume that inter-cell forces can 2006). The availability of accurate soil data, both in terms of their
be neglected, and structural finite elements are used to calculate spatial variation in type and thickness, and associated geotechnical
the local Factor of Safety (SoF). This resulting FoS, which depends and hydrological properties can be a major limiting factor. An inte-
on a soil description and soil water behavior, are successfully grated approach to model slope failure, debris flows and hydrology
applied as a prediction for slope failure in a variety of regions could only be used when sufficient input data is available. In the past
around the world (Van Beek & Van Asch, 2004; Kuriakose et al., decade however, the data problem has become less severe due to
2009a,b; Mergili et al., 2014a). Recent approaches, such as the increasing availability of detailed data. High-resolution elevation
method by Mergili et al. (2014b), are focused on providing esti- products such as Lidar DEMs have become widely available (Tarolli,
mation of failure probability by using more accurate estimates of 2014). Major improvements have been made in estimating soil data
stability based on three dimensional rotational slip surface analysis from various sources such as national soil maps and satellite data
following Xie et al. (2004). A similar approach, but instead using (Hengl et al., 2017). Furthermore, soil depth estimations based on
Bishops Simplified method, is used in Scoops3D to find approxi- statistical correlation of topographical parameters have shown
mate failure volumes for rotational slides (Reid et al., 2015). increasing accuracy. Lastly, when slope failure is not estimated an
Specialized models are thus frequently used to investigate shallow initial volume is often used for runout calculation. In this case, the
landslides and debris flows. influence of hydrology and other processes on debris flow runout is
Despite all the currently available models, combined approaches still neglected. Thus, an integrated model can provide further im-
to debris flow initiation and runout are scarce. Empirical initiation provements compared to traditional methods even when insuffi-
and runout models, such as Flow-R combine landslide failure cient data is available for prediction of slope failure.
criteria and possible flow patterns to create susceptibility maps The objective of this paper is to develop and test an integrated
(Horton et al., 2013). However, their method is semi-physical in model to analyze the influence of rainfall-triggered shallow land-
nature and doesn't include failure volumes, rheology nor entrain- slides and debris flows in a hydrological catchment model. Slope
ment. Fan et al. (2017) introduced an innovative way to approach failures are estimated by using an adaptation on the classic infinite
catchment-scale landslide modeling. For their approach, the failure slope stability method. To simulate debris flows, the two-phase
plane is static, and defined as the lithic contact. Furthermore, hy- generalized debris flow equations by Pudasaini (2012) are imple-
drological processes are still generally neglected, and the stability mented. These methods are included in the OpenLISEM model
and runout models are coupled through a one-way link. Finally, (Bout and Jetten, 2018). In order to test the performance of the
Mergili et al. (2011), provided a more in-depth analysis of the what aforementioned model, we attempt to model the impact of a
sort of approaches to integrated simulation of debris flows could be convective storm that hit the south-eastern coast of Sicily in 2009
taken. Their analysis covered wetting front based slope stability and (Lombardo et al., 2015, 2018a). Finally, several alternatives of the
failure in a small catchment, and simulation of erosion by runoff, developed modeling method are tested and a sensitivity analysis is
both leading to debris flow runout. The authors solved runout by performed.
using a two-parameter semi-deterministic frictional model routed
over the terrain. While the currently available models provide 2. Materials and methods
useful investigative and predictive tools, integrated simulations
using fully physically based descriptions of all related processes 2.1. Schematic model description
could still increase understanding and usability of numerical
simulations. Combining the discussed methods, we created a modeling
Most existing hydrological, slope stability and debris flow run- method that incorporates both hydrological processes such as
out models focus on specific processes without the possibility of rainfall, interception, infiltration, flow, and morphological pro-
interactions or feedbacks between the processes. However, in cesses such as shallow landslides, slope failure, and landslide run-
practice processes such as slope failure, flow directions, infiltration out. A simplified flow chart for the final model is shown in Fig. 1.
and flow properties such as viscosity all influence each other. When sediment components are absent, the model reduces to a
Flooding and debris flows in particular have frequent interactions fully functioning hydrological catchment model.
due to their common metrological trigger. In many cases, in-
teractions between debris flows and flooding substantially influ- 2.2. Model basis
ence the behavior of both processes. When these phenomena are
neglected, the predictive power of models is substantially limited. A In order to integrate the occurrence of shallow landslides, debris
major example of these interactions are blocked rivers or drainage flows and flash flooding within a single model, we used the Open
channels by debris flows, resulted in alternating waves of debris Source Limburg Soil Erosion Model (OpenLISEM) as a basis. Open-
flows and flooding (Tang et al., 2012; Adegbe et al., 2013; Luna et al., LISEM is a physically-based numerical model with the purpose of
2014). Debris flows can also interact with overland flow causing event-based runoff, flooding and erosion modeling on a catchment
decreased viscosity. Nguyen et al. (2013) found that dilution of a scale. LISEM is fully spatially distributed and uses a topography-
debris flow by directed overland flow caused runout over a larger following grid to solve both cell specific processes, and the differ-
area, including the streets of a nearby village. In order to increase ential equations governing flow.
the understanding of hydrology, shallow landslides and the debris The OpenLISEM model implements multiple types of infiltration
flows that are caused by these, a holistic and integrated approach models such as Smith and Parlange (1978) and the SWATRE full
should be considered. vertical soil water balance model (Bastiaanssen et al., 1996). The
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 3

Fig. 2. The input data layers for OpenLisem.

     
vhux v hu2x v hux uy
þ þ ¼ gh Sx  Sf ;x (3)
vt vx vy
 
   
vhuy v hu2y v hux uy
þ þ ¼ gh Sy  Sf ;y (4)
vt vy vx

Where h is the flow height (m), u is the flow velocity (m s1 ), R is

the rainfall (m), I is the infiltration (m), g is the gravitational ac-
celeration (m s2 ), S is the friction term (m s2 ) and Sf is the mo-
mentum source term (m s2 )
 
Fig. 1. A simplified flow chart for the new OpenLISEM model. gh2
d 2
Sx ¼   ghSfx (5)
dx
simulations in this paper use the Green & Ampt infiltration model,  2
which assumes a wetting front moving down into the soil due to d gh2
infiltrating rainfall (Green and Ampt, 1911). The resulting potential Sy ¼   ghSfy (6)
dy
infiltration is subtracted from the available surface water (Equation
(1)). The friction slope terms, which are the friction forces divided by

 the water height and the gravitational acceleration, can be calcu-
q  qi lated using the Darcy-Weisbach friction law (Chow, 1959) (Equation
fpot ¼ Ks j s þ1 (1)
F (7)).

Where fpot is the potential infiltration rate (m s1 ), F is the cu- !!

g ujuj
mulative infiltrated water (m), qs is the porosity (m3 m3 ), qi is the
Sf ¼ (7)
n2 h43
initial soil moisture content (m3 m3 ), j is the matric pressure at
the wetting front (m) and Ks is the saturated conductivity (m s1 ). Where n is the Manning's n friction coefficient (s m3 ).
1

Input data consists of soil, land surface and terrain properties

(Fig. 2). Surface properties such as buildings, roads and channels are
defined as fraction of a cell's surface. The hydrological processes that 2.4. Slope stability
are simulated within OpenLisem are extensive and include inter-
ception by vegetation, surface micro-ponding and dynamic flow. The implemented method for estimating slope stability is based
Further details on the underlying physical principles of OpenLISEM on the infinite-slope method. In this method both the local
can be found in Baartman et al. (2012) and Jetten and De Roo (2001). downslope and local resisting forces are calculated. Downslope
force is based on the assumption that the failure plane is parallel to
2.3. Flow equations the surface plane. The weight of the soil section is calculated by
using soil density, porosity and effective soil saturation. The effec-
Overland flow within OpenLISEM is a combination of runoff and tive soil saturation is the average of the saturation above the wet-
flooding. All flow computations are based on the full Saint-Venant ting front and below the wetting front (Equation (8)).
equations for shallow flow (equations (2)e(6)).
W ¼ ðg  mgw Þz (8)
 
vh vðhux Þ v huy
þ þ ¼ RI (2) Where W is the weight of the soil column (N), g is the soil density
vt vx vy
(kg m3 ), gw is the density of water (kg m3 ), z is the soil depth (m)
and m is the ratio between the depth of the saturated zone and the
4 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

soil depth. row of soil sections in pixels, the total of force capacity and demand
Force capacity is based on the Mohr-Coulomb failure criterion must therefore be taken into account. However, we assume there is
(Equation (9)) no information about subsurface structure. There is therefore, no
 0 knowledge about subsurface forces and any possible plane that lies
t ¼ c þ N tan f (9) parallel to the surface is allowed to act as failure plane. To overcome
this, we utilize another approach. For every group of unstable cells
Where t is the shear stress, c is the cohesion (kPa), N the normal at any moment, a toe can be defined between the lowest unstable
0
force acting on the soil column (N), f the internal friction angle ( ). cells, and the stable cells underneath those. Here the downstream
Within the soil description, a wetting front is assumed, together cell provides enough force capacity and the upstream cell does not.
with a homogeneous initial soil moisture content (Fig. 3). Therefore, a force equilibrium must be present between these cells.
Combining the provided equations, the factor of safety is calculated A visual representation of such as situation is provided in Fig. 4,
based on simple raster-based soil data (Equation (10)). where Fd is the force demand, and Fc is the resisting force.
At the toe location, slope failure is initiated for the upslope cell.
0  0 Based on this assumption, the safety factor equation is solved to
c þ c þ ððg  mgw Þz þ mgw zÞcosðbÞ2 tan f
SF ¼ (10) find the depth of the remaining soil that is required for stability.
ððg  mgw ÞzÞsinðbÞcosðbÞ
This is possible since we assume any plane parallel to the surface
could be the failure plane. Any material above this stable depth is
Where SF is the safety factor (), b is slope of the soil section (), c
0 indicated as slope failure. Due to the slope failure, the slope of the
is the cohesion of the soil (kPa), c is the apparent cohesion of the
0 surrounding cells changes. In particular, the upslope cells experi-
soil (kPa) and f is internal friction angle of the soil (). Here the
ence an increase in slope, leading to a decrease in stability. We thus
apparent cohesion consists of additional root cohesion and a matric
repeat the process of finding the stable depth, and initiate slope
suction term.
failure repeatedly, until no unstable cells are left. This iterative
When the safety factor decreases below 1, the soil section is
method of finding stable soil depth, provides an estimation of
unstable and failure depth will be calculated. Within the safety
failure volume.
factor equation, stability decreases when soil depth increases.
In reality, the change in stability is not caused by changing slope,
Calculations for the entire soil depth thus provide the lowest
but rather by a change in force propagation through the subsurface.
possible safety factor. Therefore, by calculation stability for the
Any cell directly upslope of a failure could experience a change in
complete soil column, failure of smaller surfaces is similarly
resisting force. Thus, the stability of the upslope cell changes. Since
detected.
there is no detailed knowledge about the forces that propagate
through the subsurface, we use the assumption that the failure
2.5. Slope failure plane must be parallel to the surface. When the local slope changes,
the local failure plane direction and thus stability similarly changes.
When a raster element with a particular slope is determined to An example of results from the iterative slope failure method is
be unstable based on equation (10), the next step is to determine shown in Fig. 5. The safety factor of the fourth cell indicates
the failure volume. Within the infinite-slope method, several as- instability. A stable soil depth is calculated and the slope is altered.
sumptions should be made to estimate failure volume. In reality, Based on the altered slope, the safety factor is recalculated for
forces propagate through the subsurface, and stability should be upslope cells. In this particular example, three upslope cells
calculated through the method of slices, using circular or irregular become unstable, and slope failure is initiated in four iterations.
sliding surfaces with methods such as Fellenius (1936), Janbu After the iterative process, the full failure volume is instanteneously
(1975) and Morgenstern and Price (1965). For the stability of a added to the flow equations.

Fig. 4. A force diagram for the landslide toe. An equilibrium point must exist imme-
Fig. 3. An example of a simplified physical description of the soil layer. diately downstream of the toe.
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 5

Fig. 5. A 2D example of the result of using iterative slope failure with a finite element factor of safety.

2.6. Debris flow equations 8 2

<
vb

1 v h2

vb 1 v2 vf
To simulate the flow dynamics and interactions of floods and Sy;f ¼ af g  ε4 Pbf þ Pbf  2 2
: vy h vy 2 vy af NR vy
nun-uniform debris flows, a set of two-phase equations is required.
!


An extensive set of two-phase debris flow equations are available v2 uf v2 v f cvf 1 v vas  
from Pudasaini (2012). This set of equations contains a physically þ þ  þ 2 v  vs
vyvx vy2 ε2 h2 af NR vy vy f
based two-phase momentum balance. Besides pressure and grav-
 39
itational forces, it includes viscous forces, non-Newtonian viscosity,
 xas uf  us =
v vas   va 
s 5
two-phase drag and a Mohr-Coulomb type friction force for the þ uf  us þ v  vs  2
solid phase (Equations (11)e(14)). Based on the current and local vy vy vx f ε af NRA h2 ;
state of flow, forces increase in magnitude. This approach allows for   
1 ! ! j1
a smooth transition between non viscous flow, hyper concentrated  CDG uf  us  u f  u s 
g
streamflow and debris flows. Furthermore, the interactions be-
tween distinct flow types are automatically solved. (14)

Where as and af are the volume fraction for solid and fluid phases
(), Pb is the pressure at the base surface (Kg m1 s2 ), b is the basal




surface of the flow (m), NR is the Reynolds number (), NRA is the
vb us   vb vh
Sx;s ¼ as g  ! tan vPbs  εPbs  εas gPbf quasi-Reynolds number (), CDG is the drag coefficient (), rf is the
vx j u sj vx vx
 j1 density of the fluid (kg m3 ), rs is the density of the solids (kg m3 ),
vb  
! !  g is the density ratio between the fluid and solid phase(), c is the
þ þ CDG uf  us  u f  u s 
vx vertical shearing of fluid velocity (m s1 ), ε is the aspect ratio of the
(11) model (), x is the vertical distribution of as (m1 ).
Within the momentum equations, the Reynolds number is used
to scale turbulent and viscous forces.
Here, we follow Pudasaini (2012) for the definitions of the




normal and quasi-Reynolds number (Equations (15) and (16).
vb vs   vb vh
Sy;s ¼ as g  ! tan vPbs  εPbs  εas gPbf
vy j u sj vy vy pffiffiffiffiffi
  gLHrf
vb   NR ¼ (15)
! ! j1
af h
þ þ CDG vf  vs  u f  u s 
vy
(12) pffiffiffiffiffi
gLH rf
NRA ¼ (16)
Ah

8 2 Where L is the length scale of the flow (m), H is the height of the
<
vb

1 v h2

vb 1 v2 uf flow (m), h is the viscosity (kg s1 m1 ) and A is the mobility of the
Sx;f ¼ af g  ε4 Pbf þ Pbf  2 2
: vx h vx 2 vx af NR vx interface ().
! To apply these two-phase equations successfully in a



v2 vf v2 u f cuf 1 v vas  catchment-based model, we replace the frictional force for the
þ þ  þ 2 u  u s water phase with the Darcy-Wiesbach equation for water flow
vyvx vy2 ε2 h2 af NR vx vx f

 3 9 friction. To complete the set of equations that govern debris flow-

 xas vf  vs =
v vas   va 
s 5
dynamics, several flow properties are estimated based on the
þ v  vs þ uf  us  2 volumetric sediment content. Viscosity is based on an empirical
vy vx f vy ε af NRA h2 ;
relation by O'Brien et al. (1993) (Equation (17)).
  
1 ! ! j1
 CDG uf  us  u f  u s  h ¼ aebas
g (17)
(13)
Where a the first viscsosity parameter () and b the second vis-
cosity parameter ().
6 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

The drag coefficient is based on the relation provided by

Pudasaini (2012) (Equations (18)e(21)).
 
r
af as 1  rfs
CDG ¼ (18)
εUT ðPFðReÞ þ ð1  PÞGðReÞÞ

rf dUT
Re ¼ (19)
hf

3
rf af
F¼ Re (20)
180 rs as

MðReÞ1
G ¼ af (21)

Where Re is the Reynollds number (), d is the median grain

diameter (), UT the settling velocity (m s1 ) and M is an empirical
parameter depending on the Reynolds number ().
Finally, the settling velocity of small (d < 100 mm) grains is
estimated by Stokes equations for a homogeneous sphere in water.
For larger grains, the equation by Zanke (1977) is used (Equation
(22)).
0v
u
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ffi 1
u h 2 Bu ðrs  rf Þ 3
0:01 rf gd C
rf Bu C
ws ¼ 10 Bt1 þ h  1C (22)
d @ r A
f

2.7. Application of the model

The described methods were applied in a test area to model the

initiation and runout of the debris flows that were caused by the 1
October, 2009 convective storm hitting the North-Eastern part of
Sicily, Italy. This event had a maximum rainfall intensity of 120 mm/
h and a return period of 1 in 30 years (Aronica et al., 2012). Rainfall
data for the S. Stefano Briga and Fiumedinisi rainfall stations are
Fig. 7. Input data for the Scaletta catchment. Location (top) Elevation model and soil
texture (middle), landslide inventory and land use map (bottom).

shown in Fig. 6, being provided by the local administration. The

event caused landslides, debris flows and flash flooding in the area
between St. Stefano Briga and Fiumedinisi. A statistical analysis of
the landslides and causal factors was carried out by Lombardo et al.
(2016) and Trigila et al. (2015). Furthermore, debris flow simula-
tions have been performed for several villages that were hit by
debris flows during this event (Cama et al., 2017; Lupiano et al.,
2016; Stancanelli et al., 2013).
For this study, the Scaletta catchment, a 4.3 km2 area between
St. Stefano Briga and Fiumedinisi has been selected as study site
(Fig. 7). A landslide inventory for the event is available, based on
areal images (Lombardo et al., 2014). A total of 395 shallow land-
slides took place during or directly after the event in the selected
area. Landslides in the neighboring catchment consisted mostly of
the layer of weathered surface material, between 0.5 and 2.5 m in
depth (Lombardo et al., 2015). The inventory does not distinguish
between source and deposition areas of each individual shallow
landslide. Instead, the full impact area is delineated by the poly-
gons. Fig. 5 illustrates the elevation, landslide inventory, land use
types and soil types. This coastal region is characterized by steep
Fig. 6. Rainfall data for the 1-10-2009 rainfall event in south-west Sicily. Temporal slopes, and three small channels that drain towards the sea.
resolution of the rainfall is 10 min. Source: (http://www.osservatorioacque.it/).
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 7

Infrastructure, in the form of a highway and railway lines are pre- case, the final statistical model took the form of equation (24).
dominantly located parallel to the coast, flanked with a string of
small settlements. The highway is partly tunneled, or built on SD ¼ Intersect þ a*dcoast þ b*dchannel þ c*s þ d*c (24)
bridges, and does not actively block any flow. Land use types consist
predominantly of shrubs/sparse forest mixed with areas of agri- Where SD is the soil depth (m), a,b,c,d are constants () derived
culture. Small patches of forest are present in the inland regions. from statistical analysis, dcoast is the distance to the coast (m),
The most common soil is of colluvial origin and consists of sandy dchannel is the distance to the nearest channel (m), s is the slope of
soils upstream, and clayey soils near the coast. These are the result the surface (m m1 ) and c is the profile curvature of the surface
of pedogenetic processes taking place on a medium to high grade (m1 ).
metamorphic parent rock primarily consisting of paragneiss, Based on the statistical correlation, the values for the constant
micaschists, phyllites and meta-arenites (Cama et al., 2015). Their are provided in Table 3.
weakening, transport and deposition in the form of unstable un- The fitted relationship between predicted and estimated soil
consolidated soils is further facilitated by a long history of depth has a slope of 0.90 and R2 of 0.69 (Fig. 8). Due to the low
compressional and more recent extensional tectonic regimes number of available soil depth measurements, the R2 is relatively
(Somma, 2006). low. However, the trend line between predicted and measured soil
depth gives a slope of 0.90, which confirms the expected trend. The
2.8. Required data resulting average soil depth of the study catchments is 1.29 m. The
spatial distribution of soil depth is shown in Fig. 8. Schiliro et al,
An overview of the input data for the model, the data sources, (2016) have used the soil depth model by Saulnier et al. (1997)
and the spatial resolution are provided in Table 1. Most of the which relates soil depth to slope, and depends on the maximum
required input parameters were estimated based on the digital and minimum soil depth as parameters. The results from these
elevation model (DEM), land use map and soil texture map calculations show similar trends as the results from equation (23).
(Table 2). The native spatial resolution of the national soil texture However, besides our higher spatial resolution, our statistical cor-
map of Sicily is coarse (Lombardo et al., 2018b). The soil texture relation uses a higher number of variables to predict soil depth and
map was therefore validated and adapted with field measurements, can therefore account for a greater variability.
collected by Costantini et al. (2014) and Schillaci et al. (2017). Initial
conditions for soil moisture content were based on the hydrological
2.10. Simulations and calibration method
modeling by Aronica et al. (2012). The landslide inventory was
interpreted from high resolution images, field checked and ras-
To compare our modeling method to existing similar method-
terized. Rainfall data, with a temporal resolution of 10 min, was
ologies, several simulations were performed (Table 4). Firstly, the
collected for two nearby station, S. Stefano and Fiumedinisi.
full method was used, simulating slope failure locations, timing and
We calculated both the leaf area index and vegetation cover
runout, together with catchment hydrology. Secondly, the landslide
based on the Normalized Differential Vegetation Index (NDVI). For
inventory was used to provide a predefined failure volume. The
vegetation cover, we linearly scale between minimum and
highest 20 percent of the inventory impacted areas was taken as
maximum NDVI value (Choudhury et al., 1994). For the leaf area
initiation area. The depth of failure was varied to be 50, 70 and 90
index, the vegetation cover is used in an empirical relationship
percent of the soil depth. This volume is introduced in the middle of
(Equation (22)) (Choudhury, 1987).
the precipitation event. Finally, both these simulations were
lnð1  vegetationcoverÞ repeated without the addition of the hydrology component. When
LAI ¼ (23) hydrology is not simulated, rainfall is neglected and soil moisture
9:1
contents are equal to their initial value. This simplifies the model to
With LAI the leaf area index (). a more traditional catchment model. However, without additional
infiltration, slopes can not become unstable. Therefore, we alter the
2.9. Spatial soil depth estimation initial moisture content to an estimation of soil moisture after 50
percent of the rainfall event has passed. This causes the slope
Unfortunately there was no spatial data concerning soil depth failure to be initiated immediately.
available within the extent of the study area, except for a limited From satellite images, reports and photos it can be observed that
number of field measurements, and estimations based on the a debris flow fan was formed during the event (Fig. 9).
interpretation of oblique photographs of landslide scars after the Finally, in order to validate the transportation of debris flow
event. Soil depth forms one of the most vital parameters in the case solids, a simulation was performed where coastal water is included
of shallow landslides, and a spatial estimate had to be made in the model. This allowed for the formation of the debris fan. A
(Kuriakose et al., 2009a,b). Spatial soil depth was estimated using section of 100 m of linearly decreasing elevation was added along
the method described by Kuriakose et al. (2009a,b). In this method, the coast line. Slope was estimated from the coastal elevation
variables such as slope, profile curvature, distance form channels, published by Casalbore et al. (2011). At the computational bound-
wetness index and distance from the coast are statistically related aries, forced water volumes were set to simulate the sea level.
to the measured soil depth. This statistical model was used to Initial velocities were set to zero, since no information was available
provide an estimation of soil depth in the entire catchment. For our considering sea currents during the event.

Table 1
The data used in the model, their sources and spatial resolutions.

Input Data Method Spatial resolution Source

Elevation Model Lidar DEM (0.22 m vertical-accuracy) 2m A.R.T.A (2008)

Land Use Map Sentinel-2 Supervised Classification 10 m Copernicus Sentinel Data, 2015
Soil Texture Map Country-wide Soil Texture Product 200 m Fierotti and Fierotti, 1988
NDVI Spot-6 Satellite Product 4m KompSat 2 Data, 2010
8 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

Table 2
Maps that are derived from the input data maps, and how they were validated.

Elevation Model Soil Texture Map Land Use Map NDVI

Slope Saturated Conductivityb Mannings N Vegetation Covera

Soil Depthb Porosity Surface roughness Leaf Area Indexa
Channel Locationa Soil Cohesion Vegetation Heighta
Channel Deptha Densityb Root Cohesion
a
Validated with areal and site images.
b
Validated with field measurements.

Table 3
p0  pe
The constants derived from the statistical
k¼ (25)
correlation of soil depth to topographical 1  pe
parameters.

intersect 2.37 Where p0 is the total accuracy () and pe is the probability of un-
correlated agreement ().
a 0.01
b 0.12 The probability of uncorrelated agreement indicated the chance
c 1.03 that simulated and predicted maps agreed due to chance. For a
d 2.73 binary system with two observers, pe can be calculated with
equation (24)

All simulations were calibrated to the runout patterns in the Nagr N1;pos *N2;pos N1;neg *N2;neg
p0 ¼ pe ¼ * (26)
landslide inventory. Calibration was performed manually by Ntotal Ntotal Ntotal
altering important input parameters. Based on theoretical consid-
erations, the soil cohesion, internal friction angle, infiltration rate Where Ntotal is the total number of measurements (), Nagr is the
and initial moisture were varied by multiplying these input maps number of agreed measurements and Ni;pos and Ni;neg are the
with a scalar value. Accuracy was calculated using Cohen's Kappa numbers of measurements where observer i is respectively posi-
(Equation (25)). This measure of model efficiency takes into ac- tive or negative ().
count both correctly predicted positives, negatives, and incorrectly Since our inventory of the event consists only of full shallow
predicted areas. landslide/debris flow impact areas and does not distinguish be-
tween source and deposition areas, the validation of failures is
performed using the 20 percent highest elevations for each

Fig. 8. The spatial distribution of predicted soil depth, and the correlation between predicted and estimated soil depth.

Table 4
Performed simulations.

Simulation Includes Hydrology Failure from Inventory Fraction of soil used in failure Coast

A1 No Yes 0.5 No
A2 No Yes 0.7 No
A3 No Yes 0.9 No
B1 Yes Yes 0.5 No
B2 Yes Yes 0.7 No
B3 Yes Yes 0.9 No
C Yes No e No
D Yes No e Yes
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 9

3.1. Simulation results with inventory-based slope failure

To further investigate the influence of added hydrology on the

runout, we carried out another simulation in which we did not
simulate the slope failures themselves in the model. Instead, we
took the slope failures from the landslide inventory and used the
volumes of these as the starting points for the debris flows. A
comparison of debris flow runout patterns and the landslide in-
ventory are shown in Fig. 12. Here, initiation of shallow landslides is
Fig. 9. Images of the Scaletta catchment before and after the 1-10-2009 event (left).
simultaneous, and based on the inventory instead of physically
Photograph of the Giampillieri catchment outlet one day after the event (right). predicted. Furthermore, Fig. 12 compares maximum debris flow
height for inventory-based landslides with and without hydrology.
The accuracies for these simulations are 0.64, 0.60, 0.60 and 0.56,
landslide impact polygon. An overlap of a failure area might not 0.58, 0.60 for simulations with and without hydrology respectively,
overlap with the actual source area of the inventory, but instead with a failure fraction of 0.5, 0.7 and 0.9, respectively. A full over-
with part of the runout. For completeness with give the correlation view of accuracies is provided in Table 6.
between failures and the inventory as a mere indication of corre-
lation between the two.
3.2. Simulations with coastal deposition

The results of the model simulations that included a coastal

3. Results section are shown in Fig. 13A. A visual comparison of these results
with the satellite images after the event are shown in Fig. 13B.
First we analyzed the runout of shallow landslides and debris
flows with the inventory. The resulting slope failure patterns and 4. Discussion
debris flow runout patterns were compared with the inventory.
However, the comparison was complicated due to several aspect. 4.1. Slope failure
The landslides in the inventory were only mapped to the point
where they reached a channel. In reality, however, the runout often The resulting slope failure patterns show a fairly good match
extended much further, and landslide materials were incorporated with the highest 20 percent of elevation for each landslide impact
in debris flows. Since the model counts any location, including a polygon. While this does not directly confirm correlation with
channel, with high solids content as a debris flow, the channel re- actual landslide source areas, it does provide confidence in the
gions have been masked to be excluded from the accuracy patterns predicted by the method. Individual landslides are in
calculation. numerous cases predicted with high accuracy. Moreover, general
The simulated slope failure depth, and a comparison between patterns and areas of high hazard are well estimated. In terms of
simulated slope failure and the landslide inventory are shown in the number of landslides, 341 out of the 395 landslides overlap
Fig. 10 respectively. A comparison of the failures with the estimated with the modelled slope failures. However, the modelled slope
failure zones shows an accuracy of 95.9 percent, with a kappa value failures are over-estimated with 591 predicted landslide locations.
of 0.27. This is an over-estimation of actual accuracy since failures But many of these may also be represented by numerous small
can overlap with the runout regions. failure locations located on the edges of the channels, which were
Fig. 11 shows maximum debris flow height during simulations not mapped in the landslide inventory. One of the possible reasons
with and without hydrology. Here, both slope failure and runout for this could be that these were subsequently eroded by the
are simulated. Furthermore, the comparison between simulated intense debris flood that also left deposits on the channel sides,
runout patterns and the landslide inventory is shown. Table 5 making it harder to map small slope failure events. Since the over-
shows the accuracy of the predicted debris flow runout and slope predicted landslides were small in volume, they added relatively
failures. The accuracy of the simulated runout patterns are 89.9 and small volumes to the total runout. Thus, while predictions for single
91.8 percent with kappa coefficients of 0.211 and 0.216 with and landslides can still increase in accuracy, general patterns of slope
without hydrology, respectively. instability in the catchment are well predicted.

Fig. 10. A comparison of simulated slope failure with the landslide inventory for the Scaletta catchment (Left). The simulated failure depth (Right).
10 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

Fig. 11. Maximum debris flow depth, both with hydrology (left) and without hydrology (right), and a comparison of debris flow runout with the landslide inventory, with channels
removed (bottom).

Table 5
Comparison between model results and the landslide inventory (TP ¼ Percentage True Positive, TN ¼ True negatives (%), FP ¼ Fasle positives (%), FN ¼ False negatives (%)).

TP TN FP FN Cohens Kappa Accuracy Predicted Nr Landslides Missed Nr

(%) (%) (%) (%) (%) Landslides

Landslide Inventory 4.9 95.1 0 0 1 100 395 0

Debris Flow Runout
With Hydrology 2.0 87.9 7.0 3.2 0.224 89.9 496 42
Without Hydrology 1.6 89.2 5.6 3.4 0.216 91.8 403 81

Many factors might be the cause of inaccuracies in predicting where sub-surface structural details are important. These details,
landslide locations. Primarily, the digital elevation model and es- often formed by the morphological history of an area, can deter-
timations of soil depth appear to influence the location of slope mine spatial patterns in slope failure. In the case of the Scaletta
failure in the performed simulations. In the case of relatively small area, knowledge on subsurface structure was not available, limiting
shallow landslides, high spatial detail is required for accurate slope the accuracy of the model.
failure estimations. For both the calculation of slope, and in esti- The method that was implemented to estimate the failure
mating soil depth, the quality of the elevation model determines fraction of slope failure provided values between 0.1 and 1.7 m.
the level of detail. In our study case, limited knowledge was Data on the actual depth of slope failure in the Scaletta catchment is
available on the spatial distribution of soil depth. Soil depth was not only available from photographs. Based on field photographs taken
varying for the various types of soils (e.g. colluvial, alluvial, residual after the event, it can be concluded that the slope failure depth
soils). This could be responsible for most of the incorrect estimation provides reasonable estimations of the depth of the failure plane
of slope failure. Furthermore, landslides are complex processes for the modelled landslides. The average failure depth within the
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 11

Fig. 12. Debris flow runout and a masked comparison of debris flow runout with the landslide inventory.

Table 6 visits several days after the event. Within the iterative slope failure
Cohens Kappa values for all performed calibrated simulations. (TP ¼ True Positive model, the slope failure depth is significantly influenced by cali-
(%), TN ¼ True negatives (%), FP ¼ Fasle positives (%), FN ¼ False negatives (%)).
bration for the stability margin parameter. While over-estimating
Inventory TP (%) TN FP FN Accuracy Cohens Kappa failure volumes caused an over-estimation of runout, decreasing
(%) (%) (%) (%) viscosity could, to a limited extent, compensate this. Due to the
4.9 95.1 0 0 100 1 nature of the model, and the large number of free parameters that
Debris Flow Impact Area With Hydrology determine the outcome, there is an unavoidable equifinality.
Failure Fraction ¼ 0.5 3.6 92.5 2.5 1.4 96.1 0.638 Because of this, runout patterns alone are not sufficient evidence in
Failure Fraction ¼ 0.7 3.3 92.5 2.4 1.8 95.8 0.597 determining accurate failure depths.
Failure Fraction ¼ 0.9 4.3 90.4 4.4 0.9 94.8 0.596
In recent years, several models have attempted to estimate both
Debris Flow Impact Area Without Hydrology
Failure Fraction ¼ 0.5 2.4 93.9 2.8 0.9 96.3 0.559
slope failure and failure volumes. Mainly, approaches are based on
Failure Fraction ¼ 0.7 3.0 93.0 2.2 1.8 96.0 0.583 random ellipsoid sampling of slopes. With this method, each slope
Failure Fraction ¼ 0.9 3.4 92.3 1.8 3.5 95.7 0.599 segment is approached by a large number of randomly shaped
elipsoids. The soil within this ellipsoid is then taken as a possible
failure volume, and stability is calculated using a limit equilibrium
Giampillieri region, which is the neighboring catchment with method. SCOOP3D (Reid et al., 2015) and the modeling method by
identical characteristics, is estimated to be near 0.75 based on field Mergili et al. (2014a,b), implement this method. This method

Fig. 13. A) Final solid and fluid height with the simulation of a part of the coast. Simulation uses predicted slope failure and debris flow runout.
12 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

provides a highly detailed simulation of subsurface stability in a estimation of debris flow susceptibility. Based on the data that is
region with the best accuracy in the field. However, since the iter- available for the catchment, a perfect prediction of slope failure and
ative method used in this research does not requires the sampling debris flows can furthermore not be expected. To further validation
of large amounts of random ellipsoids, it can be performed for each of the debris flood behavior, the maximum flow depth can be
hydrological calculation step. Compared to random ellipsoid compared to photographs from the days after the event. In the
methods, the iterative slope failure method allows for a quick calibrated simulation, the maximum flow height is equal to 2.8 m.
assessment of slope stability and possible failure volume. Thus, When compared with field photos, this shows the realistic esti-
while the assumptions in the iterative slope failure method are mation of debris flood properties by the model.
more ambiguous, calculation time strongly improves applicability. Simulation results for catchment-scale debris flow runout with
In another type of model, infiltration simulations are included in inventory-based initiation show a good correlation with the land-
slope stability models to predict event-based slope failure on a slide inventory. The values for Cohens Kappa reach 0.638 with
catchment scale. Examples of this are iCRESTRIGRS (Zhang et al., simulated hydrology and a failure fraction of 0.5. One of the causes
2016) and CRESLIDE (He et al., 2016). In both these models, the of the reached accuracy is the lack of individual calibration. Due to
methods that are used to estimate slope stability are highly similar. the complex behavior of debris flows, calibration is required to
Slope stability equations are coupled with a dynamic soil moisture achieve sufficient correlation. In the case of individual debris flows,
simulation. Accuracy results are therefore similar. The accuracy of parameters such as viscosity, yield stress and surface flow resis-
the proposed iterative slope failure method produces similar ac- tance are altered in order to match simulated behavior with mea-
curacy in estimating slope failures, when compared to other surements. In the case of catchment-scale estimations of debris
methods. All current models are forced to make similar assump- flow runout, individual calibration is too time-consuming, and re-
tions considering sub-surface structure due to a lack of available duces predictive power. However, two important properties of the
data. However, in the case of iterative slope failure, failure volumes used methodologies compensated the lack of calibration. Due to the
are provided without the increase in computation time required scalability of the two-phase debris flow equations by Pudasaini
with random ellipsoid sampling. (2012), estimations of rheological parameters was automatically
performed. Furthermore, because of the soil data provided by
4.2. Runout Schillaci et al. (2017), the automatic estimations increased in ac-
curacy. Simulation for catchment-scale debris flow runout with
Debris flow simulations were performed for two types of situ- inventory-based initiation were carried out with varying failure
ations, with predicted initiation and inventory-based initiation. It is depth. From the achieved accuracies it can be seen that the addition
important to note that the experts only mapped the landslides from of hydrology significantly improved the simulation results. Without
the 1-10-2009 event up to the channel intersection, without hydrology, a failure depth of 90 percent of the soil depth provided
including the full runout. Because of this, debris flow runout was the highest kappa value. In reality slope failures were shallow, and
masked out in channel locations. Based on reports, photographs the simulation with hydrology, which was most accurate at 50
and the resulting debris flow deposits in the area, it is known that percent of the soil depth, was more consistent with the seen slope
debris flows took place, and a substantial amount of material was failures.
transported downstream towards the coast. However, debris flow The simulated output for both inventory-based and prediction-
runout in channels is completely dependent on the channel loca- based debris flow runout was highly influenced by the incorpora-
tion. Including these areas in the runout comparison would artifi- tion of a hydrological simulation. In both cases, the accuracy of
cially increase performance of the runout simulations without simulations was significantly higher due to the addition of simu-
indicating improved performance. Another important factor in the lated hydrology. As it is shown in Fig. 11, several aspects of the
accuracy of the model comes from the small slope failures directly simulation are altered. Firstly, the runout distance of solid materials
nest to the channels. These small side-channel slope failures are not is highly increased. In the case of simulations without hydrology,
included in the landslide inventory. Because of the debris flooding, runout distances often only reach near the closest channel. How-
mapping these landslides is a difficult processes after the event. ever, when a hydrological simulation is included, more material is
However, since no data is known considering these shallow slope taken through the channel toward the coast, where a fan of debris
failures, they are indicated as inaccurate predictions. However, material was created at the channel outlet. Secondly, the further
based on images and reports, numerous of these landslides runout of solid material decreases the depth of deposits that are left
occurred in the region (Goswami et al., 2011; Ardizzone et al., on the slopes below a slope failure. Finally, debris material has a
2012). significant influence on the behavior of channel flow and flooding
In the case of prediction-based runout simulations, debris flow in the catchment. Compared to a simulation with only hydrology,
initiation is determined by predicted slope failure. Therefore, cor- flood velocities are reduced by the solid material that is taken
relation with the landslide inventory is highly bound by the accu- through the channel. In some cases, debris material nearly stag-
racy of the slope failure predictions. The runout patterns of nates channel flow, or temporarily blocks it entirely. An example of
predicted slope failure show reasonable correlation with the velocity reduction by the solid volume in the channel is shown is
landslide inventory. General patterns in runout distance and debris Fig. 15.
flow density throughout the catchment are well predicted. Exam- The addition of hydrology altered the behavior of the debris flow
ples of the types of patterns that are visible in the map are provided simulation due to several effects. Firstly, viscosity lowers when
in Fig. 14. debris material is diluted with overland flow. Secondly, solid ma-
Firstly, areas where debris flows are accurately predicted are terials experience a higher drag force with a large volumes of water.
visible in Fig. 14a. In such patterns, that form the predominant area Thus, due to the scalability and flexibility of the used two-phase
of the catchment, both runout and initiation are simulated with debris flow equations, flow properties are continuously influ-
substantial accuracy. Secondly, Fig. 14b shows areas where debris enced by the flow composition. Because of this, a variety of emer-
flow runout is overestimated due to small slope failures on the gent behaviors arise in the model.
channel sides. Thirdly, several areas of the catchment are similar to Simulations of debris flow runout with hydrological modeling
the patterns in Fig. 14c. Here debris flow correlation is low. How- show the creation of a fan of debris material at the outlet of one of
ever, general patterns are highly similar, providing a good the sub-catchments in the Scaletta region. A similar fan of material
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 13

Fig. 14. An overview of several patterns in the predicted slope failure and debris flow runout. A) mostly correct B) over-estimates small slope failure along channel sides C) initiation
along the same stream, but not correctly placed.

or debris flows would have been different when these hazardous

processes where approached separately.
In the case of the Scaletta catchment, downstream behavior of
the hazardous processes has changed due the combined approach
to debris flows and flooding. A comparison of maximum fluid/solid
flow height for an integrated simulation and a flash flood simula-
tion, is shown in Fig. 16. In the case of the integrated simulation
with slope failures, the reduced infiltration, and runout of solid
materials and soil water strongly increases the total flow depth. A
Fig. 15. One of the locations in the catchment where a shallow landslide reduced flow primary change in the risk is then caused by the exposure to the
velocities in the channel, partly blocking flow. hazard. When bridges or other infrastructure are present, a further
change in risk is created. Debris flows generally provide more
impact force on blocking objects, such as bridges, and cause more
has been described in reports, and shows on satellite images taken destruction. Details such as these can, in scenarios considering
one year after the event. Based on the areal images and photo- hazardous event, make substantial difference concerning evacua-
graphs of similar fans created during the event, the total volume of tion, and rescue travel times.
deposited solids in the fan is estimated to be around 2460 cubic Finally, the hazard and risk are influenced by the presence of
meters. During the simulation that was calibrated on debris flow sediment in (flash) floods. Generally, flood simulations exclude
runout, a volume of 3490 cubic meters was deposited in the fan. sediment-related processes. The influence of flood deposits on the
The difference in the predicted and estimated volumes can be vulnerability of elements-at-risk, is therefore often ignored. How-
explained by the fact that the satellite image was taken one year ever, in the case of the Scaletta catchment, debris floods left sub-
after the event. Thus, erosion has most likely decreased the size of stantial amounts of deposits. Deposits often are required to be
the fan, and shifted it towards the most common direction of the manually removed, which can be costly and time-consuming. More
coastal current. The implementation of debris flow runout in a costs could be caused by these additional damage to crops, building
hydrological catchment model thus provided an acceptable and infrastructure. For an accurate representation of risk in the case
approximation of the processes that caused and lead to the creation of hazardous processes such as debris floods, the inclusion of
of a fan of debris deposits near the coast. sediment in flood vulnerability curves is required.
In the case of the presented research, calibration was required to
4.3. Influence on hazard obtain adequate results. In many applications of physical-based
modeling, particularly related to hazard and risk assessment, for-
Catchment scale risk assessment often incorporates simulations ward analysis is useful. The current setup for the Scaletta area is not
of debris flow runout. In the case of multi-hazard risk assessment, tested for such an application. Because of the physically-based
separate models are used to simulate debris flows and hydrological implementation of all related processes, it seems forward pre-
hazards. Since shallow landslides, debris flows and floods share a dictions should be possible provided there is sufficient input data of
common meteorological trigger, these hazardous processes often high enough quality. In our case, particularly soil information was
occur simultaneously. Similarly, in the Scaletta catchment, haz- of low certainty, which made the soil strength parameters the
ardous processes took place simultaneously. Therefore, integrating values that required the most calibration. With all currently avail-
debris flow simulations in a hydrological catchment model has able methods, forward predictions of complex multi-hazard events
significantly altered the behavior of the simulated hazardous pro- should become increasingly accurate with more available data.
cesses. As a result, the estimation of hazard or risk for flooding and/ Further research should investigate the possibility and accurate of
14 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

Fig. 16. A comparison of maximum flow height for the scenario's with (left) and without (right) predicted shallow slope failure.

forward modeling of complex multi-hazard events in an area with and anyone can participate. Source-code of the model that is pre-
high-density measurements of soil and hydrological information. sented in this paper, and older versions of OpenLISEM are therefore
available online (https://sourceforge.net/projects/lisem/).
5. Conclusions
Future research
Several aspects of the developed, integrated, modeling method
provide added value above separated approaches to debris flows, Before integrated modeling approaches could be applied in
shallow landslides and flood simulations. During the 2009 extreme hazard and risk analysis, several research directions might require
event in the Scaletta catchment, a wide variety of processes such as further investigation. Firstly, vulnerability curves are generally
slope failures, debris flooding, channel blocking and solids depo- investigated separately for flood and debris flow processes. More-
sition occurred. Our developed modeling approach allowed for over, sediment is often completely neglected in flood vulnera-
replication of the multi-hazard impact of the event. Application of a bility, while both debris flood and flood deposits provide substan-
physically based multi-hazard model allows for simulating this tial costs and possible hazard. The inclusion of sediment and the
broad spectrum of hazardous phenomena. When separated simu- combination of debris flow and flood vulnerability should, for some
lations are used in the analysis of event such as these, critical as- regions, be investigated. Secondly, the proposed modeling
pects of the event are not simulated. This affects both the potential approach was tested predominantly with data from remote sensing
for investigative use and predictive power. Furthermore, resulting sources. While this allows for usage in data-poor regions, the model
risk and hazard could be seriously under-predicted when debris performance depends largely on the quality of the digital elevation
flows and flooding, which often take place simultaneously, are model. The effect however, of digital elevation models of separate
separately analyzed. Using a model that integrates (flash) flood sources, is still to be investigated.
behavior, shallow landslides, and debris flow runout, can thus in-
crease the accuracy of hazard and risk assessment. Acknowledgements
The purpose of the developed model in this case study is not the
accurate prediction of the temporal and spatial location of slope We wish to express our gratitude to the editorial team of the
failure. Lack of available data limits the predictions of slope failure journal and the reviewers for the thorough efforts to help improve
to mere estimations and probabilities. Furthermore, the developed this research. The presented work is part of the Earth System
model does not provide a full description of all relevant processes. Analysis group, ITC, Twente University, the Netherlands in corpo-
OpenLISEM is event-based in nature, and therefore neglects long- ration with Computational Earthquake Seismology and Extreme
term physical processes, such as lateral ground water flow and Statistic groups, King Abdullah University of Science and Technol-
evaporation, which are relevant in slope stability modeling on a ogy. Satellite images were provided by the European Space Agency
longer time-scale. It is rather meant for physically based simulation through the project titled: “A remote sensing based approach for
of behavior of slope failure, debris flow runout and interactions storm triggered debris flow hazard modeling: application in Med-
with catchment hydrology. Detailed local knowledge of landslide iterranean and tropical Pacific areas”, code: C1P.14151, PI: Luigi
risk requires both field visits by experts, and collecting subsurface Lombardo, ID: 14151. Funded by the State Key Laboratory of Geo-
structure data. While OpenLISEM is not meant for long-term slope hazard Prevention and Geoenvironment Protection Open Fund
prediction modeling, it can be used for several other important SKLGP2018K001.
applications. Applications for the Investigation into the physical
processes that lead to and cause hazardous processes in a catch-
References
ment is the first and primary usage. A second important use is the
simulation of scenarios in risk and hazard assessment. Finally, Adegbe, M., Alkema, D., Jetten, V.G., Agbor, A.T., Abdullahi, I.N., Shehu, O.U.,
investigation into the effectiveness of measures, and the assess- Unubi, A.S., 2013. Post seismic debris flow modeling using Flo-2D; case study of
ment of downstream hazards are an important functionality of Yingxiu, Sichuan Pronvince, China. J. Geogr. Geol. 5 (3), 101.
Ardizzone, F., Basile, G., Cardinali, M., Casagli, N., Del Conte, S., Del Ventisette, C.,
OpenLISEM. et al., 2012. Landslide inventory map for the Briga and the giampilieri catch-
ments, NE Sicily, Italy. J. Maps 8 (2), 176e180.
Aronica, G.T., Brigandí, G., Morey, N., 2012. Flash floods and debris flow in the city
Software availability
area of Messina, north-east part of Sicily, Italy in October 2009: the case of the
Giampilieri catchment. Nat. Hazards Earth Syst. Sci. 12 (5), 1295.
OpenLISEM is a freely available model that features a full A.R.T.A, (2008), http://www.sitr.regione.sicilia.it/geoportale/it/metadata/details/
interface and documentation (https://blog.utwente.nl/lisem/). It 502 (retrieved 10-04-2017).
Baartman, J.E., Jetten, V.G., Ritsema, C.J., Vente, J., 2012. Exploring effects of rainfall
furthermore provides tools for erosion modeling and (flash) flood intensity and duration on soil erosion at the catchment scale using openLISEM:
predictions. Development of OpenLISEM is an open-source process, Prado catchment, SE Spain. Hydrol. Process. 26 (7), 1034e1049.
 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16 15

Bartelt, P., Buehler, Y., Christen, M., Deubelbeiss, Y., Graf, C., McArdell, B.W., 2013. the 2009 storm event in Messina (Sicily, southern Italy). Nat. Hazards 79 (3),
RAMMSerapid Mass Movement Simulation, a Modeling System for Debris 1621e1648.
Flows in Research and Practice, User Manual V1. 5, Debris Flow. WSL Institute Lombardo, L., Fubelli, G., Amato, G., Bonasera, M., 2016. Presence-only approach to
for Snow and Avalanche Research SLF. assess landslide triggering-thickness susceptibility: a test for the Mili catch-
Bastiaanssen, W.G.M., Huygen, J., Schakel, J.K., Van Den Broek, B.J., 1996. Modeling ment (north-eastern Sicily, Italy). Nat. Hazards 84 (1), 565e588.
the Soil-water-crop-atmosphere System to Improve Agricultural Water Man- Lombardo, L., Opitz, T., Huser, R., 2018a. Point process-based modeling of multiple
agement in Arid Zones (SWATRE). debris flow landslides using INLA: an application to the 2009 Messina disaster.
Beguería, S., Van Asch, T.W., Malet, J.P., Gro € ndahl, S., 2009. A GIS-based numerical Stochastic Environmental Research and Risk Assessment To appear. https://
model for simulating the kinematics of mud and debris flows over complex doi.org/10.1007/s00477-018-1518-0.
terrain. Nat. Hazards Earth Syst. Sci. 9 (6), 1897e1909. Lombardo, L., Saia, S., Schillaci, C., Mai, P.M., Huser, R., 2018b. Modeling soil organic
Bout, B., Jetten, V.G., 2018. The validity of flow approximations when simulating carbon with Quantile Regression: dissecting predictors' effects on carbon
catchment-integrated flash floods. J. Hydrol. 556, 674e688. stocks. Geoderma 318, 148e159.
Cama, M., Lombardo, L., Conoscenti, C., Agnesi, V., Rotigliano, E., 2015. Predicting Luna, B.Q., Blahut, J., van Asch, T., van Westen, C., Kappes, M., 2015, April. AschFlow-
storm-triggered debris flow events: application to the 2009 Ionian Peloritan A dynamic landslide run-out model for medium scale hazard analysis. In: EGU
disaster (Sicily, Italy). Nat. Hazards Earth Syst. Sci. 15 (8), 1785e1806. General Assembly Conference Abstracts, 17, p. 13656.
Cama, M., Lombardo, L., Conoscenti, C., Rotigliano, E., 2017. Improving transferability Luna, B.Q., Blahut, J., Camera, C., van Westen, C., Apuani, T., Jetten, V., Sterlacchini, S.,
strategies for debris flow susceptibility assessment: application to the Saponara 2014. Physically based dynamic run-out modeling for quantitative debris flow
and Itala catchments (Messina, Italy). Geomorphology 288, 52e65. risk assessment: a case study in Tresenda, northern Italy. Environ. Earth Sci. 72
Casalbore, D., Chiocci, F.L., Mugnozza, G.S., Tommasi, P., Sposato, A., 2011. Flash- (3), 645e661.
flood hyperpycnal flows generating shallow-water landslides at Fiumara Lupiano, V., Machado, G.E., Crisci, G.M., Di Gregorio, S., 2016. Simulations of debris/
mouths in Western Messina Strait (Italy). Mar. Geophys. Res. 32 (1e2), 257. mud flows invading urban areas: a cellular automata approach with SCIDDICA.
Choudhury, B.J., 1987. Relationships between vegetation indices, radiation absorp- Lect. Notes Comput. Sci. 9863 LNCS, 291e302.
tion, and net photosynthesis evaluated by a sensitivity analysis. Rem. Sens. Mergili, M., Fellin, W., Moreiras, S.M., Sto € tter, J., 2011. Simulation of debris flows in
Environ. 22 (2), 209e233. the Central Andes based on Open Source GIS: possibilities, limitations, and
Choudhury, B.J., Ahmed, N.U., Idso, S.B., Reginato, R.J., Daughtry, C.S., 1994. Relations parameter sensitivity. Nat. Hazards 61 (3), 1051e1081.
between evaporation coefficients and vegetation indices studied by model Mergili, M., Marchesini, I., Alvioli, M., Metz, M., Schneider-Muntau, B., Rossi, M.,
simulations. Rem. Sens. Environ. 50 (1), 1e17. Guzzetti, F., 2014a. A strategy for GIS-based 3-D slope stability modeling over
Chow, V. Te, 1959. Open Channel Hydraulics. McGraw-Hill Book Company, Inc, New large areas. Geosci. Model Dev. (GMD) 7 (6), 2969e2982.
York. Mergili, M., Marchesini, I., Rossi, M., Guzzetti, F., Fellin, W., 2014b. Spatially
Costantini, E.A.C., Barbetti, R., Fantappie , M., L'Abate, G., Lorenzetti, R., Napoli, R., distributed three-dimensional slope stability modeling in a raster GIS. Geo-
Marchetti, A., Rivieccio, R., 2014. The soil map of Italy: a hierarchy of geo- morphology 206, 178e195.
databases, from soil regions to sub-systems. In: GlobalSoilMap: Basis of the Morgenstern, N.R., Price, V.E., 1965. The Analysis of the Stability of General Slip
Global Spatial Soil Information System - Proceedings of the 1st GlobalSoilMap Surfaces.
Conference, pp. 109e112. Nguyen, H.T., Wiatr, T., Ferna ndez-Steeger, T.M., Reicherter, K., Rodrigues, D.M.,
D'Ambrosio, D., Iovine, G., Spataro, W., Miyamoto, H., 2007. A macroscopic colli- Azzam, R., 2013. Landslide hazard and cascading effects following the extreme
sional model for debris-flows simulation. Environ. Model. Software 22 (10), rainfall event on Madeira Island (February 2010). Nat. Hazards 65 (1), 635e652.
1417e1436. Nikolopoulos, E.I., Borga, M., Creutin, J.D., Marra, F., 2015. Estimation of debris flow
Fan, L., Lehmann, P., McArdell, B., Or, D., 2017. Linking rainfall-induced landslides triggering rainfall: influence of rain gauge density and interpolation methods.
with debris flows runout patterns towards catchment scale hazard assessment. Geomorphology 243, 40e50.
Geomorphology 280, 1e15. O'Brien, J.S., 2006. FLO-2D User's Manual, Version 2006.01. FLO-2D Software. Inc.,
Fellenius, W., 1936, July. Calculation of the stability of earth dams. In: Transactions Nutrioso.
of the 2nd Congress on Large Dams, Washington, DC, vol. 4. International O'brien, J.S., Julien, P.Y., Fullerton, W.T., 1993. Two-dimensional water flood and
Commission on Large Dams (ICOLD), Paris, pp. 445e463. mudflow simulation. J. Hydraul. Eng. 119 (2), 244e261.
Fierotti, G., Fierotti, G., 1988. Carta dei suoli della Sicilia (No. 631.47458). Regione Pudasaini, S.P., 2012. A general two-phase debris flow model. J. Geophys. Res.: Earth
siciliana, Assessorato territorio ed ambiente. Surface 117 (F3).
Formetta, G., Antonello, A., Franceschi, S., David, O., Rigon, R., 2014. Hydrological Reid, M.E., Christian, S.B., Brien, D.L., Henderson, S.T., 2015. Scoops3DdSoftware to
modeling with components: a GIS-based open-source framework. Environ. Analyze Three-Dimensional Slope Stability Throughout a Digital Landscape. In
Model. Software 55, 190e200. Tech. Rep. US Geological Survey Techniques and Methods, p. 218. book 14.
Goswami, R., Mitchell, N.C., Brocklehurst, S.H., 2011. Distribution and causes of Rigon, R., Bertoldi, G., Over, T.M., 2006. GEOtop: a distributed hydrological model
landslides in the eastern peloritani of NE sicily and western aspromonte of SW with coupled water and energy budgets. J. Hydrometeorol. 7 (3), 371e388.
Calabria, Italy. Geomorphology 132 (3), 111e122. Saulnier, G.M., Beven, K., Obled, C., 1997. Including spatially variable effective soil
Green, W.H., Ampt, G.A., 1911. Studies in soil physics. J. Agric. Sci. 4, 1e24. depths in TOPMODEL. J. Hydrol. 202 (1), 158e172.
He, X., Hong, Y., Vergara, H., Zhang, K., Kirstetter, P.E., Gourley, J.J., et al., 2016. Scheidl, C., Rickenmann, D., McArdell, B.W., 2013. Runout Prediction of Debris Flows
Development of a coupled hydrological-geotechnical framework for rainfall- and Similar Mass Movements. In Landslide Science and Practice. Springer Berlin
induced landslides prediction. J. Hydrol. 543, 395e405. Heidelberg, pp. 221e229.
Hengl, T., de Jesus, J.M., Heuvelink, G.B., Gonzalez, M.R., Kilibarda, M., Blagoti c, A., , L., Montrasio, L., Scarascia Mugnozza, G., 2016. Prediction of shallow
Schiliro
et al., 2017. SoilGrids250m: global gridded soil information based on machine landslide occurrence: validation of a physically-based approach through a real
learning. PLoS One 12 (2), e0169748. case study. Sci. Total Environ. 569e570, 134e144.
Horton, P., Jaboyedoff, M., Rudaz, B., Zimmermann, M., 2013. Flow-R, a model for Schillaci, C., Lombardo, L., Saia, S., Fantappie , M., Ma €rker, M., Acutis, M., 2017.
susceptibility mapping of debris flows and other gravitational hazards at a Modeling the topsoil carbon stock of agricultural lands with the Stochastic
regional scale. Nat. Hazards Earth Syst. Sci. 13 (4), 869. Gradient Treeboost in a semi-arid Mediterranean region. Geoderma 286,
Huang, Y., Zhang, W., Xu, Q., Xie, P., Hao, L., 2012. Run-out analysis of flow-like 35e45.
landslides triggered by the Ms 8.0 2008 Wenchuan earthquake using Smith, R.E., Parlange, J.Y., 1978. A parameter-efficient hydrologic infiltration model.
smoothed particle hydrodynamics. Landslides 9 (2), 275e283. Water Resour. Res. 14 (3), 533e538.
Janbu, N., 1975, April. Slope stability computations, 1973. In: Hirschfeld, R.C., Somma, R., 2006. The south-western side of the Calabrian Arc (Peloritani Moun-
Poulos, S.J. (Eds.), Embankment-dam Engineering. Textbook. JOHN WILEY AND tains): geological, structural and AMS evidence for passive clockwise rotations.
SONS INC., PUB., NY, p. 40. In International Journal of Rock Mechanics and J. Geodyn. 41 (4), 422e439.
Mining Sciences & Geomechanics Abstracts (Vol. 12, No. 4, p. 67). Pergamon. Stancanelli, L.M., Rosatti, G., Begnudelli, L., Armanini, A., Foti, E., 2013. Single or
Jetten, V.G., de Roo, A.P., 2001. Spatial Analysis of Erosion Conservation Measures Two-phase Modeling of Debris-flow? a Systematic Comparison of the Two
with LISEM. In Landscape Erosion and Evolution Modeling. Springer US, Approaches Applied to a Real Debris Flow in Giampilieri Village (Italy). In
pp. 429e445. Landslide Science and Practice. Springer Berlin Heidelberg, pp. 277e283.
Kuriakose, S.L., Van Beek, L.P.H., Van Westen, C.J., 2009a. Parameterizing a physically Tang, C., van Asch, T.W., Chang, M., Chen, G.Q., Zhao, X.H., Huang, X.C., 2012.
based shallow landslide model in a data poor region. Earth Surf. Process. Catastrophic debris flows on 13 August 2010 in the Qingping area, south-
Landforms 34 (6), 867e881. western China: the combined effects of a strong earthquake and subsequent
Kuriakose, S.L., Devkota, S., Rossiter, D.G., Jetten, V.G., 2009b. Prediction of soil rainstorms. Geomorphology 139, 559e576.
depth using environmental variables in an anthropogenic landscape, a case Tarolli, P., 2014. High-resolution topography for understanding Earth surface pro-
study in the Western Ghats of Kerala, India. Catena 79 (1), 27e38. cesses: opportunities and challenges. Geomorphology 216, 295e312.
Kwan, J.S., Sun, H.W., 2006. An improved landslide mobility model. Can. Geotech. J. Trigila, A., Iadanza, C., Esposito, C., Scarascia-Mugnozza, G., 2015. Comparison of
43 (5), 531e539. logistic regression and random forests techniques for shallow landslide sus-
Lombardo, L., Cama, M., Maerker, M., Rotigliano, E., 2014. A test of transferability for ceptibility assessment in giampilieri (NE sicily, Italy). Geomorphology 249,
landslides susceptibility models under extreme climatic events: application to 119e136.
the Messina 2009 disaster. Nat. Hazards 74 (3), 1951e1989. Van Asch, T.W., Buma, J., Van Beek, L.P.H., 1999. A view on some hydrological trig-
Lombardo, L., Cama, M., Conoscenti, C., Ma €rker, M., Rotigliano, E., 2015. Binary lo- gering systems in landslides. Geomorphology 30 (1), 25e32.
gistic regression versus stochastic gradient boosted decision trees in assessing van Beek, L.P.H., 2003. Assessment of the Influence of Changes in Land Use and
landslide susceptibility for multiple-occurring landslide events: application to Climate on Landslide Activity in a Mediterranean Environment (Doctoral
16 B. Bout et al. / Environmental Modelling & Software 105 (2018) 1e16

dissertation). (3e4), 112e131.

Van Beek, L.P.H., Van Asch, T.W., 2004. Regional assessment of the effects of land- Xie, M., Esaki, T., Cai, M., 2004. A GIS-based method for locating the critical 3D slip
use change on landslide hazard by means of physically based modelling. Nat. surface in a slope. Comput. Geotech. 31 (4), 267e277.
Hazards 31 (1), 289e304. Zanke, U., 1977. Berechnung der sinkgeschwindigkeiten von sedimenten. EV.
Van Westen, C.J., Van Asch, T.W., Soeters, R., 2006. Landslide hazard and risk Zhang, K., Xue, X., Yang, H., Gourley, J.J., Lu, N., Wan, Z., et al., 2016. iCRESTRIGRS: a
zonationdwhy is it still so difficult? Bull. Eng. Geol. Environ. 65 (2), 167e184. coupled modeling system for cascading flood-landslide disaster forecasting.
Van Westen, C.J., Castellanos, E., Kuriakose, S.L., 2008. Spatial data for landslide Hydrol. Earth Syst. Sci. 20 (12), 5035.
susceptibility, hazard, and vulnerability assessment: an overview. Eng. Geol. 102