Site-specific nitrogen management in winter wheat supported by low-altitude remote sensing and soil data

Site-specific nitrogen (N) management in precision agriculture is used to improve nitrogen use efficiency (NUE) at the field scale. The objective of this study has been (i) to better understand the relationship between data derived from an unmanned aerial vehicle (UAV) platform and the crop temporal and spatial variability in small fields of about 2 ha, and (ii) to increase knowledge on how such data can support variable application of N fertilizer in winter wheat (Triticum aestivum). Multi-spectral images acquired with a commercially available UAV platform and soil available mineral N content (Nmin) sampled in the field were used to evaluate the in-field variability of the N-status of the crop. A plot-based field experiment was designed to compare uniform standard rate (ST) to variable rate (VR) N application. Non-fertilized (NF) and N-rich (NR) plots were placed as positive and negative N-status references and were used to calculate various indicators related to NUE. The crop was monitored throughout the season to support three split fertilizations. The data of two growing seasons (2017/2018 and 2018/2019) were used to validate the sensitivity of spectral vegetation indices (SVI) suitable for the sensor used in relation to biomass and N-status traits. Grain yield was mostly in the expected range and inconsistently higher in VR compared to ST. In contrast, N fertilizer application was reduced in the VR treatments between 5 and 40% depending on the field heterogeneity. The study showed that the methods used provided a good base to implement variable rate fertilizer application in small to medium scale agricultural systems. In the majority of the case studies, NUE was improved around 10% by redistributing and reducing the amount of N fertilizer applied. However, the prediction of the N-mineralisation in the soil and related N-uptake by the plants remains to be better understood to further optimize in-season N-fertilization.


Introduction
In Swiss agriculture, farm size is considered medium to small scale averaging around 20 ha (Swiss Federal Statistical Office FSO 2019). The scale is comparable to other small scale farming systems in the EU-28 countries, in which 85% of the farms have a size between 0.1 and 20 ha (EUROSTAT 2016), and partly to village farm systems in the North China Plain (Chen et al. 2019;Zha et al. 2019) composed of family households managing 0.3-0.5 ha. At this scale, monitoring and management of in-field variability are confronted with additional obstacles compared to large-scale precision farming operations. Specifically, farmers have a lower level of specialization and smaller investment potential. In Switzerland, the fertilizer recommendation is usually based on the "Principles of fertilization of agricultural crops in Switzerland" (PRIF) (Sinaj and Richner 2017). The PRIF recommends a total N fertilizer amount given a certain yield expectation and other soil and agronomic factors and is based on a broad long-term dataset with replicates all over Switzerland. For winter wheat, the recommendation ideally includes an initial soil N min test to account for available soil N and it recommends up to three split applications at key stages of crop development. Splitting the application has the advantage of reducing the risk of N loss by emission (Hunt et al. 2019) and providing better availability of N with a view to increased crop growth and grain quality. Nonetheless, N fertilizers are widely used and their improper application represents both an environmental risk and a cost factor (Jan et al. 2017;Spiess 2011). Therefore, precise fertilization is a key challenge for producers who need to manage their crop production systems in order to minimize N losses to air or water, while achieving high crop yields of good quality (Zebarth et al. 2009). In fact, the inaccurate calculation of fertilizer requirements is often the main cause of environmental problems such as eutrophication or pollution of water bodies and emission of the potent greenhouse gas N 2 O and thus for additional societal costs (Lassaletta et al. 2014). It is, therefore, of practical and political concern to increase nitrogen use efficiency (NUE) and reduce high N input in the agricultural sector. In 2014, the efficiency of N used in Swiss agriculture was around 30% (FSO 2019), which is lower compared to the Danish agriculture, for instance, depicting a NUE of 41% in 2012 (Hansen et al. 2017).
Fertilizers are still typically distributed in a uniform way across fields without considering in-field variability in Switzerland. This is also true for the majority of wheat production, which is the most abundant crop in Switzerland. However, the spatial and temporal distribution of available N for agricultural crops in the field varies greatly according to a multitude of factors including soil properties, climate and soil management (Kindred et al. 2015), which influence N supply, mineralization processes, plant response to N and subsequently growth, yield and quality (Samborski et al. 2009). Minimized or site-specific application of fertilizers in precision agriculture systems has the potential to mitigate leaching problems as well as the emission of greenhouse gases . In fact, it was shown that site-specific fertilization can increase the NUE at field scale (Basso et al. 2016;Cohan et al. 2018;Raun et al. 2002) which was recently confirmed for small to medium scale agriculture systems (Van Loon et al. 2018;Wang et al. 2019) similar to Swiss agriculture. Therefore, variable rate N fertilization has the potential to improve NUE in small to large-scale agricultural cropping systems (Diacono et al. 2013;Ebertseder et al. 2003;Ravier et al. 2018). Yet, previous studies have mostly taken fields of large size into account; in practice, small-scale heterogeneities within fields smaller than one ha are typically neglected.
At present, remote, aerial and ground based sensing techniques emerge as a key element for observation of in-field variability. A combination of remote and proximal sensing techniques will enable a better understanding of cropping systems, from mineralization dynamics in soils to crop growth and nutrient uptake (Gabriel et al. 2017;Khan et al. 2018;Nawar et al. 2017). The use of spectral information to detect changes in canopy structure and growth is a well-established technology for a multitude of platforms (Matese et al. 2015;Muñoz-Huerta et al. 2013;Walter et al. 2018). For N-status detection of plants, different 1 3 paths have been explored in basic research e.g. hyper-and multi-spectral image spectroscopy, which delivers responses as spectral vegetation indices (SVI), and radiative transfer models . Recently, it has been suggested that the relation to proteins is more consistent than the previously assumed relation to chlorophyll . When applying this knowledge in practical field management, many studies have made use of SVI calculated from the near-infrared (NIR) region of the light spectrum (790-840 nm) and surrounding wavelengths, as for example the red-edge (RE) region located around 730 nm. This region appears to have a consistent direct relationship with the biomass and N status of the crop-a relation shown to exist for many crops including wheat, maize and sugar beet (Cammarano et al. 2011;Li et al. 2014;Liebisch et al. 2014;Prey and Schmidhalter 2019). However, it has been suggested that reflectance in the RE region, which is higher than that in the red region and which is relatively constant throughout the vegetative season, is not more sensitive to N per se (Bean et al. 2018). Instead, indices that make use of the RE band seem to enhance the importance of the NIR reflectance. Because such spectral signals can be assessed with commercially available sensors and because SVI derived from them appear to capture the crop in-field variability in a consistent way, such SVI are typically used to display in-field variability maps for crop N-status.
Unmanned aerial vehicles (UAV) are platforms suitable for monitoring fields of small to medium size. Main advantages of these systems are their flexibility of use and their capability to deliver high spatial and temporal resolution of observations simultaneously. Main disadvantages are the initial financial and knowledge related investments and the time needed to acquire and process the remote sensing data (Hunt and Daughtry 2018). On the market, a broad selection of sensors and platforms is available ) and the quality of the obtainable data has reached levels which support precision fertilization methodologies. Still, the quantification of crop N-status and subsequent fertilization support based on remote sensing imagery is not fully standardized. Calibration with ground data usually improves the reliability of derived fertilizer application maps. However, it is still difficult to generalize the relationship of the sensed values to the crop N demand for seasons (temporal), regions (spatially) and different crop species or genotypes of the same species.
The objectives of this study were (i) to better understand the relationship between data derived from a UAV platform representing temporal and spatial variability of crops in small fields of about 2 ha, and (ii) to increase knowledge on how such data can support variable application of N fertilizer in winter wheat (Triticum aestivum). Four fields in total were treated as four case studies to show how the variability between and within fields influences the outcome of site-specific fertilization. The main hypothesis was that the implementation of site-specific N fertilization using VRA techniques would reduce average N application compared to the standard fertilization strategy without affecting yield and thus increasing NUE.

