Drone-Based Hyperspectral and Thermal Imagery for Quantifying Upland Rice Productivity and Water Use Efficiency after Biochar Application

Abstract

Miniature hyperspectral and thermal cameras onboard lightweight unmanned aerial vehicles (UAV) bring new opportunities for monitoring land surface variables at unprecedented fine spatial resolution with acceptable accuracy. This research applies hyperspectral and thermal imagery from a drone to quantify upland rice productivity and water use efficiency (WUE) after biochar application in Costa Rica. The field flights were conducted over two experimental groups with bamboo biochar (BC1) and sugarcane biochar (BC2) amendments and one control (C) group without biochar application. Rice canopy biophysical variables were estimated by inverting a canopy radiative transfer model on hyperspectral reflectance. Variations in gross primary productivity (GPP) and WUE across treatments were estimated using light-use efficiency and WUE models respectively from the normalized difference vegetation index (NDVI), canopy chlorophyll content (CCC), and evapotranspiration rate. We found that GPP was increased by 41.9 ± 3.4% in BC1 and 17.5 ± 3.4% in BC2 versus C, which may be explained by higher soil moisture after biochar application, and consequently significantly higher WUEs by 40.8 ± 3.5% in BC1 and 13.4 ± 3.5% in BC2 compared to C. This study demonstrated the use of hyperspectral and thermal imagery from a drone to quantify biochar effects on dry cropland by integrating ground measurements and physical models.

1. Introduction

Land surface variables are required for modeling of carbon, water, and energy exchanges between land and atmosphere. Miniature multispectral and thermal cameras onboard lightweight unmanned aerial vehicle (UAV) have been widely used for providing spatially explicit land surface variables and as an efficient tool supporting precision farming and crop management e.g., [1,2,3]. For example, Kandylakis et al. [4] used UAV-borne multispectral and shortwave infrared cameras to estimate LAI of vineyards, and assess water stress conditions using individual relations gained from UAV data and in-situ measurements via statistical regression. Chen et al. [5] used data from UAV-borne multispectral and thermal cameras to evaluate cotton canopy water stress. Infrared thermography has been used to diagnose crop water deficit and drought conditions [6,7,8].

Besides multispectral, hyperspectral imaging onboard an aircraft with dozens or even 100+ continuous narrow bands of high spectral resolution is also attracting research interests. However, the previously high-cost and heavyweight hyperspectral imager using motion scanning mechanics limited the hyperspectral technique to few studies on expensive airborne sensing (e.g., [9,10]). The recent advances in miniature hyperspectral snapshot imaging technology bring new opportunities for low-cost UAV-borne remote sensing [11]. These snapshot imagers use non-moving parts such as image slicer mirror or filter-on-chip technologies [12], and minimize moving artifacts during sensing to gain satisfactory spectral and spatial accuracy. The current miniature hyperspectral snapshot cameras are promising for monitoring crop growth and retrieving crop biophysical variables for precision agriculture from UAVs.

Hyper- or multiple spectral data from UAV sensing are usually analyzed in three approaches to estimate land surface properties. (1) Vegetation indices simply combining different band reflectance are often used to estimate a certain land surface variable based on empirical relations. For example, the normalized difference vegetation index (NDVI) formulated using red and near-infrared (NIR) bands is widely used to indicate land surface greenness [13]. Other narrow hyperspectral band combinations are used for leaf chlorophyll estimation [10]. (2) Machine learning (ML) techniques. ML is increasingly used to estimate land surface variables based on large field measurements, including for example partial least squared regression [14], artificial neural network [15], and random forest [16]. These empirical methods need site-specific land surface variables for model training and the trained model may not be suitable for other sites. (3) Physical model inversion. A radiative transfer model, such as PRO4SAIL [17], can simulate hyperspectral reflectance close to the field measured reflectance with optimal inputs. The inverse of the radiative transfer model based on field hyperspectral measurements can obtain leaf and canopy variables such as pigment content and leaf area index (LAI) by minimizing the difference between the modelled and the measured reflectance in all spectral wavelengths in model optimization process. Such a physically based model does not necessarily require field-measured land surface variables for parameter calibration, especially for common crops like rice, being more robust than vegetation index methods and avoiding overfitting issues in ML methods to retrieve plant and soil parameters.

These techniques can be used to assess crop-atmosphere exchanges when manipulating soil conditions, for example, using biochar for soil amendment. Biochar is an old technique used for soil amendment in the Amazon region centuries ago [18,19]. The porous charcoal has a large surface area that can bind and retain soil water and nutrients, and reduce nutrient leaching and volatilization loss, thereby helping to increase long-term plant water availability and soil fertility [20,21]. Therefore, biochar is considered a valuable soil amendment to increase crop productivity and resilience to drought [21]. Because charcoal itself is decomposition-resistant, it can contribute to soil carbon sequestration and atmospheric carbon reduction [22], and is thus regarded as a feasible negative-emission strategy to mitigate climate change [19].

Incorporating biochar into good agricultural practices requires a better understanding of its effects on crop productivity, soil water usage, and land surface energy balance. Using a conceptual model, Fischer et al. [21] showed that biochar amendment—as evident by shifting the soil water retention curve—should increase soil water availability under water-limited conditions. As a result, the long-term mean evapotranspiration (ET) rate is expected to increase after biochar amendment, which in turn should increase the gross primary productivity (GPP) due to the joint regulation of transpiration and photosynthesis by stomata. However, the impacts of biochar on plant water use efficiency (WUE = GPP/ET, the rate of GPP per unit water consumption) are not clear. Moreover, ET is not always increased in biochar-amended fields, in contrast to model predictions [21].

