Soil Physical Properties Spatial Variability under Long-Term No-Tillage Corn

: Spatial variability of soil physical and hydrological properties within or among agricultural ﬁelds could be intrinsically induced due to geologic and pedologic soil forming factors, but some of the variability may be induced by anthropogenic activities such as tillage practices. No-tillage has been gaining ground as a successful conservation practice, and quantifying spatial variability of soil physical properties induced by no-tillage practices is a prerequisite for making appropriate site-speciﬁc agricultural management decisions and / or reformulating some management practices. In particular, there remains very limited information on the spatial variability of soil physical properties under long-term no-tillage corn and tropical soil conditions. Therefore, the main objective of this study was to quantify the spatial variability of some selected soil physical properties (soil surface temperature (ST), volumetric water content ( θ v ), soil resistance (TIP), total porosity ( θ t ), bulk density ( ρ b ), organic carbon, and saturated hydraulic conductivity ( K sat )) using classical and geostatistical methods. The study site was a 2 ha ﬁeld cropped no-tillage sweet corn for nearly 10 years on Oahu, Hawaii. The ﬁeld was divided into 10 × 10 and 20 × 20 m grids. Soil samples were collected at each grid for measuring ρ b , θ t , and soil organic carbon (SOC) in the laboratory following standard methods. Saturated hydraulic conductivity, TIP at 10 and 20 cm depths, soil surface temperature, and θ v were also measured. Porosity and ρ b have low and low to moderate variability, respectively based on the relative ranking of the magnitude of variability drawn from the coe ﬃ cient of variation. Variability of the SOC, TIP, and K sat ranges from moderate to high. Based on the best-ﬁtted semivariogram model for ﬁner grid data, 9.8 m and 142.2 m are the cut o ﬀ beyond which the measured parameter does not show any spatial correlation for SOC, and TIP at 10 cm depth, respectively. Bulk density shows the highest spatial dependence (range = 226.8 m) among all measured properties. Spatial distribution of the soil properties based on kriging shows a high level of variability even though the sampled ﬁeld is relatively small.


Introduction
Knowledge of soil physical, hydrological, chemical, and biological properties is important for many purposes including precision farming, environmental management, and crop growth modeling.Soils are characterized by a high degree of spatial variability because of the combined effects of physical, chemical, and biological processes that operate with different intensities and at different scales [1].Changes in soil properties and processes, in turn, impact plant growth and the environment [2] and require accurate knowledge of these properties for efficient utilization of agricultural inputs.Reduced tillage practices in agronomic crops such as corn, soybeans, cotton, sorghum, and cereal grains were introduced in the mid-twentieth century to conserve soil and water [3].In the United States, acreage under conservation tillage has increased from 2% in 1968 to 41% in 2004 [4,5].A recent study by Claassen et al. [6] showed adoption of different conservation tillage practices varies across crops; it revealed that conservation tillage was used on a majority of wheat (67%, 2017), corn (65%, 2016), and soybeans (70%, 2012) whereas conservation tillage was used on just 40% of cotton acres (2015) [6].Crops under conservation tillage are grown with minimum soil disturbance where crop residues are left in the field after harvesting of previous crop resulting in at least 30 percent of the soil is covered by crop residues [3].Conservation tillage practices include no-till/strip-till, ridge-till, and mulch-till [7].No-tillage, also known as no-till or zero tillage protects the soil surface and reduces soil erosion.However, no-tillage practice can lead to over-compaction of the soil [8].Several studies assessed the long-term effect of no-tillage on soil compaction; they reported that long-term practice of no-tillage increases soil organic carbon (SOC) and pore connectivity in the surface layer [8,9].Some studies [10,11] reported the formation of a relatively over compacted layer under no-tillage systems whereas other studies [12,13] did not experience any compacted layer formation under no-tillage conditions.There is evidence that no-tillage practices cause significant improvements in water retention, aggregate stability, temperature, and other properties compared to conventional tillage [14,15].Dao [16] found that although tillage initially decreases ρ b ; however, soil under no-till has lower ρ b than the soil under traditional tillage practices.In contrast, Hammel [17] found that soils under no-tillage had a higher ρ b as compared to conventional tillage in the top 30 cm soil surface.It seems there is site-specific effect of the response of some soil hydraulic properties to tillage treatments [18,19].
Spatial variability of soil properties may be enhanced by tillage and other on-farm management practices [20].Soil characteristics generally show spatial dependence [21], which means that samples close to each other are more similar than samples farther apart.Classic statistical methods have generally been used to study soil variability in the field; however, geostatistical methods have been used to understand the spatial correlation structure [22] of soil properties.Geostatistics take into account the spatial autocorrelation of a given variable in predicting its value in another unsampled location.Previous studies of spatial variability of soil physical properties within or among agricultural fields have yielded contrasting results across geographical locations, crops, soils, climates, and agricultural practices particularly tillage practices [18].Although no-tillage has become a successful conservation practice, quantifying spatial variability of soil physical properties induced by no-tillage practices is a prerequisite for making appropriate site-specific agricultural management decisions or reformulating some management practices.In particular, there remains very limited information on the spatial variability of soil physical properties under long-term no-tillage corn under tropical soil conditions.Therefore, the objective of this work was to quantify some of the selected soil physical properties (ST, θ v , θ t , TIP, ρ b , SOC, and K sat ) and describe their spatial dependence using classical and geostatistical methods.

