The Value of Sentinel-2 Spectral Bands for the Assessment of Winter Wheat Growth and Development

: Leaf Area Index (LAI) and chlorophyll content are strongly related to plant development and productivity. Spatial and temporal estimates of these variables are essential for e ﬃ cient and precise crop management. The availability of open-access data from the European Space Agency’s (ESA) Sentinel-2 satellite—delivering global coverage with an average 5-day revisit frequency at a spatial resolution of up to 10 metres—could provide estimates of these variables at unprecedented (i.e., sub-ﬁeld) resolution. Using synthetic data, past research has demonstrated the potential of Sentinel-2 for estimating crop variables. Nonetheless, research involving a robust analysis of the Sentinel-2 bands for supporting agricultural applications is limited. We evaluated the potential of Sentinel-2 data for retrieving winter wheat LAI, leaf chlorophyll content (LCC) and canopy chlorophyll content (CCC). In coordination with destructive and non-destructive ground measurements, we acquired multispectral data from an Unmanned Aerial Vehicle (UAV)-mounted sensor measuring key Sentinel-2 spectral bands (443 to 865 nm). We applied Gaussian processes regression (GPR) machine learning to determine the most informative Sentinel-2 bands for retrieving each of the variables. We further evaluated the GPR model performance when propagating observation uncertainty. When applying the best-performing GPR models without propagating uncertainty, the retrievals had a high agreement with ground measurements—the mean R 2 and normalised root-mean-square error (NRMSE) were 0.89 and 8.8%, respectively. When propagating uncertainty, the mean R 2 and NRMSE were 0.82 and 11.9%, respectively. When accounting for measurement uncertainty in the estimation of LAI and CCC, the number of most informative Sentinel-2 bands was reduced from four to only two—the red-edge (705 nm) and near-infrared (865 nm) bands. This research demonstrates the value of the Sentinel-2 spectral characteristics for retrieving critical variables that can support more sustainable crop management practices. for a linear ﬁt between the SunScan and destructive sample LAI measurements. We further corrected the bias in the SunScan measurements based on the destructive data, which reduced the NRMSE from 50% to 17%. A quadratic ﬁt demonstrated a high correlation (R 2 = 0.75) between the LCC measurements and leaf N content derived from the destructive sample analysis. The variance of the measured LCC was, however, relatively high when compared to that of the LAI, with the standard deviation of SPAD measurements ranging from +/ - 3 to +/ - 16.


Introduction
Leaf Area Index (LAI) and chlorophyll content are essential indicators of crop phenological status and condition, which can be used to support a range of precision agricultural technologies. For instance, LAI is a key biophysical parameter that quantifies plant canopy structure and function. LAI is, therefore, related to canopy-scale processes, including evapotranspiration, photosynthesis, respiration and the interception of precipitation and solar radiation [1,2]. Consequently, past research has demonstrated with the exception of Clevers et al. [24] where potato leaf chlorophyll content (LCC) and CCC were estimated from a VI derived from real Sentinel-2 observations, studies involving the use of measurements that matched the spectral characteristics of Sentinel-2 are sparse.
Whilst the potential of Sentinel-2 data for supporting agricultural applications has been investigated, a thorough sensitivity analysis of each Sentinel-2 band for deriving crop variables across important crop growth stages is not well documented. This research addresses this knowledge gap by evaluating the characteristics of Sentinel-2 spectral data for retrieving key winter wheat variables-LAI, LCC and CCC. We used multi-temporal data acquired from field campaigns, which included non-destructive direct measurements of LAI and LCC at experimental winter wheat N trial plots. These measurements were also compared to data derived from analysing destructive samples. In conjunction with these ground measurements, we acquired data from an Unmanned Aerial Vehicle (UAV)-mounted multispectral camera comprising nine sensors measuring the same key central wavelengths, ranging from 443 to 865 nm, and with the same bandwidths (full width at half maximum response (FWHM)) as that of the Sentinel-2 multi-spectral instrument. Where cloud cover often limits the availability of optical Earth observation data, our use of a UAV platform ensures observations and allows us to thoroughly explore the sensitivity of the Sentinel-2 bands on the ground data. Our research objectives were, firstly, to characterise the uncertainty and correct biases in the non-destructive ground measurements based on the destructive sample data. Second, to investigate the impact of propagating the quantified uncertainty when training the Gaussian processes regression (GPR) machine learning algorithm, which was used to determine the most informative Sentinel-2 bands for retrieving LAI, LCC and CCC. We further compared the performance of the GPR model to a simpler multivariate linear fit derived from the most informative Sentinel-2 bands for estimating each crop variable.

Field Site and In Situ Crop Measurements
Multi-temporal field campaigns, carried out during a 2017/2018 winter wheat growing season, involved acquiring sets of ground measurements in coordination with UAV-mounted multispectral camera observations over an experimental field trials site.

