Sentinel 2-Based Nitrogen VRT Fertilization in Wheat: Comparison between Traditional and Simple Precision Practices

Abstract

This study aimed to compare standard and precision nitrogen (N) fertilization with variable rate technology (VRT) in winter wheat (Triticum aestivum L.) by combining data of NDVI (Normalized Difference Vegetation Index) from the Sentinel 2 satellite, grain yield mapping, and protein content. Precision N rates were calculated using simple linear models that can be easily used by non-specialists of precision agriculture, starting from widely available Sentinel 2 NDVI data. To remove the effects of not measured or unknown factors, the study area of about 14 hectares, located in Central Italy, was divided into 168 experimental units laid down in a randomized design. The first fertilization rate was the same for all experimental units (30 kg N ha⁻¹). The second one was varied according to three different treatments: 1) a standard rate of 120 kg N ha⁻¹ calculated by a common N balance; 2) a variable rate (60–120 kg N ha⁻¹) calculated from NDVI using a linear model where the maximum rate was equal to the standard rate (Var-N-low); 3) a variable rate (90–150 kg N ha⁻¹) calculated from NDVI using a linear model where the mean rate was equal to the standard rate (Var-N-high). Results indicate that differences between treatments in crop vegetation index, grain yield, and protein content were negligible and generally not significant. This evidence suggests that a low-N management approach, based on simple linear NDVI models and VRT, may considerably reduce the economic and environmental impact of N fertilization in winter wheat.

1. Introduction

In recent decades, population growth and the expanding demand of agricultural products for multiple purposes have constantly increased the environmental pressures on land and water resources. Worldwide concern of citizens and governments on environmental issues, together with the availability of improved and cost-effective methodologies and tools for spatial data acquisition, analysis, and modeling, have promoted the development of precision agriculture (PA) techniques. PA is based on innovative system approaches that use a combination of various technologies such as Geographic Information System (GIS), Global Navigation Satellite System (GNSS), computer modeling, Remote Sensing (RS), variable rate technology (VRT), yield mapping and advanced information processing for timely in-season and between season crop management [1].

Currently, a wide range of satellite data is available that varies in terms of acquisition cost, technique (active/passive), spatial resolution, spectral range, and viewing geometry [2]. A recent step forward was made thanks to the launches of Sentinel-2A (2015) and Sentinel-2B (2017) satellites by the European Space Agency (ESA), intended for Earth observation and monitoring of land surface variability. The Sentinel-2 satellites are equipped with a multispectral sensor (MSI) including 13 spectral bands, with a spatial resolution ranging from 10 m to 60 m, which provides relevant information for supporting precision agriculture [3]. Images provided by Sentinel-2 satellites are publicly available for free through the Copernicus Open Access Hub with 5 days temporal resolution averaging (2–3 days in mid-latitudes), in L1C or L2A processing levels, which make them attractive for time-series analysis and PA applications. The Level-2A product provides Bottom Of Atmosphere (BOA) reflectance images derived from the associated Level-1C product, which provides Top Of Atmosphere (TOA) reflectance images [4].

Satellite remote sensing of vegetation is mainly based on the green (495–570 nm) and red (620–750 nm) regions of the visible spectrum, the red-edge (680–730 nm), near and mid infrared bands (850–1700 nm). These bands are often combined with a range of algebraic formulas to obtain several Vegetation Indices (VIs) for assessing vegetation status and crop management (see e.g. references [5,6,7]). The Normalized Difference Vegetation Index (NDVI) [8], a normalized difference between reflectance of Red and Near-Infrared (NIR) spectral bands, is still the most used index worldwide because of its ease of calculation and interpretation. It varies between 0 and 1 in cropped areas, increasing with soil cover, LAI (Leaf Area Index), chlorophyll content, and plant N-nutritional status (see e.g. references [9,10,11]). The NDVI has been widely tested to assess wheat N-nutritional status and yield, with encouraging results [12,13,14].

