Using UAV Collected RGB and Multispectral Images to Evaluate Winter Wheat Performance across a Site Characterized by Century-Old Biochar Patches in Belgium

: Remote sensing data play a crucial role in monitoring crop dynamics in the context of precision agriculture by characterizing the spatial and temporal variability of crop traits. At present there is special interest in assessing the long-term impacts of biochar in agro-ecosystems. Despite the growing body of literature on monitoring the potential biochar e ﬀ ects on harvested crop yield and aboveground productivity, studies focusing on the detailed crop performance as a consequence of long-term biochar enrichment are still lacking. The primary objective of this research was to evaluate crop performance based on high-resolution unmanned aerial vehicle (UAV) imagery considering both crop growth and health through RGB and multispectral analysis, respectively. More speciﬁcally, this approach allowed monitoring of century-old biochar impacts on winter wheat crop performance. Seven Red-Green-Blue (RGB) and six multispectral ﬂights were executed over 11 century-old biochar patches of a cultivated ﬁeld. UAV-based RGB imagery exhibited a signiﬁcant positive impact of century-old biochar on the evolution of winter wheat canopy cover ( p -value = 0.00007). Multispectral optimized soil adjusted vegetation index indicated a better crop development over the century-old biochar plots at the beginning of the season ( p -values < 0.01), while there was no impact towards the end of the season. Plant height, derived from the RGB imagery, was slightly higher for century-old biochar plots. Crop health maps were computed based on principal component analysis and k-means clustering. To our knowledge, this is the ﬁrst attempt to quantify century-old biochar e ﬀ ects on crop performance during the entire growing period using remotely sensed data. Ground-based measurements illustrated a signiﬁcant positive impact of century-old biochar on crop growth stages ( p -value of 0.01265), whereas the harvested crop yield was not a ﬀ ected. Multispectral simpliﬁed canopy chlorophyll content index and normalized di ﬀ erence red edge index were found to be good linear estimators of harvested crop yield ( p -value (Kendall) of 0.001 and 0.0008, respectively). The present research highlights that other factors (e.g., inherent pedological variations) are of higher importance than the presence of century-old biochar in determining crop health and yield variability.


Introduction
Biochar is a carbon-rich material obtained through the pyrolysis of feedstocks to intentionally amend soil [1,2]. In the context of climate-smart agriculture, special attention has been given to biochar application owing to its supposed ability to improve soil fertility and long-term organic carbon storage [3][4][5][6]. To date, short-term biochar effects have been well addressed by studying the soil nutrient availability [7,8], soil water holding capacity [9][10][11], and plant available water content [12,13]. The aforementioned improvements of soil quality could consequently lead to an impact on the crop productivity [14]. Moreover, biochar was found to increase nutrient uptake and crop yield while being incorporated into soil matrix together with fertilizers [15]. The short-term improvement of crop productivity can be attributed to the addition of nutrients related to the biochar application [7]. However, long-term biochar effects on soil and agro-ecosystems properties are still poorly understood. Several authors have reported a higher nutrient availability in the century-old biochar sites originating from earth kiln sites [16,17]. The higher nutrient availability in long-term biochar-enriched soils can be explained by cation exchange capacity that increases over time [18]. Moreover, century-old biochar significantly affects physico-chemical properties of soils in relic hearths of agricultural fields [16].
Across the globe, many efforts have been devoted to unravel biochar effects on aboveground biomass productivity. In many of these studies, a short-term effect of biochar on increasing crop productivity was reported [7,19,20], however, Güereña et al. [21] and van Zwieten et al. [22] showed a reductive influence. Similarly, Hernandez-Soriano et al. [23] reported a long-term biochar impact on increasing aboveground productivity, whereas Mastrolonardo et al. [24] observed a decreased tree ring width in relic charcoal hearths of Belgium. Oak growth depression was also found in aged charcoal hearth patches in the United States [25].
The studies listed above all explored the biochar effects on crop productivity. However, it is highly likely that biochar could affect crop development only during the green-up phase, while the harvested crop yield could remain unaffected [26]. More specifically, long-term biochar effects on crop performance and its development over the entire growing season have received much less attention than biochar effects on total crop productivity.
In recent years, remote sensing has been promoted as a means in precision agriculture for detailed spatio-temporal monitoring of crop dynamics [27,28]. A particular focus at present is the use of satellite imagery for precision agricultural practices due to the images' fine spatial, spectral, and temporal resolution [28]. Moreover, high-resolution unmanned aerial vehicle (UAV)-based imagery can be used to address the data gap in satellite products suffering from cloud coverage [29]. Recent advances in Red-Green-Blue (RGB) and multispectral UAV remote sensing have improved the ability to accurately monitor crop performance and vegetation characteristics at the canopy level, such as canopy cover, leaf area index, and chlorophyll content, allowing for site-specific decision making in the context of precision agriculture [29][30][31][32]. In addition, due to remote UAV flights, the experimental surface is not damaged [33], thus enabling a complete crop monitoring over the entire season. Therefore, the spatio-temporal information provided by the UAV images can reveal variability in crop performance due to the presence of century-old biochar.
Every year, agricultural products suffer significant economic losses caused by multiple pests and diseases, such as yellow rust or aphid, which may lead to a shortage of food production globally [34,35] and particularly in Europe [36]. Cameras onboard UAVs are able to illustrate changes in crop dynamics as a function of a particular disease at early onset [37]. The pioneering work of Franke et al. [38] reported that pathogens reduces plant chlorophyll content in the visible and red-edge spectral bands due to necrotic and chlogotic lesions. Su et al. [37] highlighted that the red band of the spectrum has strong potential in discriminating healthy from diseased winter wheat plants. Su et al. [37] also showed that the multispectral optimized soil adjusted vegetation index (OSAVI) has a strong discriminating Remote Sens. 2020, 12, 2504 3 of 22 capability between healthy and diseased winter wheat plants. In addition, UAV imagery is able to detect changes in crop traits linked to disease incidence [39]. Therefore, high-resolution UAV remote sensing can be successfully applied for spatial explicit assessment of crop health as a consequence of the presence of century-old biochar within agricultural soils.
The primary objective of this paper is to evaluate crop performance by studying crop growth and health through the usage of high-resolution RGB and multispectral information, respectively, derived from two UAV platforms. This approach also allowed us to investigate, as a secondary objective, the impacts of century-old biochar on winter wheat crop performance. The available imagery sources were RGB and multispectral data collected using two UAV platforms. UAV remote sensing data in combination with ground-based measurements of crop traits were used to determine the influence of century-old biochar from kiln sites of an agricultural field on winter wheat crop growth. The present study evaluated, for the first time, the impact of century-old biochar on crop performance at the canopy level using UAV imagery, while paving the way to quantify the long-term effects of biochar on winter wheat crop health.

