Workflow to Establish Time-Specific Zones in Precision Agriculture by Spatiotemporal Integration of Plant and Soil Sensing Data

Management zones (MZs) are used in precision agriculture to diversify agronomic management across a field. According to current common practices, MZs are often spatially static: they are developed once and used thereafter. However, the soil–plant relationship often varies over time and space, decreasing the efficiency of static MZ designs. Therefore, we propose a novel workflow for time-specific MZ delineation based on integration of plant and soil sensing data. The workflow includes four steps: (1) geospatial sensor measurements are used to describe soil spatial variability and in-season plant growth status; (2) moving-window regression modelling is used to characterize the sub-field changes of the soil–plant relationship; (3) soil information and sub-field indicator(s) of the soil–plant relationship (i.e., the local regression slope coefficient[s]) are used to delineate time-specific MZs using fuzzy cluster analysis; and (4) MZ delineation is evaluated and interpreted. We illustrate the workflow with an idealized, yet realistic, example using synthetic data and with an experimental example from a 21-ha maize field in Italy using two years of maize growth, soil apparent electrical conductivity and normalized difference vegetation index (NDVI) data. In both examples, the MZs were characterized by unique combinations of soil properties and soil–plant relationships. The proposed approach provides an opportunity to address the spatiotemporal nature of changes in crop genetics × environment ×management interactions.