N is the main nutrient supplied to most crops, including wheat, and may cause environmental impacts on near-surface and deep aquifers [15]. Increasing the N rate generally increases crop yield since it increases the grain number and size [16]. However, increasing the N rates reduces the N uptake efficiency and increases the amount of residual N in the soil which is exposed to leaching risks [9,17]. Countless studies are available on wheat N nutrition and several of the recent ones have dealt with precision N fertilization, but most of them propose quite complex approaches that are not widely adoptable by non-specialists of PA. For example, Bourdin et al. [18] propose a complex model starting from LAI, estimated by remote sensing, and yields from previous years; Basso et al. [19] use a crop simulated model (SALUS) based on weather and yields from previous years. Another approach consists in assessing spatial variability to delineate homogeneous sub-field areas overlapping various thematic spatial maps (such as yield, soil properties) or applying a multivariate geostatistical approach of the factors that affect yield and grain quality [20]. In this context Song et al. [21] delineated management zones on the basis of soil, yield data and remote sensing information derived from Quickbird imagery. Besides these methods, to improve the use of N prescription maps in PA, various simplified approaches based on NDVI have also been developed through user-friendly web interfaces (e.g. CropSAT [22], Agrosat [23], OneSoil [24]) but these approaches can be applied only in specific areas or lack either quantitative analysis or validation.

While traditional flat-rate spreaders typically tend to over- or under-apply fertilizers, with VRT spreaders it is possible to modify the distribution rate according to on-the-go proximal sensors measuring, in real-time, plant properties or prescription maps derived through more or less advanced models including various spatial factors (e.g.: vegetation indices, soil variables, yield maps, crop nutrition status) [25]. However, the majority of precision fertilization approaches are generally map-based, because on-the-go sensors are too expensive, not sufficiently accurate, or not available [26]. VRT spreaders are subject to errors depending on the technology and the intrinsic characteristics of the fertilizer [27].

Grain yield mapping aims at providing or increasing knowledge about the spatial heterogeneity of yields. This technology is based on the recording of georeferenced yield data thanks to the integration of different subsystems mounted on the combines [28]: (i) harvest quantity measurement sensors (mass or volume); (ii) GNSS location sensors; (iii) reference area measurement sensors (working width, speed, time); (iv) data recording and processing units. To generate yield maps, the harvested quantities are georeferenced and associated with the corresponding reference area which is calculated by multiplying the working width by the area length (working speed * time interval). Yield mapping accuracy is influenced by many factors including flow sensor calibration, combined speed changes, grain flow variations, grain moisture, and the level of smoothing applied for mapping [29].

While many studies combine VIs (from various remote or proximal sensors) and VRT crop fertilization, only a small quantity of research integrates yield mapping sensors, and, to our knowledge, probably no study has ever combined VIs, VRT fertilization, and yield mapping in a PA case study and, in particular, to assess the effects of different variable N-rate treatments on winter wheat.

In this framework, this study was aimed at comparing two VRT N fertilization treatments (based on Sentinel 2) versus a standard flat N rate in terms of crop NDVI trend, grain yield, and protein content. In order to promote the adoption of PA techniques among farmers and non-specialists in PA, a simplified approach, based on free and open-source software, the widely used NDVI, and an easily-applicable linear model, was applied to calculate the VRT N fertilization rates.

2. Materials and Methods

2.1. Crop Management and Experimental Treatments

The experiment was carried out in the cropping season 2017-2018, on a 14 ha plain field of the middle Tiber valley, near Deruta, Umbria, Italy (170 m a.s.l., 42°95’07’’ N, 12°38’18’’ E), provided by the Foundation for Agricultural Education of Perugia (Figure 1). The soil was loam with increasing sand content from the west to the east side. The climate is Mediterranean, characterized by a dry season between May and September and a cold and rainy season from October-November to March-April. The cropping season 2017–2018 was unusually rainy in December and March, while temperatures were generally higher than the poly-annual trend except for a very cold end of February (Figure 2).

