Scale-Dependent Field Partition Based on Water Retention Functional Data

Abstract

Functional data are being used increasingly in recent years and in many environmental sciences, such as hydrology applied to agriculture. This means that the output, instead of a scalar variable represented by a spatial map, is given by a function. Furthermore, in site-specific management, there is a need to delineate the field into management areas depending on the agricultural procedure and on the scale of application. In this paper, an approach based on multivariate geostatistics is illustrated that uses the parameters of Dexter’s water retention model and some soil properties to arrive at a multiscale delineation of an agricultural field in Iran. One hundred geo-referenced soil samples were taken and subjected to various measurements. The volumetric water contents at the different suctions were fitted to Dexter’s model. The estimated curve parameters plus the measurements of the soil variables were transformed into standardized Gaussian variables and the values transformed were subjected to geostatistical cokriging and factorial cokriging procedures. These results show that soil properties (organic carbon, bulk density, saturated hydraulic conductivity and tensile strength of soil aggregates) influence the parameters of Dexter’s model, although to different extents. The thematic maps of both soil properties and water retention curve parameters displayed a varying degree of spatial association that allowed the identification of homogeneous areas within the field. The first regionalized factors (F1) at the scales of 508 m and 3000 m made it possible to provide different delineations of the field into homogeneous areas as a function of scale, characterized by specific physical and hydraulic properties. F1 at a short and long distance could be interpreted as “porosity indicator” and “hydraulic indicator”, respectively. Such type of field delineation proves particularly useful in sustainable irrigation management. This paper emphasizes the importance of taking the spatial scale into account in precision agriculture.

1. Introduction

Agricultural soils are very complex systems that are capable of storing and transmitting air, water, solutes and heat to microorganisms and plants. Soil consists of solid spaces (particles) and voids (pores) that determine its structural characteristics and function. In particular, water transport, which is crucial in controlling the health and productive potential of a plant, occurs through the connected pore spaces, whereas the workability of soils depends on the existence of pores in the form of microcracks. Pore spaces are determined by the size of composite particles that form a hierarchical system [1,2], consisting of various levels: groups of primary particles comprising microaggregates, microaggregate groups comprising aggregates and aggregate groups comprising clods, which are aggregates formed by tillage according to SSSA Soil glossary, or bulk soil. Similarly, the average pore size follows a hierarchical system determined by the size of the soil particles. Therefore, following [3] the total porosity of soil can be said to contain various components: residual, matrix, structural and macro porosity.

The first and last components include pores that are either too small or too large, respectively, that are generally disregarded in experiments aimed at defining optimal soil water management for irrigation purposes. In contrast, the matrix porosity represents the pore space between individual mineral particles in the soil. As groups of mineral particles can combine to form discrete microaggregates, it is also necessary to consider structural porosity, which is the pore space between the microaggregates and between the incipient aggregates. This type of porosity in most soils is largely composed of microcracks, also called ‘elongated transmission pores’ by Pagliai, Vignozzi and Pellegrini [4], which play an important role in transport processes.

Moreover, it is well known that water in pores is not in a free state, but is subject to a matrix potential whose absolute value is called suction. The relationship between suction and pore water content can be described by a function called soil water content retention curve (SWRC) [5], the shape of which depends on different physical and chemical properties of the soil, which are specific to each soil [6]. SWRC contains useful information on pore water content at different suctions and pore size distributions, which characterize the response of the soil to stress conditions, such as saturation or dryness [7]. Therefore, it is clear that knowledge of it cannot be disregarded for sustainable irrigation water management in a spatially variable agricultural field. However, there are many models of SWRC in the literature that have been compared by various authors [8,9]. In particular, Dexter’s model [3] is a double-exponential equation with five adjustable parameters and a bimodal distribution, and each mode is related to either matrix or structural porosity. Compared with the other models, it has the following advantages: all five parameters have a physical meaning related to the different components of porosity; the parameter related to residual porosity is always positive or zero; the five parameters are poorly correlated with each other and can be related to soil properties influencing its hydraulic behaviour.

Dexter’s model parameters together with other relevant properties of agricultural soil can then be treated jointly in a spatial analysis to arrive at a partitioning of the field into homogeneous areas to be subjected to differential irrigation management. The study of the spatial variability of an agricultural field generally requires soil sampling at a finite number of locations. The samples are then subjected to laboratory measurements and these measurements are processed spatially to produce thematic maps. However, in the case discussed above, the interest is not in one or more scalar variables but in a functional variable; i.e., in the modelling of a water retention curve from a set of soil moisture measurements recorded at the same location and corresponding to gradually increasing values of water suction.

Moreover, it is expected that the functional data will exhibit spatial dependence. Several alternative techniques have been developed by different authors [10] to spatially predict continuous curves instead of single-point data for each observation site. These methods first require the application of a parametric or non-parametric smoothing technique (e.g., Dexter’s model or B-spline regression) to convert discrete data into continuous functions; after that, some interpolation technique is applied to predict the whole curve in unsampled locations. These methods have also been used in several areas besides hydrology and agronomy, including meteorology [11] and marine biology.

