Kriging Versus Cokriging and Collokated Cokriging for Soil Physical-Hydraulic Attributes and their Influence on Soybean Growth

Abstract

in Brazil, management agricultural practices not currently consider the soil spatial variability as a result, crop growth can be non-uniform and yields often is low. This research aims to compare Kriging, Cokriging and Collocated cokriging using soil physical and hydraulic properties and their influences on soybean development. We hypothesized that spatial variability of physical and hydraulic properties has influence on soybean development and this variability can be better represented by Collocated Cokriging method. To test these hypotheses, we accessed the soil physical and hydraulic attributes in a field experiment under no-till system, cultivated with soybean. Geostatistical interpolators were applied to generate maps from which spatial dependence of the variables was evaluated. The experiment was conducted on a sandy clay loam Oxisol, on an experimental station located in Ponta Grossa, Paraná, Brazil. Evaluation of the soil attributes was performed: bulk density (BD), particle size distribution, saturated soil hydraulic conductivity (K _fs ), total porosity (TP), macroporosity and microporosity. The plant was plant height and stand. Data analysis were performed by geostatistical methods; the spatial dependence was established using experimental univariate and cross semivariograms with datasets. Modeling semivariograms led to the generation of attribute maps by Kriging, Cokriging and Collocated cokriging. The estimation by Cokriging and Collocated cokriging was similar from Kriging. From the semivariogram, it was possible to identify that soil and plant attributes were spatially related with each other. The soya growth was mainly changed by slope of the area and little changed by saturated hydraulic conductivity.

Keywords: soil water; spatial variability; interpolation methods

HIGHLIGHTS

The soil spatial variability is related with plant development.

There was no difference between the interpolation methods.

The hydraulic conductivity and mainly the slope change the plant growth.

INTRODUCTION

The inaccurate estimation of crop growth is partially due to spatial variability of soil properties, which is caused by the complex formation of soil process, topography, biotic factors and management (Yan and coauthors) [¹]. Spatial variability of soil can result in big variation of soil physical and hydraulic properties. In addition, these soil physical properties have a great role in crop growth (Duarte and coauthors [²]; Strudley and coauthors [³]).

The crop growth in agricultural systems can be altered by a variety of factors, such as physical quality. Poor soil physical quality occurs when soils exhibit one or more of the following symptoms: poor water infiltration, run-off of water from the surface, hard-setting, poor aeration, poor rootability, and poor workability. Good soil physical quality occurs when soils exhibit the opposite or the absence of the conditions listed above (4).

One of the most sensible parameters to detect the soil heterogeneity in the area is soil hydraulic conductivity, defined as the soil’s ability to transmit water (Klute and Dirksen, 1986) [⁵]. The saturated soil hydraulic conductivity (K _fs ), when compared with soil physical attributes, such as soil bulk density and texture is highly variable mainly due to the variation of the macroporosity - that freely drain water under the influence of gravity -contributing with more than 70 percent of the total soil water infiltration (Watson and Luxmoore, 1986 [⁶]). These attributes are central to understanding the soil-plant system from the perspective of soil water dynamics (Klute and Dirksen, 1986) [⁵].

The link between plant development within a specific management system, and the spatial dependence is possible through geostatistics. Schaffrath and coauthors [⁷] found that spatial correlations of soil physical attributes were more evident in no-tillage systems than in conventional systems, leading them to conclude that soil management have a significant influence on the spatial dependence of soil properties.

Geostatistics allows to study the spatial variability in agricultural systems. Given these capabilities, geostatistics has been successfully applied in Soil Science. The analytical procedures through which geostatistics reveals spatial patterns and defines relationships between variables relies on fitting models and the use of specific interpolators. An incorrect choice of interpolator could result in lower precision in the mapping of an area, which in turn could result in the wrong interpretation of spatial dependence of the variable, and may also lead to problems in establishing the relationships between the different mapped variables of the system (Schaffrath and coauthors, 2015 [⁷]).

In geostatistics, the examination of cross semivariograms is used when many variables are sampled in the same space (Schaffrath and coauthors, 2015 [⁷]). Interpolators like Kriging utilizes a semivariogram to estimate a variable unsampled location using weighted neighboring measured values (Nielsen and Wendroth, 2003 [⁹]), Cokriging, utilizing spatial information of two or more variables (one is the primary variable, the other is the auxiliary variable) along with spatial cross-correlation to estimate the unsampled variable and Collocated cokriging make the linear combination between main variable (Z1) and one second variable (Z2) (Xu and coauthors, 1992, [¹⁰]). According to Journel and Huijbregts (1978) [¹¹], the major contribution of cokriging is the co-estimating poorly sampled variables, and the reduction of error-variance estimation and the opportunity to estimating many attributes in the same domain (Myers, 1982 [¹²]).

Collocated cokriging join the benefit to estimate the second variable, in addition is an attempt to avoid matrix instability that comes out from the application of ordinary cokriging (Goovaerts, 1997 [¹³]). This problem is a direct consequence of a high auto-correlation between closer secondary data as opposed to lower auto-correlation and more distant primary data, resulting in negative co-estimates (Xu and coauthors 1992 [¹⁰]). According to Xu and coauthors, another reason for choosing collocated cokriging is the need of eliminating redundant information since secondary data collocated or located near unknown primary data tend to screen the influence of further away secondary data [¹⁰].

