Soil and climate factors drive spatio-temporal variability of arable crop yields under uniform management in Northern Italy

ABTRACT Soil and weather data were used to analyse spatio-temporal yield patterns of winter cereals (wheat) and spring dicots (sunflower and coriander) in an 11-ha field in Northern Italy (44.5° N, 12.2° E), during 2010–2014. Three yield stability classes (YSCs) were established over multiple years, based on spatio-temporal characteristics across the field: high yielding and stable (HYS), low yielding and stable (LYS), and unstable. The HYS class (46% of field area) staged a 122% relative yield with low temporal variability. The unstable class (24% of field area) was slightly more productive (83% relative yield) than the LYS class (30% of field area, and 80% relative yield), but less consistent over time. Crop yields evidenced negative correlations with sand content; positive correlations with silt and clay content. Soil properties were quite consistently classified among YSCs: the LYS and unstable classes were associated with higher sand content and lower cation exchange capacity, suggesting that these characteristics lead to fluctuation and depression of final yield. Establishing YSCs based on spatio-temporal yield appears a sound approach to appraise field potential. It results in strategic and tactical decisions to be taken, depending on the profile of spatial and temporal productivity in different field areas.


Introduction
Precision agriculture (PA) has a great potential to increase crop growth and final yield through the application of variable crop inputs (Basso et al. 2017). Specific crop inputs at the right time and place is highly encouraged in today's agriculture. Therefore, the focus of this study is to analyse spatial and temporal variability of field crops, in-season climate conditions, and soil nutrient status for optimizing crop productivity and sustaining the most efficient use of finite natural resources (Blackmore 2000;Maestrini and Basso 2018). Soil physical-chemical properties vary in space and time depending on their interaction with factors such as climate, topography and anthropogenic activities . Bullock and Bullock (2000) stressed the importance of adopting efficient methods to characterize soil spatial variability. Among them, apparent soil electrical conductivity (ECa) directed to soil sampling has been shown a rapid and reliable method to characterize field variability . However, ECa has not always staged consistent results with crop yield, due to its complex interactions with soil properties and external factors .
Climatic factors exert a strong influence on crop productivity under rainfed conditions (Iizumi and Ramankutty 2015;Asfaw et al. 2018), being responsible for consistently low yielding areas where insufficient moisture is the most limiting factor. Several studies demonstrate that the two main climatic factors, precipitation and temperature, significantly influence yield stability across growing seasons (Kukal and Irmak 2018;Maestrini and Basso 2018;Mohsenipour et al. 2018;Shiru et al. 2018). Precipitation is seen to be more impacting than the temperature on final crop yield (Kang et al. 2009) During the twenty-first century, it is expected that higher temperatures influence the regime of precipitation, and the ultimate availability of water (Mishra et al. 2014). Therefore, registering the weather course during the crop season stimulates farmers to think critically regarding crop management (Cuculeanu et al. 2002;Asfaw et al. 2018).
Furthermore, biotic and abiotic factors equally contribute to influence crop growth and development, and final yield. Among abiotic factors, low water availability and heat exert an influence on final crop yield (Mariani and Ferrante 2017), also depending on genotype adaptation to specific adversities (Zandalinas et al. 2018).
In this complex situation, many studies addressed different methods of delineating site-specific crop management (SSCM) zones within a field, by relating yield data with soil properties and external factors. Da Silva (2006) produced classified zones based on spatio-temporal yield maps. Lark and Stafford (1996) used an unsupervised fuzzy clustering method over multiple years' yield data. Swindell (1997) analyzed the spatial variability by using several crop harvest indices. Fraisse et al. (1999) combined topographical variables and ECa through unsupervised cluster analysis. Maestrini and Basso (2018) produced zones based on spatio-temporal yield of several crops with soil, crop reflectance and weather data during the growing seasons.
This research was aimed at identifying homogeneous areas for site-specific management, using soil and crop yield data. The following steps were carried out: i) establishment of spatio-temporal yield stability classes (YSCs) (Blackmore 2000;Panneton and Brouillard 2002;Blackmore et al. 2003), based on the yield data of a five-year crop rotation under uniform, rainfed management; ii) assessment of the spatial variability of soil properties determined in samples taken according to an ECa; iii) establishment of spatio-temporal YSCs based on soil properties; iv) analysis of the weather effects on temporal yield variability in the five crop seasons.

