Evaluation of Soil Management Effect on Crop Productivity and Vegetation Indices Accuracy in Mediterranean Cereal-Based Cropping Systems

Mostly, precision agriculture applications include the acquisition and elaboration of images, and it is fundamental to understand how farmers’ practices, such as soil management, affect those images and relate to the vegetation index. We investigated how long-term conservation agriculture practices, in comparison with conventional practices, can affect the yield components and the accuracy of five vegetation indexes. The experimental site is a part of a long-term experiment established in 1994 and is still ongoing that consists of a rainfed 2-year rotation with durum wheat and maize, where two unfertilized soil managements were repeated in the same plots every year. This study shows the superiority of no tillage over conventional tillage for both nutritional and productive aspects on durum wheat. The soil management affects the vegetation indexes’ accuracy, which is related to the nitrogen nutrition status. No-tillage management, which is characterized by a higher content of soil organic matter and nitrogen availability into the soil, allows obtaining a higher accuracy than the conventional tillage. So, the users of multispectral cameras for precision agriculture applications must take into account the soil management, organic matter, and nitrogen content.


Introduction
Providing a sufficient amount of food to satisfy the nutritional demand of the current population is the essential goal of global agriculture. By 2050, the global population is estimated to reach 2.6 billion people [1], so food production must increase by at least 70% before 2050 to support continued population growth [2]. In modern agriculture, conventional tillage (CT) techniques have allowed the adoption of crops, especially on large surfaces ensuring high yields: the mixing of surface horizons in preparing the seedbed allows the stabilization of the main crop to the detriment of the weed competitors. However, this intensification of the crops, although necessary for responding to the food needs of the growing demographic pressure, is proving unsustainable: in fact, the increment of soil erosion [3,4], the use of water, energy, and fertilizers, the disruption of soil structure, and the reduction of water use efficiency [5] will probably increase the environmental and economic pressures posed by intensified agricultural activities [6]; therefore, the negative consequences for the environment are evident [7][8][9].
To lower the pressure of pollution and costs, agricultural conservation practices are gaining worldwide popularity for their ability to optimize productivity and reduce the impact on the land's natural resources [10].
In fact, reduced tillage and even no tillage (NT) bring benefits to the environment in terms of reduction of soil erosion, leaching of nitrates, reduction in the use of agricultural machinery, as well as a lower emission of greenhouse gases and fuel costs [11].
Furthermore, the low soil disturbance, with the addition of crop residues, increases the levels of humidity [12] and nutrients in the horizons of soil explored from the roots and the soil organic carbon [13,14], and it reduces the mineralization rate of the organic matter, nitrogen losses, and the soil erosion [15], so it's possible sustain long-term crop production [16][17][18].
These economic and environmental benefits underpin the three pillars of conservation agriculture (CA) such as NT, the adoption of crop rotations, and in-situ residue conservation and permanent soil cover [19].
Conservation practices are being studied on winter cereals, which are dominant crops in the Mediterranean semi-arid climate regions [20] where climate change is putting cereal yields at risk [21][22][23], and they are often penalized by extreme events such as long periods of extreme dryness alternated with a short heavy rainfall time. In the Mediterranean area, crop production can be improved with the adoption of CA techniques [10] and with the application of the right dose of nitrogen through the site-specific application of fertilizers [24].
To understand the phenological status and the soil during crop cycles, manual measurements of agronomic characteristics are necessary, but they are so labor intensive and time consuming [10]. As a solution, a modern farming management concept that responds to such challenge is Precision Agriculture (PA) [25,26], providing spatial and temporal data on the agricultural fields in a fast and economic way [27].
In fact, its remote sensing technology offers a more efficient way to obtain the large-scale mapping of plant parameters: the development of this technology is expected to increase the effectiveness of PA [28]. In particular, studies indicated that space-borne sensors can be used to obtain spatially extensive information from landscape at the global scale [29][30][31][32][33].
The contribution of spatial information technologies [50,51] defines site-specific management units (SSMU) that are useful for understanding the spatial variations of the crop, especially in terms of yield [52]. These variations are influenced by a multitude of factors including topographic, edaphic, biological, meteorological, and anthropogenic factors [53].
Climate change, as already mentioned, contributes to influencing this variability: in fact, in the Mediterranean area, there is a decrease in rainfall which, for example, influences the food activity of the microbial component of the soil [54,55], so it will be necessary to understand through the technology offered by PA the changes that take place in crop systems.
However, several studies show that this variability makes it complicated to use precision farming tools in and so often it is rather difficult to adapt them in farms that have to make lesser decisions [56][57][58]. As a consequence, precision farming technologies require support structures to facilitate learning and the reduction of uncertainty in the implementation and adaptation process [59,60].
The uncertainties detected with the instrumentation and the climate variability [54,61,62] join the information lack related to the evaluation of the soil management (SM) effect on the crop nutritional status and productivity through multispectral imagery.
Only recently [10] was it reported that the canopy height, cover, volume, and the Normalized Difference Vegetation Index (NDVI) calculated on cotton growth under NT was statistically higher than the cotton grown under CT. This suggests that soil management can influence not only the crop growth development, but also the NDVI values. The aim of this study is to describe the effect of different SM (NT versus CT) on the unfertilized durum wheat crop parameters, nutritional status, and VIs accuracy in order to draw up vegetation maps that are useful for the correct management of soil fertility and cropping systems productivity.