Experimental fields and design
The experimental fields are located at the "Swiss Future Farm" in Tänikon, Switzerland (47.4790021° N,8.9059287° E). The farm is a concept to test and show innovative agricultural technologies, operated by the three partners AGCO (Duluth, USA), GVS Agrar 1 3 (Herblingen, Switzerland) and the cantonal education and extension service (Arenenberg, Switzerland). The climate in this region is characterized by 1170 mm annual rainfall and an average annual temperature of 8.6 °C  recorded by a local weather station from the federal office of meteorology (MeteoSwiss). In the growing season 2017/2018, the first experiment was carried out in field F1, while in the growing season 2018/2019 the fields F2, F3 and F4 were used (Fig. 1).
The soil types of the four investigated fields (F1-4) were classified during a soil survey in 1977, provided by the national soil-monitoring network (NABODAT 2019). According to the world reference base for soils (WRB, Food and Agriculture Organization FAO 2014), the soil of F1 is characterized as a Gleysol, with a more stagnating Gleysol zone richer in organic matter in the central part. F2 is characterized as a Cambisol, whereas F3 is mainly characterized as a Luvisol and in F4 two main soil types, a Gleysol (to the west) and an Alisol (to the east) are present. The heterogeneity of soil properties was high between and within the fields. Among fields, the clay content varied between 25 and 35% and the content of C org fraction calculated from organic matter content (OM) ranged between 1.5 and 3.0% (Table 1). Within fields, the standard deviation of clay was between 7 and 9% in all four fields while F1 and F4 showed a higher standard deviation (0.6-0.9%) compared to F2 and F3. The fields were sown with winter wheat (T. aestivum) of the same cultivar (Arnold, Saatzucht Donau, Austria) and covered an area of two ha on average (Table 1). Whereas wheat followed maize in the crop rotation in F1-3, sown on the 19th of October in 2017 and 9th and 12th of October 2018, respectively, F4 followed 2 years of temporal grassland and was sown 1 month later on 5th of November 2018. For F1 a randomized block design with three N-fertilization treatments: standard (ST), variable rate (VR) and non-fertilized (NF), replicated in six blocks each (n = 6) was established. The dimension of a single plot was 15 × 50 m, to match the operational range of the used pneumatic fertilizer spreader (Rauch Aero 2215, Sinsheim, Germany). The N-fertilization treatments are described in detail in the section below (Fertilization and Variable Rate method). For the second year, a non-randomized block design with two treatments ST and VR replicated in three fields F2, F3 and F4 (n = 7) plus two reference treatments NF and N-rich (NR) differing in number and distribution in the fields were established with regard to fit the different soil zones identified before the season. The dimension of a single plot was 15 × 90 m, which was chosen to match the operational range of the disc spreader (Sulky X40 + ECONOV, Sulky-Burel, Châteaubourg, France) that was mounted on the cultivation tractor (Massey Ferguson S 5713, AGCO, Duluth, USA). The tractor was equipped with a Valtra Smart Touch terminal (Valtra Inc., Suolahti, Finland), featured with "Vario Doc Pro" software (Application Status: 2.2.2, BSP: nt03_171212_134701, Kernel: 3.0.35, Based on Revision: 26565, Built on Machine: nt0x-vm, AGCO, Duluth, USA).
The fields were managed without the use of chemical plant protection and growth regulators. Mechanical weeding by means of a tined harrow (Treffler TS1520HM3) was carried out in early spring. Fertilizer was applied in the form of mineral ammonium nitrate (Agroline-Landor Fenaco, Muttenz, Switzerland, composition: 24% N, 5% Mg and 8.5% S). No additional P and K were fertilized. The fields were harvested using a combine harvester (Fendt 5275 C PLI, AGCO, Duluth, USA).