Experimental Trial Plot Description
The field experiment included a total of 50 winter wheat (Triticum aestivum L.) trial plots located approximately 3.8 km south of the village of East Saltoun, East Lothian, Scotland (55 • 52 50.5" N, 2 • 50 12.2" W; 170 m above sea level). The trial plots, with dimensions of 2 × 10 m, had a Latin square experimental design and, in order to induce variation in LAI, LCC and CCC, comprised five different levels of N application-0, 50, 100, 150 and 200 kg N ha −1 (Figure 1) -each with five replicates for two common soft group 4 winter wheat varieties. These two wheat varieties-Revelation and Leeds-were included on the 2017/2018 recommended lists published by the U.K.'s Agriculture and Horticultural Development Board (AHDB) for cereals crops [34].
Each of the wheat plots were sown on 30 September 2017, in a roughly southwest to northeast direction with a seed rate of 340 seeds/m 2 and harvested on 25 August 2018. The soil was of the Humbie soil series with a loam texture. The N fertilisation of plots was carried out as a split application-50% of the total N was applied on 22 March 2018 and the remainder was added on 26 April 2018, which corresponded to growth stages (GS) 24 and 31, respectively. A herbicide based on the active ingredients picolinafen and pendimethalin was applied at GS11 on 27 October. A robust fungicide programme based on the active ingredients triazole, chlorothalonil, cyflufenamid, proquinazid, succinate dehydrogenase inhibitor and azoxystrobin, was also applied to all plots at four growth stages (GS30 on 16 April, GS32 on 9 May, GS39 on 30 May and GS65 on 22 June) to keep all diseases to a minimum level throughout the growing season. Each of the wheat plots were sown on 30 September 2017, in a roughly southwest to northeast direction with a seed rate of 340 seeds/m 2 and harvested on 25 August 2018. The soil was of the Humbie soil series with a loam texture. The N fertilisation of plots was carried out as a split application-50% of the total N was applied on 22 March 2018 and the remainder was added on 26 April 2018, which corresponded to growth stages (GS) 24 and 31, respectively. A herbicide based on the active ingredients picolinafen and pendimethalin was applied at GS11 on 27 October. A robust fungicide programme based on the active ingredients triazole, chlorothalonil, cyflufenamid, proquinazid, succinate dehydrogenase inhibitor and azoxystrobin, was also applied to all plots at four growth stages (GS30 on 16 April, GS32 on 9 May, GS39 on 30 May and GS65 on 22 June) to keep all diseases to a minimum level throughout the growing season.

In Situ Measurements
For five different dates within the growing season (Table 1), non-destructive measurements carried out in each experimental trial plot included LAI, LCC and growth stage observations in accordance with the Zadoks decimal code [35]. Five technical replicates of LAI were taken per plot, with ten replicates of chlorophyll content on a regular grid within each plot. The replicates were then combined to give a plot average and standard deviation. The growth stage was assumed to be reached when it was observed in at least 50% of the plots. LAI was measured using a SunScan device (Delta-T Devices, Cambridge), and LCC measurements were inferred from a portable soil plant analyses development (SPAD) meter device (Konica Minolta, Japan). The CCC, expressed per unit of leaf area, was calculated as the product of the LAI and LCC [11,24,27]. Across the five observation dates and 50 trial plots, a total of 250 sets of LAI, LCC and CCC were derived from the ground measurements and were used in this study.

In Situ Measurements
For five different dates within the growing season (Table 1), non-destructive measurements carried out in each experimental trial plot included LAI, LCC and growth stage observations in accordance with the Zadoks decimal code [35]. Five technical replicates of LAI were taken per plot, with ten replicates of chlorophyll content on a regular grid within each plot. The replicates were then combined to give a plot average and standard deviation. The growth stage was assumed to be reached when it was observed in at least 50% of the plots. LAI was measured using a SunScan device (Delta-T Devices, Cambridge), and LCC measurements were inferred from a portable soil plant analyses development (SPAD) meter device (Konica Minolta, Japan). The CCC, expressed per unit of leaf area, was calculated as the product of the LAI and LCC [11,24,27]. Across the five observation dates and 50 trial plots, a total of 250 sets of LAI, LCC and CCC were derived from the ground measurements and were used in this study. Destructive analyses were carried out to sample the winter wheat vegetation at five of the revelation trial plots, covering the range of N application from 0 to 200 kg ( Figure 1). These destructive measurements, conducted on three dates (25 May,13 June and 4 July 2018), were used to correct biases and quantify the uncertainty of LAI and LCC non-destructive measurements, which were also acquired on the same day as the destructive sampling. The measurements entailed randomly placing a 0.25 m 2 quadrat within the destructively-sampled plots and removing all above-ground vegetation. The leaves were separated from the remainder of the vegetation and the LAI was then estimated by passing the collected leaves through a Li-3100C leaf area meter (Li-Cor, Nebraska, USA). The non-destructive SunScan LAI estimates could then be directly compared to the destructive measurements to firstly correct for biases. This bias correction was performed through reduced major axis linear regression, which accounts for the variance in both the destructive and non-destructive measurements [36]. Specifically, the resultant linear fit equation was used to correct the bias in the non-destructive measurements. The uncertainty of the bias corrected non-destructive measurements was then calculated as the normalised root-mean-square-error (NRMSE): where E i and O i represent the bias corrected non-destructive and destructive values, respectively, and n is the number of non-destructive and destructive comparisons. Past research has demonstrated a significant relationship between LCC and leaf N content [14]. We, therefore, quantified the uncertainty of the in situ LCC measurements by estimating the leaf N content of the destructive samples. This analysis involved drying and milling the samples for each vegetation material type (i.e. stem, leaves and ears). The samples were then weighed and analysed for percentage C and N content using a Flash 2000 elemental analyser. The SPAD LCC estimates were then compared to the average leaf N percentage for each of the five destructive sample plots.

