High spatio-temporal monitoring of century-old biochar effects on evapotranspiration through the ETLook model: a case study with UAV and satellite image fusion based on additive wavelet transform (AWT)

ABSTRACT It can be challenging to fuse remotely-sensed images with large differences in spatial resolutions. In this paper, we used additive wavelet transform (AWT) to fuse Landsat-8 (30 m) and unmanned aerial vehicle (UAV) images (7 cm and 3.7 cm for thermal and multispectral images, respectively) as one of the primary studies. AWT image fusion generated sharpened Landsat-8 (L-8) images which were significantly correlated with coarse resolution images, while also well preserving the spatial details. Surface albedo (α0), normalized difference vegetation index (NDVI), and surface temperature (ST) were computed from multispectral and thermal sensors on board of UAV and L-8 platforms. High-resolution UAV and AWT sharpened L-8 images were then used in ETLook model to estimate evapotranspiration (ET) across an agricultural farm enriched with century-old biochar. High spatio-temporal analysis demonstrated a significant decrease in α0 across the biochar patches during the early development stages of winter wheat. Moreover, biochar significantly stimulated the development of wheat canopies towards the middle of the cropping season. There were however no impacts at the end of the season due to dense wheat canopies covering the aggravated dark colour soil across the biochar patches. ST was not affected by biochar either at the beginning or towards the end of the season. Neither was there any impact of biochar on actual ET over the season. Our approach can help to develop robust techniques for fusion of UAV and satellite images in light of climate-smart agriculture, and is also applicable to other farms with any specific precision agricultural treatments. GRAPHICAL ABSTRACT


Introduction
Due to increases in global demands for water use, it has been critical to balance water supply and demand in agricultural settings. In this context, the use of hydrological models may help decisionmakers with water-management practices (Herman et al. 2018). The combination of water losses to atmosphere from soil surface (evaporation) and plants (transpiration) is called evapotranspiration (ET) (Hanson 1991). As such, knowledge of ET is essential to understand water balance and hydrological processes (Evett et al. 2012;Bastiaanssen et al. 1998a) since it massively transports lateral global energy through latent heat in the redistribution of earth's surface water (Mauser and Schädlich 1998). ET determination can provide important insights into plants interaction with available soil water content by means of water infiltration in the soil column or plant water uptake (Farquhar and Sharkey 1982).
In precision agricultural settings, biochar has been recently promoted as a means to improve crop water efficiency while also being resilient to climate variability (Fischer et al. 2019). Biochar is a solid charcoal product made by biomass pyrolysis (Sohi et al. 2010). Biochar incorporation to soil is an emerging approach in practices associated to climate-smart agriculture (Montanarella and Lugato 2013). Several previous studies have investigated the potential short-term role of biochar in enhancing soil water capacity (e.g. Farkas et al. 2020;Glaser, Lehmann, and Zech 2002), and hence, crop water use efficiency (Fischer et al. 2019). Short-term biochar amendment was found to regulate soil water capacity in some instances (e.g. Glaser, Lehmann, and Zech 2002;Abel et al. 2013;Agegnehu et al. 2016), though several studies reported divergent biochar impacts on soil water capacity (e.g. Gray et al. 2014;Schneider et al. 2020). Short-term biochar amendment can also positively affect soil properties such as bulk density and cation exchange capacity (Lychuk et al. 2015) as well as nutrients availability (Brantley et al. 2015). Such short-term effects of biochar on soils can in turn provoke crop growth (e.g. Biederman and Stanley Harpole 2013;Nair et al. 2017), and also, crop yield (e.g. Akhtar et al. 2014;Jeffery et al. 2017). Despite the aforementioned short-term biochar effects have been widely investigated, the specific short-term influence of biochar on ET across vegetable crops, due to biochar interaction with soil moisture and water retention capacity, has been studied only in few researches (e.g. Ahmed et al. 2019;EUROCHAR 2014). Furthermore, this is yet vague how long-term biochar enrichment (aged in agricultural soils more than a century ago) may affect soil-plant interactions. Few studies have shown that long-term biochar enrichment can lead to higher nutrient availability due to the increases of cation exchange capacity over time (Liang et al. 2006;Major et al. 2010). A recent study by Heidarian Dehkordi et al. (2020a) indicated a positive impact of long-term biochar enrichment on crop growth. The latter can be explained by soil physical properties, and mostly, higher available water content across the long-term biochar enriched soils (Kerré et al. 2017). However, to our knowledge, no comprehensive studies have yet evaluated long-term biochar effects on ET elements and the consequences on soil-plant system. Moreover, no research has yet been devoted to high spatio-temporal analysis of such century-old biochar effects on ET across agricultural farms at the canopy scale throughout an entire cropping season.
To fill the aforementioned gaps from the previous studies, we utilized high frequency remote sensing images offering high-resolution information on soilplant system which also remain nondestructive toward the monitoring period. Remote sensing is a rapidly developing technique for monitoring earth's surface. The unprecedented freely-available satellite images, such as Sentinel and Landsat, have been successfully utilized in many agricultural applications (e.g. Kooistra and Clevers 2016;Dutrieux et al. 2016;Clevers, Kooistra, and van Den Brande 2017;Chauhan et al. 2020). However, their coarse resolution images alongside the cloud-coverage remain problematic when addressing the precision agricultural needs for smallholder farmers (Peter et al. 2020). It is also worth drawing attention to the fact that century-old biochar patches usually compose an area of approximately 25 m diameter within the agricultural fields (Heidarian Dehkordi et al. 2020a). This underlines the necessity of deploying high spatial resolution remotely-sensed images to accurately monitor soil-plant interactions across the century-old biochar patches. High-resolution images acquired with unmanned aerial vehicles (UAVs) can achieve such necessary spatial resolution while also eliminating the drawback of clouds presence (Capolupo et al. 2015). Numerous studies have deployed UAVs for investigating a wide range of precision agricultural practices such as plant traits monitoring (e.g. Van Der Meij et al. 2017;Roosjen et al. 2018), disease management (e.g. Su et al. 2018;Heidarian Dehkordi et al. 2020c), and biomass productivity (e.g. Ten Harkel, Bartholomeus, and Kooistra 2020; Ballesteros et al. 2018). Conducting high temporal frequency of UAV imageries is however time and effort demanding in terms of image acquisition and processing.
Owing to advances in remote sensing, image fusion technique can combine images of different sensors to benefit from the associated advantages of each dataset (Zhou et al. 2021). For example, Weng, Fu, and Gao (2014) endeavored to fuse images of moderate resolution imaging spectroradiometer (MODIS) having daily temporal resolution and Landsat with finer spatial resolution to generate daily fine spatial resolution surface temperature (ST) images across Los Angeles County, California, United States. Their fusion was performed in spatio-temporal adaptive data fusion algorithm for temperature mapping (SADFAT). The authors indicated that SADFAT fusion algorithm contains a few limitations in predicting ST changes. Another study by Gao et al. (2006) fused high temporal resolution MODIS images with high spatial resolution Landsat images using spatial and temporal adaptive reflectance fusion model (STARFM). They concluded a good efficiency of STARFM in predicting daily surface reflectance. Several authors have explored various image fusion algorithms, developed based on STARFM, to successfully generate high spatio-temporal resolution vegetation images allowing for an improved precision agricultural monitoring (e.g. Rao et al. 2015;Liao et al. 2017;Chen et al. 2018). Despite space borne image fusion algorithms have been adopted to different satellites and for various intended applications and spatial scales, it remains challenging to fuse images with large differences in their spatial resolutions (Chen, Wang, and Liang 2019). Moreover, the main commonality among remote sensing image fusion techniques is that the input images are acquired by satellites with coarse spatial resolution in comparison to UAV images. To date, only few studies have investigated image fusion algorithms for fusing satellite and UAV images composing extremely different spatial resolutions from meters in the former to centimeters in the latter. For example, a recent study by Zhao et al. (2019) fused UAV and Sentinel-2A images using Gram-Schmidt (GS) transformation for crop classification. Another research by Sagan et al. (2019b) used temporal fusion of virtual constellation of UAV and WorldView-3 for crop phenotyping and stress monitoring. Gevaert et al. (2017) compared the accuracy of unmixing-based algorithm versus STARFM while fusing UAV and Formasat-2 images. They highlighted a better accuracy of unmixingbased method as compared to STARFM while retrieving vegetation indices by fusing UAV and satellite images. One important factor when fusing the remote sensing images is the preservation of both radiometric and geometric information in the fused image (Bama et al. 2013). Additive wavelet transform (AWT) is a robust image fusion method which is typically developed based on intensity-hue-saturation (IHS) transform algorithm (Nünez et al. 1999). AWT is capable of preserving the detailed information of both input images while enhancing the spectral and spatial quality of the fused image (Johnson 2014). Bama et al. (2013) reported a strong ability of AWT in preserving both spectral and geometric information while fusing multispectral (2.5 m) and panchromatic (46 cm) bands of WorldView-2 images. Their finding paves the way for testing AWT while fusing UAV and satellite images composing large differences in their spatial resolutions.
One way to estimate the ET across the agricultural settings is assimilating remote sensing images to energy balance algorithms (Mokhtari et al. 2021). The energy balance algorithms are classified into one-source or two source modeling frameworks. Two-source approaches can differentiate turbulent heat fluxes between soil and canopy, while onesource approaches use only a single resistance (Khan, Baik, and Choi 2021). Surface energy balance algorithm for land (SEBAL; Bastiaanssen et al. (1998b)), surface energy budget system (SEBS; Su et al. (2002)), and mapping ET at high resolution with internalized calibration (METRIC; Allen, Tasumi, and Trezza (2007)) are amongst the most popular one-source models to estimate ET. Two-source energy balance (TSET; Norman, Kustas, and Humes (1995)), two-source time integrated model (TSTIM; Anderson et al. (1997)), and ETLook ) are well-known two-source energy balance models. Two-source energy balance models are capable of providing more accurate estimates of energy budgets by partitioning ET into evaporation and transpiration separately (Norman, Kustas, and Humes 1995). The latter is the key advantage of two-source energy balance models since the radiometric response of vegetation pixels differs from that of soil pixels (Khan, Baik, and Choi 2021).
ETLook model makes predictions of actual ET based on the Penman-Monteith (P-M) approach (Monteith 1965) through the combination of energy balance and aerodynamic equations (as described in Allen et al. 1998). The main difference between ETLook and many other surface energy balance models is the differentiation of net available radiation and resistance in accordance with the proportion of vegetation . As such, the model eliminates the need for manual cold/hot pixel selection in the SEBAL model and solves the evaporation and transpiration components of the P-M equation (Blatchford et al. 2020).
As such, in this paper, we have attempted to fuse UAV and Landsat images as one of the first studies according to the best of our knowledge, resulting in sharpened Landsat images with the high spatial resolution of the former that is needed for high spatiotemporal monitoring. The fused images, in combination with meteorological data, were then used in ETLook model to investigate century-old biochar effects on ET across an agricultural farm planted with winter wheat. The presented approach may help to develop in the future more robust image fusion techniques for fusing UAV and satellite images in light of climate-smart agriculture. This approach is also applicable to other study sites with any specific precision agricultural treatments.

