Combining Spatial Models for Shallow Landslides and Debris-Flows Prediction

Abstract: Mass movements in Brazil are common phenomena, especially during strong rainfall events that occur frequently in the summer season. These phenomena cause losses of lives and serious damage to roads, bridges, and properties. Moreover, the illegal occupation by slums on the slopes around the cities intensifies the effect of the mass movement. This study aimed to develop a methodology that combines models of shallow landslides and debris-flows in order to create a map with landslides initiation and debris-flows volume and runout distance. The study area comprised of two catchments in Rio de Janeiro city: Quitite and Papagaio that drained side by side the west flank of the Macico da Tijuca, with an area of 5 km 2 . The method included the following steps: (a) location of the susceptible areas to landslides using SHALSTAB model; (b) determination of rheological parameters of debris-flow from the back-analysis technique; and (c) combination of SHALSTAB and FLO-2D models to delineate the areas more susceptible to mass movements. These scenarios were compared with the landslide and debris-flow event of February 1996. Many FLO-2D simulations were exhaustively


Introduction
Natural disasters constitute one of the largest socio-economic problems worldwide.Among these, earthquakes, floods and mass movements are the most striking events.According to the United Nations [1], mass movements cause major financial damage and death.In cities, they may have catastrophic effects due to the various anthropogenic changes imposed on the natural landscape [2,3].
Weather conditions and geomorphological characteristics of the Brazilian coast facilitate the occurrence of mass movements.Apart from these natural conditions, anthropogenic activity along spillways also triggers numerous mass movement events.The combination of shallow landslides with debris-flows caused catastrophic events in Rio de Janeiro, particularly after the 1960's when the low-income communities began occupying the steep slopes around the cities.
In this context, predictive models capable of identifying shallow landslide susceptible areas are fundamental tools for city planners and civil defense officials.These models offer a quick and effective way to determine areas that people should not exploit, evacuation plans, knowledge about risk zones, and regional planning.Deterministic and probabilistic approaches have been developed to assess geologic hazards [4].However, physically-based models are more frequently applied than empirical and statistical approaches (e.g., [5][6][7][8]).Physically-based models using Digital Elevation Models (DEMs) plus spatially distributed geotechnical information consider the mechanics of slope instability rather than empirical correlations among variables [9].These models generally couple a hydrologic model for the analysis of pore-water pressure regime with an infinite slope idealization for computing the safety factor (e.g., [10][11][12][13]).Thus, different physically-based models have been proposed, such as, Shallow Slope Stability Model (SHALSTAB) [10], physically-based Slope Stability Model (dSLAM) [11], Stability Index Mapping (SINMAP) [14], Transient Rainfall Infiltration and Grid-based Regional Slope-Stability Model (TRIGRS) [15], TRIGRS-unsaturated [16], SEEP/W and SLOPE/W [17], and HIRESSS [18].
The shallow landslides and debris-flows phenomena are often correlated because debris-flow starts after the occurrence of a previous shallow landslide [29].According to Carrara et al. [30], although physically-based models may be suitable for modeling the hydrological conditions leading to debris-flow initiation, they have several limitations when applied to predict the spatial distribution of debris-flows, i.e., to map debris-flow susceptibility.Consequently, for assessing debris-flow, a dynamical model can be applied considering either 1-D models [31,32], which move the flow in only one spatial dimension as a cross-section of a single pre-defined width or 2-D models [33] that move the flow in two dimensions, taking into account the topography in the plane surface and cross-section.Among the various numerical codes available in the current literature, the commercial FLO-2D is one of the most widely used [34][35][36][37].SHALSTAB is used to estimate the sediment volumes that move downstream adopted as input in the FLO-2D [38,39].
Thus, this paper aimed to develop a hazard map considering: (a) identification of areas prone to shallow landslides initiation and sediment volumes that moved downstream by SHALSTAB model and (b) determination of debris-flow reaches and final amount of volume material deposition in order to encompass all the areas that will be affected by both the processes.