UAV Platform and Multispectral Instrument
The UAV flights took place over the wheat trial plots on the corresponding ground measurement dates (Table 1) using a hexa-copter platform (DJI Matrice 600) at a height of 100 m above ground level and a constant speed of 3.1 m/s. Each of these flights were carried out close to solar noon (between 11:00 and 14:00 GMT+1) to avoid errors due to a low solar elevation angle. In order to ensure consistent and comparable multi-date coverage, the platform was equipped with a real-time kinematic Global Navigation Satellite System (GNSS) and was flown autonomously using the same flight mission, which was pre-programmed using the DJI Ground Station Pro software, for each of the measurement dates.
Multispectral imagery was acquired throughout the UAV flights using the MAIA camera system (SAL Engineering/EOPTIS), which is composed of an array of nine monochromatic sensors, each having a 1.2 Mpixel resolution. Furthermore, the nine sensors of the MAIA camera (referred to hereafter as the MAIA/Sentinel-2) have band-pass filters that have the same central wavelength and width as that of the first nine bands (i.e., bands 1 to 8A, Table 2) of the ESA Sentinel-2 multispectral instrument [37]. During the UAV flights, these sensors imaged the trial plots from a fixed nadir position with the aid of a three-axis stabilisation gimbal (DJI Ronin-MX). The pixels in the nine MAIA/Sentinel-2 sensors collected data simultaneously via global shutters, thus allowing all the nine band images to be recorded in a single acquisition [38]. The sensors had horizontal and vertical angles of view of 33.4 • and 25.5 • , respectively, and a fixed focal length of 7.5 mm, which corresponded to a ground-sampling interval of 47 mm for our flight mission at 100 m above ground level. The MAIA/Sentinel-2 system also had a standard GNSS receiver that synchronously logged the position and time at which the camera's shutter was activated. Throughout each of the UAV flights, a total of 18 MAIA/Sentinel-2 images were acquired over the trial plots and saved in a proprietary raw format with a 12-bit radiometric resolution.

Data Post-Processing
Image processing, including radial and radiometric calibration, was applied to the raw MAIA/Sentinel-2 imagery using the MultiCam Stitcher Pro software (v.1.1.8). The multiband images were first co-registered in order to correct the offsets between the nine sensors on the camera system. The radial calibration then involved a per-pixel correction for vignetting (i.e., the effect of a reduction in illumination from the centre to the edge of the image). To calculate ground-leaving reflectance, two identical white ground targets were placed at opposite ends of the UAV flight extents. Each target comprised a 1.0 m 2 panel that was made of a lambertian reflectant polyvinyl chloride (PVC)-coated material ('Odyssey' trademark material, Kayospruce Ltd.), which has previously been used in ESA remote sensing fieldwork campaigns [39]. At the beginning of each UAV flight, the spectral reflectance of these targets was measured using an Analytical Spectral Devices (ASD) FieldSpec Pro. This target reflectance was converted to absolute reflectance, following the guidelines outlined by the National Environment Research Council (NERC) Field Spectroscopy Facility [40]. These data were then convolved with the spectral response of each MAIA/Sentinel-2 band to correct the reflectance digital number (DN) recorded at pixels for each of the MAIA/Sentinel-2 spectral bands to absolute reflectance [37]: where for each MAIA/Sentinel-2 band the corrected DN, DN i , is calculated based on the ratio of the reference target spectra, Rt i , to the same value measured by the camera, Pt i . We applied corrections to the recorded position of each of the MAIA/Sentinel-2 images by matching the GNSS time-stamp to that of the UAV platform. Consequently, with the time-stamps matched, we were able to use the more precise real-time kinematic GNSS position recorded in the UAV flight log. The collected images were then loaded into Agisoft PhotoScan Professional (v.1.3.3) where a photogrammetric workflow was applied to align and produce an othomosaic of the 18 multiband images covering the trail plots with a ground sampling resolution of 0.04 m.
The multiband MAIA/Sentinel-2 orthomosaics were overlaid with vector polygons for each of the 50 winter wheat trials plots. A buffer of −0.5 m was applied to the polygon edges in order to ensure a representative coverage of the trail plots. For each plot, the vector dataset was then used to extract mean pixel values recorded in the MAIA/Sentinel-2 dataset. Since band 1 of the MAIA/Sentinel-2 data is measured by the Sentinel-2 multispectral instrument at a 60 m spatial resolution (Table 2), we considered this resolution to be too coarse for precision agricultural applications and, therefore, we omitted this band and reduced the analysis to the eight remaining bands.