These divergent effects of biochar on crop growth and soil hydraulic properties highlight our limited understanding of the effects of biochar on crop growth and soil water usage [21,23,24,25]. The apparent gap between the theory and empirical evidence on biochar effects may be also due to extrapolating documented biochar effects on soil properties from discrete in situ point measurements to the field or farm scale. Furthermore, as there are concurrent impacts of biochar on soil and plant properties, synthesizing them by combining both spatial and point measurements may provide a broader picture of crop responses after biochar application.

This study explores the potential of UAV-borne hyperspectral and thermal sensing for a spatially improved understanding of the functioning of biochar in tropical dryland. We investigate the effect of biochar in a field experiment with upland rice (Oryza sativa L.) in the North Pacific of Costa Rica to assess changes associated with energy, productivity, soil water availability, ET, and WUE. This is a drought-prone region, projected to experience more prolonged and severe droughts in the near future. A four-day field campaign was conducted using hyperspectral and thermal imaging cameras onboard a hexacopter at the end of the rice growing season. Combining the UAV snapshot images of surface variables, in situ measurements, and physical models, the study aims to (1) assess relative variations in dryland rice productivity, soil moisture content, and water use efficiency after biochar application, (2) test the applicability of drone-based hyperspectral and thermal imagery in cropland monitoring, and (3) develop a methodology combining data from UAV sensing, ground measurements, and models to support agricultural management.

2. Materials and Methods

2.1. Biochar Experiment Plots

The experiment was conducted on an about 8 × 20 m site parcel, located at the Enrique Jiménez Núñez Experimental Station in Costa Rica (10.3436° N, 85.1353° W; 17 m asl, Figure 1). The North Pacific of Costa Rica is a drought-prone region, in which severe dry events are mostly related to the warm phase of the El Niño-Southern Oscillation phenomenon [26,27]. The region features increasing aridity and rainfall reduction trends, which have been projected to exacerbate along with global warming, threatening both water resources availability and agricultural activities [26,28]. In the long-term, the region is characterized as a tropical savannah climate with a marked dry season from mid-November to April. The long-term mean annual rainfall is about 1547 mm/yr and the mean annual temperature 27.4 °C (100-year period). The soil at the experimental site has a clay loam texture. The experiment consisted of three groups in triplicate separated by 0.8 m wide corridors. Due to logistic reasons, the experimental plots were set in a fairly small area without using a randomized block design. Two treatment groups used local Guadua bamboo (Guadua angustifolia) biochar (BC1) and Taiwan sugarcane (Saccharum officinarum) filter cake biochar (BC2) respectively, and the control group (C) had no biochar addition. BC1 and BC2 biochar were incorporated in the top 20 cm of the soil on February 13 and July 18 of 2018 respectively due to their different arriving times. Upland rice seeds (variety Palmar 18) were sown on July 18, 2018, and were irrigated and fertilized to encourage germination and seedling growth following the common local farming practice. Thereafter, the water resource for rice growth was solely from precipitation (rain-fed). Harvest took place on 21 November 2018 and indicated the end of the experiment. See [29] for details of biochar experimental design.

2.2. UAV Field Campaign

The UAV field campaign was carried out using a hexacopter (DJI M600 Pro, Dajiang, China) on November 14, 16, 19, and 21, 2018 before harvest. These days were selected for the field campaign for the following reasons: (1) These days were at the start of the dry season in the region and the biochar effects on soil water retention could be easily tested. (2) It was the phenology transition period during which the rice leaves changed from green at the beginning of the campaign to yellow (senescence) at the end of the campaign, likely creating large spatial and temporal variability in rice water use. (3) The LAI and GPP at the time of the campaign are expected to reflect the accumulated growth of the whole growing period and consequently, allowing the integrated effect of biochar to be quantified. (4) A field campaign spanning a whole rice growing period was logistically unfeasible.

The Dajiang DJI M600 Pro hexacopter has a payload capacity of 6 kg and a flight time of 35 minutes under minimal payload. In the campaign, the UAV was flying at a height of about 30 m above the ground. The UAV position was recorded by the three sets of GNSS and IMU systems inbuilt with the drone. A gimbal (Gremsy T7, Gremsy, Vietnam) was used to enable the camera to consistently viewing the nadir direction. We used a hyperspectral camera (FireflEYE 185ST, Cubert, Germany) and a thermal camera (324×256 Pixels, Tau2 324, FLIR, USA) with the attained ground resolutions of 4.5 mm and 2.25 cm respectively to investigate the plant productivity, soil, and surface properties. The Cubert hyperspectral camera has 50 × 50 pixels, with each pixel spanning from 450 to 950 nm in 138 wavebands. The camera has an extra panchromatic band of 1000 × 1000 pixels to facilitate sharpening the hyperspectral bands. It has been shown that the best spatial and spectral accuracy can be achieved from Cubert panchromatic sharpening at a flying height of 30 m [11]. In our application, we adopted the default pan-sharpening function in the Cubert-pilot software. During the flight, a spectrometer (Flame S VIS-NIR, Ocean Optics, Dunedin, FL, USA) was mounted on the drone to simultaneously collect spectral irradiance.

It should be noted that the hyperspectral camera images had about 55% of side overlap and 40% of forward overlaps in the flight. The overlaps were smaller than the suggested figures (60% of side overlap and 80% of forward overlap) in [30] due to the small field-of-view (FOV) of the Cubert hyperspectral camera (13° × 13°). On the one hand, a small FOV camera with nadir viewing has a great advantage of minimizing the bi-directional reflectance distribution function (BRDF) effect [31], which otherwise needs to be adjusted carefully for quantitative remote sensing [32]. On the other hand, the small FOV creates difficulties in further camera alignment processing when the image overlap is not large enough as afore-mentioned in [30] and consequently the cloud points are not dense enough for estimating camera locations using Agisoft Metashape software. Otherwise, much more hyperspectral images should be collected, heavily increasing the field workload. To attain accurate camera alignment with small hyperspectral image overlaps and without increasing field workload, we used an extra RGB camera (Sony DSC-RX100, Sony, Japan) with a larger FOV (64.8° × 45.9°), and therefore at least 80% overlaps in RGB images along with hyperspectral capturing. We described the use of an RGB camera to help align hyperspectral camera images accurately in the following section.

