Spatial variability of soil properties and soil erodibility in the Alqueva reservoir watershed

The aim of this work is to investigate how the spatial variability of soil properties and soil erodibility (K factor) were affected by the changes in land use allowed by irrigation with water from a reservoir in a semiarid area. To this end, three areas representative of different land uses (agroforestry grassland, lucerne crop and olive orchard) were studied within a 900 ha farm. The interrelationships between variables were analyzed by multivariate techniques and extrapolated using geostatistics. The results confirmed differences between land uses for all properties analyzed, which was explained mainly by the existence of diverse management practices (tillage, fertilization and irrigation), vegetation cover and local soil characteristics. Soil organic matter, clay and nitrogen content decreased significantly, while the K factor increased with intensive cultivation. The HJ-Biplot methodology was used to represent the variation of soil erodibility properties grouped in land uses. Native grassland was the least correlated with the other land uses. The K factor demonstrated high correlation mainly with very fine sand and silt. The maps produced with geostatistics were crucial to understand the current spatial variability in the Alqueva region. Facing the intensification of land-use conversion, a sustainable management is needed to introduce protective measures to control soil erosion.


Introduction
Soil erosion is a significant economic and environmental problem worldwide as a driving force affecting landscapes (Zhao et al., 2013).It is a very dynamic and complex pro-cess, characterized by the decline of soil quality and productivity, as it causes the loss of topsoil and increases runoff (Lal, 2001;Yang et al., 2003).Furthermore, soil erosion often causes negative downstream impacts, such as sedimentation in rivers and reservoirs, decreasing their storage volume as well as lifespan (Pandey et al., 2007;Haregeweyn et al., 2013).
One of the main causes of soil loss intensification around the world is associated with land-use change (Leh et al., 2013).The relationship between different land use and soil susceptibility to erosion has attracted the interest of a variety of researchers (Yang et al., 2003;Cerdà and Doerr, 2007;Blavet et al., 2009;Biro et al., 2013;Wang and Shao, 2013), who have shown the impact of changes in vegetation cover and agricultural practices on soil properties and therefore on overland flow.Generally, cultivated lands experience the highest erosion yield (Cerdà et al., 2009;Mandal and Sharda, 2013).In the Mediterranean regions, in combination with these anthropogenic factors, climate change has amplified the concerns about soil erosion since it is expected that there will be an increase of dry periods followed by heavy storms with concentrated rainfall (Nunes et al., 2009).
Some models have been developed to predict soil loss and sediment delivery.The Revised Universal Soil Loss Equation (RUSLE) is the most used empirical equation for modeling annual soil loss from agricultural watersheds (Renard et al., 1997).The susceptibility of soil erosion and land degradation depends largely on various inherent soil properties, namely chemical, physical, biological and mineralogical properties (Cambardella et al., 1994;Pérez-Rodríguez et al., 2007).However, according to the RUSLE model only some Published by Copernicus Publications on behalf of the European Geosciences Union.
V. Ferreira et al.: Spatial variability of soil properties and soil erodibility of the soil's properties define soil erodibility (K factor), such as particle-size composition, the content of organic matter, soil structure and permeability.Therefore, the K factor is the most used and is an important index to measure soil susceptibility to erosion (Panagopoulos and Antunes, 2008).
Spatial variability in soils occurs naturally as a result of complex interactions between geology, topography and climate.Moreover the spatial variability of soil properties, which influence soil susceptibility to erosion, is highly related to anthropogenic factors particularly in cultivated lands (Paz-González et al., 2000;Wang and Shao, 2013).Thus, information on the spatial variability and the interactions between soil properties is essential for understanding the ecosystem processes and planning sustainable soil management alternatives for specific land uses (Pérez-Rodríguez et al. 2007;Ziadat and Tamimeh, 2013).
Classical statistics and geostatistics methods have been widely applied in studies about spatial distribution of soil properties (Pérez-Rodríguez et al., 2007, Tesfahunegn et al., 2011).Geostatistical techniques based on predictions and simulations have been used to describe areas where predicted information is established by a limited number of samples (Goovaerts, 1997).Geostatistics provides tools for analyzing spatial variability structure and distribution of soil properties and evaluating their dependence (Panagopoulos et al., 2014).
The biplot methodology provides an added value for analyzing spatial variability of soil properties.This multivariate statistical technique allows the graphical representation of a large data matrix (Gabriel, 1971), whereby it is possible to interpret the relations between individuals (samples) and between variables, as well as between both.Biplot can also indicate clustering of units with close characteristics, showing inter-unit distances as well as displaying variances and correlations of the variables (Gallego-Álvarez et al., 2013).The HJ-Biplot permits not only the analysis of the behavior by sample but also the determination of which variable is responsible for such behavior (Garcia-Talegon et al., 1999), allowing a visual appraisal to establish relations between soil properties and land uses.
The construction of the Alqueva dam in a semiarid area of southern Portugal created one of the largest artificial lakes in Europe.Taking advantage of water availability from the reservoir, this Mediterranean region has been subjected to land-use conversion from the native montado grassland to intensive agricultural uses.Land-use conversion from the native ecosystem to agriculture may alter physical, chemical and biological soil properties, which consequently may increase soil erosion and siltation in the reservoir.Soil erosion in the area has to be carefully evaluated in order to undertake sustainable soil management measures.Therefore, the aim of this study was to evaluate the effects of cultivation practices on some chemical and physical soil properties and on soil erodibility (K factor on RUSLE), and to characterize their spatial variability using geostatistics and HJ-Biplot methodology.