Study Area
This study was conducted at the Waimanalo Research Station (21 • 20 ′ 15 ′′ N; 157 • 43 ′ 30 ′′ W) of the College of Tropical Agriculture and Human Resources, the University of Hawaii at Mānoa (Figure 1).The 2 ha experimental field (130 m × 160 m) was continuously cropped with sweet corn crop (Zea mays var.rugosa) under no-tillage for 10 years.The climate of the study area is tropical with a mean annual air temperature of 23 • C, and an average annual precipitation of 938 mm most of which occurs between November and April.The soil at this location has 2% to 6% slope and is classified as Waialua gravelly clay variant (isohyperthermic Pachic Haplustolls).The percentage of sand, silt, and clay at the top 18 cm of soil profile are 8.7, 45.2, and 46.1, respectively [23].Corn residues were not removed and management traffic patterns were parallel to the crop rows to minimize compaction in planting zone.

Soil Sample Collection and Laboratory Analysis
The 10 m × 10 m sampling grid was established on the study field (Figure 1).Undisturbed soil core samples (radius, 2.5 cm; height, 6 cm) were collected from the top 10 cm depth of each grid point between 6 December, 2008 and 9 December, 2008 using brass cylinders.These samples were used to determine ρ b and θ t , following the methods of Grossman and Reinsch [24] and Dane and Hopmans [25], respectively.
In the laboratory, soil core samples were placed in a tray filled with water for saturation after securing the soil from the bottom using nylon mesh.Samples were left in water for approximately 24 h to ensure their complete saturation.The saturated samples were weighed and then oven dried at 105 • C for at least 24 h.
After ρ b and θ t , analysis, the soil samples of the top 10 cm soil layers were used to determine the SOC using the loss on ignition (LOI) method [26] using a muffle furnace (Model 550 Isotemp Series, Fisher Scientific, Pittsburgh, PA, USA) and following the standard procedure of Schulte and Hopkins [27].Five grams of air-dry soil samples with a diameter <2 mm were placed in 15 mL ceramic cups, oven-dried at 105 • C for 24 h, cooled in a desiccator for a couple of hours, and then weighed (M 105 ).The samples were then combusted at 550 • C for 5 h [28].The combusted samples were again transferred to an oven at a temperature of 105 • C for a couple of hours.The samples were then taken out, cooled down in desiccators, and weighed again (M 550 ).The percent LOI was calculated as follows: A regression model was established between LOI and SOC using the Hawaii soil SSURGO data from the United States Department of Agriculture Natural Resources Conservation Services (USDA-NRCS).Established regression equation (Figure 2) was used to calculate SOC based on LOI values.