Study Area
The study area comprised of the Quitite and Papagaio river basins located on the west slope of Maciço da Tijuca in Jacarepaguá neighborhood of the city of Rio de Janeiro.This region has an area of around 5 km 2 (Figure 1).In February 1996, the area experienced strong storms that triggered several shallow landslides and subsequent debris-flows in the "Baixada de Jacarepaguá" region (Figures 1 and 2).According to GEORIO [40], the study area experienced 250 mm of precipitation in 48 h between 12 and 13 February 1996.Strong rainfall occurs in this area due to the presence of an orographic barrier.The total volume mobilized by the debris-flow was around 132,000 m 3 (90,000 m 3 in Quitite River Basin and 42,000 m 3 in Papagaio River Basin).The debris-flow reached a speed of 5.3 m/s and 2.8 m/s in the Quitite and Papagaio River Basin, respectively.
Bedrock is composed of a complex combination of Precambrian high-grade metamorphic rocks with granite intrusions [40].The most frequent lithological unit is the highly foliated banded gneiss (Archer Gneiss).
With regard to land use and land cover, native forest occurs in the basin headwaters.In the middle part of the Quitite River Basin, there were pasture lands and mining activities were stopped after the mass movement events occurred.The Papagaio river basin contains small farms with typical agricultural activities.The basins contain anthropogenic elements, such as houses and streets that have been built over previous debris-flow deposits (Figure 2).

Methodology
In this paper, we used two models: SHALSTAB for shallow landslides mapping and FLO-2D for debris-flow mapping.According to Pirulli and Sorbino [41], three steps are necessary to run numerical analyses: (a) topography of the study area as a DEM; (b) initial sediment volume that moved downstream; and (c) guidelines for rheological parameter calibration.In this paper, the initial volume is calculated using two approaches: (a) photo-interpretation of the scar areas with field-work and (b) SHALSTAB model.If the SHALSTAB result is compatible with the photo-interpretation of the scar areas, this can be used to predict shallow landslides and initial volume in other areas with similar environmental conditions.
The determination of rheological parameters of debris-flow is obtained from the back-analysis technique considering the event of February 1996 and the FLO-2D model.Finally, a combination of SHALSTAB and FLO-2D models is used to delineate the areas where shallow landslides and debris-flows occurred in February 1996.

Elaboration of DEM and Photo-Interpretation of Landslide Scars and Debris-Flow Deposits
Both the SHALSTAB and FLO-2D models depend strongly on the quality of the DEM.The DEM was generated by interpolation of contour lines and drainage network from a 1:10,000 scale topographic map.A DEM with a spatial resolution of 2 m was generated for the SHALSTAB model, and a DEM with a resolution of 10 m was generated for the FLO-2D model.The low spatial resolution for the DEM of FLO-2D model was due to an extremely high computational effort.Stolz and Huggel [37] described that FLO-2D can have problems with high resolution DEMs resulting in extensive computing time or numerical instability.Moreover, most of the debris-flow events in the area were more than 10 m wide.Both DEMs were generated using the TOPOGRID module of ArcInfo software.
Guzzeti et al. [42] subdivided the methodological approach of remote sensing for landslide studies into three broad categories: (i) visual (heuristic) interpretation of optical images, including panchromatic, composite, false-color, and pan sharpened ("fused") images, (ii) analysis of multispectral images, including image classification methods and semi-automatic detection and mapping of landslides, and (iii) analysis of synthetic aperture radar (SAR) images.In Metternicht et al.'s [43] review on application of remote sensing techniques for landslide studies, photo interpretation was the most frequent remote sensing tool applied in landslide inventory mapping.The recognition of landslide characteristics has used normally aerial photographs, in addition to satellite imagery [44][45][46][47].Recently, use of satellite data has increased significantly, due to high resolution (HR) and very-high resolution (VHR) sensors, and improvements in both computer hardware and software for processing, visualization, and analysis of satellite images [42].Tofani et al. [48] confirmed this statement from a questionnaire about actual application in landslide detection and mapping in which 17 European countries demonstrated that aerial photos with visual interpretation techniques are in much use.A partial reason of the extensive use of aerial photographs is free access to large historical data with low costs for acquisition and processing.Mantovani et al. [49] verified that the most useful photographic scale is about 1:15,000 because it increases the chances of identifying landslide scar.
In this paper, landslide scars and debris-flow path and deposits were mapped by visual interpretation from aerial photographs that were taken about two months after the storm.A complete landslide inventory was elaborated according to the 2 m grid.Those scar maps served to validate both the landslide model and the debris flow simulation (Figure 3).