Introduction
Crop yields and resource use efficiency (e.g., nutrients and water) have a strong spatial component, which can be observed over a wide range of scales, from regional to subfield [1].At the field scale, yield variability in uniformly managed fields is often related to the spatial variability of soil properties and their impact on plant growth [2][3][4][5][6].This variability can be addressed using precision agriculture practices [7], such as variable rate management (VRM) [8].According to VRM principles, efficiency or crop production can be increased by varying agronomic inputs over a field according to varying soil and crop conditions [9].Within a field, areas with similar soil properties (e.g., texture, topography, water holding capacity) that can be managed uniformly are commonly called management zones (MZs) [10].To justify the subdivision of a field into MZs, there should be sizeable difference in soil properties between MZs [9,11].The use of MZs for VRM has been shown to increase productivity, decrease costs, and/or reduce environmental impacts of agronomic practices [12,13].Several authors have shown that within-MZ management should change over time [14][15][16][17][18][19][20], such an approach is often referred to as "dynamic VRM" [14].
According to current practices, MZ-designs derived from soil maps, crop sensing, and/or historical yield maps are generally designed once and then implemented year after year.Nevertheless, several authors have reported inconsistent benefits for this precision management strategy [21,22].Spatial patterns in yield tend to change from year-to-year, mostly because of changes in meteorological conditions [12,[36][37][38].In other words, the spatial patterns of most soil properties are fairly stable in time, but at different times plants may be limited in different ways at the same location because of the influence of transient factors affecting the soil-plant relationship, such as meteorological factors and agronomical management [12,38,39].
Sadler et al. [40] indicated that there is a need for accurate and inexpensive systems to delineate dynamic management zones, obtained by sensing within-field variability in real time, so that agricultural management can be controlled adaptively.Recent research strongly suggested that MZ designs should change over time, both intra-and inter-seasonally [19,20,38,[41][42][43][44].Myers [45] formally justified the need for time-specific and spatially dynamic VRM through the "fundamental theorem of precision agriculture production" where Crop Yield is a function of Genetics × Environment × Management × Space × Time interactions.According to this theorem, the spatiotemporal variability of crop performance should be addressed by adjusting the agronomic prescriptions over time and space.Several examples of protocols and data analysis workflow for static soil and/or plant-based MZ delineation are present in the scientific literature [46][47][48][49].To our best knowledge, time-specific MZ delineation based on soil and in-season plant information has not been commonly discussed in the literature.Particularly, there is a lack of protocols and analytical workflow that farm managers, agricultural consultants, and scientists can use to take advantage of free/inexpensive in-season crop information (e.g., from UAV, Sentinel 2 satellite) and high-resolution soil maps.
We aim to present a novel workflow for the selection of time-specific MZs according to in-season spatial measurements of crop growth status and its relationship with high-resolution soil spatial information.The MZ should identify areas with homogeneous and unique (within a single field) soil-plant relationships.The MZ-delineation workflow will be described in detail.We will also provide two examples on how to implement the workflow.The first example is based on synthetic data.The second example uses data from a maize (Zea mays L.) field in northeastern Italy.

Time-Specific MZ-Delineation Workflow
The time-specific MZ delineation can be implemented as follows: • STEP 1. Soil and time-specific plant spatial information acquisition, pre-processing, interpretation, and interpolation In STEP 1, high-resolution spatial measurements for target soil properties and in-season plant-canopy information, such as canopy reflectance, are acquired, pre-processed, interpreted, and interpolated.In-season plant-canopy reflectance is acquired and used as an indication of crop status.Soil spatial information is used to interpret the crop canopy measurements.
Soil sensor data acquisition should be carried out according to established protocols [50][51][52] to increase the accuracy and consistency of the survey across large areas.Attention should be paid to selecting those sensors which represent the spatial variability of soil properties known or believed to influence a crop at the site of interest [53,54].For information on sensor data pre-processing (e.g., conversion of spatial coordinates, removals of outliers), readers are referred to the first protocol step of Córdoba, Bruno, Costa, Peralta, and Balzarini [49].Subsequent to the sensor surveys, soil sampling should be carried out across the root zone (e.g., 0-1.2 m) to calibrate/interpret the sensor readings [50].Sensor-directed sampling schemes can be used to minimize the number of sampling sites [55,56].Exploratory analyses, such as correlation analysis and principal component analysis (PCA) should be carried out to investigate the relationships between soil sensor and laboratory soil analyses.The strength of the relationships between collocated soil sensor values and laboratory soil analyses should be investigated.If these relationships are moderately to very strong, the sensor data can be considered as an indicator of spatial variability of the target soil properties.Then, sensor data can be used to generate maps of the selected soil properties [57,58].Soil maps should be generated only if acceptable prediction errors [59] are obtained.Alternatively, for weak to moderately strong relationships, the soil sensor maps should be used as qualitative indicator of soil spatial variability.Córdoba, Bruno, Costa, Peralta, and Balzarini [49] describe how to practically process interpolated data to obtain a raster of desired block support (e.g., of the same resolution chosen for the MZ design) in the second step of their MZ-delineation protocol.
Free (e.g., Sentinel 2 satellite) or inexpensive (e.g., from UAV) crop measurements are available throughout the growing season with moderately-high temporal resolution.For example, the Sentinel 2 satellite from the European Space Agency provides multi-spectral canopy reflectance at the 10 × 10-m resolution, with a 5-day revisit time, free of charge.Remote sensing of crop canopy data can be used to calculate vegetation indices [34].Point measurements, such as those from tractor-mounted active spectrometers [30], should be pre-processed and interpolated similarly to soil sensing measurements.High-resolution raster data, such as that from satellite imagery, may need to be re-gridded and scaled to the selected block support.

Time-Specific Sub-Field Soil-Plant Modeling
In STEP 2, soil information from STEP 1 is used to interpret in-season (i.e., time-specific) measurements of crop status.The interpretation of crop canopy sensing measurements is not straightforward-as they are influenced by species × growth-stage × stress levels × soil background interactions [60][61][62].
Moving window spatial regression modeling, such as geographically-weighted regression (GWR) [63,64], can be used to understand the local (i.e., sub-field scale) variation in the plant-soil relationship.With GWR, a regression is run for each grid location, rather than for the whole study area [65].Soil map(s) are used as the independent (explanatory) variable(s) and the in-season crop status maps are the dependent variable.When multiple explanatory variables are available, one should consider standardizing them.The GWR allows non-stationarity of the regression equation parameters a (e.g., intercept, slope), estimating their values at each location i.For a dependent variable y the equation reads: where ε, is a random error term, a 0 is the regression intercept, and a 1 to a k are the regression coefficient(s) for each explanatory variable(s) x 1 to x k .The spatial variability of the soil-plant relationship can be described with the maps of local a 1 to a k coefficients (i.e., the regression slope[s]).Maps of GWR slope coefficient(s) show the spatial variability of the impact (i.e., sensitivity) of the explanatory variable on the regression [65].
In the GWR framework, spatial weighting is determined by incorporating all the dependent and explanatory variables falling within a geographical kernel of each target feature [65].The values of the regression parameters and goodness-of-fit of the GWR depend on how the kernel size is chosen [66].Maps of the estimated dependent variable, local coefficient of determination (R 2 ), and local Pearson correlation coefficient r can be generated with the GWR [65].The GWR is available in commercial GIS software platforms (e.g., ArcMap's Spatial Statistics package [65], version 10.5.1;ESRI, Redlands, CA, USA), and freeware (e.g., spgwr package in R, version 3.4.1; the R Foundation for Statistical Computing Platform, Vienna, Austria).

Time-Specific MZ Delineation with Cluster Analysis
In STEP 3, the MZ delineation is carried out using the soil map(s) and the time-specific GWR slope maps from STEP 2 as ancillary variables.As indicated by Córdoba, Bruno, Costa, Peralta, and Balzarini [49], fuzzy c-means (also known as "k-means") unsupervised clustering algorithms [67] can be employed to classify the data into MZs.ArcMap's Grouping Analysis tool (e.g., [68]), the Management Zone Analyst software [69] or the EZZone online tool [70] can be used to delineate MZs.Several MZ designs can be tested.The optimum number of MZs can be identified using cluster validity functions, including: the Calinski-Harabasz criterion [71], the fuzziness performance index [67], the normalized classification entropy index [67,72], and the Jenks optimization method [73].The Calinski-Harabasz criterion (CHC), also known as pseudo F-statistic, describes the ratio between within-MZ similarity and between-MZ differences.It is defined as: where N is the number of data points, MZn is the number of considered MZs, BMZSS is the between-MZ sum of squares, and WMZSS is the within group sum of squares.The larger the value of the CHC the higher are the within-MZ homogeneity and between-MZ differences.Finally, one may consider smoothing the fuzzy c-means clustering results to reduce zone fragmentation [49,70].

Time-Specific MZ-Design Quality Control and Interpretation
In STEP 4, the quality of the MZ design from STEP 3 should be checked.Each MZ should identify a unique combination of soil and soil-plant relationship characteristics.To infer differences across MZs, parametric analysis of variance or the Kruskal-Wallis (KW) rank test [74] can be used.The KW test is a nonparametric analysis assessing if samples originate from the same distribution.The test can be used as alternative to the standard analysis of variance when assumptions for parametric testing are not met.

Synthetic Data Example
The proposed workflow can be demonstrated using synthetic data for soil clay content (percentage) and normalized difference vegetation index (NDVI) [75].
A spatial random field for clay content was generated using the RFsimulate function in the RandomFields package in R. The function estimated a Gaussian random field (RPgauss function) having an exponential covariance function (RMexp function with variance = 25 and scale = 300), a nugget effect (RMnugget function with variance = 5), and a pure trend model with covariance 0 (RMtrend function with mean = 15).The simulations were carried out over a square of size 250 × 250 cells.The random field was then converted into a 250 × 250-pixel raster with pixel size = 1 m using the function raster from the raster package in R. The clay content raster was then exported from R to a text file using the writeRaster function.We refer to this raster as the true clay content (TCC) map.
Next, we mimicked a field procedure for generating a soil clay content map, using TCC as the true, unknown, high resolution clay content.Accurate approximations of true soil properties can be obtained when high-resolution covariates, such as those from proximal-soil sensing [76], are available [77].Soil sensor surveys can be calibrated to estimate the spatial variability of a target soil property by using laboratory measurements on collocated soil cores [58].Heggemann et al. [78] reported that gamma-ray sensor readings can be calibrated to predict texture values with high accuracy (up to ~95 percent of observed variance in soil texture).We mimicked a typical sensor survey (e.g., Lesch [55]) in which 1870 data points (average nearest neighbor distance = 3.2 m) were spread across 24 nearly-parallel transects (average nearest neighbor distance = 9.9 m).Values from the TCC map were extracted at the sensor survey locations.A random error having mean = 0 and variance equal to 5 percent of the extracted TCC values was added to the sensor data.This simulated a realistic sensor calibration with goodness-of-fit with a R 2 close to 0.95.The spatial autocorrelation of the calibrated sensor measurements was described with a spherical semivariogram having range = 50.3m, and nugget equal to 67% of the total sill.This spatial structure was similar to those reported by other authors in a 5-ha clay loam field in Italy [79] and in a cluster of fields with contrasting soil properties (from clay to gravelly) in Germany [80].These calibrated sensor measurements were then interpolated using simple kriging in with ArcMap's Geostatistical Analyst package.The resulting sensor-derived clay content (SCC) map was retained for further analyses.
NDVI is a vegetation index calculated from surface reflectance in the red (RED) and near-infrared (NIR) regions of the electromagnetic spectrum.It is defined as: NDVI ranges from −1 to +1.NDVI for agricultural crops usually ranges from ~0.1 to ~0.9, with lush vegetation generally having high NDVI [81,82].The NDVI information was simulated with the understanding that, often, the spatial relationship between remote sensing canopy measurement and collocated soil properties has both deterministic and spatial random components [83][84][85].The deterministic component of the relationship can be a linear model between the soil property and the available remote sensing plant information [84].The spatial random component is often equal to the field of spatially correlated residuals from the deterministic linear model [57,83].We simulated the NDVI as follow: where SpERR was a spatial error raster, S was a scaling factor = 0.01, and O was an offset coefficient = 0.4.SpERR was generated in R, using the RandomFields package.We estimated a Gaussian random field based on a model with exponential covariance function (variance = 7.5 and scale = 50) and a nugget effect (variance = 5).
A GWR with 30-m bandwidth size was used to model the spatial variability of the synthetic NDVI image using the SCC map as explanatory variable.Then, MZs were delineated with a fuzzy c-mean unsupervised cluster analysis.Clustering was carried out in ArcMap using the Grouping Analysis tool from the Spatial Statistics package.The SCC and the time-specific GWR slope maps were used for the cluster analysis.No spatial constraints were set for the cluster analysis: the covariate data points could be grouped according to their value, disregarding the values of their geographical coordinates (i.e., points did not need to be neighbors to be part of the same MZ).The classifications were carried out for scenarios outputting two, three, and four MZs.Then, the optimal number of MZs was identified using the CHC.Differences in clay (from the SCC map), NDVI, and GWR slope across MZs were investigated with the KW test using STATISTICA (version 12, StatSoft Inc., Tulsa, OK, USA).

Site Description
Another example is given using experimental data from researchers at the University of Padua, Italy [86][87][88].Here, we describe two time-specific MZ designs for a maize field based on remote sensing imagery acquired during the late kernel blister stage (R-2), in two consecutive growing seasons.
The study site (Figure 1) is a 21-ha field located at Chioggia, Venice, Italy, along the southern margin of the Venice Lagoon (45 • 10 7.0 N, 12 • 13 55.0 E).Previous studies [86][87][88][89] discussed the relationship between soil and maize yield at this site.The site lies below average sea level, and the groundwater level is kept fairly shallow (approximately between −0.6 to −1.5 m below ground level) by a pumping station [87] to promote sub-irrigation [88].The soil in the area is classified as Molli-Gleyic Cambisol [90] with two coarsely-textured paleochannels crossing it with SW-NE direction (Figure 1).The field is affected by soil salinity due to its proximity to the Venice Lagoon [86].
Here, we analyzed the two maize growing seasons (April to September), 2010 and 2011.The two growing seasons differed greatly in terms of meteorology.Compared to the average April to September rainfall from 1993 to 2012, 2010 was rainy (in the upper third quartile, 540 mm) and 2011 was a drought year (in the first quartile, 200 mm) [88].The daily average reference evapotranspiration was 4.01 mm day −1 in 2010 and 4.42 mm day −1 in 2011.For the entire growing season, the reference evapotranspiration was 569 mm in 2010 and 672 mm in 2011 [88].Agronomic management (base-dressing of 64 kg N ha −1 , and 94 kg P 2 O 5 ha −1 and urea top-dressing of 184 kg N ha −1 ) and maize hybrid (PR32P26, Pioneer Hi-Bred Italia, Gadesco-Pieve Delmona, Italy) were the same in the two years [86].
was identified using the CHC.Differences in clay (from the SCC map), NDVI, and GWR slope across MZs were investigated with the KW test using STATISTICA (version 12, StatSoft Inc., Tulsa, OK, USA).

Site Description
Another example is given using experimental data from researchers at the University of Padua, Italy [86][87][88].Here, we describe two time-specific MZ designs for a maize field based on remote sensing imagery acquired during the late kernel blister stage (R-2), in two consecutive growing seasons.
The study site (Figure 1) is a 21-ha field located at Chioggia, Venice, Italy, along the southern margin of the Venice Lagoon (45°10′7.0′′N, 12°13′55.0′′E).Previous studies [86][87][88][89] discussed the relationship between soil and maize yield at this site.The site lies below average sea level, and the groundwater level is kept fairly shallow (approximately between −0.6 to −1.5 m below ground level) by a pumping station [87] to promote sub-irrigation [88].The soil in the area is classified as Molli-Gleyic Cambisol [90] with two coarsely-textured paleochannels crossing it with SW-NE direction (Figure 1).The field is affected by soil salinity due to its proximity to the Venice Lagoon [86].
Here, we analyzed the two maize growing seasons (April to September), 2010 and 2011.The two growing seasons differed greatly in terms of meteorology.Compared to the average April to September rainfall from 1993 to 2012, 2010 was rainy (in the upper third quartile, 540 mm) and 2011 was a drought year (in the first quartile, 200 mm) [88].The daily average reference evapotranspiration was 4.01 mm day −1 in 2010 and 4.42 mm day −1 in 2011.For the entire growing season, the reference evapotranspiration was 569 mm in 2010 and 672 mm in 2011 [88].Agronomic management (basedressing of 64 kg N ha −1 , and 94 kg P2O5 ha −1 and urea top-dressing of 184 kg N ha −1 ) and maize hybrid (PR32P26, Pioneer Hi-Bred Italia, Gadesco-Pieve Delmona, Italy) were the same in the two years [86].

Soil Maps and Plant Spatiotemporal Information
Geospatial measurements of soil apparent electrical conductivity (EC a ) were used as an indicator of soil spatial variability at the study site.EC a over the 0-1.5 m soil profile (EC a Deep) was measured in April of 2011 across the site with a frequency-domain electromagnetic induction sensor (CMD-1, GF Instruments, Brno, Czech Republic) at 20,470 geo-referenced locations.Soil samples were taken at 91 locations in May 2010 (Figure 1) down to 1.2 m with 0.3-m increments.Here we discuss salinity and texture for the 0-1.2-m soil profile.Texture was measured with a Mastersizer 2000 (Malvern Instruments Ltd., Great Malvern, UK) and salinity was measured as EC 1:2 (i.e., the electrical conductivity of a 1:2 soil-water extract) [91].Principal component analysis (PCA) carried out with STATISTICA was used to visually assess similarities and differences between the EC a Deep dataset and the selected soil properties [92].The point measurements were preprocessed with a procedure comparable to that of Córdoba, Bruno, Costa, Peralta, and Balzarini [49] and spatially interpolated with ordinary kriging following Scudiero, Teatini, Corwin, Deiana, Berti, and Morari [86] on a 10 × 10-m block support-the desired support for MZ at the site.
Two remotely sensed measures of crop status were available.On 31 July 31 2010 (10:39:40 UTC) and 9 July 2011 (10:47:30 UTC) WorldView-2 (DigitalGlobe, Westminster, CO, USA) satellite scenes were acquired over the study area.For both years, acquisition dates corresponded to late R-2 [88].Pre-processing procedures, including sensor calibration, atmospheric correction, and radiometric normalization, were applied according to Vicente-Serrano et al. [93].Radiometric calibration was required to convert digital numbers to top-of-atmosphere radiances [W m −2 sr −1 µm −1 ], using the absolute radiometric calibration factors and effective bandwidths for each band, according to the satellite data provider [94].No topographic correction was applied to the images because the study area is relatively flat and the solar incident angles (28.11 • in 2010; 23.58 • in 2011) were quite similar at the two acquisition times.The top-of-atmosphere radiance was then transformed to surface reflectance through the 6S code [95].The values of aerosol optical thickness (AOT) at 550 nm collected from a nearby AERONET (Aerosol Robotic Network, https://aeronet.gsfc.nasa.gov/)[96] station (45 • 18 50.0"N, 12 • 30 29.9" E) were used as input for 6S.The AOT measurements were obtained simultaneously to the satellite overpasses.WorldView-2 reflectance has 2 × 2-m spatial resolution over eight spectral bands, at wavelengths spanning from 400 to 1040 nm.The red band (RED, 630-690 nm) and the near-infrared band at 770-895 nm where used to calculate NDVI according to Equation (3).The NDVI maps were aggregated (i.e., the coarsened pixel is the average of the pixels within the aggregated cell) to the 10 × 10-m cell size, as suggested by Córdoba, Bruno, Costa, Peralta, and Balzarini [49].

Time-Specific Spatial Soil-Plant Modeling
The time-specific relationship between soil properties and in-season NDVI were described using geographically weighted regressions (GWRs).The GWRs were carried out using the NDVI maps as the dependent variable and spatial soil information (i.e., EC a Deep) as the independent variable.GWRs were carried out with ArcMap using an adaptive kernel (i.e., moving window) of 70 neighbors.For locations with significant GWR, the observed-estimated NDVI relationship was used as indication of goodness-of-fit of the GWR models.The GWR slope map was selected as indicator of soil-plant relationship type.

Delineating MZs with Cluster Analysis
Clustering was carried out using the Grouping Analysis tool in ArcMap.Spatial soil information (i.e., EC a Deep) and the time-specific GWR slope map were used for the cluster analysis.No spatial constraints were set for the cluster analysis.The classifications were carried out for scenarios outputting three, four, and five MZs.Then, the optimal number of MZs was identified using the CHC.To compare the proposed time-specific soil-plant based MZ delineation with a more traditional, static MZ design, EC a Deep alone was used to select a static MZ design.

Evaluation and Interpretation of MZ Design
Differences across MZs in the two years were investigated using the KW test in STATISTICA.The following variables were tested for differences across MZs: EC a Deep, GWR slope, GWR local r, EC 1:2 , and sand and clay contents.Additionally, in-season potential yield estimations were used to evaluate the MZ designs.Appendix A reports on how yield maps (obtained from Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88]) and NDVI measurements were used to calculate yield prediction maps.

Synthetic Data Example
Every time the crop status is mapped during the season (e.g., every week), a new set of MZ is delineated (i.e., time-specific MZ design).Clearly, this applies only for crop growth stages when crop canopy plays a relevant role in determining reflectance readings.For soil tillage, sowing, and early vegetation stages when surface reflectance is mostly determined by soil, a static MZ delineation approach based on the spatial variability of soil properties (e.g., [9]) may be more adequate.
Figure 2 outlines the proposed analytical protocol to delineate time-specific MZs from soil information and in-season crop information using the synthetic dataset.In Step 1 of the protocol, the soil and plant information were acquired and pre-processed.The clay content (percentage) map had mean = 18.6%, standard deviation = 4.5%, minimum = 9%, and maximum = 30.1%.The NDVI map had mean = 0.58, standard deviation = 0.04, minimum = 0.45, and maximum = 0.69.The NDVI values are typical of vegetative stages for crops such as maize [97] and soybean (Glycine max (L.) Merr.) [98].A simple linear model between the two maps had R 2 = 0.61, with an intercept of 0.44 and a positive slope of 0.007.Similar goodness-of-fit values between soil properties and vegetation reflectance were reported by Gomez et al. [99] for clay and by Gomez et al. [100] for soil organic carbon.In Step 2, the clay map was used as the explanatory variable in a GWR with NDVI as dependent variable.Through the GWR analysis, clay explained 94.5% of the observed variance of NDVI.Such high goodness-of-fit can be found in real-world data: Scudiero, Corwin, Wienhold, Bosley, Shanahan, and Johnson [54] observed R 2 ranging between 0.83 and 0.94 for their GWR analyses between soil maps and remotely sensed winter wheat (Triticum aestivum L.) canopy reflectance.The GWR slope had mean = 0.009, standard deviation = 0.008, minimum = −0.018,and maximum = 0.038.In Step 3, the clay and the GWR slope maps were used to delineate MZs via unsupervised fuzzy c-means clustering.The CHC indicated that the best number of MZs was three.
In Step 4, the MZ-design was evaluated and interpreted.The three MZs identified unique combinations of clay content and GWR slope values.MZ 1 was characterized by low NDVI, low clay content, and large GWR slope values (indicating high sensitivity of NDVI to clay).MZ 2 was characterized by high NDVI, high clay, and moderately low GWR slope values.MZ 3 was characterized by low NDVI, low clay content, and the smallest GWR slope values.

Soil and Plant Information
Figure 3a shows an example of short-scale soil spatial variation at the site.Figure 3b depicts the frequency distribution of measured soil salinity.At this site, soils with EC 1:2 >0.79 dS m −1 should be considered salt affected, with EC 1:2 >1.4 dS m −1 indicating strongly saline soils [88].Figure 3c reports the observed texture at the 91 soil sampling locations.Soil texture is predominantly in the loam and sandy loam classes.Few samples are in the clay loam and loamy sand class.

Soil and Plant Information
Figure 3a shows an example of short-scale soil spatial variation at the site.Figure 3b depicts the frequency distribution of measured soil salinity.At this site, soils with EC1:2 >0.79 dS m −1 should be considered salt affected, with EC1:2 >1.4 dS m −1 indicating strongly saline soils [88].Figure 3c   The ECa Deep map is shown in Figure 4a.Descriptive statistics for the map are reported in Table 1.ECa Deep showed a Pearson correlation coefficient of 0.89 with salinity, −0.51 with sand, and 0.39 with clay.The PCA extracted two factors (with eigenvalues of 2.61 and 0.89) explaining 65.3% (first component) and 22.3% (second component) of the total variance in the dataset.The bi-plot in Figure 3d-see Abdi and Williams [92] for interpretation guidelines-indicates that the first component contrasts the positive contribution of clay content and salinity on the ECa Deep measurements with the negative correlation between sand content and ECa Deep.In Figure 3d, clay content was not clustered with salinity and ECa Deep, indicating that the three variables were not collinear, yet negatively correlated with sand content.According to the correlation and PCA analyses, high values of ECa Deep were mainly interpreted as an indication of high salinity and fine soil texture.ECa readings larger than 1-2 dS m −1 are most likely due to high soil salinity rather than other edaphic factors contributing to soil conductivity [101,102].Conversely, lower ECa Deep was interpreted as an The EC a Deep map is shown in Figure 4a.Descriptive statistics for the map are reported in Table 1.EC a Deep showed a Pearson correlation coefficient of 0.89 with salinity, −0.51 with sand, and 0.39 with clay.The PCA extracted two factors (with eigenvalues of 2.61 and 0.89) explaining 65.3% (first component) and 22.3% (second component) of the total variance in the dataset.The bi-plot in Figure 3d-see Abdi and Williams [92] for interpretation guidelines-indicates that the first component contrasts the positive contribution of clay content and salinity on the EC a Deep measurements with the negative correlation between sand content and EC a Deep.In Figure 3d, clay content was not clustered with salinity and EC a Deep, indicating that the three variables were not collinear, yet negatively correlated with sand content.According to the correlation and PCA analyses, high values of EC a Deep were mainly interpreted as an indication of high salinity and fine soil texture.EC a readings larger than 1-2 dS m −1 are most likely due to high soil salinity rather than other edaphic factors contributing to soil conductivity [101,102].Conversely, lower EC a Deep was interpreted as an indication of coarser soils with low salinity.Note that Scudiero, Teatini, Corwin, Deiana, Berti, and Morari [86] reported that the spatial patterns of EC a Deep remained stable over time at the study area.
The two in-season NDVI maps for the late R-2 growth stage are reported in Figure 4b (2010) and Figure 4e (2011).Table 1 reports the descriptive statistics for the two NDVI maps.The 2010 NDVI was more heterogeneous (mean = 0.79, standard deviation = 0.046) than that of 2011 (mean = 0.83, standard deviation = 0.018).During the R-2 growth stage, environmental stressors (e.g., drought, heat, severe nutrient deficiency) may prevent kernels from developing properly (e.g., kernel abortion at the ear tip) [103,104].indication of coarser soils with low salinity.Note that Scudiero, Teatini, Corwin, Deiana, Berti, and Morari [86] reported that the spatial patterns of ECa Deep remained stable over time at the study area.
The two in-season NDVI maps for the late R-2 growth stage are reported in Figure 4b (2010) and Figure 4e (2011).Table 1 reports the descriptive statistics for the two NDVI maps.The 2010 NDVI was more heterogeneous (mean = 0.79, standard deviation = 0.046) than that of 2011 (mean = 0.83, standard deviation = 0.018).During the R-2 growth stage, environmental stressors (e.g., drought, heat, severe nutrient deficiency) may prevent kernels from developing properly (e.g., kernel abortion at the ear tip) [103,104].4 summarizes the GWR analysis between soil ECa Deep and in-season NDVI for 2010 and 2011.For 2010, 1149 cells (62.5% of total) were characterized by a significant (p < 0.05) local Pearson correlation coefficient r.For 2011, 56.8% of the cells had a significant local r.At the locations with significant local r, the observed-estimated NDVI relationship was characterized by R 2 = 0.73 in 2010 (Figure 3d) and R 2 = 0.62 in 2011 (Figure 3g).There was a strong negative relationship in both years (r = −0.36 in 2010 and r = −0.43 in 2011) between the significant local r and ECa Deep maps.The slope of this relationship was steeper in 2011 (−0.46 with standard error = 0.027) than in 2010 (−0.38 with standard error = 0.026).This may indicate that low ECa Deep (e.g., indicating coarse texture) was a   4 summarizes the GWR analysis between soil EC a Deep and in-season NDVI for 2010 and 2011.For 2010, 1149 cells (62.5% of total) were characterized by a significant (p < 0.05) local Pearson correlation coefficient r.For 2011, 56.8% of the cells had a significant local r.At the locations with significant local r, the observed-estimated NDVI relationship was characterized by R 2 = 0.73 in 2010 (Figure 3d) and R 2 = 0.62 in 2011 (Figure 3g).There was a strong negative relationship in both years (r = −0.36 in 2010 and r = −0.43 in 2011) between the significant local r and EC a Deep maps.The slope of this relationship was steeper in 2011 (−0.46 with standard error = 0.027) than in 2010 (−0.38 with standard error = 0.026).This may indicate that low EC a Deep (e.g., indicating coarse texture) was a greater constrain in 2011 than 2010.Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] showed that water stress in areas with coarser soil texture is particularly limiting in dry years, such as 2011.
In areas where the significant local r values were consistently negative in both years, the EC a Deep averaged 1.17 dS m −1 (standard deviation = 0.39 dS m −1 ), whereas, in areas with negative local r only in 2010, it averaged 0.94 dS m −1 (standard deviation = 0.28 dS m −1 ).This suggested that only very high EC a Deep, which indicates high salinity, consistently limited crop growth in the two years.As discussed by other authors [105][106][107], the effects of high salinity on crops are generally stable throughout different growing seasons.
The GWR slope maps (Figure 4c for 2010, Figure 4f for 2011) changed remarkably between the two years, indicating a change in the sensitivity of plant NDVI to in soil properties.Changes of crop growth and yield spatial patterns are expected between seasons having widely different meteorology, as discussed by Maestrini and Basso [38] and McBratney, Whelan, Ancev, and Bouma [37].

Time-Specific MZ Delineation
The GWR slope and EC a Deep maps were used to drive the delineation of in-season site-specific MZs.The CHC index indicated that a five-MZ configuration was optimal for 2010.The CHC was 1901.Overall, 45.1% of the cells were classified in the same MZ at both times.The consistency rate was 94.3% for MZ I, 75.6% for MZ II, 20.1% for MZ III, and 48.5% for MZ IV.Conversely, changes in local GWR slope led to 54.9% of the cells being assigned to a different time-specific MZ.As reported by Maestrini and Basso [38], it is reasonable to expect portions of a field to have stable (e.g., consistently high yields) and unstable crop outputs.Management in the unstable areas should be addressed within each growing season to meet desired agronomic output goals [18,38].The proposed methodology helps characterizing time-specific changes in the crop output according to soil spatial variability.greater constrain in 2011 than 2010.Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] showed that water stress in areas with coarser soil texture is particularly limiting in dry years, such as 2011.
In areas where the significant local r values were consistently negative in both years, the ECa Deep averaged 1.17 dS m −1 (standard deviation = 0.39 dS m −1 ), whereas, in areas with negative local r only in 2010, it averaged 0.94 dS m −1 (standard deviation = 0.28 dS m −1 ).This suggested that only very high ECa Deep, which indicates high salinity, consistently limited crop growth in the two years.As discussed by other authors [105][106][107], the effects of high salinity on crops are generally stable throughout different growing seasons.
The GWR slope maps (Figure 4c for 2010, Figure 4f for 2011) changed remarkably between the two years, indicating a change in the sensitivity of plant NDVI to changes in soil properties.Changes of crop growth and yield spatial patterns are expected between seasons having widely different meteorology, as discussed by Maestrini and Basso [38] and McBratney, Whelan, Ancev, and Bouma [37].Overall, 45.1% of the cells were classified in the same MZ at both times.The consistency rate was 94.3% for MZ I, 75.6% for MZ II, 20.1% for MZ III, and 48.5% for MZ IV.Conversely, changes in local GWR slope led to 54.9% of the cells being assigned to a different time-specific MZ.As reported by Maestrini and Basso [38], it is reasonable to expect portions of a field to have stable (e.g., consistently high yields) and unstable crop outputs.Management in the unstable areas should be addressed within each growing season to meet desired agronomic output goals [18,38].The proposed methodology helps characterizing time-specific changes in the crop output according to soil spatial variability.2).Note that the local GWR r and slope, which can be used to interpret the ECa Deep-NDVI relationship, remained fairly consistent over time at the different MZs.

