Steep slopes mantled by pyroclastic deposits are favorable areas prone to generate hazardous volcaniclastic flows. In Italy, such a setting is well represented in the Campania Region, where pyroclastic deposits from the explosive activity of the Neapolitan volcanoes (Ischia, Campi Flegrei, and Somma-Vesuvius) cover the Apennine range bordering the Campanian Plain. In order to provide a useful contribution to the mitigation and prevention of these calamitous natural events, this work presents a multidisciplinary approach to improve the understanding of the volcaniclastic flow hazard zonation in an Apennine area of 340 km2 surrounding the Somma-Vesuvius volcano. The disruption proneness index (DPI) was calculated in order to identify the drainage basins potentially prone to generate volcaniclastic flows. This index is obtained by combining satellite and morphometric data in a geographic information system (GIS) environment. It is calculated for 1100 drainage basins, considering the main parameters influencing the slope stability (slope angle, basin shape factor, curvature, relative relief, aspect, and land cover). The land cover mapping is obtained from Landsat data and airborne high-resolution images, while the morphometric parameters are derived from a digital elevation model (DEM) with a cell size of 10 m. The result is a zonation map that classifies the drainage basins according to different degrees of proneness to generate volcaniclastic flows (low, moderate, high, and very high). The drainage basins falling within high and very high classes are 66%, while 28% fall in the moderate class, and the remaining 6% fall in the low proneness class.
The volcanic areas affected by pyroclastic deposits and significant hillslopes are considered zones with high proneness for triggering volcaniclastic flows. Volcaniclastic flows are variously saturated mixtures of water and sediment that flow along open slopes or drainage channels driven by the force of gravity (Smith and Fritz, 1989; Smith and Lowe, 1991). Volcaniclastic flows that occur on volcano slopes are usually termed “lahars” (Smith and Lowe, 1991; Vallance, 2000). One of the most common causes of volcaniclastic flows is the combination of heavy rainfall, steep slopes, and loose substrata (Pareschi et al., 2000; Scott et al., 2001; Pareschi et al., 2002; Zanchetta et al., 2004a, 2004b; Sulpizio et al., 2006). Other volcaniclastic flow triggers include volcanic eruptions where pyroclastic density currents cause melting of snow and ice at the volcano summit (e.g., Nevado del Ruiz, Colombia, 1985; Pierson et al., 1990; Pierson, 1995); sector collapse, debris avalanches, and landslides (e.g., Casita volcano, Nicaragua, 1998; Kerle et al., 2003), crater-lake breakouts that produce remobilization of existing deposits (e.g., Mount Pinatubo, Philippines, 1994; Newhall and Punongbayan, 1996), and earthquakes (Scott et al., 2001).
On the Italian peninsula, one of the areas most prone to generate volcaniclastic flows is the Campanian Apennine relieves bordering the southern Campanian Plain (Fig. 1), where alternating poorly to well-preserved fallout beds, volcaniclastic colluvium, and buried soils with variable amounts of bedrock material cover that mantle the carbonate bedrock (e.g., Di Crescenzo and Santo, 1999; Calcaterra et al., 1999, 2000; Zanchetta et al., 2004a, 2004b; Sulpizio et al., 2006). In particular, the fallout deposits are from explosive activity of the Somma-Vesuvius and Campi Flegrei volcanoes that occurred in the past 40 k.y. (Andronico et al., 1995; Cioni et al., 2003), which were dispersed mainly to the east and northeast by the prevailing stratospheric winds (e.g., Barberi et al., 1990; Cioni et al., 2003; Costa et al., 2009).
Volcaniclastic flows may be disastrous in terms of damage to infrastructures, private property, and human life (Lowe et al., 1986; Behncke et al., 2012). In the Campania Region (Fig. 1), the most calamitous volcaniclastic flow event occurred on 5–6 May 1998, causing the tragic death of more than 150 people and damage to several villages (Migale and Milone, 1998; Zanchetta et al., 2004a). The largest flows followed nine consecutive rainy days, with uninterrupted rainfall that lasted for 30 h. The slope failures were mainly concentrated in the vicinity of the Sarno Mountains, Cancello Hills (Fig. 1). The lack of rain gauges at the head of the drainage basins prevented precise correlation between rainfall intensity and the triggering of volcaniclastic flows. However, statistical analyses show that the total rainfall of 4–5 May (100–180 mm), although well below the maximum-recorded annual values, was seasonally exceptional (Pareschi et al., 2000; Zanchetta et al., 2004a). Many flows occurred simultaneously affecting villages of Sarno, Quindici, Episcopio, Siano, and Bracigliano (Figs. 2A–2C) disrupting buildings, bridges, and roads connecting these villages (Cascini, 2004).
The 1998 flows represent only the latest episode of these catastrophic events in the Apennine belt bordering the southern Campanian Plain, and it is worth noting they occurred in a period of volcanic quiescence (the last eruption occurred in A.D. 1944 from Somma-Vesuvius; Santacroce et al., 2008). During the previous two centuries in the Sarno area ∼40 flow episodes were recorded (Migale and Milone, 1998). Other relevant volcaniclastic flows occurred in 1954 at Villaggio di Vettica (30 casualties), in 1924 at Amalfi (30 casualties), in 1910 at Regina Major (a suburb of Maiori village) and Cetara (60 and 110 casualties, respectively), in 1841 at Mulino delle Capre (120 casualties), in 1823 at Salerno-Positano-Siano-Tramonti (120 casualties), and in 1764 at Gragnano (43 casualties) (Di Crescenzo and Santo, 2005; Bisson et al., 2007a and references therein).
These types of events are the result of complex interactions among several factors such as meteorology, geomorphology, and lithology that influence the instability of drainage basins (Pareschi et al., 2000). In a previous work (Bisson et al., 2013), 1296 drainage basins were recognized as potential source areas of volcaniclastic flows in the Sub-Apennine Vesuvian area.
In this paper, the existing classifications of the areas more prone to trigger volcaniclastic flows (Bisson et al., 2010, 2013) were improved, resulting in a more accurate hazard map. This objective was achieved by applying a multidisciplinary approach that takes into account landslide susceptibility evaluation methods (Guzzetti at al., 2006). These methods are generally based on the study of the factors that influence the slope instability, defined as causal factors (Anbalagan, 1992; Turrini and Visintainer, 1998). In the literature, the causal factors are combined using statistical, heuristic, and deterministic methods (Di Crescenzo et al., 2008; Bijukchhen et al., 2009; Vergari et al., 2011; Ghosh et al., 2012). In the heuristic method, the causal factors are weighted subjectively; the deterministic method is based on the physical laws driving landslides; the statistical approach, instead, is founded on the multivariate relationship between causal factors and past and present landslide occurrence (Vergari et al., 2011 and references therein).
Taking into account the large number of literature works about the volcaniclastic flows in the Campania Region (Migale and Milone, 1998; Pareschi et al., 2000, 2002; Calcaterra et al., 2000, 2003; De Vita and Piscopo, 2002; Calcaterra and Santo, 2004; Zanchetta et al., 2004a; Sulpizio et al., 2006; Bisson et al., 2013) and the type of available data in the investigated area, the heuristic method was here adopted to map the hazard by volcaniclastic flow triggering. Following this method, the parameters calculated for each drainage basin were weighted according to their influence on slope instability. The parameters were identified through DEM analyses and remote sensing elaborations using a GIS platform. The result was achieved in terms of disruption proneness index mapping. The GIS platform used to achieve this objective was ArcGIS 9.3 (ESRI© Environment). All data were geocoded to the WGS84 UTM Zone 33 reference coordinate system.
STUDY AREA AND GEOLOGICAL LAYOUT
The investigated area is located in the eastern portion of southern Campanian Plain and is represented by the Apennine hillslopes surrounding the Somma-Vesuvius volcanic complex (Fig. 1). This area extends from Cancello hills to the southwestern termination of the Sorrentina Peninsula (Fig. 1), covering ∼340 km2. The boundary of the study area was defined by previous identification of 1296 drainage basins (Bisson et al., 2013) where the pyroclastic deposits that alternated to colluvial material reach thicknesses sufficient to trigger volcaniclastic flows during heavy or persistent rainfall (Mazzarini and Bisson, 2008). The pyroclastic deposits covering the hillslopes derive from the explosive eruptions of the Campi Flegrei and Somma-Vesuvius volcanoes that occurred during the past 22 k.y. (Rosi and Sbrana, 1987; Santacroce, 1987) and are dispersed mainly to the east and northeast in accordance with the prevailing stratospheric winds (e.g., Barberi et al., 1990; Cioni et al., 2003; Costa et al., 2009).
In detail, the bedrock of the investigated area consists of a Meso-Cenozoic substratum represented by carbonate sedimentary sequences belonging to a paleogeographic unit deformed by Miocene orogenic phases. At present, the sedimentary sequences constitute the carbonate mountains (Lattari, Salerno, and Sarno Mountains) surrounding Neapolitan volcanoes (CARG Project). Pleistocene extensional tectonic activity is responsible for the present-day geological setting of the southern Campanian Plain, forming a semi-graben structure (De Vita and Piscopo, 2002). In the central part of this structural depression, the volcanic activity began from different eruptive centers: Ischia Island (150 k.y. B.P. ÷ A.D. 1302; De Vita and Piscopo, 2002; Brown et al., 2008), Campi Flegrei (<100 k.y. B.P. ÷ A.D. 1538; Rosi and Sbrana, 1987; Orsi et al., 1996; Di Vito et al., 1999) and Somma-Vesuvius (<39 k.y. B.P. ÷ A.D. 1944; Santacroce, 1987; Santacroce et al., 2008). The structural depression was filled by alluvial and volcanic deposits and now corresponds to the Campanian Plain (De Vita and Piscopo, 2002).
The values of the first five parameters (S, BSF, C, RR, and A) were derived from a DEM with a spatial resolution of 10 m (Tarquini et al., 2007), whereas the land cover (LC) was obtained by elaborations of remote sensing data. All the parameters are defined and discussed in the next sections of the paper.
The morphometric parameters in Equation (2) are described below.
Slope of the Landforms
For each cell, the algorithm calculates the maximum rate of change in elevation over the distance between the cell and its eight neighbors and identifies the steepest downhill descent from the cell. Slope thresholds for triggering volcaniclastic flows in the Campania area fall in the range from 25° (minimum slope following the model of Takahashi, 2000) to 45° (limit of internal slope angle of granular deposits; Calcaterra et al., 1999, 2000; De Riso et al., 1999; Pareschi et al., 2000, 2002; Zanchetta et al., 2004a).
The BSF describes the geometry of the drainage basin. It is calculated using the ratio between the upper (Au) and lower (Al) areas of the basin with respect to the average elevation line (BSF = Au/Al). According to the results of analyses of the drainage basins affected by 1998 flows of Sarno-Quindici, those characterized by a BSF greater than 1 (Au>Al) are statistically more prone to generate flows (Pareschi et al., 2000). This can be explained by considering that the greater capacity to collect water in their upper portions and the conveyance of this water toward the central channel can increase erosion at the bottom of the channel, thus favoring the destabilization of the volcaniclastic cover and consequently flow triggering (Pareschi et al., 2000).
Like the slope parameters, the curvature is calculated by matrix data obtained from the DEM, permitting the ability to distinguish between areas having convex, flat, or concave surfaces (Shaw and Johnson, 1995). It is a dimensionless numeric value (positive, null, or negative) that indicates changes in slope obtained by applying the second derivative to the elevation matrix. In detail, the curvature value is calculated by subtracting the profile value (i.e., rate of change of slope measured along the slope direction ) from the planform value (measured transverse to the slope direction). The first value governs flow acceleration and deceleration and, thereby, sediment deposition and erosion, while the second value influences flow concentration or dispersal. In a matrix cell, the positive curvature indicates that the surface is upwardly convex, whereas the negative curvature corresponds to an upwardly concave surface. A curvature value equal to zero indicates a flat or planar surface.
The curvature influences mass-wasting processes by governing the distribution of soil water. Convex surfaces (e.g., ridges) tend to disperse groundwater, impede the formation of perched water tables, and suppress the development of high water pressures in the pores between soil particles that contribute to slope instability (Sidle et al., 1985). Concave surfaces (e.g., tributary valley), in contrast, tend to develop perched water tables and concentrate groundwater, surface water, sediment, and organic material. In these depressions, the forces promoting mass movement (e.g., gravity, soil pressure, and material weight) can exceed those resisting movement (e.g., soil shear strength and buoyancy effects of soil pore water). It is, therefore, not a coincidence that the majority of flows arise in concave surface forms such as incised channels and depressions upslope of channel heads (Shaw and Johnson, 1995).
The relative relief parameter was calculated as the difference between the maximum and the minimum elevation within each drainage basin (considered the unit area). This parameter is derived from the elevation matrix and provides a measurement of fluvial channel incision generating erosive action (Vergari et al., 2011). As relative relief increases, erosive activity also increases, thus inducing disruption proneness.
The aspect parameter refers to the horizontal direction of a slope face. The aspect was calculated by matrix data derived from the DEM where each cell stores the value of the angle that defines the facing to the cardinal points. This parameter has an indirect influence on slope instability. The southern exposed facing is subjected to a greater thermal diurnal cycling with respect to the northern ones. This is due to the duration of the slope face exposure to the sun’s radiation, since a long exposure increases the slope face temperature. In the boreal hemisphere, the northern exposed facing is often shaded, while the southern one receives more solar radiation for a specific surface area. This fact also influences the amount of humidity absorbed by the exposed facing. The amount of both thermal stress and humidity produces conditions more favorable to erosional activity (Sahi et al., 1977). Moreover, southern exposed facings have a lower vegetation density as compared to northern exposed ones (Bennie et al., 2006). The growth of vegetation roots helps keep terrain material in place (De Baets et al., 2006). In the same atmospheric conditions (rains and winds), the south and southeastern facings are the most exposed to thermal stress and humidity and have less vegetation. Due to the paucity of protective vegetation cover (Pimentel and Kounang, 1998), south and southeastern exposed facings are more prone to erosional activity and therefore are potentially more unstable.
Land Cover Parameter
The land cover parameter was used to determine the type of surface and its influence with respect to the triggering of volcaniclastic flows. Indeed, it is known that changes in land use are among the most important factors influencing the occurrence of rainfall-triggered landslides (Glade, 2003). In particular, vegetation can stabilize loose deposits (Gyssels and Poesen, 2003), while the build-up area interferes in the gravity equilibrium of loose deposits (Vrieling, 2006), and the barren soil is the type of land cover that is most exposed to rain and wind, which are the two primary causes of bare land erosional activity (Pimentel and Kounang, 1998; Cantón et al., 2001).
In order to be consistent with the weight assigned to the classes of morphometric parameters in Equation (2), the satellite multispectral images acquired before 5–6 May 1998 were selected and elaborated in order to classify the surfaces covering the drainage basins. The use of remote sensing offers the advantages of multi-temporal large-scale view data at an affordable cost. In this work, the Landsat, with its continuity space mission, made possible the use of the multispectral data pre-event acquired by Landsat 5 TM. The TM sensor (http://landsat.gsfc.nasa.gov) collects radiation in six spectral bands from visible to medium infrared with a spatial resolution of 30 m and in a thermal infrared band with a spatial resolution of 120 m (Table 1).
Therefore, in the 1296 drainage basins, three masks of land cover were obtained and converted in vector maps: vegetation, barren soil, and build-up areas. Moreover, in order to assure a more correct classification, the three maps were checked using the high-resolution (1 m) airborne images acquired by the State Agency for Interventions on the Agricultural Market (AIMA) during a flight that took place in 1997–1998 at altitudes below 2000 m above sea level (asl) (Bisson et al., 2007b; Fornaciai et al., 2008). In case of discordance between the Landsat land cover mask and the corresponding AIMA coverage due to the different spatial resolution between Landsat and AIMA images, the land cover class was assigned to Landsat pixel using the prevalence criterion applied to higher resolution images.
The land cover parameter was obtained for the entire study area in 1997 excluding an area <2.6% affected by meteorological clouds. The land cover classes are shown in Figure 6. The vegetation, build-up zones, and barren soils cover an area of 288.4 km2, 14.8 km2, and 27.8 km2, respectively. The vegetation represents the prevalent land cover (84.8%), the build-up areas are scarcely present (4.4%), and the barren soil cover corresponds to 8.2% of the entire study area.
Parameter and Class Weights Assignment
In order to calculate the disruption proneness index (DPI) to solve Equation (2), the six parameters, previously defined, were weighted and ranked into several classes. In Equation (1), the weights of morphometric parameters and of the relative classes were calculated through a multivariate analysis applied to the 37 drainage basins affected by the 1998 volcaniclastic flows event, which was considered as a reference case. The values of six parameters calculated for this reference case are reported in Table 2. The multivariate analysis used in this work consists of the identification of the main changes in slope (Fig. 7) to define the classes of each parameter and the assignment of relative weight. In detail, considering the drainage basins affected by the 1998 volcaniclastic flows, the weight of each class was calculated using the class intensity, which was used to calculate the number of drainage basins falling in each class, and then normalized to 10 (maximum weight assigned to class characterized by the highest number of drainage basins). The weight assigned to each class is reported in Table 3.
In Figure 7A, the S parameter shows four classes (S < 15°; 15° ≤ S < 20°; 20° ≤ S < 27°; and S ≥ 27°), and the relative weights are 0, 2, 3, and 10, respectively (Table 3). This figure demonstrates that these 37 drainage basins do not have S < 15° and that more than 85% of basins are characterized by S ≥ 27°. This value is in good agreement with the minimum threshold values for the volcaniclastic flow triggering defined by the literature in the range 25° to 28° (Calcaterra et al., 1999, 2000; De Riso et al., 1999; Pareschi et al., 2000, 2002).
In Figure 7B, the BSF parameter shows two main slope changes individuating the same three classes defined by Pareschi et al., 2000 (BSF < 0.6; 0.6 ≤ BSF < 1.2; BSF ≥ 1.2). These classes indicate three different basin geometric shapes—drop, rhombus, and funnel. The funnel and rhombus shapes are more prone to flow generation. The weights are 2, 10, and 9, respectively (Table 3).
In Figure 7C, the C parameter shows three classes of drainage basins—concave, convex, and flat surface. In these 37 drainage basins, only two have a flat surface, whereas 18 and 17 are characterized by concave and convex topography, respectively. The weights assigned to concave (C < –0.005), convex (C > 0.005), and flat (–0.005 ≤ C ≤ 0.005) curvature are 10, 9, and 1, respectively (Table 3).
In Figure 7D, the RR parameter shows three slope changes individuating four interval classes (RR < 200 m; 200 ≤ RR < 400 m; 400 ≤ RR < 600 m; and RR ≥ 600 m), the weights for which are 0, 1, 7, and 10, respectively.
In Figure 7E, the A parameter shows the number of drainage basins falling in each exposed facing (W, NE, E, SE, S, and SW). This analysis defines six classes, the weights of which are 1, 2, 3, 7, 9, and 10, respectively. The aspect of drainage basins exposed differently from the drainage basins of the reference case was weighted as zero value (Table 3). The figure demonstrates that most of drainage basins affected by historical volcaniclastic flows have a south-southwest exposed facing and correspond to 62% of the 37 drainage basins.
For the LC parameter, the classes and their weights were obtained from previous works in the literature (Sarkar and Kanungo, 2004; Di Crescenzo et al., 2008; Bijukchhen et al., 2009; Vergari et al., 2011) and normalized to 10 with the same criterion used for the morphometric parameters. The barren soil, the vegetation, and the build-up classes were weighted as 10, 6, and 3, respectively. In particular, in the case of drainage basins affected by more classes of LC, the weights were calculated by weighting the percentage of the area of each LC class in each drainage basin (Equation 2).
In order to assign the relative weights of all six parameters, a normalization to 1 was applied to the number of basins falling in the highest classes, except the LC parameter, for which the literature data indicate a weight halfway between 0 and 1. Using this method, the weights of six parameters (S, BSF, C, RR, A, and LC) are 0.86, 0.89, 0.48, 056, 0.62, and 0.50, respectively (Table 3). As expected, the slope and the basin shape factor are the most important parameters.
DISRUPTION PRONENESS INDEX (DPI) MAPPING
The disruption proneness map was obtained by applying the DPI (Equation 1) to each drainage basin excluding those affected by cloud cover in the elaborated Landsat 5 TM scene (15% of total drainage basins). The resulting map is shown in Figure 8, where 1100 drainage basins are classified into four different disruption proneness degrees: low (DPI ≤ 17.55), moderate (17.55 < DPI ≤ 23.49), high (23.49 < DPI ≤ 28.06), and very high (DPI > 28.06). These four classes were defined applying the natural break method to the values interval ranging from the minimum to maximum DPI (6.67 and 37.97, respectively). In the DPI resulting map, 390 (35.5%) and 337 (30.5%) drainage basins fall in the very high and high disruption proneness hazard classes, respectively; 306 (27.8%) fall in the moderate hazard class, and 67 (6.2%) fall in the low hazard class.
This map highlights the fact that more than 65% of all classified drainage basins (727 of 1100) have high and very high proneness for volcaniclastic flows generation. The high and very high drainage basins are mainly located in three different zones: the northern portion of the study area (Cancello Hills and the southeastern part of the Clanio Valley); Sarno-Quindici-Bracigliano-Siano areas; and a third area that involves almost the entire Peninsula of Sorrentina excluding the southwestern portion.
A more detailed analysis of very high and high DPI spatial distribution allowed the individuation of the municipalities more potentially affected by volcaniclastic flows hazard. The histogram in Figure 9 shows the 66 municipalities characterized by high and very high hazard areas. Sarno can be considered a reference municipality for the greatest number of historical flows events (81) recorded in the past 400 years (Bisson et al., 2007a). For this reason, the percentage of area exposed to high and very high hazard in Sarno (red color in Fig. 9) can be used as the minimum threshold value to individuate the municipalities with greater probability to volcaniclastic flows triggering. Accordingly, the municipalities characterized by significant percentages of highly hazardous area (≥20%) are 34. The highest exposed municipalities are located in Sorrentina Peninsula: Corbara, Furore, Lettere, and Minori (Figs. 8 and 9). In particular, Lettere is the most exposed to potential flow triggering for the largest percentage of area (almost 75%) falling in high and very high hazard classes (Fig. 9). In addition, the population density in the 66 municipalities (Fig. 10) shows that the highest values in the 34 most exposed municipalities (Fig. 8) belong to Minori (1167 inhabitants (inh)/km2), Amalfi (985 inh/km2), Castallemare di Stabia (3928 inh/km2), Casola di Napoli (1381 inh/km2), Siano (1090 inh/km2), Gragnano (1925 inh/km2), Nocera Superiore (1536 inh/km2), and Sarno (792 inh/km2).
The multidisciplinary approach used in this work allowed more accurate identification of the drainage basins more likely to be exposed to volcaniclastic flows generation in the Apennine belt bordering the southern Campanian Plain. One thousand one hundred drainage basins were classified using the disruption proneness index (DPI) calculated by combining the parameters individuated as the most representative for influencing volcaniclastic flows generation. In the DPI equation, the parameters were weighted according to multivariate analysis and previous literature works and validated against the data of the drainage basins that generated the 1998 flows of the Sarno-Quindici disaster, which was used as a reference case.
The resulting map classifies the 1100 drainage basins, with 66% falling in areas having very high and high proneness to volcaniclastic flows generation. In detail, all drainage basins affected by the 1998 Sarno-Quindici volcaniclastic flows fall in high or very high hazard areas, with the exception of three drainage basins that have not been classified due to the presence of clouds in the Landsat image. The DPI map highlights three zones that potentially may be more affected by volcaniclastic flows hazard in the northern (Cancello Hills and Clanio Valley), central (Sarno-Quindici-Bracigliano-Siano zone), and southern portions (Sorrentina Peninsula) of the study area. In particular, the Sorrentina Peninsula shows 20 municipalities classified as high and very high hazard areas, while the Lettere municipality is the most exposed area, with more than 75% of its area classified as high and very high DPI.
Recently, during an intense rainfall on 22 January 2014, the Gragnano municipality was affected by a flow that wounded a man and caused damage to the road network. This fact confirms the high disruption proneness of this municipality individuated in Figure 10 with 30% of its area exposed to high and very high hazard.
Since the drainage basins with high and very high DPI cover an area of 204 km2 (corresponding to 75% of total area of the classified drainage basins), they represent a significant portion of the investigated territory that should be taken into account for mitigation and prevention planning of calamitous natural events such as volcaniclastic flows.
In future works, the DPI map will be further improved by updating the land cover to the present through the analysis of the new Enhanced Thematic Mapper data that will be acquired by the recently launched Landsat 8 platform.
This work was partially funded by Istituto Nazionale di Geofisica e Vulcanologia (INGV)-DPC V1 Project 2012-2103 (UR 2-WPS2/3). The European Space Agency (ESA) is acknowledged for providing the Landsat data and in particular Dr. M. Marconcini, Dr. F. Sarti, and F. Costantini. Finally, a special thanks to Dr. F. Doumaz and Dr. L. Colini for GIS suggestions. We thank Dr. Karoly Nemeth and an anonymous reviewer for their useful comments and suggestions that greatly improved the manuscript.