SHALSTAB Model Formulation
A shallow landslide is characterized by a straight slide plane [29].It is triggered during strong storms that produce high pore pressures in the zone of contact between the regolith and an impervious layer in the soil.The positive pore pressure alters the slope stability, reducing the internal shearing strength of the material and normal tension.Shallow landslides were responsible for about 38% of the mass movement events that occurred during 1962-1992 in Rio de Janeiro [21].
In this paper, we used the SHALSTAB model developed by Montgomery and Dietrich [10].The model performs an infinite-slope stability analysis assuming equilibrium (steady-state) conditions and flow parallel to the surface and uses Darcy's law to estimate the spatial distribution of pore pressures [10].This model combines a hydrological model with a slope stability model to identify susceptible landslide areas.
The hydrological model considers the model developed by O'Loughlin [50] to estimate soil saturation levels.This model calculates the saturation level based on upstream flows from a given point, slope angle, and soil transmissivity.Hydrological model can be solved based on the saturated proportion of the soil (h/z), assuming that saturated conductivity does not vary with depth [10], as follows: . ..sinθ where, h is the water-table height, z is the soil thickness, a is the upslope contributing area (m 2 ), b is the grid cell size (m), θ is the local ground slope (degrees), T is the soil transmissivity (m 2 /day), and Q is the steady-state rainfall intensity.Infinite-slope theory can be solved based on the h/z ratio, which has the following equation [10]: The combination of hydrological and slope stability model uses the h/z ratio resulting from the following equation:

Application and Validation of SHALSTAB Model
In the application of SHALSTAB model, the geomorphological parameters (contributing area and slope) were directly obtained from DEM. Soil parameters were extracted from Guimarães et al. [23], who made 125 simulations with different values of soil properties from the SHALSTAB model.Guimarães et al. [23] observed that the best-fit model using back-analysis has C/z = 2 kPa•m −1 , ϕ = 45°, and ρ s /ρ w = 1.5, which provided the best result among all possible combinations of soil parameters for the study area.In general, the best parameterization values have high ϕ, modest to low C/z, and 1.5 < ρ s /ρ w < 1.75, whereas parameterizations with low ϕ in combination with low C/z and ρ s /ρ w > 2 resulted in poor model performance.
SHALSTAB validation procedure can be done according to Dietrich et al. [19], who compared the distribution of shallow landslides modeled with that observed from photo-interpretation. If the model is successful, the mapped landslide scars and density of shallow landslides (number of landslides per unit area) should be much more common in the least stable areas.Thus, each digitized polygon representing a landslide scar was overlaid on the grid of Q/T values (Equation (3)), where the cell with lowest Q/T value within the polygon represents the least stable site and therefore controls the site stability.Instead of the Q/T values, the instability categories defined by log (Q/T) intervals can also be used.The more unstable categories in scar polygon may be closer to the actual condition at instability.