Field Measurement
The soil surface temperature, θ v , and soil resistance (TIP) were measured at each intersection point of 10 m × 10 m grid during 17 December, 2008 and 18 December, 2008.Similar measurements were also conducted on a 20 m × 20 m grid on 6 January, 2009 to quantify the temporal variability of ST, θ v , and TIP.The soil temperature was measured using the VWR ® two-channel thermometer (VWR International, LLC., Radnor, PA, USA).In situ θv was measured using the FIELDSCOUT ™ TDR 300 soil water content sensor, which has an accuracy ± 3.0% volumetric water content with electrical conductivity <2 dS m −1 of (Spectrum Technologies, Inc., Plainfield, IL, USA) in the top 10 cm depth.The soil penetration resistance was measured at 10 cm and 20 cm depths using DICKY-john ® soil compaction tester (DICKY-john ® Corporation, Auburn, IL, USA).All the above three measurements were conducted in three replicates within few cm and their mean values were used in the data analysis.
Infiltration measurements were conducted on 20 m × 20 m grid over the entire field and a 10 m × 10 m grid on the transect (160,0) using a tension infiltrometer (Soil Measurement Systems, Tucson, AZ, USA) between 6 December, 2008 and 10 December, 2008.The infiltration tests were carried out at two tensions, −15 and −20 cm H 2 O, with 5 cm tension differences.The falling head of water in the graduated water reservoir was monitored at 2 to 5 min intervals depending on the infiltration rates.The initial readings were taken frequently.At three consecutive intervals of around 5 min, if the rates did not change, steady state conditions were assumed, and the resulting data were used to calculate K sat using the following equation [29]: where Q (LT −1 ) is steady state water flux, r (m) is the radius of perforated disk, α (T −1 ) is inverse macro-porosity capillary length, ψ (L) is tension at which Q was measured.The inverse macro-porosity capillary length α was estimated using the relationship given by Hussen and Warrick [30] as follows: where ψ 1 and ψ 2 are the lower and upper bounds of the tension.

Statistical and Geostatistical Analysis
Mean, maximum, minimum, and variance were calculated for each of the properties measured during this study.Frequency distributions were also produced for each soil property and Shapiro-Wilk test [31] was used for normality test using Statistix software package (Analytical Software 2003, Tallahassee, FL, USA).
Spatial variability of ST, θ v , TIP, ρ b , θ t , SOC, and K sat was described using semivariograms assuming stationarity among similar lag that explains the average dissimilarity between the samples separated by a distance h.Semivariogram values at each lag separation γ(h) were calculated as follows: where N(h) is the number of data pairs separated by distance h, Z(x i ) is the measured value of the variable at location x i and Z(x i + h) is the measured value of the variable at a distance of h from the location x i .
Experimental semivariogram, cross-variogram, correlogram, and Moran's I were computed for each property using GS + [32] software package-version 7 (Gamma Design Software, LLC, Plainwell, Michigan, USA) and plotted with separation distance (h).During the variogram analysis, active lag distance was taken as half of the minimum extent of sampling area to avoid the bordering effect.In the estimation of the semi-variance, as the separation distance gets bigger, the number of samples separated by that distance gets smaller.In this study, all the geostatistical analysis was performed considering the isotropic condition and directional trend was ignored.The computed semivariogram and cross-variogram values were fitted with the theoretical semivariogram models.The best-fit model with the lowest value of residual sum of squares was selected for each soil property.Mathematically, the theoretical semivariogram models used in this study are as follows: Spherical: Exponential: Gaussian: Linear: where C 0 is the nugget variance and C is the structural variance.C 0 + C and A 0 are known as the sill and the theoretical range, respectively.In case of the spherical model, effective range is same as the theoretical range, whereas for the exponential and Gaussian models effective range is equal to 3A 0 and √ 3A 0 , respectively.Effective range is the distance at which the sill (C 0 + C) is within 5% of the asymptote.
Degree of spatial dependence (SD) for each measured property was calculated using the following equation: The spatial dependence of soil properties were classified according to the method detailed by Cambardella et al. [33], which represents the spatial randomness.The spatial dependence is strong if SD < 25%, moderate for SD between 26% and 75% and weak with SD > 75%.
Similar to Pearson's product moment correlation statistics, Moran's I [34] statistics for autocorrelation was used to quantify the spatial autocorrelation among the soil properties.In addition to Moran's I, correlogram was plotted for the analysis of spatial randomness in the measured data.Moran's I analysis has the index of autocorrelation among group of paired samples separated by increasing lag distance.Similar to correlogram, the values of Moran's I are in the range of −1 to 1. Positive values of Moran's I indicate that data values are spatially clustered, whereas negative values of Moran's I show that neighboring data values are dissimilar [35].The mathematical form of the index is given by: where I(h) is the Moran's I for distance interval h, Z i and Z i+h are the measured sample values at the point i and i + h, and N(h) is the number of samples separated by distance h.