In this paper, a multivariate geostatistical approach was preferred for the following main reasons. The fitting of a linear model of coregionalization (LMC) allows the modelling of the spatial relationships of Dexter’s retention curve parameters with each other and with the soil properties considered to be particularly influencing water transport. In addition, multivariate geostatistics make it possible to arrive at a delineation of the field in homogeneous areas that is a function of scale [12,13], since the processes governing soil water transport occur over different scales.

The delineation of management zones for precision irrigation is commonly based on statistical clustering, which combines the impacts of different variables operating over different spatial and temporal scales on the hydraulic behaviour of soil [14,15,16,17]. However, these techniques suffer from two main limitations: (1) they do not take into account the generally observed spatial correlation between spatial observations, and (2) they produce a single field partition at an unspecified spatial scale, even if this partition is produced by processes operating at different scales [18].

Therefore, in light of the above considerations, the objective of this study was to define a multivariate geostatistical approach based on SWRC measurements in an agricultural field to predict the whole curve of Dexter’s model at unsampled points and produce a multiscale field delineation in homogeneous areas for site-specific water management.

	Notation Box
KS	Saturated hydraulic conductivity (ms−1)
OC	Organic carbon (℅)
BD	Bulk density (g cm−3)
TS	Tensile strength of soil aggregates (kPa)
θr	Asymptote of the equation, which is residual water content (cm3cm−3)
A1	Amount of matrix pore spaces (cm3cm−3)
A2	Amount of structural pore spaces (cm3cm−3)
h1 and h2	Characteristic pore water suctions at which the matrix and structural pore spaces start to empty, respectively (bar)

2. Materials and Methods

2.1. Study Area and Experimental Site

The experiment was carried out in a 200-ha field at the Dastjerd Research Farm of the College of Agriculture, Bu-Ali Sina University, Hamedan, Iran (35°02′10.46″ N–35°03′39.72″ N and 48°28′21.25″ E–48°30′0.90″ E) in 2019. One hundred disturbed and undisturbed soil samples were collected from 0.10–0.15 m depth at the nodes of a grid with cells of 150 × 150 m² size and in some areas of 75 × 75 m² and 37 × 37 m² size [13] (Figure 1) to increase the number of pairs of samples at shorter distances and to reduce the nugget effect in variogram models, representing the lower limit of estimation precision.

The mean annual rainfall is 479.7 mmy⁻¹ and the mean annual temperature is 14.83 °C (daily temperature varying from −9 °C to 35 °C).

2.2. Dexter’s Soil Water Retention Curve (Dexter’s Model)

Soil water content was measured in the undisturbed core samples at the matric potentials of −0.01, −0.02, −0.04 and −0.06 (with a Sandbox apparatus), and −0.3, −1, −5 and −15 bar (with a pressure plate apparatus).

Dexter’s model (Dexter et al., 2008) was fitted to data:

$$\mathrm{w}={\mathsf{\theta }}_{\mathrm{r}}+{\mathrm{A}}_1\exp \left(\frac{-\mathrm{h}}{{\mathrm{h}}_1}\right)+{\mathrm{A}}_2\exp \left(\frac{-\mathrm{h}}{{\mathrm{h}}_2}\right)$$

(1)

where, w (cm³cm⁻³): volumetric moisture content of the soil; θ_r: asymptote of the equation, which is the residual water content (cm³cm⁻³); A₁ (cm³cm⁻³): amount of matrix pore spaces; A₂ (cm³cm⁻³): amount of structural pore spaces; h (bar): matric suction corresponding to w; h₁ (bar) and h₂ (bar); characteristic pore water suctions at which the matrix and structural pore spaces start to empty, respectively [3].

Model Fitting

Dexter’s model parameters were estimated using the Excel Solver (2016). To test the goodness of fitting, root mean square error (RMSE) and determination coefficient (R²) were used as follows:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{\mathrm{N}}\sum _{\mathrm{i}=1}^{\mathrm{N}}({\mathrm{y}}_{\mathrm{i}-}{\mathrm{y}}_{\mathrm{P}\mathrm{i}}{)}^2}$$

(2)

