Claudio Bosco and Graham Sander,

Loughborough University, Department of Civil and Building Engineering Loughborough, United Kingdom.

*This paper proposes a semi-quantitative integrated methodology for a robust assessment of soil erosion rates in data-poor regions affected by landslide activity. It combines heuristic, empirical and probabilistic approaches. The architecture might contribute as a novel component to simplify future integrated analyses of the potential impact of wildfires or vegetation types and distributions on sediment transport from water-induced landslides and erosion.*

**Abstract**:

Rainfall-induced landslides and soil erosion are part of a complex system of multiple interacting processes, and both are capable of significantly affecting sediment budgets. These sediment mass movements also have the potential to significantly impact on a broad network of ecosystems, in terms of their health, functionality and the services they provide. To support the integrated assessment of these processes it is necessary to develop reliable modelling architectures. This paper proposes a semi-quantitative integrated methodology for a robust assessment of soil erosion rates in data-poor regions affected by landslide activity. It combines heuristic, empirical and probabilistic approaches. This proposed methodology is based on the geospatial semantic array programming paradigm and has been implemented on a catchment scale methodology using Geographic Information Systems (GIS) spatial analysis tools and GNU Octave. The integrated data-transformation model relies on a modular architecture, where the information flow among modules is constrained by semantic checks. In order to improve computational reproducibility, the geospatial data transformations implemented in Esri ArcGis are made available in the free software GRASS GIS. The proposed modelling architecture is flexible enough for future transdisciplinary scenario analysis to be more easily designed. In particular, the architecture might contribute as a novel component to simplify future integrated analyses of the potential impact of wildfires or vegetation types and distributions on sediment transport from water-induced landslides and erosion.

## Introduction

Hillslope processes can be envisaged as a cascade where surface erosion and mass movements are visible expressions of critical instabilities in a complex system of interacting processes that control the downslope movement of material [1] in [2]. Field observations, modelling simulations and experimental studies have shown that soil erosion can vary considerably due to the changes in soil properties, vegetation cover and topography occurring after a landslide (e.g., [3, 4, 5]). Following landslide events the changes in soil erosion rates can be strong enough to deliver significant cascading impacts on ecosystems, for example due to an increased sediment yield to a stream network. This may potentially be of ecological and economic relevance not only locally (possibly driving complex changes even at the landscape-scale [6, 7]) but also off-site, whenever ecosystem services are important for service benefit areas connected through service connecting areas [8] (e.g., stream networks).

As natural resources are intrinsically entangled in complex networks there is a growing awareness of the importance of these cascades. This, in turn is driving the development of integrated risk assessment and multi-purpose use optimization of different resources to develop appropriate management policies that can reliably model the potential influence of climate change on these process cascades, and assess the resultant economic and societal consequences.

Landslide events will result in changes in topography and vegetation cover which, in turn, will alter surface erosion rates and sediment yields. There are a number of relevant models that use an integrated approach to soil erosion and landslide processes, including SHETRAN (the name derived from Système Hydrologique Européen-TRANsport) [9], TOPOG (a physically-based, distributed parameter, catchment hydrological model) [10, 11], PSIAC (Pacific Southwest Inter-Agency Committee) [12] or SIBERIA (also known as the Willgoose Catchment Evolution Model) [13]. WEPP-SLIP (Water Erosion Prediction Project – Shallow Landslide Integrated Prediction) [3] is a model that explicitly considers post-failure sediment yield. This model integrates the physical basis of the WEPP model [14], with the infinite slope stability model of Skempton and DeLory [15]. WEPP-SLIP is able to consider the post-failure changes in soil erosion rate through the changes in topography and land cover.

Physically-based models use a dynamic hydrological approach and local terrain characteristics for estimating spatial and temporal landslide probability [16]. The main limits of physically-based models are that they are often optimized for small catchments and local conditions, and that these require in depth knowledge of local soil and climatological parameters [17]. Empirical methods are mainly based on the estimation of thresholds related to precipitation patterns which result in landslide occurrence [16]. This approach generally requires high temporal resolution rainfall data, which are not often available, and do not necessarily model the right processes. In addition, it is limited to being applicable to only the same conditions under which it was developed [18, 17]. However, there is still room to improve the modelling of the interactions of these processes, for example through assessments of the changes in surface area made more susceptible to soil erosion following landslide events.

To quantify the potential changes in soil erosion due to landslide occurrence it is necessary to know where and when on the slope a landslide initiates and how it evolves. This paper aims to present a new modelling approach for data-poor regions in an attempt to improve the estimation of sediment budgets derived from rainfall induced landsliding and soil erosion. A statistical approach is proposed that is based on incorporating the frequency-area landslide distribution model of Malamud et al. [19] within the framework of a spatially distributed empirical soil erosion model.

