In order to have an efficient use of natural resources, understanding soil properties, as well as its geographical distribution, mapping soil attributes becomes an important tool in the development of an adequate plan for soil management (Dessalegn et al.,  2014).  The study  of the soil properties spatial distribution is a key issue to support sustainable land management, applying erosion control, crop choice, and possibilities of irrigation (Van de Wauw et al., 2008).

The variability of the  soil  attributes  is  a  result  of  soilformation factors and processes interaction such as, the climate, the topography, the material origin and the vegetation (Rezaei and Gilkes, 2005a; Wang et al., 2009). According to the same authors, the interaction between the factors and processes occurs both in spatial and time scale, which can be altered locally mainly through soil erosion.

According to Gonçalves et al. (2001) the study of spatial and time variability is of extreme importance concerning the attributes that influence the soil water storage, such as: soil depth, infiltration capacity, topography and the regional climate. Li et al. (2002) mentions that the structural quality of the soil has been associated with the favorable conditions for the growth of the root system, aeration, infiltration and water movement in its profile. The attributes that present variability induced by management practices are, among others, the surface layer thickness, structure, bulk density, porosity and hydraulic conductivity of the saturated soil, are more dynamic in space and time (Rezaei and Gilkes, 2005b).

In this respect, the scientific knowledge advance in spatial variability studies of soil properties started to include, in addition to the classic statistic, the geoestatistics. The geoestatistics incorporates the spatial and temporal coordinates in the data process (Goovaerts, 1999; Webster and Oliver, 2007), presenting techniques that allows precise estimations of the studied variations (Scolforo et al., 2015). Among the geoestatistics methods utilized in spatial variability studies of soil attributes, kriging methodologies can be highlighted (Sun et al., 2003; Umali et al., 2012).

It should be noted that the studies carried out about the spatial variability of the physic and hydraulic properties are mostly applied in small scale areas, being inadequate or insufficient for the use in soil management in the catchment scale (Dessalegn et al., 2014; Tesfahunegn et al., 2011; Wang and Shao, 2013).

The soil properties are connected, having an important effect in the hydrology of the watershed, determining the water infiltration capacity into the soil, the groundwater flow, affecting the surface runoff and the baseflow (Price et al., 2010; Wang and Shao, 2013). Therefore, the comprehension of the attribute of soil behavior becomes relevant to help in the natural resources conservation, sustainable uses in watersheds and in the improvement of management practices of agriculture, as their effects on the environment (Cambardella et al., 1994; Pinto et al., 2016). According to Mello et al. (2011), the greatest concern in the studies in basin scale is in the comprehension on the environmental balance and the impact of the use and management of the soil in the process of runoff and erosion.

In this context, this study was developed with the objective of evaluating the spatial variability of physical and hydraulic soil attributes in a watershed, located in the Sul-Rio-Grandense Shield in the South of Brazil through the geoestatistic analysis.  

The Pelotas Stream watershed is located in the South of the State of Rio Grande do Sul - Brazil. It is an important water source for the region, once it is the main water public source in the city of Pelotas, where it is estimated around 342.873 residents (IBGE, 2015).  The watershed is located in two geomorphologic provinces (Sul-Rio-Grandense Shield and Planície Costeira), the first is located in the higher altitude of the basin and the second in the portion near from the coast.

The study was carried out in the do Ouro watershed, a subbasin of the Pelotas watershed. The do Ouro watershed has an area of 17.17 km², located between the UTM coordinates 352243.02 and 346693.81 longitude E and 6506001.84 and 6500135.29 latitude S, with altitude ranging between 76 and 326 m (Figure 1A).

The regional climate, according to Köppen climate classification, is type “Cfa”, that is, humid temperate with hot summers. The annual average temperature, precipitation and relative humidity in the region are 17.5°C, 1276 mm and 83%, respectively.

The soils in the watershed of do Ouro Stream were classified as Argissolos and Neossolos based on the survey by Cunha et al. (1996) similar to Acrisols and Regosols as classified by World Reference Base WRB (FAO, 2014). The topography is mainly undulating (46.9% of the total area), followed by mildly undulating and strongly undulating with 21.7 and 19.4% of the area, respectively (Figure 1B). Regarding the soil use, according to Bartels (2015), native grass is de main land cover in 43.4% of the total area, being destined mainly for beef consumption and dairy cattle. The other areas as covered by native forests (34.1%), annual crop (13.9%), cultivated forests (4.2%) and orchards (3.1%). The soil samples were collected in 3 transects according to Figure 1C.