Study site
Our study site (50°30ʹN -50°31ʹN; 4°44ʹE -4°45ʹE) is a 13 ha (nearly 32 acre) agricultural farm in Gembloux district, province of Namur, Belgium (Figure 1). The Gembloux district is famous for its historical biochar patches aged in agricultural soils more than 150 years ago. The selected farm was predominantly covered with oak, hornbeam, beech, and hazel forests (Hardy et al. 2017) and had been turned into cropland since the eighteenth century resulting in biochar enrichment in several spots within the site. The field is composed of Luvisol with silt loam texture, which was planted with winter wheat (Triticum aestivum L.). The district possesses a temperate climate with accumulated monthly precipitation of 173 mm, 41 mm, 155 mm, and 94 mm from March to June 2019, respectively, recorded by an on-site automated weather station. In addition, the average monthly temperature was 8.2°C, 10.2°C, 11.2°C, and 17.9°C from March to June 2019, respectively. Eleven 10 by 10 m plots, labeled from 1 to 11, were marked inside the century-old biochar plots and their corresponding reference plots as shown by red and blue squares in Figure 1, respectively. A distance of nearly 45 m was considered between the reference and biochar plots to ensure that no biochar is present in the reference plots. It is worth stating that the two unused biochar plots in Figure 1 were instrumented to monitor nutrient cycling and water dynamics through joint experiments, and hence, were not evaluated in the present study. Agricultural inputs were broadcast uniformly over the field; the northern portion of the field was rather flat, although there was a very mild slope descending into the southern part of the field with no topographic variations (Hedarian Dehkordi et al. 2020b).

ETLook theoretical framework
Surface temperature (ST), surface albedo (α 0 ), and normalized difference vegetation index (NDVI) are the key remote sensing parameters (as contained in Figure 3) in the ETLook model . It is also important to highlight the benefit of remotely-sensed ST in the ETLook model to estimate relative root zone soil moisture (Yang et al. 2015). ETLook actual ET rates are estimated in mm/day. It is worth noting that the associated units are converted from W.m −2 using the latent heat of evaporation while being dependent on the air temperature. The primary version of ETLook is capable of modeling ET across large-scale landscapes using either coarse (≥ 1 km) or moderate (100 to 1000 m) resolution satellite products such as advanced microwave scanning radiometer-earth observing system (AMSR-E) and MODIS, respectively . The modified version of ETLook was recently well-validated by using fine resolution (10 to 100 m) satellite images (i.e. Landsat 5, 7, and 8) over Africa (Blatchford et al. 2020).
However, only few researches have attempted to evaluate such surface energy balance models with high-resolution remote sensing data. There exists a gap to examine the role of high-resolution surface energy balance modeling approaches in addressing precision agricultural applications. Lastly, a research by Heidarian Dehkordi (2017) indicated that it is potentially feasible to model high spatial distribution of ET across small-scale agricultural fields by assimilating high-resolution aerial remotely-sensed images (on the order of sub-centimeters collected by an UAV) into the ETLook model. Therefore, we tested the ETLook model on high-resolution remotely-sensed images to examine if this may bring out high spatiotemporal information on actual ET of winter wheat over century-old biochar and reference patches.  DJI Matrice 100 is a lightweight (2.3 kg) vertical take-off and landing quadcopter with 40 min hovering time using two TB48D batteries. A FLIR Vue Pro R 640 thermal camera (FLIR Systems, Wilsonville, Oregon, USA) was used to collect the thermal dataset ( Figure 2). The thermal sensor has a noise equivalent temperature difference (NETD) thermal sensitivity of 0.05°C and a thermal accuracy of ± 0.05°C (Table 1). Global positioning system (GPS) coordinates of each thermal image were recorded during the flights using a FLIR GPS data logger (a.k.a. geo-tagger) installed on the UAV platform ( Figure 2). Multispectral data were collected by a MicaSense RedEdge-M camera (MicaSense, Seattle, WA, USA) in blue, green, red, red-edge, and near-infrared spectral channels ( Table 2). The multispectral camera was mounted on the UAV platform through a gimbal damper installing at a fixed viewing zenith angle of 46°.  . Schematic workflow of data collection, various data processing steps, image fusion, and crop-water status modeling proposed in the present study. UAV, OLI, TIRS, and NDVI refer to unmanned aerial vehicle, operational land imager, thermal infrared sensor, and normalized difference vegetation index, respectively. A downwelling light sensor (DLS) was also installed on the UAV platform to adjust within-flight variable light conditions. We operated the UAV flights around noon under clear-sky conditions allowing for optimum thermal and multispectral acquisitions. A general workflow diagram of data acquisition, data processing, image fusion, and crop-water status modeling is demonstrated in Figure 3. UAV flight parameters were programmed in Pix4D capture (Pix4D S.A., Lausanne, Switzerland). The flying altitude of the UAV platform was 60 m above ground level (AGL), the image overlap was 75% in both longitudinal and transversal directions, and the flying speed was 7 ms −1 . Flying with the aforementioned flight characteristics yielded high spatial resolution images with ground sampling distance (GSD) of 7 cm and 3.7 cm for FLIR and MicaSense RedEdge-M cameras, respectively. UAV image alignments were conducted in Pix4Dmapper (Pix4D S.A., Lausanne, Switzerland) for both thermal and multispectral datasets. Five ground control points (GCPs), i.e. four in the corner and one in the middle of the experimental field (Figure 1), were used to improve the geometric accuracy of the thermal and multispectral orthomosaic images (Figure 3).