## Study Area

The study area (Fig.1) is situated in southern Italy in the Daunia Appennines of the Puglia region, within the municipal territory of Rocchetta Sant’Antonio. It covers an area of almost 10*km*. This area is highly susceptible to landslide activity [20, 21], with a consequent negative impact on the local economy [22]. The neighboring area, to the northwest of the Rocchetta Sant’Antonio territory, presents a landslide frequency exceeding 20 percent for the overall area [23, 24, 22, 25]. Soil erosion is also widespread and the severity is largely determined by the combination of tillage practices and the high erodibility of the clay-rich flysch units from which some of the local soils are derived [26]. Within the catchment it is possible to distinguish four major classes of land use (agricultural soils, woodland, pastures and grassland) and three dominant lithologies (limestone, sandstone and clay). Slope angles are on average approximately 10 degrees with peak slope angles rarely exceeding 25 to 30 degrees. An ephemeral drainage network is fed by precipitation during the autumn-winter period when some 600 to 750

^{2}*mm*of rainfall is common [22]. The area is characterized by a Mediterranean sub-humid climate.

## A New Architecture for Coupling of the Effects of Rainfall-Induced Shallow Landslides and Soil Erosion

### Geospatial Semantic Array Programming

Array programming is an approach for simplifying complex algorithm prototyping with an accurate and compact mathematical description. It originates as a means for reducing the gap between mathematical notation and its implementation within the model’s algorithms in a formalized and reproducible way. As stated by Iverson [27]: “the advantages of executability and universality found in programming languages can be effectively combined, in a single coherent language, with the advantages offered by mathematical notation.” Array programming has been used for building the architecture for our modelling approach. For mitigating the complexity of trans-disciplinary modelling and the inconsistencies between input data, parameters and output, semantic checks on the processed information and a modularization of the key parts of the model were introduced following the semantic array programming paradigm (SemAP) [28, 29, 30]. The proposed architecture (Fig. 2) exploits the geospatial capacities of GIS in order to estimate soil erosion yield (e-RUSLE model). In our approach we integrated SemAP and geospatial tools (ArcGis and GRASS GIS) through the Geospatial Semantic Array Programming paradigm (GeoSemAP). GeoSemAP exploits geospatial tools and Semantic Array Programming for splitting a complex data-transformation-model (D-TM) into logical blocks whose reliability can more easily be checked by applying geospatial and mathematical constraints.

Semantic checks are exemplified in the following paragraphs with the notation **::constraint::**. The semantic constraints were implemented within the code with a specialized module [31] of the Mastrave modelling library.

### Applied techniques

The pre- and post-failure soil loss rate was calculated by applying the low data demanding model e-RUSLE [32]. This model retains all the equations of its predecessor (RUSLE, [33]) and implements an extra factor to account for the effects of soil stoniness on soil erosion. Due to the flexibility of the modelling architecture that e-RUSLE is based on, it is possible to calibrate the model for application at different scales [32]. e-RUSLE was implemented using the ArcGIS software to first estimate the **::nonnegative::** **::matrix::** representing the soil erosion rates within the catchment without considering the influence of mass movement. The scripts applied for calculating the soil erosion losses can also be easily carried out using an Open Source Free Software such as GRASS GIS or Quantum GIS.

*∞*) algorithm of Tarboton [34] was first used to calculate the flow direction and then the flow length. Due to the geomorphological characteristics of the study area, a multiple-neighbor flow algorithm was required with the D∞ algorithm being one of the most suitable [35, 36, 37]. In GRASS GIS it is possible to apply a multiple-flow approach using the tool ’r.watershed’ [38]. The slope steepness factor was also slightly modified in comparison to the application of the e-RUSLE presented in Bosco et al. [32]. This was based on the Nearing’s [39] equation which performs best for higher slopes [40, 32]. However the Moore and Burch [41] formula is more appropriate for slopes lower than 12.73 degrees because it gives the correct limiting value of zero in absence of any steepness. A comparison of both formulas is presented in Fig. 3, where a close matching trend is observed between 0 and 12.73 degrees (or 0 – 0.22

*rad*). Consequently a merged formula can be obtained by using the Moore and Burch equation for slopes less than 12.73 degrees and then the Nearing formula for higher slopes. To calculate the slope steepness factor of the model, the tool r.slope.aspect [42] of GRASS can be used. The majority of the equations that e-RUSLE is based up have been applied using the ArcGis tool ‘Map Algebra’ that in GRASS corresponds to ‘r.mapcalc’ [43].