Time-Specific MZ-Design Quality Control and Interpretation
In both years, each time-specific MZ was characterized by unique combinations of the variables used in MZ delineation.Moreover, the time-specific MZs differed greatly in terms of EC a Deep, significant GWR r and slope, salinity, and texture (Table 2).Note that the local GWR r and slope, which can be used to interpret the EC a Deep-NDVI relationship, remained fairly consistent over time at the different MZs.In both years, MZ I consisted of soils with low EC a Deep, low salinity, and coarse soil texture.MZ I had the highest positive GWR slope values, meaning NDVI would decrease with increasing sand content.Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] monitored crop water and salt stress at five soil-plant-water monitoring stations at the site.One of their stations ("Station E") was located in our MZ I. Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88] indicated that at Station E maize was not stressed by salinity but was under water stress, particularly in 2011, when water stress was described as "severe".
In contrast, MZ II in both years had GWR slope and r distributions overlapping with 0, indicating little to no influence of soil on crop NDVI variability.This is perhaps the reason why MZ II was consistently characterized by the highest potential yield estimations (Appendix A).The very saline MZ III (highest measured average EC a Deep and EC 1:2 ) was characterized by lower potential yield estimations for both years in comparison with the other MZs (Table 2), as expected for (unmanaged) very saline portions of farmlands [106].Similarly, the moderately saline MZ IV (second highest average EC a Deep), was characterized by the second lowest yield predictions in 2010 and 2011, together with MZ I. MZ V was only selected in 2010.It grouped soils with moderately high EC a Deep (third highest MZ average) together with the strongest negative GWR slope values and low potential yield estimations.In 2011, 89.5% of the 2010 MZ V locations were classified into MZ II.Previous research focusing on spatiotemporal variability of yield patterns [38], indicates that fields generally have areas of stable and unstable yield patterns.For unstable areas, in-season NDVI is suggested as the best predictor for yield spatial distribution [38].Our proposed approach further refines the NDVI information by interpreting its spatial variability as a function of soil spatial variability.
A static MZ design based on EC a Deep only (Figure 5c) would identify four areas with significantly different soil properties (Table 3), but very heterogeneous and inconsistent in terms of GWR local r and slope (Table 4).