FLO-2D Model Formulation
Debris-flow are rapid movements in which the materials have a high viscosity mobilizing a large volume of material (including large rock blocks) in a short period and over long distances [4].Various processes start the debris-flow, including shallow landslides.Borga et al. [51] consider shallow landslides being the most common debris-triggering process.This type of mass movement is highly destructive due to the strong impact of their flow and power to reach areas with lower slopes.
We used the FLO-2D model developed by O'Brien and Julien [52], which simulates debris-flow using finite-difference routines in two dimensions.The modeling of a debris-flow is controlled by topography and is performed using numerical integration of movement and continuity equations and flow-resistance parameters.Sediment flows are simulated using hyper-concentrated sediment-flow routines with continuous flows, which enable to predict the behavior of fluid flow from a fluid matrix governed by sediment concentration.FLO-2D uses a quadratic rheological model, developed from field and laboratory mudflow data, and enables appropriate simulations of flooding conditions ranging from clear water to hyperconcentrated sediment flows [33].O'Brien and Julien [52] define the FLO-2D rheological model with the following approximation: where, S f is the total friction slope, S y is the sum of the yield slope, S v is the viscous slope, and S td is the turbulent-dispersive slope.The total friction slope can be written as follows: where, τ y is yield stress, γ m is the specific weight of the slurry, h is flow depth, K is an empirical resistance parameter, η is the fluid viscosity, V is flow velocity and n td is Manning's roughness coefficient.The yield stress and fluid viscosity are defined from [52]:  where, α 1 and β 1 are empirical coefficients defined by laboratory experiments [53].Simulation of debris-flow with the FLO-2D model requires measurement of the flow tension and viscosity [54].

Determination of Rheological Parameters of Debris-Flow
In this paper, back-analysis approach was used in order to reproduce, virtually, the rheological properties of the debris-flow event that occurred in February 1996.FLO-2D model is the one that is mostly used for debris-flow back-analysis [33][34][35][55][56][57][58][59][60].This procedure focused on the simulations of a reliable range of parameters capable of describing the occurred phenomenon (Figure 3).
In the application of FLO-2D model, the topographic attributes were obtained from DEM, while the initial material volume was defined from estimates of landslide scars.The rheological parameters were estimated from the back-analysis technique considering 150 simulations.The rheological coefficients of each simulation were back-calculated by fitting with the field data (in terms of area and depth of debris-flow).Five points taken immediately after the event were used to validate the model results (Figure 4).The best simulation for debris-flow area considers the absolute difference between the real and simulated results.The best simulation for thickness was determined by greater similarity with the real event that occurred in 1996, considering the root-mean-square deviation method of the five points.
In the literature, the range of values adopted for rheological tests using empirical back-analysis are extensive and vary depending on the study area [36,55,57].Value ranges selected for rheological properties were based on works of Macias et al. [61] and O'Brien and Julien [51], besides extreme values reported for debris-flow events.Macias et al. [61] described the mechanical behavior from back-analysis of debris-flow occurred in the Quitite and Papagaio River Basins to the event of February 1996; obtaining viscosity range from 0.092 to 3.44 kPa•s and yield stress range from 12.06 to 24.42 kPa.O'Brien and Julien [52] analyzed the main physical and mechanical properties of the debris-flow from the literature.Because of the lack of knowledge in the study area on rheological parameters, we chose to conduct an intensive empirical test evaluating broad ranges of values.Figure 5 shows the values adopted for the empirical tests, in which the viscosity varies from 0.01 to 10.0 kPa•s; yield stress ranges from 0.0002 to 48 kPa; laminar flow resistance from 0 to 4,000, and time 1-4 h.In the FLO-2D model, if K = 0, the value of K is automatically computed from Manning's n-value.