Study site description
The experimental site was an 11.07-ha field of the Agrisfera Cooperative, located near Ravenna, Italy, at N 44° 29ʹ 26", E 12° 07ʹ 44", 0 m above sea level ( Figure S1). The area falls in the Mediterranean North Environmental Zone (Metzger et al. 2005). The field was managed in a uniform rotation system with winter cereals Durum Wheat in 2010 (DW 2010) and Bread Wheat in 2012and 2014(BW 2012and 2014, and spring dicots, Sunflower in 2011 (SF 2011) and Coriander in 2013(CO 2013. Cultivation was based on the good practices for each specific crop, depending on the local conditions. The previous field history from 1976 to 2005 ( Figure S2) shows three separate parts of approximately equal length (200 m each in the north-south axis), cultivated with fruit orchard and vineyard (upper, i.e. northern, part), and arable crops (lower, i.e. southern, part). In 2006, the three fields were merged into a single arable field ( Figure S2).
Thereafter, a geostatistical analysis was performed on GY data to i) examine the degree of spatial dependence (SpD) in terms of semivariogram; ii) produce continuous grid points over the entire field before mapping; iii) combine the interpolated data intersected on the regular grid. Three main parameters describe semivariogram characteristics: i) nugget (C 0 ), the measurement error at 0 distance (h = 0); ii) sill (C 0 + C), the maximum y-axis value that increases with increasing lag distance (h), and remains constant at a higher distance; iii) range (a), the maximum distance at which data points are still correlated, i.e. the lag distance at sill value. The degree of SpD as given by Cambardella et al. (1994) explains the nugget to sill ratio (C 0 /(C 0 + C)): < 25%, indicates strong SpD; (ii) 25-75%, moderate SpD; (iii) >75%, weak SpD.
We employed the iterative cross-validation technique seeking the highest coefficient of determination (R 2 ) and minimum mean absolute error (MAE) to choose the best fitting semivariogram model among Circular, Spherical, Exponential, Gaussian, and Stable (Xiao et al. 2016). Spatial variability maps were computed by simple kriging (SK) with 10 m cell size, resulting in 24 columns and 72 rows (Moral et al. 2010;Ali et al. 2019). SK was chosen as it provides, normally, maximum R 2 and minimal error parameters (Xiao et al. 2016).
For each crop, standardized interpolated data with 1156 regular grid points were used for comparison among years, by replacing the actual GY (t/ha) with a relative GY where 100% equals field average. This allowed data from different crops to be jointly analyzed. Equation (1) is used to characterize the spatial variability maps over a single crop: where, S i =standardized yield (%) over 100% field average at point (i), y i =interpolated yield at point i (t/ha), and � y=mean interpolated yield over the entire field (t/ha). For multiple crops, a spatial variability map was produced by simply calculating the mean standardized yield, laid over the five years according to Equation (2).
Where, � S i = mean interpolated yield over 100% field average over n years, Si I = interpolated standardized yield (%) at point (i). For multiple crops, a temporal variability map was produced to assess the stability of standardized GY over the five crop years. The coefficient of variation (CV) of each grid point over the five years was calculated based on Equation (3) (Blackmore 2000).
where, CVS i = coefficient of variation of standardized yield at point (i) over n years; Si t = standardized yield (%), at point (i); S i ¼ mean standardized yield at point (i).
To define the threshold levels in spatial maps, four classes were established in both single-and multiple-year yield, based on the natural break classification method (Toshiro 2002): very low (VL), medium-low (ML), medium-high (MH) and very high (VH). Each class showed maximum difference with other classes, while the within-class variability was minimized. Likewise, four classes were defined for temporal variability map across CV ranges between 2% and 73%.