$${\mathrm{R}}^2=\left[\frac{\left[\mathrm{N}({\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{Y}}_{\mathrm{i}}{{\mathrm{Y}}_{\mathrm{P}}}_{\mathrm{i}})-{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{Y}}_{\mathrm{i}}{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{{\mathrm{Y}}_{\mathrm{P}}}_{\mathrm{i}}\right]}{\sqrt{\left[\mathrm{N}({\sum }_{\mathrm{i}=1}^{\mathrm{N}}{{\mathrm{Y}}_{\mathrm{i}}}^2)-({\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{Y}}_{\mathrm{i}}{)}^2\right]\times \left[\mathrm{N}\left({\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{Y}}_{{\mathrm{P}}_{\mathrm{i}}}^2\right)-({\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{Y}}_{{\mathrm{P}}_{\mathrm{i}}}{)}^2\right]}}\right]$$

(3)

where, N is the number of samples; Y_i and Y_pi are the measured and predicted values of the variable at location i, respectively.

2.3. Soil Property Measurements

Disturbed soil samples (2 kg) were air-dried and passed through a 2 mm sieve before physical and chemical analyses (Table 1).

Table 1. Method of measuring study variables.

Variable	Measurement Method	Unit	References
Soil organic carbon (OC)	oxidation method	(%)	[19]
Bulk density (BD)	standard core method	(gcm−3)	[20]
Saturated hydraulic conductivity (KS)	The constant head method	(ms−1)	[21]

Table 1. Method of measuring study variables.

Variable	Measurement Method	Unit	References
Soil organic carbon (OC)	oxidation method	(%)	[19]
Bulk density (BD)	standard core method	(gcm−3)	[20]
Saturated hydraulic conductivity (KS)	The constant head method	(ms−1)	[21]

The studied soils covered five textural classes according to the USDA soil classification (Figure 2) [13]. Soil samples included the following classes: clay loam (n = 6), loam (n = 10), sandy clay loam (n = 57), sandy loam (n = 22) and loamy sand (n = 5).

2.4. Aggregate Tensile Strength Measurement

According to Dexter and Kroesbergen (1985), tensile strength is probably the most useful measure of the strength of individual soil aggregates, where soil strength is intended to the highest stress that any type of soil can be exposed to before failure [22]. It is influenced by soil moisture, the degree of connection among soil aggregates and the occurrence of soil cracks [23]. An increase in soil moisture causes a decrease in the bonding strength between the organic and inorganic particles in soil aggregates and a reduction in the tensile strength of soil aggregates [24]. Therefore, given the close relationship between the soil water content and tensile strength, it was decided to include this soil property among the variables chosen for the delineation of soil in hydraulically homogeneous areas. To measure tensile strength, air-dried aggregates with a diameter of 4 to 8 mm were used (Figure 3).

The reason for using these diameters was the ease of collection and measurement with the available devices. First, twenty aggregates were randomly selected from each soil sample; second, the device described by Dexter and Watts [25] was used to measure the tensile strength (Figure 4). In this device, soil aggregates were placed between the two disks, and pressure was applied to them by adding water to a water container and using a leverage system. Before starting the test, the disk and container were levelled so that no pressure was applied to the aggregates. Pressure was then applied to the aggregates until they broke, at which point the water uptake stopped.

The force (F) in Newton (N) required to break aggregates was calculated using the following equation:

$$\mathrm{F}=\left[{\mathrm{M}}_{\mathrm{c}}\mathrm{g}\times \left(\frac{{\mathrm{X}}_1}{{\mathrm{X}}_2}\right)\right]$$

(4)

where, M_c is the mass of water (kg) required to break soil aggregates; X₁ is the length of the upper lever; X₂ is the length of the lower lever (cm); g is gravity acceleration (9.807 ms^–2).

Finally, the tensile strength of aggregates was calculated using the following equation based on the method by Rogowski [26] and Dexter and Kroesbergen [27]:

$$\mathrm{Y}=0.576\times \frac{\mathrm{F}}{{\mathrm{d}}_{\mathrm{e}\mathrm{f}\mathrm{f}}^2}$$

(5)

where, Y is tensile strength (kPa); F is the maximum force (N) required for aggregate breaking; d_eff is the effective diameter of the aggregates assumed to be of spherical shape, which was calculated using the following equation:

$${\mathrm{d}}_{\mathrm{eff}}=\frac{{\mathrm{d}}_{\max }+{\mathrm{d}}_{\min }}{2}$$

(6)

d_max and d_min are the maximum and minimum diameters of soil aggregates.

2.5. Multivariate Geostatistical Analysis

The proposed multivariate geostatistical approach jointly treats two different types of data (data fusion): point data concerning the measurements of some soil properties (KS, BD, OC and TS), which were considered to be influencing the hydraulic characteristics of soils, and functional data for the determination of the parameters of Dexter’s model. The objective of such an analysis is two-fold: (1) to determine SWRC at an unsampled point in the field; (2) to characterize the field from a hydraulic point of view to arrive at a partition into homogeneous areas to be subjected to differential irrigation water management.

The approach consists of different steps as follows:

• Fitting of Dexter’s model to determine its five parameters;
• Basic statistics of all study variables including Dexter’s model parameters and the selected soil properties;
• Normalization and standardization of all variables through Gaussian anamorphosis, if they show significant deviations from symmetrical data distribution;
• Fitting of a linear model of coregionalization (LMC) to the whole data set of the experimental direct and cross-variograms of the Gaussian transformed variables to model the spatial relationships of the curve parameters with each other and with soil properties;
• Testing the goodness of fitting through cross-validation;
• Ordinary cokriging on all Gaussian transformed variables and back transformation of the estimates of Gaussian transformed parameters using the anamorphosis Gaussian function previously determined;
• Mapping of soil properties for characterizing spatial variability, and of the curve parameters to draw SWRC in any point of the field;
• Factorial cokriging on Gaussian transformed variables and extraction of the regionalized factors at short and long ranges;
• Mapping of the retained regionalized factors at each scale with eigenvalue greater than 1, which cumulatively explains a good proportion of the spatial variability at that scale;
• Delineation of within-field zones with different hydraulic properties by using the retained regionalized factors.

With regard to the geostatistical procedures used, the interested reader can refer to one of the many published manuals [28] for a detailed description. Only the main steps are briefly described below (Figure 5).

2.5.1. Basic Statistics and Gaussian Anamorphosis

The maximum, minimum, mean and median values, standard deviation (SD), skewness and kurtosis were calculated for all study variables to characterize the type of data distribution. In the case of skewed variables, Gaussian anamorphosis was performed to determine a continuous mathematical function that transformed a standardized Gaussian variable into a random variable [29]. For such a transformation, Hermite’s ortho-normal polynomial series was truncated to the first 100 terms.

2.5.2. Spatial Interpolation

• Linear Model of Coregionalization

The LMC represents a simplified model of multiscale spatial relations in a multivariate data set [30,31,32]. It assumes that all direct and cross semi-variograms can be described by the same set of basic spatial structures (standardized semivariograms with unit sill) for each spatial scale (range). What differentiates variables or pairs of variables is only the partial sill of each spatial structure corresponding to a specific scale. Although quite crude, this model nevertheless works well under many experimental conditions [30] since it is based on the assumption that the variables are different manifestations of the same processes operating over different scales in the system under study. Cokriging and Factorial cokriging are two of the standard multivariate geostatistical techniques using the same LMC.

• Factorial Cokriging Analysis

Factorial cokriging analysis (FCA) is a technique very similar to the principal component analysis (PCA) of traditional statistics. The fundamental difference is that instead of being applied to the variance-covariance matrix, it is applied to each coregionization matrix of a given spatial scale and derived from the LMC. The latter can be considered as a variance-covariance matrix, but specific to a given spatial scale. As in PCA, orthonormal components can be extracted, which in geostatistics are called regionalized factors and are scale-dependent. From the loading coefficient values, it is then possible to attribute a particular physical meaning to each factor, which can be regarded as a synthetic indicator of spatial variability. All geostatistical analyses were performed using the Isatis geostatistical software package [33].

• Cross-validation

The quality of the LMC can be evaluated using cross-validation [34], which consists of removing one observation at a time and then producing the estimate at the point where the observation was removed with the remaining observations. The difference between the estimate and the observation represents the error. Repeating the procedure for all n observations yields n errors for which different statistics can be calculated.

In particular, the mean error (ME), root mean square error (RMSE) and root mean square standardized error (RMSSE) were calculated using the following equations:

$$ME=\frac{{\displaystyle \sum _{i=1}^n}\left(\stackrel{\wedge }{Z}\right(x_i)-z(x_i\left)\right)}{n}$$

(7)

$$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(\stackrel{\wedge }{Z}\left(x_i\right)-Z\left(x_i\right)\right)}^2}}{n}}$$

