Mapping Root-Zone Soil Moisture Using a Temperature-Vegetation Triangle Approach with an Unmanned Aerial System

: High resolution root-zone soil moisture (SM) maps are important for understanding the spatial variability of water availability in agriculture, ecosystems research and water resources management. Unmanned Aerial Systems (UAS) can ﬂexibly monitor land surfaces with thermal and optical imagery at very high spatial resolution (meter level, VHR) for most weather conditions. We modiﬁed the temperature–vegetation triangle approach to transfer it from satellite to UAS remote sensing. To consider the effects of the limited coverage of UAS mapping, theoretical dry/wet edges were introduced. The new method was tested on a bioenergy willow short rotation coppice site during growing seasons of 2016 and 2017. We demonstrated that by incorporating surface roughness parameters from the structure-from-motion in the interpretation of the measured land surface-atmosphere temperature gradients, the estimates of SM signiﬁcantly improved. The correlation coefﬁcient between estimated and measured SM increased from not signiﬁcant to 0.69 and the root mean square deviation decreased from 0.045 m 3 · m − 3 to 0.025 m 3 · m − 3 when considering temporal dynamics of surface roughness in the approach. The estimated SM correlated better with in-situ root-zone SM (15–30 cm) than with surface SM (0–5 cm) which is an important advantage over alternative remote sensing methods to estimate SM. The optimal spatial resolution of the triangle approach was found to be around 1.5 m, i.e. similar to the length scale of tree-crowns. This study highlights the importance of considering the 3-D ﬁne scale canopy structure, when addressing the links between surface temperature and SM patterns via surface energy balances. Our methodology can be applied to operationally monitor VHR root-zone SM from UAS in agricultural and natural ecosystems.


Study Site
The study site is an 11-hectare short rotation coppice (SRC) willow bioenergy plantation adjacent to the Risoe national laboratory, Roskilde, Denmark (55°41′31.95″N, 12°6′14.69″E). The two willow clones of this field are Salix schwerinii × S. viminalis × S. vim. and Salix triandra × S. viminalis and are planted in rows with 1.5 and 0.75 m distance. The site has a temperate maritime climate with a mean annual temperature of 8.5 °C and precipitation about 600 mm. The soil texture is loam. In February 2016, the willow was harvested. Afterwards, the willow grew rapidly to the height of approximately 4 m during the growing seasons of 2016 and 2017. Rapeseed (Brassica napus) was grown in the nearby field. An approximate 3 m wide path between the willow plantation and the rapeseed field was covered by grass, as shown in Figure 1.
An eddy covariance (EC) observation system (DK-RCW) has been operated from 2012 until now to continuously measure Ta, air pressure (Pa), relative humidity (RH), wind speed (WS), and land surface radiation components at a height of 10 m. Components of the land surface energy balance including incoming shortwave radiation SWin, outgoing shortwave radiation SWout, incoming longwave radiation LWin and outgoing longwave radiation LWout were measured by a CNR4 net radiometer (Kipp & Zonen, Delft, The Netherlands) on the flux tower. Photosynthetically Active Radiation (PAR) measurements were obtained from ten PAR sensors (Apogee SQ-200, Apogee Instruments Inc., Logan, UT, USA) including one sensor to measure the incident PAR above the canopy (PAR ), one sensor to measure canopy-reflected PAR (PAR ) and eight sensors to measure understory PAR (PAR ). For in-situ SM measurements, fixed dielectric probes (5TM ECH2O probes, Decagon Inc., Pullman, WA, USA) were installed to continuously measure SM in two soil profiles, as shown in Figure 1. Profile A has measurements at depths of 5, 15, 30 and 60 cm and Profile B has measurements at depths of 5, 15 and 30 cm. Moreover, a portable TDR (Field Scout TDR 300 portable moisture meter, Spectrum Technologies Inc., Plainfield, IL, USA) was used to measure SM at the layer of 0-10 cm across the willow field on 26 May 2017 (10 samples) and 18

Unmanned Aerial System (UAS) and Flight Campaigns
For the unmanned aerial vehicle, this study used a DJI Hexacopter Spreading Wings S900 (DJI S900, DJI Inc., Shenzhen, China), as shown in Figure 2. DJI S900 can carry an approximate payload weight of 1.5 kg with a flight duration of 15 min. The payload onboard included a Global Navigation Satellite System (GNSS) rover station, a Single-Board Computer (SBC) Beaglebone Black for sensor communication and data storage, and an imaging system with three types of cameras. A GNSS base station was installed close to the study area to serve as the reference for the differential carrier-phase GNSS system, with the rover antenna located on the drone. The base station was a NovAtel receiver (flexpak6) with NovAtel GPS-703-GGG pinwheel triple frequency GPS and GLONASS antenna. The accurate position of the base station was measured by a real-time kinematic (RTK) GNSS (Trimble RTK GNSS R8s, Trimble Inc., Sunnyvale, CA, USA). To compute the absolute antenna position, a carrier-phase differential based solution was computed in post-processing. Raw pseudo ranges and carrier phase measurements were stored at 1 Hz. The position solution was post-processed using Leica Geomatic Office v 8.1 in kinematic mode. In the post process, an ensemble Kalman filter was applied in both forward and backward directions for best position performance. For details, refer to [37]. Further, after each flight campaign, high accuracy ground control points (GCPs) were measured with Trimble RTK GNSS.
The imaging payload included a thermal infrared camera, a multispectral camera, and a normal Red-Green-Blue (RGB) channel camera. Due to the limitation of the UAS payload capacity, a light weight thermal camera (FLIR Tau2 324, Wilsonville, OR, USA), which has an uncooled VOx microbolometer, was used in this study. It has a focal length of 9 mm with an image dimension of 324 × 256 pixels and a field of view (FOV) of 48.5 • × 39.1 • . It records thermal radiation in the wavelength range of 7.5 to 13.5 µm and is able to measure temperatures ranging from −25 • C to 135 • C with high gain mode. The multispectral camera (MCA, Multispectral Camera Array, Tetracam, Chatsworth, CA, USA) consists of an array of six individual channels for the visible and near-infrared bands, each consisting of a CMOS sensor with a progressive shutter, an objective lens, and mountings for interchangeable band-pass filters. Each channel has a FOV of 38.3 • × 31 • and the focal length of 9.6 mm. The center wavelengths for the camera's six channels are 470, 530, 570, 670, 710 and 800 nm and the full width at half maximum (FWHM) for each channel is 10 nm. The RGB camera (Sony DSC-RX100, Corporation, Tokyo, Japan) has a focal length of 10.7 mm with FOV of 64.8 • × 45.9 • . It is mainly used to take RGB images to generate high accuracy digital surface model (DSM) to orthorectify the multispectral and thermal infrared images and also to obtain the canopy height (h c ). UAS flight campaigns were conducted in the willow site over eight days across different growing stages of willow and weather conditions during 2016 and 2017. Images were acquired with 60% side overlap and 80% forward overlap with the horizontal speed of the UAS at 3 m·s −1 . Detailed information on these flight campaigns is shown in Table 1.

Sensor Calibration
To obtain high quality UAS data, thorough laboratory geometric and radiometric calibrations of imaging sensors were firstly conducted. Geometric calibration of the RGB and multispectral cameras was conducted with standard checkerboard calibration patterns to retrieve intrinsic camera geometric parameters. Radiometric calibration for each channel of MCA was conducted with an integrating sphere (CSTM-USS-2000C, LabSphere, NH, USA) and the in-lab calibration showed that the absolute errors of the measured radiance were within ±4.8%. For details on the calibration of the multispectral camera, please refer to [38].
The accuracy of FLIR is influenced by both the target temperature and the sensor body temperature. To obtain accurate thermal imagery usable for SM monitoring, this study conducted a calibration of FLIR with a Landcal P80P black body radiation source (Land Instruments, Leicester, UK). The calibration was performed with ten different target temperatures ranging from 0 to 45 °C in a climate-controlled chamber, which provided three different ambient temperatures, 10, 20 and 30 °C. Before taking the first snapshot of a series, the thermal camera was powered on for at least 15 min to reach thermal equilibrium. In total, 560 images were taken for the calibration. Based on the uniform target temperature provided by the black body, we conducted a pixel-wise calibration to simultaneously remove the image vignetting effects and link the image digital number (DN) with temperature.
For the FLIR camera, the DN of each pixel in the thermal imagery can be assumed to be proportional to the incoming thermal radiance Bp (W•m −3 •sr −1 ) of the target with additional radiance emitted from the camera itself, as shown in Equation (1).

