Regional Analyses of Rainfall-Induced Landslide Initiation in Upper Gudbrandsdalen (South-Eastern Norway) Using TRIGRS Model

: In Norway, shallow landslides are generally triggered by intense rainfall and/or snowmelt events. However, the interaction of hydrometeorological processes (e.g., precipitation and snowmelt) acting at different time scales, and the local variations of the terrain conditions (e.g., thickness of the surﬁcial cover) are complex and often unknown. With the aim of better deﬁning the triggering conditions of shallow landslides at a regional scale we used the physically based model TRIGRS (Transient Rainfall Inﬁltration and Grid-based Regional Slope stability) in an area located in upper Gudbrandsdalen valley in South-Eastern Norway. We performed numerical simulations to recon-struct two scenarios that triggered many landslides in the study area on 10 June 2011 and 22 May 2013. A large part of the work was dedicated to the parameterization of the numerical model. The initial soil-hydraulic conditions and the spatial variation of the surﬁcial cover thickness have been evaluated applying different methods. To fully evaluate the accuracy of the model, ROC (Receiver Operating Characteristic) curves have been obtained comparing the safety factor maps with the source areas in the two periods of analysis. The results of the numerical simulations show the high susceptibility of the study area to the occurrence of shallow landslides and emphasize the importance of a proper model calibration for improving the reliability.


Introduction
The study of rainfall-induced landslides is an important research topic for the improvement of risk mitigation strategies for current and future climate scenarios. For this reason, over the years a number of solutions have been developed on this issue, often based on statistical-probabilistic techniques (e.g., [1][2][3][4]). However, even though these methods can be very effective in the identification of the areas more prone to slope failure (i.e., landslide susceptibility), the same cannot be said about the temporal prediction of such phenomena, which still represents a challenging research topic. In this respect, also considering the improvements of environmental sensing technology, different landslide early warning systems (LEWSs) have been recently drawn up, both at slope scale (see [5] and references therein) and over national/regional scale (see [6] and references therein), such as the one developed. These systems basically monitor the most significant variables related to the triggering conditions (e.g., rainfall, soil moisture, displacement) and, when given thresholds are exceeded, warning messages are generated to alert the public and authorities [7]. Every early warning system operational worldwide has a different number of thresholds and algorithms to activate a warning level. Such thresholds can be defined using different tools, including physically based models that reproduce, often via numerical simulations, the progressive onset of the instability within the potentially unstable soil during heavy rainfall events (e.g., [8][9][10][11]). However, a correct spatial evaluation of the input parameters (e.g., rainfall field, initial water table, thickness and physical-mechanical parameters of the soil cover) is a fundamental aspect for the reliable application of such models. In this respect, back-analyses of preceding landslide events is a common strategy to calibrate the inputs parameters of physically based models.
The results presented in this paper have been obtained applying the physically based model TRIGRS [12], which is a model specifically devoted to the prediction of rainfallinduced landslides. Since the model calculates the slope stability conditions in response to specific rainfall events, we have back-analysed two separate landslide events that occurred in the upper Gudbrandsdalen, Eastern Norway, on 10 June 2011 and 22 May 2013 [13]. The unstable areas obtained by the application of TRIGRS (Safety Factor lower than 1) have been compared with the source areas observed in situ. The outcome of this comparison was then evaluated through a ROC (receiver operating characteristic) curve analysis, allowing is to quantify the performance of the two predictions.