Experimental Site
The experimental site is located at the "Pasquale Rosati" experimental farm of the Polytechnic University of Marche in Agugliano, Italy (43 • 32 N,13 • 22 E, at an altitude of 100 m above sea level and a slope gradient of 10%), on a silty-clay soil classified as Calcaric Gleyic Cambisols [63] (Figure 1). Only recently [10] was it reported that the canopy height, cover, volume, and the Normalized Difference Vegetation Index (NDVI) calculated on cotton growth under NT was statistically higher than the cotton grown under CT. This suggests that soil management can influence not only the crop growth development, but also the NDVI values. The aim of this study is to describe the effect of different SM (NT versus CT) on the unfertilized durum wheat crop parameters, nutritional status, and VIs accuracy in order to draw up vegetation maps that are useful for the correct management of soil fertility and cropping systems productivity.

Experimental Site
The experimental site is located at the "Pasquale Rosati" experimental farm of the Polytechnic University of Marche in Agugliano, Italy (43°32′ N,13°22′ E, at an altitude of 100 m above sea level and a slope gradient of 10%), on a silty-clay soil classified as Calcaric Gleyic Cambisols [63] (Figure 1). The climate of the site is Mediterranean, on which was recorded a total rainfall of 801 mm between October 2017 and July 2018, while a contraction of 30% of rainfall was recorded during October 2018-July 2019 with a 560.8 mm of rainfall (Table 1).
In order to better represent the water dynamics into the soil-crop system, we estimated the monthly soil water balance (SWB) by using the following formulas (Equations (1) and (2)): The climate of the site is Mediterranean, on which was recorded a total rainfall of 801 mm between October 2017 and July 2018, while a contraction of 30% of rainfall was recorded during October 2018-July 2019 with a 560.8 mm of rainfall (Table 1). In order to better represent the water dynamics into the soil-crop system, we estimated the monthly soil water balance (SWB) by using the following formulas (Equations (1) and (2)): where P: monthly precipitation (mm); ETc: monthly crop evapotranspiration (mm); ETo: reference evapotranspiration calculated with the Hargraves formula (mm) [64]; Kc: crop coefficient [65] The soil water balance calculated during the 2017-2018 growing season was 230 mm higher than the 2018-2019 growing season (Table 1) with a marked difference in the February-March period. The average minimum air temperature was higher on October 2017-July 2018 than October 2018-July 2019 with values respectively of 10 • C and 9.7 • C. Otherwise, the average maximum air temperature was higher on October 2018-July 2019 than October 2017-July 2018 with values respectively of 18.2 • C and 17.9 • C.
Soil properties in compared experimental plots are indicated in Table 2. Soil sampling was made with a Hand Huger (mod. Suelo HA-3) immediately before sowing. From each subplot, 3 samples were taken for a total of 12 soil samples analyzed for each year.