Sensor Calibration
To obtain high quality UAS data, thorough laboratory geometric and radiometric calibrations of imaging sensors were firstly conducted. Geometric calibration of the RGB and multispectral cameras was conducted with standard checkerboard calibration patterns to retrieve intrinsic camera geometric parameters. Radiometric calibration for each channel of MCA was conducted with an integrating sphere (CSTM-USS-2000C, LabSphere, NH, USA) and the in-lab calibration showed that the absolute errors of the measured radiance were within ±4.8%. For details on the calibration of the multispectral camera, please refer to [38].
The accuracy of FLIR is influenced by both the target temperature and the sensor body temperature. To obtain accurate thermal imagery usable for SM monitoring, this study conducted a calibration of FLIR with a Landcal P80P black body radiation source (Land Instruments, Leicester, UK). The calibration was performed with ten different target temperatures ranging from 0 to 45 • C in a climate-controlled chamber, which provided three different ambient temperatures, 10, 20 and 30 • C. Before taking the first snapshot of a series, the thermal camera was powered on for at least 15 min to reach thermal equilibrium. In total, 560 images were taken for the calibration. Based on the uniform target temperature provided by the black body, we conducted a pixel-wise calibration to simultaneously remove the image vignetting effects and link the image digital number (DN) with temperature.
For the FLIR camera, the DN of each pixel in the thermal imagery can be assumed to be proportional to the incoming thermal radiance B p (W·m −3 ·sr −1 ) of the target with additional radiance emitted from the camera itself, as shown in Equation (1).
where B p is the thermal radiance from the target, which can be calculated based on the Planck's Law (Equation (2)); l (W·m −3 ·sr −1 ) is a gain factor for each pixel; O is a constant offset term; T core is the DN of the sensor core temperature; and t is a factor. Here, we assumed that the additional radiation from the sensor core is approximated as a first-order polynomial. In Planck's Law (Equation (2)), h is the Planck constant (6.626 × 10 −34 J·s), c is the speed of light (2.998 × 10 −8 m·s −1 ), k B is the Boltzmann constant (1.381 × 10 −23 J·K −1 ), and λ (m) is the emitting wavelength. Solving Equations (1) and (2) for B p , we can get the transfer function (Equation (3)) for the target temperature (T b ), based on DN in the thermal imagery and the sensor core temperature (T core ).
where parameters l, O and t were obtained by fitting the calibration dataset. For λ, the sensor reference wavelength of 10.076 µm was used.

Image Processing and Validation
The workflow for image processing to obtain DSM, normalized difference vegetation index (NDVI), T b and T s (land surface temperature after emissivity correction) for this study is shown in Figure 3. The first step of image processing was to geo-reference images with GNSS data from UAS. Then, these georeferenced images along with pre-calibrated intrinsic camera geometric parameter values were imported into Agisoft Photoscan (Agisoft LLC, St. Petersburg, Russia) to generate orthophotos. Agisoft Photoscan software uses the structure-from-motion (SfM) approach [34] to generate orthomosaic and surface elevation maps from overlapping images from different positions and orientations. The Agisoft software firstly aligned images and generated sparse 3D point clouds. After that, high accuracy GCPs measured by the Trimble RTK GNSS were added to the aligned images. Then, the GNSS information from these accurate GCPs was used to optimize the image alignment by adjusting the estimated camera locations. Further, the DSM generated from RGB images was imported into both multispectral and thermal projects to orthorectify multispectral and thermal images. Finally, the six-band reflectance was calculated from the generated multispectral orthophotos based on the radiance method. An ASD spectroradiometer (FieldSpec HandHeld 2 TM Inc., Boulder, CO, USA) was used to collect with a Spectralon panel (nominal reflectance of 99.99%) before and after each UAS flight campaign. Then, the averaged spectral radiance was used as the incoming solar spectral radiance (L in,λ ), while MCA measured the reflected spectral radiance for each band (L MCA,λ ). Due to the low flight altitude, we did not conduct the atmospheric correction. Six-band reflectance was calculated as Equation (4). The reflectance of the near infrared (800 nm) and red (670 nm) bands was used to calculate NDVI as Equation (5).
where L MCA,λ is the reflected spectral radiance measured by MCA (W·m −2 ·sr −1 ·nm −1 ), L in,λ is the incoming solar spectral radiance measured by ASD (W·m −2 ·sr −1 ·nm −1 ), λ is the wavelength (nm). ρ 800 is the reflectance at 800 nm, and ρ 670 is the reflectance at 670 nm. The observed fraction of vegetation cover ( f c ) as Equation (6) was calculated based on NDVI [39]. Further, the observed f c (Equation (7)) was used to validate the estimated f c from UAS ( f c_UAS ).
where NDVI max = 0.97 for fully vegetated surface and NDVI min = 0.24 for bare soil were used in this study based on the field investigation. f c_obs was calculated based on measurements of the incident PAR above the canopy (PAR above ), canopy-reflected PAR (PAR reflected ) and the average of eight understory PAR sensors (PAR below ). The T b from the thermal orthophotos was influenced by the longwave radiation emitted by the object, the reflected longwave radiation and the atmospheric attenuation. UAS flight campaigns were conducted with a low attitude (Table 1), thus, in this study, we did not conduct atmospheric correction. To obtain the thermal signal emitted from the object (T s ), T b was corrected with the land surface emissivity (ε s ) to exclude the signal from reflected longwave radiation as Equation (8). ε s can be approximated based on an empirical relation with NDVI as Equation (9) [40].
where σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ), T s is the surface temperature with emissivity correction (K), T b is the brightness temperature from FLIR (K), and LW in is incoming longwave radiation (W·m −2 ) and was calculated based on the Stefan-Boltzmann law using the atmospheric emissivity ε a as Equations (10)-(13) [41].
where e 0 is the actual water vapor pressure (hPa) and was calculated based on the Clausius-Clapeyron equation, L v = 2.5 × 10 6 J·kg −1 is the latent heat of vaporization, R v = 461 J·kg −1 ·K −1 is the gas constant for water vapor, and T a is the air temperature (K). Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 28 To check the accuracy of the generated DSM, the Trimble RTK GNSS was used to measure the ground elevation in different places across the site. It should be noted that the validation was mainly conducted during periods with low vegetation height (hc < 0.5). During the dense vegetation period, we did not conduct validation, due to the difficulty in measuring the tree height.
To validate the acquired surface reflectance data, the reflectance of tarpaulins with four different colors (green, blue, black and silver) was measured with the ASD spectroradiometer for each flight. Tarpaulins are acknowledged to have low anisotropic effects [42]. The measured reflectance by ASD was used to validate reflectance from MCA. Validation was shown for 670 and 800 nm bands, since, in this study, only these two bands were used for calculating NDVI.
To conduct validation of the thermal infrared imagery in field conditions, the generated Tb orthophotos were compared with the brightness temperature converted from the outgoing longwave radiation (LW _ ) from CNR4 on the tower ( Figure 1). Even though FLIR and CNR4 have different thermal wavelength responses (FLIR: 7.5 to 13.5µm, CNR4: 4.5 to 45 µm) and sensor FOVs (FLIR: 35° × 27°, CNR4: downward 150°), this comparison can provide insights into the quality of the thermal data to some extent. To compare with Tb from FLIR, LW _ was converted to brightness temperature Tb_CNR4 by applying Stefan-Boltzmann's law as Equation (14).
where σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W•m −2 •K −4 ), T _ is the brightness temperature (K), and LW _ is the recoded longwave radiation by CNR4 (W•m −2 ). To get the corresponding pixels inside the FOV of CNR4, the flux tower position, the FOV and sensor height of CNR4, and the vegetation height were used to calculate the radius of the circle To check the accuracy of the generated DSM, the Trimble RTK GNSS was used to measure the ground elevation in different places across the site. It should be noted that the validation was mainly conducted during periods with low vegetation height (h c < 0.5). During the dense vegetation period, we did not conduct validation, due to the difficulty in measuring the tree height.
To validate the acquired surface reflectance data, the reflectance of tarpaulins with four different colors (green, blue, black and silver) was measured with the ASD spectroradiometer for each flight. Tarpaulins are acknowledged to have low anisotropic effects [42]. The measured reflectance by ASD was used to validate reflectance from MCA. Validation was shown for 670 and 800 nm bands, since, in this study, only these two bands were used for calculating NDVI.
To conduct validation of the thermal infrared imagery in field conditions, the generated T b orthophotos were compared with the brightness temperature converted from the outgoing longwave radiation (LW out_CNR4 ) from CNR4 on the tower ( Figure 1). Even though FLIR and CNR4 have different thermal wavelength responses (FLIR: 7.5 to 13.5 µm, CNR4: 4.5 to 45 µm) and sensor FOVs (FLIR: 35 • × 27 • , CNR4: downward 150 • ), this comparison can provide insights into the quality of the thermal data to some extent. To compare with T b from FLIR, LW out_CNR4 was converted to brightness temperature T b_CNR4 by applying Stefan-Boltzmann's law as Equation (14).
where σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ), T b_CNR4 is the brightness temperature (K), and LW out_CNR4 is the recoded longwave radiation by CNR4 (W·m −2 ). To get the corresponding pixels inside the FOV of CNR4, the flux tower position, the FOV and sensor height of CNR4, and the vegetation height were used to calculate the radius of the circle according to Equation (15), as one example, with the vegetation height equal to 0 m, as shown in Figure 1.
where r is the radius of the projected sensor FOV on the surface (m), H is the sensor height (10 m), h c is the height of the willow plantation obtained from DSM (m), and FOV is the downward FOV of CNR4, namely 150 • .