Studies focused on the spatial variability of soil properties and comparing the performance of kriging and Cokriging had been developed. However, there are few studies including Collokated cokriging. This manuscript wanted to advance the investigation of spatial variability of soil physical and hydraulic properties with geostatistics analyzing the effects of spatial variability on soil physical and hydraulic properties and their influences on soybean plant development by Kriging, Cokriging and Collokated cokriging.

MATERIAL AND METHODS

Study area

Field data acquisition was obtained in an experimental farm located in Ponta Grossa, within the state of Paraná, Brazil (25º 24' 56" S and 50º 06' 14" W) (Figure 1). The climate of the region is described as mesothermal, subtropical and moist, being a type humid temperate climate with moderaty hot summer Cfb within the Köppen classification scheme. The local climate is characterized by mild summers, hard frosts in the winter, and the absence of a dry season. Average annual rainfall is 1545 mm and the mean temperature is 18.7 °C (IAPAR, 1978) [¹⁵].

The area selected for the study is situated in the middle third of an area with low slope; the relief is convex, with a little slope towards the South of 10 percent. The altitude varies between 1012 and 1022 meter above sea level. The soil has been classified as a Oxisol clay soil according to the USA Soil Taxonomy (Soil Survey Staff, 2014 [¹⁷]). The area has been cultivated under no tillage for approximately 20 years in accordance with the following crop rotation: maize (Zea mays L.) and soybean (Glycine max L.) during the spring/summer followed by wheat (Triticum aestivum L.) intercropped with black oat (Avena strigosa Schreb.) þ vetch (Vicia sativa L.) during the autumn/ winter. When the field experiments were established, the experimental areas were covered with crop residue after the maize harvest, and cultivated with soybean, cultivar Nidera 5909.

Sampling scheme

A regular square grid of sampling points, with intervals of 10 m ( 10 m in the X and Y directions, was set up within the experimental field over a rectangular area of 70 m by 110 m, resulting in a total of 96 sampling points. At each of the sampling points, the altitude was measured with a GPS system; altitude constituted one of the variables within the subsequent analysis of the field data (Figure 2).

In the sampling grid, points located between lines of the soybean plants, undisturbed soil samples were collected from 0-5 cm depth, utilizing a soil sampler of the Uhland type. This sampler was coupled with ring with a cross-section of 5 cm, with a volume of approximately 100 cm³.

The soil samples were saturated by capillarity, placed on tension table, and submitted to a tension of -6 kPa. After the establishment of hydraulic equilibrium, the samples were weighed, dried in an oven at 105 ºC and weighed again to determine the mass of dry soil from which the soil water content under a tension of -6 kPa was obtained by difference. The total porosity (TP) was calculated from the soil bulk density (BD) and the particle density (PD) using the value 2,65 g cm^-3 (Flint and Flint, 2002 [¹⁶]) using the following equation:

$TP\left(\%\right)=1-\frac{BD}{PD}*100$ (1)

The macroporosity was calculated as the difference between TP and from the quantity of water retained in the soil at -6 kPa), according to the methodology suggested by Lorraine and coauthors 2002 [¹⁸]. Disturbed samples were also collected from between the rows of plants at each of the 96 grid points and submitted for particle size analysis (sand, silt, and clay). The particle size analysis was performed using the hydrometer method (Gee and coauthors,2002 [¹⁹]).

The field saturated soil hydraulic conductivity (K _fs ) was determined in the field by a Simplified Falling Head method (SFH) (Bagarello and coauthors 2004 [²⁰]). The measurements were conducted between the rows of soybean plants, using a steel cylinder with dimensions of 21.2 cm (diameter) ( 9.4 cm (height); the area of the base (A) was 353 cm². The measurement site was prepared by removing any plant straw from the soil surface, the steel cylinder was inserted into the surface of the exposed soil to a depth of 5 cm. A fixed volume of water (V = 0.35 L) was evenly applied to the soil across the entire area encircled by the steel ring. The infiltration time (ta) was recorded from the application of the water to the soil to the moment at which water was no longer visible on the soil surface. The dry soil sample was collected from an area adjacent to the steel ring and a sample of saturated soil was collected from inside the ring after the complete infiltration of the water into the soil. The difference in the volumetric water content ((() between the saturated soil (sample collected from inside the ring after infiltration) and non-saturated soil (dry soil sample collected from outside the ring) was obtained from measurements of volumetric moisture. Following the acquisition of the raw data, K _fs was calculated according to the following equation provided by (Bagarello and coauthors 2004 [²⁰]):

$K_{fs}=\left[\frac{\Delta \theta }{\left(1-\Delta \theta \right)ta}\right]\times \left[\frac{D}{\Delta \theta }-\left\{\frac{D+\left(1/\alpha *\right)}{\left(1-\Delta \theta \right)}\right\}\times \ln \left\{1+\frac{\left(1-\Delta \theta \right)D}{\Delta \theta \left(D+\left(1/\alpha *\right)\right)}\right\}\right]$ (2)

In Equation 2, $D=V/A$ is the ratio of the volume of water applied to the area of the steel ring inserted into the soil. This ratio determines the depth of the water in the ring above the soil at the start of the measurement of the infiltration time ta. (* is a constant which depends on the texture and structure of the soil; (* varies between 4 and 36 m^-1 (Elrick and Reynolds, 1993 [²¹]). Within the experimental area the texture of the soil was approximately sand 74 %, silt 5 %, clay 21 %, and a value of 12 m^-1 was accordingly adopted for the constant (*.

The following plant attributes were measured at each of the sampling points: plant height, reproductive stage, and stand (number of plants per unit length along the rows). In order to gather the plant height and reproductive stage data, six plants were evaluated at each sampling point, three in the row left the sampling point and three in the row right in east direction. The height measurements were performed with a graduated ruler, recorded in centimeters, measuring from soil level up to the highest totally expanded leaf. The reproductive stage was assessed according to the procedures of Fehr and Caviness (1977) [²²], in which visual observations were made of the flowers on the last nodes on the stem with a fully developed leaf. The recorded value for each plant was the mean of the six measurements obtained at each point. The value for the stand was obtained by manually counting the number of plants per meter along a plant row.

Data analysis

The results were submitted to descriptive statistical analyses (Vieira and coauthors 2002 [²³]) which included testing of the normality of the data through the Kolmogorov Smirnov test at alpha levels of 0.01 and 0.05. In the case of the K_fs data, a transformation to the logarithm of K_fs was required to obtain the assumption of normality.

Geostatistical analyses (Vieira and coauthors 2002 [²³]) were employed to verify if the variables displayed spatial dependence, and then to perform interpolation and build maps through the methods of ordinary Kriging, Cokriging and collocated Cokriging. Semivariograms were calculated using Avario program for each variable using estimator and examined evidence of spatial dependence and trend. The data for altitude, plant height, K_fs, sand, clay, microporosity and BD, displayed no trend and presented stationarity. The experimental semivariograms were represented by a Gaussian model for the following variables: altitude, plant height, K_fs, sand, clay, and microporosity. The semivariogram for BD was described by a spherical model.

The procedures for these analyses are detailed in (Vieira and coauthors 2002 [²³]). The nugget effect ($C_0$) and the sill (C₀+C₁) for each variable were obtained after fitting the semivariance models to the experimental semivariograms and used to calculate a ratio Standard deviation ratio (SDR):

$SDR=\left(\frac{C_0}{C_0+C_1}\right)\times 100$ (3)

which is a measure of the spatial dependence of the variable and was used to classify the dependence as weak (SDR > 75 %), moderate (25 % ( SDR ( 75 %), or strong (SDR < 25 %) (Cambardella and coauthors 1994 [²⁴]). For the variables which presented a spatial dependence, the values adjacent to the sampling points were estimated by Kriging. For the variables between which there was spatial correlation, the cross semivariograms (examples are presented in Figure 5) were fit to the Gaussian model of equation, and the estimation of the values of the correlated variables at points adjacent to the sampling points could be performed by Cokriging and Collocated cokriging (Vieira and coauthors 2002 [²³])

These methods provide estimated values without biases and with minimum variance in relation to the known values at the sampling points. Once these interpolation procedures for the estimation of the values of variables at any point within the sampling grid were set up, it was possible to produce contour (isoline) maps over the sampled area for each interpolator (Kriging, Cokriging and collocated Cokriging) and variables utilizing the program Surfer 7.0 (Golden Software, 1999 [²⁵]).

The quality of the estimation methods (Kriging, Cokriging and collocated Cokriging) was evaluated through two approaches. First, from the variance of the estimation, which is extremely sensitive to the form of the semivariogram or cross semivariogram (Vieira and coauthors 2002 [²³]) Second, through the normalized root mean square error (RMSE), which measures the error of predictions from the model (Fortin and coauthors 2014 [²⁶]).

Analyses of correlation were performed to identify relationships between variables selected as dependent K _fs , plant height, BD and a set of potential explanatory variables (easily measured soil physical attributes, including texture). The significance of the explanatory variables in predicting the dependent variables was analyzed by performing an F-test utilizing the ‘Biotools’ package (Silva, 2015 [²⁷]) of the R statistical software (R Core Team, 2014 [²⁸]).

RESULTS

Descriptive Statistics

The descriptive statistics of the 96 recorded values for each of the measured variables are presented in Table 1. Of particular interest is the coefficient of variation (CV), which is a measure of the values dispersion of the vales across the different sampling points in relation to the mean value. Following the classification scheme of Warrick and Nielsen (1980) [²⁹], the CV was found to be low for altitude, plant height, sand, clay, bulk density, microporosity, and total porosity, with CV values ranging from 0.31 % (altitude) to 12 % (clay and microporosity), medium for stand, log(K _fs ) and macroporosity, varying between 17 % and 34 %, and high for K _fs (CV = 98 %). Plant height, bulk density, and total porosity presented normal distributions, while the distributions for altitude, plant stand, sand, clay, macroporosity, microporosity and log(K _fs ), showed a significant asymmetry (Table 2).

Spatial Dependence

The experimental semivariograms for each of the variables in Table 1 were examined to identify those displaying spatial dependence. In the case of the macroporosity, the total porosity, and the plant stand, the semivariograms showed a semivariance which did not vary with the distance h. Therefore, these attributes displayed only a nugget effect, indicative of a total absence of spatial dependence (Journel and Huijbregts, 1978 [¹¹]).

The semivariograms demonstrated the existence of spatial dependence, these semivariograms were fit to the Gaussian model, and spherical model to the bulk density (BD). The results obtained from the fitted models are presented in Table 2, including the conclusions from the spatial dependence ratio (SDR), defined by Equation (3). While a strong spatial dependence ratio (SDR) was found for the clay soil content, there was a moderate SDR for altitude, plant height, sand content, bulk density (BD), and microporosity, and a low SDR for K _fs . However, if we analyze the range of the K _fs (40 meters) is higher than of bulk density (27m) and sand (37m) (Table 2, a) but these two soil attributes have SDR classified as moderated, while K _fs is weak

Starting from the fitted semivariogram models of the spatially dependent variables in Table 2, a contour map of each variable was generated by ordinary Kriging. The resulting maps for altitude, plant height, and the following soil variables: saturated hydraulic conductivity (log [K _fs ]), clay content, microporosity, and bulk density (BD) are in Figure 3. Examination of the maps revealed that the saturated hydraulic conductivity was highest in areas of greater altitude (corresponding to areas of smaller slope), and plant development (height) was also varied from 71.5 to 79 cm in the areas with smaller slope, while in the areas of greater slope the plant height range was 64 to 69 cm. Although the height of the plants has low correlation with K _fs , in some areas of the map, the low K _fs are associated with low plant development. This difference in plant height of 15 cm may have an influence on subsequent crop production, while its correlation with a difference in K _fs suggests that this parameter should be considered in crop management practices.

Correlations between the variables

Spearman correlation coefficients were calculated between all pairs of the following variables: altitude, plant height, soil bulk density (BD), log (K _fs ), sand, silt, and clay content of the soil (Figure 4).

The positive correlation established between the soybean vegetative development and soil hydraulic conductivity (Spearman correlation coefficient between plant height and log (K _fs ) = 0.27, Figure 4) may be due to its macroporosity and sand content (Table 1) that are generally associated with high values for K _fs . Though the correlation between K _fs and plant height was low, the comparison of the contour maps for K _fs and plant height (Figure 3) demonstrated that in some areas the K _fs were associated with vegetative development of the soybean plants, with zones of high and low values for K _fs being associated with zones of higher and lower vegetative development, respectively. Soil bulk density have low correlation with the plant development (Figure 4)

Spatial relationships and estimator performance

An inspection of the contour maps of the attributes obtained by ordinary Kriging (such as those presented in Figure 3) suggests that many of the attributes displayed joint variations in space, which were revealed by constructing experimental cross semivariograms. The Gaussian model, were considered to give the best representations of the experimental cross semivariograms, as assessed from the variance and the root mean square error. The cross semivariograms confirmed the existence of high spatial correlations between the variables, which were already known to have high spatial dependence, with the magnitude of the nugget cross semivariance being close to zero and much less than the magnitude of the sill cross semivariance (Goovaerts, 1997 [¹³]). Four cross semivariograms are presented in Figure 5. The range of the cross semivariograms was of 50 m, indicating that there was correlation between the variables up to this distance.

When two variables present joint variation in space, the mapping of one of the variables using the data collected for the other is possible through Cokriging and collocated Cokriging, which use interpolators built from the modelled cross semivariogram between the two variables. Contour maps constructed by Cokriging and Collocated cokriging for K _fs , using the K _fs ( clay cross semivariogram, for sand, using the sand ( altitude cross semivariogram, and for clay, using the microporosity x clay cross semivariogram are presented in Figure 6, together with the ordinary Kriging contour maps for K _fs , sand, and clay. The mapping obtained utilizing the Cokriging and Collocated cokriging methods makes possible the visualization of the joint variation of the variables, from which inferences concerning the influences of management and slope may be reached.

The low influence of the soil texture on K _fs can be seen on the Figure 6 contour maps where the zones with high values for K _fs (log [K _fs / cm h^-1] 1.2 to 1.4) corresponded to high sand content (752 to 788 g kg^-1) and low clay (175 to 189 g kg^-1). However, soil texture was not the major factor in setting the zone of maximum K _fs , for within this area the sand content (716 to 752 g kg^-1) was less than the maximum and the clay content (177 to 201 g kg^-1) was greater than the minimum.

The highest clay concentration was found on the western side of the study area, in a zone where the altitude was between 1016 and 1020 m, between 2.35 and 6.19 m below the southwestern corner of the study area. In contrast, the highest sand content occurred towards the eastern side of the study area, in a zone of higher altitude (1020 to 1021 m), with a maximum altitude difference of slope from the southwestern corner of 1.68 m.

The values for the normalized root mean square error (NRMSE) and the variance by the different interpolation methods confirmed for all variables, Kriging was similar from those provided by Cokriging and Collocated cokriging (Figure 7). The use of a restricted neighborhood for multicollocated cokriging provided similar results than the utilization of kriging. Its means that a secondary information was highly satisfactory and leads to estimates that correspond closely with the natural behavior of soil physical and hydraulic properties.

Figure 1
A) Study area located in Ponta Grossa, Paraná, Brazil; B) Block diagram indicating the topography of the experimental area as a relief map: darker colors correspond to higher altitudes, lighter to lower altitudes.

Figure 2
A) Satellite image of the study area, Ponta Grossa, Paraná, Brazil. B) Sampling grid established within the study area.

Table 1
Statistical moments of the altitude, plant development attributes and physical and hydraulic properties of the soil from the measurements made at the 96 sampling points in the study area, located in Ponta Grossa, PR, Brazil.

Variable	Mean	Variance	²CV	Minimum	Maximum	³(	Skewness	Kurtosis
Altitude / m	1018	9.80	0.31	1011	1022	11	-0.72	-0.82
Plant height / cm	72.1	44.49	9.25	56.67	83.17	26.50	-0.29	-0.76
Stand / plants m-1	11.39	3.95	17.44	3	17	14	-1.29	3.99
¹log[K fs / cm h-1]	1.14	0.16	34.6	0.17	2.19	2.02	-0.35	-0.13
Sand / g kg-1	736.1	734.7	3.68	669.3	795.9	126.6	0.12	-0.38
Clay / g kg-1	206.8	604.6	11.89	149.7	280	130.3	0.18	-0.03
Bulk density/ g cm-3	1.55	0.01	5.08	1.37	1.72	0.35	-0.03	-0.42
Microporosity / %	27.06	10.49	11.97	16.73	33.96	17.23	-1.01	1.79
Total porosity / %	40.02	6.81	6.52	34.8	46.8	12.0	0.36	-0.29
Macroporosity / %	12.96	16.08	30.94	5.18	28.4	23.2	1.36	2.30

• ¹Kfs - field-saturated hydraulic conductivity; ²CV - coefficient of variation; ³( = (maximum - minimum)

Table 2
Parameters estimated from the experimental semivariograms for the following variables: altitude, plant height, field-saturated hydraulic conductivity of the soil (as log[K _fs ]), soil texture - sand and clay, soil bulk density (Ds), and the microporosity of the soil.

Variable	Model	¹C 0	²C 1	a	³r 2	4 RMSE	5 SDR
Altitude / m	Gaussian	0.02	0.05	50.73	0.79	0.18	Moderate
Plant height / cm	Gaussian	15	23.89	40	0.49	0.42	Moderate
log[K fs / cm h-1]	Gaussian	0.11	0.03	40	0.20	0.36	Weak
Sand / g kg-1	Gaussian	151.75	185.85	36.98	0.74	1.29	Moderate
Clay / g kg-1	Gaussian	148.36	520.23	45.00	0.77	1.60	Strong
Bulk density / g cm-3	Spherical	0.005	0.002	27.22	0.23	0.07	Moderate
Microporosity / %	Gaussian	6.70	5.03	44.36	0.61	0.06	Moderate

• ¹C ₀ - nugget effect; ²C ₁ - sill; a - range; r² - coefficient of determination; ⁴ RMSE - root mean square error; ⁵ SDR - spatial dependence ratio defined by equation (3) (SDR < 25 %: Strong; 25 % ( SDR ( 75 %: Moderate; SDR> 75 %: Weak).

Figure 3
Spatial variations in the altitude, soybean development (plant height), and soil properties: hydraulic conductivity (K _fs) , clay and sand content, microporosity, bulk density (BD) presented as contour maps generated from the raw data by ordinary Kriging.

Figure 4
Correlations between altitude, soil physical ad hydraulic properties, vegetative development of soybean plants. Larger circles indicate stronger correlation, red colored circle denotes a negative correlation, blue a positive correlation.

Figure 5
Cross semivariograms for sand and altitude, clay and microporosity, clay and sand, and K_fs and clay, the experimental cross semivariograms were represented by Gaussian model (equation 3). The adjusted parameter values are given with each cross semivariogram and the fitted model is plotted with experimental points.

Figure 6
Contour maps showing the spatial variation of the field-saturated hydraulic conductivity (log(K_fs )) of the soil (top row), the sand content (middle row) and the clay content (bottom row). The estimation methods which produced the maps were Kriging (left column), Cokriging (central column) and Collocated cokriging (right column). The cross semivariograms used to generate the Cokriging and collocated Cokriging maps were log(K_{fs x} clay for K_fs, sand _X altitude for sand, and clay x microporosity for clay.

Figure 7
Comparison of the performance of the Kriging (Kri), Cokriging (coKri) and collocated Cokriging (coloc) interpolators in the prediction of the selected variables at the sampling points. RMSE- root mean square error: Kfs- field saturated hydraulic conductivity.

DISCUSSION

Descriptive Statistics

The high variation occurring with K _fs , implies that this variable always changes at short distance.

Field-saturated hydraulic conductivity, K _fs , is known to be highly variable (Bagarello and Sgroi, 2004 [³⁰]; Keller and coauthors 2012 [³¹]). This arises because K _fs is strongly influenced by many factors that impinge upon the formation of the local soil including landscape and the soil biology, as well as soil management. Utilizing the same methodology for the measurement of K _fs , Keller and coauthors 2012 [³¹] has reported a CV value for log (K _fs ) of 34 % in fields where the soil texture varied from clay loam to clay. The low values of CV for sand, bulk density, and total porosity are like results previously reported by Keller and coauthors 2012 [³¹], who studied a soil like the present study.

The difference between the maximum and minimum values (( in Table 1) for the plant attributes (height and stand) and some soil properties can give an indication of how the physical and hydraulic properties influence on the plant development. Part of the influence can be attributed to altitude difference (11 m) may have contributed to the soil texture variation (((sand) = 126.6 g kg^-1; ((clay) = 130.3 g kg^-1).

Spatial Dependence

The altitude can influence the trajectories along which water drains over the ground; for this reason, local features of the landscape can reduce the spatial dependence (Sobieraj and coauthors 2002 [³²]). The weak spatial dependence (small SDR) found for log (K _fs ) can be attributed to no-tillage system. Alletto and Coquet 2009 [³³] has asserted that when a no-till system is implanted the consequent improvements in the soil attributes are not expressed homogeneously across the entire area, and as time passes, the variability in the soil may increase.

The moderate spatial dependence for BD in the present study was similar to that found in Cecília and coauthors 2017 [³⁴] in Oxisol under no tillage system, the same type of soil and system management used in our research. In addition, the low values for the nugget effect (Table 2, C₀) are an indication that the spatial variability of the soil was carefully detected, which may be describes to the small spacing (10 m) of the sampling grid.

The existence of spatial dependence on the variables, as well as the similarity of behavior between them, suggests that the spatial variability of soil properties be analyzed together, before adopting an experimental design at random, since any treatment adopted in this plot that requires homogeneity will lead to false results (Grego and Vieira, 2005).

The subareas with different levels of physical soil quality were demonstrated by the bulk density (BD) and the hydraulic conductivity (K _fs ) variable and may be the result of low slope of the area. The low slope in an area can affect the soil hydraulic properties, and the spatial and temporal dynamics of processes occurring in the soil. These processes control the movement of the chemical elements within the soil and consequently alter plant development ( Strudley and coauthors [³]).

Correlations between the variables

The positive correlation between altitude, K _fs and plant development (plant height) can be associated to the no tillage system that was adopted approximately 19 years before the study starts, it is possible that biological factors, and in particular the channels formed through the soil by ants and earthworms may have had a greater influence on the saturated conductivity of the soil (Blouin and coauthors 2013 [³⁶]), than the soil texture, which would account for the weak correlation found between log(K _fs ) and the soil textural attributes. In addition the direct association between K _fs and the soil physical quality, which arises because soil structure, macroporosity and soil structure determine the water flow (Vieira and Klein, 2007 [³⁷]).

The low correlation between BD and plant development have found by numerous studies that has proved that the increase in bulk density restrict plant development (Zarehaghi and coauthors 2017 [³⁸]). However, in the present study the highest value for the bulk density was 1.72 g cm^-3, found in two small regions within the area of higher sand content (Figure 3). According to Daddow and Warrington 1983 [³⁹], the upper limit values for bulk density in sandy loam soils, for the successful development of plants is between 1.70 and 1.80 g cm^-3. Across most of the study area, the soil was a sandy loam with an average density of 1.55 g cm^-3 (Table 1), which is below the critical value cited above, so the bulk density should not have altered the plants development.

Spatial relationships and estimator performance

The observations concerning the spatial distribution of the soil texture indicate that the altitude did influence the distribution across the study area. For example, erosive processes may have contributed to the transport of clay to lower parts of the landscape (Figure 6), as verified by Sadeghi and coauthors 2017 [⁴⁰]. In addition, it is possible that variations in the relief may have influenced the soil texture distribution across the study area, since the spatial distribution of the soil particle size distribution is related to the relief (Sobieraj and coauthors 2002 [³²]). Comparing kriging with collocated cokriging for water table mapping Xu and coauthors 1992 [¹⁰] found better representation for collocated cokriging.

Due to Cokriging and Collocated cokriging no require the determination of the second variable in all points, the use of these methods allow safe money and time in the procedure of field sample and laboratory analysis, but the statistical procedure is more complex. However, statistically kriging is the simplest method in comparison to Cokriging methods, but the numbers of samples is bigger, requiring more time and money in laboratory and field.

The mapping using the different interpolators can be compared by examining the maps within the rows of Figure 6. Although the maps obtained with the Collocated cokriging interpolator present greater detail, an established characteristic of maps generated by the Collocated cokriging method (Goovaerts, 1999 [¹³]), the behavior and distribution are like the maps generated by ordinary Kriging. Thus, it may be assumed that the variable map produced by ordinary Kriging is enough to evaluate the behavior of the variable across a study area, and offers a similar result to Cokriging and Collocated cokriging.

These results can contribute to soya farmers, especially in representative areas sloping in Brazil under No-tillage system, that like our results in these areas the soil physical and hydric can be variable in short space. The knowledge of this variability is so important to management system that which must agree with this variability.

CONCLUSION

From the semivariogram, it was possible to identify that soil and plant attributes were spatially related with each other. The slope of the area altered the soil texture and had a positive correlation with sand content. This alteration on sand content and hydraulic conductivity contributed low to the soybean development, the mainly contribution was the slope.

Differences were evident between the maps generated by the Kriging, Cokriging and Collocated cokriging interpolators, the Cokriging methods produced more detailed maps than ordinary Kriging, but this does not imply greater precision because the RMSE and variance was similar between methods. Therefore, according to objective, we can choose one of these methods. If the objective is safe time and money in analyses and sampling process, the Cokriging and Collocated cokriging interpolators are more recommendable but if the objective is facility in statistical procedure the Kriging is easier than Cokriging and Collocated cokriging.

Acknowledgments

I dedicate this study to my supervisor Alvaro Pires da Silva (in memory) for his precious guidance, support and contribution to this research. I thank the co-authors Sydney and Neyde for their precious contributions, the University of Ponta Grossa (UEPG), University of São Paulo (USP) for their support and CNPQ for their financial support.

REFERENCES

• ¹ Yan P, Peng H, Yan L, Lin K. Spatial Variability of Soil Physical Properties Based on GIS and Geo-Statistical Methods in the Red Beds of the Nanxiong Basin, China. Pol. J. Environ. Stud. 2019; 28(4): 2961-72. Available from: http://www.pjoes.com/Spatial-Variability-of-Soil-Physical-Properties-Based-non-GIS-and-Geo-Statistical,92245,0,2.html [Acessed 21th March 2021].
  » http://www.pjoes.com/Spatial-Variability-of-Soil-Physical-Properties-Based-non-GIS-and-Geo-Statistical,92245,0,2.html
• ² Duarte SJ, Glaser B, Lima RP, Cerri CEP. Chemical, physical and hydrological properties as affected by one year of Miscanthus biochar interaction with sandy and loamy tropical soils. Soil Syst. 2019; 3(24): 2-18. doi:10.3390/soilsystems3020024
  » https://doi.org/10.3390/soilsystems3020024
• ³ Strudley MW, Green TR, Ascough JC. Tillage effects on soil hydraulic properties in space and time: State of the science. Soil and Tillage Research, 2008; 99(1): 4-48. doi:10.1016/j.still.2008.01.007
  » https://doi.org/10.1016/j.still.2008.01.007
• ⁴ Dexter AR. Soil physical quality: Part I. Theory, effects of soil, texture, density, and organic matter, and effects on root growth. Geoderma. 2004; 120(3-4) 201-14. doi: 10.1016/j.geoderma.2003.09.004
  » https://doi.org/10.1016/j.geoderma.2003.09.004
• ⁵ Klute A, Dirksen C. Hydraulic conductivity and diffusivity. Laboratory methods. In: Klute A (eds), Methods of soil analysis - part 1. Physical and mineralogical methods. SSSA. 2 ed; 1986.p. 687-734.
• ⁶ Watson KW, Luxmoore RJ. Estimating macroporosity in a forest watershed by use of a tension infiltrometer. Soil Science Society of America Journal. 1986; 50(3): 578-82. doi:10.2136/sssaj1986.03615995005000030007x.
  » https://doi.org/10.2136/sssaj1986.03615995005000030007x.
• ⁷ Schaffrath VR, Andrade Goncalves AC, Sousa AJ, Tormena CA. Correlação espacial entre atributos físicos do solo e de plantas espontâneas em dois sistemas de manejo. Ver. Bras. de Ciência do Solo. 2015; 39(1): 279-88. doi:10.1590/01000683rbcs20150568.
  » https://doi.org/10.1590/01000683rbcs20150568.
• ⁸ Kim D, Hirmas DR, McEwan RW, Mueller TG, Park SJ, Šamonil P, Thompson JA, Wendroth O. Predicting the influence of multi-scale spatial autocorrelation on soil-landform modeling. Soil Science Society of America Journal. 2016; 80(2): 409-19. DOI:10.2136/sssaj2015.10.0370
  » https://doi.org/10.2136/sssaj2015.10.0370
• ⁹ Nielsen DR, Wendroth O.) Spatial and Temporal Statistics: Sampling Field Soils and Their Vegetation. 1 da ed. Catena Verlag: Reiskirchen; 2003. Geoscience. Available on: https://www.nhbs.com/spatial-and-temporal-statistics-book
  » https://www.nhbs.com/spatial-and-temporal-statistics-book
• ¹⁰ Xu W, Tran TT, Srivastava RM, Journel AG. Integrating Seismic Data in Reservoir Modelling: The Collocated cokriging Alternative. In: Annual Technical Conference of the SPE. Washington: Proceedings, paper, SPE 24741; 1992. p. 833 - 42.
• ¹¹ Journel AG, Huijbregts CJ Mining geostatistics. London: Academic Press; 1978
• ¹² Myers DE. Matrix Formulation of Cokriging Mathematical. Geology 14, 1982; 14(3): 249-257.Available on:https://www.researchgate.net/publication/232607996_1982_MyersDE_Matrix_Formulation_of_Cokriging_Mathematical_Geology_14_249-257 Accessed in: April, 2020.
  » https://www.researchgate.net/publication/232607996_1982_MyersDE_Matrix_Formulation_of_Cokriging_Mathematical_Geology_14_249-257
• ¹³ Goovaerts P. Geostatistics for natural resources evaluation. New York: Oxford; 1997
• ¹⁴ Kiani M, Hernandez-Ramirez G, Quideau SA. Spatial Variation of Soil Quality Indicators as a Function of Land-use and Topography. Canadian Journal of Soil Science. 2020; 100(4): 1500-11. doi:10.1139/cjss-2019-0163
  » https://doi.org/10.1139/cjss-2019-0163
• ¹⁵ Instituto Agronômico do Paraná-IAPAR. Cartas climáticas do estado do Paraná. Congresso e mostra de agroinformática. 1978; 1-6. Available on: https://infoagro.deinfo.uepg.br/artigos/pdf/info_102.pdf [Accessed on: April, 2020].
  » https://infoagro.deinfo.uepg.br/artigos/pdf/info_102.pdf
• ¹⁶ Flint AL, Flint LE. Particle Density. In Dane JH and Topp GC. (eds) Physical Methods, 3rd Edition, Madison: SSSA; 2002. 1555-60.
• ¹⁷ Soil Survey Staff. Keys to Soil Taxonomy, 12th ed. USDA-Natural Resources Conservation Service, Washington: DC; 2014
• ¹⁸ Lorraine E. Flint, Alan L. Flint. Porosity. In: Jacob H. Dane,G. Clarke Topp. 2002 Book Series: SSSA Book Series doi:10.2136/sssabookser5.4.c11.
  » https://doi.org/10.2136/sssabookser5.4.c11.
• ¹⁹ Gee GW, Dani OR. Particle Size Analysis. In: Dane JH, Topp GC. Methods of Soil Analisis. Part 4- Physical Methods. Madison: SSSA; 2002. p. 229-240.
• ²⁰ Bagarello V, Iovino M, Elrick D. A simplified falling-head technique for rapid determination of field-saturated hydraulic conductivity. Soil Sci. Soc. Am. J. 2004; 68(1): 66-73. doi:10.2136/sssaj2004.0066
  » https://doi.org/10.2136/sssaj2004.0066
• ²¹ Elrick DE, Reynolds WD, Methods for analyzing constant-head well permeameter data. Soil Sci. Soc. Am. J. 1993; 56 (1),320-323. doi:10.2136/sssaj1992.03615995005600010052x
  » https://doi.org/10.2136/sssaj1992.03615995005600010052x
• ²² Fehr WR, Caviness CE. Stages of soybean development. Ames: Iowa State University of Science and Technology;1977
• ²³ Vieira SR, Millete J, Topp GC. Handbook for geostatistics analysis of variability in soil and climate data Campinas: TOP; 2002
• ²⁴ Cambardella CATB, Moorman JM, Novak TB, Parkin DL, Karlen RF, Turco AE, Konopka. Field scale variability of soil properties in central Iowa soils. Soil science society of America journal. 1994; 58(5):1501-1511. doi:10.2136/sssaj1994.03615995005800050033x
  » https://doi.org/10.2136/sssaj1994.03615995005800050033x
• ²⁵ Golden Software Surfer 7.0: contouring and 3D surface mapping for scientists and engineers: user’s guide Available on: https://www.goldensoftware.com/products/surfer [Acessed on: August, 2020].
  » https://www.goldensoftware.com/products/surfer
• ²⁶ Fortin V, Abaza M, Anctil F, Turcotte R. Why Should Ensemble Spread Match the RMSE of the Ensemble Mean? Journal of Hydrometeorology. 2014;15(4): 1708-1713. doi:10.1175/jhm-d-14-0008.1
  » https://doi.org/10.1175/jhm-d-14-0008.1
• ²⁷ Silva AR. Tools for Biometry and Applied Statistics in Agricultural Science. R package version 2.2. Available on: http://CRAN.R-project.org/package=biotools [Accessed on April, 2020].
  » http://CRAN.R-project.org/package=biotools
• ²⁸ R Development Core Team. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing; 2014
• ²⁹ Warrick AW, Nielsen DR. Spatial variability of soil physical properties in the field. In: Hillel D (Ed.) Applications of Soil Physics. New York, Academic.1980; p.319-344.
• ³⁰ Bagarello V, Sgroi A. Using the single-ring infiltrometer method to detect temporal changes in surface soil field-saturated hydraulic conductivity. Soil Tillage Res. 2004; 76(1): 13-24. doi:10.1016/j.still.2003.08.008
  » https://doi.org/10.1016/j.still.2003.08.008
• ³¹ Keller T, Sutter JA, Nissen K, Rydberg T. Using field measurement of saturated soil hydraulic conductivity to detect low-yielding zones in three Swedish fields. Soil & tillage research. 2012; 124(1): 68-77. doi:10.1016/j.still.2012.05.002
  » https://doi.org/10.1016/j.still.2012.05.002
• ³² Sobieraj JA, Elsenbeer H, Cameron G Scale dependency in spatial patterns of saturated hydraulic conductivity. Catena. 2004; 55(1): 49-77. doi:10.1016/s0341-8162(03)00090-0
  » https://doi.org/10.1016/s0341-8162(03)00090-0
• ³³ Alletto L, Coquet Y. Temporal and spatial variability of soil bulk density and near-saturated hydraulic conductivity under two contrasted tillage management systems. 2009; Geoderma. 152(1): 85-94. doi:10.1016/j.geoderma.2009.05.023
  » https://doi.org/10.1016/j.geoderma.2009.05.023
• ³⁴ Cecília M, Da C, Andreotti M, Costa NR, Da G, Lima R, Pariz CM. Soil Physical Attributes and Yield of Winter Common Bean Crop Under a No - Till System in the Brazilian Cerrado. Rev. Catinga 2017; 30(1) 155-163. doi:10.1590/1983-21252017v30n117rc
  » https://doi.org/10.1590/1983-21252017v30n117rc
• ³⁵ Grego CR; Vieira SR. Spatial variability of soil physical properties on an experimental plot. Rev. Bras. Ciênc. Solo. 2005; 29(2):169-77.
• ³⁶ Blouin M, Hodson ME, Delgado EA, Baker G, Brussaard L, But KR, Brun J-J. A review of earthworm impact on soil function and ecosystem services. Eur. J. Soil Sci., 2013; 64(2):161-82. doi:10.1111/ejss.12025
  » https://doi.org/10.1111/ejss.12025
• ³⁷ Vieira ML, Klein VA, Propriedades físico-hídricas de um latossolo vermelho submetido a diferentes sistemas de manejo. Rev. Bras. Ciênc. Solo. 2007; 31(6): 1271-80.
• ³⁸ Zarehaghi D, Neyshabouri M. R, Gorji M, Hassanpour R, Bandehagh A. Growth and development of pistachio seedling root at different levels of soil moisture and compaction in greenhouse conditions. Soil and Water Research. 2017; 12(1): 60-6. doi:10.17221/146/2015-swr
  » https://doi.org/10.17221/146/2015-swr
• ³⁹ Daddow R, Warrington G. Growth-limiting soil bulk densities as influenced by soil texture. Watershed: System Development Group; 1983.
• ⁴⁰ Sadeghi M, Babaeian E, Tuller M, Jones SB The optical trapezoid model: A novel approach to remote sensing of soil moisture applied to Sentinel-2 and Landsat-8 observations. Remote Sensing of Environment. 2017; 198(1): 52-68. doi:10.1016/j.rse.2017.05.041
  » https://doi.org/10.1016/j.rse.2017.05.041

• Funding:
  This research was funded by National Council for Scientific and Technological Development (CNPQ).

Edited by

• Editor-in-Chief:
  Alexandre Rasi Aoki

• Associate Editor:
  Adriel Ferreira da Fonseca

Publication Dates

• Publication in this collection
  19 Nov 2021
• Date of issue
  2021

History

• Received
  03 Apr 2020
• Accepted
  10 June 2021