Soil sampling
The positions for soil samples were based on the procedure developed by Corwin and Lesch (2005). First, a soil ECa survey was conducted using an on-the-go sensor CMD Tiny Electromagnetic Conductivity Meter (GF Instruments, s.r.o., Brno, Czech Republic) along a 8 m transect over the field area ( Figure 1), removing outliers from raw ECa values, which left a total of 2651 data points ( Figure 1). Then, the ECa-directed Response Surface Sampling Design module in ESAP-95 version 2.01 (Lesch et al. 2000) was used to delineate the scheme for 20 soil samples to be taken ( Figure 1). Soil cores were taken at the 0-30 and 30-60 cm soil depth. The samples were air-dried at 40°C and sieved at 2 mm diameter.

Soil physico-chemical analysis, and spatial variability
The 20 soil samples (200-250 g) at 0-30 and 30-60 cm depth were subjected to determination of the following properties: particle size distribution (sand, silt, and clay content), pH, total carbonates (CaCO 3 ), total organic carbon (C), total nitrogen (N), available P (P Olsen), exchangeable cations (K, Ca, Mg, Na), cation exchange capacity (CEC), and electrical conductivity of a soil extract with a 1:2.5 (w/w) soil-to-water ratio (EC 1:2.5 ). The particle size distribution was determined by the pipette method (Gee and Bauder 1986). Soil pH was measured at 1:2.5 (w/w) soil-to-water ratio. The total carbonate content (CaCO 3 ) was volumetrically determined (Loeppert and Suarez 1996). Total organic C and total N concentrations were determined by a CHN elemental analyzer (EA 1110 Thermo Fisher, Waltham, MA, USA). The available P was extracted according to Olsen et al. (1954) and was measured by inductively coupled plasma optical emission spectrometer (ICP-OES, Ametek, Spectro Arcos, Kleve, Germany). The cation exchange capacity (CEC) and the exchangeable cations were determined according to the method proposed by Orsini and Remy (1976) and modified by Ciesielski et al. (1997), and the amounts of Co and exchangeable cations were measured by ICP-OES. Soil electrical conductivity (EC) was determined on 1:2.5 (w/w) soil-to-water ratio aqueous suspension and then reported as EC on the saturation extract (ECe).
For soil spatial variability, we produced the maps of soil properties in the 0-60 cm soil depth (average of the 0-30 and 30-60 cm layers), by using ordinary kriging with 10 m grid resolution. Kriging outperforms normally the inverse-distance weighted method in spatial soil mapping (Kravchenko and Bullock 1999;Reza et al. 2010;Daniel et al. 2017).

Relationship between spatio-temporal YSCs and soil data
Thirty m wide buffers around the 20 positions determined by the ESAP software were created for statistical correlations between spatio-temporal yield and soil properties. The values of interpolated GY and selected soil properties falling within the range of each buffer were averaged for Pearson's correlations (r) involving the 20 data points. Thereafter, it was evaluated if multi-years spatio-temporal yield could effectively be described by the differences in soil properties within YSCs. To this aim, interpolated soil data were associated with the YSCs, then the statistical differences of soil properties among YSCs were assessed in the same way as described by Li et al. (2008) and Scudiero et al. (2018).
The weather information during the five growing seasons (Hydro-meteorological Service of the Emilia-Romagna region) was used to interpret temporal yield variability. The wet and dry periods from initial to maturity stages of the surveyed crops were represented by the balance between precipitation (P) and crop evapotranspiration (ET C ), this latter determined according to Allen et al. (1998). In the supplementary materials, total precipitation and the average temperature were computed monthly according to Bagnouls and Gaussen (1953), to indicate wet and dry periods during the five crop seasons.
Map production and geostatistical data analysis were carried out with the ArcGIS software (Version 10.3, ESRI, Redlands, CA, USA) under the reference system WGS 84/UTM zone 32N. Statistical analyses were performed with the Statistica 10 software (StatSoft Corp., Tulsa, OK, USA).