The Original "DT" Triangle Approach
Numerous satellite studies have applied the temperature-vegetation triangle approach to estimate both SM and evaporative fraction [11,18,21,22,28,29]. However, the linkage between applying the triangle approach for SM and evapotranspiration estimation is seldom explained. For SM estimation, originally, the temperature vegetation triangle index or soil wetness index (SWI) was defined as Equation (16), by rescaling the observed radiometric surface temperature (T s,obs ) at each pixel between extreme dry and wet T s for a given level of f c [11,21,28] (hereinafter, "original triangle approach"). This approach assumed that the only factor modifying T s in the triangle space apart from f c is SM. The SWI is used to indicate SM conditions.
where the subscripts "dry", "wet" and "obs" indicate dry, wet and observed conditions, respectively. T s,dry and T s,wet are the minimum and maximum T s for a given f c . For the wet edge, T s,wet corresponds to the temperature of the wet land surface and numerous studies proved that T s,wet is in equilibrium with T a [22,36,43]. That means T s,wet can be replaced by T a and DT wet is equal to 0 [29]. DT dry can normally be extracted from the observed dry edge in the triangle space. However, if there are no sufficient pixels to represent dry and wet conditions, DT dry can be approximated from a linear mixture of DT bs,dry (temperature difference between the dry bare soil temperature and the air temperature) and DT wet as Equation (17) [18,44]. The triangle approach can be quantitatively explained by the surface energy balance and one source evapotranspiration model [18,29], as shown in Figure 4. As shown in Figure 4b, the relationship between SM and Non-Evaporative Fraction (NEF = H/AE) for any given pixel within a fixed f c can be described in a linear model approach Equation (18), due to the land-atmosphere coupling [1]. In Figure 4a, for a fixed level of f c , the increases in T s are due to increases in sensible heat flux (H), which are determined by decreases in SM.
where AE is the available energy and is equal to the sum of the sensible heat flux (H) and the latent heat flux (LE), θ is the volumetric SM (m 3 ·m −3 ), θ fc is the field capacity (m 3 ·m −3 ), and θ wt is the wilting point of soil (m 3 ·m −3 ), representing the upper and lower limits of available SM [11]. (a) Conceptual sketch of the temperature-vegetation triangle approach. The x-axis is the fractional vegetation cover (fc). The y-axis is the sensible heat flux (H) and it can be simplified to be proportional to DT by canceling air density (ρ), the specific heat of air at constant pressure (cp) and aerodynamic resistance for heat transfer (ra). DT represents the difference between air and surface temperatures. DTdry and DTwet are for the dry and wet edges, while DTobs is the observed temperature difference for a given pixel and DTbs_dry is the temperature difference for the dry bare soil. ra,wet, ra,dry, ra,obs and ra,bs,dry are aerodynamic resistance for heat transfer for wet, dry, observed and driest bare soil conditions, respectively. The upper red line and lower blue line represent the dry and wet edges, respectively. (b) Diagram (adapted from [1]) representing the relationship between non-evaporative fraction (NEF) and SM. In the dry conditions (θ < θ ), the available energy (AE) converts into sensible heat. In the wet conditions (θ ≥ θ ), the available energy (AE) converts into latent heat flux. For transitional conditions, the NEF has a linear relationship with SM. θ is the wilting point and θ is the field capacity in this study. These two represent the lower and upper limits of available SM, respectively.
For NEF, it can be further expressed as a ratio between H and Hdry as Equation (19). Hdry is the sensible heat flux at the dry edge, where LE is assumed to be negligible (Equation (20)). With the bulk sensible heat flux equation, H can be expressed as a ratio between DT and aerodynamic resistance (ra). Further, if we assume that meteorological forcing and aerodynamic resistance (ra) is constant within the area, NEF can further be simplified into a ratio between DTobs and DTdry by canceling the air density (ρ), the specific heat of air at constant pressure (cp) and aerodynamic resistance for heat transfer (ra).
where DTobs is the temperature difference between the surface and air for each pixel, DTdry is the temperature difference for the dry pixel corresponding to a certain fc level (Figure 4a), ρ is the air density (kg•m −3 ), cp is the specific heat of air at constant pressure (J•kg −1 •K −1 ), ra is the aerodynamic . The y-axis is the sensible heat flux (H) and it can be simplified to be proportional to DT by canceling air density (ρ), the specific heat of air at constant pressure (c p ) and aerodynamic resistance for heat transfer (r a ). DT represents the difference between air and surface temperatures. DT dry and DT wet are for the dry and wet edges, while DT obs is the observed temperature difference for a given pixel and DT bs_dry is the temperature difference for the dry bare soil. r a,wet , r a,dry , r a,obs and r a,bs,dry are aerodynamic resistance for heat transfer for wet, dry, observed and driest bare soil conditions, respectively. The upper red line and lower blue line represent the dry and wet edges, respectively. (b) Diagram (adapted from [1]) representing the relationship between non-evaporative fraction (NEF) and SM. In the dry conditions (θ < θ wt ), the available energy (AE) converts into sensible heat. In the wet conditions (θ ≥ θ fc ), the available energy (AE) converts into latent heat flux. For transitional conditions, the NEF has a linear relationship with SM. θ wt is the wilting point and θ fc is the field capacity in this study. These two represent the lower and upper limits of available SM, respectively.
For NEF, it can be further expressed as a ratio between H and H dry as Equation (19). H dry is the sensible heat flux at the dry edge, where LE is assumed to be negligible (Equation (20)). With the bulk sensible heat flux equation, H can be expressed as a ratio between DT and aerodynamic resistance (r a ). Further, if we assume that meteorological forcing and aerodynamic resistance (r a ) is constant within the area, NEF can further be simplified into a ratio between DT obs and DT dry by canceling the air density (ρ), the specific heat of air at constant pressure (c p ) and aerodynamic resistance for heat transfer (r a ).
where DT obs is the temperature difference between the surface and air for each pixel, DT dry is the temperature difference for the dry pixel corresponding to a certain f c level (Figure 4a), ρ is the air density (kg·m −3 ), c p is the specific heat of air at constant pressure (J·kg −1 ·K −1 ), r a is the aerodynamic resistance for heat transfer (s·m −1 ), r a_obs is the aerodynamic resistance for the observed pixel, and r a_dry is the aerodynamic resistance for the dry conditions. In this way, we can explain the theories for using the triangle approach for both SM and evapotranspiration fraction (approximately 1-NEF). In this study, SWI (Equation (16)) was converted to the volumetric SM with a piecewise approach [11]. For energy limited evapotranspiration regime (θ ≥ θ fc ), SWI is equal to 0. For water limited, dry condition and no evapotranspiration conditions (θ ≤ θ wt ), SWI is equal to 1. For water limited evapotranspiration regime (transitional conditions in Figure 4b) (θ wt < θ < θ fc ), the SWI can be expressed as Equation (21).
3.4.2. The Modified "DT/r a " Triangle Approach In a real-world condition, the assumption of the triangle approach may not be valid. For instance, if r a is not homogeneous for a given f c and causes variations in DT unrelated with SM, DT should be normalized with r a , especially with fast surface elevation changes. In this study, the willow bioenergy plantation grows fast during the growing season and the change of h c can alter the surface roughness to further affect DT. To exclude the influence of the change of surface roughness, we need to normalize DT with r a for heat transfer. Thus, in this study (Equation (16)), we modified SWI to more accurately represent NEF as Equation (22) instead of Equation (16) To apply a modified triangle approach, we calculated an aerodynamic resistance and DT dry . The aerodynamic resistance for heat transfer (r a ) was calculated based on the algorithm proposed by Brutsaert (1982) as Equation (23) [45].
where z is the height of the reference wind velocity (m), i.e., the sensor height (10 m). To calculate r a,obs in the real conditions, d is the zero-plane displacement height (m) and was chosen equal to 2/3 of h c . z 0m is the surface roughness length for momentum transport and was equal to 0.1 of h c . h c is the effective averaged h c for the willow plantation. h c was obtained by subtracting UAS-based DSM, which was obtained from each UAS flight campaign, with the digital elevation model (DEM). The DEM was obtained from UAS flight campaigns before the willow grew and only contained the elevation information of the ground surface (DSM without trees). z 0h is the roughness lengths for heat transport and was calculated based on Equation (24). kB −1 is a parameter to account for the local boundary resistance for the heat transfer and a constant value of 2.3 was adopted in this study [46]. k is von Karman constant (0.4). u is the horizontal wind speed at the reference height (m·s −1 ). For the dry bare soil conditions, d = 0 and z 0m = 0.005 were used to calculate r a,dry in Equation (23) [47]. The dry edge, represented by DT dry , cannot be extracted from the observed dry edge in the triangle space as UAS imagery does not have sufficient pixels to represent dry and wet conditions. Thus, we used ground meteorological observations to calculate theoretical dry edges. The dry edge DT dry /r a,dry was obtained from a linear mixture of DT bs,dry /r a,bs,dry (extreme dry bare soil conditions) and DT wet /r a,wet (fully vegetated surface conditions). DT wet /r a,wet can be approximated to 0, since the sensible heat flux for fully vegetated conditions is equal to 0 [18,44].
DT dry r a,dry = (1 − f c ) DT bs,dry r a,bs,dry (25) For extremely dry bare soil conditions, LE is equal to 0. DT bs,dry can be derived based on the surface energy balance principle and the heat transfer equation as Equation (26) [44,47].
where SW in is the shortwave incoming radiation (W·m −2 ) and, in this study, it was obtained from the meteorological tower. α is the albedo of the dry bare soil and was adopted as 0.2 in this study according to the field investigation. ε ss is the land surface emissivity for the dry bare soil (NDVI = 0.24) and the value of 0.94 was adopted according to Equation (10). ε a is the emissivity for the atmosphere and was calculated based on Equation (10). σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W·m −2 ·K −4 ). LE was assumed to be 0 for the dry bare soil conditions. c is a ratio between the ground heat flux and net radiation. The value of 0.3 was adopted in this study [48]. ρ is the air density (kg·m −3 ). C p is the specific heat of air at constant pressure (J·kg −1 ·K −1 ).