Band Analysis and Model Evaluation Approaches
The analysis of the eight MAIA/Sentinel-2 bands was carried out using the machine learning algorithm of Gaussian processes regression (GPR; [41]). Specifically, we applied GPR in order to determine the most informative MAIA/Sentinel-2 bands for the retrieval of LAI, LCC and CCC from the mean of the spectral data extracted at each trial plot. GPR provides a non-parametric and probabilistic modelling approach to establishing relationships between inputs (i.e. MAIA/Sentinel-2 bands) and outputs (vegetation variables), allowing for both the predictive mean and variance to be obtained [20,41,42]. Studies have demonstrated the calibration of GPR for the estimation of biophysical variables from satellite and airborne sensors [43,44] and has been shown to perform favourably in comparison to alternative machine learning algorithms [45].
Using the ARTMO (automated radiative transfer models operator) machine learning regression algorithms (MLRA) toolbox [46], we applied GPR as a scaled Gaussian kernel function (Equation (2); [20]): where for a given covariance function relating two observations, k x i , x j , the GPR model hyper parameters include the scaling factor, v, a standard deviation describing the variance of the estimates, σ n , and the length-scale, σ b , for each of the MAIA/Sentinel-2 bands, b. These hyper parameters, along with the model weight, are automatically optimised by maximising the marginal likelihood when training the GPR model using the MAIA/Sentinel-2 dataset and corresponding ground measurements for LAI, LCC and CCC. The inverse of σ b represents the importance of each spectral band on k; accordingly, a higher σ −1 b value indicates a higher information content when developing a GPR model for estimating a variable of interest.
To quantify the sensitivity of the MAIA/Sentinel-2 bands on the GPR model estimates, we trained the model using all 250 ground measurements recorded during the 2018 field campaign at the experimental trail plots and applied it within the ARTMO MLRA sequential backward band removal algorithm. This band removal algorithm entails an iterative procedure, whereby a GPR model is first developed using all eight MAIA/Sentinel-2 input bands. The least informative band (i.e., lowest σ −1 b ) is then removed and a new GPR model is developed with the remaining bands, with this process repeating until the single most sensitive band remains. At each iteration of the backward band removal, we used the same input data for training and validating the GPR model. However, in order to ensure a robust analysis of each band, we applied a three-fold cross-validation sampling scheme. The GPR band analysis results, therefore, included the mean and standard deviation of the cross-validation statistics, including the coefficient of determination (R 2 ) and the NRMSE.
From the statistical outputs of the sequential backward band removal procedure, we analysed the GPR models for each of the crop variables based on Akaike information criterion (AIC) [47]. AIC is an approach for model selection based on relative performance and works by balancing the trade-offs Remote Sens. 2019, 11,2050 8 of 18 between model complexity (i.e., number of MAIA/Sentinel-2 bands) and goodness-of-fit against the validation data. We calculated the AIC based on the R 2 value for each GPR model [48]: where the AIC of a GPR model, AIC gpr_n , is calculated based on the maximum likelihood of the parameter vector, L θ , comprised of a number of bands, K. For each variable, the GPR model with the lowest AIC value was selected for further analysis. In order to determine the impact of uncertainty on GPR model development, this band analysis procedure was repeated both with and without the observational uncertainty. Our best performing GPR models for estimating LAI, LCC and CCC from the Sentinel-2 spectral data, with the propagated observational uncertainty, were further made available for the remote sensing community in a format that can be imported into the ARTMO MLRA toolbox (see the Supplementary Material).
In addition to the band analysis using the three-fold cross-validation statistics, we performed an independent validation of the GPR model performance by re-training the model with only half of the observations (i.e., with remaining observations used for validation). We further compared this independent model evaluation to a multivariate linear model, using ordinary least squares regression, which was developed using the most explanatory MAIA/Sentinel-2 bands and calibrated and validated using the same observations as the GPR model. We also quantified the spectral responses of these individual bands to the variables that were measured directly (i.e. LAI and LCC) and, in doing so, determined the extent to which these variables can be retrieved from the single bands. We thus compared the performance of the GPR approach to simple parametric models.