Site Description and Data Acquisition
The data were collected in a farm of approximately 13 hectares cultivated with winter wheat located in Isnes, Belgium ( Figure 1). Winter wheat (Triticum aestivum L.) seeds were sown on 6 December 2018. The region is of temperate climate with dry and hot summers. The soil type is a Luvisol characterized by a silt loam texture [40]. Homogenous agricultural treatments were applied throughout the field. The detailed description is available in Heidarian Dehkordi et al. [26], in which the same experimental configuration was followed by eleven 10 × 10 m century-old biochar plots (red squares in Figure 1) together with the corresponding reference plots (blue squares in Figure 1). Ground control point (GCP) targets, with known centroid coordinates through a real-time kinematic GPS, were deployed throughout the farm (Figure 1) to help with accurate geo-referencing of the data (i.e., 3 cm horizontal accuracy).
Two remote sensing datasets (i.e., RGB and multispectral) were acquired during the 2019 growing season (Table 1) using high-resolution sensors onboard two UAV platforms to monitor the crop status over time ( Figure 2). RGB images were collected using a DJI Phantom 4 Pro and the multispectral images were obtained using a MicaSense RedEdge-M and downwelling light sensor, onboard a DJI Matrice 100 platform specifically designed by Heidarian Dehkordi et al. [26], capturing blue (475 nm central wavelength), green (560 nm), red (668 nm), red-edge (717 nm), and near-infrared (840 nm) spectral channels with bandwidths of 20, 20, 10, 10, and 40 nm, respectively. The UAV platforms were programmed to fly at 65 and 60 m above ground altitude (AGL) with an airspeed of 5 and 7 ms −1 for RGB and multispectral sensors, respectively. The forward and sideway image overlaps were 85% and 75% for RGB, and 75% and 75% for multispectral acquisitions, reaching a ground sampling distance (GSD) of 1.8 and 3.7 cm for RGB and multispectral images, respectively.  Two remote sensing datasets (i.e., RGB and multispectral) were acquired during the 2019 growing season (Table 1) using high-resolution sensors onboard two UAV platforms to monitor the crop status over time ( Figure 2). RGB images were collected using a DJI Phantom 4 Pro and the multispectral images were obtained using a MicaSense RedEdge-M and downwelling light sensor, onboard a DJI Matrice 100 platform specifically designed by Heidarian Dehkordi et al. [26], capturing blue (475 nm central wavelength), green (560 nm), red (668 nm), red-edge (717 nm), and near-infrared (840 nm) spectral channels with bandwidths of 20, 20, 10, 10, and 40 nm, respectively. The UAV platforms were programmed to fly at 65 and 60 m above ground altitude (AGL) with an airspeed of 5 and 7 ms −1 for RGB and multispectral sensors, respectively. The forward and sideway image overlaps were 85% and 75% for RGB, and 75% and 75% for multispectral acquisitions, reaching a ground sampling distance (GSD) of 1.8 and 3.7 cm for RGB and multispectral images, respectively.   The photogrammetric processing of the acquired UAV images was conducted in Pix4Dmapper photogrammetric software (Pix4D S.A., Lausanne, Switzerland). After the aerial triangulation of the UAV images, the geometric processing was optimized in Pix4Dmapper using the coordinates of the GCPs, and dense point clouds and final orthomosaic images were then created. For the multispectral data, the radiometric calibration was applied using MicaSense calibrated reference panels with known reflectance values (MicaSense, Seattle, WA, USA) taken immediately before and after each flight to convert the raw image digital number (DN) to the calibrated reflectance values. In addition, sensor's sensitivity to light, expressed as ISO, values recorded by the downwelling light sensor were used to aid the radiometric calibration. Both RGB and multispectral datasets were converted to GeoTIFF format with the projected geographic coordinates system WGS 1984 UTM Zone 31N.
In addition to the aforementioned remotely-sensed datasets, ground based measurements of crop traits were performed ( Figure 2). Crop growth stages were determined on 28 May 2019 according to the stabiologische bundesanstalt, bundessortenamt and chemical industry (BBCH) scale for winter wheat [41]. For each plot, a regular non-destructive sampling method of 5 plants, i.e., 4 in the corners and 1 in the middle to capture the intra-plot variability, was followed. The BBCH code of each plant was visually determined and the median value of the five plants represented the plot growth stage.
Yield measurements were conducted on 20 July 2019. For each experimental plot, an area of approximately 20 m 2 was harvested to determine the yield.