Descriptive Statistics
Daily precipitation, minimum, and maximum temperatures during the sampling period are shown in Figure 3.The average value of K sat (0.91 cm min −1 or 1.52 × 10 −4 m s −1 ) is lower than the average values reported in similar studies of the long-term effects of no-tillage, for example, K sat = 3.5 × 10 −7 m s −1 for soils with 25% clay reported in Ferreras et al. [36], 4.2 × 10 −6 m s −1 for soils with 19.3% clay in Bhattacharyya et al. [37], and 3.54 × 10 −6 m s −1 for soils with 26% clay in Azooz and Arshad [38].Surface temperature has the lowest variability (CV = 1.87) whereas K sat has the highest (CV = 98.64) among the measured properties (Table 1).In previous works, Tsegaye and Hill [39], and Warrick and Nielsen [40] showed that K sat variability is normally high and CV values can reach up to 190%; whereas Mulla and McBratney [41] reported the values of CV for K sat in the range of 48%-352%.The average θ t value of 0.61 cm 3 cm −3 (Table 1) is relatively high as compared to the θ t values reported by Fares et al. [42] for Waialua gravelly clay variant amended.These high values of θ t can be attributed to the lower compaction of the field and the high plant residue and the root activity after harvesting.No-tillage also maintains large pores [43] inside soil aggregates as compared to conventional tillage which destroys soil aggregates as a result of soil tilth.Hill [44] showed greater soil strength for no-till soils compared to conventionally tilled soils at three different locations.Karlen et al. [45] reported an increase in soil aggregate stability, total carbon, microbial activity, and earthworm populations after eliminating tillage in the U.S. Corn Belt region over a 12-year period.
The average soil bulk density of 1.13 g cm 3 (Table 1) is slightly lower than that reported in the literature (1.20 g cm −3 ) for the top 38 cm soil depth of the Waialua gravelly clay variant This low bulk density value could be partially due to its shallow sampling depth as bulk density tends to increase with depth [17,20].Hill [44] reported statistically significant difference in bulk density under no-tillage at three locations when compared with conventional tillage.However, Karlen et al., [45] reported no significant difference in bulk density in top 50 cm depth under chisel plow and no-tillage treatments but there was significantly higher total carbon (almost twice as compared to plow and chisel) under no-tillage in the top 10 cm depth.As presented in Table 1, soil penetration resistance increases with depth and the resistance in the deeper depth (20 cm; TIP-20) showed higher spatial variability (CV = 52.54)as compared to the top 10 cm (TIP-10).Compaction in the lower depth tends to be more uniform due to fewer disturbances as compared to the top soil.
Coefficients of variation for bulk density and porosity were lower than the values of these soil physical properties reported by Mulla and McBratney [41].According to the relative ranking of magnitude of variability, porosity, and bulk density show low and low to moderate variability, respectively [41].However, variability of soil organic carbon, soil resistance, and K sat ranges from moderate to high.Frequency distributions of the measured properties show highly skewed and no-normal distribution for majority of soil parameters (Figure 4).Based on Shapiro-Wilk normality test, K sat data follows a log-normal distribution; this concurs with the findings of other studies, e.g., Vieira et al., [46]; Wierenga et al. [47]; Russo and Bouton [48].The volumetric water content and soil resistance at 10 cm depth data follow a normal distribution (Table 1).The volumetric water content and soil resistance data distribution varied as a function of the sampling date; this could be due to temporal effects.However, porosity, bulk density, soil surface temperature, and soil organic carbon seem to follow neither a normal nor a log-normal distribution.Soil properties collected based on the 20 m × 20 m grid on 6 January, 2009 were also skewed except for the volumetric water content and soil resistance at 20 cm (Table 1).Parkin and Robinson [49] reported that many soil properties follow skewed distributions are log-normally distributed.Data of soil properties, which did not exhibit normal distributions, according to Shapiro-Wilk test, were log transformed; however, only K sat data met the stationarity assumption and showed a normal distribution.Cambardella et al. [33] reported similar results for soil organic carbon in central Iowa soils; the log-normal transformation was only able to reduce skewness, data distribution remained kurtotic.