Sensitivity Test of the Modified Triangle Approach to the Canopy Height
To investigate the sensitivity of the proposed "DT/r a " triangle approach to h c , a sensitivity test to model inputs (h c , DT and f c ) was conducted. The atmospheric forcing was based on 26 May 2017 and there were two experiments in the test. The first experiment was to control f c and quantified the relationship between h c and SWI with different DT. A pixel was assumed to be f c equal to 0.5, h c from 0.1 to 3 m with an interval of 0.01 m, and DT from 0 to 5 • C with an interval of 1 • C. The second one was to control DT, which was assumed to be equal to 1 • C, and explored the relationship between h c and SWI with different levels of f c . A pixel was assumed to have h c from 0.1 to 3 m with a 0.01 m interval and f c from 0 to 1 with a step of 0.1. Totally, 3492 simulations were conducted. Then, these results were used to analyze the sensitivity of the "DT/r a " approach to h c , DT and f c . This sensitivity analysis quantified the sensitivity of model inputs. Further, with the information of input uncertainty from the validation of UAS observations, this analysis could also provide insights into the uncertainty of the simulated SWI.

Validation of Soil Moisture Estimates
The estimated SWI was converted to the volumetric soil moisture (θ) as Equation (22) with θ fc and θ wt equal to 0.31 and 0.15, respectively. The values of θ fc and θ wt were determined according to the soil texture (loam) and long-term in-situ SM observations from 2012 to 2017. Here, we use the "DT" and "DT/r a " to denote the original and modified triangle approaches, respectively. In-situ SM observations from the portable and fixed SM sensors were used to validate the spatial and temporal dynamics of estimated SM from these two approaches. The Root Mean Square Deviation (RMSD, Equation (27)), Correlation Coefficient (R, Equation (28)), Relative Error (RE, Equation (29)), unbiased Root Mean Square Deviation (ubRMSD, Equation (30)) and Standard Deviation (STD, Equation (31)) were used as statistical indices to evaluate the performance of estimated SM. The Taylor diagram (Taylor, 2001), which presents the complementary statistics among R, Normalized STD (NSTD, as Equation (32)) and Normalized ubRMSD (NubRMSD), was used. In the Taylor diagram, these three statistical indices have a triangle-cosine-law-like relationship, as Equation (33). The radial distance of the diagram stands for the NSTD and the angle in the polar plot represents R. The reference point located on the x-axis with R = 1, NSTD = 1, and NubRMSD = 0 is the observation. The distance from the simulation point to the reference point means the NubRMSD of simulations and it is the integrated performance of the simulation.
NubRMSD 2 obs,sim = NSTD 2 obs + NSTD 2 sim − 2NSTD obs NSTD sim cos CC obs,sim (33) where sim is the simulation, obs is the observation, i refers to the i th simulation or observation, N is the total number of data points, sim is the average of the simulation, and obs is the average of the observation. A scale discrepancy may exist between pixel-based UAS imagery and point-based in-situ SM measurements. In the spatial validation, to compare with in-situ measurements, we used the mean values from buffer zones (circles) of different radiuses around the sample location in UAS imagery. In the temporal validation, we compared the estimated SM and in-situ measurements at different soil depths.