Conclusions
Characterizing the temporal variability of the soil-plant relationship within the growing season and over years is very relevant to developing best crop management practices with VRM.Soil influences plants according to complex interactions between factors that may change in time, such as meteorological conditions, nutrient availability, and water content.Static MZs are not ideal when spatial patterns of soil-plant relationship change in time-because of changing meteorology and/or other transient factors.We presented a time-specific MZ-delineation workflow to address such spatiotemporal variability.
The proposed approach for time-specific MZs should be further tested and evaluated in the field.Future work should focus on multi-field and multi-year comparisons of time-specific and static MZs and their relative impact on crop yields and resource use efficiency (e.g., water, nutrients).

Figure 1 .
Figure 1.(a) Geographical location of the study area relatively to the Venice Lagoon, Italy, and (b) locations of soil sampling locations and coarsely textured paleochannels.Modified after Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88].

Figure 1 .
Figure 1.(a) Geographical location of the study area relatively to the Venice Lagoon, Italy, and (b) locations of soil sampling locations and coarsely textured paleochannels.Modified after Scudiero, Teatini, Corwin, Dal Ferro, Simonetti, and Morari [88].

Figure 2 .
Figure 2. Diagram outlining the proposed analytical protocol to delineate time-specific management zones (MZs) from soil and (in-season) crop information.In Step 3, MZn refers to the number of considered management zones.