Study area
Located in the semiarid Alentejo region of Portugal, at the Guadiana River, the Alqueva reservoir (8 • 30 W, 38 • 30 N) covers an area of 250 km 2 , and the capacity of the reservoir is 4.15 km 3 .The main arguments for the implementation of what is considered the largest artificial lake in Europe were based on the need to combat the growing effects of desertification and to prevent the annual and monthly fluctuations in precipitation.One of the main goals of the Alqueva Multipurpose Project was the implementation of 120 000 ha of newly irrigated land in the Alentejo.The Alentejo region, covering an area of 27 000 km 2 , is considered one of the most depressed regions of the European Union and characterized by a Mediterranean climate with very hot and dry summers and mild winters.The average temperature ranges from 24 to 28 • C in hot months (July/August) and from 8 to 11 • C in cold months (December/January).The average annual precipitation at the nearest meteorological station, for the last 30 years, is 517.2 mm.The region is affected by intense dry periods followed by heavy, erosive rains concentrated in the autumn season.
The study experimental site (farm "Herdade dos Gregos"), located in the surrounding area of the reservoir (Fig. 1), is a private property with 900 ha.The landscape is characterized by its hilly topography with significant altitude variations (mainly between 100 and 250 m).The bedrock of the study area is rocky, and, according to the World Reference Base for Soil Resources (FAO, 2006), the two types of soil in this area are Haplic Luvisols (LVha) and Lithic Leptosols (LPli).This farm was selected to include a diversity of land uses, including native montado grassland and more intensive land uses, with irrigation, namely olive tree orchard and lucerne cultivation.Direct pumping from Alqueva reservoir is done on this private property since it is near the reservoir.
The typical landscape in the Alentejo region is the montado native grassland, an agrosilvopastoral system characterized by savannah-like, low-density woodlands with evergreen holm oaks (Quercus ilex).For that reason, an area of the montado grassland (20.7 ha), used as a permanent pasture for the cattle, was selected for this study.This small area is located in the high altitudes of the "Herdade dos Gregos" (from 200 to 240 m) with a slope that varies from 1.4 to 20.9 %.Tillage (at about 15 cm depths) was done only once every 10 years to decrease shrub competition (the most recent one was 4 years before the study implementation), and the soil is not subjected to any fertilizer.Four years before the study implementation, there was a fire on this agrosilvopastoral area of the farm.
Taking advantage of the water availability, another land use (with 33.5 ha) is an irrigation area (pivot sprinkler irrigation system) on which lucerne (Medicago sativa) is sown four times a year.Lucerne, once dried, is nutritional for cat- tle, and it incorporates nitrogen in the soil.In this area, conventional tillage is used, involving multiple aspects: plough (about 20 cm depth) in fall, fallowing cultivator (about 15 cm depths) and disc harrow (about 10 cm depths) subsequent to soil tillage.Inorganic fertilizers were applied to the cultivated field at a rate of 100 kg NPK ha −1 .This land use is placed in the midland (194-220 m), and the slope varies from 0 to 9 %.
Another irrigated land use consists of an olive tree plantation (57.5 ha), which is done in strips.This cultivation has a drip irrigation system, is fertilized once every 2 years and is ploughed once a year to decrease weed competition.The olive orchard is located in the low elevations of the farm (150-186 m), and it is on the side of the reservoir (Fig. 1).The slope varies from 0 to 14.2 %.