Statistical analysis
Crop yields and soil data were subjected to descriptive statistics. Pearson's correlation (r) was used to evaluate the relationships of soil properties and spatio-temporal relative yield in single and multiple crops. One way analysis of variance (ANOVA) was run to assess the differences in soil and yield traits among the three YSCs. The least significant difference (LSD) at P ≤ 0.05 was used to indicate significantly different levels. Table 1a summarizes the characteristics (mean, minimum, maximum, SD, kurtosis, and skewness) of standardized GY data in the five years. Crop yield varied greatly across the field. The widest min.-max. range (183) was found in BW 2012, whereas the tightest range (143) was shown in DW 2010. Standardized GY variability was generally high, as indicated by SD ranging from 29% for DW 2010 to 38% (SF 2011 andCO 2013).

Geostatistics of crop yields
The spatial patterns of crop yields were evaluated in terms of semivariograms and the respective model fittings (Table 1b). DW 2014 showed a zero nugget effect, followed by SF 2011 and DW 2010 with very low nugget values. All crops exhibited a quite similar total variance (sill variance (C 0 + C) ranging from 0.92 to 1.17), whereas the range (a) varied noticeably between 38 and 121 m. A high/ low range indicates high/low continuity, respectively, within the dataset. Based on the degree of SpD (Cambardella et al. 1994), crop data showed a 'strong' continuity in their SpD in all cases except BW 2012. The results of semivariogram model fitting (R 2 and MAE) confirmed the good performance of the stable and exponential variogram, depending on years, over the empirical data (Xiao et al. 2016;Bhunia et al. 2018).  (27) and lower maximum (148) relative GY, resulting in a narrower range (121) compared to single crops. Nevertheless, spatial variability maps in single vs. multiple crops were quite consistent, i.e. areas at high or low GY tended to repeat in the same position. The upper field portion (4.12 ha) always showed low and below-average yield, whereas the middle and lower field portions (6.95 ha) always featured above average and high yield.
( Figure 2g). Lastly, the three yield stability classes (HYS, LYS and unstable) depicted the features of spatio-temporal yield over multiple years (Figure 2h).

Yield data distribution within spatio-temporal yield stability classes
Data distribution within the classes of spatio-temporal maps is reported in Table 2.
For spatial variability classes, in DW 2010 GY data distribution was more skewed towards high yielding classes, meaning that more grid points belonged to the MH and VH classes. The same occurred, to a lesser extent, in all the other years except SF 2011. Also the multiple crops combined showed a slight prevalence of the MH and VH classes. The differences between the relative GY values Figure 2. Consolidated spatio-temporal yield maps: VL, very low; ML, medium-low; MH, medium-high; VH, very high; 2a, spatial variability map of DW 2010; 2b, spatial variability map of SF 2011; 2 c, spatial variability map of BW 2012; 2d, spatial variability map of CO 2013; 2e, spatial variability map of BW 2014; 2 f, spatial variability map over 5 years' multiple crops; 2 g, temporal variability map over 5 years' multiple crops; 2 h, yield stability classes over spatio-temporal variability over 5 years' multiple crops.
in single and multiple crop maps indicate that class limits and width are not the same between single and multiple datasets.
For yield stability classes, 527 data points (46% or 5.3 ha) were found in the HYS class, 352 (30% or 3.5 ha) in the LYS class, and 277 points (24% or 2.8 ha) in the unstable class.

Spatial variability of soil properties
The descriptive analysis of soil physico-chemical properties in the average of the two depths (0-60 cm) is reported in the supplementary material (Table S2). The soil was loamy, moderately alkaline, rich in carbonates, poor of organic carbon (<10 g/kg), and with a low C:N ratio (<10). Available P and exchangeable K were also quite low (<10 mg/kg and <0.3 cmol + /kg, respectively). Exchangeable Ca was relatively high (almost 80% of the CEC). ECe denoted a negligible salinity across the field. Lastly, the three particle size classes (sand, silt, and clay) were more heterogeneous (higher SD in proportion to mean data) than the rest of the soil properties.

Spatial soil maps
Spatial maps of soil properties interpolated with ordinary kriging are reported in Figure

