Simpliﬁed and Advanced Sentinel-2-Based Precision Nitrogen Management of Wheat

: This study compares simpliﬁed and advanced precision nitrogen (N) fertilization approaches for winter wheat relying on Sentinel-2 NDVI, grain yield maps, and protein content. Five N fertilization treatments were compared: (1) a standard rate, calculated by a typical N balance (Flat-N); (2) a variable rate calculated using a simpliﬁed linear model, adopting a proportional strategy (NDVI directly related) (Var-N-dir); (3) a variable rate calculated using a simpliﬁed linear model, adopting a compensative strategy (NDVI inversely related) (Var-N-inv); (4) a variable rate calculated using the AgroSat model (Var-N-Agrosat); and (5) a variable rate calculated applying the Agricolus model (Var-N-Agricolus). The study was carried out in four ﬁelds over two cropping seasons with a randomized blocks design. Results indicate that the weather remains the main factor inﬂuencing yield, as it typically happens in a rainfed crop. No substantial differences in crop yield were observed among the N fertilization models within each year and experimental location. However, in the more favorable season, the low-input direct model (Var-N-dir) resulted as the best choice, providing the higher NUE (nitrogen use efﬁciency) value. In the less favorable season, results showed a better performance of the advanced models (Var-N-Agricolus and Var-N-Agrosat), which limited yield losses and reduced intra-ﬁeld variability, with relevant importance given to the increasing frequency of abnormal climate phenomena. In general, all these VRT approaches allowed reduction of the excess of fertilizers, preservation of the environment, and saving money.


Introduction
Reducing inputs like irrigation water, fertilizers, pesticides, and herbicides, while increasing the efficiency of their use, represents modern agriculture's primary challenge in terms of limiting environmental impact [1] and keeping pace with world population growth [2]. This can be pursued by applying precision agriculture (PA) approaches [3], which aim at providing the various inputs where they are necessary, at the correct rate (limiting growth deficiencies, yield loss, surpluses, money loss, and pollution) and at the right time (synchronizing application with crop needs).
On a field scale, PA measures require monitoring of the crop with remote sensing. The European Space Agency (ESA)'s Sentinel-2 satellites, aimed at Earth observation and monitoring land surface variability, can be a powerful tool to provide relevant information in support of PA [4][5][6][7]. The Sentinel-2 satellites are equipped with multispectral sensors (MSIs), including 13 spectral bands, with a spatial resolution ranging from 10 m to 60 m. Images provided by Sentinel-2 satellites are publicly available for free through Copernicus mean value expressed as integers of standard deviations as a classification criterion per pixel. The user should input the algorithm by selecting the appropriate index, the number of homogeneous zones, and the distribution method. The selection of these inputs is strongly dependent on the user's experience and field-by-field evaluation.
The use of spreaders equipped with variable-rate technology (VRT) is one of the most important aspects of precision agriculture technology, allowing the modification of the distribution rate during the fertilization operation as needed, unlike traditional flat-rate spreaders that often imply over-or under-application of fertilizers. VRT spreaders can be map-based or proximal-sensor-based (on-the-go technology) [42], but the first approach is the most used, because on-the-go sensors are usually either too expensive, insufficiently accurate, or unavailable [43]. Several studies have investigated VRT technology's accuracy with encouraging results, reporting a very good match between the prescription map and the as-applied map [6,44].
Another widely used PA technique is based on grain yield mapping systems. This technology, useful for increasing knowledge about the spatial heterogeneity of yields, is based on the recording of georeferenced yield data recorded using different sensors mounted on the harvesters [45]: (1) harvest quantity measurement sensors (mass or volume); (2) GNSS sensors; (3) reference area measurement sensors (working width, speed, time); and (4) data recording and processing units. Yield mapping systems are constantly improving, but to date, their accuracy is influenced by many factors, which include flow sensor calibration, combine speed changes, grain flow variations, grain moisture, and the level of smoothing applied for mapping [46]. As a consequence, the outputs often require post-processing correction [47].
While many studies combine vegetation indices (from various remote or proximal sensors) and VRT crop fertilization, only a few integrate yield mapping sensors. Moreover, to our knowledge, very few studies have combined vegetation indices, VRT fertilization, and yield mapping in a PA case study designed explicitly for evaluating variable-N-rate treatments on winter wheat. In this context, two preliminary studies reported data from a two-year experiment aimed at comparing a standard rate calculated by a typical N balance and two variable rates defined by simplified NDVI-based models, in order to define the N rate to be supplied to winter wheat in the second application (i.e., in late winter, just before shoot elongation) [48,49]. The authors admitted that the models used to define the variable rates were very simplistic, and that more advanced approaches might have further improved the efficiency of fertilization.
In this context, the present work follows up by adding two advanced NDVI-based models, so that in four fields and two experimental years, a total of five approaches are evaluated: a flat N fertilization rate, calculated by a typical N balance, and four variable rates, derived from simplified and advanced approaches, all based on Sentinel-2 NDVI. Variable approaches include both proportional (NDVI directly related) and compensative strategies (NDVI inversely related). The overall aim is to compare all of these approaches in terms of NDVI, NUE (nitrogen use efficiency), grain yield, and protein content, and to identify the most responsive methods in the cultivation conditions of the experimental fields.