The experimental crop was a rainfed winter wheat (Triticum aestivum L., cv. PRR58) sown on 10 November 2017 at a nominal density of 450 viable seeds m^-2. The previous crop had been pea (Pisum sativum Asch et Gr). The crop was managed according to ordinary practices, while weeds and diseases were controlled chemically. The N rate was split in two application times in order to increase N uptake efficiency, support the formation of yield components and limit N leaching by fall-spring rainfall. The first N application (as urea) occurred on 18 January 2018 with 30 kg N ha⁻¹, while the second N fertilization (as urea), occurred on 26 March 2018 and was managed according to three experimental treatments: 1) a standard rate of 120 kg N ha⁻¹ (Flat-N) derived by an N balance approach (the relatively high rate is justified by the very rainy winter); 2) a variable rate of 60 to 120 kg N ha⁻¹, based on NDVI, where the maximum rate was equal to the standard rate (Var-N-low); 3) a variable rate of 90 to 150 kg N ha⁻¹, based on NDVI, where the medium rate was equal to the standard rate (Var-N-high). An inverse linear relationship between NDVI and VRT N-rates was adopted in Var-N-low and Var-N-high on the assumption that NDVI and other correlated VIs (e.g. NIR/Red simple ratio), before the second N fertilization (7–8 Feekes’ stage), are directly related to crop N nutritional status [10,13,30]. Thus, the Var-N-high was defined with the aim of keeping the average N-rate equal to the flat N-rate while optimizing the N fertilization according to the supposed N nutritional status. The range adopted in Var-N-high represents a variation of ±25% of the flat N rate, and a ±20% of the total N supplied (including the 30 kg N ha⁻¹ of the first N application), which was considered as being appropriate to see appreciable effects on wheat yield. The Var-N-low was defined with the aim of testing a sensible N reduction, the results of which were very important in those areas where a more eco-compatible approach is required. So, the flat N-rate was considered as the maximum N level and the −50% of the flat N rate (−40% of the total rate, equal to 90 Kg/ha including the first N application) was adopted as a lowest reference level based on the evidence from a previous study in the same location [16].

To account for the possible not measured or unknown factors, the 14 ha study area was divided into 168 plots of about 700 square meters each (35 m long, 21 m wide) grouped in 28 zones. The three treatments were laid down according to a randomized design with 2 replicates per treatment in each zone for a total of 56 replicates per treatment in the whole study area (Figure 1). All the precision on-field operations were performed using a tractor equipped with a Topcon GNSS automatic guide device connected to a Real Time Kinematic (RTK) network. The variable rate treatments were performed using a Sulky 40+ (ECONOV) VRT fertilizer spreader connected through ISOBUS (a widely used software protocol complaint to ISO 11783 standard) to the Topcon system console. The console, after each VRT fertilization, releases the output of each fertilization map showing the estimated distributed rates. This output was compared to the prescription map to investigate the accuracy of the N distribution.

2.2. Calculation of Rates for VRT Treatments

A Level 2A Sentinel-2 image collected on 22 March 2018 (i.e., four days before the second fertilization), georeferenced according to WGS84-UTM32, was downloaded from The Copernicus Open Access Hub. NDVI, calculated from bands 4 and 8 using QGIS Raster calculator (Figure 3a), was used for the N rates calculation of the two variable rate treatments. In both cases, a linear relationship was imposed between the average NDVI value calculated for all experimental units and fertilizer-N rates where the 5° percentile of NDVI value corresponded to the minimum fertilizer-N rates (60 or 90 kg N ha⁻¹) and the 95° percentile of NDVI value corresponded to the maximum fertilizer-N rate (120 or 150 kg N ha⁻¹) (Figure 3b). All the data processing and analysis were carried out using QGIS software, version 2.18.12 64bit [31] and MS Excel 2016. Average NDVI values were calculated using the SAGA grid statistics for the polygons algorithm included in the QGIS processing framework. To simplify this step, pixels overlapping the border of the experimental plots were not excluded from the calculation. To assess the effect of these pixels, average NDVI values were compared with those obtained while considering only non-overlapping pixels.

2.3. Determination of NDVI