Remote sensing of the crop
The low-altitude remote sensing UAV platform was composed of a multispectral camera "Parrot Sequoia" (Parrot, Paris, France) mounted on a quadcopter "P4P" (DJI, Shenzhen, China) flown over the field on a weekly to bi-weekly basis. Spectral information of the crop was recorded in four bands of the light spectrum, namely green (G) centred at 550 ± 40 nm, red (R) at 660 ± 40 nm, red-edge (RE) at 735 ± 10 nm and near-infrared (NIR) at 790 ± 40 nm. The output of the camera consisted of one separate image for each band. The images were captured while the camera was flown over the field in a single grid automatic flight plan generated by the software Pix4D Capture (Pix4D, Lausanne, Switzerland) by Table 1 Soil and field properties for the four experimental fields F1 (2018) and F2-F4 (2019) For average values reported: n = 6 a CO 2 release reaction to hydrochloric acid (+) indicates the presence of carbonate b C org = OM/1.724 (Howard 1965) 3.2 ± 0.9 2.4 ± 0.1 1.5 ± 0.3 2.8 ± 0.6 P (mg kg −1 ) c 1.9 ± 0.6 1.3 ± 0.5 2.7 ± 1.1 1.3 ± 0.46 K (mg kg −1 ) c 38.5 ± 14 17.9 ± 5.8 32.6 ± 9.5 18.7 ± 6.4 Mg (mg kg −1 ) c 381. the drone at 5 m s − 1 speed and 1.9 s time interval between the single images at 50-80 m height. The optimal flight parameters were calculated with the help of the PhenoFly Planning Tool (Roth et al. 2018). The ground sampling distance (GSD) of the singleband images was between 4.5 and 8.4 cm pixel − 1 depending on flight altitude. The raw images were processed in the photogrammetry software Pix4D Mapper (Pix4D, Lausanne, Switzerland, version 4.5). First, a radiometric correction was applied by using an Airinov reflectance target (Parrot Airinov, Paris, France), then the images were transformed to an orthomosaic to create a reflectance map and finally, they were georeferenced using ground control points (GCP). GCPs were annually established within and around each field before generating final index maps of the field. The normalized difference red-edge index (NDRE = (NIR-RE)/(NIR + RE); Barnes et al. 2000) was selected as N-status indicator. This spectral vegetation index has been chosen as it showed consistent relationship with N uptake of the plants (Argento et al. 2019;Basso et al. 2016;Li et al. 2018). Additionally, the normalized difference vegetation index (NDVI = (NIR-R)/(NIR + R); Rouse et al. 1974) was calculated, because it correlates closely to canopy cover and biomass (Liebisch et al. 2017;Tucker 1979;Tucker et al. 1980) and is therefore often used to provide decision support for variable rate N fertilization Walsh et al. 2013).