UAS Sensor Calibration and Data Validation
Results of pixel-wise calibration in the laboratory are shown in Figure 5. With the pixel-wise calibration, the RMSDs of the thermal imagery were less than 0.57 • C. The nearly uniform pattern indicates that, through the calibration, the typical camera image errors, e.g., vignetting effects, were removed. The orthophotos of T b , NDVI, and DSM obtained from UAS flight campaigns are shown in Figure 6. It can be found that the first three flights of 2016 were obtained at the beginning of the growing season with NDVI < 0.5 and DSM < 15.  To check the quality of these orthophotos, Tb, reflectance and DSM were compared with in-situ observations (Figure 7). Comparison between the UAS Tb and CNR4 Tb showed a close correspondence between the brightness temperatures with R 2 equal to 0.97 and RMSD of 0.93 °C. Figure 7b indicates a good accuracy of the reflectance at 670 and 800 nm with R 2 more than 0.98 and RMSD less than 3%. Figure 7c shows the accuracy of the generated DSM from the RGB images. The RMSD of the DSM is 0.06 m from three campaigns in May 2016 when the willow plants had a small growth height (hc < 0.5 m). Validation of fc (Figure 7d) revealed a good correspondence with the observations with R 2 equal to 0.99 and RMSD of 0.04.  To check the quality of these orthophotos, T b , reflectance and DSM were compared with in-situ observations (Figure 7). Comparison between the UAS T b and CNR4 T b showed a close correspondence between the brightness temperatures with R 2 equal to 0.97 and RMSD of 0.93 • C. Figure 7b indicates a good accuracy of the reflectance at 670 and 800 nm with R 2 more than 0.98 and RMSD less than 3%. Figure   To check the quality of these orthophotos, Tb, reflectance and DSM were compared with in-situ observations (Figure 7). Comparison between the UAS Tb and CNR4 Tb showed a close correspondence between the brightness temperatures with R 2 equal to 0.97 and RMSD of 0.93 °C. Figure 7b indicates a good accuracy of the reflectance at 670 and 800 nm with R 2 more than 0.98 and RMSD less than 3%. Figure 7c shows the accuracy of the generated DSM from the RGB images. The RMSD of the DSM is 0.06 m from three campaigns in May 2016 when the willow plants had a small growth height (hc < 0.5 m). Validation of fc (Figure 7d) revealed a good correspondence with the observations with R 2 equal to 0.99 and RMSD of 0.04.   Figure 8 shows the modified triangle approach is sensitive to the change of hc, DT and fc. In general, in the same conditions of the observed DT and fc, the higher hc corresponds to higher SWI. Figure 8a shows that, with the same level of SWI and fc, the higher DT corresponds to the higher hc. This agrees with the fact that, for pixels with the same SM, fc and atmospheric conditions, pixels with a higher hc tend to have higher surface roughness and lower ra. Similarly, as shown in Figure 8b, with the same levels of fc and DT, the pixels with higher hc correspond to higher SWI. Different from the sensitivity of SWI to DT with the same level of hc, which is approximately linear, the sensitivity of SWI to fc with the same hc is highly non-linear. In low fc conditions, the sensitivity of SWI to fc is low, while it is much higher in the high fc conditions. Within the same levels of DT and fc, the sensitivity of SWI to hc is closer to linear.

Sensitivity of the Modified Triangle Approach to the Vegetation Height
Further, Figure 8 also provides insights into simulation uncertainty of this study. As shown in the validation of UAS orthophotos (Figure 7), with the pixel-wise calibration approach, UAS-based Ts can have an uncertainty around 0.93 °C in the field. The RMSD of reflectance is around 3% and the RMSD of hc is about 0.06 m. It should be noticed that RMSD of hc is assessed when hc < 0.5 m and errors for high hc conditions is likely to be higher. Here, we use 0.5 m to consider the uncertainty of hc. In Figure 8a, for instance, with the condition DT = 1 °C and fc = 0.5, the uncertainty of 0.5 m in hc can contribute to the uncertainty of SWI around 0.02. With the condition fc = 0.5 and hc = 2 m, the 1 °C uncertainty of Ts can contribute to the uncertainty of SWI about 0.15. As for fc in Figure 8b, with the condition DT = 1 °C and hc = 1 m, it can contribute to the uncertainty of SWI around 0.05 in low fc conditions (fc < 0.5), while the uncertainty of SWI is much higher in high fc conditions. Therefore, in this approach, the accuracy of SWI is highly sensitive to DT and, in high fc conditions, it is also highly sensitive to fc.  Figure 8 shows the modified triangle approach is sensitive to the change of h c , DT and f c . In general, in the same conditions of the observed DT and f c , the higher h c corresponds to higher SWI. Figure 8a shows that, with the same level of SWI and f c , the higher DT corresponds to the higher h c . This agrees with the fact that, for pixels with the same SM, f c and atmospheric conditions, pixels with a higher h c tend to have higher surface roughness and lower r a . Similarly, as shown in Figure 8b, with the same levels of f c and DT, the pixels with higher h c correspond to higher SWI. Different from the sensitivity of SWI to DT with the same level of h c , which is approximately linear, the sensitivity of SWI to f c with the same h c is highly non-linear. In low f c conditions, the sensitivity of SWI to f c is low, while it is much higher in the high f c conditions. Within the same levels of DT and f c , the sensitivity of SWI to h c is closer to linear.

Sensitivity of the Modified Triangle Approach to the Vegetation Height
Further, Figure 8 also provides insights into simulation uncertainty of this study. As shown in the validation of UAS orthophotos (Figure 7), with the pixel-wise calibration approach, UAS-based T s can have an uncertainty around 0.93 • C in the field. The RMSD of reflectance is around 3% and the RMSD of h c is about 0.06 m. It should be noticed that RMSD of h c is assessed when h c < 0.5 m and errors for high h c conditions is likely to be higher. Here, we use 0.5 m to consider the uncertainty of h c . In Figure 8a, for instance, with the condition DT = 1 • C and f c = 0.5, the uncertainty of 0.5 m in h c can contribute to the uncertainty of SWI around 0.02. With the condition f c = 0.5 and h c = 2 m, the 1 • C uncertainty of T s can contribute to the uncertainty of SWI about 0. 15. As for f c in Figure 8b, with the condition DT = 1 • C and h c = 1 m, it can contribute to the uncertainty of SWI around 0.05 in low f c conditions (f c < 0.5), while the uncertainty of SWI is much higher in high f c conditions. Therefore, in this approach, the accuracy of SWI is highly sensitive to DT and, in high f c conditions, it is also highly sensitive to f c .

Spatial Validation of UAS Estimated Soil Moisture
The spatial patterns of the estimated SWI from the two triangle approaches ("DT" and "DT/ra") are shown in Figure 9(a1-a8) and (b1-b8), respectively. It can be seen that before the willow grew (the first three flights of 2016), there was no difference in spatial patterns and histogram distributions between "DT" and "DT/ra" as shown in Figure 6(a1-a3) and (b1-b3). However, when the willow grew up (the last five flights) and the canopy grew higher, the spatial patterns of SWI from these two approaches became different. In the "DT/ra" approach, the normalization of ra tends to stretch SWI of willow plantation to a higher value. From the histogram distributions (Figure 9(c4-c8)), it can also be seen that, after normalization, more pixels have higher values of SWI.

Spatial Validation of UAS Estimated Soil Moisture
The spatial patterns of the estimated SWI from the two triangle approaches ("DT" and "DT/r a ") are shown in Figure 9(a1-a8) and (b1-b8), respectively. It can be seen that before the willow grew (the first three flights of 2016), there was no difference in spatial patterns and histogram distributions between "DT" and "DT/r a " as shown in Figure 6(a1-a3) and (b1-b3). However, when the willow grew up (the last five flights) and the canopy grew higher, the spatial patterns of SWI from these two approaches became different. In the "DT/r a " approach, the normalization of r a tends to stretch SWI of willow plantation to a higher value. From the histogram distributions (Figure 9(c4-c8)), it can also be seen that, after normalization, more pixels have higher values of SWI.