To monitor and compare the NDVI of experimental treatments, all relevant level 2A Sentinel-2 images with no cloud cover on the study area, were collected from 22 March to the beginning of the senescence. Average NDVI values were calculated for each plot using the SAGA grid statistics for the polygons algorithm included in the QGIS processing framework. Six to nine S2 pixels fell within each experimental unit. To highlight possible differences between the effects of the experimental treatments while considering different crop vigor status before the second N fertilization, the NDVI time series analysis was performed by classifying plots according to three NDVI classes (NDVI ≤ 0.68; 0.68 < NDVI ≤ 0.79; NDVI > 0.79), as revealed by the S2 image from 22 March.

2.4. Determination of Grain Yield and Protein Content

Harvest was carried out on 26 June 2018 using a combined harvester Claas Lexion 780–740 equipped with a Topcon YieldTrakk system (processing data from the optical sensor measuring grain mass flow and moisture sensors), which produced a georeferenced yield map as an ESRI polygon shapefile. The combine harvester had a cutting width of 7.50 m and was equipped with a tilt sensor to correct the effect of slope on the sensor readings. An initial on-field calibration was performed on the combine to adjust for the actual working width and measure the unit weight of grain, which was used by the system to convert the measured mass flow (l/s) to Mg.

The output shapefile containing the georeferenced yield data was intersected with the experimental units to calculate the average yield expressed in Mg ha⁻¹ which was used for the subsequent analysis. To produce an easy-to-read yield map, the shapefile was converted to a 5 m resolution raster. Subsequent gaussian kernel smoothing, based on a 50 m search radius, produced the final yield map.

Protein content was measured on 20 June 2018 (i.e., six days before final harvest). To reduce the number of samples, the sampling scheme was defined according to the three classes of NDVI previously defined for a total of 36 field samples, then merged two by two to obtain 18 lab samples, which were analyzed according to the official Kjeldahl method, a widely used chemical procedure for the quantitative determination of protein content in food, feed, feed ingredients and beverages [32].

3. Results

The NDVI values on 22 March show within-field differences mainly associated with the west-east textural gradient (Figure 3a). Average NDVI values including overlapping pixels with plot borders differ only slightly from those obtained while excluding these pixels (0.015 in average absolute value, 0.011 St. Dev.). According to the linear model used for N-rates calculation, these differences correspond only to 2.6 kg/ha (1.8 St. Dev.) for Var-N-low treatment and 2.6 kg/ha (2.1 St. Dev.) for the Var-N-high treatment.

A comparison between the prescription map and the Sulky output (Figure 4a) shows that the correlation between the prescribed rate and the rate actually supplied by the VRT machine within each plot was very good (R² = 0.91) (Figure 4b), with only a little bias which resulted in smoothing the extremes, raising the lower rates and lowering the higher ones (

Table 1. Amounts of N prescribed and supplied with the second N application in the three treatments: Flat-N (standard rate of 120 kg N ha⁻¹), Var-N-low (variable rate from 60 to 120 kg N ha⁻¹), Var-N-high (variable rate from 90 to 150 kg N ha⁻¹).

Treatments	Average N prescribed (Kg N ha−1)	Average N supply (Kg N ha−1)	Average Δ (%)	Δ St. Dev. (Kg N ha−1)
Flat-N	120.0	116.5	−1.3%	0.03
Var-N-low	90.2	95.4	5.4%	0.06
Var-N-high	120.5	117.0	−3.0%	0.05

).

After the second fertilization, the trend of NDVI, considering the three vigor classes previously defined, was not affected by fertilization treatments (Figure 5). In the portion of the field where the NDVI was already high on 22 March (class III—about 0.85), all treatments showed a further moderate increase up to over 0.9, while in the portion of the field where NDVI was low (class I—just over 0.6), all treatments showed an increase up to 0.85 in about one month, and nearly 0.9 later on. As expected, an intermediate trend can be observed for treatments with NDVI falling within class II. Thus, whatever the crop vigor status was on March 22nd, all treatments reached a high NDVI one month after the second N application (Figure 5).

Grain Yield and Quality