Figure 2 .
Figure 2. Diagram outlining the proposed analytical protocol to delineate time-specific management zones (MZs) from soil and (in-season) crop information.In Step 3, MZn refers to the number of considered management zones.
Figure3ashows an example of short-scale soil spatial variation at the site.Figure3bdepicts the frequency distribution of measured soil salinity.At this site, soils with EC1:2 >0.79 dS m −1 should be considered salt affected, with EC1:2 >1.4 dS m −1 indicating strongly saline soils[88].Figure3creports the observed texture at the 91 soil sampling locations.Soil texture is predominantly in the loam and sandy loam classes.Few samples are in the clay loam and loamy sand class.

Figure 3 .
Figure 3. (a) Photograph showing the sharp textural change at the site (light vs. dark colors); (b) histogram for measured soil salinity (EC1:2); (c) soil textural triangle; and (d) bi-plot of selected variables (soil apparent electrical conductivity [ECa Deep], salinity, clay, and sand) on the two larger factors in the principal components analysis.

Figure 3 .
Figure 3. (a) Photograph showing the sharp textural change at the site (light vs. dark colors); (b) histogram for measured soil salinity (EC 1:2 ); (c) soil textural triangle; and (d) bi-plot of selected variables (soil apparent electrical conductivity [EC a Deep], salinity, clay, and sand) on the two larger factors in the principal components analysis.
3 for three MZs, 1979.7 for four MZs, and 2002.7 for five-MZs.In 2011, the CHC index indicated that a four-MZ configuration was optimal.CHC was 1248.9 for three MZs, 1383.4 for four MZs, and 1374.0 for five MZs.The two MZ delineations are shown in Figure 5a (for 2010) and Figure 5b (for 2011).