UAV-Based Crop Monitoring
Canopy cover was estimated from the acquired high-resolution RGB images ( Figure 2) based on the classification method by Heidarian Dehkordi et al. [26]. This algorithm classifies the pixels in the orthomosaic to vegetation and non-vegetation classes according to the spectral differences between the green band and the blue and red bands as below: where Blue, Green, and Red show the digital number (DN) of the corresponding spectral channel of the visible spectrum, α and β are the classification threshold values with α representing the spectral difference between the green and red bands and β revealing the spectral difference between the green and red spectral channels. Similar to the finding of Heidarian Dehkordi et al. [26], a threshold value of greater than 20 for both α and β was also reached in the present study through the visual inspection for the best vegetation classification. Canopy cover estimates were used as a potential discriminative parameter between the reference and century-old biochar plots in each acquisition date [26]. Plant height, as an important phenotyping crop trait, has received no attention in the previous studies exploring the impacts of century-old biochar on crop performance. In this study, plant height was estimated from the RGB images ( Figure 3) as below: where Height is the computed wheat height, DSM and DTM are the UAV-based digital surface model and digital terrain model, respectively, with index t representing the acquisition date. DSMs were generated from the RGB flights for each acquisition date. The RGB flight on 22 February 2019 (Table 1) was executed when the field was completely bare soil and, hence, the UAV-based DSM, revealing the pure terrain surface topography, was used as the DTM for this study. GCP_2 (Figure 1) was used as the reference levelling point to calibrate the generated UAV-based DSMs at each acquisition date, bringing the elevation values onto the same scale. The created RGB-based vegetation mask (Equation (1)), with GSD of 1.8 cm, was used to retrieve the DSM from the true canopy pixels for the height computation. In addition to the RGB derivatives listed above, narrow-band vegetation indices given in Table  2 were computed from the multispectral dataset ( Figure 2) as a proxy to crop traits allowing for demonstrating spectral changes related to crop structure [39]. Furthermore, the multispectral information derived from the vegetation indices were analyzed to study the impacts of century-old biochar on winter wheat crop spectral information. The normalized difference vegetation index (NDVI) is a well-known proxy of photosynthetic activity and vegetation greenness. Canopy cover and leaf area index (LAI) were accurately predicted based on the weighted difference vegetation index (WDVI) in several studies [42][43][44]. The normalized difference red edge index (NDRE) was shown as a good proxy of wheat nitrogen concentration [45]. The optimized soil adjusted vegetation index (OSAVI) was found to be a robust indicator of leaf chlorophyll content [46], and the same conclusion was reached for the chlorophyll vegetation index (CVI). The enhanced vegetation index (EVI) was shown to be sensitive to canopy structural parameters such as LAI while adjusting soil background effects and atmospheric resistance using the blue band of the spectrum [47]. The chlorophyll index red (CI-red) and the simplified canopy chlorophyll content index (s-CCCI) are known alternatives of structural indices which are often used for estimating canopy chlorophyll and nitrogen contents [48,49].
The aforementioned vegetation indices were used for identifying the crop spectral differences between reference and century-old biochar plots throughout the growing season. First, WDVI was used to generate a multispectral-based vegetation mask. A threshold of WDVI higher than 0.30 was used to retrieve the so-called true canopy pixels because this threshold yielded the best results through the visual inspection. This method was based on the approach adopted by Heidarian Dehkordi et al. [26]. Finally, the created WDVI-based mask was applied to the vegetation indices (Table 2) to remove the soil pixels from the analysis.
Unsupervised k-means clustering [50] was applied to the multispectral OSAVI imagery in order to generate a MSP crop health map with "good", "moderate", and "poor" classes ( Figure 2). The kmeans algorithm identifies the k number of cluster centroids (i.e., 3 in the present study) within the image and allocates each pixel to the nearest cluster by minimizing its Euclidean distance from the cluster centroid; this iterative algorithm applies until a steady state with the optimized centroid In addition to the RGB derivatives listed above, narrow-band vegetation indices given in Table 2 were computed from the multispectral dataset ( Figure 2) as a proxy to crop traits allowing for demonstrating spectral changes related to crop structure [39]. Furthermore, the multispectral information derived from the vegetation indices were analyzed to study the impacts of century-old biochar on winter wheat crop spectral information. The normalized difference vegetation index (NDVI) is a well-known proxy of photosynthetic activity and vegetation greenness. Canopy cover and leaf area index (LAI) were accurately predicted based on the weighted difference vegetation index (WDVI) in several studies [42][43][44]. The normalized difference red edge index (NDRE) was shown as a good proxy of wheat nitrogen concentration [45]. The optimized soil adjusted vegetation index (OSAVI) was found to be a robust indicator of leaf chlorophyll content [46], and the same conclusion was reached for the chlorophyll vegetation index (CVI). The enhanced vegetation index (EVI) was shown to be sensitive to canopy structural parameters such as LAI while adjusting soil background effects and atmospheric resistance using the blue band of the spectrum [47]. The chlorophyll index red (CI-red) and the simplified canopy chlorophyll content index (s-CCCI) are known alternatives of structural indices which are often used for estimating canopy chlorophyll and nitrogen contents [48,49].
The aforementioned vegetation indices were used for identifying the crop spectral differences between reference and century-old biochar plots throughout the growing season. First, WDVI was used to generate a multispectral-based vegetation mask. A threshold of WDVI higher than 0.30 was used to retrieve the so-called true canopy pixels because this threshold yielded the best results through the visual inspection. This method was based on the approach adopted by Heidarian Dehkordi et al. [26]. Finally, the created WDVI-based mask was applied to the vegetation indices (Table 2) to remove the soil pixels from the analysis.
Unsupervised k-means clustering [50] was applied to the multispectral OSAVI imagery in order to generate a MSP crop health map with "good", "moderate", and "poor" classes ( Figure 2). The k-means algorithm identifies the k number of cluster centroids (i.e., 3 in the present study) within the image and allocates each pixel to the nearest cluster by minimizing its Euclidean distance from the cluster centroid; this iterative algorithm applies until a steady state with the optimized centroid positions is reached [50]. Moreover, principal component analysis [51] was performed on the red and green spectral channels retrieved from the RGB imagery ( Figure 2). The first component (PC1) of the principal component analysis (PCA), which accounts for the maximum proportion of variance [52] from the red and green Remote Sens. 2020, 12, 2504 7 of 22 bands of the RGB orthomosaic (i.e. 99.06% in the present study), was used in the unsupervised k-means clustering to create a RGB crop health map ( Figure 2). More precisely, as recommended by Heidarian Dehkordi et al. [26], the need to explore the impact of century-old biochar on crop health during the senescence phase was addressed in this study.
Moreover, the use of multispectral vegetation indices for the remote sensing-based estimation of crop yield at the field-scale, as suggested by [53], was evaluated in this study ( Figure 2). In addition, the impact of century-old biochar on the predicted multispectral crop yield was studied.