Each transect presents the same distance of 1500 m between each other, containing spaced points every of 150 m, totaling 67 samples. Thus, the sampled points represented the multiple uses present in the basin, as well as its different physiographical characteristics, such as slope and altitude. The soil samples were collected in the surface layer (0-0.15 m) and subsurface layer (0.15-0.25 m). Deformed samples were collected for particle size distribution analysis. Samples with preserved structure were also collected for determining the macroporosity (Ma), microporosity (Mi), total porosity (TP), soil bulk density (BD), hydraulic conductivity (K) and Field Capacity (qFC). This soil sampling was performed with Uhland soil sampler in cylinders with a diameter and height of 7.6 cm (approximately 344.1 cm³).

Afterward, soil samples were taken to Laboratory of Soils and Hydro-sedimentology of the Water Resources Engineering course - UFPel, in plastic bags previously identified. The particle size distribution was determined by the pipette method, according to what was proposed by Embrapa (1997), using chemical dispersion with NaOH solution 0.1 mol L^-1 and mechanical agitation with high rotation for 15 min. The samples with preserved structures collected in the cylinders were saturated for 48 h through gradual water depth elevation.

The total porosity, macroporosity, microporosity and soil bulk density were determined according to methodology proposed in Embrapa (1997), with the TP determined by the difference between the saturated soil mass and the dried soil at 105°C for 24 h. For determining the microporosity, the tension table method was used, in which a -6 kPa tension was applied. The macroporosity was calculated by the difference between the total porosity and the microporosity.

The soil water contents at field capacity (θFc) were taken as those corresponding to potential at -10 kPa (Reichardt and Timm, 2012). After sample dried out at 105°C for 24 h, the soil bulk density was determined by the relation between the mass of dry soil and the cylinder volume. The hydraulic conductivity of the saturated soil was determined through a constant  head  permeameter,  according to Embrapa (1997). Soil water content at permanent wilting point (qPWP) corresponding to potentials of -1500 kPa (Reichardt and Timm, 2012), was determined in previously air-dried and sieved (2 mm sieve) soil samples, with a psychrometer (Gubiani et al., 2013).

Firstly, the obtained data were submitted to a descriptive analysis using the statistic software Assistat (Silva and Azevedo, 2009), where average, minimum, maximum, standard deviation (S), asymmetry, kurtosis, coefficient of variation (CV) and normality test, were obtained. In addition to this, boxplot graphics were analyzed for removing outliers. The coefficient of variation (CV) was classified, according to Wilding and Drees (1983) as: CV ≤ 15%; 15% < CV ≤ 35%; CV > 35%, as low, medium and high, respectively. For testing the hypothesis of normality of data distribution, the Shapiro and Wilk (1965) test was applied at a level of 5%.

With Geoest package (Vieira et al., 2002), the geostatistical analysis was carried out to characterize the spatial variability of the physical and hydraulic properties. Therefore, the calculation of an experimental and theoretical semivariogram and its respective adjustment parameters was performed (Equation 1).

Wherein: N(h) is the number of observation pairs Z(xi) and Z(xi+h) separated by a distance h. A theoretical model was adjusted  to  the experimental semivariogram, where C0 = nugget effect; C0+C1 = sill; A = range.

All the models of semivariograms were submitted to validation by the “Jack-Knifing” method (Vieira et al., 2002). The spatial dependence degree (SDD) was classified according to Zimback (2001), as: spatial dependence ≤ 25%; 25% < spatial dependence ≤ 75% and spatial dependence > 75%, in low, moderate and high, respectively.

The values of asymmetry show an asymmetric distribution of the soil attributes (Table 1). Values of positive asymmetry show a tendency of concentration of values below the observed average, being this tendency more significant the higher the obtained value. An inverse situation occurs for negative values (Neves Neto et al., 2013). However the asymmetry in the data is not significant, once the asymmetry coefficients are inside the ±1 limit proposed by Kerry and Oliver (2007).

According to the classification of the coefficient of variation (CV) proposed by Wilding and Drees (1983), sand, BD and TP presented low variability (CV < 15%) for the two depths.  It can be noted in Table 1 that the K was the attribute that showed the highest variability on both analyzed layers, followed by macroporosity. The high values for CV represent heterogeneity in relation to the average. The other attributes presented a moderate variability, with CV between 15 and 35%.

