Optimizing Top Dressing Nitrogen Fertilization Using VEN µ S and Sentinel-2 L1 Data

: Environmental and economic constraints are forcing farmers to be more precise in the rates and timing of nitrogen (N) fertilizer application to wheat. In practice, N is frequently applied without knowledge of the precise amount needed or the likelihood of signiﬁcant protein enhancement. The objective of this study was to help farmers optimize top dress N application by adopting the use of within-ﬁeld reference N strips. We developed an assisting app on the Google Earth Engine (GEE) platform to map the spatial variability of four different vegetation indices (VIs) in each ﬁeld by calculating the mean VI, masking extreme values (three standard deviations, 3 σ ) of each ﬁeld, and presenting the anomaly as a deviation of ± σ and ± 2 σ or deviation of percentage. VIs based on red-edge bands (REIP, NDRE, ICCI) were very useful for the detection of wheat above ground N uptake and in-ﬁeld anomalies. VEN µ S high temporal and spatial resolutions provide advantages over Sentinel-2 in monitoring agricultural ﬁelds during the growing season, representing the within-ﬁeld variations and for decision making, but the spatial coverage and accessibility of Sentinel-2 data are much better. Sentinel-2 data is already available on the GEE platform and was found to be of much help for the farmers in optimizing topdressing N application to wheat, applying it only where it will increase grain yield and/or grain quality. Therefore, the GEE anomaly app can be used for top-N dressing application decisions. Nevertheless, there are some issues that must be tested, and more research is required. To conclude, satellite images can be used in the GEE platform for anomaly detection, rendering results within a few seconds. The ability to use L1 VEN µ S or Sentinel-2 data without atmospheric correction through GEE opens the opportunity to use these data for several applications by farmers and others.