Quality control of spatio-temporal YSCs
Significant correlations were evidenced between soil properties and the spatio-temporal yield data (Table 3)  Classification of stable soil physico-chemical properties within the three spatio-temporal yields classes depicted statistical differences of soil properties (Table 4).
The lowest mean value of sand (47.9%) was associated with the HYS class, resulting in maximum yield over multiple crops (YSCs). Higher sand content was found in LYS (56.8%) and unstable class (57.4%), which featured a similar but statistically different yield (80% and 83%, respectively). Conversely, silt, clay, and CEC had higher values in the HYS class, compared to LYS and unstable class.

Ambient conditions during the five cropping seasons
The balance of ambient moisture in the five crop seasons is reported in Table 5.
During DW 2010, the crop growing period from tillering to heading (initial to mid-season) received a surplus of 206 mm as P-ETc difference and was considered a wet period, whereas lateseason (ripening stage) staged a 68 mm deficit. BW 2012 and 2014 also received enough precipitation from tillering to stem elongation (28 and 212 mm surplus in the two respective years), whereas a dry period occurred from heading to ripening in both years. It resulted in a respective deficit of 298 and 249 mm. Compared to winter cereals, spring dicots SF 2011 and CO 2013 suffered an increasing water deficit across growth stages. At ripening, a cumulated deficit of 433 and 219 mm was attained in the two respective crops.
The representation of temperature and precipitation during the five growth seasons according to Bagnouls and Gaussen (1953) exhibits a similar picture (Figures S3 and S4).  The stronger drought experienced by the spring dicots vs. autumn cereals reflected in a stronger variation of yield data (Table 1a). Standard deviation was 38% of the mean GY in both SF 2011 and CO 2013, whereas it was 31%, averagely, in DW 2010, BW 2012, and BW 2014.