For quantifying the effect of size, position and number of landslides affecting this catchment the frequency-size distribution model proposed by Malamud et al

*.*[19] was adopted. They found that landslide data from three quite different locations around the world (Italy, Guatemala and the United States) could be described quite well with the inverse gamma distribution.

In (1),

*p*= probability density (

*km*), Γ is the gamma function,

^{-2}*A*= the landslide area (

_{L}*km*

^{2}),

*ϱ*(-) is a parameter which controls the power law decay for medium and large landslide areas,

*a*(

*km*

^{2}) determines the position of the maximum in the probability distribution and

*s*(

*km*

^{2}) is a parameter which fits the exponential decay behaviour for small landslide areas. Parameter values of

*ϱ=1.4, a=*1.28

*10*and

^{-3 }*km*^{2}*s*= -1.32

*10*were shown to provide a good fit to the measured data. A dataset of more than 400 reported landslides that affected the catchment in 2006 was made available and published by Dr. Janusz Wasowski of CNR-IRPI, Bari [22, 25]. For obtaining the landslide inventory, high resolution IKONOS satellite imagery was used. To make the interpretation easier, the satellite images were orthorectified and pansharpened. This dataset is not freely available but the IFFI (Inventario dei Fenomeni Franosi in Italia) database [44] is a valuable alternative to apply our modelling approach whenever enough data are available.

^{-4 }km^{2}Overall a reasonable correlation between the inverse-gamma distribution of Malamud et al. [19] with the above parameter values and the frequency-size distribution of the landslide database was found (Fig. 4). The fit is very good for landslide areas greater than or equal to the peak in the distribution. For smaller landslide areas to the left of the peak the agreement is not as good, though modifications to parameters

*a*and

*s*could be made to improve this section. However the distribution of Malamud et al. [19] and the parameter values they used were shown to work over a wide range of landslide sizes from various countries around the world. It was found that these same parameter values also provided a similar fit to the data from our field site, suggesting the possibility of universality in the parameter values and therefore removing the need for calibrating the distribution for local applications.

On this basis, we kept the original Malamud parameter values and wanted to see how well they would perform against data from the Rocchetta catchment. The data for the smaller landslides do have a greater degree of uncertainty since their collection could easily have led to either an over- or underestimation of number of landslides. This could occur through either medium landslides being classified as smaller due to being covered by larger landslides, or through the smaller landslides being covered by larger ones and therefore missed completely. The main point of this exercise was not to match exactly the landslide-area probability distribution, but to have a physically realistic distribution on which to base our modelling. To predict when and where a landslide will occur is one of the main challenges for calculating post-failure soil loss in data-poor regions. We exploited the correlation between the measured data and Malamud’s distribution using Monte Carlo simulation to analyze the effects of mass movements on soil erosion by water.

Assuming the validity of the proposed inverse-gamma function for calculating the probability distribution of landslide areas we implemented a simple script (based on SemAP) in MATLAB language. Starting from a

**::scalar_positive::**number to represent the number of landslides that occurred in the catchment, we then calculate the number of landslides

*δN*in the

_{L}(h)*h-th*class of landslides. Each class is a

**::categorical-interval::**

**,**which includes all the landslides with an area from

*A*to

_{L}(h)*A*The classes thus form a partition of

_{L}(h+1).**::contiguous_interval::**

*s*in

*[0, A*] whose values are found from:

_{L}(hmax)In order to evaluate the effect of the post-failure changes on the soil erosion rates in the catchment, we applied the Monte Carlo method twice. Once to randomly determine the location of a landslide and a second time to sample the Malamud distribution to assign its size. The Monte Carlo simulation was also implemented in the MATLAB language following the SemAP paradigm and exploiting the potentiality offered by the Mastrave Library [29] whose tools were largely used within the code.

To be more explicit: considering

*Y*as a random variable distributed according to a given probability distribution, it is possible to generate n pseudo-random instances

*Y*with the same distribution. This may be accomplished with a classical Monte Carlo extraction. Let us define

_{1},…, Y_{n }*f*(

*⋅*) as a certain function of

*Y*which is implemented, within the SemAP paradigm, as a D-TM transforming an instance of

*Y*into the desired output data. Suppose we are interested in computing the integral

*A*of

*f*(

*⋅*) over a given domain

*Ω*. This implies considering the probability density function

*π*(⋅) of

*Y*over

*Ω*:

Numerically, it is possible to approximately estimate

*A*by exploiting the

*n*Monte Carlo instances

*Y*as

_{1},…, Y_{n }where

