Spatial Variability of Soil Physical and Hydraulic Properties in a Durum Wheat Field: An Assessment by the BEST-Procedure

Spatial variability of soil properties at the field scale can determine the extent of agricultural yields and specific research in this area is needed. The general objective of this study was to investigate the relationships between soil physical and hydraulic properties and wheat yield at the field scale and test the BEST-procedure for the spatialization of soil hydraulic properties. A simplified version of the BEST-procedure, to estimate some capacitive indicators from the soil water retention curve (air capacity, ACe, relative field capacity, RFCe, plant available water capacity, PAWCe), was applied and coupled to estimates of structure stability index (SSI), determinations of soil texture and measurements of bulk density (BD), soil organic carbon (TOC) and saturated hydraulic conductivity (Ks). Variables under study were spatialized to investigate correlations with observed medium-high levels of wheat yields. Soil physical quality assessment and correlations analysis highlighted some inconsistencies (i.e., a negative correlation between PAWCe and crop yield), and only five variables (i.e., clay + silt fraction, BD, TOC, SSI and PAWCe) were spatially structured. Therefore, for the soil–crop system studied, application of the simplified BEST-procedure did not return completely reliable results. Results highlighted that (i) BD was the only variable selected by stepwise analysis as a function of crop yield, (ii) BD showed a spatial distribution in agreement with that detected for crop yield, and (iii) the cross-correlation analysis showed a significant positive relationship between BD and wheat yield up to a distance of approximately 25 m. Such results have implications for Mediterranean agro-environments management. In any case, the reliability of simplified measurement methods for estimating soil hydraulic properties needs to be further verified by adopting denser measurements grids in order to better capture the soil spatial variability. In addition, the temporal stability of observed spatial relationships, i.e., between BD or soil texture and crop yields, needs to be investigated along a larger time interval in order to properly use this information for improving agronomic management.