Introduction
Wheat is one of the key cereal crops grown worldwide, providing the primary caloric and nutritional source for millions of people. The climatic conditions of the Israeli fields, typical of Mediterranean regions, are highly variable. These conditions cause substantial variability in wheat (Triticum aestivum L.) grain production and quality, a matter of great concern for both producers and bakers. Hard spring wheat (HSW) is the most important crop in Israel. HSW has a high protein concentration (>11.5%), which is also associated with greater gluten content thus making HSW superior for milling and bread-making purposes. In Mediterranean rain-fed conditions, the amount of rainfall and its distribution during the growing season have a marked influence on wheat yield and quality [1][2][3][4][5][6]. Drought restricts yields and test weight but can increase grain protein, a result of shriveled kernels. In years with good moisture, wheat yields are less restricted, but this can lead to low grain protein. In conditions of abundant moisture, low grain quality has been attributed to a dilution of grain protein by larger grain biomass [1,7]. Temperature, rainfall, moisture stress before anthesis, solar radiation during grain filling, soil nitrogen (N), and rate and timing of supplemental N application can influence wheat protein composition through sink/source effects on N during early growth stages [1,5,[7][8][9][10].
Efficient use of N fertilizer is important for economical wheat production. Environmental and economic constraints are forcing farmers to be precise in the rates and timing of N fertilizer applications to wheat. Insufficient N reduces wheat yield and profit, while excessive N results in wheat that is susceptible to water deficiency, disease, and lodging, with consequently reduced quantity and quality of yield. N fertilizer is typically applied in Israel in one of the following options: (i) All N fertilizer is applied in the fall before sowing; (ii) the same process but followed by possible N topdressing application according to growing conditions; (iii) Some N fertilizer is applied in the fall, followed by a mid-winter or early spring top dressing. Although pre-sowing N fertilizer applications decrease the potential for N deficiencies in the early stages of growth, the presence of plant-available residual soil N from the previous season is at risk of being exported to the environment (leaching or denitrification) and/or the creation of growth problems. Excessive pre-sowing applications of N encourage vegetative growth and, therefore, the crop can utilize too much soil water during that growth stage [11], leaving insufficient water during the grain filling stage [2,12]. Late-season N applications, combined with a better temporal distribution of N during the wheat cycle, significantly improved grain protein over yield and enhanced the bread-making quality of the flour [7,8,13,14]. In contrast, late-season moisture stress is essential for increasing protein up to acceptable levels, as moisture stress was found to increase protein content and reduce grain weight [9,[15][16][17][18].
Acquiring information about crops from remote platforms, such as unmanned aerial vehicles, aircraft, or satellites, has great potential for decision-making facilitation in season fertilizer management practices [19]. At present, most application algorithms use within-field reference N strips, regardless of the sensor [20][21][22][23][24]. Bonfil [25] showed that the use of a proximal-sensing technique allows for rapid and accurate crop monitoring and yield estimation but limits future use of remote sensing techniques as well. Agricultural monitoring is conducted typically with vegetation indices (VIs) [26][27][28][29]. Remote and proximal sensing were previously used to forecast crop grain yield and protein content [30][31][32][33][34][35][36]. Spatial resolution is an important factor determining the accuracy of crop yield and quality assessments, with satellites having sensors that range from high to coarse spatial resolutions.
Google Earth Engine (GEE) is a planetary-scale geospatial analysis cloud-based platform that processes a large number of freely available scientific datasets such as Sentinel-2 (L1 and L2). GEE leverages its cloud computing service for analysis capabilities of time series data, which are otherwise heavy consumers of time and computation resources, and it provides predefined processes and functions in two APIs (JavaScript and Python). GEE also helps researchers to easily disseminate their products to other users and researchers, which enabled us to provide code and a web application, for the benefit of farmers and the science community. GEE apps are designed to provide non-expert stakeholders (farmers) with different tools, e.g., for optimizing crop management [37,38]. These may be used for different agricultural purposes. For example, while the coarse-resolution advanced very high-resolution radiometer (AVHRR;~1 km) enables a stable and accurate yield estimation at the county-to-regional levels [39], high spatial resolution data from sensors on satellites such as IKONOS, VENµS, and Sentinel-2 (~4-10 m) may be used for precision farming applications [40][41][42][43][44]. Even though wheat is the main field crop sown globally, there is still a lack of knowledge regarding the optimal practice that should be exerted to properly manage this crop. VENµS and Sentinel-2 satellites with their special features (high resolutions: spatial, super-spectral, and temporal) will hopefully create new abilities in detecting within-field reference N strips and improve wheat crop management through continuous monitoring (in terms of time and space) of wheat fields.
The objectives of this study are (i) to develop algorithms for within-field variation calculation (SD or percentage) on the GEE platform using Sentinel-2 and VENµS data; (ii) display the variations with respect to the reference strips by GEE apps; (iii) to evaluate the outcomes of farmers practice in using these algorithms. In the following, we describe the wheat fields, reference strips (Section 2.1), ground observations (i.e, biophysical properties of crop, Section 2.2), Sentinel-2-and VENµS-based Vegetation Indices (Section 2.3), and the GEE platform (Section 2.4). In the Results Section, we explore the correlations between VIs and biophysical properties of crops (Section 3.1) and use them to interpret the relationship between the VI anomaly maps (presented in the GEE app, Section 3.2), and the grain yield and quality maps (presented in Section 3.3). Finally, we discuss the potential, limitations, and existing knowledge gaps in Section 4.

Wheat Fields and Reference Strips
The data for this study were collected during three growing seasons 2018-2020 at individual producers' fields across Israel, one or more fields in each farm ( Figure 1). Farmers (with their extension advisers) established reference strips in many ways in the first season. At the second and third seasons, farmers were asked to perform the following tasks: (1) keeping a minimum strip size of 30 × 30 m (to include few pure VENµS pixels), preferably 60 × 60 m or larger; (2) keeping a buffer of at least 10 m from the field edge (to avoid edge problems, e.g., double stand, weeds, roads) if the strip does not cross the whole field; (3) applying higher and/or lower amount (±50 kg N ha −1 ) base and/or top N vs. the commercial (standard) application; (4) maintaining a north-south or west-east positioning, which are preferred, as the satellite data are disseminated in UTM grid; (5) marking clearly the strip edges (to ease identification).

Ground Observations
Ground observations including extensive field observations of crop progress were conducted along the growing season. Samples of biomass, N uptake, and leaf area index (LAI) were collected from commercial and experimental fields along the main development stages of active vegetative growth, spanning from two leaves (Zadoks 12) up to heading (Zadoks 59), reflecting large differences in canopy height and vegetation cover. Four subsamples (0.3 m 2 ) were taken from an area of 10 × 10 m in each sampled location (treatment*field*day). A total of 121 mean values of biophysical variables (biomass, N content) were calculated based on 960 subsamples, and 53 mean LAI were calculated based on 229 records. The samples were dried for 48 h at 70 • C and weighted (dry weight, DW), and the above-ground plant was ground and mixed. N concentration was obtained using NIRS (FOSS 6500) following the micro-Kjeldahl method [45]. N content was calculated by multiplying the DW by N concentration. These data were used to assess the relationship between VIs (see Section 2.2) and biophysical properties of crops in Section 3.1.
Yield samples were collected from fields that were not harvested for hay or silage production. In total, 685 yield samples were collected from Gilat experiments and 187 from commercial fields. Yield maps were available for most fields, and data were cleaned up prior to analysis. Saad fields were harvested by a combine equipped with a grain protein sensor (CropScan 3300H, Next Instruments, Condell Park, NSW, Australia). More than 300 samples (each season) from few fields were manually sampled during harvest to test the accuracy of the grain protein sensor.

Sentinel-2-and VENµS-Based Vegetation Indices
Sentinel-2 (A and B) level 1C and VENµS level 1 (L1) top-of-atmosphere (TOA) reflectance were processed to assess variation in wheat biomass, grain yield, and grain quality. No atmospheric correction was used, enabling to run the algorithm as an application in the GEE platform. VENµS Super Spectral Camera and Sentinel-2 MultiSpectral Instrument have many overlapping spectral bands with small differences [46]. Overall, 42 VIs (Supplementary Materials) were computed, as a potential tool for monitoring the crop's nutritional contents by quantitatively retrieving the nutrient content in the crop. The bands that were used for the four selected VIs (NDVI, REIP, NDRE, ICCI) and the equations are listed in Tables 1 and 2. The spatial resolution for all VENµS VIs is 5 m, while for the Sentinel-2 Vis, it is 20 m except for NDVI, which is 10 m. The temporal resolution of VENµS (revisit interval of two days) is also better than that of Sentinel-2 (five days). The tradeoff is the global coverage of Sentinel-2 vs. local coverage of VENµS. The higher temporal resolution of VENµS is translated into much more data along the season, and (more important) more "clear sky" images. Figure 2 shows the NDVI phenology of one field, emphasizing less high cloud mask conditions by VENµS (red circles). Table 1. Vegetation indices used in the study: normalized difference vegetation index (NDVI), Normalized difference red-edge (NDRE), red-edge inflection point (REIP), and red-edge canopy chlorophyll index (ICCI). Band numbers are indicated in Table 2.

Google Earth Engine (GEE)
The present study provides a GEE web app that enables optimization of yield by providing the farmer with the ability to monitor the spatial heterogeneity of VIs anomaly (with respect to the mean) while using a reference strip with known N application to interpret the retrieved anomaly and apply N based on this information. VENµS images and selected polygons of wheat fields in Israel were uploaded to the GEE platform. Figure 3 describes in orange the users' inputs before calculating the anomalies and, in green, the anomaly calculations performed (in the background) on each polygon (field). Figure 4 presents the results (i.e., the NDVI anomaly) for 15 Jan 2019, displayed in a browser as a deviation of ±σ and ±2σ (or percentage) from the mean using predefined symbology. Five wheat fields in Saad and references strips polygons are displayed, while 3σ extreme values are masked.

Correlation between VIs and Biophysical Properties of Crop
Non-linear equations provided higher correlations (r) between the wheat parameters and VIs ( Figure 5). Wheat biophysical variables (biomass, nitrogen content, and LAI) were regressed against the best three VIs (REIP, NDRE, ICCI) and NDVI from the pixel of the sampling location (samples cannot represent spatial distribution) ( Table 3).  The correlations for the averaged (per field*day) samples were higher than those for the individual sampling points for both VENµS and Sentinel-2. REIP, ICCI, and NDRE exhibit a high correlation with the biophysical variables. For these parameters, the VENµS higher resolution (5 m pixel) does not exhibit significant improvement over Sentinel-2.

Correlation between VIs and Biophysical Properties of Crop
Non-linear equations provided higher correlations (r) between the wheat parameters and VIs ( Figure 5). Wheat biophysical variables (biomass, nitrogen content, and LAI) were regressed against the best three VIs (REIP, NDRE, ICCI) and NDVI from the pixel of the sampling location (samples cannot represent spatial distribution) ( Table 3).  The correlations for the averaged (per field*day) samples were higher than those for the individual sampling points for both VENµS and Sentinel-2. REIP, ICCI, and NDRE exhibit a high correlation with the biophysical variables. For these parameters, the VENµS higher resolution (5 m pixel) does not exhibit significant improvement over Sentinel-2. The correlations with canopy N concentration were low for most indices. The correlations were much higher for canopy N content, similar to those for biomass, as was also observed over rice [50]. The correlations were highest for REIP, r = 0.81 (VENµS) and 0.87 (Sentinel-2), and steady along the three growing seasons ( Figure 5). VENµS VI's values were in high correlation with Sentinel-2 VI's ( Figure 6), although several sampling fields/dates can be considered outliers probably because of abnormal atmospheric conditions.  Figure 7 shows VENµS and Sentinel-2 NDVI and REIP anomalies for five fields in Saad on 15 January 2019, as seen in the GEE anomaly app. The NDVI within field patterns for VENµS and Sentinel-2 were similar. VENµS-based REIP exhibited much more variation than Sentinel-2, probably due to the coarser spatial resolution of Sentinel-2.

Evaluation of Grain Yield and Quality
Reference strips should be used for topdressing nitrogen application (i.e., if a significant difference is identified between the reference strip and the whole field, a topdressing is recommended such as that performed with proximal sensing) and evaluated with respect to the grain quantity and quality. Figure 8 shows row data of grain yield and grain protein from the five fields that are presented in Figure 7. In field 101N, herbicide drift from the nearby field (west to east) decreased grain yield and increased grain protein; this effect is shown nicely by the NDVI-based anomaly app from 3 weeks after herbicide application. In field 105S, the anomaly app ( Figure 7) shows a significant difference between the west and east parts regarding soil management (previous crop), which is reflected in the yield map ( Figure 8a). Field 105N exhibits a higher yield in the west south area, similar to the NDVI-based anomaly app (Figure 7). Yield reduction is seen at two small areas, in the east and in the south of field 106S (probably weeds and water, respectively). Yield production is more uniform than grain protein in field 106N (Figure 8a). The variation shown by the anomaly app (Figure 7) for this field can be attributed to the lower grain protein content in the south of field 106N with respect to the north of that field (Figure 8b, control and reference, respectively), which is negatively correlated to yield and NDVI.
Farmers and extension advisors used the first version of the GEE app during the 2020 season. Time series (along the growing season) of red-edge-based VIs were very useful for the detection of infield anomalies. Figure 9 shows time series from one field in Beit Guvrin (Field 109) to demonstrate the anomaly detection by the reference strips. N treatments (base and top) are described in Table 4. On average, during tillering, the values of the high-reference strip (refH) were similar to the commercial (comm) management but were higher than the neighboring area. The values of the VIs in the lower-reference strip (refL) were significantly lower, starting about 30 days after emergence, thus forcing the farmer to decide whether to apply topdressed nitrogen. In this example, the eastern part of the field received the topdressed fertilization (topN), while the areas near the reference strips (and the strips) were left without this topN application (comm). The topdressed N indeed had a positive effect on the crop, reaching higher values of VIs (since application) and grain yield (Figures 9 and 10, Table 4). More important, grain quality was improved by the topN application, showing an increase in protein and gluten content, reaching the quality standard (protein >11% and wet gluten > 25%, Table 4). Split nitrogen application (commercial + topN) yielded the same grain quantity as the base application (refH) but with higher quality, emphasizing the importance of splitting nitrogen application in this growing region (Table 4). This advantage of the split application is significantly shown by all VIs, both VENµS, and Sentinel-2.   Figure 10). Average VI values are plotted for commercial (comm), reference strips (refH, refL), and topdressed N (topN) managements. VENµS operation caused missing period January-February data.  Figure 10. Wheat grain yield (raw data, t ha −1 ) map of Beit Guvrin, field 109, 2020 season. N management is described in Table 4. Figure 11 shows that the refH (triangles) and refL (circles) values are only slightly higher/lower (up to 5%) than the control (except ICCI-yellow), whereas topN (squares) shows a much higher difference by all VIs. Since refL did not receive the topN (after 16 January), all VIs exhibit the deficiency (circles are always lowest). ICCI showed the highest difference between managements (for both VENµS and Sentinel-2). Figure 11. Reference strip or topdressed N managements: VENµS-and Sentinel-2-based NDVI, NDRE, REIP, and ICCI difference from the commercial management along the growing season. Beit Guvrin 2020 (Field 109). REIP difference was calculated based on REIP-700. VENµS operation caused missing period January-February data.

Discussion and Conclusions
Manivasagam et al. [46] developed transformation coefficients from VENµS to Sentinel-2 surface reflectance for their overlapping spectral bands (i.e., visible, red edge, and near infrared), enabling the images from Sentinel-2 and VENµS to be combined into one time series that facilitates continuous, temporally dense vegetation monitoring. However, the VI values per se are not important for anomaly detection and/or the difference in the reference strip identification. Although retrieving absolute values is preferable, the high correlation between the sensors-based VIs ( Figure 6) enabled using an algorithm that is based on VIs regardless of the sensor and detecting anomalies even without atmospheric corrections. The correlations with NDVI were lower than all red-edge-based VIs, for biomass canopy N content ( Table 3). The main advantage of using NDVI is the higher spatial resolution in many satellites' sensors (including Sentinel-2).
The potential to retrieve biophysical parameters using VENµS and Sentinel-2 L1 images was demonstrated. Although the correlations between satellite-based VIs and wheat parameters were lower than those from proximal-sensing-based (multi-and hyperspectral) VIs [49,[51][52][53][54][55], only a small reduction is shown for VI correlations with nitrogen uptake. Moreover, crop phenology (growth stage) affects these correlations [51,54,56], lowering the correlations where samples cover many growth stages in different environments, as was the case in this study. REIP exhibited the highest correlations with the biophysical variables. WDRVI) [53] reached high correlations but lower than REIP (data not shown). The fits based on proximal sensing [25] were better than those based on satellites data (Table 3). This can be attributed to timing, as proximal sensing was taken at the same time as the samples, while satellite remote sensing data were retrieved mostly not at the same sampling day, even with a long period between them ( Figure 6).
VENµS spatial, spectral, and temporal resolutions are higher than other available satellites. The combination of 5 m pixels for the red-edge bands, and the two days revisit time, provides an advantage over Sentinel-2 data. The tradeoff is the specific small regions that are covered by VENµS versus the global coverage of Sentinel-2. The ability to use L1 VENµS or Sentinel-2 data through GEE opens the opportunity to use these data for several applications by farmers and others. An obvious opportunity is to use satellite images through GEE for anomaly detection. Nowadays, most farmers have a polygon layer containing the borders of all their fields. After uploading this layer to the GEE platform, the farmer can look for anomalies within each field along the growing season. The only limitation is a clear sky that enables acquiring a good image. An anomaly can be a result of many abiotic or biotic factors such as water leakage/blockage from/in an irrigation system, fertilization level, soil compaction, weed/disease/insect infestation/infection, fallow, salinity, etc. Anomaly identification does not point to the cause of the anomaly, but it can assist farmers to monitor the problematic points in the field.
We based the GEE algorithm on NDVI since it is sensitive to almost every biotic and abiotic (significant) anomaly [36] and is provided at the highest spatial resolution. Although Sentinel-2 red-edge VI (REIP) is provided at a lower resolution (20 m), it was included in the GEE anomaly app, as it is sensitive to wheat crop conditions even in 24 m wide reference strips ( Figure 10). These traits emphasize the possibility to use the GEE anomaly app in the "real world" immediately.
Reference fertilization strips are based on a lower or higher level of base/to-applied nitrogen. Where the lower N level has a significant effect, the N deficiency symptom should be shown by a decrease in biomass production and/or reduce N content in wheat dry matter. Conversely, significant superfluous N levels should be shown by an increase in biomass production and/or increase N content in wheat dry matter. In both cases, wheat above ground N uptake (biomass multiplied by N concentration) represent the real situation and can be best estimated by REIP ( Figure 5). Therefore, the REIP GEE algorithm is suggested for N reference strip estimations. The satellite-based correlations were in some cases lower than proximal sensing data. Nevertheless, Figure 5 shows that the correlations of the N uptake estimations are high enough to be used (with the S2 L1 REIP) for making decisions based on reference strips. This is in agreement with the use of red-edge-based VIs for reference strip or nitrogen uptake in wheat, rice, and corn [22,51,57,58].
Several algorithms have been developed for fertilization and yield estimation based on proximal sensing sensors such as RapidScan, GreenSeaker, Yara, and other similar sensors based on NDVI and/or NDRE. When proximal and remote sensing occur within 4 days difference, high correlations exist between RapidScan NDVI and NDRE, and L1 VENµS and Sentinel-2 VIs (NDVI 0.62 and 0.7, NDRE 0.83 and 0.77 for VE and S2, respectively). Hence, algorithms that have been developed for topdressing nitrogen application based on proximal sensing should work well also with remote sensing data. There are two main limitations: (1) image availability: the revisit time of Sentinel-2 is 5 days; hence, real applications could face limiting cloudy conditions. In some cases, it can continue for 2-3 or more image dates. Therefore, reducing the revisit time, such as the two days VENµS revisit resolution, is very important; (2) reference strip: Since the reference strip should be identified by the sensor, it should include pure pixels. Sentinel-2 NDVI and REIP/NDRE/ICCI spatial resolution is 10 and 20 m, requiring larger reference strips. VENµS spatial resolution is 5 m for both VIs, enabling the use of smaller reference strips. Hence, the high temporal and spatial resolution of VENµS enables monitoring agricultural fields during the growing season better than the lower resolution Sentinel-2.
It is important to stress that, in many cases, yield and or yield quality did not vary between the reference area and the whole field. In some cases, other biotic/abiotic stress canceled the effect of the reference nitrogen strip. For example, in field 101N (Figure 8), the reference strip "effect" was canceled by herbicide drift from a nearby field that affected grain yield and quality. In some cases, excessive N base application in the reference strip encouraged vegetative growth, which was shown by the Anomaly Explorer app but the crop utilized too much soil water during that growth stage, leaving insufficient water during the grain filling stage and even causing grain yield reduction. Hence, many fields exhibited lower correlations than field #109 (Figures 9-11). Therefore, the GEE anomaly app can be used for decision making as regards topdressing N application, but there are some major issues that must be addressed, and more research is required. the most important points among them are highlighted as follows: • Anomaly: Which calculation will best support decision making-deviation of variance/STD or deviation of percentage from the mean? • Threshold: What kind of threshold will best support a decision to apply topdressed N?Is it a fixed threshold or should it be affected by parameter(s), such as the VI used, phenology stage, cultivar, fallow, precipitation? • Reference strip: It is obvious that more pure pixels would enable better identification, but the minimum strip size should be determined, as well as its shape (square/rectangle), distance from the field edge, and orientation (facing north/ignore direction). • Cloud/Shade: Clear sky provides optimal conditions, while complete cloud cover prevents using the image. However, in many images, part cloudy condition prevails, causing complications such as areas covered by clouds that cannot be used, and areas shaded by clouds (or semitransparent clouds) that are questionable to be used for reference strip decisions. There are some confusing preliminary data showing that anomalies can or cannot be detected under these problematic conditions, which requires further studies to determine the limits.
To conclude, satellite images can be used in the GEE platform for anomaly detection, rendering results within a few seconds. The ability to use VENµS or Sentinel-2 data through GEE opens the opportunity to use these data for several applications by farmers and others.