*Y*is the

_{run}*run-th*instance of

*Y*corresponding to the

*run-th*Monte Carlo iteration. From the law of large numbers, if n⇒∞,Â

_{n}⇒A.

In our particular application,

*Â*is the average over n runs of simulated landslides; in each of them the total erosion by water

_{n}*f*(⋅) is computed for the particular array of landslides

*Y*. The

_{run}*n*arrays of simulated landslides are the basis for

*f*(⋅) to estimate the corresponding post-landslide soil erosion. Each landslide occurring in the

*run-th*simulation has an area distributed according to

*̄p*(⋅). This defines

*π*(⋅) as the probability density function with which each

*run-th*array of landslides is distributed.

The Monte Carlo simulation was iterated 1,000 times. For each of the iterations the post-failure changes in soil erosion were calculated and compared with the pre-failure estimates.

The

**::matrix::**representing the cover management factor of the e-RUSLE model was calculated using a 5-by-5 meter resolution land cover map of the study site, produced by CNR-IRPI of Bari using ASTER satellite multi-spectral imagery and published in [22]. The map is not freely available but the CLC [45] is a valid open access alternative. The post-failure changes in vegetation cover were used within the model for estimating the effect of mass movement on soil erosion. Because of the modular modelling architecture (Fig. 2), the module that calculates the pre-failure C factor can be used as a link among our model and other approaches for measuring different land disturbance effects, in order to measure their effects on soil erosion.

The post-failure vegetation cover results were only partially altered by the slow mass movements that characterize this catchment (see Fig. 1). As locally the slide surface may also remain unchanged, we introduced into the model a value representing the post-failure percentage of bare soil. By analyzing the landslide dataset, the available pictures, satellite images and accounting for all the information collected during a field survey carried out within the study area, the percentage of the post-failure bare soil cover was estimated to be not less than 20 percent of the landslide area. For each of the pixels of the modelled landslides in each of the 1,000 Monte Carlo iterations, the

**::scalar_positive::**

**::proportion::**of bare soil was therefore randomly determined in the range 0.2 – 1.

## Results and Discussion

Table 1 shows the results of the Monte Carlo simulations. We replaced the mean values obtained by applying equation 4, with the median, because it is more stable in that it is only marginally affected by extreme values. By analyzing the median on 1,000 simulations of the cumulated pre-failure and post-failure soil erosion, an increase of 20 percent of the total soil loss was estimated. The post-failure soil erosion rate in areas where landslides occurred is, on average, around 3.5 times the pre-failure value.A bootstrap analysis based on 10,000 runs was performed in order to assess uncertainty. The analysis of the changes in the rate of soil erosion due to landslide occurrence shows post-failure increases in soil loss of approximately 1,700 tons per year (bootstrap p ≤ 0.05). This corresponds to an increase of around 22 percent of the total soil erosion. We also analyzed the extension of the area affected by slope instability. The bootstrap analysis shows that in each simulation at least 76 hectares, corresponding to around 8.5 percent of the catchment, are affected by landslide activity (bootstrap p ≤ 0.05). By comparing this value with the area that presented slope instability in 2006 (around 55 hectares), the applied methodology seems to result in a slight overestimate. The graph in Figure 3 shows that Malamud’s distribution seems to underestimate the number of small landslides (< 300

*m*). Nevertheless, the probability density distribution for the Rocchetta landslides from 2006 is in line with those reported by Malamud et al. [19] for precipitation-triggered landslides that took place in Guatemala in 1998. The model is in its early developmental phase and fine tuning the fit of the Malamud distribution to small landslides should help to improve the model predictions. However, for better evaluating the limits or the robustness of the proposed inverse-gamma distribution or of a modified version, further data would be necessary.

^{2}The bootstrap analysis, with 10,000 runs, performed on the measured data (Fig. 4) shows the uncertainty associated with a single year landslide dataset is too high to extrapolate different parameter values. A more detailed analysis based on datasets covering a longer time interval would help to improve the applied methodology. An additional source of error contributing to the predictions, which needs further investigation, arises from the selection of the model for estimating soil erosion and its running with limited data: thus there is considerable scope for errors in prediction to be strongly linked to this simplification.

Because the capacity to estimate the changes in soil erosion from landslide activity is largely dependent on the quality of the available datasets, the applied methodology broadens the possibility of a quantitative assessment of these effects in data-poor regions. The obtained results, even considering a possible overestimation, confirm the important role of mass movements on soil erosion and the consequent necessity to better integrate these processes into soil erosion modelling.

## Conclusions