The yield map used for the quantitative analysis is shown in Figure 6a, while the generalized version used only for visual analysis purposes is shown in Figure 6b. The three treatments did not differ significantly for total yield (

Table 2. Grain yield and protein content for the three N fertilization treatments: Flat-N (standard rate of 120 kg N ha⁻¹), Var-N-low (variable rate from 60 to 120 kg N ha⁻¹), Var-N-high (variable rate from 90 to 150 kg N ha⁻¹).

Treatment	Average Yield (Mg ha−1)	St. Dev. Yield (Mg ha−1)	Average Protein Content (%)	St. Dev. Protein Content (%)
Flat-N	6.74	0.36	9.4	0.31
Var-N-low	6.73	0.39	8.9	0.37
Var-N-high	6.76	0.38	9.2	0.48

) and no relationship was found between the N rate and yield (R² = 0.087). Similarly, no correlation was found between N treatments and protein content (R² = 0.001). However, the slight difference in protein content observed between Flat-N and Var-N-low (Table 2) was statistically significant (p = 0.02) even though it was not agronomically relevant. Finally, grain yield was weakly correlated to NDVI at any time of NDVI measurements (

Table 3. Correlation between yield and NDVI at any time of NDVI measurement plotted over all of the 168 plots.

Date	R2
22 Mar	0.29
6 Apr	0.26
21 Apr	0.24
29 Apr	0.19
11 May	0.19
10 May	0.18
26 May	0.27
31 May	0.26

).

4. Discussion

Our study investigated the differences between a flat N fertilization rate and two variable rates defined through a simplified method based on Sentinel-2 NDVI collected a few days before the second fertilization. One of the two variable rates was aimed at drastically reducing the overall N input (Var-N-low), the other one was aimed at maintaining the same overall input (Var-N-high) as in the flat rate (Flat-N), i.e., 120 kg N ha⁻¹ while optimizing the N fertilization according to the supposed N nutritional status. Our results show that a VRT approach with a lower overall N rate may be more efficient, giving same grain yield and quality (Table 2) with a lower N input (Table 1). This result is consistent with the evidence of Raun et al. [33] which reports that the VRT method improves the NUE (Nitrogen Use Efficiency) by 15% compared to a flat rate. In vulnerable contexts, as in the study area, the N-rate reduction results environmentally and economically very relevant since it could reduce water pollution (still a critical issue in Umbria and all over the world) and decrease expenditures on N. However, differently from our case study, a Var-N-low treatment could determine a lower protein grain content. In this regard, the economic trade-off between lower N costs and reduced yield value due to lower grain protein content could be calculated by the farmers using their local costs and prices. Indeed, the treatment Var-N-low differed from the other two treatments by only about 20 kg N ha⁻¹ on average, which was supposed to be not so important compared to the total amount available (i.e., about 160 vs 180 kg ha⁻¹), accounting for the residual N left by the previous pea crop (likely around 30 kg N ha⁻¹) and the amount of N supplied with the first mineral N application (30 kg N ha⁻¹). It is worth noting that about 180 kg ha⁻¹ of available N was in line with the usual practice adopted by farmers in the Tiber valley of Umbria and adequate, if not limiting, in view of the high yielding potential of this area (6–8 tons grains ha⁻¹). Nonetheless, the rainy fall-spring of 2017–2018 caused a relevant N loss from the soil, thus, in this case, a difference of 20 kg N ha⁻¹ was expected to make a difference. Several reasons can be addressed to explain the lack of difference among treatments. The main one is that crop development was not associated to N availability but to the soil textural gradient, with NDVI values tendentially as lower as higher appeared the sand content in the soil. In such a condition, and considering the rainy season, forcing the N rate in plots with low NDVI (i.e., with higher sand content) resulted in losing N (mainly by leaching) without affecting yield. As reported by Asseng & Turner [34], NUE is mainly influenced by soil and rain. Thus, in our case, the N rate increase did not translate in a NUE increase as this was probably limited by the other factors. This would justify that, in some cases, the N rate should be directly related to the NDVI [35], i.e., increased where NDVI is higher and decreased where NDVI is lower, on the assumption that lower NDVI values could prove that other soil characteristics, besides N availability, are not suitable for allowing high yields. For example, a soil with low water retention is not suitable in case of dry grain filling periods (as it is in Mediterranean areas) due to water shortage and thus lower yields whatever the N availability [36,37,38]. The similar time course of NDVI observed in all the three treatments, distinguishing three different vigor status before the second fertilization (Figure 5), also proves that N was not the limiting factor, or at least, the 20 kg N ha⁻¹ of difference among treatments was not relevant for crop NDVI.