Spatial Validation of UAS Estimated Soil Moisture
The spatial patterns of the estimated SWI from the two triangle approaches ("DT" and "DT/ra") are shown in Figure 9(a1-a8) and (b1-b8), respectively. It can be seen that before the willow grew (the first three flights of 2016), there was no difference in spatial patterns and histogram distributions between "DT" and "DT/ra" as shown in Figure 6(a1-a3) and (b1-b3). However, when the willow grew up (the last five flights) and the canopy grew higher, the spatial patterns of SWI from these two approaches became different. In the "DT/ra" approach, the normalization of ra tends to stretch SWI of willow plantation to a higher value. From the histogram distributions (Figure 9(c4-c8)), it can also be seen that, after normalization, more pixels have higher values of SWI. Figure 9. (a,b) Spatial patterns of the estimated SWI. "DT" (the first and third columns) represents the scheme without normalizing surface roughness and "DT/ra" (the second and fourth columns) Figure 9. (a,b) Spatial patterns of the estimated SWI. "DT" (the first and third columns) represents the scheme without normalizing surface roughness and "DT/r a " (the second and fourth columns) indicates the scheme with normalizing aerodynamic resistance. (c) Histograms of SWI. The blue color is the "DT" scheme and the red color is the "DT/r a " scheme. The magenta color is the overlapping of the histograms from these two schemes.
The corresponding T s -f c scatterplots for these UAS flights are shown in Figure 10. From the scatterplots, it is clear that there is a large mismatch between the observed dry and wet edges and the theoretical ones, since UAS campaigns do not have a large coverage of heterogeneous landscapes to represent the dry and wet edges. This confirms that it is necessary to calculate the theoretical dry and wet edges for the triangle approach. It should also be noticed that, in the scatterplots, there are outliers in the feature space enclosed by the theoretical dry and wet edges. This can be attributed to the difference between aerodynamic and radiometric temperatures [49] or errors in the retrieval of T s from the UAS sensors. Due to the lower r a of the willow plantation than the pixels of grass or bare soil, after normalization, points in the "DT/r a " approach shifted toward the dry edge as compared to the "DT" approach ( Figure 10). This also agrees with the change of the spatial pattern and the histogram distribution in Figure 9.
indicates the scheme with normalizing aerodynamic resistance. (c) Histograms of SWI. The blue color is the "DT" scheme and the red color is the "DT/ra" scheme. The magenta color is the overlapping of the histograms from these two schemes.
The corresponding Ts-fc scatterplots for these UAS flights are shown in Figure 10. From the scatterplots, it is clear that there is a large mismatch between the observed dry and wet edges and the theoretical ones, since UAS campaigns do not have a large coverage of heterogeneous landscapes to represent the dry and wet edges. This confirms that it is necessary to calculate the theoretical dry and wet edges for the triangle approach. It should also be noticed that, in the scatterplots, there are outliers in the feature space enclosed by the theoretical dry and wet edges. This can be attributed to the difference between aerodynamic and radiometric temperatures [49] or errors in the retrieval of Ts from the UAS sensors. Due to the lower ra of the willow plantation than the pixels of grass or bare soil, after normalization, points in the "DT/ra" approach shifted toward the dry edge as compared to the "DT" approach ( Figure 10). This also agrees with the change of the spatial pattern and the histogram distribution in Figure 9.  Figure 11 shows the spatial validation of "DT" and "DT/ra" approaches with different buffer zone sizes. Because the hc was calculated based on the averaged hc, the normalization did not change the correlation coefficient R. However, it can also be seen that there was a significant improvement in RMSD and relative errors (RE) with normalized ra than without. With a buffer zone of 1.5 m radius, as shown in Figure 11g, on 26 May 2017, the RE reduced from 14.02% to −1.51%. Similarly, on 18 June  Figure 11 shows the spatial validation of "DT" and "DT/r a " approaches with different buffer zone sizes. Because the h c was calculated based on the averaged h c , the normalization did not change the correlation coefficient R. However, it can also be seen that there was a significant improvement in RMSD and relative errors (RE) with normalized r a than without. With a buffer zone of 1.5 m radius, as shown in Figure 11g, on 26 May 2017, the RE reduced from 14.02% to −1.51%. Similarly, on 18 June 2017, as shown in Figure 11h, the RMSD decreased from 0.04 to 0.02 m 3 ·m −3 and the RE reduced from 15.07% to −2.53%. Further, the comparison between different sizes of buffer zones indicates that the best match of R between SWI and in-situ SM measurements was with buffer zones of 1.5 m radius. 2017, as shown in Figure 11h, the RMSD decreased from 0.04 to 0.02 m 3 •m −3 and the RE reduced from 15.07% to −2.53%. Further, the comparison between different sizes of buffer zones indicates that the best match of R between SWI and in-situ SM measurements was with buffer zones of 1.5 m radius. Figure 11. Validation of the spatial variability of SM. The first and third columns (a,c,e,g,i,k) are on 26 May 2017 and the second and fourth columns (b,d,f,h,j,l) are on 16 June 2017. The blue color represents the "DT" approach, which considered theoretical dry and wet edges but not ra normalization. The red colors indicate the "DT/ra" approach, which considered ra normalization and theoretical dry/wet edges. The validation was conducted with buffer zones of the different radius: 0.15 m (initial spatial resolution) (a,b); 0.5 m (c-e); 1 m (e,f); 1.5 m (g,h); 2 m (i,j); and 2.5 m (k,l). Figure 12 shows the temporal variations of in-situ and UAS derived SM estimates for Soil Profile A (Figure 1) during the growing seasons of 2016 and 2017. In the non-growing season (November-April), this site had high SM and was energy limited, while, during the growing season (May-October), the condition switched to SM limited ( Figure 4). This supports the assumption of this study to apply NEF to estimate SM during the growing season. By comparing the "DT" and "DT/ra" approaches, we found that SM estimates from the "DT/ra" scheme had a close correspondence with the temporal variations of in-situ SM measurements at depths of 15 and 30 cm. Figure 11. Validation of the spatial variability of SM. The first and third columns (a,c,e,g,i,k) are on 26 May 2017 and the second and fourth columns (b,d,f,h,j,l) are on 16 June 2017. The blue color represents the "DT" approach, which considered theoretical dry and wet edges but not r a normalization. The red colors indicate the "DT/r a " approach, which considered r a normalization and theoretical dry/wet edges. The validation was conducted with buffer zones of the different radius: 0.15 m (initial spatial resolution) (a,b); 0.5 m (c-e); 1 m (e,f); 1.5 m (g,h); 2 m (i,j); and 2.5 m (k,l). Figure 12 shows the temporal variations of in-situ and UAS derived SM estimates for Soil Profile A (Figure 1) during the growing seasons of 2016 and 2017. In the non-growing season (November-April), this site had high SM and was energy limited, while, during the growing season (May-October), the condition switched to SM limited ( Figure 4). This supports the assumption of this study to apply NEF to estimate SM during the growing season. By comparing the "DT" and "DT/r a " approaches, we found that SM estimates from the "DT/r a " scheme had a close correspondence with the temporal variations of in-situ SM measurements at depths of 15 and 30 cm. The Taylor diagram and the scatterplot (Figure 13) shows the temporal validation of the "DT" and "DT/ra" approaches with SM at different depths (5, 15, 30 and 60 cm) for Soil Profiles A and B (Figure 1) with 1.5 m buffer zones. Significant improvements in R by normalizing ra is evident than without normalization. In the Taylor diagram (Figure 13a), all markers representing "DT/ra" approach (solid colors) had closer correspondence to the observation than markers representing the "DT" approach. Compared to the "DT" approach, R of the "DT/ra" approach increased from −0.25-0.1 to approximately 0.4-0.75. The "DT/ra" approach had the best matches with the in-situ SM at 15 and 30 cm depth, which corresponds to the root-zone depth of the willow plantation [50]. The R between the SM estimates and observations at the depth of 15 and 30 cm improved from the value around 0 (not significant correlation) to 0.7 (<0.05 significance level). Similarly, in Figure 13b, the comparison shows the "DT/ra" approach had the lowest RMSD around 0.025 m 3 •m −3 and RE about 3% when compared with the SM at 15 and 30 cm. However, in the "DT" approach, RMSDs were around 0.045 m 3 •m −3 and REs were about 10%. Generally, the improvement for the "DT/ra" approach in the temporal domain was substantially larger than the comparison in the spatial domain. This is because, in the temporal domain, hc changed more significantly from 0 to approximately 4 m. Further, regarding SM measurements at different depths, it can be found that the estimated SM from the triangle approach has better correlation with measurements at 15-30 cm depths, as shown by the green and light blue symbols in Figure 13. These results suggest that the estimated SM from the triangle approach reflects SM at the root zone. The Taylor diagram and the scatterplot (Figure 13) shows the temporal validation of the "DT" and "DT/r a " approaches with SM at different depths (5, 15, 30 and 60 cm) for Soil Profiles A and B (Figure 1) with 1.5 m buffer zones. Significant improvements in R by normalizing r a is evident than without normalization. In the Taylor diagram (Figure 13a), all markers representing "DT/r a " approach (solid colors) had closer correspondence to the observation than markers representing the "DT" approach. Compared to the "DT" approach, R of the "DT/r a " approach increased from −0.25-0.1 to approximately 0.4-0.75. The "DT/r a " approach had the best matches with the in-situ SM at 15 and 30 cm depth, which corresponds to the root-zone depth of the willow plantation [50]. The R between the SM estimates and observations at the depth of 15 and 30 cm improved from the value around 0 (not significant correlation) to 0.7 (<0.05 significance level). Similarly, in Figure 13b, the comparison shows the "DT/r a " approach had the lowest RMSD around 0.025 m 3 ·m −3 and RE about 3% when compared with the SM at 15 and 30 cm. However, in the "DT" approach, RMSDs were around 0.045 m 3 ·m −3 and REs were about 10%. Generally, the improvement for the "DT/r a " approach in the temporal domain was substantially larger than the comparison in the spatial domain. This is because, in the temporal domain, h c changed more significantly from 0 to approximately 4 m. Further, regarding SM measurements at different depths, it can be found that the estimated SM from the triangle approach has better correlation with measurements at 15-30 cm depths, as shown by the green and light blue symbols in Figure 13. These results suggest that the estimated SM from the triangle approach reflects SM at the root zone. Remote Sens. 2018, 10, x FOR PEER REVIEW 21 of 28