A new method for empirically estimating the importance and extent of landslides on soil erosion losses in data-poor regions has been developed. This has been achieved by sampling the frequency-size landslide distribution proposed by Malamud et al. [19], and stochastically distributing the landslide location across the catchment. Given the increasing threat of soil erosion all over the world and the implications this has on future food security and soil and water quality, an in-depth understanding of the rate and extent of soil erosion processes is crucial.

Each year, on average, 8.5- 10 percent of the catchment shows evidence of landslide activity that is responsible for a mean increase in the total soil erosion rate between 22 and 26 percent above the pre-failure estimate. These results confirm the potential importance of integrating the landslide contribution into soil erosion modelling. While this approach clearly has limitations the proposed approach can be seen as a first attempt to assess the landslide-erosion interaction in areas with limited data.

The proposed modelling approach is also suitable to be applied in applications having a wider spatial extent and to be potentially implemented in a transdisciplinary context. For example, the relevant effect of wildfires on soil erosion and landslide susceptibility [46, 47] could be modelled with a higher reliability integrating the proposed approach. As stated in de Rigo et al. [47], wildfires can considerably increase soil erosion by water and landslide susceptibility. The changes in landslide susceptibility may in turn affect soil erosion. In general, considering the modelling architecture (Fig. 2), if the module that calculates the pre-failure C factor value would provide the layer altered by a different disturbance (e.g., wildfires or outbreak of pests), the presented modelling architecture could be applied for estimating the indirect effect of these disturbances on soil erosion, provided a new landslide susceptibility map, that considers the altered vegetation cover, is produced.

Although the preliminary results are promising, further research is required before this method can be applied by the scientific community and relevant authorities with any level of confidence. Consideration of, and integrating within the model, post-failure changes in topography and soil characteristics (e.g., soil armouring [48]) is fundamental for increasing the predictive capacity of the model. Also a better estimation of the bare soil exposed within a landslide is fundamental for improving our model. It would also be worthwhile to improve the fit of the Malamud distribution to the data, but at present this is not possible due to the limited availability of measured data. In order to obtain more reliable results, and more robust estimates of the effects of landslides on soil and vegetation cover, it will be also necessary to focus attention on producing a less uncertain zonation of the spatial probability of the landslide susceptibility in areas characterized by low data availability [49].

### Acknowledgements

We would like to thank Dr. Tom Dijkstra for his valuable comments on the manuscript. We also would like to thank Dr. Janusz Wasowski and Dr. Caterina Lamanna for providing the landslide data and Dr. Wasowski for his fundamental support during fieldwork. This paper is published with the support of the Maieutike Research Initiative.

### Authors Bio

Claudio Bosco graduated in 2002 from the University of Milan with a degree in natural sciences. His more recent research activities are focused on natural hazards and their link with climate change, combining research into quantitative, robust modelling approaches with expert-driven understanding of environmental processes. His research interests also cover quantitative geomorphology, spatial analysis (GIS-based) and wildfire effects on soil degradation processes.

Graham Sander is professor of hydrology in the School of Civil and Building Engineering at Loughborough University, U.K. His research interests cover sediment transport and soil erosion modelling, shallow overland flow, unsaturated subsurface water and contaminant transport.

References

[1] W.J. van Asch, “Water erosion on slopes and landsliding in a Mediterranean landscape,” Ph.D dissertation, Utrecht Univ., 1980.

[2] P.H. Van Beek, “Assessment of the Influence of Changes in Landuse and Climate on Landslide Activity in a Mediterranean Environment,” Ph.D dissertation, Utrecht Univ., Utrecht, 2002

[3] A. Cochrane, and Acharya, G., “Changes in sediment delivery from hillslopes affected by shallow landslides and soil armouring,” *Journal of Hydrology (New Zealand)*, vol. 50, no. 1, pp. 5-18, 2011.

[4] Acharya, T.A. Cochrane, T. Davies, E. Bowman, “The influence of shallow landslides on sediment supply: A flume-based investigation using sandy soil,” *Engineering Geology*, vol. 109, pp. 161-169, 2009. doi: http://dx.doi.org/10.1016/j.enggeo.2009.06.008.

[5] Acharya, T. Cochrane, T. Davies, and E. Bowman, “Quantifying and modeling post-failure sediment yields from laboratory-scale soil erosion and shallow landslide experiments with silty loess,” *Geomorphology*, vol. 129, no. 1-2, pp. 49-58, 2011. doi: http://dx.doi.org/10.1016/j.geomorph.2011.01.012.