Experimental Design and Crop Management
The experimental site is a part of a long-term experiment established in 1994 and is still ongoing [66] consisting of a rainfed 2 years rotation with durum wheat (Triticum turgidum L. var. durum cv. Grazia, ISEA) in rotation with maize (Zea Mays L., DK440 hybrid Dekalb Monsanto, FAO Class 300) [67].
Within each field, two soil management techniques (main plot, 1500 m 2 ) were repeated in the same plots every year and arranged according to a split plot experimental design with two replications. The conventional tillage (CT), which is representative of the business as usual tillage practice in the study area, was ploughed along the maximum slope every year by a moldboard (with 2 plows) at a depth of 40 cm in autumn. The seedbed was prepared with harrowing before the sowing date. The no-tillage (NT) soil was left undisturbed and was sprayed with herbicides before sowing prior to direct seed drilling. In this study, we will examine the unfertilized plots in order to describe the effect of different soil management techniques on the durum wheat crop parameters and on the crop nutritional status through the vegetation indices (VIs) computation. The dates (dd/mm/yy) of all the agronomic practices are reported in Table 3. Table 3. Agronomic management practices adopted during the two-year experimental period.

Measurements
At stem elongation and anthesis phenological stages (ZS35 and ZS60 respectively were ZS = Zadoks Scale [68]), we have measured crop parameters such as dry matter (g) and nitrogen (N) content (% and g m −2 ), and we have acquired multispectral images (MAIA S-2 multispectral camera) by using a UAV platform (DJI Matrice 600 pro) in order to compute the VIs algorithm. At crop maturity (ZS92), we measured the typical agronomic measurements, number of kernels per spike (KS), thousand kernel weight (TKW), and the grain yield (t ha −1 ) for both years under analysis in order to characterize the yield of the different soil management techniques.