Soil and plant analysis
In both seasons, soil samples were collected in each field in early spring from six sampling locations based on the soil type mapping (approximately three sampling locations/ha). For each sampling location, six to eight soil cores were taken and mixed separately for the three depths 0-30 cm, 30-60 cm and 60-90 cm, respectively. Out of these samples N min available in the soil (Table 2) was measured as NH 4 (ammonium test) and NO 3 (Aluminium-sulphate extraction combined with a nitrate electrode) and converted in kg N ha − 1 using the reference method by the federal agriculture research centre (Agroscope 1995). For the depth of 0-20 cm, phosphorus (P), potassium (K) and magnesium (Mg) content, organic carbon (C org ), texture, pH and carbonate content were measured additionally, following the guidelines for the mandatory soil samples that farmers perform to be eligible for direct payments (ÖLN, 2018).
Plant biomass samples were collected at various growth stages during the season and the respective growth stages were recorded using the BBCH decimal scale (Meier et al. 2009). In the first year, the samples were collected two times at BBCH 32 (second knot) and at BBCH 84 shortly before harvest in two subsamples per plot each covering an area of 50 × 60 cm (four rows at 15 cm spacing) for a total of 36 samples. The BBCH 32 samples were dried at 60 °C for 48 h to quantify dry matter (DM) and afterwards N concentration (N c ) was measured with a Flash EA 1112 Series elemental analyser (Thermo Italy, Rodano, Italy) coupled to a Finnigan MAT Delta plus XP isotope ratio mass spectrometer (Finnigan MAT, Bremen, Germany). The samples at harvest (BBCH 84) were collected in two subsamples per plot as described for BBCH 32 above. Plants were cut at 10 cm over the soil and air-dried for 2 weeks. The grains were mechanically separated from the straw to differentiate between weights of grains and straw that together make up the plant dry weight. N concentration was measured in both the milled grains and straw at 7% water content. In the second year, biomass samples were collected four times at BBCH 25 (end of tillering), BBCH 31 (beginning of stem elongation), BBCH 39 (flag leaf) and BBCH 84 (pre-harvest, 15th July 2019). From each plot, three subsamples were collected, each covering an area of 50 × 60 cm (four rows at 15 cm spacing), for a total of 72 sampling locations across three fields. The samples were treated and analysed as described above for the previous season. From a subsample of the grains, one thousand grains were counted with a seed counter (Contador, Pfeuffer GmbH, Kitzingen, Germany) and weighed to determine the thousand-grain weight (TGW, g). Another subsample of the grains was used to determine the protein content (GVS Agrar AG, Herblingen, Switzerland). To investigate wheat N status, three plant traits i.e. N concentration (N c ), total N in the aboveground biomass (N up , Eq. 1) and the nitrogen nutrition index (NNI, Eq. 2) were chosen. These indicators are often used to assess or calibrate crop N fertilizer demand. The N uptake (N up , kg N ha − 1 ) was calculated by multiplying the dry matter biomass (DM, kg ha − 1 ) with the corresponding N c (%) of the plant sample: The NNI is an index based on the principle that N c decreases with growing biomass. The experimental dilution curve reflecting the relationship between DM and N c values can be used to relate the measured N c to a critical N crit i.e. the lowest N c necessary to obtain maximum biomass at a given growth stage (Justes et al. 1994;Prost and Jeuffroy 2007). NNI values < 1 indicate an N deficiency while NNI values > 1 indicate an N surplus. (1)

Fertilization and variable rate method
The selected fertilization strategy is based on the N min method suggested in the PRIF (Sinaj and Richner 2017) as well as the standard application of the farm manager at the experimental site. Based on the recommendations, N-fertilization was divided in three split applications targeting the growth stages: end of tillering, beginning of stem elongation and flag leaf, respectively. On the field F1, the mean 0-90 cm N min was used as a reference for the calculation of the first split application to be applied as shown in Table 3. For the ST treatment, an average field N min value was used as farmers usually do. For the VR, the specific N min of each plot was subtracted from the total. According to this procedure, e.g. plot 2 with N min = 30 kg N ha − 1 received 80 − 30 = 50 kg N ha − 1 fertilizer with the first split. For the second and third split, the values were qualitatively adjusted based on the NDRE index map produced few days before the fertilization. That means e.g. plot 10 with NDRE mean value of 0.25 (20% higher than the field average) received 20% less fertilizer i.e. 40 kg N ha − 1 in the second split.
In the second year, a different approach was selected. The mean 0-30 cm N min , considered more relevant for the initial phase of plant growth, was used as a reference for the calculation of the first split application for VR only, as shown in Table 3. The ST in this season was based on the farm managers' standard application (hence without applying N min correction). For the VR treatment, the N min values from the six locations were interpolated to create an N min map of the field. The fields were then divided in N min management zones (Mattei et al. 2020) and for each VR plot, the corresponding N min value was used to adjust the first split. For example plot 2 of field F2 was in between two zones with 18 and 24 kg N ha − 1 , respectively, therefore it received two variable applications of 80-18 = 62 kg N ha − 1 and 80-24 = 56 kg N ha − 1 in the first split. The second and third split applications were adjusted based on an application map produced by applying Eq. 3, experimentally derived from the first season, to an NDRE index map produced few days before fertilization. For the optimization of the second and third fertilization split in the VR treatment, the NDRE values were used to adjust the standard amount of fertilization (N ST , kg N ha − 1 ). Assuming that higher reflectance intensity corresponds to a better plant N status, N fertilizer (N fert,i , kg N ha − 1 ) values for VR were calculated by applying Eq. 3 to each pixel.
(3) N fert,i = N ST − N corr,i 1 3 The correction factor at pixel i in the field (N corr,i , kg N ha − 1 ) was calculated as in Eq. 4 by adjusting the N ST according to the relative difference of reflectance intensity at a point i (x i ) and the mean of the reflectance intensity over the whole field (x m ).
The fertilizer was spread in variable amounts over the field, according to the final prescription map, which was created using the software (NEXT Farming AG Office, version 1.8.1.14, Pfarrkirchen, Germany) and uploaded via the Vario Doc Server onto the tractor terminal.