Discussion
The geostatistical analysis of the five crop yields featured quite similar parameters, despite using two different semivariogram models (Table 1b). BW 2012 represents the only partial exception, having a far more extended range than the rest of crops, associated with less spatial dependence. However, negligible nugget or barely exceeding 25% of the total sill, as in the case of BW 2012, indicates high spatial continuity between data points. This is a circumstance strengthening the value of the spatial variability maps obtained through kriging interpolation.
In these maps, the differences among the four yield classes that are evidenced in individual crops ( Figure 2f). Therefore, the multiple crops play a buffering role vs. single crops, meaning that operating with the former data is as a sounder basis for crop management decisions to be taken.
The three YSCs proposed by Blackmore (2000) set themselves one step beyond spatial variability maps, as they combine spatial and temporal variability into a single indicator. The unstable class is that deserving most attention, as it is an area with a potential for improving crop yields. In our case, this area covers almost one-fourth of the total field surface (Table 2). Additionally, the unstable class has a patchy distribution across the field, whereas the two stable classes, HYS and LYS, have a more consistent shape and distribution (Figure 2h). Lastly, the unstable and LYS classes denote an increase in sand content and a parallel decrease in silt and clay content, and CEC (Table 4). Hence it is sensed that sharper values in soil properties lead to fluctuation and depression of final yield, whereas more balanced values conduct to consistently higher crop yields (Table 4).
The ECa survey directs soil sampling towards field areas more prone to indicate variations in soil properties, in contrast to regular grid sampling . In our case, it is perceived that a higher density of sampling points was placed in the upper field portion (Figure 1), where the multiple-year data indicate systematic yield loss vs. the rest of the field (Figure 2f). Therefore, the ECa survey was shown able to predict soil constraints for plant growth. However, sampling and analysis were needed to detect the underlying causes, as premise for taking decisions to deal with constraints.
Overall, the mean values of spatial soil properties across the three YSCs show considerable differences and align with the multiple year yield map (Figure 2h). The considerable variations of crop yield across the whole field were quite well correlated with the variations of soil properties (Table 3); SF 2011 was a partial exception, but CO 2013, the other spring sown crop, behaved as the three winter wheat crops (DW 2010, BW 2012and 2014. The general good correlations between crop yield and soil properties are in accordance with the findings of . It is evinced from their work and ours how much it is important to understand the causes of yield variation through soil factors that are expected to contribute to crop productivity. Temporal stability is seen as a relevant property in multiple crops across multiple years. The unstable is made a class of its own in the YSC system (Blackmore 2000), to account for fluctuations which are due to crop type, weather, and undefined factors interacting with them ( Figure 2g). We believe that this field portion that gave a slightly higher yield (83%) than the LYS class (80%) (Table  4), averagely, may require separate cultivation practices, depending on in-season weather conditions during the specific crop season. The unstable class was characterized by similar values of soil properties as the LYS class (Table 4), indicating that these levels of soil properties are prone to reduce crop yield. Therefore, the unstabilizing effect associated with these properties is responsible for reducing crop yield in areas that could potentially give a high yield. Separate cultivation practices aimed at reducing the negative effects of soil constraints in such areas should provide a gain in crop yields.
The weather pattern during the five growing seasons (Table 5, Figures S3 and S4) provides some clues to understanding why some field areas gave higher yield than others, and why some other areas behaved differently over time. The five rainfed crops were exposed to irregular weather, and the erratic pattern of water availability is acknowledged as one of the main determinants of crop yield and its variability (Kang et al. 2009;Kukal and Irmak 2018). Yield losses consequential to drought are commonly reported in the literature for wheat (Karim et al. 2000;Mirzaei et al. 2011), as well as sunflower (Nel et al. 2001) and coriander (Unlukara et al. 2016). However, the effect of ambient moisture that we noticed on yield spatial variability cannot be ascertained in small plot experiments and is relatively novel in the literature.
Additionally, our study suggests that, under favorable weather in a specific year, unstable field zones could be managed as high yielding ones, i.e. supplying a non-limiting amount of inputs to harness the favourable conditions conducive to high yield. Conversely, under unfavourable weather, savings could be made to avoid inefficient use of crop inputs.
In other words, while the stable yield zones of a field should be managed by strategic planning, the unstable zones shall better be managed by tactical approach, e.g. based on crop growth status and soil moisture conditions, which are key factors for final yield in many agricultural areas around the world. Therefore, it is of paramount importance for the farmers to monitor their crop and receive timely we ather information for alternative decisions to be taken (Basso et al. 2011).
One last point concerns the management of fields that become larger and larger by merging previously separated fields, under the urge to increase the efficiency in agricultural practice. In the surveyed case, three consecutive fields, each approximately 200 m long in the north-south axis, were combined into a single field approximately 600 m long in 2006 ( Figure S2). The result is that the northern third, which was planted with deep rooted tree crops, is scarcely productive, once converted to annual crops: in the multiple year average, a relative yield of 76% can be calculated for the upper third, compared to 104% and 124% for the central and lower third.
There is no univocal answer to the dilemma whether to pursue the enlargement of crop fields, to the expenses of crop advocacy, or save advocacy, to the expenses of efficiency. In the former case, the assessment of the causes for low productivity in a field portion provides the grounds for applying the most suited crop husbandry in a site specific manner.

Conclusions
Site-specific zones are the basis of precision agriculture by understanding where variable crop inputs are needed, based on the spatial and temporal variability of field characteristics. This paper defines the concept of classified zones by delineating the potential yield stability classes based on spatiotemporal maps over a five-year series of yields obtained with different crops.
It is evinced that the field areas featuring unstable yield across years should be managed by considering the in-season weather information to predict whether the unstable field part will behave as high yielding or low yielding in a specific year. This will provide farmers with valuable support to decide the appropriate level of crop intensity, e.g. fertilizers or water supply, in a site-specific way.
Our work concludes that multi-year yield stability classes are a more practical and cost-effective approach than uniform management, setting the premise for variable inputs to optimize crop productivity. Based on yield stability classes, strategic and tactical decisions must be taken in different field areas, depending on the spatial and temporal profile of productivity owned by these areas.
However, despite encouraging results based on a five-year data from a mixed cropping system, it is quite likely that the surrounding conditions play a relevant role in each specific case. This makes the approach described in this work reproducible, not simply generally valid, in different crop conditions.