Statistical Analysis
Statistical paired and global t-tests [59] were used to determine whether the contrasts between reference and century-old biochar plots were statistically significant. Several levels of significance were considered to investigate the observed differences. Throughout the manuscript, asterisks *, **, ***, and **** indicate the significance levels of p-values less than 0.05, 0.01, 0.001, and 0.0001, respectively.
Relationships between UAV-based derivatives and ground-based measurements were evaluated using the statistical F-test. Subsequently, the significance of the obtained correlations was evaluated through the non-parametric Kendall's rank correlation tau (τ) test which is less sensitive to the small statistical sample size, i.e., 11 experimental plots in the present study [60].

Canopy Cover
The temporal pattern of canopy cover shows vegetation development from 20 March to 29 April (Figure 4). At the beginning of May the canopy cover was stable around the maximum until mid-May. From this time onwards, canopy cover slightly decreased towards the end of the season since leaves lost their greenness and began yellowing. The canopy cover clearly showed larger values in biochar plots than the reference plots for most of the season (Figure 4). The Area under the curve (AUC) of the canopy cover was significantly higher for the century-old biochar plots than the reference plots over the season (p-value of 0.00007 and 0.00002 for paired and global t-tests, respectively). This finding is in accordance with the positive impact of century-old biochar patches on the evolution of chicory canopy cover in the same experimental farm [26]. Another study by Kerré et al. [6] also identified this positive impact of century-old biochar from relict charcoal hearths on maize crop growth in Belgium.
The greater canopy cover development over the century-old biochar plots could be explained by improved soil physical properties [6] such as water holding capacity (WHC), and soil porosity and infiltration rates, leading to higher water content available to plants that, in turn, increased the vegetation development in the biochar plots. Moreover, higher soil chemical properties (i.e., organic carbon and nitrogen content) were observed over the same century-old biochar plots [26] characterized by a positive relationship with the evolution of canopy cover in the 2018 growing season.
(AUC) of the canopy cover was significantly higher for the century-old biochar plots than the reference plots over the season (p-value of 0.00007 and 0.00002 for paired and global t-tests, respectively). This finding is in accordance with the positive impact of century-old biochar patches on the evolution of chicory canopy cover in the same experimental farm [26]. Another study by Kerré et al. [6] also identified this positive impact of century-old biochar from relict charcoal hearths on maize crop growth in Belgium. The greater canopy cover development over the century-old biochar plots could be explained by improved soil physical properties [6] such as water holding capacity (WHC), and soil porosity and infiltration rates, leading to higher water content available to plants that, in turn, increased the vegetation development in the biochar plots. Moreover, higher soil chemical properties (i.e., organic carbon and nitrogen content) were observed over the same century-old biochar plots [26] characterized by a positive relationship with the evolution of canopy cover in the 2018 growing season.
In addition, the earlier germination of winter wheat plants over the century-old biochar plots on the first acquisition (Figure 4) could be attributed to the dark background soils increasing the soil In addition, the earlier germination of winter wheat plants over the century-old biochar plots on the first acquisition (Figure 4) could be attributed to the dark background soils increasing the soil temperature and subsequently also plant temperature, leading to a faster crop phenological development. Hence, future research using high-resolution thermal imagery should investigate the impact of century-old biochar on soil-plant water exchanges and their relationship with crop development. Figure 5 illustrates the outcome of the height analysis. Winter wheat height increased from 20 March to 29 April ( Figure 5) which is in accordance with our earlier finding of canopy cover evolution (Figure 4). The height estimates were stable in May and June around the maximum. Wheat height was rather accurately derived with an accuracy (SD/mean) of 4.92 cm counting the entire dataset ( Figure 5), while the previous work by ten Harkel et al. [61] reached a higher accuracy of 3.4 cm for wheat height estimates using UAV-based light detection and ranging (LiDAR) data. The usage of UAV-based LiDAR enables a proper selection of the top of the canopy pixels to be used in height computations [61] increasing the modeling accuracy, which was not possible in the present study. The highest variation in height estimates ( Figure 5) was observed on 28 March (SD of 10.73 and 11.25 cm for the century-old biochar and reference plots, respectively), which could be attributed to the erectophile and open structure of winter wheat at the beginning of the season, lowering the number of true canopy pixels used for height estimates. However, the observed variations were much lower during the maturity phase (with a maximum SD of 5.45 cm on 13 May) because wheat canopies are rather dense, increasing the number of true canopy pixels in height computations. Figure 5 illustrates the outcome of the height analysis. Winter wheat height increased from 20 March to 29 April ( Figure 5) which is in accordance with our earlier finding of canopy cover evolution (Figure 4). The height estimates were stable in May and June around the maximum. Wheat height was rather accurately derived with an accuracy (SD/mean) of 4.92 cm counting the entire dataset ( Figure 5), while the previous work by ten Harkel et al. [61] reached a higher accuracy of 3.4 cm for wheat height estimates using UAV-based light detection and ranging (LiDAR) data. The usage of UAV-based LiDAR enables a proper selection of the top of the canopy pixels to be used in height computations [61] increasing the modeling accuracy, which was not possible in the present study. The highest variation in height estimates ( Figure 5) was observed on 28 March (SD of 10.73 and 11.25 cm for the century-old biochar and reference plots, respectively), which could be attributed to the erectophile and open structure of winter wheat at the beginning of the season, lowering the number of true canopy pixels used for height estimates. However, the observed variations were much lower during the maturity phase (with a maximum SD of 5.45 cm on 13 May) because wheat canopies are rather dense, increasing the number of true canopy pixels in height computations. Figure 5. Comparison of the plant height of all the 11 century-old biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median height, respectively. The error-bar presents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum height, respectively. Outside of the plot, asterisks *, **, ***, and **** reveal the statistical levels of significance; the acronym NS stands for nonsignificant. Figure 5. Comparison of the plant height of all the 11 century-old biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median height, respectively. The error-bar presents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum height, respectively. Outside of the plot, asterisks *, **, ***, and **** reveal the statistical levels of significance; the acronym NS stands for non-significant.