(8)

$$\begin{array}{l}\\ RMSSE=\sqrt{{\frac{{\displaystyle \sum _{i=1}^n\left[\left(\stackrel{\wedge }{Z}\left(x_i\right)-z\left(x_i\right)\right)/\stackrel{\wedge }{\sigma }\left(x_i\right)\right]}}{n}}^2}\end{array}$$

(9)

where, n is the number of samples; $\stackrel{\wedge }{Z}$ (x_i) and z (x_i) are the predicted and measured values of the variable at point x_i, respectively; $\stackrel{\wedge }{\sigma }\left(x_i\right)$ is the root mean squared prediction error. The optimal values for ME and RMSE should be 0, whereas RMSSE should be close to one within the heuristic tolerance interval $\left[1-3\sqrt{\frac{2}{n}},1+3\surd \frac{2}{n}\right]$ [35].

3. Results and Discussion

3.1. Basic Statistics of Soil Properties

The basic statistics of the soil properties are presented in

Table 2. Basic statistics of soil properties.

	Min	Max	Mean	Std. Dev	Skewness	Kurtosis
KS (ms−1)	0.33	34.59	8.73	5.93	1.58	7.04
OC (%)	0.08	0.66	0.40	0.12	−0.51	3.20
BD (gcm−3)	1.21	1.73	1.42	0.11	0.50	2.98
TS (kPa)	51.47	124.21	85.71	17.05	0.21	2.32