Linking Soil Moisture, Surface Heat Flux and the Triangle Approach
The temperature-vegetation triangle space is defined by the scatterplot of surface temperature (assuming a constant air temperature) and vegetation cover, where the triangle boundaries correspond to dry and wet endmembers of SM ( Figure 4). In this 2D space, the relative position of an observation of surface temperature between the wet/dry extremes has already been extensively applied to map SM using satellite data and at spatial scales of 60 m to 1 km [11,21,28]. Using high resolution imagery from UASs provides opportunity to match the typical spatial scale of SM variations (meter level) [51], identify areas of low SM levels where water is limiting for plants and will be obscured in coarser grids [52], and also interpret how fluxes aggregate within a footprint of the EC system.
In our study, we mapped changes in SM over a willow forest at canopy scales by considering a simple piecewise SM relation with sensible heat flux normalized by available energy at the scale of tree crowns (Figure 4) [1,53] together with observations of Ts and fc. The linkage between soil moisture and surface heat flux within the triangle approach shows that using a triangle approach for mapping SM requires considering any environmental factors susceptible to modify Ts (vertical axis of Figure  4) other than SM [11,21,28,29]. As it happens in our study site, it is possible that pixels with the same fractional cover have different aerodynamic properties due to tree height variations. In this framework, finding the absolute or "exact" value of the aerodynamic resistance is less important than accounting for relative variations among pixels with the similar vegetation cover, not being as critical as in SEB models that explicitly estimate sensible heat flux using aerodynamic resistances [54]. Additionally, as the aerodynamic properties of soil and vegetated surfaces are quite different, just accounting for those in the triangle for the full cover and bare soil cases corrected most of the original bias (spatial validation in Figure 11 and temporal validation in Figure 13). To do that, we used the formulation of bulk aerodynamic resistance under neutral conditions from Brutsaert with roughness effective parameters similar to Moran et al. (1994) and Hoffman et al. (2016) [16,18]. In their studies, those parameters were derived from static prescribed crop height for the entire area, while we used estimates of tree height derived from SfM photogrammetry techniques, which is a more accessible alternative to LiDAR to retrieve canopy structure information [55]. To have effective parameters of

Linking Soil Moisture, Surface Heat Flux and the Triangle Approach
The temperature-vegetation triangle space is defined by the scatterplot of surface temperature (assuming a constant air temperature) and vegetation cover, where the triangle boundaries correspond to dry and wet endmembers of SM ( Figure 4). In this 2D space, the relative position of an observation of surface temperature between the wet/dry extremes has already been extensively applied to map SM using satellite data and at spatial scales of 60 m to 1 km [11,21,28]. Using high resolution imagery from UASs provides opportunity to match the typical spatial scale of SM variations (meter level) [51], identify areas of low SM levels where water is limiting for plants and will be obscured in coarser grids [52], and also interpret how fluxes aggregate within a footprint of the EC system.
In our study, we mapped changes in SM over a willow forest at canopy scales by considering a simple piecewise SM relation with sensible heat flux normalized by available energy at the scale of tree crowns ( Figure 4) [1,53] together with observations of T s and f c . The linkage between soil moisture and surface heat flux within the triangle approach shows that using a triangle approach for mapping SM requires considering any environmental factors susceptible to modify T s (vertical axis of Figure 4) other than SM [11,21,28,29]. As it happens in our study site, it is possible that pixels with the same fractional cover have different aerodynamic properties due to tree height variations. In this framework, finding the absolute or "exact" value of the aerodynamic resistance is less important than accounting for relative variations among pixels with the similar vegetation cover, not being as critical as in SEB models that explicitly estimate sensible heat flux using aerodynamic resistances [54]. Additionally, as the aerodynamic properties of soil and vegetated surfaces are quite different, just accounting for those in the triangle for the full cover and bare soil cases corrected most of the original bias (spatial validation in Figure 11 and temporal validation in Figure 13). To do that, we used the formulation of bulk aerodynamic resistance under neutral conditions from Brutsaert with roughness effective parameters similar to Moran et al. (1994) and Hoffman et al. (2016) [16,18]. In their studies, those parameters were derived from static prescribed crop height for the entire area, while we used estimates of tree height derived from SfM photogrammetry techniques, which is a more accessible alternative to LiDAR to retrieve canopy structure information [55]. To have effective parameters of roughness length at blending height for momentum and displacement height, we averaged tree height (e.g., roughness elements) over the footprint area in each date. The SM retrieved showed larger improvements in accuracy over time more than over space. This is logical, as we did not consider the effect of fine scale spatial variability (canopy level) and because the tree height variation is very high (0-4 m) over time in this willow forest.