Combining the SHALSTAB and FLO-2D Models
According to Gentile [38], an integrated approach to debris-flow risk analysis using SHALSTAB and FLO-2D allows a good determination of the real risk conditions.Thus, new debris-flow simulations were performed considering the best-fit values for the rheological properties (obtained in the back-analysis stage) and the initial material volume estimated from SHALSTAB results.The delimitation of the susceptible areas to shallow landslides considers values of log Q/T lower than −3.1.The coupling models (SHALSTAB and FLO-2D) enable to delimitate susceptible area to mass movements (landslide and debris-flow) in the study area or in other areas along the Serra do Mar characterized by analogous environmental conditions.Result of the coupled models (SHALSTAB/ FLO-2D) was validated by comparison between the simulated and observed results (thickness and area of the debris-flow deposits).

Results of the SHALSTAB Model
The SHALSTAB map shows seven classes that vary from stable (tan θ ≤ tan ϕ [1 − ρ w /ρ s ]) to unstable (tan θ > tan ϕ) (Figure 6(a)).The stable class is not susceptible to shallow landslides even with strong rainfall.In contrast, the unstable class is very susceptible to landslides even without rainfall.Notably, unstable areas are predominantly composed of bedrock outcrops.As the ratio Q/T was a small number, a log function was applied.The five intermediate classes demonstrated degrees of instability.To better understand how those class intervals were defined, we considered transmissivity to be around 65 m 2 /day, log Q/T = −3.1 represents a rainfall amount of 51 mm/day, and log Q/T = −2.2represents a rainfall amount of 404 mm/day.Therefore, the intermediate classes interval is −0.3, that corresponds a factor 2, i.e., the effect increases twice for precipitation.Thus, the log Q/T can be calculated between real, reasonable and natural values.In addition, each landslide was classified as associated with the minimum log (Q/T) class within its boundary.The comparison of model output with the landslide scars mapping shows that almost all of the 89 scars were detected (Figure 6(a)).The model placed only one scar in the stable class.In our approach, the unstable class occupied only 2% and the <−3.1 class only 3% of the total area (Figure 6(b)).Thus, the high correspondence between observed shallow landslides and potentially unstable ground attest to the quality of SHALSTAB results.Due to the accuracy of the results obtained, tests with other models (such as TRIGRS) are not necessary.

Results of the Back-Analysis for Rheological Properties
In the back-analysis for rheological properties, the initial material volume in FLO-2D is calculated from scars map.Among the 150 simulations, only five simulations showed area and thickness similar to the debris-flow event of February 1996.Table 1 shows the values of rheological properties obtained with the back-analysis.The top five simulations had the same time (two hours).The best-fit simulation shows viscosity of 0.092 kPa•s, flow tension of 0.02 kPa, K equal to zero, and simulation time of two hours.The only change in the second best-fit simulation was the yield stress (0.002 kPa).Figure 7 shows the graph that represents the Square Meters per Final Order Positioning (FOP) with the best simulations to the area of debris-flows, which are compared with the measured area by photo interpretation (731,539 m 2 , red dashed line).The two best simulations exhibit absolute differences of debris-flows area less than 2,000 m 2 .Therefore, the error introduced by the model is of an extremely low percentage corresponding to 0.27 of the total area of debris-flows, showing a high conformity with reality.FOP5 simulation achieved an area of 739,600 m 2 , which when compared with a true area shows an absolute difference of approximately 8,000 m 2 ; although it is four times the FOP1, it also has a low percentage considering the debris-flow total area of 1%.The thickness analysis of each simulation was performed by the root-mean-square deviation among the actual and simulated values for the five points of fieldwork.Figure 8 shows the best-fit simulations for flow thickness, which are coincident with those obtained for the area.The two best-fit models obtained root-mean-square deviations equal to 3.2 m, i.e., comparable values with the event of February 1996.These rheological parameters are close to the best result using back-analysis obtained by Armento et al. [55] for the events and debris-flow in Acquabona K19 (Boite Valley, near Cortina), wherein the viscosity is equal to 0.097 kPa•s and yield stress equals 0.004177 kPa.