. Related to variations in soil texture [13] and structure (Table 2), 40% of the soil samples had a very high KS (>0.12 m/h), 23% had high KS (0.06–0.12 m/h), 29% had moderate KS (0.02–0.06 m/h), 4% had low KS (0.01–0.02 m/h) and 4% had very low KS (<0.01 m/h) [36]. For low KS, runoff is quite irregular and soil may be suitable for irrigation [36], whereas for very high KS, (about 40% of soil samples), if soil is used for effluent or waste disposal, possible problems might arise due to potential groundwater contamination. Only 8.9% of soil samples had high BD (1.6–1.7 gcm⁻³), whereas the others had medium BD (1.3–1.6 gcm⁻³) [26]; therefore, the permeability of the soil does not cause water runoff or irrigation problems, at least in most circumstances. The box plots reveal the existence of some outliers in the soil properties, with the exception of TS (Figure 6). Strong positive skewness and kurtosis > 3 were observed for KS; however, all variables showed non-symmetrical data distribution.

3.2. Basic Statistics of Dexter’s Water Retention Curve Model (SWRC)

The goodness of fitting for the SWRC function to Dexter’s model was satisfactory at any location with an average RMSE of 0.010 (range: 0.001–0.08) and R² varying between 0.96 and 0.99. The values of model parameters showed general wide variability due to the heterogeneity of soils belonging to different texture classes (Figure 2). High skewness in absolute values (>0.5) was observed in A₂ and h_2, indicating markedly non-symmetrical data distribution of these parameters (

Table 3. Basic statistics of Dexter’s model parameters.

	Min	Max	Mean	Std. Dev	Skewness	Kurtosis
θr (cm3cm−3)	0.00	0.13	0.05	0.05	0.22	1.39
A1 (cm3cm−3)	0.06	0.25	0.15	0.05	0.10	2.07
A2 (cm3cm−3)	0.10	0.34	0.24	0.05	−0.72	3.32
h1 (bar)	0.27	43.51	14.74	10.91	0.38	2.00
h2 (bar)	0.01	0.19	0.06	0.03	1.46	10.04

).

The box plots (Figure 7) show the presence of outliers only for A₂ and h₂, which confirms the high asymmetry of the data distribution of these parameters as already indicated by the skewness coefficient.

Although the geostatistics estimator is not parametric and, in theory, can be applied to any type of data distribution, variogram modelling is facilitated on Gaussian transformed data. Furthermore, it is also possible to calculate the uncertainty of strongly skewed data [37]. Therefore, before using multivariate analysis, it was preferred to normalize and standardize to zero mean and unit variance for all data (including soil properties and model parameters) due to the presence of skewed variables and outliers.

3.3. Correlation Analysis

In order to detect significant linear correlations between the model parameters and the investigated soil properties, Pearson coefficients were calculated for all variables (

Table 4. Correlation coefficients (r) between the study variables.

	θr	A1	A2	h1	h2	KS	OC	BD	TS
θr	1.00	−0.72 **	−0.06	−0.91 **	−0.27 **	0.18	−0.11	−0.18	−0.01
A1		1.00	−0.01	0.60 **	0.07	−0.11	0.29 **	−0.06	0.08
A2			1.00	0.11	0.19	0.16	0.29 **	−0.38 **	0.12
h1				1.00	0.41 **	−0.22 *	0.11	0.14	0.10
h2					1.00	−0.18	0.00	−0.05	0.01
KS						1.00	0.10	−0.29 **	−0.15
OC							1.00	−0.29 **	0.23 *
BD								1.00	−0.06
TS									1.00

* and ** Correlations are significant at the 0.05 and 0.01 probability levels, respectively.

).

A₁ and A₂ showed significant and positive correlations with OC (0.29 **), while A₂ had a negative correlation with BD (−0.38 **). These results are expected because organic carbon is effective in increasing the porosity at the pore sizes investigated and decreasing the bulk density. However, BD seems to be more influenced by structural porosity than by matrix porosity, at least according to Dexter’s model. h₁ was negatively correlated with KS (r = −0.22 **), which reflects the fact that high matrix pore water suction has a negative impact on KS. With regard to the relationships between the model parameters, it is worth emphasizing the negative correlation of θ_r with A₁ and h₁ and partly with h₂, which shows that soils with a high proportion of residual water have a low proportion of matrix pores. Additionally, other significant observations are as follows: the positive correlation between A₁ and h₁ (0.60 **); the higher the pore water suction at which the matrix pore spaces starts to empty, the higher the proportion of matrix pores; the positive correlation (0.40 **) between the two suctions related to the two sizes of pores. Concerning the significant correlations between soil properties, the positive correlation (0.23 *) of OC with TS is explained because OC increases the soil’s resistance to being crushed. Saturated hydraulic conductivity and OC were negatively correlated with BD (−0.29 **), which was as expected because OC increases porosity, improves KS and reduces BD [38].

It can be deduced from the above results that, although some correlations are significant, there are no strong relationships between the investigated soil properties and the parameters of Dexter’s model. However, these are linear correlations that do not take into account the spatial correlations of the observations. The spatial analysis could confirm these results.

