Mapping soil erodibility from composed data set in Sele River Basin, Italy

Nat Hazards DOI 10.1007/s11069-010-9679-2 ORIGINAL PAPER Mapping soil erodibility from composed data set in Sele River Basin, Italy Nazzareno Diodato • Massimo Fagnano • Ines Alberico Giovanni Battista Chirico • Received: 26 February 2010 / Accepted: 24 November 2010  Springer Science+Business Media B.V. 2010 Abstract Evaluation of soil erodibility is an important task for Mediterranean lands, in which fertility and crop yield are significantly affected by soil erosion. The soil physicochemical parameters affecting soil erodibility are highly variable in space and, as for many other environmental variables, sample measurements are generally not enough for assessing its spatial variability with an acceptable level of uncertainty at the scales of practical interest. This study illustrates the procedure applied for estimating the pattern of soil erodibility across the Sele Basin (Southern Italy), where soil properties have been measured on a limited number of sparse samples. Sampled data were integrated with other sparse data estimated by local regression functions, which relate soil erodibility to auxiliary variables, such as terrain attributes and land system class memberships. Sampled and estimated data were merged in a composed data set to assess the spatial pattern of soil erodibility by ordinary kriging. The proposed approach offers effective spatial predictions, and it is exportable to regions where financial costs for soil sampling are not feasible. Keywords Soil erodibility  Geospatial mapping  Topotransfer function  Southern Italy 1 Introduction The high variability of soil properties and their unpredictability from a deterministic perspective have led researchers to consider spatial variability as a key characteristic of N. Diodato MetEROBS—Met European Research Observatory, GEWEX-CEOP Network, World Climate Research Programme, 82100 Benevento, Italy M. Fagnano  G. B. Chirico (&) Dipartimento di Ingegneria Agraria e Agronomia del Territorio—University of Naples Federico II, Via Università 100, 80055 Portici, Italy e-mail: gchirico@unina.it I. Alberico CIRAM—Centro Interdipartimentale di Ricerca Ambiente—University of Naples Federico II, 80134 Naples, Italy 123 Nat Hazards soils to be exploited for generating prediction maps within a stochastic framework (e.g. Castrignanò and Lopez 2000). Predictive soil mapping became a feasible task after 1970s, when automatic spatial interpolation techniques were implemented within geocomputational tools (Burgess and Webster 1980; Scull et al. 2003). Soil-landscape mapping studies have generally focused on physicochemical, biogeophysical and hydraulic properties (e.g. Sinowski and Auerswald 1999; Diodato and Ceccarelli 2004; Hengl et al. 2004; Herbst et al. 2006; Yugang et al. 2007; Bourennane et al. 2007; Sigua and Hudnall 2008), while rarely on soil erodibility (e.g. Zucca et al. 1999; Veihe 2002; Ozcan et al. 2008; Castrignanò et al. 2008), although soil erodibility is a significant problem for environmental, social and economic fields. Soil erodibility is difficult to be determined, since it is not a directly measurable physical property. Soil erodibility represents the susceptibility of the soil to be eroded by detachment and transport processes triggered by erosive forces. Thus, soil erodibility is a synthetic description of soil erosion for given reference conditions, and it is influenced by several soil properties controlling different subprocesses (e.g. Salvador Sanchis et al. 2008). At subregional and local scales, soil erodibility is a very important parameter for soil erosion models aiming at the definition of best planning management practices for protecting fragile ecosystems (Park and Egbert 2005; Bayramin et al. 2008). For instance, an interesting study was carried out by Baskan and Dengiz (2008) to evaluate and compare the relationship between soil erodibility maps computed by traditional and geostatistical methods for a small basin in Turkey. However, upon large regions characterized by erosion problems, soil property data are generally too scarce for an effective direct evaluation of soil erodibility patterns. So that, different multivariate geostatistical procedures have been suggested to improve the knowledge of the spatial structure of the primary targeted variable, i.e. soil erodibility, by exploiting auxiliary variables consistently correlated with the former (e.g. Goovaerts 2000; Monestiez et al. 2001; Wang et al. 2003; Diodato 2005, 2006; Alison et al. 2005; Simbahan et al. 2006; Casa and Castrignanò 2008). A drawback of multivariate geostatistical techniques is that inference demands a high number of samples of primary and auxiliary variables at appropriate spatial scales in order to gather a detailed specification of spatial cross-covariance structure (e.g. Rodriguez-Iturbe et al. 1998). This is a common issue in several regions characterized by erosion problems, where primary and auxiliary variables of soil features affecting soil erodibility are scarce. The aim of this work is to propose a new procedure, based on both statistical analysis and expert knowledge, to define soil erodibility patterns of areas characterized by data scarcity, as in the Sele Basin (Southern Italy). 2 Material and method In this paper, we evaluate the soil erodibility according to the RUSLE erodibility K–factor, which is defined for a specific soil as the rate of soil loss per rainfall erosivity unit, in a clean-tilled fallow condition of a plot with 9% slope (Renard et al. 1997). It is one of the factors employed in the RUSLE equation to evaluate the average soil loss rate per year in the bare Wischmeier reference plot (Renard et al. 1997). The flow chart in Fig. 1 illustrates the procedure applied in this study for generating a soil erodibility map of the Sele Basin, in comparison with the application of ordinary kriging. The proposed procedure is similar to the one proposed by Abbaspour et al. (1998). Soil erodibility data (hereafter referred to as sampled data) calculated by directly measured soil properties in a limited number of sample locations are integrated with soil erodibility data estimated from 123 Nat Hazards Fig. 1 Flow chart of composed ordinary kriging procedure adopted in this work and its comparison with the ordinary kriging method site-specific regression equations (hereafter referred to as estimated data). These regression equations, hereafter referred to as Local Topotransfer Functions (LTFs), provide estimates of soil erodibility from other auxiliary variables, such as terrain attributes and specific information regarding land system class memberships, which show significant correlation with soil erodibility sampled data. The combination of estimated and sampled data generates a larger data set (hereafter referred to as composed data set) usable as input for composeddata ordinary kriging (CO_OK) interpolation method. This procedure can be used as an alternative to other types of multivariate geostatistical techniques in areas characterized by a limited number of sampled data, but where it is possible to obtain estimated data from other auxiliary variables. LTFs can be calibrated with primary and auxiliary data pairs available at same locations and they can be then applied to estimate values of the primary data at locations for which auxiliary data, but not primary data are available. As evidenced by Abbaspour et al. (1998), CO_OK, similarly to other geostatistical techniques, is also able to treat geostatistical parameters (i.e. mean, variance, nugget, range and shape of the semivariogram) as uncertain random variables and therefore allows analysis of parameter uncertainty, which is inherently associated with semivariogram modelling. 2.1 Study area The Sele Basin is located across the Southern Campania and Western Basilicata regions, in the western side of Southern Italy (Fig. 2a, b). The Sele basin covers an area of 3,236 km2, which includes large alluvial plains and a large sector of the Campania Apennines (highest elevation 1,899 m a.s.l.) characterized by a Mesozoic bedrock, mainly consisting of limestone and subordinate dolostones partly covered by pyroclastic products derived from the explosive eruptions of the Vesuvius. 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 2,000 mm, with average 123 Nat Hazards Fig. 2 Geographical setting (a), peninsula of South Italy (b), land-use map of Sele River Basin in bold line (c). Land-cover source in (c) map was arranged from Vector Layers Processing (Y. Farhat) http://hydis.eng.uci.edu/gwadi/ 1,180 mm and standard deviation 367 mm. The original forest covers are fragmented by agricultural lands, which characterize the land use mostly across the plains and the hills of this region (Fig. 2c). 2.2 Soil sampling and erodibility assessment A set of 114 soil samples was collected from the top layer (0–10 cm) and subjected to laboratory measurements to determine particle size distribution with hydrometer method (Bouyoucos 1962). Organic carbon in soil was determined with the dichromate method (Walkley and Black 1934), whereas organic matter content was calculated by multiplying the organic carbon content by 1.724. The RUSLE K-factor (Mg h MJ-1 mm-1) has been computed according to the equation suggested by Torri et al. (1997): K ¼ 0:0293ð0:65  Dg þ 0:24Dg2 Þ " #   OM OM 2 2  0:00037  exp 0:021 4:02CL þ 1:72CL CL CL ð1Þ where OM is the organic matter content expressed as a percentage, CL is the clay content expressed as a fraction, Dg is the decimal logarithm of the geometric mean of soil particle size, which is calculated as follows (Shirazi et al. 1988): pffiffiffiffiffiffiffiffiffiffiffiffi X fi  log di di1 Dg ¼ ð2Þ i where fi is the fraction of the texture classes (clay, silt and sand) with corresponding maximum di and minimum di-1 sizes. The main statistics of the measured soil properties and estimated RUSLE K-factor are illustrated in Table 1. Figure 3 shows the location of the sample data. The soil erodibility values computed according to Eqs. (1) and (2) are refereed to as sampled data, to emphasize that these erodibility values are directly derived by measured soil properties. We evaluated the possibility to apply an updated method for estimating soil erodibility, such as the one suggested by Borselli et al. (2009), using the KUERY software package, which has been shown to have better performance than the former method by Torri et al. (1997). However, we had to discard the option of using the KUERY software, 123 Nat Hazards Table 1 Main statistics of the soil samples in the study area OM (%) CL (-) Dg (log10 mm) K (Mg h MJ-1 mm-1) Average 2.52 0.28 -1.79 0.028 Median 1.81 0.27 -1.76 0.028 Standard deviation 1.94 0.10 0.36 0.005 Minimum 0.52 0.04 -2.61 0.010 Maximum 12.78 0.53 -0.93 0.044 OM organic matter content, CL clay content, Dg decimal logarithm of the geometric mean of soil particle size, K erodibility factor of the RUSLE equation Fig. 3 a Coarse soil erodibility map from JRC–European Soil Bureau, http://eusoils.jrc.ec.europa.eu/) with soil sampling locations (white circles); b shaded terrain map with supplementary locations (white circles) where soil erodibility has been estimated by Local Topotransfer Functions for generating a ‘‘composed data set’’ (black?white circles) as this model requires some data, such as rock content, which was not available for the entire data set across the whole study area. We preferred therefore to use the approach suggested formerly (Torri et al. 1997), which requires only data (organic matter and soil particle size distribution) which are traditionally measured in soil sample analysis executed for agricultural practices. Moreover, in the limited sample points where a complete data set was available, the erodibility values estimated with the approach proposed by Borselli et al. (2009) were not significantly correlated with the organic matter (R2 = 0.02; n.s.). This result was unexpected since, according to our experience, soil organic matter should contribute significantly in reducing soil erodibility. The procedure suggested by Torri et al. (1997), instead, gave erodibility values well correlated with soil organic matter (R2 = 0.56; P = 0.01) 2.3 Supplementary data and regression analyses When the number of sampled data is low for reliable spatial estimates of the primary variable, auxiliary variables such as terrain attributes and soil class membership can be 123 Nat Hazards used as potential predictors to estimate the primary variable in supplementary locations by using Local Topotransfer Functions (Salski 2006; Wessolek et al. 2008). A prerequisite for the development of these functions is the availability of auxiliary variables in the same locations, covering representative value ranges in which these variables are expected to vary within the study area (Illian et al. 2008). Similarly to the pedotransfer functions (Matula and Špongrová 2007), LTFs can be linear or non-linear regression equations. In this case study, we explored the possibility to use two types of auxiliary variables: • class memberships according to the land system classification of Campania Region edited by Di Gennaro (2002) following an integrated approach suggested by FAO (1995) and Dalal-Clayton and Dent (2001); • terrain attributes, such as elevation, terrain slope and aspect. Sampled data have been first classified according to the land system class membership (see Table 2). A correlation analysis has been performed between terrain attributes and sampled data, within each land system class. Terrain attributes are the environmental variables most commonly used as auxiliary variables of soil data, since topography is an important pedogenetic factor and it is strongly related with other pedogenetic factors (Jenny 1980). Elevation has been selected among other terrain attributes, being most suitable for inferring soil erodibility K-factor, based on a Pearson correlation test. The following LTFs have been identified by linear regression analyses: KeðE22Þ ¼ 0:02842 þ 4  105  ele; r ¼ 0:79; n ¼ 5 ðP ¼ 0:112Þ ð3Þ KeðE23Þ ¼ 0:02965 þ 4  105  ele; r ¼ 0:77; n ¼ 8 ðP ¼ 0:025Þ ð4Þ KeðB11þD34Þ ¼ 0:01610  0:9  105  ele; r ¼ 0:73; n ¼ 12 ðP ¼ 0:007Þ ð5Þ KeðD13Þ ¼ 0:03680  2:5:9  105  ele; r ¼ 0:68; n ¼ 12 ðP ¼ 0:015Þ ð6Þ where Ke in Mg h MJ-1 mm-1 is the estimated soil erodibility factor, ele is the elevation in m a.s.l., r is the correlation coefficient, n is the number of samples employed in the regression analysis, P is the significance level. Ke subindex meaning (land class) is specified in Table 2. The above reported LTFs have been applied to 109 new locations. These are locations whose land system class memberships could be identified based on expert knowledge. Sampled (114 samples) and estimated (109) soil erodibility values resulted in a composed data set of 223 samples, as represented in Fig. 2b. Some internal areas and the eastern side of the basin resulted still unvisited, due both to lack of information on land system class Table 2 Land system classes explored with the available data set 123 Land system class No. of samples Description B11 6 Internal calcareous relieves with ash sediments D13 12 Clayey hills of Cilento D34 6 Marl-calcareous and marl-arenaceous high hills of Irpinia and of upper Sele basin E22 5 Marl-arenaceous coastal hills of Cilento E23 8 Clayey coastal hills of Cilento Nat Hazards Fig. 4 Voronoi maps highlighting the local outliers distribution of the sampled data (a), and for the composed data set (b) membership and to scarce soil sampling data, which prevented to explore the correlation between soil erodibility and other auxiliary variables. It is interesting to compare the pattern of the sampled erodibility values to the pattern of the composed data set by evaluating the corresponding Voronoi maps constructed with Thiessen polygons around each location (Fig. 4a, b). In Fig. 4a, the sampled erodibility data appear spatially random, and they can hardly be interpolated to draw a map across the study catchment. In Fig. 4b, a declustered pattern of the composed data set of 223 samples is represented. The second pattern appears more suitable for interpolation. 2.4 Kriging Kriging is a generic name, adopted by the geostatisticians for a family of generalized leastsquares regression algorithms. The basic idea is to estimate the unknown attribute value at the unsampled location so as a linear combination of the neighbouring observations. A covariogram model, representative of the spatial structure of the targeted variable, is fitted to the sampled set and is used to determine weights for the neighbouring samples considered in the inference process. Remembering that the covariance function measures the average degree of similarity between an unsampled value z(so) and the nearby data values, a kriged estimate of the soil variable at the location so is given by: Z  ðso Þ ¼ n X ki  zðsi Þ ð7Þ i¼1 where z is the n-vector of si observed primary data (in our case they correspond to K– factor), selected in the so neighbourhood; ki is the weight vector associated with the distance ho(i) (separating so from si) and calculated by solving of the kriging simultaneous equation system (Johnston et al. 2001). In this study, ordinary kriging is applied to both the composed data set and the sampled data set only. A model of regionalization was fitted by using an iterative procedure developed by Johnston et al. (2001) and consisted of two stages. 123 Nat Hazards Stage 1 begins by assuming an isotropic and spherical model and computing the k k 1 empirical covariance functions on the scaled data Z^k j ðsi Þ ¼ Zj ðsi Þ  sk , where Zj ðsi Þ is used to denote the jth measurement of variable type k at the ith spatial locations si, and sk is the sampled standard deviation. At this first stage, it was also accepted that soil erodibility collected close to one another is more similar than sampling further apart. Hence, property values lie on a continuum between two extremes and will exhibit a relationship between spatial dependence and distance (Webster and Oliver 1992). In stage 2, the parameters of the model are interactively calibrated, such as the number of lags, the lag size h, nugget and range. Once the isotropic assumption was verified, the covariance functions were modelled as a combination of two distinct spatial structures: nugget variance and a K–Bessel model for the composed data set (Fig. 5a), and nugget plus exponential model for the sampled data set (Fig. 5b). In particular, K-Bessel model was set following the expression (Johnston et al. 2001): 2 3  Xhk khk hk   h X h k k r 6 7 hk ð8Þ Kh Cðh; hÞ ¼ hs 4 hk 1 5 for all h hr 2 Cðhk Þ k The partial sill hs C 0 was estimated equal to 0.00003765, while range hr C 0 is estimated equal to 24,000 m. Xhk is a scaling parameter assumed equal to 10, Khk(•) is the modified Bessel function of the second kind of order hk (Abramowitz and Stegun 1965), while C(hk) is the gamma function equal to: CðyÞ ¼ Z1 xðx1Þ expðxÞdx ð9Þ 0 In Fig. 5, it is possible to detect that the experimental covariogram fits well to a nonstandard function, such as the K-Bessel model, similar to a Gaussian one. Both modelcurves are characterized by a finite covariance. As referred by Krasilnikov (2008), it is important to notice that such a form of covariogram model might also result from an abrupt change of values (e.g. a change in soil classes), particularly when we deal with a covariogram computed in directions perpendicular to the boundary between two distinct regions. It should be also noted that the covariogram of the composed data set shows a smaller nugget, with a reduction of about 30% when compared with the covariogram of the sampled data set (Fig. 5b). The nugget represents the unexplained or random variance that may be caused by measurement errors and/or property variability which cannot be detected at the sampling scale. Fig. 5 Experimental covariogram (dots) and fitted model (curves) for erodibility composed data set (a), and for only erodibility sampled data (b) 123 Nat Hazards 3 Results Figure 6a shows the erodibility maps defined by using ordinary kriging based on composed data set approach and the corresponding kriging error (Fig. 6b) at a 0.5 km 9 0.5 km grid resolution. The estimates of the areal average over the whole region and the standard deviation are 0.026 ± 0.0049 Mg h MJ-1mm-1, respectively, with a median of 0.028 Mg h MJ-1 mm-1. The computed K values span from 0.010 to 0.045 Mg h MJ-1 mm-1, thus within the range estimated by the European Soil Bureau at a larger scale (Van der Knijff et al. 2000). Kriging mapping process is not only used to describe spatial structures, but can be also used to understand or begin to explore the underlying processes that are responsible for soil property variation (Trangmar et al. 1985). For instance, the 24,000 m–range of the Bessel– model suggests a moderate–low spatial variability in soil erodibility, thus producing a mosaic pattern of source and sink areas. From the kriging uncertainty map (Fig. 6b), it is possible to observe negligible errors (around to 0.004 Mg h MJ-1 mm-1). Just slightly higher errors are expected around the Eastern watershed of the basin (around to 0.006 Mg h MJ-1 mm-1). Significant errors (around to 0.010 Mg h MJ-1 mm-1) have been found only in the Eastern lands out of Sele basin. The performances of the two approaches, OK applied to the composed data set (CO_OK) and OK applied to the sampled data only, were assessed and compared by using cross-validation. Cross-validation procedures use all the data to estimate the autocorrelation model. Each K datum location is removed, one at a time, and the associated value is predicted. The interpolated and actual values are then compared. The following crossvalidation statistics have been calculated: Mean Errors (ME), Root Mean Square Error (RMSE), Average Standard Error (ASE), Mean Standard Error (MSE) and Root-MeanSquare Standardized Error (RMSSE). ME represents the overall bias. MAE, RMSE and Fig. 6 Spatial variability by CO_OK of topsoil erodibility (a), and kriging standard error map (b), across and around the Sele River Basin 123 Nat Hazards Table 3 Cross-validation statistics of soil erodibility maps (Mg h MJ-1 mm-1) derived by composed data set ordinary kriging (CO_OK) and sampled data set ordinary kriging (OK) Error statistics CO_OK OK N. samples 223 114 ME -0.0000193 -0.0000277 RMSE 0.0046610 0.0049190 ASE 0.0045560 0.0047420 MSE -0.0007913 0.0062010 RMSSE 1.046 1.044 ME mean errors, RMSE root mean square error, ASE average standard error, MSE mean standard error, RMSSE root-mean-square standardized error ASE are indices that represent the variability prediction, while RMSSE compares the error variance with the theoretical variance, i.e. the kriging variance. These cross-validation statistics are reported in Table 3. The cross-validation statistics show that the usage of supplementary information, derived from local topotransfer functions, can help to compensate for the lack of local erodibility data from which any interpolation method, such as kriging with only sampled data, might suffer. This is confirmed also by the scatterplots drawn during the crossvalidation stage (Fig. 7a, b). Afterwards, in order to simulate possible undersampling, a subset of erodibility data was randomly subtracted from the complete pattern and it was successively used at validation stage. The results are depicted in Fig. 7a1, which also shows that the proposed approach can provide satisfactory prediction performance. 4 Conclusions Predicted erodibility (x 100) Although human judgment involved in information construction is an additional source of uncertainty, combining more types of data under a geostatistical approach can be a successful strategy for improving soil erodibility mapping in undersampled regions. This work evaluated the possibility to estimate patterns of soil erodibility at catchment scale by interpolating both data obtained from soil samples and data estimated with auxiliary terrain data and land system class memberships. This approach has been applied to the Sele Basin, 4.40 4.40 a b 4.40 3.73 3.73 3.73 3.07 3.07 3.07 2.40 2.40 2.40 1.73 1.73 1.07 y = 0.587x + 1.015 2 R = 0.61 (N = 223) 0.40 0.40 1.07 1.73 2.40 3.07 3.73 4.40 1.07 a1 1.73 y = 0.174x + 2.287 2 R = 0.14 (N = 114) 0.40 0.40 1.07 1.73 2.40 3.07 3.73 4.40 1.07 y = 0.645x + 0.847 2 R = 0.69 (N = 114) 0.40 0.40 1.07 1.73 2.40 3.07 3.73 4.40 Actual erodibility (x 100) Fig. 7 Scatter diagrams among sampled and estimated soil erodibility values according to the two kriging approaches: CO_OK (a) and OK (b) at cross-validation stages, and CO_OK at validation stage (a1) 123 Nat Hazards Southern Italy, where measurements of soil properties were too scarce to allow an effective description of the spatial variability of soil erodibility. Estimates of soil erodibility at unsampled locations have been obtained with Local Topotransfer Functions. In this basin, elevation was the best explanatory variable for soil erodibility likely because it is characterized by a large elevation range (0–1,899 m a.s.l.). Of course, in other study cases with less contrasting altitudes, some other terrain attributes could better explain soil erodibility. Estimated and sampled data have been merged to compose a larger data set and interpolated with ordinary kriging (OK). The covariogram of the composed data set allowed a 30% reduction in the nugget effect (random variability), with a significant increase in the accuracy of the predicted soil erodibility map when compared with the map predicted with OK applied to the sampled data only. The proposed approach offers effective spatial predictions, and it is exportable to regions where financial costs for soil sampling are not feasible. Acknowledgments Annamaria Castrignanò, Fabio Terribile, Amedeo D’Antonio and Antonio Di Gennaro are gratefully acknowledged for the long and very useful discussions and for their valuable suggestions. The KUERY software package is available for free available for free download at www.irpi.fi.cnr.it/software.html. We are grateful to two anonymous reviewers for useful comments. This study was financially supported by VECTOR Project (line 2 VULCOST- chief, Bruno D’Argenio). References Abbaspour KC, Schulin R, van Genuchten MT, Schläppi E (1998) An alternative to cokriging for situations with small sample sizes. Math Geol 30:259–274 Abramowitz M, Stegun IA (1965) Handbook of mathematical functions. Dover, NY, USA Alison BT, Moore KJ, Bullock DG, Dixon PM, Burras CL (2005) Improving map accuracy of soil variables using soil electrical conductivity as a covariate. Precis Agricolture 3:255–270 Baskan O, Dengiz O (2008) Comparison of traditional and geostatistical methods to estimate soil erodibility factor. Arid Land Res Manage 22:29–45 Bayramin İ, Basaran M, Erpul G, Cang MR (2008) Assessing the effects of land use changes on soil sensitivity to erosion in a highland ecosystem of semi-arid Turkey. Environ Monit Assess 140:249–265 Borselli L, Cassi P, Salvador Sanchis MP (2009) Soil erodibility assessment for application at watershed scale. In: Costantini EAC (ed) Manual of methods for soil and land evaluation. Science Publisher Inc, Enfield, NH, p 600 Bourennane H, King D, Coututirer A, Nicoullaud B, Mary B, Richard G (2007) Uncertainty assessment of soil water content spatial patterns using geostatistical simulation: an empirical comparison of a simulation accounting for single attribute and a simulation accounting for secondary information. Ecol Modell 205:323–335 Bouyoucos GJ (1962) Hydrometer method improved for making particle-size analyses of soils. Agron J 54:464–465 Burgess TM, Webster R (1980) Optimal interpolation and isarithmic mapping of soil properties: the semivariogram and punctual kriging. J Soil Sci 31:315–331 Casa R, Castrignanò A (2008) Analysis of spatial relationships between soil and crop variables in a durum wheat field using a multivariate geostatistical approach. Eur J Agron 28:331–342 Castrignanò A, Lopez N (2000) La geostatistica nella scienza del suolo: stato dell’arte e prospettive future. In: Stefano Bocchi (ed), Proceeding geostatistica per lo studio e la gestione della variabilità. Applicazioni nelle scienze fisiche, ambientali e agronomiche, University of Milan, pp 61–86 Castrignanò A, Buttafuoco G, Canu A, Zucca C, Madrau S (2008) Modelling spatial uncertainty of soil erodibility factor using joint stochastic simulation. Land degrad dev 19:198–213 Dalal-Clayton B, Dent D (2001) Knowledge of the land. Lad resources information and its use in rural development. Oxford University Press, New York, p 448 Di Gennaro A (2002) I sistemi di terre della Campania. Carta 1:250.000 e Legenda. Selca, Firenze Diodato N (2005) The influence of topographic co-variables on the spatial variability of precipitation over little regions of complex terrain. Int J Climatol 25:351–363 123 Nat Hazards Diodato N (2006) Idroinformatics system for pollutant potential leaching spatial uncertainty assessment. Environ Geosci 13:227–238 Diodato N, Ceccarelli M (2004) Multivariate indicator kriging approach using a GIS to classify soil degradation over Mediterranean agricultural lands. Ecol Indic 4:177–187 Diodato N, Fagnano M (2010) A simple geospatial model climate-based for designing erosive rainfall pattern. In: El Nemr A (ed) The environmental pollution and its relation to climate change. Nova Science Publishers, Hauppauge, NY. ISBN 978-1-61122-144-2 FAO (1995) Planning for sustainable use of land resources. Toward a new approach. Land and Water Buletin, 2, Rome Goovaerts P (2000) Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. J Hydrol 228:113–129 Hengl T, Heuvelink GBM, Stein A (2004) A generic framework for spatial prediction of soil variables based on regression-kriging. Geoderma 120:75–93 Herbst M, Diekkrüger B, Vereecken H (2006) Geostatistical co-regionalization of soil hydraulic properties in a micro-scale catchment using terrain attributes. Geoderma 132:206–221 Illian J, Penttinen A, Stoyan H, Stoyan D (2008) Statistical analysis and modelling of spatial poin patterns. John Wiley and Sons Ltd, Chichester, pp 264–265 Jenny H (1980) The soil resource, origin and behaviour, Springer-890 Verlag, New York Johnston K, Ver Hoef JM, Krivoruchko K, Lucas N (2001) Using ArcGISTM geostatistical analyst. ESRI, Redlands, p 300 Krasilnikov P (2008) Variography of discrete soil properties. In: Krasilnikov P, Carré F, Montanarella L (eds), Soil geography and geostatistics, European commission—JRC, Report EUR 23290 EN, pp 12–25 Matula S, Špongrová K (2007) Pedotransfer function application for estimation of soil hydrophysical properties using parametric methods. Plant Soil Environ 53:149–157 Monestiez P, Courault D, Allard D, Ruget F (2001) Spatial interpolation of air temperature using environmental context: application to a crop model. Environ Ecol Stat 8:297–309 Ozcan AU, Erpul G, Basaran M, Erdogan HE (2008) Use of USLE/GIS technology integrated with geostatistics to assess soil erosion risk in different land uses of Indagi Mountain Pass—Çankırı, Turkey. Env Geol 53:1731–1741 Park S, Egbert SL (2005) Assessment of soil erodibility indices for conservation reserve program lands in Southwestern Kansas using satellite imagery and GIS techniques. Environ Manage 36:886–898 Renard KG, Foster GR, Weesies GA, McCool DK, Yoder DC (1997) Predicting soil erosion by water: a guide to conservation planning with the revised universal soil loss equation (RUSLE). Agricultural Handbook, United States Department of Agriculture, Washington, DC Rodriguez-Iturbe I, Marani M, D’Odorico P, Rinaldo A (1998) On space-time scaling of cumulated rainfall fields. Water Resour Res 34:3461–3469 Salski A (2006) Ecological applications of fuzzy logic. In: Recknagel F (ed) Ecological informatics. Springer, Berlin, pp 3–14 Salvador Sanchis MP, Torri D, Borselli L, Poesen J (2008) Climate effects on soil erodibility. Earth Surf Proc Land 33:1082–1097 Scull P, Franklin J, Chadwick OA, McArthur D (2003) Predictive soil mapping: a review. Prog Phys Geogr 27(2):171–197 Shirazi MA, Boersma L, Hart W (1988) A unifying analysis of soil texture: improvement of precision and extension of scale. Soil Sci Soc Am 52:181–190 Sigua GC, Hudnall WH (2008) Kriging analysis of soil properties: implication to landscape management and productivity improvement. J Soils Sediments 8:193–202 Simbahan GC, Dobermann A, Goovaerts P, Ping J, Haddix ML (2006) Fine-resolution mapping of soil organic carbon based on multivariate secondary data. Geoderma 132:471–489 Sinowski W, Auerswald K (1999) Using relief parameters in a discriminant analysis to stratify geological areas with different spatial variability of soil properties. Geoderma 89:113–128 Torri D, Poesen J, Borselli L (1997) Predictability and uncertainty of the soil erodibility factor using a global dataset. Catena 31:1–22 Trangmar BB, Yost RS, Uehara G (1985) Application of geostatistics to spatial studies of soil properties. Adv Agron 38:45–94 Van der Knijff JM, Jones RJA, Montanarella L (2000) Soil erosion risk assessment in Europe. EUR 19044 EN, European soil Bureau, space application institute. JRC, Ispra, p 34 Veihe A (2002) The spatial variability of erodibility and its relation to soil types: a study from northern Ghana. Geoderma 106:101–120 123 Nat Hazards Walkley A, Black IA (1934) An examination of the Degtjareff method for determining soil organic matter and a proposed modification of the chromic acid titration method. Soil Sci 37:29–38 Wang G, Gertner G, Fang S, Anderson AB (2003) Mapping multiple variables for predicting soil loss by geostatistical methods with TM images and a slope map. Photogramm Eng Remote Sensing 69:889–898 Webster R, Oliver MA (1992) Sample adequately to estimate variograms of soil properties. J Soil Sci 43:177–192 Wessolek G, Duijnisveld WHM, Trinks S (2008) Hydro-pedotransfer functions (HPTFs) for predicting annual percolation rate on a regional scale. J Hydrol 356:17–27 Yugang W, Duning X, Yan L (2007) Temporal–spatial change in soil degradation and its relationship with landscape types in a desert–oasis ecotone: a case study in the Fubei region of Xinjiang Province, China. Env Geol 51:1019–1028 Zucca C, Canu A, Madrau S, Castrignanò A (1999) Utilizzo di tecniche GIS e geostatistiche per la produzione di mappe del rischio erosivo. Proceedings of the 3rd national conference ASITA, vol II, Napoli, pp 1269–1274 123