Spatial variability of soil physical and hydraulic properties in the southern Brazil small watershed

In order to provide an adequate usage to agricultural soils it is crucial to comprehend the soil physical attributes and their spatial variability. The variability of soil attributes is dependent on several formation processes and interactions. To properly describe such complex variability, the statistical methods used must incorporate the spatial and temporal influences. In this context, the present study was developed with the objective of evaluating the spatial variability of physical and hydraulic soil attributes, such as sand, clay and silt contents, soil macroporosity (Ma), soil microporosity (Mi), soil total porosity (TP), soil bulk density (BD), hydraulic conductivity (K), water contents at field capacity (tension of 10 kPa) and permanent wilting point (tension of 1500 kPa) (θFc and θPWP), in a watershed, located in the Sul-RioGrandense Shield, in the South of Brazil, through the geoestatistics analysis based on Ordinary Kriging. It was found a better performance of exponential model for the semivariograms of most physical and hydraulic soil attribute in the 0-0.15 and 0.15 to 0.25 m layers. Spatial dependence was observed for all the soil physical and hydraulic attributes studied, classifying as high for the variables macroporosity (0 to 0.15 m) and clay (0.15 to 0.25 m), with the other variables classified as moderate spatial dependence. Besides, 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).

formation 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.

MATERIALS AND METHODS
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.02and 346693.81longitude 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.15m) 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 (FC).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 (PWP) 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).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.

RESULTS AND DISCUSSION
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.15m).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.
The lower concentrations of sand are located in the southern part of the basin in both evaluated layers (Figure 2).The opposite occurs with the spatial distribution of clay contents.The spatial distribution of the variable soil bulk density in the layer 0.15 to 0.25 m is more heterogeneous when compared to its distribution in the surface layer.The higher values of BD can be considered normal throughout do Ouro watershed for the sand contents in the watershed.Reinert et al. (2008) analyzed the BD effect in the root growth in a soil with grain size similar and the observed in the study area.A normal root growth for bulk densities up to 1.75 kg dm -3 in the experimental area, with the highest values found in the northern part where the higher concentration of sand was obtained.The total porosity and macroporosity maps show a opposite behavior in the 0.15 to 0.25 m layer, in other words, areas that present the lowest BD values have the highest TP values and Ma.The compaction caused by animal trampling can be linked to the increase in the soil BD.This basin is mainly covered by native grass which is destined to dairy cattle, which according to Menezes et al. (2016) is the main cause for the increase in soil BD followed by a decrease in TP in Oxisols.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.

Conclusion
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.15m) 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.), TP (m 3 m -3 ), Ma (m 3 m -3 ), Mi (m 3 m -3 ), BD (kg dm -3 ), K (mm h -1 ), θFc (m 3 m -3 ) and θPWP (m 3 m -3 ), in the 0-0.15 and 0.15-0.25 m layers, in the do Ouro watershed.fellowship to the second author, and for the postgraduate scholarship to the first author and CNPq for the scholarship to the fourth author.

Figure 1 .
Figure 1.Location of watershed in do Ouro Stream (A), map of slope (B) and soil sampling points (C).
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.

Table 1 .
Parameters of descriptive statistic for physical and hydraulic attributes of the soil.
Sand, Clay and Silt, sand, clay and silt contents; Ma, soil macroporosity; Mi, soil microporosity; TP, soil total porosity; BD, soil bulk density; K, hydraulic conductivity; θFc and θPWP, soil water contents at field capacity (tension of 10 kPa) and permanent wilting point (tension of 1500 kPa), respectively; SD, Sample Standard Deviation, CV, Coefficient of Variation (%), Cs, Coefficient of Asymmetry, Ck, Coefficient of Kurtosis, SW, Shapiro and Wilk Test, Significance of 5%, N, Follows the Normal Distribution NN, Does not follow the Normal Distribution; * Surface layer; ** Subsurface layer.

Table 2 .
Theoretical models of semivariograms and respective parameters of adjustment of physical and hydraulic soil attribute in the 0-0.15 and 0.15-0.25 m layers.