Uncertainty Analysis of In Situ Measurements
An overall high agreement existed between the in situ non-destructive measurements and destructive sample analysis ( Figure 2). The R 2 was 0.80 for a linear fit between the SunScan and destructive sample LAI measurements. We further corrected the bias in the SunScan measurements based on the destructive data, which reduced the NRMSE from 50% to 17%. A quadratic fit demonstrated a high correlation (R 2 = 0.75) between the LCC measurements and leaf N content derived from the destructive sample analysis. The variance of the measured LCC was, however, relatively high when compared to that of the LAI, with the standard deviation of SPAD measurements ranging from +/-3 to +/-16.

Sentinel-2 Band Analysis and Responses
Overall, the GPR model produced a high agreement with the LAI, LCC and CCC observations for all MAIA/Sentinel-2 band combinations ( Figure 3). For the most explanatory band combinations, the mean R 2 was 0.89 and 0.83 without and with the propagated uncertainty, respectively. The model performance was greatest when estimating the CCC with the R 2 ranging from 0.92 (without uncertainty) to 0.86 (with uncertainty). On the basis of the standard deviations of R 2 and RMSE values, with the sequential band removal, the estimation uncertainty for each of the variables remained relatively stable until the successive removal of bands from the most explanatory band combinations, where error and uncertainty of these estimates sharply increased when using less than the optimum number of spectral bands.

Sentinel-2 Band Analysis and Responses
Overall, the GPR model produced a high agreement with the LAI, LCC and CCC observations for all MAIA/Sentinel-2 band combinations ( Figure 3). For the most explanatory band combinations, the mean R 2 was 0.89 and 0.83 without and with the propagated uncertainty, respectively. The model performance was greatest when estimating the CCC with the R 2 ranging from 0.92 (without uncertainty) to 0.86 (with uncertainty). On the basis of the standard deviations of R 2 and RMSE values, with the sequential band removal, the estimation uncertainty for each of the variables remained relatively stable until the successive removal of bands from the most explanatory band combinations, where error and uncertainty of these estimates sharply increased when using less than the optimum number of spectral bands.
In comparing the GPR band selection with and without propagating observational uncertainty, it was found that without uncertainty, the number of bands included in the optimum bands was four for each variable, whereas the number of key bands varied from two to three when including uncertainty (Figure 3). The MAIA/Sentinel-2 red-edge band at 705 nm and near-infrared 783 nm band were frequently included in the most explanatory band combinations. The selected band combinations for the estimation of LAI and CCC were identical-comprising two red-edge and two near-infrared bands (705, 740, 783 and 865 nm) for the GPR models without the observational uncertainty. For the models including uncertainty, however, the two most explanatory bands selected for estimating LAI and CCC were that of the red-edge (705 nm) and near-infrared 865 nm bands only.
For the most explanatory bands selected from the GPR framework developed with the propagated uncertainty, we analysed the sensitivity of these individual bands for estimating LAI and LCC based on the spectral responses to these two variables ( Figure 4 and Table 3). The red-edge band at 705 nm showed a general exponential decrease in reflectance with increasing LAI (R 2 = 0.61), whereas the near-infrared band (865 nm) was characterised by a linear increase in reflectance with increasing LAI (R 2 = 0.67). For individual band responses to LCC, both the blue (490 nm) and green (560 nm) band reflectance exhibited a weak linear correlation (R 2 was 0.46 and 0.55 for the blue and green bands, respectively) with the reflectance in these bands decreasing with increasing LCC. Reflectance in the red-edge band (783 nm), however, showed a reasonable linear positive correlation (R 2 = 0.61) with increasing LCC. In comparing the GPR band selection with and without propagating observational uncertainty, it was found that without uncertainty, the number of bands included in the optimum bands was four for each variable, whereas the number of key bands varied from two to three when including uncertainty (Figure 3). The MAIA/Sentinel-2 red-edge band at 705 nm and near-infrared 783 nm band were frequently included in the most explanatory band combinations. The selected band combinations for the estimation of LAI and CCC were identical-comprising two red-edge and two near-infrared bands (705, 740, 783 and 865 nm) for the GPR models without the observational uncertainty. For the models including uncertainty, however, the two most explanatory bands selected for estimating LAI and CCC were that of the red-edge (705 nm) and near-infrared 865 nm bands only.
For the most explanatory bands selected from the GPR framework developed with the propagated uncertainty, we analysed the sensitivity of these individual bands for estimating LAI and LCC based on the spectral responses to these two variables ( Figure 4 and Table 3). The red-edge band at 705 nm showed a general exponential decrease in reflectance with increasing LAI (R 2 = 0.61), whereas the near-infrared band (865 nm) was characterised by a linear increase in reflectance with increasing LAI (R 2 = 0.67). For individual band responses to LCC, both the blue (490 nm) and green (560 nm) band reflectance exhibited a weak linear correlation (R 2 was 0.46 and 0.55 for the blue and green bands, respectively) with the reflectance in these bands decreasing with increasing LCC. Reflectance in the red-edge band (783 nm), however, showed a reasonable linear positive correlation (R 2 = 0.61) with increasing LCC.  Figure 3. Sentinel-2 band analysis using Gaussian processes regression (GPR) modelling trained using multi-date wheat observation, both without (left) and with (right) accounting for uncertainty in observations, for deriving Leaf Area Index (LAI), leaf chlorophyll content (LCC) and canopy chlorophyll content (CCC). Statics include the mean and standard deviations of the coefficient of determination (R 2 ) and normalised root-mean-square-error (NRMSE) from a three-fold cross-validation of the corresponding GPR models. The best performing GPR model, selected based on the lowest Akaike information criterion (AIC) value, is shown in boldface. Table 3. Summary of regression analysis statistics, including the coefficient of determination (R 2 ) and normalised root-mean-square-error (NRMSE), from comparing single-band and multivariate linear regression and Gaussian processes regression (GPR) modelling for retrieving ground measurements of LAI and LCC.