Evaluation of N use efficiency
The NUE was evaluated by means of three different indicators: the apparent fertilizer recovery (AFR) (Kindred et al. 2015), also known as recovery efficiency, was used to quantify the N fertilizer recovered by the crop (Eq. 5). Therefore, N up from the NF plot (as a measure of the soil N supply) was subtracted from the N uptake of a treatment plot (VR, ST, or NR) and divided for the N fertilizer applied. AFR was calculated for both the total N up (AFR tot) and for the N up in the grains (AFR grain).
The partial factor productivity (PFP) was calculated according to Wang et al. (2019) as a measure of the relationship between the grain yield in a fertilized plot (Y fert , kg ha − 1 ) and the N applied in that plot (kg N ha − 1 ) in Eq. 6.
A simplified measure for the economic income, defined as the marginal return of N fertilization (E, CHF ha − 1 ) was calculated in Eq. 7 according to Wang et al. (2019).
where Y is the grain yield (kg ha − 1 ), P Y is the grain price (CHF kg − 1 ), N is the N fertilizer applied (kg N ha − 1 ), P N is the N fertilizer price (CHF kg − 1 ). The prices of wheat grain and N fertilizer were 0.52 and 0.45 CHF kg − 1 , respectively, in Switzerland in 2019. Currency exchange values are set at 1 CHF = 1.05 USD = 0.93 EUR (UBS 2020).

Statistics
The spectral data were extracted from the orthomosaics and analysed in a pipeline, which combined the GIS software QGIS (QGIS Development Team 2019, version 2.18) and the software RStudio Desktop (RStudio Team 2016, version 1.0.143) with R version 3.4.1 "Single Candle" (R Core Team 2017). The packages "sp", "raster", "gstat", "rgdal" and "dplyr" were used in the pipeline. The statistical analysis including the one-way ANOVA and Fisher-LSD test (package "agricolae") for the grain yield and efficiency parameters were also performed in Rstudio. In order to select the most sensitive vegetation index, twenty-three vegetation indices from the literature were correlated to selected plant traits namely DM, N c , N up and NNI using quadratic linear regression (y = x 2 ). R 2 and RMSE values were calculated and evaluated. The data preparation and quadratic linear regression for the sensitivity analysis between vegetation indices and selected plant traits were performed in Rstudio. The data preparation and graphics for yield and spectral data were done with the packages "dplyr", "reshape2" and "ggplot2".

Grain yield
The grain yield in all fields showed no significant difference between the VR treatment and the ST treatment (Fig. 2). The fertilized treatments in fields F1, F2 and F3 were in the same range with no significant differences. The NF plots had significantly lower yield than the other treatments, except in field F4, where no significant differences were observed and yield was generally lower. Additionally, the NF treatment was not significantly different from the ST treatment in F2, where ST also tended to be lower than VR and NR. Average yields of fertilized plots were 6.9 t ha − 1 for both F2 and F3 and 3.9 t ha − 1 for F4. The nonfertilized plots reached 5.7, 5.3 and 4.3 t ha − 1 for F2, F3 and F4, respectively.
The thousand-grain weight (TGW) after harvest (Table 4) showed no significant differences between ST and VR (43.74 and 43.73 g) in F1 while both were above the NF weight (41.9 g). In 2019, the pattern was inverted; in fact, the grains in the NF plots were in all three fields about 10% heavier than the average of the fertilized plots. Protein content measured in 2019 was on average 15% for the fertilized plots with no significant differences between VR, ST and NR treatments while in the NF plots was on average 13.7%. The protein content in F4, which yielded less grain in total, was higher compared to F2 and F3. Error bars indicate the standard error. Letters denote significant differences between treatments according to Fisher-LSD test at p < .05. The values inside the columns represent the average N rate (kg N ha − 1 ) applied in that treatment and field (Color figure online)