Soil sampling and laboratory analysis
Since the objective was to study the relation between soil properties and the K factor from RUSLE, the soil samples were collected from 0 to 20 cm depth, according to Renard et al. (1997).In order to predict variations in short distances, 25, 27 and 52 soil samples were randomly collected respectively in montado, lucerne and the olive orchard (see Fig. 1).
Samples were air-dried and then dried for about 6 h at 40 • C on a ventilated oven, and they were passed through a 2 mm sieve to remove rocks and gravel.The particle-size distribution was determined by the Bouyoucos hydrometer method (Bouyoucos, 1936).Soil organic matter content was determined using the Walkley and Black (1934) method, a wet oxidation procedure.The soil's total nitrogen content was determined according to Kjeldhal digestion, distillation and the titration method (Bremmer and Mulvaney, 1982).Soil pH and electrical conductivity were measured with a glass electrode in a 1 : 2.5 soil-water suspension (Watson and Brown, 2011).

Soil erodibility factor
Soil erodibility factor (K) (Mg ha h ha −1 MJ −1 mm −1 ) was estimated using soil property values -such as particle-size composition, content of organic matter, soil structure and permeability -in the 104 sample points described above.This factor represents the soil-loss rate per erosion index unit for a specified soil as measured on a standard plot (Renard et al., 1997).An algebraic approximation of the nomograph (Eq. 1) was used to estimate the K factor (Renard et al., 1997): .14(1) where OM is the percentage of organic matter, s is soil structure class, p is permeability class and M is the product of the percentage of modified silt (silt particles and very fine sand) or the 0.002-0.1 mm size fraction and the sum of the percentage of silt and percentage of sand.K is expressed in SI units of Mg ha h ha −1 MJ −1 mm −1 .To estimate the permeability, the field-saturated hydraulic conductivity was measured in the field using a double-ring infiltrometer (six site measurements per land use, each one with five repetitions).Permeability class and soil structure class were defined in accordance with Renard et al. (1997).

Statistical and geostatistical analysis
Data were subjected to classical analysis using SPSS 17.0 software to obtain descriptive statistics, namely the mean, minimum and maximum; standard deviation (SD); coefficient of variation (CV); and skewness of each parameter.Soil data were introduced in the ArcGIS environment, and geostatistical analyses were performed using Geostatistical Analyst tool, in other to examine spatial distribution of soil properties.Prior to using geostatistics to obtain prediction maps, a preliminary analysis of data was done to check data normality and global directional trends.Skewness is the most common statistical parameter to identify a normal distribution that is confirmed with skewness values varying form −1 to +1.Data transformation to normal distribution was necessary for some soil properties, and geostatistical analysis tools were used (log or Box-Cox method).Trend analysis was performed to examine the presence of any global directional trend in our data, an overriding process that affects all measurements in a deterministic way (nonrandom).So, when necessary, the trend removal was done using geostatistical analysis tools to more accurately model the variation (Panagopoulos et al., 2006).
The geostatistical methodology is based on the creation of a semivariogram (SV), a graphical representation (Eq.2) that describes how samples are related to each other in space, and it is based on where γ (h) is the variance (the most related samples have lower values of variance), N (h) is the number of samples that can be grouped using vector h, Z i represents the value of the sample and Z i + h is the value of another sample located at a distance h from the initial sample Z i (Chiles and Delfiner, 1999).Ordinary Kriging (OK) was selected as a geostatistical method.OK is considered one of the most accurate interpolation techniques which assumes that variables close in space tend to be more similar than those further away (Goovaerts, 1999).
Using the Geostatistical Analyst tool (ArcGIS) and selecting the OK methods, a semivariogram was created for each measured property.In the Kriging method different semivariogram models can be used (e.g., spherical, exponential) and the selection is usually performed by employing the crossvalidation technique, which permits the evaluation of the prediction accuracy.Cross validation was executed to investigate the prediction performances through the statistical values, as the mean error (ME) or root-mean-square standardized error (RMSSE), which results from comparing the estimated semivariogram values and real observed values.Additional semivariogram parameters were analyzed to better understand the spatial structure and dependence of each variable.Nugget is the variance at distance zero and reflects the sampling error.Sill is the semivariance value at which the semivariogram reaches the upper bound and flattens out after its initial increase; it is the variance in which the samples are no longer spatially related to the study area.
Once the cross-validation process was completed, interpolation maps of spatial distribution, for each soil variable, were produced according to the semivariogram model selected, in the ArcGIS software.