Plant Height
The impact of century-old biochar on winter wheat height considering all of the 11 biochar plots versus the reference plots is shown in Figure 5 for each acquisition date. The mean height (mean value of the 11 plots) was higher in the century-old biochar plots than in the reference plots for each acquisition date except 28 March ( Figure 5). It is worth noting that the observed mean height differences between reference and century-old biochar plots were statistically non-significant for all the monitoring dates ( Figure 5). This slight promotion of wheat height as a consequence of biochar presence was previously reported for lettuce-cabbage-lettuce [62] and oat [63] fields. The greater plant height over the century-old biochar patches could be attributed to the higher soil organic carbon and nitrogen contents as previously found by Schulz et al. [63]. Figure 6 shows the outcome of the height analysis for each experimental plot. Most of the experimental plots showed greater height values in century-old biochar plots compared to their corresponding reference plots throughout the season, except for plots 8 and 11 where an inverse trend was observed ( Figure 6). This can be explained by the location of the experimental plots 8 and 11 which were situated in a part of the field with a high soil moisture content (images not shown). The observed difference of plant height between the century-old biochar plots and the corresponding reference plots at each acquisition date was significant for all of the experimental plots ( Figure 6) except for plots 5, 8, and 10 on 28 March and plot 10 on 16 April. The observed height differences between reference and century-old biochar plots were significant for most of the experimental pairs separately (as shown in Figure 6), whereas the differences were not of particular concern when counting the entire dataset ( Figure 5). The latter could be explained by the within-field spatial variability of plant height that accounts for the spatial variation in field topographic derivatives such as elevation, groundwater flow, or surface wetness. Validation with the ground-based measurements of plant height should be performed in future studies to explore the robustness of the proposed plant height derivation method in this research.  Figure 7 illustrates the temporal profiles for OSAVI. There is a slight vegetation development from 20 March until the end of that month. The OSAVI then increased gradually until the end of April and remained constant at the maximum towards mid-May. From this time onwards, the OSAVI started to decrease as the leaves started to lose their greenness. The overall OSAVI pattern (Figure 7) is completely in line with the trend observed from the canopy cover evolution in Figure 4.  Figure 7 illustrates the temporal profiles for OSAVI. There is a slight vegetation development from 20 March until the end of that month. The OSAVI then increased gradually until the end of April and remained constant at the maximum towards mid-May. From this time onwards, the OSAVI started to decrease as the leaves started to lose their greenness. The overall OSAVI pattern (Figure 7) is completely in line with the trend observed from the canopy cover evolution in Figure 4. The impact of century-old biochar on OSAVI considering all of the 11 century-old biochar plots versus the reference plots for each acquisition date is shown in Figure 7. The mean OSAVI (mean value of the 11 plots) was higher in the century-old biochar plots than the reference plots at each acquisition date except on 24 June when the mean OSAVI was 0.73 for both reference and centuryold biochar plots (Figure 7). This finding is in accordance with the decreased plant greenness at the end of growing season for chicory in century-old biochar patches of the same experimental farm [26]. The observed OSAVI difference was non-significant between the reference and century-old biochar plots towards the end of the season, while it was highly significant from the beginning of the season until mid-April (Figure 7). This could be partially explained by the sparseness of the wheat canopies at the beginning of the season, increasing the interactions between plants and dark background soils in the century-old biochar plots.