Independent Model Evaluation
The GPR model estimates when using the most explanatory band combinations demonstrated a high agreement to observations when validated using an independent dataset-the mean R 2 was 0.77 and the NRMSE was 12% ( Figure 5 and Table 3). The GPR model performance was very similar when deriving the LAI and CCC estimates, with the R 2 being 0.85 for both variables and the NRMSE being 9 and 10% for LAI and CCC, respectively. In comparison to the other two variables, the model performance was weaker when estimating LCC, with R 2 and NRMSE being 0.60 and 18%, respectively.
When compared to the GPR model estimates, a multivariate linear regression model using the same dataset for calibration and validation showed a generally weaker performance in the estimation of all three variables, with the R 2 ranging from 0.67 to 0.73 ( Figure 6). The LCC estimates by the linear regression model, however, had a higher agreement to the observations, with the R 2 and NRMSE being 0.67 and 13%, respectively.
The GPR model estimates when using the most explanatory band combinations demonstrated a high agreement to observations when validated using an independent dataset-the mean R 2 was 0.77 and the NRMSE was 12% ( Figure 5 and Table 3). The GPR model performance was very similar when deriving the LAI and CCC estimates, with the R 2 being 0.85 for both variables and the NRMSE being 9 and 10% for LAI and CCC, respectively. In comparison to the other two variables, the model performance was weaker when estimating LCC, with R 2 and NRMSE being 0.60 and 18%, respectively. When compared to the GPR model estimates, a multivariate linear regression model using the same dataset for calibration and validation showed a generally weaker performance in the estimation of all three variables, with the R 2 ranging from 0.67 to 0.73 ( Figure 6). The LCC estimates by the linear regression model, however, had a higher agreement to the observations, with the R 2 and NRMSE being 0.67 and 13%, respectively.

Ground Measurement Analysis and Uncertainty Characterisation
We quantified the uncertainty and bias of non-destructive measurements of LAI and LCC by comparisons to data derived from the destructive sub-sample analysis (Figure 2). The linear relationship established between the SunScan LAI and destructive measurements was biased (NRMSE = 62%), which was likely to have been due to the SunScan measurements over-estimating the LAI by including the contribution of stems and branches [49]. A non-linear relationship was respectively. When compared to the GPR model estimates, a multivariate linear regression model using the same dataset for calibration and validation showed a generally weaker performance in the estimation of all three variables, with the R 2 ranging from 0.67 to 0.73 ( Figure 6). The LCC estimates by the linear regression model, however, had a higher agreement to the observations, with the R 2 and NRMSE being 0.67 and 13%, respectively.

Ground Measurement Analysis and Uncertainty Characterisation
We quantified the uncertainty and bias of non-destructive measurements of LAI and LCC by comparisons to data derived from the destructive sub-sample analysis ( Figure 2). The linear relationship established between the SunScan LAI and destructive measurements was biased (NRMSE = 62%), which was likely to have been due to the SunScan measurements over-estimating the LAI by including the contribution of stems and branches [49]. A non-linear relationship was

Ground Measurement Analysis and Uncertainty Characterisation
We quantified the uncertainty and bias of non-destructive measurements of LAI and LCC by comparisons to data derived from the destructive sub-sample analysis ( Figure 2). The linear relationship established between the SunScan LAI and destructive measurements was biased (NRMSE = 62%), which was likely to have been due to the SunScan measurements over-estimating the LAI by including the contribution of stems and branches [49]. A non-linear relationship was observed between the LCC (SPAD) with increasing N content. This result is in agreement with research by Rostami et al. [50], where a quadratic plateau was observed between SPAD data with higher levels of N content in maize. In comparison to the LAI measurements, the LCC measurements made at each plot and sampling date were less consistent, which could be attributed to a high variability in the distribution of chlorophyll at both the leaf and canopy level [11,51].