HJ-Biplot
HJ-Biplot represents a matrix, without assumptions related to its probabilistic distribution, permitting a graphic representation of the geometric data structure, representing the data set (samples and variables) variability.The prefix "bi" is due to a simultaneous representation of the matrix rows and columns, searching for the maximum representation quality possible, at the same scale (Martín-Rodríguez et al., 2002;González-Cabrera et al., 2006;Gallego-Álvarez et al., 2013).
A data matrix X suffers a factorization to reduce its dimensionality through single-value decomposition, the algebraic base of biplot representation (Gabriel, 1971): where (r×r) is a diagonal (λ 1 , λ 2 , . .., λ r ) corresponding to the r eigenvalues of XX' or XX , U (n×r) is an orthogonal matrix whose columns are the eigenvectors of XX or X X and U (n×r) is an orthogonal matrix whose columns are the eigenvectors of XX .
With the MultBiplot software, developed by the University of Salamanca (Vicente Villardón, 2014), an HJ-Biplot was used to determine the relation between soil properties and land uses, and the correlations between both (soil properties and land uses), thereby defining patterns and clustering the samples in groups.
On the HJ-Biplot graphic representation, the points represent individuals (samples) and the vectors represent variables (in this case, chemical and physical soil properties).To interpret and discuss the graphs obtained with this methodology, it is essential to be aware of the following (Gallego-Álvarez et al., 2013): -The distance between points represents the variability and can be interpreted as similarity or dissimilarity; i.e., the close samples have similar behaviors.
-The angle formed by variable vectors is interpreted as correlation; i.e., small angles between variables represent similar behaviors with high positive correlations, and the obtuse angles that are almost a straight angle are associated with variables with high negative correlations; i.e., the cosine value of the angles represents the correlation between variables.
-The proximity of individual points and variable vectors means high preponderance; in other words the closer a point is to a variable vector, the more important this sample is to explain this variable.
-The length of the vector represents the variable's variability; the longer the vector, the higher this variability.

Descriptive statistics
The descriptive statistics of soil properties are given in the first part of Table 1.All measured parameters varied considerably within the areas (different land uses) as indicated by the coefficient of variation (varies from 4.2 to 70.2 %).Nitrogen (N) and organic matter (OM) show the highest variation values, especially for cultivated fields (lucerne cultivation and olive orchard), which can be explained with the lack of homogeneous fertilization or tillage practices applied to soil in these areas.
The skewness results, which vary from −1.48 to 3.54 in this study, indicated that some soil properties of the different uses were not normally distributed, especially OM and N. The principal reason for some soil properties having nonnormally distributions may be related to soil management practices (Tesfahunegn et al., 2011).As already mentioned, data were transformed to normal distribution when necessary (see Table 1).
These mean results show significant differences between land uses for all the properties analyzed.From the particle size distribution reported in Table 1, the soils are mostly sandy loam, formed mainly of sand, followed by silt and low quantities of clay.However, there are some differences between land use areas that can be explained by soil type.The LPli soils are characterized by a thin layer (about 10 cm), in that case upon a schist rock, justifying the higher clay content at the montado grassland.The LVha soils in the lucerne cultivation and the olive orchard are characterized by a loam or sandy loam layer (first 20 cm) with good drainage over clayenriched subsoil (upon a basic crystalline rock), explaining the lower values of clay and fine sand, especially in the olive orchard.Despite the same soil type, soil texture is different between lucerne and the olive orchard, which can be justified by land use.The lucerne is a more intensive cultivation (intensive irrigation, tillage and continuous cultivation; fertilizers and lime application), involving conditions that promote changes in the soil weathering and moisture and, consequently, in soil texture (Yimer et al., 2008).On the other hand the soil between olive trees is kept without vegetation for most of the year and can explain the clay drainage to a sub-layer.