UAV image acquisition and image processing
A reference reflectance target (MicaSense calibration panel, Seattle, WA, USA) was sensed with MicaSense RedEdge-M camera, over a flat ground with no shadows, before and after the UAV flights to perform the radiometric calibration of the multispectral images (Figure 3) in Pix4Dmapper. Moreover, the information on in-flight sensor's sensitivity to light, acquired by the DLS during each UAV flight, was used to optimize the radiometric calibration.
A hand-held FLIR TG167 thermal imager (FLIR Systems, Wilsonville, Oregon, USA) was used to collect reference ground-based thermal data ( Figure 3) during the UAV flights. We used twelve reference targets of 2.5 × 2.5 m with three different colors according to Santesteban et al. (2017), and as such, four black, four blue, and four white panels were deployed within the experimental site. ST of each target was measured with the hand-held thermal imager at five replicates with 10 s delay, taken both at the beginning and the end of each UAV flight. The average of the ten measurements made over each target was used to describe the ground temperature of each target. Subsequently, sensor temperature of each target, measured by the FLIR Vue Pro R camera on board of the UAV, was extracted from the thermal orthomosaic image in Quantum GIS software (QGIS, Open Source Geospatial Foundation, Chicago, IL, USA). For each flight, a linear regression line was modeled between the surface temperature and sensor temperature of the twelve reference targets similar to the methodology proposed by Sagan et al. (2019a). Finally, we used the attained regression fit of each flight in order to convert UAV sensor temperature orthomosaic image into UAV surface temperature orthomosaic image ( Figure 3).

UAV-based inputs for ETLook
As stated before in Section 2.3.1, ST values were calculated from FLIR Vue Pro R camera, mounted on the UAV, based on the ground-based temperature data collected with hand-held FLIR TG167 thermal imager.
We used the approach of Liang et al. (2003), which has been initially designed for Landsat-7 enhanced thematic mapper plus (ETM+), by adapting the surface reflectance spectra to the multispectral MicaSense RedEdge-M sensor (as proposed by Heidarian Dehkordi (2017)) in order to retrieve the surface albedo (α 0 ) as follows: where R blue , R red , and R NIR are the reflectance values at the subscripted spectral bands of the MicaSense RedEdge-M sensor (Table 2). In addition, NDVI was calculated by taking the ratio between the difference and sum of the spectral responses in the red and nearinfrared spectral channels (Rouse et al. 1974) of MicaSense RedEdge-M contained in Table 2, as below:

Landsat-8 data
Cloud-free Landsat-8 (L-8) images (Appendix A) of the study site during the 2019 cropping season were available on 24 April, 10 May, and 27 June. L-8 operational land imager (OLI) atmospherically corrected surface reflectance images of path 198 and row 25 at worldwide reference system (WRS) were downloaded from the USGS earth resources observation and science platform (https://landsat. usgs.gov/landsat-surface-reflectance-highlevel-data -products) as on-demand products. OLI reflectance images were co-registered to MicaSense RedEdge-M reflectance images ( Figure 3) using the image-to -image registration tool in QGIS in order to optimize the geometric accuracy needed for the fusion of OLI and MicaSense RedEdge-M reflectance images (Section 2.4.2). The latter avoids pixel shifting and geometric distortions while conducting the image fusion (Zhou et al. 2021). We used Eq.
where R blue , R red , R NIR , R SWIR1 , and R SWIR2 are the reflectance values at the subscripted spectral bands of Lansat-8 OLI sensor. L-8 thermal infrared sensor (TIRS) sensor captures the earth ST using two thermal infrared channels (Appendix A; Table A1). The corresponding L-8 TIRS level-1 thermal products (WRS path 198 and row 25) were collected from the USGS archives (https://earthexplorer.usgs.gov). In the present study, we only investigated TIRS band 10 since its central wavelength, 10,895 nm (Appendix A; Table A1), is closer to FLIR Vue Pro R central wavelength, 10,250 nm (Table 1), as compared to TIRS band 11.
We retrieved surface temperature (ST) from TIRS band 10 ( Figure 3) following the temperature emissivity separation (TES) algorithm (Zareie, Khosravi, and Nasiri 2016). Similar to the multispectral dataset, we first co-registered TIRS images to FLIR Vue Pro R thermal images ( Figure 3) in QGIS through the image-to-image registration function for the optimization of the geometric accuracy (Zhou et al. 2021) required for the fusion of TIRS and FLIR Vue Pro R thermal images (Section 2.4.2). Calibrated level-1 pixel values, a.k.a. digital numbers (DN cal ), recorded by TIRS band 10 were then converted to top-of-atmosphere (TOA) radiance (L λ ) values based on radiance multiplicative (M L ) and additive (A L ) scaling factors, stored in the metadata file of L-8, as follows: Next, TOA radiance values (L λ ) were translated to TOA sensor brightness temperature (T TOA ) values ( Figure 3) considering thermal conversion constants (K1 and K2) which are reported in the metadata file of L-8 as: To estimate ST from T TOA , the atmospheric effects should be corrected (Vanhellemont 2020). We utilized the theory of Zareie, Khosravi, and Nasiri (2016) to determine surface emissivity (SE) from land use type, surface roughness, and vegetation coverage. According to the relationships between emissivity and NDVI of natural surfaces reported in Van De Griend and Owe (1991), surface emissivity values of 0.98 and 0.91 for vegetation and soil pixels, respectively were considered in the present study. Plausible discrepancies in surface roughness were omitted from the computation of surface emissivity since geometric distribution effects and internal reflections are relatively small within the homogenous winter wheat canopies throughout the experimental site that is without any topographic variation (Section 2.1). Hence, we can write the surface emissivity as a linear combination of NDVI with soil and vegetation emissivity as: where SE v , and SE s are the surface emissivity of vegetation and soil pixels, respectively. NDVI max and NDVI min were extracted from L-8 OLI NDVI for fully vegetated and contaminated bare soil pixels, respectively. It is worth noting that the coarse spatial resolution of OLI spectral bands used to compute NDVI (30 m as contained in Appendix A; Table A1) does not allow the selection of pure bare soil pixels for determination of NDVI min in Eq. (6). Finally, we computed ST from T TOA and SE so that: where λ is the wavelength of the emitted radiance in TIRS band 10 (Appendix A; Table A1), σ and h are Boltzmann and Planck constants (Zareie, Khosravi, and Nasiri 2016), respectively, and c is the velocity of light (2.99 × 10 8 ms −1 ). Eq. (7) allows us to retrieve ST in Kelvin, we hence translated ST values to Celsius to be compatible with those of UAV FLIR Vue Pro R.

Image fusion
With high-resolution UAV remote sensing it is common to postulate a modeling resolution in the order of a few centimeters (Gevaert et al. 2015). However, surface energy balance models can never achieve such subcentimeter resolution using satellite constellation. Image fusion is a possible path to enhance the modeling resolution while using multi-sensor satellite remote sensing images (Guzinski et al. 2020). Unfortunately, very little information is yet available for fusing Satellite and UAV imageries (Maimaitijiang et al. 2020). To fill this resolution gap, we hence convolved all co-registered L-8 OLI and TIRS derivatives to those of UAV MicaSense RedEdge-M and FLIR Vue Pro R (Figure 3), respectively based on AWT. To perform AWT, we first used a B 3 cubic spline filter (Nünez et al. 1999) in order to add the wavelet coefficients of the high-resolution image (collected by the UAV sensors) to the co-registered coarse resolution image acquired by L-8 sensors. Then we computed the sharpened image similar to fast IHS (FIHS)-like image fusion algorithm (Tu et al. 2001) as follows: where Landsat-8 high is the sharpened image of L-8 as the output of AWT fusion, Landsat-8 low is the original co-registered coarse resolution L-8 image which has been resampled to the spatial resolution of the original UAV image without any spectral enhancement, UAV is the original high-resolution image captured by the UAV sensor, and UAV low is the degraded UAV image using the B 3 cubic spline function. It is worth mentioning that AWT was applied separately to each remotelysensed input needed for the ETLook model (

ETLook meteorological inputs
Meteorological input datasets comprised of air temperature, wind speed, solar radiation, and relative humidity are needed to describe the meteorological condition in the ETLook model ) as shown in Figure 3. To apply the ETLook model on the study site, meteorological datasets should be obtained at both sensor overpass time (a.k. a. instant inputs) and daily averages (periodic inputs). As mentioned earlier (Section 2.1), meteorological data were recorded as 30 min averages by an on-site automated weather station. This dataset is presented in Table 3 for our nine acquisition dates over the year 2019. The closest measurement of the weather station to sensors (i.e. UAV or satellite) overpass time was used in the ETLook model for each acquisition. It is worth mentioning that UAV overpass time was considered as the middle of the start-end UAV flight times for each acquisition date. Satellite overpass time included in the metadata file of L-8 was used for each monitoring date.

Spatiotemporal analysis
Data analysis was carried out using R (R Core Team, Foundation for Statistical Computing, Vienna, Austria) through its fundamental packages and dependencies libraries. Moreover, we utilized violin plot visualization to express the density of the intended ET elements. We used UAV and (AWTsharpened) satellite derived maps of ST, α 0 , NDVI, and actual ET to assess whether there may be an impact of century-old biochar enrichment on cropwater status indicators. Statistical t-test (Kim 2015) was undertaken to determine the significance of the differences across the study site when comparing biochar and reference plots over the 2019 cropping season.
In addition, correlation analysis was undertaken to test the correlations between the coarse resolution L-8 images (used as inputs to the AWT image fusion with the UAV images) and the corresponding AWTsharpened L-8 images. For this spatial correlation, each coarse resolution L-8 pixel (GSD of 30 m) was compared with the average of AWT-sharpened L-8 pixel values (GSD of 7 cm and 3.7 cm for thermal and multispectral images, respectively) which are covering that coarse resolution L-8 pixel. Given the different range of datasets in the present study (i.e. α 0 , ST, and NDVI), we computed the normalized root mean square error (nRMSE), normalized by range (Finger et al. 2021), using Eq. (9), allowing us to evaluate AWT image fusion performance. Furthermore, we examined the plausible AWT persistent errors (i.e. over and under predictions) through percent bias (PBIAS) in Eq. (10) as suggested by Abolafia-Rosenzweig et al. (2021).
nRMSEð%Þ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi where i identifies each coarse resolution L-8 pixel, n represents the total number of pixels, X i is the coarse resolution L-8 model input value, and Y i is the corresponding AWT sharpened L-8 model output value.

Image fusion performance assessment
AWT image fusion yielded sharpened L-8, images which were all linearly correlated (p-value < 0.001) with their corresponding coarse resolution images for each fusion date (Figure 4, Table 4). As it stands, the strongest correlations were obtained in NDVI image fusion (Figure 4g-i, Table 4). For ST image fusion, a small decrease of the correlation coefficients was observed, although r values were yet all greater than 0.6 (Figure 4d-f, Table 4). When α 0 images were used for AWT fusion, correlation coefficients were found to be the weakest (Figure 4a-c, Table 4). The latter could be caused by shortwave infrared 1 and 2 bands (alongside their associated weights) that are to be included in the determination of L-8 surface albedo (Eq. (3)) as compared to UAV surface albedo (Eq. (1)). Liang et al. (2003) particularly designed these weights based on the correlation between TOA reflectance and ground-based albedo measurements while considering only the blue, red, NIR, SWIR1, and SWIR2 spectral bands of Landsat ETM+. This approach allows a more accurate retrieval of α 0 in comparison to unweighted average approach, as for instance by Smit (2010), that only considers an average of reflectance values over the spectrum to compute the α 0 .
It should be stated that such overall good agreements in AWT results, when fusing UAV and L-8 images, can be partially attributed to the importance of optimized image co-registration as was initially raised by Zhou et al. (2021). AWT was also pointed out by Chen, Wang, and Liang (2019) as a robust image fusion technique for fusing satellite images with large differences in spatial resolutions such as L-8 (30 m) and Gaofen-2 (4 m) or MODIS (250 m) and SPOT (10 m) images. It is worth noting that the homogeneity of the pixel's content between UAV and L-8 images is also a major factor in obtaining such good agreements in AWT results. Future research should examine AWT image fusion over the heterogeneity conditions. Moreover, the proposed approach should be further evaluated using high-resolution satellites such as Landsat/Sentinel together with coarse resolution images Like MODIS products. Despite the generated α 0 and NDVI values by AWT were flattened around the 1:1 line on 24 April, most of the values were overestimated on 10 May, while many others were underestimated on 27 June (Figure 4). On the contrary, AWT generated ST values coherently positioned along the 1:1 line for the three image fusion dates (see Figure 4d-f). This indicates that AWT thermal image fusion was not adversely affected by the radiometric inconsistency between FLIR Vue Pro R (central wavelength of 10,250 nm; Table 1) and L-8 TIRS (central wavelength of 10,895 nm; Appendix A; Table A1) thermal sensors. Our image fusion analysis suggests that AWT may play a key role in the development of well-known spatiotemporal image fusion confidence intervals) are presented in red and dashed-black lines, respectively. AWT and NDVI refer to additive wavelet transform and normalized difference vegetation index, respectively. N is the total number of Landsat-8 pixels throughout the study site used for image fusion with UAV images; 8 pixels are missing for surface temperature image fusion on 27 June (f) due to malfunction of the thermal sensor during a part of the UAV flight. Asterisk * indicates the significance level of p-value less than 0.001. The statistical results are presented in Table 4   Table 4. Statistics comparing the coarse resolution and sharpened Landsat-8 images based on AWT fusion with their corresponding unmanned aerial vehicle (UAV) images for surface albedo (α0), surface temperature (ST), and NDVI on 24 April, 10 May, and 27 June 2019 as contained in Figure 4. The correlation coefficient (r), normalized root mean square error (nRMSE), and percent bias (PBIAS) are presented. Asterisk * indicates the significance level of p-value less than 0.001.

Index
Surface algorithms in space-borne remote sensing domain, such as regression model fitting, spatial filtering, and residual compensation (Fit-FC), which are yet highly sensitive to the radiometric inconsistency between the input images (Zhou et al. 2021). In addition, less attained PBIAS values in AWT surface temperature fusion as compared to multispectral image fusions (α 0 and NDVI) is also expected since the ratio of spatial resolution of the input images decreases (Chen, Wang, and Liang 2019) from 30 m/3.7 cm in multispectral images to 30 m/7 cm for thermal images. Zhou et al. (2021) also indicated that several image fusion algorithms, such as STARFM or one-pair dictionary-learning (OPDL), can be heavily affected by increases in spatial resolution ratio of the input images. Future work may also consider including radar data (e.g. Sentinel-1) to enhance the prediction of ST by taking the advantage of soil water variability in the thermal images as raised by Amazirh, Merlin, and Er-Raki (2019). Our primary objective of implementing AWT image fusion was to generate high spatiotemporal resolution maps throughout the cropping season, with an overarching goal of comparing biochar and reference plots across the study site. It is worth mentioning that the acquisition times were not identical between L-8 and UAV images (Table 3). Although ST values can be different due to various acquisition times, this has, however, no impact on our spatiotemporal analysis because ST discrepancies within the field (i.e. the relative ST differences between biochar and reference plots) are essential to provide an indication of biochar effects on ST. Moreover, to further assess the strength of the implemented image fusion approach, the recorded air temperature values by the on-site weather station were plotted against the average remotely sensed ST (averaged across the field) values retrieved from either L-8 or UAV images for all the nine acquisition dates (Appendix B; Figure B1). Remotely sensed ST values exhibited a strong correlation (r = 0.95, p-value < 0.001; Appendix B; Figure B1) with the air temperature values. This result is confirming the validity of implementing an AWT image fusion approach using input images with different acquisition times.
It is noteworthy that conclusions regarding overall AWT performance for UAV and satellite image fusion should be drawn with care because only three fusion dates were investigated in the present study. However, these promising results may help to develop in the future more robust image fusion techniques for the merging of UAV and satellite images. We recommend repeating the implemented methodology over the other dates and datasets to draw more strong conclusions regarding AWT fusion of UAV and satellite images. Moreover, future investigation should be aimed at improving IHS-like image fusion techniques with less over and under predictions while merging UAV and satellite images.
Figure 5a-f shows L-8 and UAV images, acquired on 10 May and 13 May, respectively, that were used as inputs to the AWT image fusion in order to generate sharpened L-8 images for 10 May (Figure 5g-i). The corresponding graphs for the other fusion dates are depicted in Appendix C; Figure C1 and Figure C2. α 0 exhibits similar patterns for L-8 ( Figure 5a) and UAV (Figure 5d) images in which the eastern part of the field revealed less values for α 0 while slightly higher values in the middle portion of the site. Neither was there any footprint of biochar patches (see Figure 1 for their locations) in the α 0 images acquired by L-8 ( Figure 5a) or UAV (Figure 5d) in May. It should be that the high proportion of vegetation, as is clearly seen from the NDVI images (Figure 5c,f) with pixel values mostly greater than 0.5, had entirely covered the background dark soil pixels within the biochar patches. Within field α 0 , variations observed in L-8 image ( Figure 5a) were yet preserved in the AWTsharpened image (Figure 5g). In addition, AWTsharpened L-8 α 0 image (Figure 5g) matched well with the high-resolution UAV α 0 image (Figure 5d) in light of the spatial discrepancies. More precisely, both images reveal an area with relatively higher α 0 in the southern part of the site while a zone with less values of α 0 in the north-eastern part of the site. Furthermore, microscale spatial surface details (such as tractor driving paths as well as the roads along the field borders with relatively less α 0 ) were clearly discerned in both images. It is also encouraging that the AWT sharpened α 0 image (Figure 5g) contained most of the spatial surface details (e.g. tractor driving paths), which were clearly discerned in the highresolution UAV image (Figure 5d) but not in the coarse resolution L-8 image (Figure 5a). Chen, Wang, and Liang (2019) also reported that spatial details can be well preserved while fusing satellite images with extremely different spatial resolutions through AWT.
As expected, L-8 and UAV images show similar overall pattern of NDVI across the site (Figure 5c,f). High NDVI values accounted for the majority of the vegetation pixels throughout the field, while there were some soil pixels in the eastern and northern borders of the site that appeared darker in the NDVI images. What was noteworthy is the impact of contaminated (i.e. soil and vegetation) pixels in the coarse resolution L-8 image, with a GSD of 30 m, declining the upper range of NDVI values from 1 to a NDVI of 0.7 (Figure 5c). However, high-resolution UAV derived NDVI image (GSD of 3.7 cm) was able to capture pure vegetation pixels (NDVI > 0.80) over the site (Figure 5f). In general, AWT-sharpened L-8 NDVI image (Figure 5i) revealed similar patterns to both L-8 ( Figure 5c) and UAV (Figure 5f) NDVI images. Interestingly, AWT was again capable of preserving microscale spatial surface circumstances while generating sharpened L-8 NDVI image (Figure 5i) that was the same for the aforementioned α 0 image fusion.
In contrast to α 0 as well as NDVI image fusion chains, ST images were not entirely identical between L-8 ( Figure 5b) and UAV (Figure 5e) imageries. ST values were found to be of various ranges in L-8 image of 10 May (Figure 5b) as compared to the UAV image on 13 May (Figure 5e), which is likely ascribable to differences in sensor overflights, i.e. 10:33 for L-8 versus 13:44 for UAV (Table 3). UAV overpass was the most ideal for revealing thermal discrepancies across our agricultural farm since it was operated when the surface had been entirely warmed up till noon. As shown in Figure 5e, ST values were ranging from 16°C to 25°C in the UAV image, while L-8 indicated ST with pixel values ranging solely from 14.5°C to 17.5°C (Figure 5b). This small range of temperature values in L-8 image might be also explained by the contaminated pixels increasing ST of the pure vegetation pixels and vice versa. More specifically, the minimum ST values of L-8 and UAV differed in only 1.5°C, whereas there was a difference of 7.5°C in the maximum ST values. This indicates the significance not only of an adequate time for conducting the thermal imagery, but also of sensor spatial resolution recording the ST. The latter was also reported by Weng, Fu, and Gao (2014) Table). This is in agreement with Chen, Wang, and Liang (2019) who also documented that AWT serves a small sensitivity to the differences in acquisition time of the input images.
Finally, the actual ET maps (modeled based on ETLook) are shown in Figure 5j,k for pure UAV case of 13 May along with the fused L-8 and UAV scenario on 10 May, respectively. The mid-northern part of the field clearly revealed low daily ET rates (less than 0.1 mm) on both 10 and 13 May. The western portion of the site exhibited a gradual increase in actual ET rates from 10 May (Figure 5j) to 13 May (Figure 5k), while the middle part of the field underwent a sharp decline in actual ET values over this time. Similar patterns of actual ET were observed on 10 and 13 May for the eastern portion of the site, indicating high (greater than 2 mm/day) and low (less than 0.4 mm/day) ET rates in its left and right sides, respectively. More precisely, the average actual ET across the site increased from 1.32 ± 0.79 mm/day on 10 May to 1.97 ± 0.97 mm/day on 13 May, which is highly attributed to 6.25 mm rainfall of 11 May as recorded by the on-site weather station. Nevertheless, gradual increases of solar radiation from 10 to 13 May (Table 3) made also a prominent contribution to the actual ET variety between these two days. The latter was also observed by Xu and Yu (2020), who identified solar radiation as one of the primary causative parameters in ET estimates. Last and not least was to provide an indication of the modeling inputs on the ET estimates through the ETLook model. This is apparent from Figure 5 that the spatial patterns of ST images were significantly preserved in the estimated actual ET maps. By applying correlation analysis to remotely sensed ETLook inputs, we verified that ST was, indeed, significantly correlated with the actual ET estimates (r > 0.88, p-value < 0.001; graphs not shown). This reveals that majority of surface energy balance models are yet highly dependent on the remotely sensed estimates of ST (e.g. Sun et al. 2016;Filgueiras et al. 2019;Ait Hssaine et al. 2020), although several models such as simplified surface energy balance index (S-SEBI) seem to be less sensitive to ST (Da Rocha et al. 2020). In contrast, the actual ET estimates showed no substantial correlation with either α 0 (r < 0.35, p-value > 0.05; graphs not shown) or NDVI (r < 0.31, p-value > 0.05; graphs not shown) values. This confirms the visual inspection of the overall patterns presented in Figure 5, indicating that spatial patterns of neither α 0 nor NDVI were remarkably preserved in the estimated actual ET maps. Our finding demonstrates that most of the ET variations can be derived from the observed discrepancies in ST. The latter is expected because of the well-developed wheat canopies in which higher ST values exhibit greater soil moisture stress and, hence, lower ET rates. Therefore, investigating the sensitivity of ETLook model to its various parameters should be a research priority.

Surface albedo
Century-old biochar enrichment in patches across the site (see Figure 1 for locations) caused notable decreases in surface albedo for most of the plots (Figure 6). The dark background soils across the biochar plots led to less α 0 values (Verheijen et al. 2013) compared to those of the reference plots. The observed statistical significance of differences is presented in Appendix D for each plot individually. All the 11 plots indicated significantly less α 0 across biochar plots from 20 March till 16 April (Appendix D; Table D1), except for plot 10 on 20 March when the observed difference was not significant. On 24 April, there has yet been a significant decline in α 0 for the majority of the biochar plots, except in plots 1, 8, 9, and 10. During the second half of the cropping season, lasting from the end of April until late June, there were no significant differences in α 0 between biochar and reference plots (Appendix D; Table C1). The very dense coverage of wheat canopies over this period, which had entirely covered the background dark soils across the biochar patches, seemed to be the main cause of this. Furthermore, the plausible local shadows due to this very dense wheat canopy structures had somewhat covered the dark background soil across the biochar plots (Zhang et al. 2017). α 0 averaged across the 11 biochar plots within the site was significantly different from the reference plots until 16 April (Appendix D; Table D1). Interestingly, the observed significance of differences diminished from March to April with p-values being less than 0.001 and 0.05, respectively. This was accompanied by less average surface albedo values across the biochar plots; the average α 0 across the biochar plots (in comparison to the reference ones) was 0.09 ± 0.01 (0.12 ± 0.00) on 20 March, 0.11 ± 0.01 (0.14 ± 0.00) on 28 March, and 0.17 ± 0.02 (0.19 ± 0.01) on 16 April, respectively. From 16 April, onwards, the average α 0 did not appear to be significantly affected by centuryold biochar enrichment due to dense wheat canopies. A similar conclusion was made by Genesio et al. (2012) who reported a notable decrease in α 0 due to biochar enrichment in soils of Tuscany, Italy, planted with winter durum wheat. They also showed that this impact seem to be less effective by the development of durum wheat plants toward the growing season. Zhang et al. (2017) also indicated a significant decrease in α 0 of biochar-amended soils over the early development stages of maize and wheat crops.
Century-old biochar and reference plots indicated rather similar distributions throughout the monitoring period (Figure 6). At the beginning of the growing season, the distribution was quite uniform in both biochar and reference plots, with values mostly from minimum to third quartile composing the largest α 0 density for the majority of the plots. This seems to be associated with the surface reflectance coming from a large number of soil pixels within the sparse wheat canopies over this period. In addition, a small number of emerged wheat plants within the plots exhibited α 0 values above the identified maximum. Subsequently, the distribution shapes were becoming much more uniform with evenly distributed α 0 values for most of the plots in April (Figure 6). This uniform distribution  Table D1.
is most likely attributable to the reflectance of the well-emerged wheat canopies in April. It should also be noted that the limited occurrence of few pure soil pixels within the wheat canopies caused negligible dispersed values in the observed distribution shapes. Toward the end of the cropping season, α 0 values mainly from maximum toward first quartile composed the substantial portion of α 0 density in both biochar and reference plots ( Figure 6). The latter may relate to the pronounced reflectance of a large amount of yellow leaves among the wheat canopies, dominantly contributing to the α 0 retrieval. Furthermore, several soil pixels appeared among the falling wheat leaves and exhibited α 0 values below the identified minimum.
A clear similarity was also observed in the temporal evolution patterns of surface albedo between biochar and reference plots ( Figure 6). α 0 started to slightly increase over March, which is expected because of the emergence of wheat plants increasing the reflectance in α 0 retrieval. Following March, development of wheat plants sharply increased α 0 values over April, reaching the peak at the end of April. α 0 did not appear to be stable toward the end of the cropping season and underwent a slight decrease. This decline in α 0 values is, in fact, indicative of yellow leaves in wheat canopies decreasing the reflectance for α 0 retrieval.

Surface temperature
ST did not appear to be functionally affected by century-old biochar enrichment throughout the cropping season (Figure 7). This result contrasts with our finding in Figure 6 where α 0 decreased across the biochar plots. Despite a limited number of plots, for instance plot 3, indicated a notable effect of biochar on ST at some acquisition dates (Appendix D; Table D1), there were no significant impacts for the majority of the experimental plots across the site. This was accompanied by no differences in average ST (averaged across the eleven plots within the site) between the biochar and reference plots during the monitoring period (Appendix D; Table D1). The sole exception was 29 April when the average ST was found to be affected (p-value < 0.05) by biochar enrichment, indicating an average ST of 20.46 ± 1.67°C as compared to 20.91 ± 1.32°C for the reference plots. This decrease in ST across the biochar patches may relate to the higher soil moisture content in biochar enriched soils. Several previous studies have shown that biochar enrichment increases ST because of its aggravated dark color (e.g. Oguntunde et al. 2008;Ventura et al. 2012). However, our findings reveal that biochar enrichment in soils may not be the sole factor in changing the ST since a set of interactions between α 0 , soil thermal conductivity, soil bulk density, and soil moisture relate to biochar enrichment are needed to influence the ST (Zhang et al. 2013).
ST distribution shape varied among the experimental plots across the site (Figure 7). In addition, distributions were not consistent toward the cropping season. Interestingly, ST distributions were much more clustered around the median at the beginning as well as the end of the cropping season, composing the largest ST density over March and June, respectively. In contrast, distribution shapes were fairly dispersed in April and June. We clearly noticed that plot 5 had a broad range of ST values especially on 24 and 29 April (Figure 7), which can be explained by the presence of a large number of pure soil pixels within plot 5 that had a higher range of ST values. This is somewhat consistent with the α 0 distribution shape of plot 5 on the two aforementioned dates (see Figure 6) at which those pure soil pixels with less reflectance values exhibited low α 0 values. What was noteworthy is the impact of dark pure soil pixels within the biochar patches on 10 and 13 May, where the dark soil of biochar plots remarkably showed higher range of ST values in the distribution shapes (plots 5, 6, and 8; Figure 7). However, there was no difference in ST distributions as a result of pure soil pixels within the reference plots on 10 May, nor was on 13 May.
As expected, the overall temporal pattern of the remotely sensed retrieved ST (Figure 7) matched well with the air temperature values recorded by the onsite weather station (Table 3). For example, the highest range of ST values was observed on 24 June with values ranged from 28°C to 35°C (Figure 7). Similarly, the highest air temperature within the acquisition dates was also recorded on 24 June with an instant air temperature of 30.1°C at the UAV overpass time and a periodic (daily average) air temperature of 23.6°C (Table 3). Such good agreements somewhat validates the utilized approaches in this study to retrieve ST from UAV and L-8 brightness temperature images, respectively, based on the reference groundbased thermal panels (Section 2.3.1) and TES algorithm (Eq. 4-7).

NDVI
NDVI was found to be markedly affected by centuryold biochar enrichment in patches across the site (Figure 8). Interestingly, biochar caused an earlier greening up of wheat plants on the first monitoring date in all the experimental plots. This was accompanied by an average NDVI of 0.45 ± 0.05 across the biochar patches as compared to 0.39 ± 0.05 in the reference plots. This finding corroborates Heidarian Dehkordi et al. (2020a;their Figure 5) where an earlier greening up of chicory plants was observed over the century-old biochar patches of the same agricultural farm. Their UAV-based red-green-blue (RGB) imagery illustrated an average canopy cover of 4.97 ± 2.34% across the biochar plots as compared to 3.46 ± 3.81% in the reference plots on 22 May when the chicory plants started to emerge. Moreover, biochar significantly provoked the development of wheat plants from 20 March until 24 April (Figure 8). Over this period, all the biochar plots indicated significantly higher NDVI values than the reference plots (Appendix D; Table D1), except plot 2 on 16 April as well as plots 8 and 9 on 24 April. Moreover, the average NDVI (averaged across the 11 plots within the site) was significantly higher in the biochar plots as compared to the reference plots from 20 March to 24 April. More precisely, the observed significance of differences declined from March toward April, with p-values being less than 0.001 on 20 March, 28 March, and 16 April and less than 0.01 on 24 April, respectively (Appendix D; Table D1). The average NDVI across the biochar plots (in comparison to the reference ones) was 0.45 ± 0.05 (0.39 ± 0.05) on 20 March (as already stated above), 0.61 ± 0.04 (0.53 ± 0.04) on 28 March, 0.82 ± 0.04 (0.77 ± 0.04) on 16 April, and 0.89 ± 0.02 (0.88 ± 0.02) on 24 April, respectively. The supposed ability of biochar to increase soil water availability (Kerré et al. 2017) as well as nitrogen uptake (Shi et al. 2020) could be the causes of stimulating the crop development (Kerré et al. 2017). It is also worth stating that the increases in NDVI across the biochar patches were mostly due to the significantly higher NIR reflectance values (see Eq. (2)) instead of less red reflectance values, which is more attributable to the soil composition. Following April, NDVI remained steady around the peak till 10 May with wheat canopies being well-developed (Figure 8). From this time onwards, NDVI started to smoothly decrease as wheat leaves began yellowing toward the end of the cropping season. It is worth stating that biochar effect on average NDVI was still rather significant (p-value < 0.05) from the end of April till mid-May (Appendix D; Table D1). The average NDVI of biochar plots (in comparison to the reference ones) was 0.92 ± 0.02 (0.91 ± 0.02) on 29 April, 0.96 ± 0.01 (0.95 ± 0.02) on 10 May, and 0.95 ± 0.01 (0.94 ± 0.01) on 13 May, respectively; there was, however, no impact in June (Appendix D; Table D1).
As expected, the overall temporal evolution of NDVI (Figure 8) was quite consistent with the observed temporal patterns of α 0 (see Figure 6) throughout the cropping season. α 0 increased toward the end of April with the development of wheat canopies increasing the reflectance for α 0 retrieval. Subsequently, α 0 decreased toward the end of the season when wheat leaves began yellowing with less NDVI values. In the study by Zhang et al. (2017), the authors highlighted a strong correlation between leaf area index (as an indicator of vegetation coverage) and α 0 measurements over the biochar-treated soils.
Although NDVI distribution shapes were not consistent over the cropping season (Figure 8), biochar and reference plots indicated rather similar distributions with each other. NDVI distributions were uniform with values between first and third quartiles composing the largest NDVI density at the beginning of the growing season. Over this time, the range of NDVI values was fairly dispersed for all the experimental plots except in plot 9 where the values were more clustered around the median. Such dispersed NDVI distributions are most probably attributable to the occurrence of several pure soil pixels as well as emerged wheat plants within the plots, respectively, revealing NDVI values below the minimum and above the maximum. During the middle of the cropping season, the distribution shapes were getting much more clustered between first quartile and maximum. Surprisingly over this period, well-emerged wheat plants showed no NDVI values above the identified maximum, which is most likely due to the saturation of ratio-based vegetation indices such as the utilized NDVI in the present study (Aklilu Tesfaye and Gessesse Awoke 2021) representing similar range of values for well-developed vegetation pixels. This pattern is in contradiction with the dispersed α 0 values that were observed above the maximum over the same period (see Figure 6). The use of orthogonalbased vegetation indices, which are less sensitive to the saturation effects, such as the weighted difference vegetation index (Clevers 1989), may overcome this  Table D1. issue. Evaluating ET estimates through ETLook by using orthogonal-based vegetation indices as a proxy of vegetation coverage should be a topic for future research, which falls behind the scope of the present study. It should also be stated that the occurrence of few pure soil pixels within the wheat canopies led to dispersed NDVI values below the identified minimum over the middle of the cropping season ( Figure 8). Toward the end of the season, NDVI distribution shapes for most of the plots were becoming rather uniform once again, with values mainly from maximum to first quartile composing the substantial portion of NDVI density. This could be attributed to the yellow leaves among the wheat canopies that had less reflectance and, hence, a lower range of NDVI values. In addition, some soil pixels that are visible among the falling wheat leaves showed NDVI values below the identified minimum. The latter is in line with the observed patterns in α 0 distribution shapes (see Figure 6) at the end of the cropping season.

Biochar effects on actual evapotranspiration
In general, century-old biochar enrichment in patches across the site showed no significant change in actual ET rates either at the beginning or at the middle of the cropping season (Figure 9). Neither was there any impact toward the end of the season. Over the first five monitoring dates, all the individual plots exhibited no significant difference in actual ET as a result of biochar enrichment (Appendix D; Table D1), except plots 1 and 8 only on 29 April where the observed differences were rather significant (p-value < 0.005). Following April, the majority of the experimental plots indicated significant differences in actual ET rates between biochar and reference plots on 10 May (plots 1-6 and 11; Appendix D; Table D1). From these plots, there was, however, only a significant impact in plots 5 and 11 two days later on 13 May. The rest of these plots might be highly influenced by (i) the large amount of rainfall (6.25 mm as recorded by the on-site weather station) on 11 May as well as (ii) the sharp increases of solar radiation from 10 to 13 May (see Table 3) that had jointly overcome the impact of biochar enrichment on actual ET. No experimental plot showed any significant effect of biochar on actual ET on either 24 June or 27 June except plot 5 on 27 June (Appendix D; Table D1).
In addition, there were no differences in average actual ET (averaged across the 11 plots within the site) between the biochar and reference plots throughout the monitoring period except on 10 May (Appendix D; Table D1). The latter was accompanied by an average actual ET of 1.66 ± 0.44 mm/day across the biochar patches as compared to 1.68 ± 0.45 mm/day in the reference plots; the observed difference was statistically significant (p-value < 0.005; Appendix D; Table  D1). Although several studies have shown that biochar amendment can (nonsignificantly) increase ET rates in vegetable crops due to increases in soil water retention capacity (e.g. EUROCHAR 2014; Ahmed et al. 2019;Fischer et al. 2019), our remotely sensed analysis exhibited no significant influence of biochar enrichment on actual ET estimates. The latter is most likely attributable to the rich soil quality in the agricultural farm of the present study (Hedarian Dehkordi et al. 2020b;Burgeon et al. 2021), overcoming the impact of biochar enrichment on actual ET rates. Similarly, another study by Alburquerque et al. (2013) reported no significant influence of biochar on ET across durum wheat planted on a loamy sand soil. This corroborates the recent finding of Baiamonte, Minacapilli, and Crescimanno (2020)where biochar application had no impact on ET across a sandy soil planted with wheat.
Distribution shapes of actual ET varied throughout the cropping season and also among the experimental plots across the site (Figure 9). Actual ET distributions were clustered around the median, close to zero mm/day, on the first monitoring date, though plots 5 and 9 exhibited small rates of actual ET (< 1 mm/day). Distributions were yet rather clustered around the median (< 1 mm/ day) composing the largest density of actual ET for most of the plots on 28 March and 16 April, although plots 1 and 2 indicated fairly dispersed distributions. Such clustered actual ET distributions at the beginning of the season (Figure 9) is consistent with the clustered distributions of ST observed from March till mid-April (see Figure 7). Actual ET rates then sharply increased toward 24 April (Figure 9), which is in accordance with the high development of wheat canopies (see Figure 8) that also revealed high values of α 0 (see Figure 6). This sharp increase of actual ET occurred because the air temperature was extremely warmer on 24 April with a daily average of 15.4°C as compared to 11.6°C on 16 April (Table 3). As such, wheat plants were opening-up the stomata to release water vapor that increases the transpiration rates (Allen 1999). It is worth mentioning that actual ET rates on 24 April (Figure 9) were mostly clustered between first quartile and maximum, and exhibited no values above the identified maximum, which we postulate to be related to the saturation of NDVI (Section 3.2.3; Figure 8). As such, this would be of particular interest to examine ETLook ET estimates by using orthogonal-based vegetation indices instead. There was a gradual decrease in actual ET rates on 29 April, which is expected because of the recorded rainfalls of 3.25 and 8.5 mm on 25 and 28 April, respectively. Actual ET distributions were becoming rather uniform with values between first and third quartiles composing the largest density (Figure 9) from the end of April till mid-May. Actual ET then reached the peak on 24 June (Figure 9), which is consistent with the observed patterns in ST (Figure 7) and recorded air temperature values (Table 3), indicating 24 June as the warmest monitoring date. There was a slight decrease in actual ET at the last acquisition (Figure 9), which is in line with the observed patterns in the three examined remotely sensed ET elements as α 0 (Figure 6), ST (Figure 7), and NDVI ( Figure 8).
To further assess the impacts of century-old biochar enrichment, correlation analysis was performed to investigate any potential relationship between biochar and ET. For this, scatterplots of soil organic carbon (as representative of biochar) for each experimental plot, available from Heidarian Dehkordi et al. (2020a), against their ETLook actual ET estimates (this study) as well as grain humidity, available from Hedarian Dehkordi et al. (2020b), were generated (Appendix E; Figure E1). There was a weak correlation between soil organic carbon and ET estimates, exhibiting correlation coefficient (resp. p-value) of 0.43 (0.04) as shown in Appendix E; Figure E1. This corroborates the identified relationship between soil organic carbon and ET rates across the Southern Great Plains grassland, United States by Homann, Kapchinske, and Boyce (2007). There was, however, no correlation between soil organic carbon and grain humidity (r = 0.13, p-value = 0.57; Appendix E; Figure   Figure 9. Violin plots comparing the actual evapotranspiration (ET) between century-old biochar and reference plots. Inside boxplots indicate the percentiles of the distribution values at 25% (i.e. first quartile), 50% (i.e. median, represented by circles), and 75% (i.e. third quartile) surrounded by minimum and maximum values. Width of the violin plots represents the distribution of actual ET values. p-value results of each plot (comparing biochar vs. reference) for varying dates are contained in Appendix D; Table D1. E1). A similar evaluation was made by Heidarian Dehkordi et al. (2020a) across the same agricultural farm, who revealed that biochar enrichment did not cause substantial changes in chicory dry matter content. Hence, there exists a need for further work to study, in particular, such elusive relationships concerning biochar effects on ET.
In this paper, we focused on high spatiotemporal monitoring of century-old biochar effects on ET as one of the primary studies to the best of our knowledge by fusing remotely sensed UAV and L-8 images. We demonstrated a strong ability of AWT to fuse UAV and satellite images with large differences in the spatial resolutions. However, more effort is still needed to improve the associated persistent errors in AWT while fusing UAV and satellite images. We retrieved ST from UAV and L-8 brightness temperature images, respectively, based on the reference ground-based thermal panels and TES algorithm. α 0 and NDVI were computed based on the reflectance spectra of multispectral sensors on board of the UAV and L-8 platforms. More work is yet required to enhance α 0 retrievals from the multispectral sensors on board of UAVs. In the present study, we chose only to focus on modeling ET by assimilating remote sensing images into the ETLook model. Future researches may consider verifying remote-sensing-based ET estimates with in-field techniques such as eddy covariance (Parent and Anctil 2012). We applied the presented approach to an agricultural farm enriched with multiple biochar patches. The fact that in the present study we did only consider one particular field characterized by century-old biochar patches does not allow us to provide generalized findings. Since the agricultural fields with century-old biochar patches are distributed rather disperse across Belgium, future researches should be aimed at extending the scale of the study across multiple sites using manned aircrafts and high-resolution satellite to generalize the findings. In light of climate-smart agriculture, this approach is applicable to any study site characterized with other specific precision agricultural practices and, also, many other satellite images given their easy access across the globe.

Conclusion
It may be challenging to fuse remotely sensed images with large differences in the spatial resolutions. In this study, we illustrated a strong ability of AWT to fuse L-8 (GSD of 30 m) and UAV images (GSD of 7 and 3.7 cm for thermal and multispectral images, respectively) as one of the primary studies to the best of our knowledge. We computed α 0 and NDVI based on the reflectance spectra of MicaSense RedEdge-M on board of the UAV and L-8 OLI multispectral sensors. We also retrieved ST images from FLIR Vue Pro R on board of the UAV and L-8 TIRS brightness temperature images, respectively, based on the reference groundbased thermal panels and TES algorithm. Highresolution UAV derivatives together with AWTsharpened L-8 images, and in combination with meteorological data, were used in the ETLook model to estimate ET rates across an agricultural farm enriched with century-old biochar. Fusion of L-8 and UAV images generated sharpened L-8 images that were all significantly correlated (p-value < 0.001) with their corresponding coarse resolution images at each fusion date, though small amount of persistent errors were observed. We also demonstrated that spatial details can be well preserved when fusing UAV and satellite images with extremely different spatial resolutions through AWT. Such promising results can help to establish in the future more robust image fusion techniques for fusing UAV and satellite images in light of climate-smart agriculture. Our approach is also applicable to investigate other farms with specific precision agricultural treatments.
Our high spatiotemporal analysis indicated a significant decrease in α 0 across the biochar patches especially over the early development stages of winter wheat. In addition, biochar enrichment significantly (i) provoked an earlier greening up of wheat plants and (ii) stimulated the development of wheat canopies toward the middle of the cropping season. Aforementioned biochar impacts on α 0 and crop development seem to be vanished toward the end of the cropping season, which was most likely attributable to the dense wheat canopies covering the aggravated dark color soil across the biochar patches. In contrast, ST did not appear to be affected by biochar enrichment either at the beginning or toward the end of the cropping season. Neither was there any impact of biochar enrichment on actual ET estimates throughout the cropping season. Thus, there exists a need for future studies to verify such remotely sensed findings in light of biochar effects on ET elements with infield techniques.

Highlights
• Additive wavelet transform allows the fusion of UAV and satellite images • Century-old biochar has no impact on the modeled actual evapotranspiration • Biochar significantly decreased surface albedo at the beginning of cropping season • Biochar showed no change in the remotely sensed surface temperature • Biochar significantly stimulated plant development at the beginning of the season TIRS thermal bands 1 and 2 are initially designed at 100 m but have been resampled to 30 m in order to match OLI multispectral bands. Figure B1. Scatterplots of air temperature recorded by the on-site weather station against the average remotely -sensed surface temperature (averaged across the field) retrieved from either UAV or Landsat-8 (L-8) images for all the acquisition dates. Blue and green dots refer to UAV and L-8, respectively.

Appendix D
Output of statistical t-test examining whether there are significant differences in evapotranspiration elements between century-old biochar and reference plots; p-values are color-coded with light, moderate, and dark blues indicating low (< 0.05), moderate (< 0.01), and high (< 0.001) levels of significance, respectively. α 0 , ST, NDVI, and Actual ET refer to surface albedo, surface temperature, normalized difference vegetation index, and actual evapotranspiration, respectively.