2.3. Ground Measurements

Hyperspectral reflectance signatures of rice canopy and soil background were measured respectively using a spectroradiometer (FieldSpec HandHeld 2, ASD, Boulder, CO, USA) and a white Spectralon reference panel (Labsphere, North Sutton, NH, USA). Canopy reflectance was measured on each plot three times. The average spectral reflectance was used to verify Cubert hyperspectral camera measurements.

The leaf gravimetric water content was measured using a destructive method. Three samples of whole rice leave (including leaf blade and sheath, and the stem under the sheath) were taken in each plot on 19 November 2018. Each sample was weighed fresh, cut into small pieces, and heated in an oven at 70 °C for 24 hours in a laboratory. The leaf water mass was determined as the difference between the fresh and the dry leaf mass. The leaf gravimetric water content was defined as the ratio of leaf water mass to fresh mass (g/g). Leaf nutrient content was also determined in the laboratory using the total digestion method on leaf samples. The macro-elements of nitrogen (N), phosphorous (P), and potassium (K) concentration are determined as a percentage (%) of the dry matter. One soil sample was collected from each plot. The total nitrogen concentration in each sample was tested using a CN 628 Dumas analyzer (LECO, St. Joseph, MI, USA) in the laboratory in National Agricultural Technology Institute (INTA), Costa Rica.

An automatic meteorological station (Vaisala WT520, Vaisala, Finland) was installed at a height of 1.5 m above the ground to continuously monitor precipitation, wind speed and direction, air temperature, relative humidity, and atmospheric pressure (Figure 1). Each plot was instrumented with two in situ sensors at 15 cm below the surface (Figure 1), including one sensor (Decagon GS3, METER Group, Pullman, WA, USA) measuring volumetric soil water content, soil electrical conductivity, and soil temperature, and another sensor (Decagon MPS6, METER Group, USA) measuring soil matric potential and soil temperature. The sensors were connected to a data logger (CR1000, Campbell Scientific, Logan, UT, USA) to collect the data at 30-minute intervals from 18 July to 21 November 2018.

2.4. Data Processing

2.4.1. Radiometric Correction: Digital Number (DN) to Physical Values

The DN values of pan-sharpened high-resolution hyperspectral image were transformed to absolute radiance (W/m²/nm/sr) values using per-band per-pixel sensitivity factors that were determined in a photonics laboratory of the Technical University of Denmark before the campaign. The spectral irradiance (W/m²/nm) measured by the on-flight spectrometer was converted to per-band per-pixel irradiance value for the Cubert camera using the spectral response function of each pixel that was determined in the same laboratory. See [33] for details about the Cubert camera calibration. The reflectance of each pixel was calculated as:

$$\rho =\pi \frac{radiance}{irradiance},$$

(1)

The brightness temperature was calculated using calibration parameters and the inverse of Planck’s Law following Köppl [34]:

$$T_b={\log }^{-1}\left(\frac{k_1}{DN+k_2+k_3{T}_c}+1\right),$$

(2)

where k₁, k₂, and k₃ are calibration parameters, DN is the digital number of FLIR camera images, and T_c is the camera detector core temperature.

2.4.2. Image Orthorectification and Mosaicking

After transforming all the images into hyperspectral reflectance and land surface temperature, the Agisoft Metashape software (Agisoft, St. Petersburg, Russia) was used to generate orthorectified and mosaicked (orthomosaic) images that covering the whole experiment rice field. In hyperspectral image processing, we pooled both RGB and hyperspectral camera images for camera alignment to guarantee cloud points dense enough for an accurate reconstruction of Cubert hyperspectral camera locations. To avoid processing of large volume of hyperspectral data iteratively, only panchromatic gray images were used in the alignment process. This allowed reducing the computational time on a desktop computer with a 16-core 3.40 GHz CPU and two GPUs (GTX1080 Ti, Nvidia, Santa Clara, CA, USA). The DEM and orthomosaic hyperspectral images were generated. The DEM was used to estimate rice canopy heights by subtracting the average surrounding soil background elevation from the canopy elevation. The hyperspectral images were used in the subsequent analysis.

2.4.3. NDVI and Variations in Gross Primary Productivity (GPP)