Vegetation Indices
Contrast between reference and century-old biochar plots for each vegetation index is presented in Table 3 for all the acquisition dates. Century-old biochar significantly affected the crop spectral information derived from NDVI, NDRE, OSAVI, CI-red, and s-CCCI on 20 and 28 March. However, WDVI, CVI, and EVI did not exhibit significant differences in crop spectral status between reference and century-old biochar plots during March. All of the examined vegetation indices showed a significant impact of century-old biochar on crop spectral status on 16 April except WDVI, EVI, and s-CCCI (Table 3). On 29 April only NDVI, NDRE, and CI-red, and on 13 May only NDVI and CI-red, Figure 7. Comparison of the optimized soil adjusted vegetation index (OSAVI) of all of the 11 century-old biochar plots (red) versus the 11 reference plots for each acquisition date. The black circle and cross represent the mean and median OSAVI, respectively. The error-bar represents the corresponding standard deviation. The bottom and top black pluses (+) indicate the minimum and maximum OSAVI, respectively. Outside of the plot, asterisks *, **, ***, and **** reveal the statistical levels of significance and the acronym NS stands for non-significant.
The impact of century-old biochar on OSAVI considering all of the 11 century-old biochar plots versus the reference plots for each acquisition date is shown in Figure 7. The mean OSAVI (mean value of the 11 plots) was higher in the century-old biochar plots than the reference plots at each acquisition date except on 24 June when the mean OSAVI was 0.73 for both reference and century-old biochar plots (Figure 7). This finding is in accordance with the decreased plant greenness at the end of growing season for chicory in century-old biochar patches of the same experimental farm [26]. The observed OSAVI difference was non-significant between the reference and century-old biochar plots towards the end of the season, while it was highly significant from the beginning of the season until mid-April (Figure 7). This could be partially explained by the sparseness of the wheat canopies at the beginning of the season, increasing the interactions between plants and dark background soils in the century-old biochar plots.
Contrast between reference and century-old biochar plots for each vegetation index is presented in Table 3 for all the acquisition dates. Century-old biochar significantly affected the crop spectral information derived from NDVI, NDRE, OSAVI, CI-red, and s-CCCI on 20 and 28 March. However, WDVI, CVI, and EVI did not exhibit significant differences in crop spectral status between reference and century-old biochar plots during March. All of the examined vegetation indices showed a significant impact of century-old biochar on crop spectral status on 16 April except WDVI, EVI, and s-CCCI (Table 3). On 29 April only NDVI, NDRE, and CI-red, and on 13 May only NDVI and CI-red, exhibited a significant impact of century-old biochar on crop spectral status. All of the examined vegetation indices exhibited non-significant differences between reference and century-old biochar plots on 24 June (Table 3). This finding is in accordance with [26] who reported no significant difference of the multispectral vegetation indices between the century-old biochar and reference plots for chicory crop at the beginning of July. Table 3. Mean vegetation index values over the study area comparing the century-old biochar and reference plots during the growing season of 2019. The statistical significance of the difference is expressed as the p-value of the paired t-test. Next to the p-values asterisks, *, **, ***, and **** reveal the statistical levels of significance and acronym NS stands for non-significant. It is worth mentioning that the WDVI and the EVI performed the worst in discriminating spectral difference between century-old biochar and reference plots throughout the season. The latter could be attributed to the use of the blue band in EVI's computation which is not of particular importance in the vegetation spectral signature [64]. Figure 8 shows the outcome of the OSAVI analysis for each experimental plot. Most of the experimental plots exhibited higher OSAVI values in century-old biochar plots compared to their adjacent reference plots for most of the season (Figure 8). Graphs of the other examined multispectral vegetation indices (Table 2) are not presented here. No clear trend was observed in terms of the OSAVI difference between reference and century-old biochar plots (Figure 8). The contrast between reference and century-old biochar plots was highly significant at several acquisitions (e.g., for plot 1 on 29 April), whereas some acquisitions exhibited a non-significant difference (e.g., plot 1 on 28 March).