V. Ferreira et al.: Spatial variability of soil properties and soil erodibility
Montado shows the highest content of OM (5.22 %), whereas lucerne and olive fields show the lowest values (2.08 and 2.10 %, respectively).Other studies suggest that OM is higher in no-tillage soils compared to minimum tillage that increases aeration (Celik, 2005).Tillage mixes the subsoil with topsoil; after soil erosion, the nutrients are easily leached and the surface becomes poor in nutrients (Al-Kaisi and Licht, 2005).As for OM, the highest value of N nutrient occurs in the montado (0.19 %) and the lowest values in lucerne (0.11 %) and the olive orchard (0.10 %), which is related to the tillage practice that is frequently employed in these last two land uses, while in the montado grassland the cattle enriches the soil.
Soil EC values (Table 1) were similar when comparing the montado grassland (0.100 dS cm −1 ) and the lucerne field (0.107 dS cm −1 ); they were slightly higher in the olive orchard (0.182 dS cm −1 ) but not enough to raise salinity problems.Usually, the addition of fertilizers (that happens on lucerne and the olive orchard) can cause high EC due to the percentage of the salts which are leached by water irrigation (higher in the lucerne field).
The soil pH was significantly higher in the lucerne cultivation land (7.1) compared to the montado grassland (5.9) or in the olive tree orchard (5.5) (Table 1).The soil pH in the lucerne was greater due to lime application to increment the soil pH in that area.Lucerne's optimum pH for production is between 6.5 and 7.2, and lime application has been found to produce a significant improvement in nodulation of lucerne (both number and dry weight of nodules per plant) (Grewal and Williams, 2001).
Saturated hydraulic conductivity (HC) values were greater in the lucerne area (5.95 cm h −1 ), slightly lower in the montado grassland (4.56 cm h −1 ) and lowest in the olive orchard (2.60 cm h −1 ).The lower permeability in the olive orchard can be explained by the clay-enriched subsoil or soil crust problems, and it may explain the higher values of EC, i.e., the greater concentration of salts.Also it can be explained by the frequency of tillage in the different land uses because aggregate stability and water infiltration rate are higher in soils subjected to limited tillage systems (Alvarez and Steinbach, 2009).
As a result, the K factor was different for the typical land use, montado grassland, compared to the lucerne cultivation and the olive orchard.The values increased with the intensification of the cultivation field, with the lowest values for montado grassland (0.021 Mg ha h ha −1 MJ −1 mm −1 ) and the highest for the lucerne cultivation (0.039 Mg ha h ha −1 MJ −1 mm −1 ) and the olive orchard (0.038 Mg ha h ha −1 MJ −1 mm −1 ).Other studies had similar results, showing that the removal of permanent vegetation, the loss of OM and the reduction of aggregation, caused by intensive cultivation, contribute to decrease the K factor (Celik, 2005).

Spatial dependence of soil properties
Model selection for each soil property was based on the nugget, sill, ME and the RMSSE presented in the second part of Table 1

(Geostatistics).
Nugget is low in most soil properties studied, implying strong spatial dependence.The nugget-to-sill ratio is used to define spatial dependence of soil properties: if the ratio is < 0.25, there is strong spatial dependence; if it is 0.25 to 0.75, there is moderate spatial dependence; and if the ratio is > 0.75, spatial dependence is weak (Cambardella et al., 1994).As shown in Table 1 the ratio values indicate the presence of high to moderate spatial dependence for all soil parameters (values between 0 and 0.64).In general, there is stronger spatial dependence in montado (low nugget-to-sill ratio), which can be explained with the non-existence of extrinsic factors, such as management cultivation practices, that influence soil properties, and soil is left as it is for permanent pasture.
Cross validation facilitated the selection of the best-fit semivariogram for an interpolation map, which could provide the most accurate predictions.Closer values of the ME to 0, and closer values of the RMSS to 1, suggested that the prediction values were close to measured values (Wackkernagel, 1995).Most of the soil properties were best fitted with an exponential model, particularly in the montado area and olive orchard, whereas in lucerne the Gaussian, circular and stable semivariogram models were used.

