Hydrological modelling is described as the art to reconstruct the architecture of a largely unknown system, where the complex interaction between water and soil and the exchange of mass and energy between the functional components of a basin are not observable (Savenije, 2009). In this context, evaluation of hydrologically damaged land can be conceived as an interdisciplinary research, where geospatial techniques can transfer theoretical knowledge in physics and mathematics toward applied disciplines such as bioclimatology, hydrogeomorphology and agriculture (Diodato and Ceccarelli, 2004).
Among the hydrological damaging events, water erosion is one of the most pervasive, but often the less visible form of land degradation. Its increase was related to prosperity decline in the Mediterranean Europe of the past (Bintliff, 2002). The impacts of erosion on soil productivity are particularly significant because soil losses are generally cumulative and, on a human timescale, permanent (Blaschke et al., 2000). Particularly at watershed scale, soil water erosion generates important environmental and socio-economic consequences such as: decrease in soil fertility and productivity, transformation of land into fallow land not suitable for reforestation, irreversible reduction in arable soil, increase in flooding events or diffuse pollution of river networks (Vezina et al., 2006).
Responding to the soil erosion hazards of the modern era, interest groups and policy-driving organizations worldwide fostered the development and use of techniques for landscape restoration. For instance, the European Commission adopted a communication on soil protection, endorsed by the Council of Ministers in June 2002 (Gobin et al., 2003). Erosion data collected in Southern Europe show that the amount of rainfall has a crucial effect on soil erosion. Generally, in hilly Mediterranean shrublands, increasing runoff and sediment loss were observed as rainfall decreased especially in regions with rainfall higher than 300-500 mm year–1 (Kosmas et al., 1997; Zuazo and Pleguezuelo, 2008). Consequently, rainfall aggressiveness – estimated by the erosivity index – has been often used to assess the soil erosion risk (Renschler et al., 1999; Le Bissonnais et al., 2002). In long-term soil loss evaluation, the universal soil loss equation (USLE; Wischmeier and Smith, 1978) is used. However, in USLE, as well as in its revised version RUSLE (Renard et al., 1997), the spatial variability of rainfall erosivity and soil erodibility are the major sources of uncertainty (Boardman, 1993; Wang et al., 2002).
In this paper, a geosaptial modeling of soils at basin scale is proposed to improve the evaluation of some RUSLE factors. This model was based on the potentiality of spatial analysis of geographic information system (GIS) framework and on geostatistical analysis. Many papers published in the last decade prove that GIS softwares are useful to assess soil erosion (Panagopoulos and Antunes, 2008; Terranova et al., 2009) even if its prediction remains difficult because of the different spatial scales characterizing the RUSLE factors and the date availability. Despite these difficulties, it remains necessary to define soil erosion and to map the erosional classes. This is the aim of the present work, for the Sele River basin (Campania, southern Italy) since this area is subjected to stormwater, especially in the agricultural and river-torrential areas (Camarasa Belmonte and Segura Beltrán, 2001; Diodato, 2005a), where it leads to erosional soil degradation processes, landslide and flash-flood events (Thornes and Alcantára- Ayala, 1998; Ramos and Mulligan, 2003).
The Sele River basin (3236 km2) is located in the western (i.e. Tyrrhenian) side of southern Italy (Figure 1) and includes a large alluvial plain, that derives from the Pliocene to Quaternary aggradation of a structural depression located along the western margin of the southern Apennines, known as Salerno Gulf-Sele Plain graben (Ippolito et al., 1973; Bartole et al., 1984; Barra et al, 1999). The plain is about 400 km2 wide and shows a triangular shape delimited seaward by a straight sandy coast, stretching between the towns of Salerno and Agropoli, encircled landwards by a range of calcareous mountains (Lattari, Picentini, Alburni, and Soprano-Sottano mounts) as well as by the siliciclastic Cilento hills with a mean elevation of about 1000 and 400 m a.s.l., respectively.
The eastern most portion of this structural embayment is characterized by hills formed by very thick clastic wedge, known as Eboli Conglomerate lower-middle Pleistocene in age (Cinque et al., 1988). The southernmost portion of the plain is characterized by the presence of thick successions of travertine deposits known as Travertini di Paestum (Violante and D’Argenio, 2000; Amato et al., 2009). Seawards and laterally to the travertines, there is a narrow coastal plain characterized by Tyrrhenian beach-dune ridges interfingered with lagoonal and fluvio-palustrine deposits (Brancaccio et al., 1986, 1988; Barra et al., 1998). A younger coastal sector, including a dunal sandy system elevated 1 to 5 m a.s.l., is interposed between the Tyrrhenian sandy-coastal ridge and the present shoreline.
The soil thickness in mountains and hill areas ranges from 50 to 150 cm; andosols, luvisols cambisols and calcisols are the prevailing soils types. The andosols are developed on pyroclastic fall deposits and reached this area from the neapolitan volcanoes. The luvisols are developed on alluvial fan deposits and the cambisols- calcisols on calcareous- marly or clayey deposits. The highest thickness (50 to >150 cm) characterizes the alluvial and coastal plains, where cambisols, luvisols and arenosols are prevailing.
The climate in Sele basin is of Mediterranean type, with important spatial variation of both erosive rainfall and temperature, according to the elevation and the distance from the coast (Diodato and Fagnano, 2010). The mean annual precipitation, measured at 62 raingauges distributed across the entire catchment, ranges from 700 to 2000 mm, with average 1180 mm and standard deviation 367 mm.
The 52% of the area is covered by natural vegetation, mainly represented by decidous trees (36%) and natural grasslands (9%). The agricultural areas (45%) are mainly represented by non irrigated arable crops (19%), mixed cropping systems (12%) and olive groves (7%). The urbanized areas cover only 2% of the total Sele River basin (SRB) area.
Methods and spatial modeling
The average annual soil loss (E in Mg ha–1) due to water eosion per unit area and per year (30-years mean upon 1957-1996 period) is estimated by RUSLE (Renard et al., 1997):
R is the rainfall-runoff erosivity factor [MJ mm (ha h)–1];
K is the soil erodibility factor [Mg ha h (ha MJ mm)–1];
LS2D is the two dimensional topographic factor, with L = length and S the slope steepness;
C is the cover and management factor;
P is the conservation support-practices factor.
The LS, C, and P values are dimensionless. The accurate modelling of RUSLE factors is the crucial condition to obtain reliable results (Renschler and Harbor, 2002), especially in basins with low-resolution data (Diodato, 2005b) and scarce pedological samplings. The potentiality of spatial analysis, available in GIS environment, was used to define the C and LS factors as grid data, while climate (R) and soil factors (K) have been spatialized by using the kriging interpolaztion method (Burrough, 2001). Kriging provide statistically unbiased estimates of spatialized values starting from a set of sample points. The basic idea of kriging is to estimate the unknown attribute value at the unsampled location s0 as a wighted average of the neighboring observations (si):
z is a vector of the observed data (erosivity and erodibility) selected in the neighborhood location;
λ is weight vector associated with the distance between so and si;
Z is the value at (s0) prediction location (Johnston et al., 2001).
By multipliying the grids of all RUSLE factors, the soil erosion map of Sele Basin was drawn.
The available time series of meteorological variables don’t record sub-hourly rainfall, limiting the application of USLE\RUSLE to obtain erosivity values. Diodato (2004), predicted annual erosivity from rainfall data of Southern Italy:
with r2 =0.867 significant at P=0.01, and where
R is average annual rainfall erosivity;
Pr is the average annual rainfall (mm);
d is the annual mean max daily rainfall (mm);
h is the annual one hour max rainfall (mm).
Soil texture and organic matter data of soil samples were used to estimate soil erodibility (K) by using the equation from Torri et al. (1997):
K is the RUSLE erodibility factor (Mg h MJ–1 mm–1) of the sampled locations;
C is the % of total clay;
OM is the % of organic matter;
Dg is a grain size parameter calculated by the equation of Shirazi et al. (1988):
fi is the fraction of particles in the range of diameters di and di–1 (mm). For each particle size class (clay, silt, sand), di is the max diameter (mm), di–1 is the min diameter and fi is the corresponding mass fraction. Considering that only the percentages of sand (S), loam (L) and clay (C) were available and that the factor log10(√didi–1) is -3.5 for C, - 2.0 for L and -0.5 for S, the following simplified formula is used (Borselli, 2006): Dg = 0.01 (-3.5C -2.0L -0.5S).
Slope and length factors
The slope (S) and length (L) factors represent the topographic constraints and estimate the effects of slope angle and slope length on the sheet and rill erosion. To capture the effect of flow convergence, the LS factor is computed with the upslope contributing area per unit lenght (Desmet and Govers, 1996; Mitasova et al., 1996; 1998) by using the equation of Mitasova et al. (1996):
A is the upslope contributing area;
b is the slope angle,
a0 and b0 are the length (22.1 m) and slope (0.09) standard RUSLE values;
m and n are parameters related to the prevailing type of flow and soil condition.
Using the hydrological extension of ArcGIS release 9.1, the LS factor is estimated from a digital elevation model (DEM) of Sele River basin with a grid size of 20 m. This spatial resolution minimizes the effect of grid size on the soil loss value assessment (Lee and Lee, 2006; Rojas et al., 2008).
To this aim, four main steps were followed to define the upslope contributing area: i) the depressionless DEM was assessed by filling the sinks recognized in the DEM (Hutchinson, 1989); ii) the depressionless DEM was used to determine the flow direction; iii) the flow direction was used to determine the flow accumulation; iv) a threshold value of 125 cells, based on the comparison of flow accumulation and stream network at scale 1:5.000 (Regione Campania, 1998), was applied to the flow accumulation grid to identify the area where runoff is active. The slope angle was calculated in degree by using the slope algorithm of ESRI-Spatial Analyst. The parameter m, according to Wischmeier and Smith (1978), was defined by using a grid with variable exponent of 0.2, 0.3, 0.4 and 0.5 for slope gradient <1.1, to 3.3, to 4.5 and >5%, respectively. The parameter n was defined as 1.3 (Moore et al., 1993).
Vegetation cover and land management factors
The C-factor value describes the effect of vegetation cover on rillinterrill erosion. Based on the Corine Land Use map of Campania on a 1:50.000 scale (SESIRCA, 2004), 27 land use classes were attributed by a specific C factor (Zanchi and Giordani, 1995; Angeli, 2004; Bakker et al., 2008; Märker et al., 2008). To urban areas and other land uses where soil is not present (i.e. bare rocks, water courses) we assigned a C factor value of 0 for a a priori assumption (Bakker et al., 2008). To coniferous, mixed and broad-leaved forest we assigned slightly higher values (from 0.001 to 0.003), since in the study area the broad leaved forest are mainly represented by chestnut trees that are deciduous and therefore reduce winter protection of soil in comparison with coniferous (Angeli, 2004). Pastures was assimilated to alfalfa stand (0.02), while natural grassland was assimilated to unmanaged grassland (0.05). Moors and heathland and Sclerophyllous vegetation were assimilated to shrublands, but because of missing data for distinguish dense from sparse shrubs, we used their C factor average value (0.05).
Complex cultivation patterns and Agro-forestry areas were assimilated to arable dense tree cover, while annual crops associated with permanent crops were assimilated to arable medium tree cover (0.20 and 0.25 respectively). For fruit trees and berry plantations, olive groves and vineyards, the C factor values 0.30, 0.30 and 0.45, respectively were used (Angeli, 2004; Märker et al., 2008), since they were adopted in Tuscany Region that has similar conditions to southern Italy.
Considering that in SRB arable land is a mixture (both in space and in time) of different crops, in some cases also both irrigated and not irrigated, a long-term average value was used (0.30 from Bakker et al., 2008). Land mainly occupied by agriculture, with significant areas of natural vegetation was assimilated to arable sparse tree cover (0.30), while it was assigned a higher value (0.36 from Angeli, 2004) to the sparsely vegetated areas. For beaches, dunes and sands we assigned the max value (1.00) of fallow (Bakker et al., 2008), because these lands are not protected by vegetation. C factor values were corrected considering the state of the canopy cover during the year (Table 1), by multiplying the C factors reported in Table 2 per the monthly rainfall erosivity index (Rm = ratio between monthly erosivity and yearly R). The average of the new monthly C factors allowed the computation of the yearly C factors corrected on the basis of the annual distribution of ground cover and rainfall erosivity (Morgan, 1995). Since the vegetation cover protects the soil from raindrop impacts, a high value of the adjusted C factor indicates low ground cover during heavy rain periods (Vezina et al., 2006). Spatial and temporal variability of the erosion control practice factor (P) were not considered in this study, due to insufficient data about changes in management through time.
Soil loss tolerance
A tolerable soil loss is the max annual amount of soil that can be removed before the long term natural soil productivity will be adversely affected. The impact of erosion on a given soil type, and hence the tolerance level, depends on the type and depth of soil. In this way, the tolerance level (T) for most soils in SRB reaches values of 10-15 Mg ha–1 y–1 for soil depth higher than 100 cm (Schertz, 1983). For economic reasons, soil loss in optimal farm management is often considered higher than the soil loss tolerance: soil loss tolerance for economic planning (TEP) = 2T (Shi et al., 2004).
Results and Discussion
Spatial pattern of erosivity and soil erodibility
Rainfall erosivity (R) over SRB ranges from 600 to 4000 MJ mm (ha h yr)–1, with mean and standard deviation of 2000±883 MJ mm (ha h yr)–1. Figure 2a shows the kriging map of annual long-term R-factor with a 0.5 km grid resolution (Diodato and Fagnano, 2010), indicating that the annual rainfall erosivity has a moderate spatial variability in SRB. The soil erodibility map (Figure 2b) was derived from a new approach used for improving the spatial variability estimates reported in Diodato et al. (2010). Soil erodibility increases from the mountainous areas, where Mesozoic carbonates largely crop out (Platano, Melandro and Calore sub-basins), to the Sele alluvial plain, where erodibility values are very high [>0,03 Mg h(MJ mm)–1].
Slope-lenght and crop grid maps
The highest LS factor values in the SRB characterize the zones with steeper and longer slopes (Figure 3a). By applying the C factor according to the different land uses of the study area (Table 2), also the landcover map was drawn (Figure 3b).
From this map, it is clear that the vegetation covers with the lowest protection of soil are the arable lands and the olive groves and vineyard, mainly concentrated in the plains and in the hilly areas.
Mapping and classifying sediment source areas
Although the C factor-map in Figure 3b shows the the highest part of SRB is covered by vegetation with low C values (forest ~50%), in Figure 4a it can be seen that SRB is a region where soil erosion by water is an actual problem. In this map, the amount of the soil losses per year from each grid-cell is calculated on the basis of the rainstorm events averaged upon a long-term annual period. Their values range from 0 to 1000 Mg ha–1 y–1 with mean rates of soil losses of 53±43(SD) Mg ha–1 y–1. Using the soil delivery ratio (SDR), as defined in Diodato et al., (2009), about the 80% of the eroded soils were trapped in the depressions and valleys during the transport process via drainage network. This is in agreement with simulation results of WEPP model of the National Soil Erosion Research Laboratory (USA) realized in the Diano Valley, (http://topsoil.nserl.purdue.edu/nserlweb/weppmain/) along a typical hill slope profile (Figure 4b1) .
The WEPP model predicted the erosional and depositional trends, along the selected profile (Figure 4b2). After the first phase of soil detachment, there was a strong increase of the erosion rate (until 90 Mg ha–1) along the higher slopes (33%). Afterwards, at about 300 m along slope-length, the material eroded upstream is settled down, and only a small part of it was conveyed in the main hydrographic network.
In order to depict a best visualization of areas characterized by critical erosion rates, a classification with a tolerable-and-no tolerable soil loss is linked to the legend of Figure 4a. The tolerance threshold for economic planning TEP = 20 Mg ha–1 y–1, for soil depth >100 cm is considered. In this way, 32% of the basin was subjected to no-tolerable soil losses, and about 7% was affected by catastrophic erosion (e.g., rates >80 Mg ha–1 y–1). The remaining 68% of the basin was protected by important erosive processes, however, about 9% of its territory still remains vulnerable because of the low soil profile depth, local shallow hilly-lands, lower weathering rates and eventual salinitasion that can exacerbate soil degradation.
Geovisualization of erosion processes
The foothills and the greatest part and hilly areas are affected by erosion rates exceeding the tolerable threshold of soil loss. Furthermore, comparing topography and land-cover (Figure 3a,b) with soil loss of Figure 5, it is evident that soil erosion rates exceed the acceptable thresholds on a wide range of landforms throughout the basin with medium-low slopes and medium-high C and K values. This is more easily perceptible from the perspective view of Figure 5, where the red areas Cilento siliciclastic hills located between the carbonate ridges of Mt. Alburno and Mt. Cervati, made of marls, sandstone, clay and conglomerate, to the clay and other terrigenous rocks of the upper Sele Valley as well as to the Pliocene-Pleistocene clastic successions of the Auletta basin (Gioia and Schiattarella, 2010). Wide areas of the Vallo di Diano basin and Sele plain, are classified with negligible soil loss (green areas) because of the very low slope. Increases of soil erosion, but still below TEP, typify the low-mountain areas, with estimated rates around 10 Mg ha–1 y–1. The results show that 30%, is subjected to uncontrolled erosion, while in about 10% of area, the erosion could be controlled with appropriate soil management practices.
These results are derived from a soil erosion map based on long-term average erosivity pattern. However, geomorphological processes can be dominated by few severe and stochastically rainstorms, that are not represented by RUSLE model. It is interesting to note that, although runoff, and the consequent sediment transport, experienced a common decrease over recent decades, extreme rainfalls and erosivity are growing in an erratic way, leading to a greater and chaotic mobilization of soil across the Sele basin, as reported by Diodato et al. (2009). This evidently creates not only socio-economical costs but hazardous environmental impacts because of repeated dispersion in rivers of chemical nutrient, previously lost with the soil erosion.
Modelling assumption and validation
As mentioned above, the spatial variability of rainfall erosivity (Diodato and Fagnano 2010) and soil erodibility (Diodato et al., 2010) are the major sources of uncertainty in RUSLE-GIS approach. Also the accuracy of DEM can affect up-slope contributing area and slope that are used to calculate LS factor. Mitasova et al. (1996) model the propagation of uncertainty in the prediction of I by interpolating a DEM to finer spatial resolution, suggested that a DEM at a spatial resolution coarser than 20 m is insufficient to calculate the up-slope contributing areas. Gertner et al. (2002) supported this assertion. Conversely, Awa et al. (1987) studied the effects of three spatial resolutions (0.25, 1 and 4 ha) on the prediction of the RUSLE using three non-parametric tests and did not find statistically significant differences in the predicted sediment loading due to cell size. Molnár and Julien (1998) instead showed that, although large grid sizes tend to underestimate soil losses in RUSLE model, a correction factor must be included only when DEM-grid sizes exceeds 100×100 m. Rojas et al. (2008), found that at grid sizes coarser than 150 m, the sediment source areas became less appropriately depicted by RUSLE and that therefore the best results to simulate soil erosion can be obtained at grid sizes smaller than 150 m. Therefore a cell resolution of 20 m adopted in this work could be considered appropriate for the geomorphological feature of the SRB.
The RUSLE model validation was not possible, because direct measurements of soil erosion were not available in the SRB. However, the estimates here calculated (53 Mg ha–1 y–1), are within the 95th confidence interval of 53-93 Mg ha–1 per year, found in the same basin by previous simulations made with the CliFEM-approach calibrated on a close basin (Diodato et al., 2009).
By integrating geographic information system and the revised universal soil loss equation, the soil erosion map of Sele River basin was drawn. The spatial and geostatistical analysis, available in GIS framework, allowed a spatial assessment of the erosion hazard over the study area.
The results suggest that soil erosion rates exceeding tolerable thresholds affect about 50% of SRB area involving different types of soils and land uses, that mainly correspond to the hilly areas located between the carbonate mountains bounding the study area.
In these areas, showing high-very high soil erosion hazard, a strategic program for a more detailed erosion mapping is necessary to define in advance the impacts of land use change. This could allow the identification of the agricultural management practices more suitable for soil protection in the different cropping systems, such as perennial cover crops in fruit trees, olive groves and vineyard or conservative agronomic techniques (i.e. sod-seeding, minimum tillage) in arable lands.