Variogram Analysis
Semivariogram parameters (range, nugget, and sill) for each soil property with the best-fitted model are shown in Figure 5, and their spatial dependencies are presented in Table 2. Values of semivariogram fitting parameters and models (Figure 5) varied with sampling date and grid size (Table 2).However, irrespective of sampling date and grid size, soil properties show moderate to strong spatial dependence based on the definition of Cambardella et al. [33].Nugget defines the micro-scale variability and measurement error, whereas partial sill indicates the amount of variation which can be defined by spatial correlation structure.Out of the total variation, the average nugget component is 35%; it ranges between 15% and 50%, which could be due to high micro-scale variation and/or measurement error.Nugget components for soil bulk density (48%), porosity (43%), and K sat (38%) were higher than the average nugget component (35%).Theoretically, the value of the range can be considered as the scale of distance beyond which the fitted parameters do not show any spatial correlation [46].For the 10 m × 10 m grid, value of the range varies between 9.8 m for soil organic carbon and 142.2 m for soil resistance at 10 cm depth.The volumetric water content measured at the 20 m × 20 m grid shows the highest spatial dependence (226.8 m) among all measured properties.Higher values of spatial dependence in volumetric water content can be attributed to the rainfall events occurred just before the field measurements (Figure 3).The goodness of fit measured as indicated by R 2 and the residual sum of squares (RSS) between measured and modeled semivariogram for the 10 m × 10 m sampling grid were slightly better than those for the 20 m × 20 m grid (Table 2).Among the different semivariogram models tested, the Gaussian model is found as the best fitting model.The spherical model best-fitted volumetric water content and soil resistance data for the 20 m × 20 m grid and volumetric water content and soil resistance combined samples collected on 17 December, 2008 and 18 December, 2008.The spatial correlation of soil organic carbon, porosity, and ln(K sat ) data followed exponential model.Cokriging is often recommended for spatial interpolation of measured data due to its smaller estimation variance compared to kriging [50] and its ability to accurately predict the property that is difficult or expensive to measure using well-behaved cross-semivariograms with another more easily sampled property [51].To assess the feasibility of cokriging, cross-semivariograms were established between measured properties.Before calculating the cross-semivariance, the degree of association between the data of every two studied properties was assessed using Pearson's correlation matrix (Table 3).There is a strong negative correlation (R = −0.72) between bulk density and porosity, and a strong positive correlation (R = 0.84) between soil resistance at 10 cm and 20 cm depth.The cross-semivariance between (i) porosity and bulk density, (ii) bulk density and soil organic carbon, (iii) K sat and bulk density, and (iv) K sat and soil resistance were negative.However, the cross-semivariance between soil resistance and bulk density, K sat, and porosity, K sat and soil organic carbon were all positive.For the volumetric water content, soil resistance, and soil temperature, the cross-semivariance was sometimes positive and sometimes negative depending on the date and/or the grid size of sampling (Table 4).

Analysis of Correlogram
Mulla [52] reported a negative cross-semivariance between gravimetric water content and surface temperature, and gravimetric water content and soil resistance with significant variation in range and partial sill.Variation in cross-semivariance within a small field is not surprising since both surface temperature and soil resistance depend on several properties including water content, soil albedo, solar radiation and soil bulk density, organic matter content, and caliche layer.The inconsistent correlation between soil resistance, soil temperature, and volumetric water content could be partly due to the wet condition of the field during the sampling time (Figure 3).In this study, soil temperature was found to be highly correlated with water content as compared to soil resistance and can be used as secondary variable for cokriging.The modeled variogram was poorly correlated with cross-semivariance between soil resistance and volumetric water content at 20 m × 20 m gird (Table 4) indicating its inappropriate scale for explaining the semivariance.
Moran's I correlogram analysis (Figure 6) for spatial dependence shows no autocorrelation at any lag distance for K sat .However, soil organic carbon was positively autocorrelated up to lag of 20 m.Porosity, volumetric water content, bulk density, soil resistance, and soil temperature show linear decreasing trend with point of inflection between 40 and 60 m.With lag distance higher than 60 m, the autocorrelation became negative except for porosity.The highest autocorrelation was associated with the 10 m (sampling grid size) lag distance.Most of these properties did not meet the assumption of stationarity, which induces some uncertainty in these results.However, these results can be used as primary information for designing future sampling plan in similar conditions.