The Case Study
Gudbrandsdalen is a major north-south trending valley in the mountainous part of the former Oppland county (now Innlandet county) in South-Eastern Norway ( Figure 1). Specifically, the study area (Veikledalen) is in the upper part of this region, near the village of Kvam. The area is characterized by hill slopes mainly consisting of green or grey phyllite, within which layers of grey sandstone, feldspathic sandstone or quartzite can be found [14]. These slopes have been largely carved by glaciers, whose erosive action is also testified by the presence of till deposits of variable thickness. In general, the thinner deposits (0.5 m or less) are in the upper part of the slopes, while in the lower parts, the thickness can reach several tens of meters. The valley floor is mainly filled by glacifluvial and fluvial deposits, which can be locally covered by Holocene colluvial fans generated at the toe of the numerous gullies cutting the till deposits [15]. In this respect, the area has experienced different landslide events recently, mainly channelized debris flows generally initiated as shallow slides within the soil cover (i.e., soil slips). The most recent events occurred in May 2008, June 2011 and May 2013. In particular, the last two occurred in very similar hydrometeorological conditions that also triggered flood events along the Glomma river.
In the period 6-12 June 2011, significant rainfall occurred, and the simultaneous high temperature led to strong snow-melting process. Due to such a combined effect, several shallow landslides were triggered in the night between the 9 and 10 June, causing the closure of numerous roads in Gudbrandsdalen. More than 270 people were evacuated, and the overall damages cost is estimated at 800 million Norwegian Kroner (NOK), ca. 100 million USD. On 22 May 2013, the combination of intense rainfall and prolongated, intense snow melting caused widespread flooding and many shallow landslides in the same area [13]. The worst damage was recorded in the town of Kvam, where more than 200 residents were evacuated. As in 2011, the two main roads, the National Highway Rv3 and the European highway E6, were closed at several locations. The damage cost is estimated at a total of 1500 million NOK, ca. 170 million EUR [16]. In this work, the over 100 landslides triggered during both events have been identified through the analysis of high-resolution aerial orthophotos by distinguishing, for each landslide, the source area from the propagation and deposition zones ( Figure 2).

The TRIGRS Model
TRIGRS is a widely used numerical model in the framework of landslide susceptibility analyses (e.g., [17][18][19]). It can reproduce, in a simplified form, the triggering process of

The TRIGRS Model
TRIGRS is a widely used numerical model in the framework of landslide susceptibility analyses (e.g., [17][18][19]). It can reproduce, in a simplified form, the triggering process of

The TRIGRS Model
TRIGRS is a widely used numerical model in the framework of landslide susceptibility analyses (e.g., [17][18][19]). It can reproduce, in a simplified form, the triggering process of rainfall-induced shallow landslides considering relatively few and simple input parameters. Specifically, it performs a slope stability analysis over large areas through the infinite slope limit equilibrium method, considering a transient 1-D analytic solution for pore water pressure induced by rainfall infiltration based on [20] linearized solution of the Richards' equation. In other words, it is possible to consider a space-time variability of rainfall input (I Z ) that, together with a simple runoff routing scheme, allows calculating the pressure head response ψ(Z, t) for each grid cell starting from initial saturated or unsaturated conditions. In this latter case, the model uses a simple analytic solution proposed by [21] to linearize the Richards' equation, considering four hydrodynamic parameters: the saturated (θ s ) and residual (θ r ) water content, the saturated hydraulic conductivity (K s ) and a specific parameter linked to the pore size distribution (α G ). If the amount of infiltrating water reaching the water table exceeds the maximum amount that can be drained by gravity, TRIGRS simulates the water-table rise, comparing the exceeding water quantity to the available pore space directly above the water table or capillary fringe ( Figure 3) and then, for each time step, it applies the water weight at the initial top of the saturated zone to compute the new pressure head [22]. For the calculation of the safety factor in the unsaturated configuration, the pressure head is multiplied by the effective stress parameter χ, as suggested by [23]: where θ is the soil water content. This approximation represents a simplified form of the suction-stress characteristic curve [24,25]. The Safety Factor at a specific depth Z is then calculated according to: where δ is the slope angle, c' is the effective soil cohesion, ϕ' is the effective friction angle of the soil, γ w is the unit weight of groundwater and γ is the total unit weight of soil. A more detail explanation of the assumptions, fundamental concepts and equation derivation of the TRIGRS model can be found in [12,22].     To calculate the Safety Factor in response to specific rainfall events, TRIGRS uses a series of input parameters (Table 1). Two different detailed (5 × 5 m) digital elevation models (DEM) are available for the study area. The first one was dated 2010, i.e., before the 2011 landslide event, and the second one was dated 2012. Considering this aspect, we decided to use the 2010 DEM for deriving the slope map of the area in the case of the 2011 event simulation, and the 2012 DEM for the 2013 event.
For the estimation of the spatial variation of the thickness of the surficial cover, a modified version of the model proposed by [26] has been applied to the study area. The original model correlates the soil depth to the local slope angle according to the following equation: where α i is the slope value at pixel i, while α max and α min are the maximum and minimum slope values encountered in the study area, h i is the soil thickness computed at pixel i, h max and h min are the maximum and minimum soil thickness values measured in the area. However, Equation (3) was modified in such a way that a series of real thickness data deriving from different sources could be used, then improving the reliability of the model itself. The available data are summarized in Table 2.  Table 2. Thickness data of the surficial cover available for the study area.