Efficiency assessment
The efficiency analysis showed a consistent trend whereby VR treatments performed better than ST treatments (Table 4). There were no differences in total N up between ST and VR in grains, straw and total crop biomass; however, both were significantly higher than the NF in three of the four cases (F1-3). The soil N supply estimated from the NF plots in 2018 was 84 kg N ha − 1 , whereas it increased in the second year ranging from 127 to 143 kg N ha − 1 . The N c measured in the grain after harvest ranged from 1.9 to 3.1% of DM and from 0.3 to 1.0% of DM in the straw. Generally, a decrease of N c from the beginning of the season (mean 4%) until harvest (mean 1.5%) could be observed over all four fields. The NNI values in the grain at harvest ranged from 0.6 (NF plots) to 1.2 (NR plots) with values of VR and ST plots around 0.9 on average in the four fields.
Compared to ST, the reduction of average applied N in the VR treatments resulted in a trend of higher AFR (Table 5). However, these differences were not significant neither for the total crop N nor for the grain N uptake (N up ) only. The efficiency of grain production in relation to the N applied (PFP) showed a better performance of VR compared to ST in all four fields over 2 years. The marginal returns of VR also showed a trend of improved financial gain, when compared to the ST treatment. The differences ranged between 31 and 335 CHF ha − 1 due to the reduction in applied fertilizer and, in two cases, the increase in grain yield. The NR had higher marginal returns in two out of three cases. However, these differences were statistically not significant.

Sensitivity analysis
The sensitivity analysis between selected vegetation indices and plant traits (Table 6) performed on the joint dataset over the 2 years (n = 254) showed that by using a quadratic regression, the SVI with the highest sensitivity to both dry matter yield and N up are those making use of the red-edge channel. The NDRE shows the best correlation and lowest error with DM (R 2 = 0.72), N up (R 2 = 0.80) and NNI (R 2 = 0.75). Very similar results are achieved with the simple ratio NIR/RE with DM (R 2 = 0.70), N up (R 2 = 0.80) and NNI (R 2 = 0.75). The N c , however, showed generally low correlations and the best (0.40) was achieved with the combined index MCARI/MTVI2. As mentioned above, NDRE was the best SVI, showing a quadratic relationship to dry matter (DM), nitrogen uptake (N up ), nitrogen nutrition index (NNI) (Fig. 3). The N c curve instead, showed a quadratic decrescent curve with higher index values corresponding to lower N c values. The best index representing this relationship was the combined index MCARI/MTVI2, which, however, results in a very low correlation. Data scattering is particularly enhanced in the second biomass sampling of 2019.

Seasonal dynamic of NDRE compared to NDVI
The seasonal development of NDRE and NDVI of the studied fields showed different characteristics over time. The NDRE data extracted from each plot and averaged per treatment from field F1 to F2 showed a steady increase from the time of the first split fertilization (beginning of stem elongation), to its peak at the beginning of the spike emergence (Fig. 4, left). After flowering and during senescence the intensity of the signal decreased. The N input treatments were clearly separable during the assumed fertilization period from stem elongation to spike appearance. In F1, the fertilized treatments VR and ST were significantly different from the NF treatment starting from DAS 200. In F2, NR was significantly higher than NF from DAS 195 on. AT DAS 203, there was a significant difference between treatments: NR > VR, ST > NF. The fertilized treatments were then in the same range and separable from NF until late senescence. In contrast to NDRE, NDVI (Fig. 4, right) saturated at the beginning of the season shortly after stem elongation before decreasing rapidly during the senescence phase. The fertilizer treatments generally did not show distinct differences between each other. In F1, the differences between NF treatment and the fertilized treatments VR and ST were significant later in the season from DAS 214 on. In F2, there were no significant differences between treatments. In comparison, with NDRE the relative difference between fertilized and non-fertilized plots was clearly visible. The fields F3 and F4 showed a similar behaviour (Supplementary Fig. 1).

Discussion
Yield analysis confirmed the hypothesis that site-specific fertilization is able to reduce total N application without a loss of yield in the small-to-medium-sized agricultural system investigated in Switzerland being in line with the findings of Stamatiadis et al. (2018) and Wang et al. (2019). In three out of four cases (F1, F2 and F3), grain yield (6.5-7.5 t ha − 1 ) and protein content (14-15%) were in the range of the Swiss average for top quality wheat varieties (Levy and Brabant 2016;Sinaj and Richner 2017;Hofer 2019). In F4, the lower yield and lack of differences between all treatments could be explained by the combination of the late sowing date (1 month later than F2 and F3) and excessive damage provoked by the mechanical weeding with the tined harrow combined with the insufficient establishment of the wheat plants in early spring. Consequently, the field showed very heterogeneous patches of crop development (data not shown).

