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

Abstract

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.

1. 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].

Using multispectral images collected by satellite, traditional aircraft, and unmanned aerial vehicles (UAVs), several studies [34,35,36,37] have examined vegetative conditions in agriculture.

In precision farming, UAVs are very widespread and are provided with multispectral cameras that measure different wavelength bands within visible and near infrared regions of the spectrum, which allow the formulation of a wide range of vegetation indices (VIs) informing on biomass [38,39], leaf area index [40,41], pigment content [42,43,44], nitrogen content [45,46], photosynthetic efficiency [47], water status [47,48], and cover (ground and residue) [49].

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.

2. Materials and Methods

2.1. 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. Thermo-pluviometric trend related to the durum wheat biological cycle during the experimental period.

Months	November	December	January	February	March	April	May	June	July	2017–2018
Rainfall (mm)										Total
2017–2018	124	96	29	173	143	37	95	48	57	802
2018–2019	42	61	70	22	36	59	165	1	105	561
Δ Rain	82	35	−41	151	107	−22	−70	47	−48	241
Soil Water Balance (mm)										Total
2017–2018	20	88	20	163	119	−44	−42	−8	−4	312
2018–2019	0	53	62	9	−4	−27	49	−51	−8	82
Δ Soil Water Balance	21	35	−42	154	123	−17	−91	43	4	229
Min air T (°C)										Average
2017–2018	7.9	4.1	5.2	2	5.7	11.8	14.9	17.7	20.5	10
2018–2019	9.3	3.7	2.4	4.5	7.4	9.1	11.7	19.1	20.3	9.7
Δ Min air T	−1.4	0.4	2.8	−2.5	−1.7	2.7	3.2	−1.4	0.2	0.3
Max air T (°C)										Average
2017–2018	11.1	11.9	12.8	8.3	13.3	21.5	23.7	27.8	30.8	17.9
2018–2019	11.9	11.1	9	13.2	17.7	18.6	20.3	30.5	31.5	18.2
Δ Max air T	−0.8	0.8	3.8	−4.9	−4.4	2.9	3.4	−2.7	−0.7	−0.3

).

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)):

$$SWB=P-ETc$$

(1)

$$ETc=ETo \left(Hargraves\right)\ast Kc (F\mathrm{AO})$$

(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 properties of the 0–20 cm layer in the conventional tillage (CT) and no tillage (NT) unfertilized plots in 2019.

Soil Proprieties	SM 1
	NT 2	CT 3
Sand (g kg−1)	127 (±21) a	120 (±19) a
Silt (g kg−1)	410 (±30) a	397 (±19) a
Clay (g kg−1)	463 (±36) a	483 (±22) a
SOM 4 (g kg−1)	18.0 (±2.8) a	13.2 (±2.1) b
Total nitrogen (g kg−1)	1.30 (±0.11) a	0.98 (±0.03) b

¹ SM: soil management; ² NT: no-tillage; ³ CT: conventional tillage; ⁴ SOM: soil organic matter. Within the same factor of variation, means that are followed by the same letter (a,b) are not significantly different at p < 0.05%.

. 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.

2.2. 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²) 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. Agronomic management practices adopted during the two-year experimental period.

Agro-Technique	SM 1	2017–2018	2018–2019
Ploughing (40 cm)	CT 2	02/10/2017	26/09/2018
Weed control: Glyphosate 3	NT 4	30/10/2017	28/09/2018
Harrowing and seed bed preparation	CT	20/11/2017	01/10/2018
Sowing 5	All	21/11/2017	30/11/2018
Weed control: Pinoxaden 6	CT	28/03/2018	08/03/2019
Pest control: Azoxystrobin, Cyproconazole 7	All	24/04/2018	22/04/2019
Harvest	All	06/07/2018	07/07/2019

¹ SM: soil management; ² CT: conventional tillage; ³ dose: 2.25 kg ha⁻¹ of active ingredient; ⁴ NT: no tillage; ⁵ Seed rate: 220 kg ha⁻¹; row spacing: 0.17 m; ⁶ 30 g ha⁻¹ of active ingredient; ⁷ dose: 0.16 l ha⁻¹ of active ingredient.

.

2.3. 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⁻²), 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⁻¹) for both years under analysis in order to characterize the yield of the different soil management techniques.

2.3.1. 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⁻²) by using the following formula (Equation (3)):

$$gN m^{-2}=\frac{\frac{N\% \ast dry matter \left(g\right)}{100}}{0.085 m^2}.$$

(3)

2.3.2. 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⁻¹).