Crop Parameters
For each plot, we have randomly selected three test areas (Figure 1). At each test area, we have manually cut and collected fresh plants biomass in a georeferenced 0.5 m long-row using the GNSS HiPer HR receiver (Topcon, Ancona, Italy) for a total of 48 ground control points (GCPs).
The fresh plant biomass was oven-dried at 80 • C for 48 h and then, we weighed the dry biomass (g). Before analyzing for total N content, we ground the dry biomass to pass a 0.5 mm.
The N content (%) was determined by automated combustion analysis Dumas method [69,70] in an oxygen-enriched atmosphere at a high temperature (EA 1110 LECO CHNS-0 analyzer, Leco Corporation, St. Joseph, MI) in order to ensure complete combustion of the whole sample.
Starting from the N content (%) results, we calculated the N content (g m −2 ) by using the following formula (Equation (3)

Yield Components
In order to characterize the yield obtained by the compared treatments, we measured at crop maturity (ZS92) the number of KS, the thousand kernels weight (TKW), and the grain yield (t ha −1 ).
The KS and the TKW were estimated on 30 spikes randomly collected per plot. The grain yield (t ha −1 ) expressed in dry matter was measured by using a laboratory thresher for the three test areas (1 m long-row) per plot.

Image Acquisition Processing
To generate the orthomosaic reflectance maps, we followed a process consisting of three steps: alignment and mosaicking of raw multispectral images, point cloud and mesh generation, and orthomosaic map export. For the first and third steps, we used the Pix4Dmapper (Pix4D, Lausanne, Switzerland), which is based on the structure from motion (SfM) algorithm [71]. This allows us to generate the orthomosaic reflectance map from the raw multispectral images acquired by each flight. For the second step, we used the geographical reference recorded by the D-RTK GNSS module equipped on the UAV platform. The newly generated orthomosaic reflectance map has been imported in QGis 3.4.8, an open source Geographic Information System, which was the software we used to complete the remaining two main steps of the image processing.
To complete the second main step, we inserted on QGis the GCPs by using a csv file format with the data source manager tool, and then we created for each GCP a polygon shape file of 0.085 m 2 , which corresponded to the sampling surface.
While in order to select the most relevant vegetation index (VI) calculated starting from multispectral imagery for precision agriculture application in a conservation agriculture context, we compered five vegetation index categories according with Xue and Su [72]. The VIs analyzed in this study are reported in the following Table 4. The VIs calculation was carried out by a "Raster calculator" of QGis 3.4.8, which allows performing calculations on the basis of existing raster pixel values, and the results are written to a new raster layer with a GDAL supported format. The extraction of the VIs values was performed by using the "zonal statistics plugin" of QGis 3.4.8 by using the polygon shape file created for each GCP.

Statistical Analysis
All statistical analysis was performed with R. To highlight the significant effect of soil management (SM), year (Y), and the SMxY factorial combination to all the crop parameters analyzed, we performed an analysis of variance (ANOVA) to a linear model generated by using the generalized least squares approach.
Before performing any statistical analysis to identify a significant difference between the two soil managements in analysis, we performed a Shapiro-Wilk W test to evaluate the normality of distribution. When the P-value of the Shapiro-Wilk W test was below 0.05, we assumed that the data are not normally distributed; otherwise, the data are considered normally distributed.
When data were normally distributed, we performed the Bartlett test, which is used to test if k samples are from populations with equal variances or not. If the p value of the Bartlett output test was below 0.05, we assumed that the k samples are not from populations with equal variances, and so we performed the Welch One-Way ANOVA to identify a significant difference between the treatments under study. When the p value of the Bartlett output test was greater than 0.05, we assumed that the k samples are from populations with equal variances, and so we performed the t-test independent samples (p value = 0.05) to identify significant differences between soil managements.
When data were not normally distributed, we performed the Levene test, which is used to check that variances are equal for all samples when your data come from a non-normal distribution. If the p value of the Levene test was below 0.05, we performed the Friedman Test to highlight the significant difference between the treatments under study. When the p value of the Levene test was higher than 0.05, we performed the Kruskal-Wallis test to identify a significant difference between the soil management techniques.
To evaluate if the soil management can affect the relationships between VIs and N content (g m −2 ), we performed a linear regression analysis that is used to identify the existence of significant relationships (*: p ≤ 0.05; **: p ≤ 0.01; ***: p ≤ 0.001). In addition, we reported the coefficient of determination (R2) and relative root mean square error (RMSE) for each relationship.

Crop Parameters
The ANOVA shows that the year (Y) factor has significantly affected all the crop parameters analyzed, while the soil management (SM) factor has significantly affected the nitrogen (N) content variables (% and g m −2 ).
For the N content (%), the ANOVA shows a significant effect of the interaction of year per soil management (Y x SM) ( Table 5). The 2019 year showed a significantly higher mean value of dry matter (DM) (g) and both N content variables (% and g m −2 ) than 2018 (Table 6), with a difference of 9.60 g, 0.74 and 2.70 for DM and N content (% and g m −2 ) respectively.
The no tillage (NT) showed a significantly higher N content (% and g m −2 ) than conventional tillage (CT) for both years (Table 6), which was equal to +0.63 for N content (%) and +0.76 for N content (g m −2 ) in 2018, and equal to +0.17 for N content (%) and +1.34 for N content (g m −2 ) in 2019.

Yield Components
The ANOVA shows for both year (Y) and soil management (SM) factors a significant effect on the yield components. In detail, the Y factor significantly affects the number of kernels per spike (KS) and the thousand kernel weight (TKW); the SM factor significantly affects the KS and grain yield (t ha −1 ) ( Table 7). No significant effect of Y x SM interaction was observed. The 2019 year showed a significantly higher value on the KS (+7) and a significantly lower value on the TKW (−7.7 g) than 2018, while no significant difference was observed for the grain yield (t ha −1 ) in the two-year comparison (Table 8).  The NT leads to a significantly higher value of the KS and grain yield (t ha −1 ) than CT in both the years under study (Table 8). In 2018, the NT obtained higher values of approximately 46% and 48% respectively for KS and grain yield than CT. While in 2019, the NT obtained higher values of approximately 35% and 35% respectively for KS and grain yield than CT.

Relationship between Vis and N Content (g m −2 )
In the growing season of 2018, the NT system showed an R 2 value of 0.81 on average and root mean square error (RMSE) of 0.57 on average, while the CT system showed an R 2 value of 0.31 on average and an RMSE of 0.58 on average. During the growing season of 2019, the NT system showed an R 2 value of 0.69 on average and RMSE of 1.44 on average; the CT system showed an R 2 value of 0.45 on average and an RMSE of 1.35 (Table 9). Table 9. Coefficient of determination (R2) and root mean square error (RMSE) between the calculated vegetation indices and the nitrogen content (g m −2 ) within variation.
The previous discussion can also be extended to each individual VI analyzed; in fact, the values of R 2 are always higher in NT than in CT in both growing seasons (Table 9).
Considering the 2018 year, we observed that Modified Soil-adjusted Vegetation Index (MSAVI2) is the most accurate VI, which reported a R 2 on the NT of 0.96 while for CT, the R 2 was 0.70. For 2019, we observed that the Normalized Difference Red Edge Index (NDRE) was the most accurate VI, which reported an R 2 on the NT system of 0.95 and an R 2 of 0.76 on the CT.
The NDRE and MSAVI2 are the only VIs that show a significant relationship with N content (g m −2 ) for both soil managements in each year.
By evaluating the average R 2 obtained for all the VIs analyzed for each year and soil managements, we reported that NDRE is the most accurate VI to be related with the N content (g m −2 ) with a mean R 2 value of 0.80 (Table 9).  The year 2019 showed a higher greenness than 2018 in each phenological stage; this is due to a significantly higher value of the DM (g) and N (% and g m −2 ) content (Table 6).

Vegetation Index Maps
Within the same phenological stage, in the comparison between different years, NT showed significantly greater levels of greenness attributable, as previously mentioned, to the greater content of N (% and g m −2 ), KS, and grain yield (t ha −1 ) in both years under study (Figures 2 and 3). The year 2019 showed a higher greenness than 2018 in each phenological stage; this is due to a significantly higher value of the DM (g) and N (% and g m −2 ) content (Table 6).
Within the same phenological stage, in the comparison between different years, NT showed significantly greater levels of greenness attributable, as previously mentioned, to the greater content of N (% and g m −2 ), KS, and grain yield (t ha −1 ) in both years under study (Figures 2 and 3).

Discussion
The year (Y) factor showed a significant impact on DM (g), KS (n) and TKW (g) as reported on the same experimental site by Seddaiu et al., 2016 [67] and on both N content variables (%N and g m −2 ). These results show, as described from several authors [78,79], that the development of durum wheat during the season is strongly influenced by the climatic trend; in fact, the rainfall recorded in 2017-2018 growing season was 30% higher than the rainfall observed during the 2018-2019 (Table 1) season, and this probably led to a higher N leaching, which implies a reduction in the availability of N for the crop [80].
The probable N leaching occurring during the 2017-2018 growing season is confirmed by the monthly-estimated soil water balance (Table 1), which showed a difference of 230 mm with respect to the 2018-2019 growing season.
The annual difference is especially concentrated in the February-March period (154 mm and 111 mm respectively), so this indicates that during this period, some of the nitrogen that was made available for soil organic matter mineralization may have been leached. All these consequences are much more accentuated in the CT because it has a greater porosity of the soil than NT where there is an increased number of soil micropores that facilitate the storage of soil moisture [81][82][83], a lower soil organic matter than NT ( Table 2) that plays a key role in water [84][85][86] and nutrient [87][88][89] retention also thanks to the mulching effect of the straw [88], as well as having no crop residues on the topsoil during the season due to the soil tillage, which involves a re-mixing of the horizons and consequently a dilution of the crop residues [90][91][92].
The year (Y) factor didn't have a significant impact on the grain yield (t ha −1 ), this result could be induced by its two intrinsic variables such as KS and TKW (g), where we observed a dynamic balance.
In the 2018 growing season, the KS showed a lower value than the 2019 growing season, which implies a lower nutritional availability, due to the higher rainfall recorded, and therefore less fertility of the spike.
For TKW, we observed an inverse behavior; in fact, lower KS values correspond to higher TKW values as described also by Mohammadi et al., 2013 [93], who reported a significant Pearson correlation value of −0.52 between KS and TKW.
During June 2019, the period in which the milk and dough kernel development is occurring, the maximum and minimum air temperature were higher than June 2018 (2.4 • C and 1.7 • C respectively) ( Table 1), this may have contributed to a greater loss of water from the caryopses with a consequent effect on the TKW reduction [94].
The NT involve a number of other advantages with respect to CT, such as reduction of the management costs of the company [97][98][99], increased fertility of the soil, and positive effects on soil biochemical properties and biomass microbial [92,[100][101][102], and this implies a stabilization of production in the medium to long term [103].
In contrast, the SM factor did not significantly affect the DM (g) and the TKW (g), confirming reports by De Vita et al., 2007 [104], according to which durum wheat grown at Vasto (Italy) did not show any significant difference in DM (g) and TKW for the years 2000 and 2001 for the Ct versus NT soil management analyzed.
The factorial combination of year and soil management (Y × SM) showed a significant effect (p ≤ 0.05%) on N content (%) as reported also by López-Bellido et al., 2013 [104].
Regarding the relationships between VIs and N content (g m −2 ), soil management shows a significant effect, as reported by Orsini et al. 2019a [66]. This may probably due to the greater amount of crop residues present on the NT system, which covers the soil surface, reducing soil disturbance [78] in the calculation of VIs starting from multispectral images.
Moreover, since the NT system is not disturbed by plowing, the residues of previous crops substantially increase water retention and consequently there is a greater availability of this element, thus determining greater crop development [5].
This dynamic is also confirmed by Ashapure et al. (2019) [10], who in cotton observed that the NT system, compared to the CT system, allows a significant increase on the NDVI (basic vegetation index category) in comparison with the CT system.
By evaluating the performance of the VIs to be related with the crop N content (g m −2 ), we suggest the use of NDRE and MSAVI2 to provide to farmers the vegetation index maps and the prescriptions maps for precision agriculture application.

Conclusions
The thermo-pluviometry trend strongly influences the development of durum wheat, both in yield and chemical composition.
This study shows the superiority of conservative agriculture over conventional agriculture for both nutritional and productive aspects on durum wheat.
We reported a dynamic balance on the yield components, in which KS and TKW are inversely proportional.
In addition, we confirmed again that the accuracy of VIs are related with the nitrogen nutrition status of durum wheat, and they also depend on the soil management. All the VIs analyzed obtained a higher accuracy in the NT system than in the CT system in both the years analyzed, which is due to the soil is not being disturbed by plowing and cultivation, previous crop residue substantially increasing water retention, and soil organic matter content contributing to higher plant growth and performance.
So, we advise to the potential users of multispectral images for precision agriculture application to take into account the soil management and related organic matter and nitrogen content into the soil.
In addition, we suggest the use of NDRE and MSAVI2 indices for durum wheat grown under a conservative agriculture context to provide vegetation maps and related prescription maps for the optimal monitoring of the nutritional status of durum wheat in Mediterranean agricultural contexts.