Remote sensing as a base for variable rate application
A literature survey reflected that a broad range of spectral indices is used for the prediction of in-field N status of crops (Heege et al. 2008;Heege 2013;Li et al. 2014;Prey and Schmidhalter 2019;Zhao et al. 2018). To confirm the applicability of the selected SVI for the given camera setting, a sensitivity analysis was conducted investigating twenty-three SVI from the literature survey applicable to the camera band combination. Most of them were shown to be sensitive to either canopy structure, biomass or N status (Basso et al. 2016;Bendig et al. 2015;Chen 2015;Cilia et al. 2014;Schmidhalter et al. 2003). The study clearly showed that indices combining RE and NIR like NDRE and NIR/RE ratio indeed had the best correlation with the N-status traits N up and NNI in this study. The correlation value (R 2 = 0.80 and R 2 = 0.75 for the N up and NNI, respectively) confirmed the assumption that NDRE can be used to monitor the N status of the wheat field. Bean et al. (2018) reported similar observations for corn N side dress also pointing to a superiority of RE-NIR based spectral indices over red-NIR based ones like NDVI. The seasonal monitoring patterns of NDVI and NDRE reflected climate and field management affecting the plants' growth. The NDVI curves saturate at the time of the second split fertilizer application, earlier than the NDRE. In Fig. 4, NDVI curves for F1 and F2 showed a high saturation already at 190 days after sowing corresponding to BBCH 29-30, between canopy closure and beginning of stem elongation (Baret and Guyot 1991). At this stage, NDRE is still increasing with values around 0.3-0.4. Because of the early saturation, NDVI is less useful for the derivation of N fertilizer prescription maps for the second and third fertilizer application. However, NDVI seems to be a viable N status indicator for a first N application when the canopy is not yet closed. The NDRE development was linear until the stage of spike emergence, which takes place after the third and last fertilizer application in winter wheat and therefore NDRE and NIR/RE are a better base for the creation of fertilizer prescription maps than NDVI or the other investigated SVIs.
Although NDRE allows quantifying N up , the relatively large variability suggests a qualitative assessment, linking NDRE to the general fertilization strategy. In this study, the N fertilizer prescription for VR was based on the NDRE index map combined with the Swiss standard fertilization recommendation indicating the average N demand subsequently adjusted for higher or lower NDRE values. The choice to use the average value as reference was supported by the presence of non-fertilized and surplus-fertilized plots in the trial acting as minimum and maximum reference values, respectively. However, other approaches to obtain fertilizer prescription maps were suggested. Holland and Schepers (2013) and Stamatiadis et al. (2018) used the 95th percentile of the histogram produced with spectral index values of the whole field. In other studies, the yield response to N-application in combination with spectral indices were used to predict the N-application rate (Bean et al. 2018;Franzen et al. 2016;Holland and Schepers 2010;Raun et al. 2005). This study indicates that the use of NDVI is less sensitive for in-season fertilization support, despite being a commonly used solution in other studies Walsh et al. 2013;Wood et al. 2003). Nevertheless, for early season fertilizer application NDVI seems to provide a reliable approximation of biomass. In-depth studies are needed to clarify the sensitivity change of spectral indices in relation to crop development features like canopy coverage to suggest the optimal sensing strategy to support in-season fertilizer application.
For the determination of the field variability by imaging, the use of a UAV system was superior to satellite data in this study (data not shown), because the plot size of the experiments was relatively small and delineation of the plots with satellite pixels such as Sentinel-2 was not robust, creating interference in the plot measurement precision. However, if considering the whole field level, the identification of N-status and extraction of fertilization prescription maps would be feasible from such satellite data. In fact, the information contained in the UAV images needed to be down-sampled to match the resolution of the application map for the tractor terminal of a minimum 7 × 7 m (derived by testing). In general, satellite images with a resolution of 10 × 10 m such as for Sentinel-2 are therefore a viable basis for creating prescription maps. Nonetheless, the drone-based approach allowed to study the vegetation development in more detail and to obtain information even on cloudy days, which can be an advantage in climates and regions in which the first and the second N fertilization split is usually performed during cloudy spring periods like in Switzerland. Therefore, the outlined approaches might lead to a more favourable use of small-scale equipment for the management of small-to-medium-sized agricultural systems in future. Based on such approaches, it might become possible to return to smaller sized tractors or even to use smaller machinery such as robots or small autonomous tractors equipped with pneumatic spreader systems, which can distribute fertilizer at a finer scale and thus may use prescription maps with higher resolution.
The characterization of the soil variability and crop growth variation is an important step to understand, whether variable rate application methodology is viable for a certain field or not (Griepentrog et al. 2007;Heege 2013). The four field case studies presented in this study showed various degrees of variability mainly due to underlying soil properties but also due to weather differences. Field F3 was relatively homogeneous and therefore the average saving of N fertilizer was only about 5%. Nonetheless, the redistribution of fertilizer had a positive effect on yield and NUE. The influence of climate was crucial for F1, as 2018 was a very dry season (475 mm cumulative rainfall and mean T 10.5 °C in January-July 2018). However, for field F1 the combination of the dry season and a commonly water-stagnating soil was beneficial. In the central zone of the field, which had a higher level of organic matter and water holding capacity, likely supporting a higher N mineralisation, it was possible to reduce the fertilizer amount down to 30-35% in two VR plots, without loss in yield. The season 2019 instead was a wet season with generally beneficial weather (657 mm and T 9.5 °C January-July 2019) which was reflected in the higher soil N supply in the NF compared to the previous year. The favourable weather conditions combined with the high N-mineralisation rate of the soil during the vegetation period likely led to a good grain filling and higher TGW in the NF-plots. Ultimately, a deeper understanding of N mineralization potential, related to spatially and temporally varying environmental conditions, is necessary to define the most efficient fertilization strategy. The use of models or sensor-based information could help to better understand the dynamics of N throughout the season (Yin et al. 2020).