The KS and the TKW were estimated on 30 spikes randomly collected per plot. The grain yield (t ha⁻¹) expressed in dry matter was measured by using a laboratory thresher for the three test areas (1 m long-row) per plot.

2.3.3. 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², 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. Agronomic management practices adopted during the two-year experimental period.

Vegetation Indices	Formula	References
ARVI 1	$\mathrm{ARVI}=\frac{\mathrm{NIR}- \mathrm{RB}}{\mathrm{NIR}+\mathrm{RB}}$ Where: $\mathrm{RB}= \mathrm{Red}-\mathrm{y}\left(\mathrm{Blu}-\mathrm{Red}\right)$	Korhonen et al. (2015) [73]
MSAVI2 2	$\mathrm{MSAVI}2=\frac{2\times \mathrm{NIR}+1-\sqrt{{\left(2\times \mathrm{NIR}+1\right)}^2-8\left(\mathrm{NIR}-\mathrm{Red}\right)}}{2}$	Leprieur et al. (2000) [74]
NDRE 3	$NDRE=\frac{NIR-Red Edge}{NIR+Red Edge}$	Barnes et al. (2000) [75]
VDVI 4	$\mathrm{VDVI}=\frac{2\times \mathrm{Green}-\mathrm{Red}-\mathrm{Blue}}{2\times \mathrm{Green}+\mathrm{Red}+\mathrm{Blue}}$	Wang et al. (2015) [76]
WDRVI 5	$\mathrm{WDRVI}=\frac{\mathrm{a}\times \mathrm{NIR}-\mathrm{Red}}{\mathrm{a}\times \mathrm{NIR}+\mathrm{Red}}$ Where:$\mathrm{a}=0.2$	Gitelson (2004) [77]

¹ ARVI: Atmospherically Resistant Vegetation Index; ² MSAVI2: Modified Soil-adjusted Vegetation Index; ³ NDRE: Normalized Difference Red Edge Index; ⁴ VDVI: Visible-Band Difference Vegetation Index; ⁵ WDRVI: Wide Dynamic Range Vegetation Index.

.

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.

2.4. 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⁻²), 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.

3. Results

3.1. 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⁻²).

For the N content (%), the ANOVA shows a significant effect of the interaction of year per soil management (Y x SM) (

Table 5. Results of the ANOVA applied to a linear model using generalized least squares for durum wheat.

Factor of Variation	df 1	DM 2	N Content
		g	%	g m−2
Y 3	20	**	***	***
SM 4	20	n.s.	***	*
Y × SM	20	n.s.	*	n.s.

¹ df: degree of freedom; ² DM: Dry Matter; ³ Y: Year; ⁴ SM: Soil management; *: Significant at p < 0.05%; **: Significant at p < 0.01%; ***: Significant at p < 0.001%; n.s.: not significant.

).

The 2019 year showed a significantly higher mean value of dry matter (DM) (g) and both N content variables (% and g m⁻²) than 2018 (

Table 6. Durum wheat crop parameters analyzed during the growing seasons 2018 and 2019.

Year	SM 1	DM 2	N Content
		g	%	g m−2
	NT 3	13.71 (±9.47) a	1.43 (±0.28) a	2.10 (±1.29) a
	CT 4	14.34 (±7.18) a	0.80 (±0.12) b	1.34 (±0.67) b
2018		14.02 (±8.23) B	1.11 (±0.38) B	1.72 (±1.08) B
	NT	26.92 (±18.95) a	1.93 (±0.51) a	5.09 (±2.71) a
	CT	20.32 (±13.53) a	1.76 (±0.35) b	3.75 (±2.02) b
2019		23.62 (±16.45) A	1.85 (±0.44) A	4.42 (±2.44) A

¹ SM: soil management; ² DM: dry matter; ³ NT: no-tillage; ⁴ CT: conventional tillage; means within columns that are followed by the same letter (lowercase letters for SM (a,b); uppercase letters for year (A,B)) are not significantly different at p < 0.05.

), with a difference of 9.60 g, 0.74 and 2.70 for DM and N content (% and g m⁻²) respectively.

The no tillage (NT) showed a significantly higher N content (% and g m⁻²) 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⁻²) in 2018, and equal to +0.17 for N content (%) and +1.34 for N content (g m⁻²) in 2019.