Type of Data Points Source
Field data-GPS position 3 The first dataset [29] is composed by measurements of the surficial cover thickness for which the GPS position is available from field surveys (Figure 4), whereas the second series of data are measurements whose position has been identified by using the information reported in [30]. Regarding the third dataset, we first identified from high-resolution aerial orthophotos 36 points within the footprint areas where the bedrock was clearly exposed after the landslide events. Then, we made the difference between two high-resolution (1 × 1 m) DEMs derived from LIDAR data ( Figure 5). Since the DEMs had been released before (i.e., 2000) and after (i.e., 2012) the 2011 landslide event, it was possible to measure the effective thickness of the surficial cover in each of the 36 points. Afterwards, all the 43 thickness values (and the corresponding slope values measured in each point) were plotted in a graph "thickness vs. tangent of slope", and a least squares linear fit of the points was used as the soil thickness model ( Figure 6). It is important to note that for the numerical simulations we decided to neglect the areas with thickness values higher than 2 m in the lower part of the valley, since the real thickness of the surficial cover is considerably higher than the expected, and it could not be simplified with a linear trend. On the other end, we imposed a constant value of 0.4 m (i.e., the lowest thickness value measured on field) for the steepest areas. However, the thickness of the surficial cover also depends on the geological and geomorphological features of the area, which are those typical of glacial environments. In this respect, the map of the Quaternary deposits available at the Geological Survey of Norway (NGU) shows the presence of "shallow" glacial deposits in the upper part of the slopes, whereas the "thick" glacial deposits can be found only in the lower part of the valleys. Considering this last aspect, we decided to impose a constant value of 0.4 m in the areas characterized by the presence of shallow glacial deposits, and a null value for the areas where the bedrock outcrops. were plotted in a graph "thickness vs. tangent of slope", and a least squares linear fit of the points was used as the soil thickness model ( Figure 6). It is important to note that for the numerical simulations we decided to neglect the areas with thickness values higher than 2 m in the lower part of the valley, since the real thickness of the surficial cover is considerably higher than the expected, and it could not be simplified with a linear trend.
On the other end, we imposed a constant value of 0.4 m (i.e., the lowest thickness value measured on field) for the steepest areas. However, the thickness of the surficial cover also depends on the geological and geomorphological features of the area, which are those typical of glacial environments. In this respect, the map of the Quaternary deposits available at the Geological Survey of Norway (NGU) shows the presence of "shallow" glacial deposits in the upper part of the slopes, whereas the "thick" glacial deposits can be found only in the lower part of the valleys. Considering this last aspect, we decided to impose a constant value of 0.4 m in the areas characterized by the presence of shallow glacial deposits, and a null value for the areas where the bedrock outcrops.  As far as the geotechnical and hydraulic parameters of the debris cover are concerned, unfortunately no data deriving from laboratory or field tests are available for the study area, as well as in many other areas in Norway. Therefore, as a first step we decided to use the values reported in [27]. In this paper, the authors calibrated the soil parameter values through Monte Carlo simulations using as constraint a series of debris flows located near Otta, northwest of the study area. Specifically, starting from the spatial distribution of the Quaternary deposits, deriving from two maps at 1:20,000 and 1:50,000 scale, they perform numerical simulations by imposing as target the maximization and minimization of areas classified as unstable for a 1000-year return period rainfall event and in absence of precipitation, respectively. The resulting values related to the glacial deposits are reported in Table 1.
For the estimation of the hydrodynamic parameters employed by TRIGRS in the unsaturated modality, we used the ROSETTA Lite module [31] available within the HYDRUS-1D software [32]. This module uses a database of measured water retention and other properties for a wide variety of media. Specifically, for a given grain-size distribution the model estimates a soil water retention curve (SWRC) based on the hydraulic model proposed by [33] and with good statistical comparability to known retention curves of other media with similar physical properties [34]. In this work, an average grain size distribution has been considered according to four samples collected in the case study area [35]. However, the α G parameter corresponds to a fitting parameter of the SWRC based on a different ( [36] hydraulic model. This parameter was then estimated using the conversion formula from [28], which defines a correspondence between the above-mentioned hydraulic models through the capillary length approach [37].  The initial water-table depth and the background rainfall rate are, instead, closely related to the preceding hydrological conditions, which can highly affect the spatio-temporal occurrence of rainfall-induced landslides [38]. Such parameters were estimated from the information reported in Varsom Xgeo (www.xgeo.no). It is an open access decision-support tool developed from the cooperation between the Norwegian Water Resources and Energy Directorate (NVE) and other governmental agencies such as the Norwegian Meteorological Institute, the National Public Road Administration, the National Rail Administration and the National Mapping Authority. Specifically, it is composed of a web portal and several modules that store different types of data such as: gridded maps (spatial resolution: 1 km 2 ) of meteorological data and simulated hydrological, glaciological and hydrogeological conditions; real-time data recorded at automatic stations; field observations (i.e., snow conditions, landslide events); thematic maps; warnings maps; and satellite images. The tool is daily employed for national forecasting and warning services (e.g., floods, landslides and snow avalanches) carried out by NVE over the entire national territory [39].   [26]. The three different types of thickness data used as constraints of the model are also reported.  [26]. The three different types of thickness data used as constraints of the model are also reported.
The hydrological conditions reported in Varsom Xgeo derive from simulations performed at the daily scale with HBV, which is a semi-distributed conceptual rainfall-runoff model [40,41]. Specifically, two versions of the model are daily employed at NVE: a river-catchment-type version that is calibrated over 150 river catchments and a spatially distributed-type version (also known as GWB-Gridded Water Balance). The latter version considers each of 385,000 cells over the entire Norwegian territory as a separate basin, for which a specific water balance is performed. The model runs 24 h a day and uses precipitation and air mean temperature as inputs. The main structure of the HBV model comprises four storage components: snow, soil moisture, and upper and lower runoff zone. The model can reproduce not only runoff, but also infiltration and snow melting, then estimating the corresponding variations of groundwater level, soil saturation, soil frost, etc. In particular, the groundwater-level map provides information on the groundwater level compared to the average for the same day for the reference period 1981-2010. The degree of soil-saturation map (Figure 7) shows the percentage of soil saturation, which describes the relationship between current soil-water storage compared to the maximum soil-water storage simulated with the HBV model in the reference period 1981-2010. In this work, the degree of soil saturation was used for the estimation of the background rainfall rate (I ZLT ). Such a parameter represents the steady (initial) surface flux, which can usually be approximated by the average precipitation in recent weeks or months that is needed to maintain the initial conditions [22]. However, information about methods for estimating I ZLT is rare in the literature. If the soil is saturated, I ZLT can be the same as hydraulic conductivity, while it can be zero for dry soil. For this reason, several authors [42][43][44][45] calculated this parameter as a fraction of the hydraulic conductivity considering the long-term weather conditions prior the investigated period. Thus, we spatially varied the background rainfall rate as a function of the degree of soil-saturation map available for both the events: the higher the degree of saturation, the higher the I ZLT value, with the maximum value equal to Ks for Sr = 100%. In this respect, the degree of saturation simulated by the HBV model accounts for the contribution of all the hydrologic processes on the slopes, including the snowmelt, which had a relevant role in developing the instability conditions during the 2013 event. Finally, the rainfall field related to the 2011 and 2013 events was obtained from hourly meteorological radar measurements operated by the Norwegian Meteorological Institute (Figure 8).
long-term weather conditions prior the investigated period. Thus, we spatially varied the background rainfall rate as a function of the degree of soil-saturation map available for both the events: the higher the degree of saturation, the higher the IZLT value, with the maximum value equal to Ks for Sr = 100%. In this respect, the degree of saturation simulated by the HBV model accounts for the contribution of all the hydrologic processes on the slopes, including the snowmelt, which had a relevant role in developing the instability conditions during the 2013 event. Finally, the rainfall field related to the 2011 and 2013 events was obtained from hourly meteorological radar measurements operated by the Norwegian Meteorological Institute (Figure 8).   45] calculated this parameter as a fraction of the hydraulic conductivity considering the long-term weather conditions prior the investigated period. Thus, we spatially varied the background rainfall rate as a function of the degree of soil-saturation map available for both the events: the higher the degree of saturation, the higher the IZLT value, with the maximum value equal to Ks for Sr = 100%. In this respect, the degree of saturation simulated by the HBV model accounts for the contribution of all the hydrologic processes on the slopes, including the snowmelt, which had a relevant role in developing the instability conditions during the 2013 event. Finally, the rainfall field related to the 2011 and 2013 events was obtained from hourly meteorological radar measurements operated by the Norwegian Meteorological Institute (Figure 8).