An isotropic LMC was fitted to the matrix of both the direct and cross-experimental variograms of all Gaussian transformed variables, including the following basic spatial structures: (S0) nugget effect, (S1) spherical model (range = 508.64 m) and (S2) spherical model (range = 3000 m) (Figure 8a,b). The proportions of the total spatial variance explained by each spatial component were as follows: S0 50%, S1 13% and S2 37%, which shows that most spatial variation is due to the nugget effect, i.e., spatially uncorrelated error. One way to reduce this spatial component might be to intensify sampling because at the scale at which it was performed, it proved incapable of resolving micro-variability due to both intrinsic soil properties and anthropogenic activities.

From Figure 8a of the LMC, some considerations can be made as follows: all direct variograms show a high nugget, as already noted in general above, which can be attributed to the adopted sampling scale, which should be reduced, as well as the measurement error associated with the measurement of soil properties and parameter estimation error. In addition, all direct variograms clearly show a short-range structure (around 500 m), which can be attributed to intrinsic soil variation within the field and the effects of anthropogenic operations, and a long-range structure (around 3000 m), causing the sill to exceed the sampling variance, which can be attributed to the study area belonging to soil map units extending beyond the field size [13]. Deviating from the characteristics stated above, there is a KS variogram, which is similar to the pure nugget effect, most likely determined by the low precision of its measurement marked by high stochasticity; and that of BD, in which the long-range component prevails since BD is generally closely linked to the type of map unit to which the soil belongs at a scale greater than that of the field (Figure 8a).

Regarding the cross-variograms between the model parameters and soil properties, it can be said that in general the spatial relationships are quite low. This can be ascertained from the large distance of the model curve of the cross-variograms from the dashed line (Figure 8b) representing the intrinsic correlation, i.e., the maximum obtainable spatial correlation between the two variables considered [32]. The only soil properties that show spatial correlations of some significance with the model parameters are OC, which is positive with A₁ and A₂ due to the positive effect already noted for OC on both matrix and structural porosity; and BD, which is negative with θ_r and A₂ but positive with h₁, meaning that an increase in density causes a reduction in structural porosity and the amount of residual water, but an increase in the suction at which matrix pores become empty. Regarding the spatial relationships between the various parameters of the model, from the cross-variograms, it can be seen that θ_r is negatively correlated with all but A₂; A₁ is positively correlated with h₁ and, to a lesser extent, with h₂. It is noteworthy that A₁ and A₂ are not significantly correlated; that is, there is no significant spatial relationship between the two types of porosity (Figure 8b). This is probably attributable to the different origins of the two types of porosity; A1 (matric porosity) is intrinsically linked to the granulometric characteristics of the soil (texture), whereas A2 (structural porosity) is mostly related to the type of tillage carried out in the field. Moreover, A2, similarly to the other parameters of Dexter’s model, shows a high degree of stochasticity (nugget effect) and is the one that appears least correlated with the other parameters (Figure 8b).

There is a positive relationship between h₁ and h₂ and between h₁ and A₁; therefore, higher porosity, at least matrix porosity, is associated with higher suction, at which pores (of both types) begin to empty.

Finally, with regard to the spatial relationships among the soil properties investigated, KS appears to be poorly correlated with all of them, mostly due to its essentially stochastic nature. It shows a slight negative correlation with BD, which is quite expected since the densest soils are those with the lowest porosity and, therefore, with low hydraulic conductivity. TS is spatially uncorrelated with the other soil properties, except for a somewhat positive correlation with OC. Therefore, it can be said that both of these properties (OC and TS) contributed to improving the soil structure and consequently its hydraulic characteristics.

Summarising the previous results concerning (1) the high nugget effect for direct variograms, and (2) the general low sill of cross-variograms, denoting low spatial correlation between the curve parameters and between these and the soil variables, it can be stated that this is attributable to two main causes. The first is related to the insufficient sampling density, as already noted; the second is to the uncertainty inherent in the measurements of the soil variables (especially KS) and to the fact that the parameters of Dexter’s model are not directly measurable physical variables but are estimated parameters with an estimation error that may also be high. However, despite the imprecision associated with the measurement and estimation procedures, we believe that it is possible to disclose interesting spatial relationships among variables related to the soil’s capability of transmitting water.

The goodness of fit of the LMC for any Gaussian transformed variable was evaluated through cross-validation, and the results are shown in

Table 5. Results of cross-validation of the studied variables.

	gθr *	gA1	gA2	gh1	gh2	gKS	gOC	gBD	gTS
ME	−0.03	−0.03	0.03	0.01	0.01	0.03	0.03	0.002	0.05
RMSE	−0.06	−0.06	0.04	0.03	0.01	0.03	0.04	0.004	0.07
RMSSE	0.95	1.00	0.85	0.99	0.98	0.98	0.85	1.08	1.10

* g before the name of variable stands for Gaussian transformed variable.

. The ME and RMSE were quite close to 0, and the RMSSE fell within the tolerance interval [0.58–1.42] for all variables.