Spatial Map from Kriging
Spatial maps of measured soil properties prepared using the ordinary kriging are presented in Figure 7.The spatial structure analysis based on kriging indicates a spatial variability across the field for volumetric water content, bulk density, porosity, soil organic carbon, soil resistance, and K sat .Soil temperature did not show a significant spatial variability as it is expected due to small amplitude of diurnal air temperature at the study site and the Hawaiian Islands in general.The volumetric water content is higher in the downslope areas than in the upslope areas.Values of the soil bulk density and porosity were lower and higher, respectively in the lower portion of the study site.This shows the presence of restricting subsurface layer in the northern part of the field that could be due to the extensive continuous corn cultivation in these areas as compared to lower part of the field.We did not observe significant difference in soil resistance at surface (0-10 cm) and subsurface (10-20 cm) layers.

Summary and Conclusions
Soil physical and hydrological properties were measured on grid basis from an agricultural field cropped with sweet corn under a long term no-tillage practice.Classic statistics and geostatistics were used to assess the spatial dependence in several soil physical properties (soil surface temperature, volumetric water content, soil resistance, porosity, bulk density, soil organic carbon, and saturated hydraulic conductivity).Surface temperature has the lowest spatial variability (CV = 1.87) whereas K sat has the highest spatial variability (CV = 98.64).Values of soil porosity and bulk density spatial variability indicators were relatively higher and lower, respectively, as compared to those of similar soils as reported in the literature.Soil resistance increased with depth indicating compaction of the subsurface layer.The degree of spatial dependence varied from strong to moderate for all tested properties.The variability of volumetric water content, soil temperature, soil resistance, bulk density, porosity, soil organic carbon, and K sat exhibited spatial dependence of varying range that can be described using different semivariogram models.The range varied from 9.8 m for soil organic carbon to 142.2 m for soil resistance at 10 cm depth.However, volumetric water content measured at 20 m × 20 m grid had the highest range value (226.8 m).Spatial distribution of the soil properties based on the results of the kriging analysis has shown high level of spatial variability across a relatively small area.These observations suggest that soil physical properties are spatially variable and that spatial variability should not be ignored.Thus, assuming uniform soil properties can lead to mismanagement of natural resources.Furthermore, results of this work revealed that since the scale of spatial variability is not constant across the different soil hydrological and physical properties, it is crucial to use the appropriate scale for the specific practical application.For instance, the degree of the spatial dependency for the volumetric water content is about 227 m.Based on this, it is practical to design field irrigation units for this tropical soil at that level where it is safe to apply constant irrigation amounts at that scale.However, for other agricultural practices such as fertilizer's applications, it might be prudent to use a finer scale, e.g., 10 m, that corresponds to extend of the maximum of spatial dependency of soil organic carbon content which is an important factor that should be considered in any crop fertility program.Further work is needed to fine-tune some of these practical implications of soil hydrological and physical properties spatial variabilities before recommending them for farmers of the region.

Figure 1 .
Figure 1.A diagram and map of the study area at the University of Hawaii at Mānoa Waimanalo Research Station; it shows sampling grids (Empty circles-10 m × 10 m sampling grids, circle with red dots-20 m × 20 m sampling grids).

Figure 2 .
Figure 2. Relationship between loss on ignition (LOI) and soil organic carbon for the Waialua gravelly clay variant soil in Waimanalo, HI.

Figure 4 .
Figure 4. Frequency distribution of selected soil properties sampled at the 10 m × 10 m grid on 17-18 December 2008.

Table 1 .
Descriptive Figure 3. Measured values of daily precipitation, and minimum (T min ) and maximum (T max ) air temperatures at the study site between 1 December 2008 and 31 January 2009.statistics and Shapiro-Wilk normality test results (grouped by sampling date and grid size) for measured soil physical properties: soil temperature (ST), volumetric water content (θv), soil resistance at 10 cm depth (TIP-10), soil resistance at 20 cm depth (TIP-20), bulk density (ρ b ), porosity (θ t ), soil organic carbon (SOC), and saturated hydraulic conductivity (Ksat).* Log transformed data, Var = variance, CV = coefficient of variation, Max = maximum, Min = minimum, N = number of samples.
Significant at the 0.05 level, ** Significant at the 0.01 level. *

Table 4 .
Summary of cross-variogram model.