The orthomosaic hyperspectral reflectance values were corrected using an empirical line correction method by selecting two pseudo-invariant features in the scene: an unshaded fixed white instrument box and a patch of dark bare soil that were identified from the composite true-color hyperspectral images. Then the LAI ([-]) and canopy chlorophyll content (CCC, [g/m² ground area]) were retrieved from inversion of PRO4SAIL model (available at http://teledetection.ipgp.jussieu.fr/prosail/, accessed on 8 April 2021) with the corrected hyperspectral reflectance, pixel by pixel at 2.25 cm resolution (resampled from 4.5 mm to match thermal images). The NDVI was calculated using the NIR band (843 nm) and the red band (664 nm): NDVI = (NIR − red)/(NIR + red). The NDVI has been shown to have a linear relationship with the fraction of photosynthetically active radiation (fAPAR) absorbed by plant canopy [35,36]. The rice plant variables GPP, LAI, NDVI, and CCC were used as rice growth indicators in this study. The GPP can be simulated using a light use efficiency model [37] with the photosynthetically active radiation (PAR) absorbed by canopy chlorophyll [38]. Based on these relations and assuming other factors constant across biochar treatments, we inferred the relative GPP variation ($\mathsf{\Delta }GPP/GPP$) after biochar application from the variation in NDVI and CCC using a propagation approach [39] (see Supplementary Material for further details):

$$\frac{\mathsf{\Delta }GPP}{GPP}=\frac{\mathsf{\Delta }NDVI}{NDVI}+\frac{\mathsf{\Delta }CCC}{CCC}.$$

(3)

2.4.4. Estimation of Shortwave Surface Albedo and Directional Emissivity

The shortwave surface albedo was estimated from the PRO4SAIL model results following the method in [40]:

$$\alpha =\frac{{\alpha }_{model}}{{\overline{R}}_{model} } \overline{\mathrm{R}},$$

(4)

where $\overline{\mathrm{R}}$ and ${\overline{R}}_{model}$ are the Cubert-measured and PRO4SAIL-modelled mean blue-sky directional reflectance from wavelength 450 nm to 950 nm respectively; ${\alpha }_{model}$ is PRO4SAIL-modelled shortwave surface albedo.

The directional land surface emissivity was calculated from the canopy gap fraction following François [41]:

$${\epsilon }_s\left({\theta }_V\right)=1-b\left({\theta }_V\right)\left(1-{\sigma }_f\right)\left(1-{\epsilon }_g\right)-c\left[1-b\left({\theta }_V\right)\left(1-{\sigma }_f\right)\right]\left(1-{\epsilon }_v\right),$$

(5)

where ${\theta }_V$ is the thermal camera viewing direction, with ${\theta }_V=0$ for nadir-viewing; ${\epsilon }_v$ and ${\epsilon }_g$ are leaf emissivity and bare ground emissivity respectively, and in this work ${\epsilon }_v=0.98$ and ${\epsilon }_g=0.94$ were used; c is a cavity factor, accounting for the volumetric multiple scattering inside rice canopy, and c = 0.3 was used for nadir viewing; b(${\theta }_V$) is the directional gap fraction at viewing direction ${\theta }_V$; ${\sigma }_f$ is called shading factor and $1-{\sigma }_f$ is the hemispherical transmittance. Both b(${\theta }_V$) and ${\sigma }_f$ were estimated from LAI. See Supplementary material for further details. The estimation of hemispherical-directional emissivity in Equation (7) has been shown to have high computation efficiency and high accuracy [41].

2.4.5. Estimation of Land Surface Temperature

The land surface temperature (T_s) was estimated from the brightness temperature with a correction for the reflected longwave radiation by land surface:

$$T_s={\left[\frac{T_b{}^4-\left(1-{\epsilon }_s\right){\epsilon }_a{T}_a{}^4}{{\epsilon }_s}\right]}^{1/4},$$

(6)

where T_a is the air temperature from the weather station, [K]; ε_s is the land surface emissivity estimated from Equation (5) for ${\theta }_V=0$; ε_a is the atmosphere emissivity, estimated from air temperature and relative humidity following Prata [42]. For details see [43].

2.4.6. Energy Components and Variations in Water Use Efficiency (WUE)

Net radiation $R_n$ (W/m²), the difference between incident and outgoing radiation energy, was calculated as:

$$R_n=\left(1-\alpha \right)R_s+\sigma {\epsilon }_s{\epsilon }_a{T}_a^4-\sigma {\epsilon }_s{T}_s^4,$$

(7)

where $R_s$ is the incoming solar radiation measured by the weather station at the site, W/m²; $\alpha$ is the land surface albedo derived from Equation (4); σ = 5.67 × 10⁻⁸ W/m²/K⁴, Stefan–Boltzmann constant. $T_s$, $T_a$, ${\epsilon }_s$, and ${\epsilon }_a$ have the same definitions as in Equation (6).

Latent heat flux ($\lambda ET$, (W/m²), the ET in energy form) was estimated as the residual component of the surface energy balance:

$$\lambda ET=R_n-H-G,$$

(8)

where $R_n$ is the net radiation from Equation (7), H [W/m²] is sensible heat flux, and G [W/m²] is the heat storage flux in soil and plants. G is estimated from net radiation $R_n$ and fractional vegetation cover $f_c$ using an empirical equation in [44]:

$$G=R_n\left[0.315\left(1-f_c\right)+0.05f_c\right],$$

(9)

where $f_c$ is calculated from the aforementioned directional gap fraction at the nadir-viewing direction (see Supplementary material): $f_c=1-b\left(0\right)=1-e^{-0.245LAI}$.

The sensible heat flux H was estimated from the temperature difference between the land surface and the air using the bulk transfer equation:

$$H=\rho C_p\frac{T_s-T_a}{r_a+r_{ex}},$$

(10)

where $\rho$ is the air density, $C_p$ is the specific heat of air at constant pressure, and $r_a$ is aerodynamic resistance to heat transfer, estimated from wind speed and the average height of rice canopy following [45]; $r_{ex}$ is an extra resistance added to $r_a$ for correcting the difference between radiometric temperature and aerodynamic temperature in sensible heat flux calculation. $r_{ex}$ is approximated by $5.75/u_{*}$ [46,47], where $u_{*}$ is friction velocity. Note that Equation (10) uses the surface temperature as a proxy of aerodynamic temperature, even though $r_{ex}$ is introduced to account for this, underestimates on $r_{ex}$ may result in an overestimation of the sensible heat flux for bare soil area from an overestimated aerodynamic temperature, and subsequently result in negative evapotranspiration in Equation (10). When this happened, we assigned a zero value for evapotranspiration (reasonable for dry bare soil) and recalculated sensible heat flux H using Equation (8). The available energy (AE, (W/m²)) was estimated as the difference between net radiation $R_n$ and the ground heat flux G, i.e., $AE=R_n-G$. Evaporative fraction (EF, [-]) is the latent heat flux normalized by the available energy: $EF=\lambda ET/AE$.

The variations in WUE after biochar application were inferred from differences between the GPP variation and the ET variation using the variation propagation method again (see Supplementary Material for further details):

$$\frac{\mathsf{\Delta }WUE}{WUE}=\frac{\mathsf{\Delta }GPP}{GPP}-\frac{\mathsf{\Delta }ET}{ET}.$$

(11)

2.4.7. Soil Moisture Content and Soil Matric Potential Estimation

We estimated volumetric soil water content from UAV data using a temperature-vegetation dryness index (TVDI [48], a triangle method by identifying the pixel dryness from the space of temperature variations between surface and bulk atmosphere ($\Delta T=T_s-T_a$) and fractional vegetation cover $f_c$, and then related TVDI to soil moisture content in the root zone as shown in [43]:

$$TVDI=\frac{\Delta T-\Delta T_{min}}{\Delta T_{max}-\Delta T_{min}}=1-\frac{\theta -{\theta }_{WP}}{{\theta }_{FC}-{\theta }_{WP}},$$

(12)

where $\mathsf{\Delta }{T}_{max}$ and $\mathsf{\Delta }{T}_{min}$ are the dry and wet edges of the ΔT-f_c triangle, estimated following [49]; and $\theta$ (% in m³/m³) is the soil moisture content to be estimated. ${\theta }_{FC}$ and ${\theta }_{WP}$ are soil moisture content at field capacity and wilting point respectively, estimated from soil water retention curve at soil matric potential of −0.05 and −1.5 MPa respectively (see [29] for details of water retention curve estimation using in situ measurements). The spatially explicit soil matric potential (ψ, (MPa)) was then estimated from the UAV-derived soil moisture content using the Van Genuchten model [50]:

$$\psi =\frac{1}{{\alpha }_v}{\left(\frac{1}{{\mathsf{\Theta }}^{1/m}}-1\right)}^{1/n},$$

(13)

where $\mathsf{\Theta }=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$, $m=1-\frac{1}{n}$, and ${\theta }_r$, ${\theta }_s$, ${\alpha }_v$, and n are the water retention curve parameters reported in [29].

2.4.8. Statistical Analysis

The biophysical and land surface variables derived from UAV hyperspectral and thermal sensing were summarized for each plot to assess their relative changes in biochar amendment vs. control plots. The sunlit soil, shaded soil, and rice canopy pixels were identified using an unsupervised classification [51]. A k-means clustering method was used in the classification. After testing with different numbers of classes and bands, we found that using 5 clusters and 3 bands (Band 12, 29, and 104, for the false composite image in Figure 2) can effectively separate the aforementioned three types of pixels. Sunlit and shaded leaves were not separable probably due to their mixed spectral signatures and the foliage transparency (Supplementary Figure S1).

Sunlit soil pixels at plot corridors were extracted for determining bare soil albedo and surface temperature and comparing their differences between the biochar groups the control group. Rice canopy pixels were extracted for analyzing the biochar effects on root zone soil moisture, canopy biophysical variables, and canopy surface energy components. To reduce edge effects [52], the pixels within 15 cm to borders of each plot, roughly corresponding to the edge rows of rice plants, were discarded. The valid pixels extracted from each plot were aggregated to three sub-plots (Figure 1) for further statistical analysis, to avoid spatial autocorrelation from fine resolution.

We analyzed four days’ repeated UAV results over three treatment groups with triplicate using repeated measures analysis of variance (ANOVA) [53] by the free statistic software JASP (available at https://jasp-stats.org/, accessd on 8 April 2021). The Greenhouse-Geisser correction was used if sphericity had been violated in the repeated measures. The marginal mean and standard error of each variable were calculated for each group of three plots and four days considering the degrees of freedom within- and between-subjects. The mean variable differences between the treatment groups and the control group were analyzed using Tukey’s honest significant difference test [54]. For the one-time measurement of leaf and soil samples, one-way ANOVA was used to analyze if their mean values were different. The four days’ repeated UAV results were compared with the in-situ measured soil water content and soil matric potential using violin plots and descriptive plots with mean and 95% confidence interval. A simple Pearson correlation was used to test the relationships between key rice growth variables and soil water availability.

3. Results

3.1. Hyperspectral and Thermal Imagery

We obtained orthomosaic Cubert hyperspectral images of 4.5 mm resolution in this field campaign. Figure 2a shows an example of a false-color composite image taken on 14 November 2018. Individual rice leaf blades are visible. Figure 2b shows the comparison of hyperspectral reflectance curves measured by Cubert camera (after the empirical line correction) and ASD spectroradiometer on rice canopy. The two curves share a similar spectra pattern and show the reliability of UAV-based hyperspectral sensing in the campaign. Note that the last few bands of Cubert camera have large variations and were discarded from further analysis following the manufacturer’s recommendation. In four days of measurements, the assumed static dark soil and white metal box objects had up to 30% deviations from day to day (Figure 2c,d), indicating the necessity for correcting the systematic day-to-day reflectance drift to obtain comparable results from the hyperspectral camera. The PRO4SAIL model using normalized hyperspectral reflectance resulted in an average high R² of 0.96 and a low root-mean-squared-difference of 0.02 between the UAV-measured and the Pro4SAIL-modeled hyperspectral reflectance (Supplementary Figure S2).

The land surface radiometric temperature (LST) map (Figure 3) shows warmer bare soil and cooler green rice canopy during the local noontime flight on 16 November 2018, matching well with the land covers in Figure 2 false-color image. The air temperature was 26.7 °C at the flight time and 99% of the LST in the area was within 29.1 to 42.4 °C. The LST histogram of the whole area was bimodal, with two peak temperatures at 31.6 and 36.5 °C, corresponding to vegetation and bare soil respectively. The LST histogram of rice plots showed only one peak at 30.9 °C. Since all the parameters in the PRO4SAIL model were set for the rice canopy, the output LAI and subsequently estimations of albedo and emissivity may not fit other vegetation types in the research area. Therefore, to avoid possible large errors, we only analyzed and presented the results of rice plots in the following sections.

Combining the temperature and hyperspectral reflectance from UAV, and the auxiliary data from ground measurements, including weather station data, leaf water content, concentrations of leaf nutrients and soil total nitrogen, and soil water retention curve, we mmapped soil moisture content, soil matric potential, leaf and canopy properties, and energy components over the three experimental groups. The average changes of key variables after biochar application are summarized in

Table 1. Average changes (%) in key variables of biochar groups compared with the control (marginal mean ± SE) derived from UAV-derived data. Bold fonts indicate the variations significant at p < 0.05.

Variables		Bamboo Biochar BC1 (%)	Sugarcane Biochar BC2 (%)
Soil variables	Soil moisture content from UAV	17.7 ± 0.1	10.8 ± 0.1
	Soil matric potential from UAV 1	44.8 ± 0.7	−66.9 ± 0.7
	Soil albedo	0.5 ± 0.1	1.4 ± 0.1
	Soil surface temperature	−0.7 ± 0.2	−1.4 ± 0.2
Leaf and canopy variables	Gross primary productivity	41.9 ± 3.4	17.5 ± 3.4
	Normalized difference vegetation index	10.0 ± 1.4	7.4 ± 1.4
	Canopy chlorophyll content	32.0 ± 3.0	10.1 ± 3.0
	Water use efficiency	40.8 ± 3.5	13.4 ± 3.5
	Leaf area index	6.0 ± 1.3	−4.9 ± 1.3
Land surface energy components	Net radiation	−0.2 ± 0.3	2.3 ± 0.3
	Latent heat flux (evapotranspiration)	1.1 ± 1.0	4.0 ± 1.0
	Evaporative fraction	3.7 ± 1.0	− 0.2 ± 1.0
	Ground heat flux	−1.8 ± 0.6	3.5 ± 0.6
	Sensible heat flux	−0.4 ± 0.7	−0.4 ± 0.7

¹ From the UAV-derived soil moisture content and the in situ sensor-derived soil water retention curve

. Detailed information is given in Supplementary Table S1.

3.2. Variations in Soil Variables after Biochar Application

Figure 4a shows the map of soil moisture content in the three groups BC1, BC2, and C. The UAV-derived average soil moisture content was significantly higher in biochar amended soils BC1 by 17.7 ± 0.1% and BC2 by 10.8 ± 0.1% compared with the control (Table 1). The soil moisture content measured by in-situ sensors also showed overall higher values in biochar groups than in the control, but the spatial variations within the group were large (Figure 4b,c, Supplementary Figure S3) and those changes were not statistically significant (Supplementary Table S1). UAV-derived soil matric potential was significantly higher in BC1 biochar by 44.8 ± 0.7% and lower in BC2 biochar by −66.9 ± 0.7% than the control (Table 1, Figure 5, Supplementary Figure S4). The soil matric potential measured by in-situ sensors did not vary across treatments (Supplementary Table S1) and the range of spatial variation was also large, from −1.4 to −0.2 MPa (Supplementary Figure S4), contrasting to the UAV-derived soil matric potential of about −0.45 to −0.06 MPa. The fitted soil water retention curves of biochar groups shifted towards higher soil water content at a given water potential compared to the control (Supplementary Figure S5).

The average soil surface albedo values were significantly higher in the two biochar groups, and the soil surface temperature decreased by −0.7 ± 0.2% in BC1 and −1.4 ± 0.2% in BC2. Measurements of soil samples revealed non-significant changes in soil total nitrogen concentration (Supplementary Table S1).

3.3. Variations in Rice Leaf and Canopy Variables after Biochar Application

There were significant increases in canopy GPP, NDVI, chlorophyll content, and WUE in the biochar groups compared to the control (Table 1 and Table S1; Figure 6). Opposite changes were observed in LAI, fractional vegetation cover (f_c), and canopy LST in two biochar groups. The LAI and f_c were significantly higher in BC1 but lower in BC2 than in the control. The canopy LST significantly decreased by −0.6 ± 0.2% in BC1 and increased by 0.6 ± 0.2% in BC2, along with a significant decrease in canopy albedo only in BC2 (−3.1 ± 0.5%). Leaf sample analyses on November 19, 2018 revealed non-significant changes in leaf water and leaf P concentration, whereas the leaf K concentration was significantly higher in both biochar groups, and the leaf N concentration was significantly higher only in the BC1 group (Supplementary Table S1) compared with the control.

3.4. Variations in Evapotranspiration and Land Surface Energy Components

The net radiation, latent heat flux (or evapotranspiration), and ground flux were significantly higher in BC2 by 2.3 ± 0.3%, 4.0 ± 1.0%, and 3.5 ± 0.6% respectively than in the control. The ground heat flux significantly decreased by −1.8 ± 0.6% in BC1, and the evaporative fraction was 3.7 ± 1.0% higher in this group than in the control. There were no significant differences in sensible heat flux among treatments (Table 1, Figure 7).

3.5. Comparison of Key Rice Biophysical Variables with Soil Matric Potential.

Over the experimental plots, the canopy chlorophyll content, NDVI, and LAI were correlated with the soil matric potential with R² of 0.16, 0.12, and 0.32, respectively (Figure 8). However, within-group variations in these parameters were larger than variations across treatments and could not be explained by soil matric potential. The canopy chlorophyll content, NDVI, and LAI had correlations with soil moisture content with relatively lower R² of 0.12, 0.09, and 0.12 respectively (Supplementary Figure S6).

3.6. Comparison of Evaporative Fraction (EF) with Soil Moisture Content and Matric Potential

EF had strong correlations with soil moisture content (θ) in individual groups (R² = 0.87, 0.78, and 0.82 for group C, BC1, and BC2 respectively, Figure 9). The ET-θ fitting lines of two biochar groups BC1 and BC2 shifted toward higher θ compared with the control group C. The slopes of the fitting lines were 0.29, 0.15, and 0.30 for group C, BC1, and BC2 respectively. The evaporative fraction had strong exponential relationships with soil matric potential (ψ) in individual groups (Figure 9b), with the EF-ψ fitting line shifting toward higher ψ after BC1 biochar application and toward lower ψ after BC2 biochar application.

4. Discussion

This study combined UAV-based hyperspectral and thermal imagery with ground measurements and physical models to provide a comprehensive assessment of relative changes in upland rice growth, soil water availability, and WUE after biochar application. While our campaign was relatively short, it allowed monitoring the integrated effects of biochar during most of the rice growing season, from sowing to early senescence. During the four-day campaign, we found that in the two biochar groups, rice GPP, NDVI, chlorophyll content, soil moisture content, and WUE were significantly higher than the control group. The increments were larger in the BC1 group that also had higher leaf nitrogen concentration. In contrast, we found opposite behaviors between the two biochar treatments for soil matric potential, leaf area index, and fractional vegetation cover, which were higher in BC1, but lower in BC2, compared to the control (Table 1). Looking at surface energy fluxes, the net radiation, latent heat flux (evapotranspiration), and ground heat flux significantly increased in BC2, but not in BC1, which translated to a significant increase in the evaporative fraction in BC1, but no change in BC2 compared to the control plots.

4.1. Biochar Effects on Bare Soil Albedo

From the exposed soil surfaces between plots, we did not observe lower soil surface albedo after the biochar application at a rate of 0.4% (weight/weight) within 20 cm topsoil in our experiment. The biochar had been added to the soil for eleven months in BC1 and five months in BC2, so the charcoal particles were likely already incorporated into the soil by the time the measurement campaign was conducted, thereby limiting any direct soil albedo effect from biochar. Our albedo results from the field biochar experiment contrast sharply with other studies that claimed that biochar addition could decrease soil surface albedo and increase radiative forcing feedback to climate based on dark surface albedo either from satellite measurements at largescale or laboratory experiments (e.g., [55,56]).

4.2. Biochar Effects on Soil Water Availability

Biochar addition has been shown to alter soil moisture content, and soil physical and hydraulic properties [20,21]. We found from in situ sensor measurements (one soil moisture probe in each plot), that the soil moisture content (θ) was higher in both biochar plots but differences across treatments were not significant due to the large within-group variability. Using the UAV-derived θ that covered the whole area of each plot to better account for the variability within each plot, we found significant increases in θ in biochar-amended plots (Figure 4). This demonstrated the advantage of high-resolution drone-based remote sensing approach over ground point measurements in cropland and agricultural assessment.

The significantly higher θ in biochar groups than the control was probably due to the increase in water holding capacity in biochar amended soil. The higher water holding capacity is likely linked to the porous biochar structure that increases the particle surface area, so that more water can be retained by the soil. The increased water retention with biochar addition can also reduce nutrient leaching [57] and result in more efficient plant nutrient uptake [24]. This might partly explain the higher leaf nitrogen concentration observed in BC1.

It should be noted that the θ derived from temperature-vegetation relationships reflects the water availability in the whole root zone [43], and the rice rooting depth can reach 0.5–1.0 m [45], far deeper than the depth of 15 cm at which the in-situ sensors were placed in this study. The shallower the measurement location, the stronger the influence from atmospheric forcing and consequently the larger the fluctuations of the measurements, explaining the higher variability of the in-situ measurements compared to the remote sensing estimates. Besides, the groundwater level was reported 0.8 m below the surface during the dry period [29] and the capillary rise in clay loam soil of the site might dampen the fluctuations of root zone soil moisture. The mismatch of in-situ sensor-measured and UAV-derived θ makes the two sets of data not readily comparable.

The soil matric potential presented a different picture for the two biochar applications, increasing in BC1 and decreasing in BC2, contrary to the θ that increased in both biochar applications. The lower soil matric potential found in BC2 despite higher soil moisture is due to the different water retention properties in BC1 and BC2. In fact, at given soil moisture, the matric potential in BC1 is lower than in BC2 (Supplementary Figure S5). Along with the decreased soil matric potential in BC2, the BC2 plots also exhibited lower leaf nitrogen concentration, LAI, GPP, and WUE. Apart from the biochar feedstock difference between BC1 and BC2, [29] attributed the lower plant performance to the shorter time of BC2 settlement in soil. The BC1 had been applied six months earlier and might be better established in soil.

4.3. Biochar Effects on Rice Growth Indicators

We found that rice GPP, NDVI, and CCC were higher in both biochar applications than in the control group (Figure 6). A similar positive effect of biochar addition on maize was reported by Agegnehu, et al. [58] under Australian tropical dry conditions. We indirectly inferred the changes in plant growth from the changes in NDVI, a proxy of the fraction of photosynthetic energy absorbed by the canopy [35,36], and the changes in canopy chlorophyll content that determines the plant photosynthetic capacity and therefore maximum potential light use efficiency. The soil background biases on NDVI values were likely minimal due to the dark-colored soil and relatively dense rice canopy [59] in our study. Ultimately, all these positive changes in canopy properties might be traced back to the aforementioned improvement in soil water availability and subsequently the resilience of plant to dry spells after biochar application in the arid cropland, which was also supported by the results of isotopic composition of plant water [29].

A further correlation analysis showed that soil matric potential was a slightly better predictor of variations in CCC, NDVI, and LAI than the θ (Figure 8 and Figure S5). Vascular plants absorb soil water and transport it through the xylem to the leaf surface thanks to the gradients of water potential along the soil-plant-atmosphere continuum SPAC [60,61]. Assuming a steady-state flow along the SPAC, higher water availability as indicated by higher soil water potential will result in more water transported to the leaves, which guarantees higher leaf water potential and gas exchanges [60]. In turn, this improved water status is expected to allow for larger leaf area and photosynthetic capacity, thereby explaining the observed correlations.

4.4. Biochar Effects on Evaporative Fraction and Plant WUE

We found that the rate of change of the evaporative fraction (EF; i.e., the fraction of net radiation energy converted to latent heat) with soil moisture content in BC1 was about half of that in the other groups, indicating that the rice in BC1 was less sensitive to soil moisture changes at a relatively higher moisture level. Extrapolating the EF-θ fitting curve to the line of EF = EF_max defines hypothetical soil moisture contents at which the available energy is maximally converted to latent heat, indicating the maximum transpiration and minimum sensible heat of the net radiation energy partition. These soil moisture values are higher in the two biochar amendment groups compared to the control (Figure 9), suggesting that the stomata of rice growing in biochar amended soil started closing (thus reducing the evaporative fraction) at higher soil water content than in the control. Moreover, the evaporative fraction remained lower in biochar amended soil at any given soil water level compared to the control, all the way to the wilting point. This finding provides empirical support to the model results in ([21] Figure 4B). This pattern can also be interpreted in a different way—that biochar-amended soils hold more water, but trigger stomatal closure and reduced evapotranspiration at relatively higher soil moisture compared to control plots. As a consequence, despite holding more water, the amount of plant-available water in biochar-amended soils was comparable to control plots.

We found that EF increased nonlinearly with soil matric potential (ψ) in each group (Figure 9b) and may be fitted using sigmoidal curves [62]. The EF-ψ fitting curve of BC2 group shifted toward lower ψ compared to the control, suggesting a higher evaporative fraction at a given ψ in BC2. Indeed, BC2 had the highest evapotranspiration among all three experimental groups (Supplementary Table S1). This curvilinear relation between EF and ψ is a direct consequence of the nonlinear shape of the soil water retention curves [50].

However, the average EF in BC2 was slightly lower than the other two groups during the UAV campaign period, indicating that the rice in this group was relatively water-stressed. Correspondingly, the rice canopy surface temperature, an indicator for instantaneous plant water stress [63], was slightly higher in BC2 than the other two groups (Supplementary Table S1). Moreover, the leaf water content, an indicator for sustained water stress [14], was the lowest in BC2 among the three groups. Although these differences between the biochar groups and the control were generally small in magnitude during the campaign days, if the mild plant stress had been accumulated over time, it could lead to significantly inferior rice growth, explaining the lower LAI and canopy cover in BC2 compared with the control.

How could the biochar amendment affect plant WUE? The answer to the question is very critical to crops in water-limited areas where the experiment locates. In the review by Fischer et al. [21], few studies reported a positive effect of biochar on crop WUE based on yield. Plant WUE is expected to increase if the biochar promotes plant productivity without a significant increase in evapotranspiration. This is the case for BC1, where the gross primary productivity (GPP) increased significantly and the evapotranspiration did not change. If the GPP increases at a significant cost of evaporative water loss, the plant WUE will change depending on which (GPP or evapotranspiration) increases more. This is the case for BC2, where the GPP increased along with the significant increase in evapotranspiration. The GPP increment largely surpassed the increment of evapotranspiration and the consequent WUE was increased in BC2 as well. This result is key as it provides observational evidence on the influence of biochar for WUE, suggesting the amendment would potentially benefit climate-smart agriculture by adding resilience to crops in regions with limited water resources [21,22].

5. Conclusions

This study demonstrated that the joint use of hyperspectral and thermal imagery from a UAV is promising to provide a wide range of indicators for quantifying upland rice productivity and water usage after biochar application. We showed that the indicators derived from the imagery using physical models were useful in exploring biophysical information behind indirect measurements, and in explaining the mechanisms from the observed changes after biochar application. The integration of direct and indirect measurements allows for the identification of relative changes in plant biophysical variables, soil water availability, and surface energy fluxes. Drone-based remote sensing provided spatially explicit information without disturbing the measurement target, densifying ground point measurements from sensors or laboratory samples, enabling a wide range of applications that can benefit precision agriculture and crop management in general.

The results show relatively large increases in GPP, chlorophyll content, and soil moisture content after the biochar amendments. In addition, non-significant or significant but small changes in energy fluxes and evapotranspiration were revealed, suggesting increases in WUE after biochar application. However, the biochar feedstock and application management might have influenced the outcomes as in our case the increment in bamboo biochar amended soil was larger than in sugarcane biochar amended soil.

This study motivates a larger scale experiment so that other aspects related to the impact of the amendment for WUE, GPP, and soil fertility can be tested at a full crop scale. The results demonstrate the potential of biochar for agricultural water management to withstand dry spells by allowing rice plants to efficiently use water and nutrients use. A noteworthy result is that potential climate feedbacks derived from high temperatures and dark soils might not be as significant as previously thought, based on the detected changes in the surface energy balance components affecting radiative forcing such as the land surface temperature, albedos, or sensible heat fluxes after limited biochar application in our experiment. Biochar amendment is therefore promising as part of the strategies to improve the resilience of crops under rainfall reduction scenarios associated with the impacts of climate variability and change. The implementation of biochar to improve WUE, soil quality, and likely carbon uptake can also represent an innovative solution to the treatment of agricultural residues and foster good agricultural practices.