The maps of soil attributes display that higher KS and OC values were localized in the central parts of the field, whereas lower values were in the south and south eastern corner (Figure 9a,b). Thus, the OC and KS showed a positive degree of association, and their higher values were associated with higher clay contents as reported elsewhere [13].

Conversely, BD was higher in the southeastern part of the field, where the soils were mostly sandy (Figure 9c) [13]. In the north and central parts, the presence of more OC enhanced the formation of larger soil pores, which reduced BD and improved permeability. Very low KS (0.05–1 cm/h) in the south part of the field, which might cause regular runoff under rainfall and irrigation, is likely to be inefficient. In contrast, in the central part of the field and the north, with hydraulic conductivity greater than 6 mm per hour, superficial water flow would rarely occur, and the soil might be too permeable for irrigation [36].

The aggregate tensile strength map showed a different trend in the field, which is more similar to that of BD (Figure 9d). The lowest values of TS were observed in the central and southeastern parts of the field, probably due to an increase in OC, which might decrease the tensile strength because it might increase the size of the aggregates [39]. Smaller aggregates are stronger because they contain fewer surface cracks [40] and are denser [2]. However, the relationship between OC and TS is not always so direct as it appears to be reversed in some southern parts of the field. This would induce one to think that the resistance of aggregates to the disintegrating action of water actually depended on multiple interrelated factors.

The increase in sand in the southeastern part [13] reduced the tensile strength of aggregates, probably because sand acted as the main nucleus to connect the particles to each other and form larger aggregates, which were more unstable, thereby causing soil fragility to increase [41].

Greater θ_r values were localized in the central part of the field, whereas smaller values of these parameters were observed in the north and south (Figure 10a). Moisture retention at greater suctions (θ_r) strongly depends on the particle size distribution [42]; therefore, a finer texture in the central part of the field might have caused an increase in θ_r.

A₁ showed the lowest values in the southern part of the field, whereas the highest values were observed in the central and northern parts of the field, probably due to the silt content, which is mostly responsible for the pores of the soil matrix [43]. A₂ map shows a very similar pattern to A₁ map with the exception of the northern part, where the structural porosity is reduced (Figure 10c) due, as already noted, to the higher silt content in soils to the north. However, throughout the field, the proportion of structural pores prevails over that of matrix pores; that is, the larger pores between the soil aggregates can favour water flow and avoid water logging.

However, a correct proportion between macro- and micro-porosity is fundamental for rational water management in an agricultural field because, while the former facilitates water flow and avoids waterlogging, the latter increases water availability to crops. This critical balance can be achieved by adopting, for example, conservation agriculture techniques [44] that improve the sustainable use of water resources. The higher OC content in the central part has a positive effect, especially on structural pores (A₂, Figure 10c), which improves the hydraulic properties and, in particular, promotes the water flow near soil-saturated conditions (KS, Figure 9a). Therefore, from the above discussion, the central part seems to have water characteristics that make it suitable for irrigation. In contrast, the greater values of sand and BD (smaller values of total porosity) in the southern part of the field would have caused a loss of structural and matrix porosity, which explains the reduction in A₁ and A₂. Therefore, the plants in this part of the field might be more frequently susceptible to water stress than elsewhere.

The maps of h₁ and h₂ (Figure 10d,e) show a somewhat opposite trend to that of A₁ and A₂, whereby greater (both structural and matrix) porosity generally corresponds to lower suctions. The highest values of h₁ and h₂ are observed in the south, with larger contents of sand. However, the relationship between porosity and suction is not so direct and is greatly influenced by the variations in texture and pore size variation spectra that occur in the field. For example, loam sandy soils (sands and sandstones) located in the northeast and southwest have a very narrow distribution of pore sizes; therefore, they empty over a very narrow range of suctions [45]. Furthermore, tillage management can affect porosity and suction at which pores become empty because soils under conventional tillage produce higher structural porosity (large and medium pores) and lower matrix porosity (small pores) than soils under no tillage or conservative agriculture [46,47], thus reducing the risk of water stress. Therefore, more conservative and less energy-demanding techniques should be promoted, as the field was only subjected to conventional tillage at least until the sampling date.

The maps in Figure 9 can also be used to derive SWRC curves at each point in the field. As an example, the curves at three locations in the field (Figure 9a) that clearly differed in their hydraulic properties, as discussed above, are shown (Figure 9a).

The hydraulic differences among the three zones are also well shown by the SWRC curves; in fact, at equal suction, the highest values of water content occur in the middle zone, the lowest in the southeast zone and the intermediate in the north zone (Figure 11a). These variations, as discussed above, can clearly be interpreted in terms of textural variations (higher content of silt in the north, clay in the center and sand in the south-east) and OC (higher content in the central zone), which have an impact on bulk density and porosity and, thus, on water availability to plants.

3.4. Factorial Cokriging