Spatial distribution
The interpolation maps obtained with geostatistics are useful to better understand spatial variability and its influences.The variability of spatial soil properties can be influenced by natural factors (such as particle-size composition and topography) and anthropogenic factors (such as land cover or management practices) (Tesfahunegn et al., 2011).Sometimes, the effect of some factors is at least 1 order of magnitude greater (as topography or soil type) than the land use.So, as mentioned trend analysis was performed to study the existence of directional trends caused by these factors with large scale of variation, and it is shown in Fig. 2. Global trend exists if a curve that is not flat (i.e., a polynomial equation) can be fitted to the data (for example for total N in montado or very fine sand (VFS) in the olive orchard).These trends were identified for part of the soil properties and for different land uses (Fig. 2).The strongest influence of a directional trend was identified from southeast to the northwest, which could be associated with the topography (Fig. 1) since the altitudes increase in accordance with these direction.So, trend removal is crucial to create more accurate prediction maps in order to justify an assumption of normality.
The interpolation maps for some studied soil properties are shown in Fig. 3. Through looking at the VFS distribution, it was noticed that the higher fractions of these particles (Fig. 3) were measured at low altitudes or on flat slopes such as the valley (see elevation in Fig. 1).This can be explained by erosion-deposition processes because these particles are easily detached and transported by water.
The highest percentages of N and OM were found on montado, as discussed previously.These two properties present similar distributions for all land uses.The nitrogen existing in the soil is mostly organic, and the inorganic forms (ammonium and nitrate) are easily leached or assimilated by plants.So, when OM breaks down due to mineralization, the N fraction decreases (Varennes, 2003).There were higher values in montado because the soil is not frequently tilled as it is in the other land uses.In the lucerne cultivation and the olive orchard, the variation of OM and N can be explained by inadequate management practices (e.g., inadequate fertilization rates, tillage, irrigation rates, seed rates).
Figure 3 illustrates the interpolation map for the K factor which was estimated through the Wischmeier nomograph (Eq.1).The values vary from 0.006 to 0.061 Mg ha h ha −1 MJ −1 mm −1 , and the prediction map shows the highest values for lucerne and the olive orchard, especially where the soils have more silt and VFS, along with less OM and N (see HJ-Biplot).In the surrounding area of the reservoir, the types of soil differ with the topography and land use; therefore, knowledge of soil properties is fundamental when facing the intensification of cultivation that could increase the K factor.These intensive practices decrease OM in soils, making them poor and vulnerable to the soil erosion process.When looking for natural vs. anthropogenic impact on the K factor, for each land use, it is evident that in the montado the spatial variability is mainly associated with natural (intrinsic) factors (as texture), with soil properties and erodibility distribution being more homogenous.In the lucerne and olive orchard the spatial variability is more dependent on inhomogeneous anthropogenic causes such as fertilization and irrigation rates and tillage/plough processes.