The results of the Shapiro and Wilk normality test (p≤0.05) indicate that only the variables K in the surface layer, and clay, Ma, K and (θPWP) in the subsurface layer, do not present normality. According to Paz-Gonzales et al. (2001), the kriging presents better results when the data normality is observed. However, the set of data normality is not a mandatory requirement for applying geoestatistics (Cressie, 1993).

The mathematical models of semivariograms are presented, as well as, its respective parameters of adjustment for the physical and hydraulic soil attribute in the 0 to 0.15 and 0.15 to 0.25 m layers (Table 2). On both layers, the variables sand, silt, Ma, Mi, TP, K and θPWP showed a better adjustment to the exponential model and the variable clay to the spherical model. The exponential model also adjusted the best to the variables BD and Field Capacity (θFc) in the 0.15 to 0.25 m layer. These same variables had a best fit to the spherical model in the surface layer (0-0.15 m). Mcbratney and Webster (1986) affirm that the spherical and exponential models are the ones most frequently adjusted to the soil attributes. 

It was also observed that the range of spatial dependence found in 0 to 0.15 m depth presented the following values for the soil physical and hydraulic attributes: Sand, 700 m, clay, 600 m, Ma, 800 m, Mi, 850 m, TP, 548 m; BD, 480 m; K, 650 m, Fc, 560 and PWP, 670 m. For 0.15 to 0.25 m depth, however, the range (A) presented the following values: Sand, 870 m; silt, 840 m; Ma,  885 m;  TP,  820 m,  BD,  700 m,   θFc   700 m   and θPWP, 680 m. According to Cora et al. (2004) the range values can influence the quality of estimates, as it determines the number of values used in the interpolation. Therefore, the estimates done with interpolation by ordinary kriging using values of higher range tend to be more reliable, presenting maps that better represent reality.

Physical and hydraulic soil attributes studied presented spatial dependence, as none has pure nugget effect (Table 2). The level of spatial dependence, according to the classification proposed by Zimback (2001), indicates that the variables sand, silt, Mi, TP, BD, K, θFc e θPWP in both layers presented spatial dependence classified as moderate (25% < spatial dependence ≤ 75%). The values of spatial dependence were high (spatial dependence > 75%) for Ma in the 0 to 0.15 m layer and clay in the 0.15 to 0.25 m layer (Table 2). The variations observed in the spatial dependence level of the physical and hydraulic soil attributes can be influenced by intrinsic factors, that is, soil formation factors, source material, relief, climate and organisms, and by extrinsic factors, that is, management practices such as fertilizing, plowing, harrowing, liming, etc. Figure 2 presents the spatial distribution maps of the attributes studied  in  the  surface and subsurface layer.

Hydraulic conductivity of saturated soil maps present very similar pattern values with the maps of the variables macroporosity and sand, in both layers. Thus, it can be inferred that soils with higher sand content present higher quantity of macropores and as a consequence increase in water flow in the soil profile. The same cannot be stated about microporosity and soil density as they present an inverse relation to the saturated hydraulic conductivity, in both layers.

It can be concluded about the spatial variability of soil physical and hydraulic properites in the do Ouro watershed that: (i) The variables K and macroporosity presented a high heterogeneity in relation to the average, with K in the surface layer and the variables clay, Ma, K and θ_PWP in the subsurface layer not presenting normality according to Shapiro and Wilk test (p≤0.05). (ii) The theoretical model of the exponential semivariogram is the one that best describes the physical and hydraulic attributes in spatial variability structure. (iii) Spatial dependence was observed for all the soil physical and hydraulic attributes studied, classifying as high for the variables macroporosity (0-0.15 m) and clay (0.15-0.25 m), with the other variables classified as moderate spatial dependence. With the maps provided by the present study, it is possible to distinguish areas in the basin according to their soil and water attributes, enabling the development of conservation practices. However, some improvements in the map creation must be highlighted. Future studies may consider more sample points in order to improve the interpolation analysis. Moreover, a higher spatial density of soil samples is also recommended.

The authors have not declared any conflict of interests.

The authors wish to  thank  CAPES  for  the  postdoctoral fellowship to the second author, and for the postgraduate scholarship to the first author and CNPq for the scholarship to the fourth author.