3.2. 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⁻¹) (

Table 7. Results of the ANOVA applied to a linear model using generalized least squares for durum wheat.

Factor of Variation	df 1	KS 2	TKW 3	Grain Yield
		n.	g	t ha−1
Y 4	20	***	***	n.s.
SM 5	20	***	n.s.	***
Y × SM	20	n.s.	n.s.	n.s.

¹ df: degree of freedom; ² KS: number of kernels per spike; ³ TKW: Thousand kernel weight; ⁴ Y: Year; ⁵ SM: Soil management; ***: Significant at p < 0.001%; n.s: not significant.

). 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⁻¹) in the two-year comparison (

Table 8. Crop yield parameter measured at crop maturity on the 2018 and 2019 years.

Year	SM 1	KS 2	TKW 3	Grain Yield
			g	t ha−1
	NT 4	13 (±2) a	52.2 (±0.9) a	2.5 (±0.2) a
	CT 5	7 (±1) b	52.8 (±1.1) a	1.3 (±0.2) b
2018		10 (±3) B	52.5 (±1.0) A	1.9 (±0.7) A
	NT	20 (±2) a	44.9 (±1.1) a	2.3 (±0.4) a
	CT	13 (±1) b	44.7 (±1.5) a	1.5 (±0.6) b
2019		17 (±4) A	44.8 (±1.3) B	1.9 (±0.7) A

¹ SM: soil management; ² KS: number of kernels per spike; ³ TKW: thousand kernel weight; ⁴ NT: no-tillage; ⁵ CT: conventional tillage; means within columns that are followed by the same letter (lowercase letters for SM (a,b); uppercase letters for year (A,B)) are not significantly different at p < 0.05.

).

The NT leads to a significantly higher value of the KS and grain yield (t ha⁻¹) 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.

3.3. Relationship between Vis and N Content (g m⁻²)

In the growing season of 2018, the NT system showed an R² value of 0.81 on average and root mean square error (RMSE) of 0.57 on average, while the CT system showed an R² 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² value of 0.69 on average and RMSE of 1.44 on average; the CT system showed an R² value of 0.45 on average and an RMSE of 1.35 (

Table 9. Coefficient of determination (R2) and root mean square error (RMSE) between the calculated vegetation indices and the nitrogen content (g m⁻²) within variation.

Vegetation Index	Year	Soil Management	N Content
			g m−2
			R2	RMSE 1
ARVI	2018	NT 2	0.80 ***	0.61
		CT 3	0.08	0.68
	2019	NT	0.73 **	1.48
		CT	0.48 *	1.53
MSAVI2	2018	NT	0.96 ***	0.28
		CT	0.70 **	0.39
	2019	NT	0.84 ***	1.15
		CT	0.42 *	1.61
NDRE	2018	NT	0.88 ***	0.47
		CT	0.59 **	0.45
	2019	NT	0.95 ***	0.62
		CT	0.76 **	0.04
VDVI	2018	NT	0.61 **	0.84
		CT	0.15	0.65
	2019	NT	0.13	1.98
		CT	0.11	2.67
WDRVI	2018	NT	0.78 **	0.64
		CT	0.01	0.71
	2019	NT	0.80 ***	1.28
		CT	0.44 *	1.59
Mean	2018	NT	0.81	0.57
		CT	0.31	0.58
	2019	NT	0.69	1.44
		CT	0.45	1.35

¹ RMSE: root mean square error; ² NT: no tillage; ³ CT: conventional tillage; *: significant at p < 0.05%; **: significant at p < 0.01%; ***: Significant at p < 0.001%.

).

The previous discussion can also be extended to each individual VI analyzed; in fact, the values of R² 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² on the NT of 0.96 while for CT, the R² was 0.70. For 2019, we observed that the Normalized Difference Red Edge Index (NDRE) was the most accurate VI, which reported an R² on the NT system of 0.95 and an R² of 0.76 on the CT.

The NDRE and MSAVI2 are the only VIs that show a significant relationship with N content (g m⁻²) for both soil managements in each year.

By evaluating the average R² 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⁻²) with a mean R² value of 0.80 (Table 9).

3.4. Vegetation Index Maps

Figure 2 and Figure 3 show the NDRE vegetation maps corresponding to stem elongation (ZS 35 phenological stage) and anthesis (ZS 60 phenological stage) for both growing seasons (2018–2019) when the durum wheat reaches the maximum vegetative development.

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⁻²) 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⁻²), KS, and grain yield (t ha⁻¹) in both years under study (Figure 2 and Figure 3).

4. 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⁻²). 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⁻¹), 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 soil management (SM) factor affected both N content variables (% and g m⁻²), KS, and grain yield (t ha⁻¹) as reported also by Orsini et al., 2019b [95] and Fiorentini et al., 2019 [96].

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⁻²), 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⁻²), we suggest the use of NDRE and MSAVI2 to provide to farmers the vegetation index maps and the prescriptions maps for precision agriculture application.

5. 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.