Results of the Integration of Models SHALSTAB and FLO-2D
In the integration of models, output data from SHALSTAB are used as input data in FLO-2D.The initial material volume in FLO-2D is calculated from the area delimited by SHALSTAB result (log Q/T low than −3.1) multiplied by the average height of debris-flows described in GEORIO (1 m) [40] (Figure 9).This class of Q/T less than −3.1 was chosen because it represents the highest risk areas defined by the model, since the unstable class corresponds only to the outcrop areas.The rheological properties are derived from the two best-fit models of the back-analysis.Both simulations obtained debris-flow areas (749,600 m 2 and 753,800 m 2 ) compatible with the event of February 1996 (731,539 m 2 ) and thickness with low values of root-mean-square deviation of around 3.2 m. Figure 10 shows a perspective view between the debris-flow area bounded just after the February 1996 event and the best-fit scenario of debris-flow event simulated with rheological parameters obtained through back-analysis.

Conclusion
The accurate prediction of mass movements can reduce damages and guide future planning and zoning.Several mathematical models have been developed so far to describe the mass movement events.In this paper, the landslide model (SHALSTAB) was integrated with the debris-flow simulation model (FLO-2D).In the FLO-2D model, the sediment volumes that moved downstream were estimate from the SHALSTAB model.
Both models provided good results to predict the mass movement events (landslide susceptibility, runout distances and deposits thickness).The empirical method used enabled us to identify landslide susceptibility zones in the landscape and to obtain topographic profiles of slides and probable deposition sites.An accurate replication of the event that occurred in February 1996 was obtained in terms of area and thickness of the debris-flow.
However, obtained data from back-analysis must be carefully evaluated due to the following reasons: (a) the model does not consider the natural changes of the rheological parameters during the event [62,63], (b) the results are very sensitive to the wide variability of rheological parameters [35,62], and (c) the set of best-fit parameters can diverge from reality.Research studies emphasize the discrepancy between simulated rheological parameters and that obtained in the laboratory [35,58,64].Prochaska et al. [64] developed a database of rheological properties (viscosity and yield strength) from literature data for both procedures' back-analysis [35,[65][66][67][68][69][70] and laboratory tests [70][71][72][73], and found significant differences between the two methods.Therefore, rheological parameters are inevitably affected by uncertainties considering either empirical or experimental procedures.
Despite the uncertainties of the numerical model, the empirical procedure adopted in this study allowed to recreate different scenarios that reproduced the mass movement event.Thus, SHALSTAB-FLO-2D

Figure 1 .
Figure 1.Location map of the Quitite and Papagaio basins.

Figure 2 .
Figure 2. Overview of various mass movements that reach the study area.

Figure 3 .
Figure 3. Map of the landslide scars (with 89 scars mapped) and debris-flow paths (with 731,539 m 2 area mapped) in the study area. and

Figure 4 .
Figure 4. Points used for checking the thickness depth from debris-flow of 1996.

Figure 5 .
Figure 5. Values of rheological parameters used in the back-analysis.

Figure 6 .
Figure 6.(a) Comparison between SHALSTAB map and the landslide scars; (b) distribution of the landslide susceptibility classes for the scars area.

Figure 7 .
Figure 7. Graph comparing the five best-fit models to the debris-flows area (Final Order Positioning-FOP).The dashed red line represents the actual area (731,539 m 2 ).

Figure 8 .
Figure 8. Graph of root-mean-square deviations of maximum debris-flows thicknesses.

Figure 9 .
Figure 9. Example of initial material volume (5,700 m 3 ) used by FLO-2D model delineated by SHALSTAB map.The volume is calculated considering a depth of 1 m and estimated area by SHALSTAB.

Figure 10 .
Figure 10.Perspective views: (a) debris-flow area bounded just after the February 1996 event; and (b) simulation result from the coupled models (SHALSTAB/FLO-2D) considering the set of rheological parameters with best-fit in the back-analysis.

Table 1 .
Ranking positioning of the five best-fit models with combinations of parameters simulated where FOP is the Final Order Positioning.