[6] M. Bakker, G. Govers, C. Kosmas, V. Vanacker, K. Oost, and M. Rounsevell, “Soil erosion as a driver of land-use change,” *Agriculture, Ecosystems & Environment*, vol. 105, no. 3, pp. 467-481, 2005. doi: http://dx.doi.org/10.1016/j.agee.2004.07.009.

[7] Geertsema and J.J. Pojar, “Influence of landslides on biophysical diversity – a perspective from british columbia,” *Geomorphology*, vol. 89, no. 1-2, pp. 55-69, 2007. doi: http://dx.doi.org/10.1016/j.geomorph.2006.07.019.

[8] -U. Syrbe and U. Walz, “Spatial indicators for the assessment of ecosystem services: Providing, benefiting and connecting areas and landscape metrics,” *Ecological Indicators*, vol. 21, pp. 80-88, 2012. doi: http://dx.doi.org/10.1016/j.ecolind.2012.02.013.

[9] Ewen, G. Parkin, P.E. OConnell, “SHETRAN: Distributed River Basin Flow and Transport Modeling System,” *Journal of Hydrologic Engineering*, vol. 5, no. 3, pp. 250-258, 2000. doi: http://dx.doi.org/10.1061/(ASCE)1084-0699(2000)5:3(250).

[10] E. M. O’Loughlin, “Prediction of surface saturation zones in natural catchments by topographic analysis,” Water Resour. Res., vol. 22, no. 5, pp. 794-804, 1986. doi: http://dx.doi.org/10.1029/WR022i005p00794.

[ 11] CSIRO. (2010). *TOPOG Home page*. [Online]. Available: http://www.per.clw.csiro.au/topog/intro/intro.html (archived at: http://www.webcitation.org/6UNGXlsOl)

[12] Pacific Southwest Inter-Agency Committee, “Report of the water management subcommittee on factors affecting sediment yield in the Pacific southwest area and selection and evaluation of measures for reduction of erosion and sediment yield,” 1968.

[13] Willgoose, and S.R. Riley, “Application of a catchment evolution model to the prediction of long-term erosion on the spoil heap at Ranger uranium mine: Initial analysis,” Supervising Scientist Report, 132, Supervising Scientist, Canberra, 1998. [Online]. Available: http://www.environment.gov.au/resource/application-catchment-evolution-model-production-long-term-erosion-spoil-heap-ranger (archived at: http://www.webcitation.org/6UNHHqeOk).

[14] J. M. Laflen, L.J. Lane, G.R. Foster, “WEPP: A new generation of erosion prediction technology,” *Journal of Soil and Water Conservation*, vol. 46, pp. 34-38, 1991.

[15] A. W. Skempton, and F.A. DeLory, “Stability of natural slopes in London clay,” *ASCE Journal*, vol. 2, pp. 378-381, 1957.

[16] P. Jaiswal, and C.J. van Westen, “Estimating temporal probability for landslide initiation along transportation routes based on rainfall thresholds,” *Geomorphology*, vol. 112, pp. 96-105, 2009. doi: http://dx.doi.org/10.1016/j.geomorph.2009.05.008.

[17] J. de Vente, J. Poesen, G. Verstraeten, G. Govers, M. Vanmaercke, A. Van Rompaey, M. Arabkhedri, and C. Boix-Fayos, “Predicting soil erosion and sediment yield at regional scales: Where do we stand? ,” *Earth-Sci. Rev.*, vol. 127, pp. 16-29, 2013. doi: http://dx.doi.org/10.1016/j.earscirev.2013.08.014.

[18] R. Hessel, “Modelling soil erosion in a small catchment on the Chinese Loess Plateau,” Ph.D. dissertation, Utrecht Univ., Netherlands, 307 pp., 2002.

[19] B.D. Malamud, D.L. Turcotte, F. Guzzetti, and P. Reichenbach, “Landslide inventories and their statistical properties,” Earth Surface Processes and Landforms, vol. 29, no. 6, pp. 687-711, 2004. doi: http://dx.doi.org/10.1002/esp.1064.

[20] G. Iovine, M. Parise, E. Crescenzi, “Analisi della franosita’ nel settore centrale dell’Appennino Dauno,” Memorie della societa’ Geologica Italiana, vol. 51, pp. 633-641, 1996.

[21] P. Magliulo, A. Di Lisio, F. Russo, and A. Zelano, “Geomorphology and landslide susceptibility assessment using GIS and bivariate statistics: A case study in southernItaly,” Natural Hazards, vol. 47, no. 3, pp. 411-435, 2008. doi: http://dx.doi.org/10.1007/s11069-008-9230-x.