The weak correlation between grain yield and NDVIs at any date of monitoring (R² always below 0.3) appears in contrast with the evidence of Sultana et al. [12] and Liu et al. [39] which show high correlation between NDVI and yield, especially for NDVI recorded during the milky-grain vegetative stage. However, the NDVI recorded in our experiment from 21 April onwards was quite high for all treatments, whereas it is known that the correlation with yield stands only when a wide range of NDVI is considered [16].

As discussed above, the textural gradient likely affected yield through N and water availability, more than the N fertilization treatments. This assumption is supported by the case study of Basso et al. [19] who observed a high correlation between yield and pedological conditions, in particular soil water capacity. Also, Zhao et al. [40] reported that fertilization, as well as the availability of water, significantly affected the content of wheat proteins, even though the latter was more influent, especially for those cultivars characterized by an intermediate protein content. Probably, further splitting the rate, with a third application at the end of shooting, would have prevented some N loss and increased the grain protein content [41], which was quite low compared to the standard of our cultivar. Again, the low protein content can be ascribed to the rainy season, as it was widespread for the harvest of 2018 in Umbria and all-over Central Italy.

Concerning the prescription map used in the study, including the overlapping pixels with plot borders in average NDVI calculation simplifies the GIS procedures without generating agronomically meaningful N-rates differences. Because of the high variability within field zones and among plots, the map could be considered as a kind of stress test for the VRT device. The very little differences between prescribed and distributed rates (Table 2) suggest that this technology is suitable for precision fertilization even with highly detailed prescription maps. The highlighted differences are clearly due to small plot size and random allotment of treatments. As a result, the VRT machine was subjected to abrupt rate variations due to the short time available to adapt the opening of the distribution valve.

The simple linear approach to calculate a N prescription map proposed in this research, even though requiring average GIS skills could be effective for VRT adoption by non-specialist farmers. Coherently with other similar experiences [22,23,24], to improve the overall usability of the method by non-specialists in PA and allow further validations and tuning, all the RS data management and GIS calculation could be implemented in an user-friendly web-GIS application where the user could upload (or digitize on aerial data) his own fields, choose the S2 NDVI reference date, and decide the most suitable approach for defining the N rate. NDVI S2 data could be conveniently accessed through the Sentinel Hub platform API [42].

Concerning the integration of NDVI from Sentinel 2, VRT fertilization, and yield mapping in PA, our experience confirms that these technologies, thanks also to both the web-based applications for calculating satellite-based prescription maps [22,23,24] and the user-friendly interfaces of system consoles, could be ordinarily used even in small or medium-sized farms. However, a preliminary cost-benefit analysis and a start-up training and support by a PA expert is often necessary.

5. Conclusions

Our results regarding crop vegetation index, grain yield, and protein content indicate that the adoption of a low-N management approach, based on simple linear models and VRT, may considerably reduce the economic and environmental impact of nitrogen fertilization in winter wheat. However, the rainfed nature of winter wheat in Mediterranean environments may cause unpredictable yield and quality variations depending on climatic trends and soil properties due to their effects on both soil nitrogen and water availability. In this view, VRT nitrogen fertilization can only partially mitigate the heterogeneity of production determined by such environmental factors. The alternative approach of providing a nitrogen supply proportional to the crop NDVI deserves to be considered when factors other than N fertilizer rate come into play, as it is with sandy soils where NDVI and yield may be limited by low N and water retention. Despite the known limitations of predictions and prescriptions based on remoted sensed vegetation indices, their use provides relevant information about within-field and between-fields variability. This information can support the implementation of crop fertilization management approaches based on GNSS/RTK and VRT technologies to replace the traditional flat-N rate, which provides results that are neither economically nor environmentally sustainable.