Crop Health
Both canopy cover and OSAVI were characterized by a decreasing trend on 24 June (Figures 4  and 7, respectively). The impact of century-old biochar on crop health during the senescence phase is presented in Figure 9. In addition, Figure 9 presents the two different crop health maps, obtained from the two different sensors, including the associated general methodological workflows as represented in Figure 9a,b for the multispectral sensor and Figure 9c,d for the RGB sensor. An unsupervised k-means clustered crop health map of OSAVI including "good", "moderate", and "poor" classes, corresponding to 24 June, is presented in Figure 9b. In addition, an unsupervised kmeans clustered crop health map derived from the first component (PC1) from the principal component analysis of the red and green bands of the RGB orthomosaic image corresponding to 24 June is shown in Figure 9d, clustering the study area into "good", "moderate", and "poor" classes. It is worth mentioning that 99.06% of the observed variance was explained by the first component (PC1) and only 0.94% of the variance retained by the second component (PC2) (Figure 9c). Moreover, there was a positive significant correlation between the red and green bands (Kendall's τ coefficient = 0.89 and p-value(Kendall) < 0.0001). The experimental field was completely covered by vegetation on 24 June enabling a proper PCA analysis over the vegetation pixels, without any soil background effects, to retrieve an accurate crop health map. Outside of the plot, asterisks *, **, ***, and **** reveal the statistical levels of significance and the acronym NS stands for non-significant. No spectral information is available for reference plot 1 on 20 March because of a flight planning constraint.

Crop Health
Both canopy cover and OSAVI were characterized by a decreasing trend on 24 June (Figures 4  and 7, respectively). The impact of century-old biochar on crop health during the senescence phase is presented in Figure 9. In addition, Figure 9 presents the two different crop health maps, obtained from the two different sensors, including the associated general methodological workflows as represented in Figure 9a,b for the multispectral sensor and Figure 9c,d for the RGB sensor. An unsupervised k-means clustered crop health map of OSAVI including "good", "moderate", and "poor" classes, corresponding to 24 June, is presented in Figure 9b. In addition, an unsupervised k-means clustered crop health map derived from the first component (PC1) from the principal component analysis of the red and green bands of the RGB orthomosaic image corresponding to 24 June is shown in Figure 9d, clustering the study area into "good", "moderate", and "poor" classes. It is worth mentioning that 99.06% of the observed variance was explained by the first component (PC1) and only 0.94% of the variance retained by the second component (PC2) (Figure 9c). Moreover, there was a positive significant correlation between the red and green bands (Kendall's τ coefficient = 0.89 and p-value (Kendall) < 0.0001). The experimental field was completely covered by vegetation on 24 June enabling a proper PCA analysis over the vegetation pixels, without any soil background effects, to retrieve an accurate crop health map. Categorical classification of the experimental plots into the crop health clusters is presented in Table 4 for both of the multispectral and RGB crop health maps (Figure 9). Based on the multispectral crop health map, the average health was 42.5% for the century-old biochar plots compared to 39.4% in the reference plots. Similarly, the average health based on the RGB crop health map was higher for the century-old biochar plots (54.6%) than the reference plots (51.6%). This finding shows that there is only a small positive impact of century-old biochar on crop health status during the senescence phase. The present study paves the way for further research focusing on the effects of biochar on crop health during the senescence phase. Categorical classification of the experimental plots into the crop health clusters is presented in Table 4 for both of the multispectral and RGB crop health maps (Figure 9). Based on the multispectral crop health map, the average health was 42.5% for the century-old biochar plots compared to 39.4% in the reference plots. Similarly, the average health based on the RGB crop health map was higher for the century-old biochar plots (54.6%) than the reference plots (51.6%). This finding shows that there is only a small positive impact of century-old biochar on crop health status during the senescence phase. The present study paves the way for further research focusing on the effects of biochar on crop health during the senescence phase. The cross-tabulation analysis of the computed crop health maps is presented in Table 5 comparing the number of all of the pixels within the study area classified into good, moderate, and poor health classes. The result yielded a clustering agreement of 74.82% ( Table 5). Out of a total of 564,760 pixels across the entire study area, 36.06% were classified as good, 31.54% as moderate, and 7.21% as poor by both methods (Table 5). Ground-based optical imagery, such as from the Fabry-Perot interferometer (FPI) camera system [39], is an essential step for future research focusing on the evaluation of crop health because of its improved spatial and spectral resolutions, enabling an accurate monitoring of the development of various plant diseases at the leaf level.

Ground-Based Crop Evaluation and Harvested Crop Yield
Ground-based measurements of crop traits are shown in Figure 10. Estimates of the general crop growth stages taken on 28 May 2019, scaled according to the BBCH, reported a more advanced crop stage for the century-old biochar plots (median BBCH = 43) compared to the reference plots (median BBCH of 40.5) as shown in Figure 10a. The observed difference in crop growth stages (BBCH) was statistically significant (p-value = 0.01265) between the reference and century-old biochar plots. This could be explained by the dark color of the century-old biochar plots, resulting in higher soil and plant temperatures stimulating plant development. The mean harvested crop yield (taken on 20 July 2019) value was 8.69 ± 0.69 t.ha −1 for the centuryold biochar plots compared to 8.84 ± 0.66 t.ha −1 in the reference plots; the observed contrast was not significant with a p-value of 0.42489 (Figure 10b). This finding contrasts with the positive promotion of the aboveground productivity for long-term biochar enriched soils [6,23]. This can be explained by the fact that our study area is situated in a quite rich soil type and hence the influence of century-old biochar on the harvested crop yield is limited. However, other studies considering poorer soil types have shown that the difference in yield due to the long-term biochar enrichment is more substantial [6].

UAV-Based Yield Map Generation
Relationships between harvested crop yields on 20 July (Section 3.2) and the UAV-based multispectral vegetation indices on the last acquisition date (i.e., 24 June which was the closest to the harvest date) are presented in Table 6. The s-CCCI, NDRE, NDVI, and CI-red exhibited a good relationship with the harvested crop yield as shown in Table 6. The obtained relationships were statistically significant with p-value(F-test) of 0.002, 0.002, 0.008, and 0.008 for s-CCCI, NDRE, NDVI, and CI-red, respectively. The aforementioned indices showed a strong correlation with the harvested crop yield with Kendall's τ coefficient of 0.60, 0.61, 0.54, and 0.54 for the s-CCCI, NDRE, NDVI, and CI-red, respectively ( Table 6). The obtained correlations were statistically significant with pvalue(Kendall) of 0.001, 0.0008, 0.004, and 0.004 for the s-CCCI, NDRE, NDVI, and CI-red, respectively (Table 6)   The mean harvested crop yield (taken on 20 July 2019) value was 8.69 ± 0.69 t.ha −1 for the century-old biochar plots compared to 8.84 ± 0.66 t.ha −1 in the reference plots; the observed contrast was not significant with a p-value of 0.42489 (Figure 10b). This finding contrasts with the positive promotion of the aboveground productivity for long-term biochar enriched soils [6,23]. This can be explained by the fact that our study area is situated in a quite rich soil type and hence the influence of century-old biochar on the harvested crop yield is limited. However, other studies considering poorer soil types have shown that the difference in yield due to the long-term biochar enrichment is more substantial [6].

UAV-Based Yield Map Generation
Relationships between harvested crop yields on 20 July (Section 3.2) and the UAV-based multispectral vegetation indices on the last acquisition date (i.e., 24 June which was the closest to the harvest date) are presented in Table 6. The s-CCCI, NDRE, NDVI, and CI-red exhibited a good relationship with the harvested crop yield as shown in Table 6. The obtained relationships were statistically significant with p-value (F-test) of 0.002, 0.002, 0.008, and 0.008 for s-CCCI, NDRE, NDVI, and CI-red, respectively. The aforementioned indices showed a strong correlation with the harvested crop yield with Kendall's τ coefficient of 0.60, 0.61, 0.54, and 0.54 for the s-CCCI, NDRE, NDVI, and CI-red, respectively ( Table 6). The obtained correlations were statistically significant with p-value (Kendall) of 0.001, 0.0008, 0.004, and 0.004 for the s-CCCI, NDRE, NDVI, and CI-red, respectively (Table 6) which is in line with the previous finding of Das and Singh (2012). Table 6. Results of correlation analysis of the multispectral vegetation indices, contained in Table 2, on 24 June and their relationships to the harvested crop yield on 20 July. The coefficient of determination, p-value of the statistical F-test, non-parametric Kendall's rank correlation tau, and its significance are expressed as R 2 , p-value (F-test) , τ, p-value (Kendall) , respectively. The predicted UAV-based multispectral crop yield map based on the identified linear regression fit between the s-CCCI and crop harvested yield, which reported the strongest relationship (Table 6), is shown in Figure 11a. The result of the predicted UAV-based MSP crop yield is shown in Figure 11b. The western part of the experimental field clearly shows lower yield values. However, the eastern part of the field exhibited rather high yield values (Figure 11b). The impact of century-old biochar on the predicted UAV-based MSP crop yield is shown in Figure 11c. The mean predicted MSP crop yield (mean value of the 11 plots) was not statistically different (p-value of 0.437) between the century-old biochar plots (8.82 ± 0.39 t.ha −1 ) and the reference plots (8.74 ± 0.42 t.ha −1 ) as shown in Figure 11c. This finding balances the previous works reporting a higher aboveground productivity as a consequence of biochar presence [6,23]. However, our results showed a highly significant difference (p < 0.0001) of the predicted MSP crop yield between the good and moderate classes derived from the MSP sensor (Figure 11d). Similarly, the observed difference of the predicted MSP crop yield was highly significant (p < 0.0001) between the moderate and poor classes based on the MSP sensor (Figure 9d). The mean predicted MSP crop yields were 8.96 ± 0.85, 8.50 ± 0.79, and 7.84 ± 0.83 t.ha −1 for the good, moderate, and poor classes derived from the MSP sensor (Figure 11d). In addition, the observed difference of the predicted MSP crop yield was highly significant (p < 0.0001) between the good and moderate classes, and between the moderate and poor classes, derived from the RGB sensor ( Figure 11e). The mean predicted MSP crop yields were 8.76 ± 0.91, 8.40 ± 0.86, and 7.54 ± 0.89 t.ha −1 for the good, moderate, and poor classes based on the RGB sensor ( Figure 11e). This result indicates that other factors (such as fertilization inputs, agricultural practices, and inherent possible pedological variations) are of higher importance than the presence of century-old biochar in determining the crop health and yield variability at the within-field scale.

Index
This paper provides a proof of concept of using high-resolution UAV images in assessing the impacts of century-old biochar on crop performance at the canopy scale. As such, the results of this study reveal important findings regarding the impact of century-old biochar on crop health and yield (in particular over the senescence period). However, since the approach described in this paper is only suitable for small-sized fields of a few hectares, further research should consider extending this analysis across a wider area, including a range of crops under contrasting weather conditions, using available fine resolution optical (e.g., Sentinel-2) and radar (e.g., Sentinel-1 Synthetic-aperture radar (SAR)) satellite images.