[22] J. Wasowski, C. Lamanna, and D. Casarano, “Influence of land-use change and precipitation patterns on landslide activity in the Daunia Appennines, Italy,” Quarterly Journal of Engineering Geology and Hydrogeology, vol. 43, no. 4, pp. 387-401, 2010. doi: http://dx.doi.org/10.1144/1470-9236/08-101.

[23] Mossa S., Capolongo D., Pennetta L., Wasowski J., “A GIS-based assessment of landsliding in the Daunia Apennines, southern Italy,” in Proceedings of the International Conference ‘Mass Movement Hazard in Various Environments’, M. Graniczny, M. Czarnogorska, *et al.*, Eds. Warsaw, Poland: Polish Geological Institute, 2005, pp. 86-91.

[24] J. Wasowski, D. Casarano, and C. Lamanna, “Is the current landslide activity in the Daunia region (Italy) controlled by climate or land use change,” in Proceedings of the International Conference on Landslides and Climate Change, R. McInnes, J. Jakeways, H. Fairbank, E. Mathie, Eds. London, UK: Taylor & Francis, 2007, pp. 41-49.

[25] J. Wasowski, C. Lamanna, G. Gigante, and D. Casarano, “High resolution satellite imagery analysis for inferring surface-subsurface water relationship in unstable slopes,” Remote sensing of Environment, vol. 124, pp. 135-148, 2012.

[26] C. Lamanna, D. Casarano, G. Gigante, and J. Wasowski, “Mappatura e studio dei fenomeni franosi nel Subappennino dauno con immagini satellitari ad alta risoluzione,” in Proc. 13a conferenza Nazionale ASITA, Bari, 2009.

[27] K.E. Iverson, “Notation as a tool of thought,” Commun. ACM, vol. 23, pp. 444-465, 1980. doi: http://dx.doi.org/10.1145/358896.358899.

[28] D. de Rigo, P. Corti, G. Caudullo, D. McInerney, M. Di Leo, and J. San-Miguel-Ayanz, “Toward open science at the European scale: Geospatial Semantic Array Programming for integrated environmental modelling,” Geophys. Res. Abstr., vol. 15, 13245+, 2013. doi: http://dx.doi.org/10.6084/m9.figshare.155703.

[29] D. de Rigo, “Semantic array programming for environmental modelling: Application of the Mastrave library,” in International Environmental Modelling and Software Society (iEMSs) 2012 International Congress on Environmental Modelling and Software. Managing Resources of a Limited Planet: Pathways and Visions under Uncertainty, Sixth Biennial Meeting, R. Seppelt, A.A. Voinov, S. Lange, and D. Bankamp, Eds., 2012, pp. 1167-1176. [Online]. Available: http://mastrave.org/bib/de_Rigo_iEMSs2012.pdf.

[30] D. de Rigo, Semantic Array Programming with Mastrave – Introduction to Semantic Computational Modelling, The Mastrave Project, 2012. [Online]. Available: http://mastrave.org/doc/MTV-1.012-1/.

[31] D. de Rigo, “Applying semantic constraints to array programming: the module “check_is” of the Mastrave modelling library,” in: Semantic Array Programming with Mastrave – Introduction to Semantic Computational Modelling, 2012. Available: http://mastrave.org/doc/mtv_m/check_is.

[32] C. Bosco, D. de Rigo, J. Poesen, O. Dewitte, and P. Panagos, “Modelling Soil Erosion at European Scale: Towards Harmonization and Reproducibility,” Nat. Hazards Earth Syst. Sci. Discuss, vol. 2, pp. 2639-2680, 2014. doi: http://dx.doi.org/10.5194/nhessd-2-2639-2014.

[33] K. G. Renard, G.R. Foster, G.A. Weesies, D.K. McCool, and D.C. Yoder, “Predicting Soil Erosion by Water: A Guide to Conservation Planning with the Revised Universal Soil Loss Equation (RUSLE),” USDA Agr. Handb., no. 703, 1997.

[34] D.G. Tarboton, “A new method for the determination of flow directions and upslope areas in grid digital elevation models,” Water Resources Research, vol.33, no.2, pp. 309-319, 1997.

[35] S. Gruber, and S. Peckham, “Land-Surface Parameters and Objects in Hydrology” in Geomorphometry: Concepts, Software, Applications, T. Hengl and H.I. Reuter, Eds., Developments in Soil Science, Elsevier, vol.33, 2009, pp. 171-194. doi: http://dx.doi.org/10.1016/S0166-2481(08)00007-X.

[36] G.B. Chirico, A.W. Western, R.B. Grayson, and G. Bloschl, “On the definition of the flow width for calculating specific catchment area patterns from gridded elevation data,” Hydrol.Process., vol. 19, pp. 2539-2556, 2005. doi: http://dx.doi.org/10.1002/hyp.5730.