HJ-Biplot
The HJ-Biplot representation matrix of soil properties is showed in Fig. 4. It was observed that the dominant axis (axis 1) takes 35.83 % of the total inertia (information) of the system.With both dimensions, an accumulative inertia of 61.04 % was achieved.Regarding this graphic representation, it was observed that samples were grouped according to the land use.The montado samples were close to OM, N and clay vectors, showing their preponderance to be a characterization of these variables.The lucerne samples were important to describe the pH and silt content.On the other hand the olive samples were more disperse but related to EC, permeability class, sand, VFS and K.
The variables demonstrating a more positive correlation were OM and N, as previously noticed.Clay and silt were also positively correlated, but they were negatively correlated with sand, as expected, because soils with more sand have less clay and/or silt.
Through the matrix representation it was detected that soils with more sand have higher EC (olive orchard), although EC normally increases with the percentage of clay.This may be explained by the addition of fertilizers, as previously discussed, that can contribute to an EC increase.These results for EC show low variability between land uses, revealing a low cation exchange capacity of these soils.This is fre-  et al., 2000).
Permeability class increases as the HC sat decreases, as defined by Renard et al. (1997).So, contrary to what was expected, for this study the soils with more sand (occurring in the olive orchard) have less hydraulic conductivity (highpermeability class).It can be explained by a clay-enriched sub-layer under the sandy loam layer or/and by the soil compaction/degradation processes.The soil compaction and degradation can be related to repeated plow operations to reduce shrubs between olive rows and irrigation (Pagliai et al., 2004).This permeability decrease in the olive orchard was correlated with the increase of the K factor.
Nevertheless, the properties more positively correlated with the K factor were the VFS and silt; this is due to the susceptibility of these particles to erosion since they can be easily detached and transported by water (Morgan, 2005).The OM and N content were negatively correlated with K and permeability.The higher OM reduces the susceptibility of the soil to detachment and increases infiltration (Bronick and Lal, 2005).The N content is not used to estimate K; however, especially for soils without fertilization, the existent N is mostly associated with OM.Nevertheless, nutrients decrease in soils that are more erodible, according to the literature (Tesfahunegn et al., 2011).The clay content also shows a negative correlation with K factor, as expected (Renard et al., 1997).
Figure 5 shows the hierarchical cluster representation.Using HJ-Biplot methodology and the aggregation tool ward, three clusters were obtained.The samples were grouped by land uses (that were already detected by the matrix representation; see Fig. 4).Cluster 1 is represented by a majority of samples from lucerne, Cluster 2 by samples from montado and Cluster 3 by samples from the olive orchard.This was explained by the effect of different management practices, vegetation cover and local soil characteristics, as discussed.Some samples in each land use had different values (higher or lower than the majority) and were grouped in a different cluster.Identifying the location of the sample, the cause of displacement can be studied and can help to improve land management practices.
Therefore, the cluster analysis is convenient to identify the effect of different land use and management on soil properties and consequently on soil erosion.On the other hand, the cluster analysis could support the delineation of zones according to soil properties, and subsequently according to erosion susceptibility, which could be used for site-specific soil management recommendations.

Conclusions
This study demonstrated that the variability of soil properties and the K factor is associated with land use, cultural practices (tillage type, fertilizer rates, conservation measures, etc.) and local conditions (complex topographic landscape, soil type, etc.).The K factor showed high correlation especially with organic matter, nitrogen, silt and very fine sand.Soils with intensively cultivated land use, and consequently with more tillage and irrigation, had lower organic matter and lower nitrogen content.This translates into a lower cation exchange capacity producing lower aggregate stability and, consequently, an increase of the K factor.
Therefore, in the surrounding area of the Alqueva reservoir, the ongoing change in land use and soil management practices can have a significant effect on chemical and physical soil properties.As a result, this affects the soil erodibility index, intensifying the risk of erosion.The increase of soil loss in the watershed might have a significant impact on a reservoir's ability to store water, reducing its lifespan.
Knowledge of soil spatial variability is fundamental for environment management and can help in the sustainable use of the resource soil.The prediction maps produced with geostatistics are an important monitoring tool, showing the exact position in the field of the specific soil properties.The HJ-Biplot methodology was demonstrated to be useful in gaining a better understanding of how soils properties were correlated and allowed not only a determination of the behavior by sample but also a conclusion as to which variable is responsible for such behavior.The simultaneous use of HJ-Biplot with geostatistics allows this information to be found on the map, which has important theoretical and practical significance for precision agriculture.Facing the intensification of cultivation in the surrounding area of the reservoir, site-specific soil management and careful land use planning are needed to take into account the spatial variability of soil properties, delineating management zones, variable fertilization management, irrigation scheduling, conservation practices and other efforts.

Figure 1 .
Figure 1.Location of the study area at the Alqueva dam watershed in Portugal.

Figure 2 .
Figure 2. Three-dimensional perspective of the trends in the input data sets.

Figure 4 .
Figure 4.The HJ-Biplot representation matrix of soil samples and studied variables.

Figure 5 .
Figure 5. Hierarchical clusters representation of soil samples and studied variables.

Table 1 .
Descriptive statistics of soil properties and parameters of the fitted variogram models and the cross-validation results.