Environmental and economic efficiency
A general trend towards improved efficiency was observed in three of the four cases studies. The total wheat crop N up was similar between fertilized treatments, indicating that the impact of the reduction in applied N on N up from the VR treatments was compensated by the improved spatial distribution. The reduction of applied N over all plots and the 4 years ranged from 5 to 40% of the standard N fertilization. PFP was the only efficiency indicator that resulted in significant differences between treatments. In this sense, the yield production in relation to the N applied was improved with VR application compared to ST. Whereas in 2018 a very high AFR of around 90% was obtained, in 2019 lower values around 24-69% were observed. The reason was that the soil in 2018 only supplied 84 kg N ha − 1 (NF plots) due to the very dry climate. In 2019, the N replenishment from the soil was much higher (127-143 kg N ha − 1 ). This might mean that the subsequent soil supply largely covered the N requirement of the plants and thus led to a very low N fertilizer utilization. Although it is not possible to know which part of the applied N or soil N was taken up by the plants, this is a clear indication of lower fertilizer NUE and of an increased risk of N loss. In 2019, F2 had higher N min (Table 2) and yield, therefore it is conceivable that the lower first application resulted in higher efficiency compared to F3 and F4. The use of image or sensor data is not sufficient to identify seasonal total N needs of the crop with high precision but arguably, the set of calculation applied to the data might have been a limitation.
In this study, it was not possible to directly estimate the emission of N via leaching and denitrification. However, the active N pool in the system can be estimated by summing the N up from the NF treatments, as measure for soil N supply during the winter wheat season, to the N applied in the fertilized treatments. Subtracting N up in the fertilized treatments from this active N pool reflects the approximate quantity of N left in the system after harvest prone to be lost or immobilized. In both years, the VR treatments showed a lower N loss risk. Whereas in 2018 the values were generally low and risk for N loss was 50% lower for VR than for ST (10 kg N ha − 1 for VR and 21 kg N ha − 1 for ST in F1), in 2019 the values were higher and VR was on average 30% lower than ST and NR (63 kg N ha − 1 for VR and 90 kg N ha − 1 for ST and NR in F2-F3) (Supplementary Table 2). Both VR application of fertilizers and considering soil N analysis for optimal N supply are viable methods to reduce the risk for N losses. The combination of both very likely creates a synergy in the reduction for field and farm management of N.
Marginal returns offer a simplified economic balance between the cost of fertilizer and the gain from the sale of the grain to the mill. The improved gain of VR when compared to the ST ranged from 1.6% corresponding to 31 CHF ha − 1 to 9.3% corresponding to 335 CHF ha − 1 . This evaluation does not take parameters such as the investment costs to obtain prescription maps or cost for the technology and machinery into account. Furthermore, in small-scale farming, a higher marginal return per ha might not be fully sufficient to sustain the required investments. Considering a scenario in which sufficient knowledge is available to the point that higher returns for a certain threshold of variability in the field can be guaranteed, cantonal or governmental authorities could consider measures to support the transition to VR technology for farmers. Another strategy could be to encourage farm contractors providing VR fertilization services.

Conclusions
The methods applied were suitable to characterize in-field variability and were able to increase the nitrogen use efficiency. Also under Swiss conditions, NDRE showed a better differentiation of the N-status of the plants, than NDVI. The potential for further improvements lies primarily in the extent of the variability of N availability in the soil and resulting crop growth within the field, but also in a better understanding of mineralization processes. Easily collectable visual data that cover fields with a high spatial resolution can help to close this gap. Ultimately, better quantification of the variability in terms of N status and plant growth is necessary to set a threshold for deciding if variable rate application is worth being implemented in a designated field. However, VRA of N has the potential to reduce fertilizer inputs while keeping yields at current levels, in particular when combined with soil N min information. Thus, it represents a viable tool among others to improve NUE in cropping systems and reduce N losses to the environment. The methodology may support increasing sustainability of small to medium scale agriculture by increasing financial return and decreasing the environmental footprint of arable cropping systems.