[37] R.H. Erskine, T.R. Green, J.A. Ramirez, L.H. MacDonald, “Comparison of grid-based algorithms for computing upslope contributing area,” Water Resour. Res., vol. 42, W09416, 2006. doi: http://dx.doi.org/10.1029/2005WR00464.

[38] C. Ehlschlaeger, “GRASS GIS manual: r.watershed,” in: GRASS Development Team, 2014. GRASS GIS 7.0.0svn Reference Manual. Open Source Geospatial Foundation, 2006. [Online]. Available: http://grass.osgeo.org/grass65/manuals/r.watershed.html .

[39] M.A. Nearing, “A single, continuous function for slope steepness influence on soil loss,” Soil Science Society of America Journal, vol. 61, no. 3, pp. 917-919, 1997. [Online]. Available: http://handle.nal.usda.gov/10113/6603.

[40] C. Bosco, E. Rusco, L. Montanarella, and S. Oliveri, “Soil erosion risk assessment in the alpine area according to the IPCC scenarios,” in Threats to Soil Quality in Europe, G. Toth, L. Montanarella, and E. Rusco, Eds., EUR 23438 EN, pp. 47-58, 2008.

[41] I. Moore, and G. Burch, “Physical basis of the length-slope factor in the universal soil loss equation,” Soil Science Society of America Journal, vol. 50, no. 5, pp. 1294-1298, 1986. doi: http://dx.doi.org/10.2136/sssaj1986.03615995005000050042x.

[42] M. Shapiro, and O. Waupotitsch, “GRASS GIS manual: r.slope.aspect,” in: GRASS Development Team, 2014. GRASS GIS 7.0.0svn Reference Manual. Open Source Geospatial Foundation, 2006. [Online]. Available: http://grass.osgeo.org/grass65/manuals/r.slope.aspect.html .

[43] M. Shapiro, and G. Clements, “GRASS GIS manual: r.mapcalc,” in: GRASS Development Team, 2014. GRASS GIS 7.0.0svn Reference Manual. Open Source Geospatial Foundation, 2006. [Online]. Available: http://grass.osgeo.org/grass64/manuals/r.mapcalc.html.

[44] V. Agnesi, L. Arziello, P. Aucelli, A. Baglioni, C. Bettucci, *et al.*, “Rapporto sulle frane in Italia. Il Progetto IFFI – Metodologia, risultati e rapporti regionali, APAT, Rapporti 782007,” p. 681, 2007. [Online]. Available: http://www.progettoiffi.isprambiente.it/cartanetiffi/carto3.asp?cat=40\&lang=IT#

[45] European Environment Agency, “Corine Land Cover 2006 raster data – version 15 (08/2011),” 2011. [Online]. Available: http://www.eea.europa.eu/data-and-maps/data/corine-land-cover-2006-raster-1

[46] Di Leo, D. de Rigo, D. Rodriguez-Aseretto, C. Bosco, T. Petroliagkis, A. Camia, and J. San-Miguel-Ayanz, “Dynamic data driven ensemble for wildfire behaviour assessment: A case study,” IFIP Adv. Inf. Commun. Technol., vol. 413, pp. 11-22, 2013. doi: http://dx.doi.org/10.1007/978-3-642-41151-9_2.

[47] D. de Rigo, D. Rodriguez-Aseretto, C. Bosco, M. Di Leo, and J. San-Miguel-Ayanz, “An architecture for adaptive robust modelling of wildfire behaviour under deep uncertainty,” IFIP Adv. Inf. Commun. Technol., vol. 413, pp. 367-380, 2013. doi: http://dx.doi.org/10.1007/978-3-642-41151-9_35.

[48] G. Acharya, and T.A. Cochrane, “Rainfall induced shallow landslides on sandy soil and impacts on sediment discharge: A flume based investigation,” in The 12th International Conference of International Association for Computer Methods and Advances in Geomechanics (IACMAG), Goa, India, 2008. [Online]. Available: http://ir.canterbury.ac.nz/bitstream/10092/3134/1/12612367_Paper-IACMAG2008.pdf.

[49] C. Bosco, D. de Rigo, T. Dijkstra, G. Sander, and J. Wasowski, “Multi-Scale robust modelling of landslide susceptibility: Regional rapid assessment and catchment robust fuzzy ensemble,” IFIP Adv. Inf. Commun. Technol., vol. 413, pp. 321-335, 2013. doi: http://dx.doi.org/10.1007/978-3-642-41151-9_31