Influence of the 3D Canopy Structure to the Triangle Approach at Fine Resolution
The main improvement of SM estimates is on the temporal domain, in which h c changes from 0 to around 4 m. In the spatial domain, there was an improvement in RMSD and RE, but no changes in the correlation coefficient. This is due to that the averaged h c from the whole willow patch was used to calculate the bulk aerodynamic roughness and a constant kB −1 of 2.3 was used to account for the local boundary resistance for the heat transfer. An important aspect to further improve the triangle approach is the proper interpretation of the spatial patterns of T s detected from UAS at fine scales of 1 m [56]. These fine spatial T s patterns vary in response to the water status of tree crowns, as trees modify root-zone SM through transpiration [52]. However, the spatial patterns are also influenced by the micro-scale structural heterogeneity of forest canopies due to gaps or tree height changes [57]. Historically, it has been considered that turbulence mixing "wipes" out such fine scale heterogeneities, but it has been shown that the eddy mixing length inside and just above the canopy sublayer is at the same scale as the tree-crown diameter and, therefore, crown-scale heterogeneity leads to persistent effects [57,58]. In fact, the spatial arrangement of soil and vegetation patches and plant architecture produces near-surface turbulent changes clearly modifying the surface fluxes and not accounted for in turbulent diffusion theory formulations [56,59].
Typically, the effect of tree scale heterogeneity when calculating turbulent fluxes is accounted via length scales using a scalar quantity (zero plane displacement and aerodynamic roughness length) that aim to account for the 3D effects of the canopy structure on the momentum transfer, using effective patch-averaged profiles when the goal is to estimate the fluxes at the blending height as we also did over time [57]. However, in this study, a zero-dimensional representation of 3D canopy effects using uniform scalar properties removes the signature of heterogeneity on fluxes and temperature when heterogeneity is relevant as in the case of mapping fine scale SM patterns.
One simple possibility to incorporate the effect of canopy heterogeneity in the triangle space would be to consider an additional boundary layer resistance to heat transport in series with the aerodynamic resistance, to account for variations in heat transfer due changes in tree height that would modify heat roughness length as in Equations (34) and (35). That means that, in Equations (34) and (35), the averaged h c can be used to estimate z om , while the pixel h c can be used to calculate z oh . Thus, the effective spatially-averaged h c is used for the aerodynamic resistance between the reference height and the displacement height and the local variation of h c with respect to the average h c is used to account for the local boundary resistance for Z oh . In practice, this is equivalent to calculating a kB −1 resistance that is dependent on tree height.
where h c is the averaged canopy height for the willow planation and ·h c is the local variation of the canopy height.
To explore if accounting for the spatial variability of h c had an impact on the retrieval of SM, we compared the local h c formulation of Equation (34) with the two prior approaches (no resistance DT and effective parameterized resistance DT/r a with average h c ) on 18 June 2017 as it had more TDR measurements available. The results are shown in Table 2. It can be seen that, when including local variations of h c at the pixel level, there is a slight improvement. This indicates the variation of h c inside the plantation is also important, even though the spatial variations are not that large. It can be expected that, for other sites with larger spatial variations of h c , more significant improvement could be achieved. The aerodynamic roughness length in this study was expressed as a ratio of h c [45]. This is valid for homogeneous surfaces as the eddy covariance flux site of this study [33]. However, for the heterogeneous landscapes, additional consideration of the effective averaged obstacle height and frontal area index can help to obtain more accuracy [58]. Moreover, h c can also be more accurately acquired from UAS LiDAR techniques. For instance, the comparison between UAS-based LiDAR and SfM techniques shows that LiDAR is relatively more accurate in the forest height and 3D structure detection, while SfM is an adequate low-cost alternative [60].
Another aspect related with SM retrieval from the relation with surface heat fluxes is due to the fact that every tree has a different root system and root depth; therefore, T s pattern will respond to slightly different SM levels from tree to tree. That might explain why we get the better correlations in the layers with higher root density. It can be expected that, in this study, during the progress of the growing season from bare soil to full canopy cover, UAS-based SWI is associated with SM in different depths. For the bare soil period, SWI indicates SM more superficial, while for high canopy covers, SWI was related with SM at deeper depths. Since this study has more UAS observations with high canopy covers (f c > 0.7) over Soil Profiles A and B, it explained better correlations with SM at layer around 15-30 cm depths. This finding agrees with other studies that showed SM estimates from the triangle approach to be more related to the root-zone SM [19,61]. For instance, Wang et al. (2018) compared the SM estimates from UAS with in-situ SM at different depths at the experimental plots with spring wheat and found the best match between estimates and observations to be at 10-20 cm depth, which corresponds to the root zone of the wheat [19].
The comparison between different sizes of buffer zones indicates that the best match of R was with buffer zones of 1.5 m radius. Several reasons can justify the best match at this size of buffer zones. The triangle approach infers SM through canopy temperature for high f c levels. Although UAS imagery provided sub-meter level VHR imagery of the canopy, there may be a spatial discrepancy between the canopy status and SM conditions underneath the canopy. This means the detected canopy status may reflect the SM conditions representing the few meters around the canopy instead of the SM condition just below the canopy and, in fact, the root of the willow plantation can spread 1-3 m around the tree [62]. Further, a complete canopy is composed of sunlit and shaded leaves and a single pixel among the canopy cannot reflect the whole status of the canopy and SM conditions in the root zone [63]. According to the DSM and visual inspection of the RGB photos, as shown in Figure 1, the radiuses of the canopy are around 1-1.5 m in dense vegetation conditions. This approximately agrees with the size of the buffer zones. Additionally, there are uncertainties in the georeferencing of the UAS image processing and the Trimble RTK GNSS, which was used to measure the location of TDR measurements and GCPs, also has errors. Therefore, errors were reduced with spatial aggregation.

Comparison with Other Studies
While numerous studies have been conducted with thermal inertia, the triangle approach and microwave to estimate SM from satellite and manned airborne remote sensing [11,36,63,64], only few studies have been conducted with UAS observations [4,17,19]. For instance, Sobrino et al. (2012) used NDVI and T s with a polynomial formulation to predict SM and this approach can estimate SM with RMSDs of 0.05 m 3 ·m −3 from observations from a manned airborne system [64]. Fan et al. (2015) used vegetation indices and T s to estimate SM with RMSDs around 0.023 m 3 ·m −3 from a manned airborne system [63].  demonstrated that UAS-based optical and thermal remote sensing data can be utilized with the water deficit index approach to estimate SM with R 2 reaching to 0.63 and RMSDs less than 0.1 m 3 ·m −3 [19]. Compared to these studies, this study achieved an accuracy to map spatial and temporal variations of SM with R around 0.58-0.69 and RMSDs around 0.025 m 3 ·m −3 . Further, this study was based on the triangle approach to utilize the optical and thermal remote sensing information. It has advantages to be applied for bare soil, sparse and dense vegetation conditions. Moreover, this approach can be applied to monitor the root-zone SM, which has added significance for predicting near-future vegetation anomalies [65]. Further, the approach of this study does not need parameter calibration and most inputs were available from UAS observations. This UAS-based SM monitoring approach can be applied in the data-scarcity regions. However, the shortage of this method is that, in the dense vegetation conditions, the accuracy of this method is very sensitive to the accuracy of the optical and thermal data, as shown in Figure 8. For the miniaturized UAS thermal sensors, it is hard to obtain surface temperature with very high accuracy (<1 • C). Thus, the precise calibration of the thermal sensor, as performed in this study, is recommended. Additionally, novel trapezoid methods utilizing the visible and shortwave infrared data may be another alternative [66,67].

Conclusions
In this study, an operational methodology to map soil moisture (SM) in the root-zone at high spatial resolution (m level) from Unmanned Aerial Systems (UAS) is proposed. We improved the state-of-the-art temperature-vegetation triangle methodology by normalizing the land-atmosphere temperature gradients with a bulk aerodynamic resistance that was derived from structure-from-motion based surface elevation maps. In this way, the effects of the tree growth driven surface roughness dynamics on surface temperatures were considered. Due to the limited coverage of the UAS, the hydrological contrast of the UAS maps was reduced, which lead us to model theoretical dry and wet edges rather than extracting them from the images. Evaluation with SM measurements in a Danish bioenergy willow short rotation coppice in two consecutive summer seasons led us to the following conclusions: (i) Accounting for changes in surface roughness considerably increased the performance of the temperature-vegetation triangle approach to estimate SM variation. (ii) Under densely vegetated conditions, the estimated SM from the triangle approach correlated better with the root-zone SM than with the surface SM (< 10 cm). This shows that the modified triangle approach is better suited to monitor water availability for vegetation compared to thermal inertia-or microwave-based approaches, which detect surface rather than deeper root-zone SM. (iii) The optimum spatial resolution or aggregation for the modified triangle approach is in the same order of magnitude as the typical length scale of the tree crowns. (iv) Given the high sensitivity of SM model to errors in surface temperature estimates, we proposed a pixel-wise sensor calibration method to improve the accuracy of the uncooled UAS thermal sensor to be around 0.55 • C for the laboratory conditions and 0.95 • C in the field. (v) The proposed methodology mainly relies on UAS observations and requires limited in-situ measurements, e.g., solar incoming radiation, air temperature and humidity, wind speed, and atmospheric pressure, and can be operationally used for routine SM monitoring in both agricultural and natural ecosystems.
This study mainly addresses the influence of the temporal changes of the canopy roughness on the triangle approach. Future research should consider the effects of strong spatial variation of the canopy height and the role of canopy gaps for the fine scale patterns of the surface temperature observed from UAS, to avoid contamination of the estimated SM maps by micro-scale variability of structural induced turbulence patterns and radiation budgets.