Study Areas and Crop Management
The experiments were carried out over two consecutive cropping seasons (2018-2019 and 2019-2020) in four different fields located in Umbria, Italy ( Figure 1). Two of the fields belong to the Fondazione per l'Istruzione Agraria (Fia, i.e., Foundation for Agricultural Education) of the University of Perugia, located near Deruta (Province of Perugia), in the middle Tiber Valley, while the other two cropping fields belong to the Sodalizio San Martino farm (Sod, i.e., Sodality of San Martino), located near Mugnano (Province of Perugia). All relevant field data are reported in Table 1. The crops were managed according to ordinary practices; weeds and diseases were controlled chemically. All of the precision on-field operations were performed using a tractor equipped with a Topcon GNSS automatic guide device connected to the regional real-time kinematic (RTK) network. The variable rate treatments were performed using a VRT fertilizer spreader connected through ISOBUS (a widely used software protocol compliant with ISO 11,783 standards) to a Topcon system console. Table 2 shows the main characteristics of the winter wheat varieties used in the study.     The crops were managed according to ordinary practices; weeds and diseases were controlled chemically. All of the precision on-field operations were performed using a tractor equipped with a Topcon GNSS automatic guide device connected to the regional real-time kinematic (RTK) network. The variable rate treatments were performed using a VRT fertilizer spreader connected through ISOBUS (a widely used software protocol compliant with ISO 11,783 standards) to a Topcon system console. Table 2 shows the main characteristics of the winter wheat varieties used in the study.

Soil Data
All of the fields were plain and characterized by a soil gradient related to their proximity to streams (located to the east of Fia19 and Fia20, and to the west of Sod19 and Sod20). Soil sampling, including 30 samples (0-40 cm depths) for each field, distributed according to a rectangular grid (about 1 sample/0.7 hectares), was performed on 14 November 2018 for Fia19, and on 17 December 2019 for Fia20. Each sample was obtained by mixing three bulks around the chosen point in order to improve the soil sampling reliability. The following parameters were analyzed: sand, loam, clay, organic matter, total nitrogen, and soil moisture. The particle size analysis was carried out following the modified pipette procedure [55]; for organic matter determination, the Walkley-Black method was used [56]; for N total determination, the Kjeldahl method [57] was used. Soil maps for clay and organic matter were produced through a B-spline interpolation procedure developed in QGIS (a free and open-source piece of geographic information system software [58]) version 3.10, using the B-spline approximation tool available in the processing plugin. According to the USDA soil classification [59], the soil of Fia19 is Fluventic Haplustepts, and belongs to the textural group silty clay loam, with increasing sand content from the west side to the east side. The soil of Fia20 is also Fluventic Haplustepts, but with a more heterogeneous texture falling in three main classes: loam, clay loam, and silty clay loam, with increasing sand content from the west side to the east side. Table 3 shows the main soil sampling statistics for the Fia19 and Fia20 fields.

Weather
Being an inner region of the Italian peninsula, Umbria is characterized by a humid subtropical climate (Cfa) according to Köppen-Geiger classification [60]. This region differs from a typical Mediterranean area in having quite cold winters (frost may occur several times between December and March) and high annual rainfall (700-800 mm), although late spring and summer are often dry. The cropping season 2018-2019 was unusually mild in winter and rainy in winter and spring ( Figure 2a). Conversely, the cropping season 2019-2020 was very rainy in November, which caused the delay of sowing in both the Fia and Sod fields, but then it was mild and very dry in winter and spring ( Figure 2b). Weather data was provided by the "FieldLab" of the Research Unit of Agronomy and Field Crops-Department of Agricultural, Food, and Environmental Sciences of the University of Perugia.