3. 2 . 3 .
Time-Specific MZ Delineation The GWR slope and ECa Deep maps were used to drive the delineation of in-season site-specific MZs.The CHC index indicated that a five-MZ configuration was optimal for 2010.The CHC was 1901.3 for three MZs, 1979.7 for four MZs, and 2002.7 for five-MZs.In 2011, the CHC index indicated that a four-MZ configuration was optimal.CHC was 1248.9 for three MZs, 1383.4 for four MZs, and 1374.0 for five MZs.The two MZ delineations are shown in Figure 5a (for 2010) and Figure 5b (for 2011).

Figure 5 .
Figure 5. In-season delineation of time-specific management zones (MZs) for maize late kernel blister stage of (a) 2010 and (b) 2011.(c) Static MZ designs using soil apparent electrical conductivity for the 0-1.5 m soil profile (ECa Deep) only.

3. 2 . 4 .
Time-Specific MZ-Design Quality Control and Interpretation In both years, each time-specific MZ was characterized by unique combinations of the variables used in MZ delineation.Moreover, the time-specific MZs differed greatly in terms of ECa Deep, significant GWR r and slope, salinity, and texture (Table

Figure 5 .
Figure 5. In-season delineation of time-specific management zones (MZs) for maize late kernel blister stage of (a) 2010 and (b) 2011.(c) Static MZ designs using soil apparent electrical conductivity for the 0-1.5 m soil profile (EC a Deep) only.

Figure A1 .
Figure A1.(a) Scatterplot between yield and WorldView-2 NDVI for both years.The modeled inseason potential yield is depicted with a solid red line.Potential yield maps calculated for (b) 2010 and (c) 2011.

Figure A1 .
Figure A1.(a) Scatterplot between yield and WorldView-2 NDVI for both years.The modeled in-season potential yield is depicted with a solid red line.Potential yield maps calculated for (b) 2010 and (c) 2011.

Table 1 .
Descriptive statistics for the soil apparent electrical conductivity map (ECa Deep), and the 2010 and 2011 Normalized Difference Vegetation Index (NDVI) maps.

Table 1 .
Descriptive statistics for the soil apparent electrical conductivity map (EC a Deep), and the 2010 and 2011 Normalized Difference Vegetation Index (NDVI) maps.

Table 2 .
Count (n), mean, and standard deviation (Std.Dev.) of apparent electrical conductivity for the 0-1.5-m soil profile (EC a Deep), local regression slope and Pearson correlation coefficient (r) obtained from significant geographically weighted regressions (GWR), soil salinity (EC 1:2 ), sand and clay contents, and potential yield estimations at each time-specific management zone (MZ) for 2010 and 2011 (n.a.= not available; NA = not assessed).
1Different letters are significantly different between MZs at the p < 0.05 level according to the Kruskal-Wallis (KW) test.2Percentage of cells having significant GWR within the MZs.

Table 3 .
Count (n), mean, and standard deviation (Std.Dev.) of apparent electrical conductivity for the 0-1.5-m soil profile (EC a Deep), soil salinity (EC 1:2 ), and sand and clay contents for the static management zones (MZs) delineated using EC a Deep only.

Table 4 .
Mean, standard deviation (Std.Dev.), and count (n) for the correlation coefficient (r) and slope obtained from significant geographically weighted regressions (GWRs) for the static management zones (MZs) delineated using soil apparent electrical conductivity for the 0-1.5-m soil profile (EC a Deep) only.
1Percentage of cells having significant GWR within the MZ.