To synthesize the spatial relationships shown above between Dexter’s model parameters and soil properties, and simultaneously arrive at a delineation of the field into hydraulically homogeneous areas, a factor analysis was carried out. Only the first factor (F1 short) at the shorter range (508.64 m) and the first and second factors (F1 long and F2 long) at a longer range (3000 m) were retained, with eigenvalues greater than 1, thereby explaining a larger percentage of the variance than that explained by each variable standardized to variance 1. The percentages of spatial variance at the corresponding scale explained by F1 short, F1 long and F2 long were 73.72%, 68.79% and 24%, respectively (

Table 6. Composition of the retained regionalized factors at the scales of 508.64 m and 3000 m.

	gθr*	gA1	gA2	gh1	gh2	gKS	gOC	gBD	gTS	Eigen val.	Var. perc.
Spherical	Model Range	508.64 m
F1short	0.30	−0.27	−0.20	−0.56	−0.49	0.06	−0.30	0.13	−0.37	1.08	73.72
Spherical	Model Range	3000 m
F1 long 1	0.35	−0.04	0.36	−0.36	−0.18	0.31	0.40	−0.57	−0.05	2.95	68.79
F2 long	−0.24	0.36	0.36	0.22	0.30	0.19	−0.01	−0.03	−0.71	1.03	24.00

* g before the symbol of variable means Gaussian transformed variable.

).

On F1 short, the Gaussian variables were as follows: gh₁, gh₂, gTS and gOC, and to a less extent gA1 and gA₂ weigh negatively, whereas gθ_r weighs positively. Therefore, this factor may be considered an indirect “porosity indicator” associated with higher values of OC and TS. It is worth underlining that, at this scale, the porosities and suctions of both types are positively correlated. Therefore, low values of F1 short are an indicator of soils with high total porosity, greater content of OC, more stable aggregates and low content of residual water. On F1 long, gOC, gA₂, gKS and gθ_r weigh more and positively,_, whereas BD weighs negatively, as expected. Therefore, at this scale, the structural porosity and the content of OC significantly affect the flow of water in the soil. This factor can then be considered a direct “hydraulic indicator” of the properties of the soil. F2 long, although it explains only 24% of the variance at a 3000 m scale, it was nevertheless retained because it can be assumed as an indirect “aggregate stability indicator”, as TS weighs much more negatively than the other variables (Table 6).

Therefore, from the composition of the retained factors, it is evident how the relationships between the variables change depending on the scale, as they are linked to processes that occur in the soil over different scales [12]. It is, therefore, expected that the maps of the three factors also have a different pattern and, thus, produce a different partition of the field into homogeneous areas. To stress the differences among the structures of spatial dependence, three classes of equal frequency, i.e., of the same spatial extension, were used, which could be defined as follows: low, medium and high classes of total porosity, water flow and aggregate stability for F1 short, F1 long and F2 long, respectively, taking into account the particular composition of these regionalized factors (Figure 12).

The F1 short map (Figure 12a) appears to be divided into several micro-areas characterized by different porosity and consequent permeability, and aggregate stability, to be associated with local variations in texture and OC, but without a clear overall trend. The high areas are those particularly suited to mechanical tillage and irrigation (especially the wide blue coloured area in the south).

The geometric pattern changes markedly in the map of F1 long, in which one can recognize a large central zone, which also extends partly to the south in a narrower form, characterized by good water flow conditions favored by high structural porosity. These properties, determined by pedological characteristics and higher OC content at a longer scale, should reduce the risk of flooding or waterlogging in the occurrence of extreme weather events.

The map of F2 long (Figure 12c) also appears rather fragmented, as shown in Figure 11a, with numerous areas characterized by higher aggregate stability, which improves resistance to the destructive action of water, extending from the center of the field towards the south (blue coloured areas). Summarizing the above, it is evident that there is no single partition of the field into homogeneous areas based on the hydraulic properties, but this depends significantly on the scale. Nevertheless, it can be said that there is a large central zone, which is more suitable for irrigation, whereas the northern and southern zones are more variable. In particular, the southern zone with a coarser texture might cause problems for irrigation, so it would be advisable to intervene with smaller but more frequent irrigation volumes and more substantial organic fertilization.

4. Conclusions

In this paper, a particular approach based on multivariate geostatistics is proposed in order to jointly process (data fusion) two types of data: soil properties by point-based measurements and functional data to determine the parameters of Dexter’s model of the SWRC curve. This model was chosen among many because its parameters have physical meaning to be associated with the hydraulic properties of the soil. This approach proved effective in interpolating functional data and in rendering scale-dependent field delineation. However, this study highlights two issues. First, the scale of sampling should be appropriate for the observed variability. An excessively high nugget, which prevails over the other structured spatial components, is indicative of the need to intensify sampling by reducing its scale.

Second, there is no single partitioning of the field into homogeneous areas to be subjected to differential management; instead, it is a function of scale. Therefore, it is critical to determine which agronomic procedure is to be optimized and to investigate the scale of the processes involved. Unfortunately, there are few scientific works that have produced a partitioning of an agricultural field according to scale, whereas this crucial topic would require a more widespread and detailed treatment in precision agriculture.