Sentinel-2 Bands and GPR Modelling for Parameter Retrievals
Using UAV-based Sentinel-2 band observations and ground measurements, we assessed the capacity of the Sentinel-2 spectral bands for the retrieval of winter wheat LAI, LCC and CCC. We acknowledge that our GPR models for retrieving each of these variables were limited to one winter wheat season at our experimental field site. The calibration of the models was, however, carried out across multiple dates and N treatments; thus, covering a range of developmental stages and nutrient stresses. We would, therefore, expect the models to be broadly applicable when applied to retrieve the key variables from similar winter wheat varieties at alternative sites.
We exploited the band-ranking capabilities of a GPR machine learning algorithm and, in doing so, we objectively identified the most sensitive Sentinel-2 bands for retrieving each of these wheat variables. When we included the uncertainty of the ground measurements used in the training of the GPR models, it was found that two to three bands were sufficient for estimating the variables (Figure 3). Specifically, for the estimation of LAI and CCC without propagating the uncertainty, the optimum GPR model included both the red-edge bands (705 and 740 nm) and two of the near infra-red bands (783 and 865 nm). When propagating the uncertainty, the most informative bands included the red-edge band at 705 nm and near infra-red band at 865 nm only. This result suggests overfitting of the GPR model when variance in the training data is not accounted for and, thus, highlights the importance of propagating uncertainty in ground measurements when developing remote sensing retrieval algorithms. Being able to provide reliable estimates of the variables from a small number of well-defined Sentinel-2 bands also has the advantage of a reduction in processing time when generating high-resolution empirical retrieval maps of the variables over large areas [20]. More broadly, the ability to estimate multiple variables of a specific crop type using the essential spectral bands also has implications for the development of compact and lightweight UAV multispectral cameras [52].
Where the red-edge and near infra-red bands were favoured for the estimation of LAI and CCC, for the retrieval of LCC when the GPR model was developed with the propagated uncertainty, the optimum bands were composed of the visible blue (490 nm), green (560 nm) and near-infrared bands (783 nm). The sensitivity of the green and infrared bands for the estimation of LCC has previously been demonstrated [9,27]. Specifically, the Green Chlorophyll Vegetation Index (near infra-red/green; [31]) has successfully been applied in several studies for deriving crop chlorophyll content [24,[53][54][55]. Research by Wang et al. [56] involved the analysis of winter wheat spectral reflectance under different N applications and demonstrated that bands centred around the green and near-infrared spectral regions were sensitive to the treatments, whereas the blue band was comparatively less sensitive. In our analysis, the inclusion of the blue band in the most informative band combination for LCC was, therefore, unexpected, but did appear to improve the performance of the GPR model when compared to using the green and near-infrared bands only. Our analysis also demonstrated a relationship between blue reflectance and LCC measurements (Figure 4).
The traditional approaches for retrieving biophysical parameters from Earth observation data have involved VIs, whereas our GPR machine learning approach allows a more thorough exploration of the sensitivity of each band to the measured variable without making assumptions. From applying the GPR model within the band analysis framework, the most informative Sentinel-2 bands that we identified were also situated around spectral regions of known reflectance and absorption and had a reasonable relationship to the ground measurements when evaluated individually ( Figure 4). The sensitivity of these spectra to the biophysical variables is also broadly in agreement with past research [11,25,32,33]. We acknowledge, however, that our calibrated GPR modelling approaches are based on near-Earth (i.e., UAV) observations and, when applied to satellite-based observations, the retrieval accuracy would be subject to atmospheric effects. Research in Clevers et al. [24] has demonstrated a good agreement between atmospherically-corrected Sentinel-2 reflectance and ground-based radiometer measurements for potato crops. Nonetheless, we recommend additional testing of the robustness and a quantification of the uncertainty of the retrieval algorithms when applied to data acquired from the Sentinel-2 platform. Using an independent dataset for calibration and validation and the most informative bands, we also demonstrated that a simpler model approach of multivariate linear regression could provide generally comparable results to the GPR model ( Figures 5 and 6). Verrelst et al. [57] discusses the computation limitations of applying GPR models to large datasets and, therefore, the simpler regression approach may be more practical when retrieving pixel-level estimates of the variables across regional and country scales.

Potential of Sentinel-2 for Supporting Agricultural Management
The recent availability of open-access data from the Sentinel-2 dual-satellite constellation provides opportunities for the development and improvement of spatial data products that can support precision agriculture. Where past Earth observation satellite sensors (including SPOT-6/7 and Landsat-7/8) have relied on observations in the visible and near-infrared wavebands, this research has demonstrated that the two Sentinel-2 red-edge bands were frequently included within the most sensitive band combinations for retrieving key crop variables.
We have demonstrated the potential of the Sentinel-2 bands for providing more accurate estimates of LAI, which would be of value for improving the efficiency of crop model-data assimilation approaches, particularly when spatially upscaling model estimates from fields to regional extents [3,58]. In this research, we have shown a high correlation (R 2 = 0.86) between LCC and leaf N along with a high retrieval accuracy of LCC and CCC using Sentinel-2 bands. Being able to provide reliable estimates of LAI and chlorophyll content, i.e., as a proxy for N, is particularly useful for farmers when deciding on mid-season N fertilisation applications [18]. The traditional uniform approaches to fertiliser N applications are economically and environmentally inefficient, since they inherently ignore spatial heterogeneities in topography and soil properties [59,60]. The availability of timely and spatially explicit estimates of crop N content could be a crucial input for variable rate applications concerned with optimum N usage. Because of links between leaf N content and crop yields [50,61], within-field estimates of N could also be used to estimate crop yields. For use on an operational basis, the overall extent to which Sentinel-2 data can reliably support precision agricultural applications is, of course, dependent on the frequency of available cloud-free imagery. For our field site and growing season, we identified 19 cloud-free images. Since the period of time between successive observations ranged from 7 to 35 days, with an average of 17 days, this would not meet the biweekly observations recommended in previous research for tracking the temporal dynamics of crop growth [22,23]. However, these cloud-free images were within five days of each of the five growth stages targeted in this research and would be expected to be of value, particularly if the retrieved variables were used to update daily process-based crop models estimates.
The traditional approaches for variable retrieval use the visible red and near-infrared bands, which correspond to Sentinel-2 bands 4 and 8 that have a spatial resolution of 10 m (Table 2). On the other hand, the Sentinel-2 red-edge bands have a spatial resolution of 20 m. Although it is clear from this research that the red-edge bands improve estimates of LAI, LCC and CCC, for practical applications, the use of these bands would be at the expense of a reduction in spatial resolution. Löw and Duveiller [62] sampled from a continuum of increasingly coarser pixel sizes in order to investigate spatial resolution requirements for the image classification of different crop types. This research demonstrated that pixel sizes of around 117 m were sufficient to identify winter wheat, but this was found to be dependent on landscape heterogeneity (i.e., field size and shape) and timing within the growing season. Further research by Colombo et al. [63] showed stability in an LAI-VI relationship when spatially aggregating IKONOS satellite imagery pixels from 12 to 36 m. In our research we would, therefore, not anticipate any substantial increases in uncertainty when retrieving the crop variables at 20 m; the increase in performance when using the Sentinel-2 red-edge bands are likely to outweigh any negative impacts of a reduction in spatial resolution. Nonetheless, where our study focuses on the spectral characteristics of Sentinel-2, we would recommend future research related to Sentinel-2 spatial resolution that includes tracking the propagation of uncertainty when estimating variables using Sentinel-2 data at 10 and 20 m resolution for specific crop types.

Conclusions
This study has evaluated the Sentinel-2 satellite multispectral instrument spectral bands for the estimation of winter wheat variables-LAI, LCC and CCC-required for supporting precision agricultural technologies. Where past research has often used synthetic Sentinel-2 data within the growing season, here, we used data from a UAV-mounted multispectral camera with sensors matching the key Sentinel-2 wavebands. The acquisition of UAV multispectral data was carried out in coordination with ground measurements. These measurements, comprising destructive and non-destructive sample analysis data, were used to calibrate and validate the performance of a GPR machine-learning algorithm we applied to identify the most informative spectral bands for estimating each of the biophysical variables. The ground measurements were also used to quantify uncertainty that was propagated into the GPR model training data.
Overall, we have demonstrated a high retrieval accuracy of the variables when using the most informative Sentinel-2 bands (mean R 2 = 0.86). The Sentinel-2 red-edge and near-infrared bands were identified as being the most informative, particularly for LAI and CCC. The propagation of uncertainty in the ground measurement reduced the number of most informative bands, indicating an overfitting of the GPR model when uncertainty was not properly accounted for.
In comparison to previous satellite missions, the results we present highlight the potential of Sentinel-2 spectral data within an operational farm-scale decision support system. Future research should include testing the robustness and characterizing the uncertainty of the GPR modelling approach when applied to data acquired from the Sentinel-2 platform, including the uncertainty linked to the spatial scaling of these estimates to the resolution of Sentinel-2 multi-spectral instrument (i.e., 10 and 20 m). Furthermore, the use of GPR modelling, which provides a predictive mean and variance, would be an ideal approach for investigating the propagation of uncertainty from ground measurements to the scale of the Sentinel-2 sensor.