Introduction
Soil physical and hydraulic properties can change drastically over space [1] and time [2] and their evaluation is essential for a rational soil management and, therefore, for increasing crop yields performance [3]. Moreover, they dynamically affect water balance components and crop yields by relating soil hydraulic functioning to climate patterns and crop water requirements [4,5].
Soil properties, such as saturated hydraulic conductivity, total organic carbon content, structure stability index, dry bulk density, as well as capacitive indicators obtained from the soil water retention curve (i.e., plant available water capacity, relative field capacity, air capacity) were widely and successfully applied to investigate soil management effects on soil physical and hydraulic properties [3,[6][7][8]. Keller et al. [3], for example, have investigated the relationships between crop yield and soil structure in three Swedish fields, applying the simplified falling head (SFH) technique [9] to evaluate whether field-saturated hydraulic conductivity (Ks) could be used as a simple and quickly measurable indicator of crop yield. The main findings of their investigation showed that Ks may be a good indicator of low yielding zones, and degradation of soil structure has been indicated as the main reason for low yield. However, in other studies, spatial and temporal variation of soil water holding capacity was suggested to be a factor partly responsible for crop yield variation, regardless of the amount and timing of the rain contributions [10,11]. Consequently, investigations addressed at evaluating new experimental procedures for soil hydraulic characterization, and at establishing the usefulness and sensitivity of soil properties as predictors for crop yield, are needed. This is a current issue in Mediterranean agro-environments, where water resources are scarce and need to be optimally managed [12].
Soil hydraulic properties assessment, i.e., water retention curve and hydraulic conductivity function (θ(h) and K(h) relationships, respectively), is expensive and time consuming, since standard methods need specific skills and their spatialization, i.e., prediction on a distributed large scale, which can be burdensome or wholly inapplicable [13].
Several quick methods are available to obtain θ(h) and K(h) (or K(θ)) relationships. Pedotransfer functions (PTFs), for example, allow for estimating soil hydraulic properties starting with basic soil variables such as soil texture, bulk density or organic matter or hydraulic conductivity [14]. However, PTFs are not able to accurately quantify the effects of different agronomic options for soil management, unless a site-specific calibration is performed [15]. On the other hand, the most widely applied laboratory methods as hanging water column apparatus [16], pressure cells [17] or evaporation method [18], may require up to a week (or more) to obtain a single soil water retention curve [19,20]. As a consequence, although reliable and accurate, they are not easily applicable for large-scale research and simplified techniques should be chosen for these purposes [21,22].
The BEST-procedure by Lassabatère et al. [23] allows for estimating hydraulic functions repeatedly over space and time with substantially limited experimental burdens. Basically, the procedure requires three sets of experimental information: (i) cumulative infiltration by a simple infiltrometric experiment (i.e., ring infiltration test of Beerkan type), (ii) soil bulk density and volumetric soil water content at the time of experiment by sampling few (generally two) undisturbed soil cores, and (iii) soil particle-size distribution or, alternatively, clay, silt and sand percentages according to the USDA classification. Specifically, the procedure makes use of well-known analytical solutions for θ(h) and K(h) relationships and estimates (i) shape parameters, which are texture dependent, from particle-size analysis by physical-empirical PTFs, while (ii) scale parameters, which are structure dependent, by a three-dimensional field infiltration experiment at zero pressure head [23]. Therefore, BEST can be considered an adequate compromise between estimation accuracy and economic-experimental load. For example, some studies have applied BEST to establish the effects of droplet impact on soil sealing and crust formation [24,25], to carry out integrated soil physical quality assessment [26,27] or to identify the effects of tillage on some soil properties (i.e., Ks) under drip irrigation [28]. A main advantage of BEST is that it can be adopted when a large number of hydraulic measurements must be obtained at the field or at irrigation district scale, for example, for precision agriculture purposes. For instance, Mubarak et al. [29] assessed the temporal stability of both Ks and spatial structure of hydraulic properties of a loamy soil. Specifically, they compared Ks-BEST data with those estimated seventeen years earlier by applying the guelph permeameter method, under relatively comparable soil and agronomic conditions. Results showed that Ks changed significantly, but observed discrepancies were not higher than a factor of three or four. This suggests that BEST can represent an easy, robust, and inexpensive way for characterizing soil hydraulic behavior and its spatial [30] and temporal [29] variability at the field scale. The availability of a large number of georeferenced hydraulic measurements would allow delineating homogeneous sub-areas within the crop field to be submitted to uniform agronomic management [31,32]. This can lead to an increase in agricultural resources optimization [31]. In addition, dense spatial data can be used as auxiliary/covariate information in mixed effects models to improve the estimates of the target variables [33] or to improve the estimation of treatment significance reducing the risk of misleading or erroneous inferences in analysis of variance [34,35]. In other words, the potential application advantages seem attractive, but BEST has not yet been tested for soil hydraulic properties spatialization and the actual reliability for the mentioned purposes must be proven.
The general objective of this study was to test the BEST-procedure for the spatialization of soil hydraulic properties. In particular, the spatial distribution of the measured-estimated by BEST variables (soil texture, bulk density, saturated hydraulic conductivity, plant available water capacity, relative field capacity, air capacity) and ancillary soil properties (total organic carbon content, structure stability index) was investigated in a typical wheat cropping system in Southern Italy. Cross-correlation analysis was applied to establish strength and the extent of the spatial relationships between selected soil properties and crop yields.

Study Area
The research was carried out in the spring-summer period of 2016 at the experimental farm "Menichella" of Council for Agricultural Research and Economics, located near Foggia (41 • 27 02 N, 15 • 30 36 E), Southern Italy ( Figure 1). The study was conducted in a field of about 5 ha (170 m × 250 m) conventionally cultivated with durum wheat (Figure 1). The field is located in a flat area characterized by Mediterranean climatic conditions and the soil was clay according to USDA classification. Additional information on the experimental site can be found in Cavallo et al. [31]. measurements would allow delineating homogeneous sub-areas within the crop field to be submitted to uniform agronomic management [31,32]. This can lead to an increase in agricultural resources optimization [31]. In addition, dense spatial data can be used as auxiliary/covariate information in mixed effects models to improve the estimates of the target variables [33] or to improve the estimation of treatment significance reducing the risk of misleading or erroneous inferences in analysis of variance [34,35]. In other words, the potential application advantages seem attractive, but BEST has not yet been tested for soil hydraulic properties spatialization and the actual reliability for the mentioned purposes must be proven.
The general objective of this study was to test the BEST-procedure for the spatialization of soil hydraulic properties. In particular, the spatial distribution of the measured-estimated by BEST variables (soil texture, bulk density, saturated hydraulic conductivity, plant available water capacity, relative field capacity, air capacity) and ancillary soil properties (total organic carbon content, structure stability index) was investigated in a typical wheat cropping system in Southern Italy. Cross-correlation analysis was applied to establish strength and the extent of the spatial relationships between selected soil properties and crop yields.

Study Area
The research was carried out in the spring-summer period of 2016 at the experimental farm "Menichella" of Council for Agricultural Research and Economics, located near Foggia (41°27′02" N, 15°30′36" E), Southern Italy ( Figure 1). The study was conducted in a field of about 5 ha (170 m × 250 m) conventionally cultivated with durum wheat (Figure 1). The field is located in a flat area characterized by Mediterranean climatic conditions and the soil was clay according to USDA classification. Additional information on the experimental site can be found in Cavallo et al. [31].
Fifty-two geo-referenced sampling points, located at the nodes of a regular grid (175 m × 250 m) with a mesh of 20 m × 40 m, were considered in this investigation ( Figure 1). The amount of observations was checked to ensure a sufficient coverage of the study area, i.e., a number of pairs for each distance class larger than the minimum required for an effective spatial analysis, as reported by Myers [36] (minimum threshold equal to 25 pairs; actual pairs observed ranging between 66 and 239). At each selected point, the wheat yield and some physical, hydraulic and chemical soil properties (e.g., soil texture, dry bulk density, saturated hydraulic conductivity, total organic carbon content) were determined to establish possible spatial structures and correlations among variables. More details on the dataset composition are given in Section 2.2. Fifty-two geo-referenced sampling points, located at the nodes of a regular grid (175 m × 250 m) with a mesh of 20 m × 40 m, were considered in this investigation ( Figure 1). The amount of observations was checked to ensure a sufficient coverage of the study area, i.e., a number of pairs for each distance class larger than the minimum required for an effective spatial analysis, as reported by Myers [36] (minimum threshold equal to 25 pairs; actual pairs observed ranging between 66 and 239). At each selected point, the wheat yield and some physical, hydraulic and chemical soil properties (e.g., soil texture, dry bulk density, saturated hydraulic conductivity, total organic carbon content) were determined to establish possible spatial structures and correlations among variables. More details on the dataset composition are given in Section 2.2.

Lab and Field Measurements
For each of the selected fifty-two sampling points, a simplified version of the BEST-procedure by Lassabatère et al. [23] was applied. The saturated hydraulic conductivity (Ks) and soil bulk density (BD) were then measured, and three capacitive indicators were estimated from soil water retention curve of BEST, namely plant available water capacity (PAWCe), air capacity (ACe) and relative field capacity (RFCe) [8]. The subscript e is used for indicating variables estimated by BEST. In detail, as required for BEST application, a falling head infiltration experiment of Beerkan type was carried out in the third decade of April at each sampling point using a metal ring (15 cm inner diameter), and the cumulative infiltration curve as function of the time (I(t)) was obtained; eighteen water volumes of 200 mL each were used in this investigation and the BEST-steady algorithm by Bagarello et al. [37] was applied to estimate the aforementioned soil properties. A relatively higher number of water volumes than usual (eighteen instead of fifteen) was chosen to be sure of sampling a representative soil volume and be more confident to reach steady-state conditions of water flow. Two undisturbed soil cores (0.05 m in height by 0.05 m in diameter) were collected at the 0 to 0.05 m and 0.05 to 0.10 m depths close to the ring to determine the soil water content at the beginning of infiltration experiments, θ i , and BD. A disturbed soil sample (0-0.10 m depth) was collected at each sampling point to quantify both the soil texture, i.e., clay (Cl), silt (Si) and sand (Sa) percentages, according to the USDA classification, and the soil total organic carbon content (TOC). Soil texture was obtained with the standard pipette method [38], and TOC was quantified through dry combustion using a TOC Vario Select analyzer [39]. Soil texture fractions, BD, θ i and I(t) were used to run BEST-steady and the infiltration constants, β and γ, were fixed at the reference values of literature [23,40]. More detail concerning the BEST-procedure application can be found in Castellini et al. [8]. At harvesting (end of June 2016), wheat grain yield data were collected on the fifty-two geo-referenced locations (i.e., nodes of a 20 m × 40 m grid). Yield data were recorded on sampling areas of 1 m 2 and normalized to 13% moisture content of grain.

Physical and Hydraulic Soil Properties
The saturated hydraulic conductivity (Ks) is the soil's ability to absorb and transmit soil water to the root zone, as well as drain excess water out of the root zone [41]. Since Ks is mainly controlled by soil structure and texture, e.g., [42], in the same soil or soil class it may be used as a measure of structural status of agricultural soils [3]. References of literature suggest optimal Ks values for agriculture soils within the range 0.005-0.05 mm s −1 to promote a rapid infiltration and redistribution of plant available water [41]. However, in order to select soil properties that directly (or indirectly) account for soil structure, we also considered (i) dry bulk density (BD), (ii) total organic carbon content (TOC) and (iii) structural stability index (SSI) by Pieri [43]. SSI is calculated from TOC and fine soil texture (SSI = 1.724·TOC%/(silt% + clay%)·100) components [27]. For these soil indicators the following critical limits were considered: optimal BD values within the range 0.9-1.2 g cm −3 ; optimal or poor TOC values, general for plant husbandry, were equal to 30-50 or <23 g kg −1 , respectively, although good values for agricultural fine-textured soils were reported in the range of 15-22 g kg −1 by Sequi and Nobili [44]. Therefore, a comparison between these two references has been made. SSI values >7% or ≤7% are representative of low or high risk of structural degradation, respectively [6].
Conversely, a second set of indicators that gives an account of the proportion between water and air into the soil was considered. Plant available water capacity (PAWCe) (cm 3 cm −3 ), in fact, is the amount of water held in the soil and available for crop growth and obtained as the difference between the water contents at field capacity (at h = −100 cm) and at permanent wilting point (at h = −15,300 cm) [45]; according to the literature [6], the following PAWCe limits were considered in this investigation: PAWCe ≥ 0.20 ideal; 0.15 < PAWCe < 0.20 good; 0.10 < PAWCe < 0.15 limited; PAWCe < 0.10 poor. Air capacity (ACe) (cm 3 cm −3 ) provided information on the soil ability to store and transmit air (Reynolds et al. [6]). According to Castellini et al. [7], an optimal ACe value falls in the range 0.10-0.26 cm 3 cm −3 , while higher or lower values represent inadequate soil aeration conditions; this interval was optimal for a clay soil bordering with that studied [7]. Finally, the relative field capacity (RFCe), obtained as the ratio between the water contents at field capacity and at water saturation, was considered, since it partially combines ACe and PAWCe, thus expressing the soil capacity to store air and water relative to the soil's total pore volume [46]. Optimal values for RFCe were suggested within the range 0.6-0.7; consequently, lower or higher values are representative respectively of "water limited" or "aeration limited" (i.e., for a given soil texture, too porous or too compact) soil conditions [6].

Preliminary Statistical Analysis
Descriptive statistics were computed in order to summarize the main features of data distribution for yield and the soil variables under study: bulk density (BD), capacitive indicators from the estimated soil water retention curve (ACe, PAWCe, RFCe), saturated soil hydraulic conductivity (Ks), soil total organic carbon content (TOC), fine soil texture components (Cl + Si) and soil structural stability index (SSI). In addition, hypothesis of normality was tested using the Kolmogorov-Smirnov test [47].

Correlation and Spatial Analysis
Relationships among studied variables were investigated using parametric correlation analysis, computing Pearson correlation coefficients.
The predisposition of the considered variables to be spatialized was investigated using Moran statistics. Moran's autocorrelation coefficient (often denoted as I) is an extension of the Pearson product-moment correlation coefficient to a univariate series [48]. Specifically, Moran statistics computes a weighted Pearson product-moment correlation of a variable against itself, where the weighting relates to the variable's spatial arrangement [49]. Moran's I allows for the investigation of correlation within a single variable due to the spatial relationship amongst its observations. The weights (w ij ) are a function of the distance between each pair of observations of the variable under study (x i ; x j ). In its simplest form, weights will take values 1 for close neighbours, and 0 otherwise; we also set w ii = 0. These weights are sometimes referred to as a neighbouring function.
Moran's I statistics is represented by the following equations: where z i is the deviation of an attribute from its mean (x i − x) and S 0 is the aggregate of all the spatial weights (w ij ).

Geostatistical Analysis
The geostatistical analysis is aimed at evaluating spatial heterogeneity of the studied variables by means of structural analysis (variography) and at producing maps of the collected data by means of spatial prediction approaches (kriging).

Variography
Under the name of (semi)variogram, two different types of functions related to geographically referenced data are usually denoted, namely the experimental variogram and the variogram model. The former is discrete and represents the half of the average squared difference between points separated by the distance h (the red-dotted line in Figure 2). The latter is parametric, continuous and a conditionally negative definite [50] (the solid line in Figure 2). The most important parameters for the model are the partial sill (σ 2 ) indicating the structured component of the variance, the nugget (σ 2 0 ), indicating the random uncorrelated component, and the range (α); this latter can be interpreted as the distance beyond which the spatial correlation becomes negligible ( Figure 2). A list of the main variogram models is summarized in Figure 2.
Water 2019, 11, x FOR PEER REVIEW 6 of 19 separated by the distance h (the red-dotted line in Figure 2). The latter is parametric, continuous and a conditionally negative definite [50] (the solid line in Figure 2). The most important parameters for the model are the partial sill (σ 2 ) indicating the structured component of the variance, the nugget ( ), indicating the random uncorrelated component, and the range (α); this latter can be interpreted as the distance beyond which the spatial correlation becomes negligible ( Figure 2). A list of the main variogram models is summarized in Figure 2. Variogram models were fitted to the experimental variograms of the variables under study. By virtue of the parsimony principle, the preferred spatial model should be isotropic. Regardless, a test has been applied to test for anisotropy. In detail, isotropic and anisotropic models were compared by means of a Likelihood Based Parameter Estimation for Gaussian Random Fields test, using REML, in order to estimate the two parameters characterizing anisotropy, namely the anisotropy angle and the ratio between the two ellipse axes [51].
Leave-one-out cross-validation was carried out and Pearson correlation coefficients (r) between predicted and observed data were used to quantify the goodness of model adaptation to the experimental variogram.

Spatial Interpolation
All the variables were interpolated using a univariate approach, the ordinary kriging (OK). Such a predictor is one of the most basic forms of kriging in which the unknown value z(x0) of a given realization of Z(x0) is predicted from the known values z(xi) i = 1, 2, ..., N, at the support points xi. The ordinary kriging predictor can be written as: where λi are weights associated with the N sampling points. The weights are chosen in such a way that the predictor is unbiased, the values are continuous and the estimation error is minimized: This ensures that kriging is an exact interpolator because the estimated values are identical to the observed values when a kriged location coincides with a sample location [33,52,53].

Model
Variogram Parameter Restrictions Variogram models were fitted to the experimental variograms of the variables under study. By virtue of the parsimony principle, the preferred spatial model should be isotropic. Regardless, a test has been applied to test for anisotropy. In detail, isotropic and anisotropic models were compared by means of a Likelihood Based Parameter Estimation for Gaussian Random Fields test, using REML, in order to estimate the two parameters characterizing anisotropy, namely the anisotropy angle and the ratio between the two ellipse axes [51].
Leave-one-out cross-validation was carried out and Pearson correlation coefficients (r) between predicted and observed data were used to quantify the goodness of model adaptation to the experimental variogram.

Spatial Interpolation
All the variables were interpolated using a univariate approach, the ordinary kriging (OK). Such a predictor is one of the most basic forms of kriging in which the unknown value z(x 0 ) of a given realization of Z(x 0 ) is predicted from the known values z(x i ) i = 1, 2, . . . , N, at the support points x i . The ordinary kriging predictor can be written as: where λ i are weights associated with the N sampling points. The weights are chosen in such a way that the predictor is unbiased, the values are continuous and the estimation error is minimized: This ensures that kriging is an exact interpolator because the estimated values are identical to the observed values when a kriged location coincides with a sample location [33,52,53].

Cross-Correlogram
In order to compare prediction maps of different variables, cross-correlograms were computed. In detail, given the prediction maps of two different variables, M A and M B , there are several methods with which to compare them [33,54,55]. To account for the inherently spatial characteristics of map representation, the cross-correlogram, which measures the correlation as a function of the distance between observations, is particularly well-suited. The analytical formulation of the cross-correlogram is the following: where z i,jA and z i ,j B represent the values at locations (i, j) and (i , j ) of the two maps separated by the h distance, respectively; h = (i − i ) + ( j − j ) represents the distance between the two locations, E denotes the mathematical expectation, m A and m B represent the populations means and s A and s B represent the populations standard deviations. If patterns are completely similar, apart from a constant, ρ A,B (0) should be equal to 1. To estimate ρ A,B (h) from the available data the following equation can be used: To compute r A,B (h), the procedure is as follows: from both the maps, all the couples whose locations are separated by the distance h are collected. Indicesm A andm B andŝ A andŝ B represent the mean and the standard deviation of mapped z i,jA and z i ,j B , respectively. N(h) is the total number of these pairs. If the methods being compared produce similar results, a decreasing cross-correlogram for increasing values of h is expected.
For summarizing the results and comparing different outcomes, results of cross-correlation analysis were reported in form of tables containing correlations at specified lag distances (0 m, 25 m, 50 m, 75 m, 100 m); in this way, practical information on the strength and extent of the spatial relationships between soil properties and crop yields was provided.

Preliminary Statistical Analysis
The preliminary statistical analysis on the studied variables highlighted some slight departures from normal distribution (Table 1). In particular, Ks and TOC showed positive and negative longer tails due to a few larger and smaller observations. The Kolmogorov-Smirnov test results indicated a significant departure from normality only for ACe and Ks, although not extremely significant (P = 0.0291; P = 0.0209; Table 2). For this reason, it was not deemed necessary to transform the original data.

Physical and Hydraulic Soil Properties
Physical and hydraulic properties of the investigated clay soil were summarized in Figure 3. Overall, according to the suggested criteria to detect relatively good or poor soil conditions of BD, TOC, SSI, Ks, PAWCe, ACe, RFCe and to obtain acceptable crop yields, results showed relatively good findings as only three of the seven soil properties indicated non optimal conditions. In further detail, the following was observed: (i) relatively low levels of organic carbon content (TOC) for threshold values reported by Reynolds [6], but good levels for ranges defined by Sequi and Nobili [44], (ii) risks of structural instability (SSI) and (iii) potential conditions of soil compaction (RFCe). In agreement with these findings, although ACe was within the suggested optimal range, it showed relatively low air capacity (Figure 3). Optimal hydrodynamic soil properties were also recognized for Ks, as the entire range of variation of the measurements fell within the limits defined by Reynolds et al. [41]. However, two clarifications should be provided because: (i) observed low TOC values are not entirely unexpected as they are quite common for a monoculture of wheat, in a typical Mediterranean environment, and (ii) BD evaluation showed optimal conditions, i.e., a soil compaction was not recognized, indicating a different assessment when compared with RFCe or ACe (Figure 3). In other words, soil physical and hydraulic assessment carried out using available references of literature has suggested optimal, or near optimal, conditions of investigated soil. This finding was in agreement with results of wheat yields, as obtained mean values equal to 3.86 t ha −1 may be considered as medium-high production levels for the investigated agro-environment. However, as TOC limits used do not appear to be entirely in line with the case under study, a specific discussion has been made in Section 4. Water 2019, 11, x FOR PEER REVIEW 9 of 19 Figure 3. Box plots of dry bulk density (BD), total organic carbon content (TOC), structure stability index (SSI), saturated hydraulic conductivity (Ks), plant available water content (PAWC), air capacity (AC) and relative field capacity (RFC). The thick-black line within each box represents the mean value (the fine-black line, the median). Open circles represent outliers (closed circle, extreme outlier for TOC). Dot green-line or dashed red-line represents optimal or critical values, respectively, according to section 2.3. For TOC, a further classification by Sequi and Nobili [43] was specifically reported (solid lines, which identify four areas, from poor to very good) for fine textured soils (clay, clay-loam, silty-clay and silty-clay-loam, according to USDA).

Correlation and Spatial Analysis
Results of correlation analysis are reported in Figure 4. Significant correlations were observed between crop yield and BD and PAWCe, with a positive and negative relationship respectively (r = 0.381, P = 0.005 and r = −0.400, P = 0.003); in addition, a positive correlation, although not significant (P = 0.063), was found with Ks. Overall, soil variables showed interesting relationships. As concerns capacitive indicators derived from the soil water retention curve of BEST, a strong negative correlation was found between RFCe and ACe (r = −0.922, P < 0.0001) as observed also in previous studies [7]. PAWCe was strongly negatively related to BD (r = −0.780, P < 0.0001). An interesting relationship was observed between Ks and the three capacitive indicators (P = < 0.0001, 0.001 and 0.009, respectively for RFCe, ACe and PAWCe). TOC was negatively correlated to BD (P < 0.006) and showed a weak relation with PAWCe (P < 0.065). Finally, SSI, synthetizing the information deriving from TOC and fine texture components (Cl + Si), showed to be an effective summary indicator highlighting good correlations also with BD (P = 0.005), PAWCe (P = 0.006) and Ks (P = 0.059). A more in-depth discussion on the correlations obtained is reported in the Discussion section.  Figure 3. Box plots of dry bulk density (BD), total organic carbon content (TOC), structure stability index (SSI), saturated hydraulic conductivity (Ks), plant available water content (PAWC), air capacity (AC) and relative field capacity (RFC). The thick-black line within each box represents the mean value (the fine-black line, the median). Open circles represent outliers (closed circle, extreme outlier for TOC). Dot green-line or dashed red-line represents optimal or critical values, respectively, according to Section 2.3. For TOC, a further classification by Sequi and Nobili [43] was specifically reported (solid lines, which identify four areas, from poor to very good) for fine textured soils (clay, clay-loam, silty-clay and silty-clay-loam, according to USDA).

Correlation and Spatial Analysis
Results of correlation analysis are reported in Figure 4. Significant correlations were observed between crop yield and BD and PAWCe, with a positive and negative relationship respectively (r = 0.381, P = 0.005 and r = −0.400, P = 0.003); in addition, a positive correlation, although not significant (P = 0.063), was found with Ks. Overall, soil variables showed interesting relationships. As concerns capacitive indicators derived from the soil water retention curve of BEST, a strong negative correlation was found between RFCe and ACe (r = −0.922, P < 0.0001) as observed also in previous studies [7]. PAWCe was strongly negatively related to BD (r = −0.780, P < 0.0001). An interesting relationship was observed between Ks and the three capacitive indicators (P = < 0.0001, 0.001 and 0.009, respectively for RFCe, ACe and PAWCe). TOC was negatively correlated to BD (P < 0.006) and showed a weak relation with PAWCe (P < 0.065). Finally, SSI, synthetizing the information deriving from TOC and fine texture components (Cl + Si), showed to be an effective summary indicator highlighting good correlations also with BD (P = 0.005), PAWCe (P = 0.006) and Ks (P = 0.059). A more in-depth discussion on the correlations obtained is reported in the Discussion section. Ks is the saturated soil hydraulic conductivity (mm s −1 ), TOC is the soil total organic carbon (g kg −1 ), SSI is the structure stability index (%), PAWCe is the plant available water capacity (cm 3 cm −3 ), ACe is the air capacity (cm 3 cm −3 ) and RFCe is the relative field capacity (dimensionless).
Moran's I spatial autocorrelation statistics gave a hint about the predisposition of the single variables to be spatialized (Table 3). For the case at hand, Moran values ranged from −0.2232 for RFCe, indicating a non-significant spatial correlation, to 0.3920, 0.4474, 0.4564 and 0.5889, for BD, SSI, PAWCe and Cl + Si respectively, indicating highly significant overall spatial dependence. Yield and TOC showed also a significant spatial correlation (Table 3). It is noteworthy that the variables correlated each other, showing also a similar overall spatial structure. Theoretical models were fitted to the experimental variograms of the variables under study, consisting generally of these nested models: a nugget effect and a spatial covariance function. Since the anisotropy ratio and the anisotropy angle parameters were not significant according to the Likelihood Based Parameter Estimation test, isotropic models were selected. In Table 4, the fitted variogram models are reported and described by the model name and the following parameters, namely nugget, partial sill and range. Actually, the Matérn model demonstrated to be the best suited theoretical function to describe spatial structure of the considered variables except for BD and PAWCe, since it was more flexible with respect to the classical models, having an additional shape parameter (k) (Figure 2). Following Cambardella et al. [56], the nugget-to-sill ratio can be a useful means to describe the spatial structure of the studied variables. In particular, the nugget semivariance expressed as a percentage of the total semivariance enables comparison of the relative size of the nugget effect among studied variables, with ratios < 25% indicating strong spatial dependence, ratios between 25 and 75% moderate spatial dependence, ratios > 75% weak spatial Red dots highlight significant correlations at P < 0.05. BD is the bulk density (g cm −3 ), Cl + Si is the sum of clay and silt content (%), Ks is the saturated soil hydraulic conductivity (mm s −1 ), TOC is the soil total organic carbon (g kg −1 ), SSI is the structure stability index (%), PAWCe is the plant available water capacity (cm 3 cm −3 ), ACe is the air capacity (cm 3 cm −3 ) and RFCe is the relative field capacity (dimensionless).
Moran's I spatial autocorrelation statistics gave a hint about the predisposition of the single variables to be spatialized (Table 3). For the case at hand, Moran values ranged from −0.2232 for RFCe, indicating a non-significant spatial correlation, to 0.3920, 0.4474, 0.4564 and 0.5889, for BD, SSI, PAWCe and Cl + Si respectively, indicating highly significant overall spatial dependence. Yield and TOC showed also a significant spatial correlation (Table 3). It is noteworthy that the variables correlated each other, showing also a similar overall spatial structure. Theoretical models were fitted to the experimental variograms of the variables under study, consisting generally of these nested models: a nugget effect and a spatial covariance function. Since the anisotropy ratio and the anisotropy angle parameters were not significant according to the Likelihood Based Parameter Estimation test, isotropic models were selected. In Table 4, the fitted variogram models are reported and described by the model name and the following parameters, namely nugget, partial sill and range. Actually, the Matérn model demonstrated to be the best suited theoretical function to describe spatial structure of the considered variables except for BD and PAWCe, since it was more flexible with respect to the classical models, having an additional shape parameter (k) (Figure 2). Following Cambardella et al. [56], the nugget-to-sill ratio can be a useful means to describe the spatial structure of the studied variables. In particular, the nugget semivariance expressed as a percentage of the total semivariance enables comparison of the relative size of the nugget effect among studied variables, with ratios <25% indicating strong spatial dependence, ratios between 25 and 75% moderate spatial dependence, ratios >75% weak spatial dependence. The analysis of the nugget-to-sill ratios (Table 4) indicated almost strong spatial dependence for TOC, Cl + Si and SSI (with values of 0.264, 0.272 and 0.263, respectively); these results were confirmed by the high predicted vs. observed correlations (r = 0.42, 0.63 and 0.74). RFCe, together with Ks, showed a weak spatial dependence, with nugget to sill ratios of 0.847 and 0.910, results confirmed also by not significant cross-validation outcomes (predicted vs. observed correlations) and by the findings of overall spatial correlation analysis (Moran statistics, Table 3). Yield fell between the moderate to weak classes, with a nugget to sill ratio of 0.761. PAWCe and BD presented a moderate adaptation as showed by the predicted vs. observed correlation (r = 0.49 and 0.46); the nugget-to-sill ratio was not computed because the nugget was not significant different from zero. The variogram estimated parameters for ACe appeared to be not physically based (an estimated value of 773 m was observed for the range). All the remaining range parameters values were consistent and below the maximum distance between observations (Table 4).  Spatial estimates of yield and spatially structured soil variables are reported in Figure 5. By visually inspecting the maps, a similar spatial behavior emerged for yield and BD with lower values in the left part of the experimental area and larger values in the lower right part. An inverse relationship was instead observed between yield and PAWCe. These results confirm the significant overall correlations already reported between yield and BD and PAWCe (Figure 4). TOC varied in a quite narrow range, apart from a few outliers, and the highest values were observed in the upper central part of the area. As expected, SSI map resembled TOC spatial behavior; at the same time, SSI brought the information related to the fine texture component (Cl + Si), showing as a consequence spatial similarity with the mentioned physical indicators maps.
To further deepen the SSI behavior observed, cross-correlograms were computed and cross-correlation coefficients at specific lags (0 m, 25 m, 50 m, 75 m, 100 m) were extracted (Table 5). SSI was strongly correlated with TOC and fine texture components, positively and negatively respectively, up to about 25 m distance for TOC and up to 50 m for Cl + Si. In addition, a moderate correlation was observed with PAWCe (positive) and BD (negative) up to 25 m, indicating the contribution of texture in driving this relationship. A weak correlation was observed between SSI and yield at lag 0, confirming that the behavior arose by the Pearson's overall correlation analysis ( Figure 4). Finally, the cross-correlation between SSI and RFCe was negligible, as expected.  Agronomic results of this investigation showed medium-high wheat yields on average (Table  1). However, according with the hypothesis that good conditions of soil physical quality should correspond to good yields [3,59], findings showed non optimal soil conditions in terms of TOC, SSI, and RFCe (Figure 3), i.e., in 43% of cases, and some considerations are due, particularly for optimal TOC ranges. For instance, the suggested lower critical TOC limits for agricultural soils (i.e., about 20 g kg −1 ) by Reynolds et al. [6] were obtained for urban soils, namely for sustainable establishment of plants in constructed landscaping soils used in urban parks, playing fields, curbside plantings, etc. [6,60]. Despite this critical threshold being applied in the literature e.g., [46,61,62], different limits should be considered for agricultural soils, and specifically for crops grown in Mediterranean environments. For example, TOC mean values observed in this investigation (i.e., about 16 g kg −1 ) were relatively high as compared to those reported in Table 1 by Ventrella et al. [63], for eight experimental farms and soils with different texture, located from north to south of Italy (5.2 ≤ TOC ≤ 15.3 g kg −1 ). In addition, in other investigations where a variety of plant species was considered, optimal TOC levels, as suggested in the literature (i.e., 30 < TOC < 50 g kg −1 ), were seldom reached. For instance, for approximately two-hundred soil cores collected in three areas and 18 spot sites of agricultural and forest Sicilian environments, Castellini and Iovino ([14]; Table 2) reported TOC values equal to 10.8 g kg −1 on average (0.9-37 g kg −1 as a measured interval). To provide more
By considering the role of BD as key soil indicator, cross-correlograms between the maps of BD and the other indicators were also computed ( Table 6). Strong negative correlations were observed with PAWCe up to about a 50 m distance. Moderate positive correlations were found with RFCe (up to 50 m) and yield (up to 25 m), whereas negative with TOC and SSI, as already observed, up to 25 m distance.

Discussion
Overall, the observed correlations showed physically plausible relationships among soil properties ( Figure 4). In particular, inverse relationships between ACe, TOC, SSI and BD were detected, as well as between ACe, Ks and RFCe. As expected, Ks was positively correlated with ACe, as increasing values of soil aeration should match with increasing saturated hydraulic conductivity (Ks); a reasonable positive relationship between Cl + Si fraction and Ks was also detected. Expected correlations have been verified between SSI and TOC (positive) or between SSI and Cl + Si (negative), since TOC and Cl + Si terms appear, respectively, in the numerator and in the denominator of Pieri's equation [6]. However, uncertain results were obtained for the significant correlations of PAWCe, because indirect relationships between PAWCe and Cl + Si, as well as BD, were detected. In fact, although the latter was quite verified in the literature for medium-high BD values [6,57,58], the former relationship does not seem entirely plausible because higher PAWCe values should correspond to higher contents of fine particles of the soil. A similar explanation can also be given for the relationship (both general and spatial) between PAWCe and crop yield.
Agronomic results of this investigation showed medium-high wheat yields on average (Table 1). However, according with the hypothesis that good conditions of soil physical quality should correspond to good yields [3,59], findings showed non optimal soil conditions in terms of TOC, SSI, and RFCe (Figure 3), i.e., in 43% of cases, and some considerations are due, particularly for optimal TOC ranges. For instance, the suggested lower critical TOC limits for agricultural soils (i.e., about 20 g kg −1 ) by Reynolds et al. [6] were obtained for urban soils, namely for sustainable establishment of plants in constructed landscaping soils used in urban parks, playing fields, curbside plantings, etc. [6,60]. Despite this critical threshold being applied in the literature e.g., [46,61,62], different limits should be considered for agricultural soils, and specifically for crops grown in Mediterranean environments. For example, TOC mean values observed in this investigation (i.e., about 16 g kg −1 ) were relatively high as compared to those reported in Table 1 by Ventrella et al. [63], for eight experimental farms and soils with different texture, located from north to south of Italy (5.2 ≤ TOC ≤ 15.3 g kg −1 ). In addition, in other investigations where a variety of plant species was considered, optimal TOC levels, as suggested in the literature (i.e., 30 < TOC < 50 g kg −1 ), were seldom reached. For instance, for approximately two-hundred soil cores collected in three areas and 18 spot sites of agricultural and forest Sicilian environments, Castellini and Iovino ([14]; Table 2) reported TOC values equal to 10.8 g kg −1 on average (0.9-37 g kg −1 as a measured interval). To provide more concrete examples, Table 7 summarizes other results of some investigations carried out in the Mediterranean environments of Apulia, Sicily and Sardinia (south-central Italy). With reference to such investigations, listed TOC values were equal, at most, to 16.7 g kg −1 for conservation soil management of durum wheat, equal to 22.6 g kg −1 for a citrus groove or slightly higher for a grassland, suggesting that it is not common to find considerably higher TOC levels in Mediterranean agricultural environments. On the other hand, TOC values close to the upper limit suggested by Reynolds et al. [6], 50 g kg −1 , were only detected for a sandy loam soil under high maquis (holm oak), i.e., when organic matter accumulation and probably slow organic matter mineralization rate can lead to such high levels ( Table 7). As an example, Figure 3 shows how the TOC classification by Sequi and Nobili [44], specifically referred to fine textured soils (clay, clay-loam, silty-clay and silty-clay-loam), is more suitable for the specific characteristics of Mediterranean agricultural soils. As a consequence, more accurate optimal or critical values should be provided for specific agro-environments, for example performing ad-hoc new research or by reviewing the available literature (e.g., through meta-analysis). Similar considerations can be made for SSI, as it is obtained from TOC and texture, and is not related directly to the soil structure, but is rather related to the soil resilience [6]. On the contrary, RFCe limits may be considered relatively more reliable as they were successfully verified in reference to other soil physical properties (i.e., BD, TOC, ACe, macroporosity) by Reynolds et al. [64] for a clay loam soil, as well as in comparison with literature guidelines. Among those investigated, only five soil physical properties were found spatially structured, i.e., Cl + Si, BD, TOC, SSI and PAWCe, but their relationships with grain yield did not appear to be always convincing ( Figure 5). In detail, in comparison to the yield map (i.e., lower or higher yielding zones), relatively coherent spatial distributions were identified in terms of fine soil texture components (Cl + Si) and BD; this suggests that, especially bulk density, can represent a reliable physical indicator to manage within-field sub-areas with different productivity. For this soil indicator, the map also shows that the sub-areas with relatively lower productivity correspond to those with very low BD values (as specified by the critical lower value of 0.9 g cm −3 of Figure 3), suggesting that this limit seems quite realistic for the case under study. On the contrary, conflicting information between yield and PAWCe maps were obtained (i.e., an inverse spatial structure), as already shown by the overall correlation. These results can be attributed to uncertainties of the PAWCe estimates obtained by BEST which would require further deepening. As expected, a relatively similar spatial distribution was observed for TOC and SSI, but not consistent with that of crop yield. In other words, findings provided by variables directly measured, both for overall correlation and for spatial analysis, were more convincing than those derived by variables estimated by BEST, for which further investigations are probably needed to quantify the accuracy degree of estimated soil water retention curve.
To deepen the aforementioned results and corroborate the relationships between yield and considered soil properties, a stepwise analysis was carried out to identify the variables most affecting wheat yield among those directly quantified (TOC, BD, fine textural components) or derived from laboratory measurements (SSI). Results of stepwise analysis showed that BD was the only variable selected (P = 0.0054), thus confirming the key role of soil porosity and compaction in affecting crop yield. Consequently, for investigated soil variables, only the BD map seems utilisable for agronomic management, i.e., for precision agriculture applications, because several factors may have contributed to produce unexpected or uncertain results for the other variables. Among these factors, particular relevance can be attributed to: (i) the uncertainties of the PAWCe estimates obtained by BEST; (ii) the relatively poor correspondence between observed physical-chemical fertility of the soil and agricultural yields in the specific conditions investigated; (iii) a possible spatial variability at a scale smaller than that experimentally measured (i.e., lower than about the mesh side). About the last statement, as Ks showed no significant spatial structure, we could make a plausible conjecture suggesting that spatial variability of Ks occurred at a smaller scale than that investigated [67]. However, no significant relationship between Ks and crop yield (or BD) was detected. Overall, Ks is reported to be a better soil structure indicator [68] than BD, as the latter does not provide any information on soil pore distribution, i.e., architecture, connectivity and tortuosity of soil porosity [3]. In our case study, a very good spatial correspondence between crop yield and soil bulk density was detected, and the cross-correlation analysis also showed a significant positive relationship with the crop yield for the lower distance, of about 25 m. Although this result could be considered questionable from the perspective of precision agriculture application (i.e., for implementation on platforms with on-the-go soil sensors), it is quite significant as BD measurements can be obtained more easily if compared to Ks, and it can be easily related to penetrometric soil measurements.
As a final remark, it should be noticed that SSI showed promising characteristics to be used as a representative indicator to identify homogeneous within-field areas, because of its tight relationships with chemical (TOC) and hydrological (fine texture components, and in turn PAWCe and BD) indicators and its predisposition to spatialization, thus suggesting a possible significant correlation with yield under a denser spatial sampling scheme. If this were the case, SSI would candidate itself as a key indicator for precision agriculture applications.

Conclusions
In this investigation a set of eight physical and hydraulic soil properties directly measured or estimated by BEST was obtained and spatialized to investigate correlations and identify intervals corresponding to medium-high levels of wheat yields, in order to provide useful information for site-specific agronomic management.
According to the guidelines of literature, a soil physical quality evaluation highlighted that soil under study had optimal bulk density, plant available water capacity, air capacity and saturated hydraulic conductivity values, but that the total organic carbon content, structure stability index and relative field capacity suggested too low levels of organic carbon or excessive soil compaction. However, both a literature analysis for different types of Mediterranean vegetation cover and the correlation analysis (overall and spatial correlation) suggested that a review of the optimal or critical TOC values for typical crops of Mediterranean environments should be made, as TOC value around 20 g kg −1 is hardly achievable even under conservation agriculture systems, and literature values are probably not entirely realistic for Southern Italy's cereal crops. In addition, capacitive indicators estimated from the soil water retention curve of BEST provided both expected and uncertain correlations, especially with regard to the inverse relationship between plant available water capacity and wheat yield. Therefore, for the soil-crop system studied, application of the simplified BEST-procedure did not return completely reliable results and further investigations are needed to quantify the accuracy and reliability of estimated soil water retention curve. These main results open up new research perspectives to improve our knowledge on this topic.
Among measured soil properties, BD showed a spatial distribution in agreement with that detected for crop yield, and the cross-correlation analysis also showed a significant positive relationship only for short lags. Finally, SSI showed promising characteristics suggesting a possible significant correlation with yield under a denser spatial sampling scheme, and a potentiality as a key indicator for precision agriculture applications.
Further research on this topic is needed for Mediterranean agro-environments, by deepening: (i) the reliability of available measurement methods for accurately estimating representative physical and hydraulic soil properties, and (ii) the temporal stability of observed spatial relationships between soil properties (soil bulk density or soil texture) and crop yields along a larger time interval.