VRT Fertilization Models, Rate Calculation, and Experimental Design
According to the usual practice, the total N rate, applied as urea, was split into three applications in order to increase N uptake efficiency and limit the risk of N leaching. The VRT fertilization was applied only for the second fertilization, at initial shoot elongation, since it is widely reported in the literature that this is the application that is crucial for crop growth and yield [61][62][63]; in fact, it often represents the greatest part of the total rate. The first application, normally performed in order to promote tillering, is generally small because the crop's need is still low, and N would be more exposed to the risk of leaching; the third application, not always scheduled, is also a tiny part of the total rate, performed close to ear emergence and mainly aimed at increasing the grain protein content, while it has little effect on yield. Hence, this research focused on the second N fertilization, comparing the typical, flat N fertilization rate calculated by a typical N balance (Flat-N) with four variable rates, derived from simplified (Var-N-dir, Var-N-inv) and advanced approaches (Var-N-Agrosat, Var-N-Agricolus), all based on Sentinel-2 NDVI.
In Var-N-dir, the direct relationship between the NDVI and VRT N-rates was adopted according to a proportional strategy, on the assumption that the higher the NDVI of the crop before the second N fertilization (Feekes stages 7-8), the higher the ability of the crop to valorize the N supply [16,23,64]. The ranges for N rates in the two years were defined assuming a maximum reduction in yield of 30% in disadvantaged areas:

VRT Fertilization Models, Rate Calculation, and Experimental Design
According to the usual practice, the total N rate, applied as urea, was split into three applications in order to increase N uptake efficiency and limit the risk of N leaching. The VRT fertilization was applied only for the second fertilization, at initial shoot elongation, since it is widely reported in the literature that this is the application that is crucial for crop growth and yield [61][62][63]; in fact, it often represents the greatest part of the total rate. The first application, normally performed in order to promote tillering, is generally small because the crop's need is still low, and N would be more exposed to the risk of leaching; the third application, not always scheduled, is also a tiny part of the total rate, performed close to ear emergence and mainly aimed at increasing the grain protein content, while it has little effect on yield. Hence, this research focused on the second N fertilization, comparing the typical, flat N fertilization rate calculated by a typical N balance (Flat-N) with four variable rates, derived from simplified (Var-N-dir, Var-N-inv) and advanced approaches (Var-N-Agrosat, Var-N-Agricolus), all based on Sentinel-2 NDVI.
In Var-N-dir, the direct relationship between the NDVI and VRT N-rates was adopted according to a proportional strategy, on the assumption that the higher the NDVI of the crop before the second N fertilization (Feekes stages 7-8), the higher the ability of the crop to valorize the N supply [16,23,64]. The ranges for N rates in the two years were defined assuming a maximum reduction in yield of 30% in disadvantaged areas: therefore, the maximum N rate was associated with the maximum potential yield, achievable in the portions of the field with maximum NDVI; the minimum rate was associated with the minimum potential yield, revealed by the minimum NDVI value; all intermediate values were defined by a linear model between the two rates, according to the NDVI gradient.
Concerning Var-N-inv, the inverse relationship between the NDVI and VRT N-rates was adopted according to a compensative strategy, on the assumption that the lower the NDVI of the crop before the second N fertilization (Feekes stages 7-8), the higher the N need of the crop in order to recover the gap in growth so as to obtain the maximum yield-whereas the crop with high NDVI is already well N-fed, and additional high N supply would cause luxury consumption. The ranges for N rates in the two years were defined by fixing a reduction of 30% of the maximum rate where the NDVI was maximum, based on the evidence from a previous study in the same location [26]. Intermediate N rates were defined by a linear model between the two rates, according to the NDVI gradient.
The Var-N-Agrosat treatment is based on the AgroSat model, which provides a nitrogen prescription map based on the user's intentions in terms of nitrogen units to apply to the crop and the % equivalent nitrogen content of the fertilizer. If the user is using the tool without AgroSat registration or any information about the crop sown and sowing date, the nitrogen prescription map tool calculates the prescription map based on the latest available NDVI map for the user's field; on the other hand, if the user is registered and the crop and sowing date data are available, the tool calculates the nitrogen prescription map considering the NDVI trend from the sowing date to the current analysis date. The output is available at 10 m of spatial resolution, and can be reclassified in 2-3-4 classes. The model behind this tool is an exponential function (NDVI vs. nitrogen rate) as a derivation of a more complex model built with NDVI, NNI (nitrogen nutrition index), chlorophyll red-edge index, and data acquired in the field (SPAD, Dualex), in the laboratory (leaf N content), and from literature [38,[65][66][67][68][69]. The main purpose of AgroSat is not to advise how many units to apply (this is defined by the user), but how to redistribute them throughout the field according to the dictates of precision agriculture. The aim is to homogenize the quantity and quality of crop production. The average of the distribution is always less than or equal to the value indicated by the user, while the variations within the field are limited to a maximum of ±30% of the value indicated by the user.
The Var-N-Agricolus model strategy is based on the last available NDVI calculated from Sentinel-2, and on the average N needs of the field calculated by the N balance model. N provided to the field in the first and third fertilization interventions was subtracted from the model results. Satellite indices and texture data (percentage of sand, silt, and clay) resulting from soil analysis were compared, in both years, in order to spatialize the resulting N rate. Since, in 2019, the water indices indicated no substantial influence on vigor, the field variability in terms of vigor was attributable to a delay of emergence in some areas of the field. The NDVI from different dates was compared in order to identify these areas; then, the seasonal balance model was run, assuming that in the areas where the NDVI was higher, a higher potential yield-and so a higher N uptake-would have occurred. In 2020, the following observations were considered: (1) The spatial trend of the NDVI was inverse to the spatial tends of NDMI (normalized difference moisture index) and TCARI/OSAVI (chlorophyll index); and (2) in the areas where the NDVI was higher, the soil analysis indicated a higher percentage of sand. In these areas, water stress and chlorophyll were lower. These factors were contrasting, as sandy soils and low chlorophyll usually require a higher N supply, while in the areas with high vigor and low water stress, low N supply is usually enough. However, during the application of the N balance model, considering the soil texture gradient, the higher N availability from organic matter mineralization induced us to decrease the nitrogen dose slightly (±5% from the mean) in places where the NDVI was higher.
The above-described treatments were combined in the four experimental fields during the two years ( Table 4). The number of treatments compared in the Fia and Sod fields was constrained by the fields' size and the unavailability of a combine equipped with a yield tracking system for the Sod fields. These constraints, as explained later, also influenced the experimental units' size, the number of replications of treatments in each field, and the harvest methods. The NDVI data used for the Var-N-dir strategy were based on Sentinel-2 L2A images from 22 March 2019 for Fia19, and from 19 March 2020 for Fia20 ( Figure 3). The NDVI data used for the Var-N-inv strategy were based on Sentinel-2 L2A images from 2 March 2019 for Sod19, and from 19 March 2020 for Sod20 ( Figure 3). All of these images were cloud-free and the closest to the fertilizer treatments.
The total N rates based on the N balance used for the Flat-N thesis were 180 kg ha −1 for Fia19, 190 kg ha −1 for Fia20, 170 kg ha −1 for Sod19, and 160 kg ha −1 for Sod20. Table 5 reports the amounts supplied in the second N application based on the different abovedescribed approaches. general, the experimental units' size and the relatively low speed of the tractor (about 10 km/h) during the operation limited the uncertainty related to the VRT spreader's possible delays occurring while changing the N doses. The prescription maps used for the second fertilization, containing the experimental units and the associated N doses expressed as urea (Kg/ha), were built in QGIS, in shapefile format, composing the outputs from the above-described approaches (Figure 3). The shapefile was uploaded on the Topcon console before executing the second N application (Table 1). The fields were divided into several plots (our experimental units) in order to exclude the effects of unknown or unmeasured factors, and the treatments were laid down according to a randomized block design. Fia19 was divided into 52 plots of about 2400 m 2 each (60 m long, 40 m wide), and the 4 treatments were replicated in 13 blocks. Fia20 was divided into 60 plots of about 2800 m 2 each (70 m long, 40 m wide), and the 4 treatments were replicated in 15 blocks. Sod19 was divided into 12 plots of about 4400 m 2 each (110 m long, 40 m wide), and the 3 treatments were replicated in 4 blocks. Sod20 was divided into 18 plots of about 4320 m 2 each (120 m long, 36 m wide), and the 3 treatments were replicated in 6 blocks. For both Sod19 and Sod20, given the experimental units' size, the VRT treatments were variated within the units following a subgrid of about 20 × 20 m. In general, the experimental units' size and the relatively low speed of the tractor (about 10 km/h) during the operation limited the uncertainty related to the VRT spreader's possible delays occurring while changing the N doses. The prescription maps used for the second fertilization, containing the experimental units and the associated N doses expressed as urea (Kg/ha), were built in QGIS, in shapefile format, composing the outputs from the above-described approaches (Figure 3). The shapefile was uploaded on the Topcon console before executing the second N application (Table 1).

NDVI Seasonal Trend
All of the level-2A Sentinel-2 images with no cloud cover on the study areas were collected from the second fertilization to the harvest, in order to monitor and compare the NDVI of experimental treatments (15 images for Fia19, 16 images for Fia20, 13 for Sod19, 19 for Sod20). Average NDVI values were calculated for each experimental unit using the SAGA "raster statistics for polygons" tool included in the QGIS processing framework. This tool allows the user to calculate statistics about raster cell values contained in zones defined by a vector layer of polygons.
Another analysis was performed by classifying plots according to three classes, defined by dividing the NDVI range revealed by the S2 image acquired before the second fertilization into three equal parts. In addition, total nitrogen distribution and yield data were retrieved and compared for each of these groupings.
A final Pearson correlation analysis between the NDVI data and the soil variables (for Fia fields) and yield variables (for Fia and Sod fields) was developed in order to explore the possibility of predicting these variables using this vegetation index.

Grain Yield
For Fia19, the harvest was carried out on 12 and 14 July 2019, and for Fia20 on 9 July 2020, using a Claas Lexion 630 combine harvester equipped with an automatic guide device connected to the regional RTK network and a Topcon YieldTrakk system (based on grain mass flow and moisture sensors). After the input of the specific grain weight by the operator, the system produced a georeferenced yield map for the whole field as an ESRI polygon shapefile. The resulting polygons' size depends on the feed rate, the collected mass, and the cutter bar's working width (in this case, 7 m). The system also corrects the base width of the polygons when the combine passes over already collected field portions.
A data cleaning technique was applied to yield data according to the following steps: (1) all values outside the general average ± 2.5 * standard deviation range were eliminated, statistically considered outliers [70]; (2) all values collected close to the drainage system were eliminated-this is because the system in those areas records the entire cutter bar as operational, while the portion above the drainage system does not collect anything; and (3) all values falling in the 10 m portion before and after the experimental units' transitions were eliminated, in order to avoid the values of one experimental unit being counted in the adjacent one due to the delay time. The delay time of the yield tracking system is the time required for the harvested grain to reach the reading sensor positioned on the hopper. This parameter is necessary in order to associate the harvested product with the correct geographical position instead of the position in which the combine is when the product reaches the reading sensor. The system sets the default delay time to 9 seconds, but the actual delay time between harvesting and reading depends on many factors, including threshing rotation speed, crop density, the state of maturation, and the moisture content of the wheat [71]. The resulting cleaned data were intersected with the experimental units in order to calculate the average yield expressed in Mg ha −1 , which was used for the comparative analysis between treatments.
For the Sod fields, the harvest, carried out on 8 July 2019 and 23 July 2020, was performed by collecting and weighing each experimental unit's yield. Although very timeconsuming, this procedure provided very accurate yield data for each experimental unit and treatment. For Sod19, the experimental unit "1" (Var-N-Dir) was successively excluded from the experiment due to water stagnation located in that area, which compromised the normal comparison conditions.

Protein Content
For all fields, grain sampling was performed a few days before harvest in order to determine the grain protein content (Fia19: 27 June 2019; Fia20: 9 July 2020; Sod19: 26 June 2019; Sod20: 10 July 2020). For Fia19 and Fia20, a sample from each experimental unit was collected, for a total of 52 samples and 60 samples, respectively. For Sod19, 2 samples per experimental unit were collected (considering their wider area), for a total of 24 samples, while for Sod20, a sample per experimental unit was collected, for a total of 18 samples. Each sample was made up of four elementary subsamples randomly taken in the plot. All samples were analyzed according to the official Kjeldahl method [57].

Statistical Analysis
The inferential statistical analysis in this study was conducted separately for the four fields, considering them as independent experiments. All analyses were conducted considering the mean values of the variables calculated for each experimental unit. The analysis of variance (ANOVA) and Tukey's HSD test (0.05 level of significance) were performed for the main variables in order to test the significance of differences between the treatments within each field. Pearson correlation was carried out in order to test the correlation between selected variables. All of the statistical analysis was performed using the statistics software R [72] (Version 4.0.3) and the "agricolae" package [73] (Version 1.3.3).

Results
The seasonal trend of crop NDVI variation for all treatments in each of the four fields is reported in Figure 4. Differences among treatments in NDVI values (average at field scale) were practically null in Fia fields and little in Sod fields, despite the N supplied by the second N application differing by up to three rates of kg N ha −1 (Table 5). Instead, the interannual difference was macroscopic in terms of within-crop-season trends and peak values.

Protein Content
For all fields, grain sampling was performed a few days before harvest in order to determine the grain protein content (Fia19: 27 June 2019; Fia20: 9 July 2020; Sod19: 26 June 2019; Sod20: 10 July 2020). For Fia19 and Fia20, a sample from each experimental unit was collected, for a total of 52 samples and 60 samples, respectively. For Sod19, 2 samples per experimental unit were collected (considering their wider area), for a total of 24 samples, while for Sod20, a sample per experimental unit was collected, for a total of 18 samples. Each sample was made up of four elementary subsamples randomly taken in the plot. All samples were analyzed according to the official Kjeldahl method [57].

Statistical Analysis
The inferential statistical analysis in this study was conducted separately for the four fields, considering them as independent experiments. All analyses were conducted considering the mean values of the variables calculated for each experimental unit. The analysis of variance (ANOVA) and Tukey's HSD test (0.05 level of significance) were performed for the main variables in order to test the significance of differences between the treatments within each field. Pearson correlation was carried out in order to test the correlation between selected variables. All of the statistical analysis was performed using the statistics software R [72] (Version 4.0.3) and the "agricolae" package [73] (Version 1.3.3).

Results
The seasonal trend of crop NDVI variation for all treatments in each of the four fields is reported in Figure 4. Differences among treatments in NDVI values (average at field scale) were practically null in Fia fields and little in Sod fields, despite the N supplied by the second N application differing by up to three rates of kg N ha −1 (Table 5). Instead, the interannual difference was macroscopic in terms of within-crop-season trends and peak values.  As far as grain yield is concerned, differences were always very little and never significant (p-value > 0.10) (Table 6, Figure 5). Consequently, based on the N supplied by the second N application, the NUE (nitrogen use efficiency-derived from the ratio between kg of harvested grain and kg of nitrogen input) was higher the lower the N supplied. Var-N-inv and Var-N-dir showed the highest NUE (p-value < 0.05), followed by Var-N-Agricolus and then by Flat-N and Var-N-Agrosat, which did not differ substantially. The high variability in Var-N-Agricolus for Fia19 accounted for the lack of significant differences between Var-N-Agrosat and Flat-N. No significant differences were recorded for protein content, whose range for all treatments within the same field was in the order of 0.5% or less.  As far as grain yield is concerned, differences were always very little and never significant (p-value > 0.10) (Table 6, Figure 5). Consequently, based on the N supplied by the second N application, the NUE (nitrogen use efficiency-derived from the ratio between kg of harvested grain and kg of nitrogen input) was higher the lower the N supplied. Var-N-inv and Var-N-dir showed the highest NUE (p-value < 0.05), followed by Var-N-Agricolus and then by Flat-N and Var-N-Agrosat, which did not differ substantially. The high variability in Var-N-Agricolus for Fia19 accounted for the lack of significant differences between Var-N-Agrosat and Flat-N. No significant differences were recorded for protein content, whose range for all treatments within the same field was in the order of 0.5% or less.  For Fia19, the NDVI time series analysis was performed by classifying plots according to three NDVI classes (NDVI ≤ 0.70; 0.70 < NDVI ≤ 0.77; NDVI > 0.77), as revealed by the S2 image from 22 March 2019 (Table 7, Figure 6). For Fia20, the NDVI time series analysis was performed by classifying plots according to three NDVI classes (NDVI ≤ 0.40; 0.40 < NDVI ≤ 0.44; NDVI > 0.44), as revealed by the S2 image from 19 March 2020 (Table 7, Figure 6). The differences in terms of yield were not significant in either field, while for the Fia20 field, the May 23 NDVI values of Var-N-Agricolus and Var-N-Agrosat were significantly different from Var-N-dir, but not from Flat-N. For Fia19, the NDVI time series analysis was performed by classifying plots according to three NDVI classes (NDVI ≤ 0.70; 0.70 < NDVI ≤ 0.77; NDVI > 0.77), as revealed by the S2 image from 22 March 2019 (Table 7, Figure 6). For Fia20, the NDVI time series analysis was performed by classifying plots according to three NDVI classes (NDVI ≤ 0.40; 0.40 < NDVI ≤ 0.44; NDVI > 0.44), as revealed by the S2 image from 19 March 2020 (Table 7, Figure 6). The differences in terms of yield were not significant in either field, while for the Fia20 field, the May 23 NDVI values of Var-N-Agricolus and Var-N-Agrosat were significantly different from Var-N-dir, but not from Flat-N.

Correlation Analysis between NDVI, Yield, and Soil Data
The results of the Pearson correlation analysis between the NDVI data and soil variables (for Fia fields) and yield variables (for Fia and Sod fields) are reported in Tables 8-11. The correlations are based on the mean values from the individual plots.
Our results indicate relevant correlations only in few cases. For Fia20, the NDVI showed a high correlation with the yield between the last weeks of May and 1 June. Moreover, other fields showed an increase in correlation during the same period, which dropped quickly after that-except for Fia19, which showed a shifted peak towards mid-June. Table 8. Pearson correlation analysis and significance levels (. = 0.10; * = 0.05; ** = 0.01; *** = 0.001) between NDVI, grain yield, protein content, and soil parameters for the Fia19 field.

Fia19
Grain   Table 9. Pearson correlation analysis and significance levels (. = 0.10; * = 0.05; ** = 0.01) between NDVI, grain yield, and protein content for the Sod19 field.  Concerning grain protein content, for the 2018/19 season, there seems to be a correlation peak for both fields in the period around the end of March, and a subsequent second peak, with Fia19 reaching a correlation of 0.4 in early June and Sod19 reaching a peak of 0.8 between the second and third weeks of June. Fia20 shows its peak towards the last week of May, with a high average correlation with the NDVI of the whole season, while Sod20 starts out having a negative correlation for the whole first part of the season, becoming positive from June onwards.

Sod19 Grain Protein Content Grain Yield
For Fia19 and Fia20, a Pearson correlation index between grain yield and the soil parameters (clay, organic matter, and total nitrogen content) was also explored. A significative inverse correlation between yield and organic matter was found in the Fia19field.

Discussion
Our results show a substantial difference at the field scale, mainly linked to the different climatic conditions in the two cropping years, with respect to the effects of the different fertilization approaches. Our results show that a "Var" approach with a lower overall N rate (Var-N-inv, Var-N-dir, Var-N-Agricolus) gave the same grain yield with lower N input in all fields (Table 6). In fact, the NUE calculated with respect to the grain yield was significantly (p-value < 0.05) higher compared to Flat-N, Var-N-Agrosat, and Var-N-Agricolus in all cases. This result is consistent with the work of Raun et al. [74], who reported that the VRT method improved the NUE (nitrogen use efficiency) by 15% compared to a flat rate. The results are also consistent with the work of Vizzari et al. [6], who reported no decrease in grain yield using "Var" approaches while reducing the rate by about 20 kg N ha −1 . This evidence is economically and environmentally relevant, since the surplus N would involve luxury consumption or leaching risks, depending on whether or not it would be absorbed by the wheat's roots. Thanks to the optimal climatic conditions of the 2018-2019 cropping season, it is evident that even the most conservative approach in terms of nitrogen intake worked in saturation conditions anyway. Moreover, an increase in grain quality may not be advocated to justify the N surplus, since the increase in the grain protein content of the flat-N and Var-N-Agrosat treatments was negligible and never significant.
The lack of relevant average differences in grain yield among treatments is in line with the NDVI values recorded after the second N application, since the differences in NDVI between treatments were either null (2019) or tiny and non-critical (2020). This occurrence may happen, especially for high NDVI values [26], when other variables aside from N nutritional status come into play, such as the soil's hydrological properties [31,75] or the rainfall regime.
In this regard, both experimental years were abnormal. The first was exceptionally rainy, and the second exceptionally dry, which masked the N treatments' effect. Thus, for example, the negative correlation found between grain yield and the clay and organic matter content of the soil in the 2018-2019 season was likely independent of the soil's nutrient content, being more likely due to soil water saturation caused by the abundant rainfall. On the other hand, in the 2019-2020 season, the rainfall in winter and spring was so scarce that the wheat growth (after shoot elongation) and the grain filling were limited by water rather than N availability. The small gradient in clay content could not have counteracted the severe water deficiency. Interestingly, in this year, both the grain protein content and the yield were lower than in the previous year. This result was unexpected since, normally, when the crop is well N-fed (as it was in our case), biomass production is inversely correlated to the concentration of N compounds-including protein production [76]. In this regard, we can only speculate that the very low water availability of 2020 somehow impaired N absorption or relocation to grains [77].
Regarding the correlation between yield and NDVI (Tables 8 and 9), the strong correlation observed in the Fia20 field could be related to the unfavorable season, which did not cause typical NDVI saturation. Consequently, the intra-field heterogeneity was indicative of the crop's phenological growth and the LAI, which was reflected in the final yield. Sod19 and Sod20 showed correlation peaks in the same period as Fia20, but with a lower Pearson value of around 0.4. This can be explained by the fact that the experimental units in these fields were very large compared to the NDVI gradient, and with a higher yield data resolution, the correlation would likely have been stronger. These results are in line with the literature, where the existence of a correlation between NDVI and yield is widely reported [22,37,78]. Fia19 showed a low correlation between NDVI and grain yield, in contrast to other fields and most of the relevant literature. This result could be explained by the optimal conditions of the season 2018-19 having translated into high crop vigor, resulting in a low intra-field variability and, above all, a saturation of the NDVI, no longer able to detect the differences in vigor above a certain threshold. This flattening of the variability could result in the low correlation found with the yield data.
Regarding the correlations between protein content and NDVI (Tables 10 and 11), although in this case the correlations are on average high for all fields, it is difficult to find a causality that explains these data, starting with the Sod20 field which, in contrast to all of the others, seems to show a negative correlation with NDVI. The significative inverse correlation between yield and organic matter in the Fia19 field is atypical, and in contrast to literature, but could be explained due to the overlapping of the clay and organic matter spatial trends along the field. The strong clayey character of the soil can cause compaction, with negative consequences on yield [79].
Any evaluation on the suitability of different Var-N-approaches (i.e., -dir, -inv, -Agrosat, -Agricolus) appears meaningless if not linked to the climatic conditions character-izing the cropping seasons and to a careful evaluation of the initial and final conditions in terms of NDVI and final yield. During the first cropping season, the optimal climatic conditions meant that reducing the overall N rate (i.e., Var-N-dir, Var-N-inv, Var-N-Agricolus) was the most suitable strategy for increasing the nitrogen use efficiency without compromising grain yield and quality. Conversely, during the second cropping season, a direct approach (Var-N-dir) did not lead to the same results. In fact, interesting evidence in support of this conclusion could be found in the most disadvantaged areas of the Fia20 field; what emerges is that, when the crop season is unfavorable, well-studied management may make the difference in limiting N losses. Analyzing the NDVI class "I" (including all of the field pixels with an NDVI value ≤ 0.40 on 19 March 2020), it emerges that, even if the differences are slight and not statistically significant in terms of yield, and additional experiments would be necessary, the advanced model strategies (Var-N-Agricolus and Var-N-Agrosat) seem to give the best results both in terms of NDVI (Figure 6, bottom) and in terms of yield (Table 7), also resulting in a more spatial homogeneous production across the field. Of the two advanced models, Agricolus resulted in the highest NUE ( Table 6). The probable reading of this result is that a full nitrogen dose was not the best choice for that year and that area, because the crop could not exploit it given the weather and climatic limits. At the same time, limiting the nitrogen rate to a minimum with a direct approach did not turn out to be a good strategy, because the crop would have been able to exploit a slightly higher dose, even if not the maximum.

Conclusions
Our study, carried out in two fields and two cropping seasons, investigated the differences between a flat N fertilization rate and four variable rates defined through simplified and complex methods, based on Sentinel-2 NDVI, applied to define the second fertilization. Three of the four variable rates were aimed at reducing the overall N input (Var-N-inv, Var-N-dir, and Var-N-Agricolus), while he other was aimed at maintaining the same overall input (Var-N-Agrosat) as in the flat rate (Flat-N) while optimizing the N fertilization according to the supposed N requirement as revealed by the NDVI data.
Our results show that the variable approaches with lower overall N rates (Var-N-inv, Var-N-dir, and Var-N-Agricolus) gave the same grain yield as the others (flat-N, Var-N-Agosat) while obtaining a higher NUE. However, the various approaches can be correctly evaluated only when contextualized to the respective climatic season. In fact, the strong effect of the opposite climatic trend across the two seasons is an explicit confirmation that weather is the main factor influencing the yield of a rainfed crop, followed by soil texture, which typically influences the vigor gradient of the crop. This evidence suggests that VRT nitrogen fertilization may only partially mitigate the heterogeneity of production determined by such environmental factors.
In general, in favorable seasons, a low-input nitrogen fertilization driven by the NDVI direct model (Var-N-dir) may be the best choice in contexts where the vigor crop gradient is a direct reflection of the soil gradient. This approach provides nitrogen where the crop can use it and saves it where the crop is limited by other factors, limiting leaching in deep water, preserving the environment, and saving money. On the other hand, in unfavorable seasons, when decision support could be more helpful, the more advanced data-driven models (Agricolus and AgroSat) may be more effective in limiting losses and maintaining more homogeneous intra-field production. This evidence takes on greater importance for rainfed wheat given the increasing frequency of anomalous climate phenomena.