Results and Discussions
The results of the numerical simulations are presented in Figure 9 in terms of temporal evolution of slope instability at the catchment scale for the 2011 event. It is important to note that just before rainfall occurs (12 p.m.), no cell is indicated as unstable by the model. After a progressive increase of the instability level in the following hours, the most critical phase verifies between 12 a.m. and 1 a.m., when the number of unstable cells jumps from 10,878 to 30,262 in response to a high rainfall peak (Table 3). Afterwards, only a slight increase in the number of unstable cells is predicted until the end of the rainfall event (11 a.m.). In the case of the 2013 event, the TRIGRS simulations describe a different temporal evolution of the event. In particular, the onset of the most critical conditions occurs since the beginning of rainfall, and the increase of unstable cells is more continuous than in the 2011 case (Table 3). At the end of the event, 52,415 cells result as unstable, against 33,692 simulated for the 2011 event. This result emphasizes how the higher soil-moisture content in 2013 represents a more critical initial condition, which leads to a higher number of unstable cells against a lower rainfall intensity. However, in both cases the timing of the occurrence of the landslides substantially agrees with the actual sequence of the events. scale is not straightforward: in our opinion, this point represents the main drawback of physically based models such as TRIGRS. Another limitation of the model concerns the difficulty to simulate complex triggering mechanisms. However, as it is shown in this work, an in-depth analysis of the parameterization process allows consideration of different hydrologic processes occurring on the slopes, including the snowmelt. Furthermore, TRIGRS can be applied over areas with incomplete landslide inventories. The latter aspect represents a clear advantage over other methods for landslide susceptibility studies (such as heuristic and statistical ones), even though it would still be better to validate the resulting susceptibility map through the comparison with actual landslide scenarios, as shown in this work.   After the reconstruction of the instability scenarios over time, we evaluated the spatial predictability of the failure process through a ROC (receiver operating characteristic) curve analysis. Thus, we compared the two FS maps corresponding to the final state of both events (i.e., 11 a.m. and 12 a.m., for 2011 and 2013 events, respectively) with the two inventories of the source areas, in order to verify if the model was able to catch the actually unstable areas (TP-True Positive) and, at the same time, to minimize the number of erroneous failure predictions over stable areas (FP-False Positive). Specifically, the curve relates, for different threshold values, the True Positive Rate (TPR) with the False Positive Rate (FPR), that can be calculated as: where TN and FN are, respectively, the number of correct and wrong nonfailure predictions. The ROC curve allows quantifying the predictive capability of the model: the larger the area under curve (AUC), the better the prediction (Figure 10). In this respect, AUC is equal to 0.940 for the 2011 event, while it is slightly higher for the 2013 event (0.947). Furthermore, in both cases the percentage of TPs, with respect to the total landslide areas, is quite high (i.e., 74.6% and 86% for 2011 and 2013 event, respectively). Although such an outcome appears very satisfactory, it is fair to highlight that a large part of this result is related to the higher proportion of stable observations with respect to unstable ones, which turns also into very high TN values for both events (Tables 4 and 5). As regards the "unstable" areas, even the graphical comparison between the two final FS maps and the spatial distribution of the source areas shows a general overestimation of the actual slope instability, in particular for the 2013 event ( Figure 11). This overestimation can be partially explained by the poor information available on several input parameters. In this respect, we performed a sensitivity analysis for evaluating the weight of each parameter on the back-analysis results. Specifically, the mean values derived from [27] for five input parameters (Table 1) were varied according to the corresponding standard deviations. We also assumed single soil-thickness values for "shallow" and "thick" glacial deposits, always on the basis of the values reported in [27]. Considering that the 2013 event is characterized by the most critical initial soil-moisture conditions, we decided to reproduce only this event for emphasizing the effect of the parameter variation on the slope stability conditions. The results of the analysis (Table 6) show that the soil mechanical parameters (friction angle and cohesion) greatly affect the simulation outcome in terms of number of unstable cells, since even small variations can induce conditions of complete stability or instability. On the contrary, soil unit weight has a significantly lower impact, which is comparable with that of the hydraulic conductivity, while the hydraulic diffusivity does not induce any variation in the final number of unstable cells (at least within the considered range of values). As regards the soil thickness, a not negligible number of slope failures can be noticed from the very beginning of the event, while at the end an underestimation and overestimation of the instability processes results, respectively, with the lowest and highest thickness values. The latter result then suggests that the use of single soil-thickness values over large areas represent an oversimplification that may induce a distorted estimation of the slope instabilities. However, it is important to stress that estimating the spatial variation of the soil thickness and other input parameters at the catchment scale is not straightforward: in our opinion, this point represents the main drawback of physically based models such as TRIGRS. Another limitation of the model concerns the difficulty to simulate complex triggering mechanisms. However, as it is shown in this work, an in-depth analysis of the parameterization process allows consideration of different hydrologic processes occurring on the slopes, including the snowmelt. Furthermore, TRIGRS can be applied over areas with incomplete landslide inventories. The latter aspect represents a clear advantage over other methods for landslide susceptibility studies (such as heuristic and statistical ones), even though it would still be better to validate the resulting susceptibility map through the comparison with actual landslide scenarios, as shown in this work.        Table 3 Final n. Table 3 Friction angle (φ') 32 3.

Conclusions
In this work we applied a well-known numerical model for rainfall-induced landslide prediction (TRIGRS) in a context where different hydrologic processes occur on the slopes, including snowmelt. In this sense, we put great efforts into exploiting the existing data available over large areas, trying to define a reliable procedure for the definition of the input parameters required by the model. The results of the performed numerical sim- Figure 11. Graphical comparison between the spatial distribution of the 2011 (left) and 2013 (right) landslide source areas and the areas predicted as "unstable" (i.e., FS ≤ 1) at the end of the event by TRIGRS.

Conclusions
In this work we applied a well-known numerical model for rainfall-induced landslide prediction (TRIGRS) in a context where different hydrologic processes occur on the slopes, including snowmelt. In this sense, we put great efforts into exploiting the existing data available over large areas, trying to define a reliable procedure for the definition of the input parameters required by the model. The results of the performed numerical simulations indicate that TRIGRS is able to reproduce the 2011 and 2013 landslide events, at least in terms of temporal evolution of the instability process. In this respect, it is worth noting the role of the initial conditions in the timing of occurrence of shallow landslides. In fact, the instability arises quite early in the 2013 case due to the wetter initial conditions. Furthermore, the number of unstable cells at the end of the event is higher than that of 2011 event, despite the lower rainfall intensity. Regarding the spatial-predictive capability of the model, the ROC analyses indicate high values of AUC in both cases, although TRIGRS tends to overestimate the unstable areas, probably due to a degree of uncertainty around the estimation of input values. In this sense, a sensitivity analysis performed on several input parameters suggests that the soil mechanical parameters and the soil thickness have the greatest impact on the slope-stability conditions. Future efforts should then focus on innovative methods and techniques for the definition of input parameters over large areas, which can enhance the reliability of landslide scenarios depicted by physically based models such as TRIGRS.