How High to Fly? Mapping Evapotranspiration from Remotely Piloted Aircrafts at Different Elevations

: Recent advancements in remotely piloted aircrafts (RPAs) have made frequent, low-ﬂying imagery collection more economical and feasible than ever before. The goal of this work was to create, compare, and quantify uncertainty associated with evapotranspiration (ET) maps generated from different conditions and image capture elevations. We collected optical and thermal data from a commercially irrigated potato ( Solanum tuberosum ) ﬁeld in the Wisconsin Central Sands using a quadcopter RPA system and combined multispectral/thermal camera. We conducted eight mission sets (24 total missions) during the 2019 growing season. Each mission set included ﬂights at 90, 60, and 30 m above ground level. Ground reference measurements of surface temperature and soil moisture were collected throughout the domain within 15 min of each RPA mission set. Evapotranspiration values were modeled from the ﬂight data using the High-Resolution Mapping of Evapotranspiration (HRMET) model. We compared HRMET-derived ET estimates to an Eddy Covariance system within the ﬂight domain. Additionally, we assessed uncertainty for each ﬂight using a Monte Carlo approach. Results indicate that the primary source of uncertainty in ET estimates was the optical and thermal data. Despite some additional detectable features at low elevation, we conclude that the tradeoff in resources and computation does not currently justify low elevation ﬂights for annual vegetable crop management in the Midwest USA.


Introduction
Irrigated agriculture uses approximately 162 billion cubic meters of water per year in the United States and accounts for 37% of total water use [1]. Of the 162 billion cubic meters per year, 70% or approximately 113 billion cubic meters enter the atmosphere as evapotranspiration (ET). The remaining 30% either runs off the surface or recharges aquifers. Fertilizer, herbicide, and pesticide contamination in runoff and groundwater have been shown to have negative impacts on ecosystems [2,3]. Precision irrigation practices have gained popularity for their ability to maximize yield while minimizing environmental impacts by spatially optimizing water and nutrient application [4,5]. By applying water precisely, overwatering is reduced and nitrogen can remain in the rootzone instead of leaching into groundwater [6][7][8]. Spatial data or maps are required to optimize spatially variable irrigation application. These maps may include parameters that are spatially variable (soil texture, surface elevation, electric conductivity) [9][10][11] and/or temporally variable (crop height, Leaf Area Index, crop water use, surface temperature) [12][13][14]. Collecting these observations at high spatial and temporal resolution over farms remains a challenge.
These different types of spatial data have been collected remotely for decades, leading to a plethora of remote sensing platforms available today. Satellite and airplane mounted cameras have been the most common platforms for agriculture over the last few decades, but can be cost and skill prohibitive, difficult to schedule frequently, prone to cloud interference (especially in the Midwest United States), and do not always have enough resolution to adequately represent smaller fields or their features [9,15]. The growing market in remotely piloted aircraft (RPA) technology has reduced the cost of entry for custom imagery collection. Remotely piloted aircrafts are defined as 'an unmanned aircraft that is piloted from a remote pilot station' [16]. We refer to these systems as 'drones' throughout the rest of the paper.
Frequently and rapidly generated maps of ET can be used as decision support tools [9,17] by showing how much water is being used in a given area and helping to determine how much water should be applied. High-resolution maps of ET and crop water stress can be created from remote sensing platforms with multispectral and thermal imagery along with meteorological data [18][19][20]. Vegetation indices, created from multispectral imagery, as well as surface temperature imagery are the main inputs for many energy balance models that represent soil and vegetation as two separate sources in a pixel [21,22] or as a single source in a pixel [23,24]. While data for vegetation indices are easier to collect, they can become 'green saturated' in dense canopies such as potato [25]. Thermal-based indices and models have a stronger relationship to stomatal conductance [26] and have been shown to detect plant water stress at an earlier, actionable time frame compared to vegetation indices [27]. Remotely sensed ET models for new crops, regions, and platforms benefit from comparison to ground-based methods of ET measurement. This can be done with physical water budget methods (lysimeters) [28,29] or micrometeorological methods (Eddy covariance or surface renewal) [19,30]. A main challenge with this is making sure that both methods are measuring the same fluxes or changes. Uncertainty maps are an often overlooked, yet insightful tool for evaluating remote sensing ET models. Uncertainty maps provide a visual layout of where the majority of ET uncertainty lies in a field, outlining the areas in which the model performs well and where it is less reliable. Barring direct application to agricultural fields where the time-consuming nature would be unreasonable, uncertainty maps have valuable insights for model development that should not be overlooked. While ET mapping has been used for decades, the majority of the tools were developed for semi-arid regions; and still today there are limited ET mapping tools developed for the Midwest United States [31].
Recently, there has been an influx of ET maps generated by drone platforms [32][33][34]. The lower elevation, when compared to satellites and airplanes, allows for higher resolution imagery with similar cameras, and drone platforms have the benefit of low experience entry, adaptable flight scheduling, lower price thresholds, and high-resolution imagery [15]. However, flying drones low over agricultural fields to collect imagery at maximum resolution increases flight time, battery, and data processing requirements. Though recent work has focused on optimizing spatial grid sizes for ET mapping algorithms [32,35], there is limited research exploring whether increasingly higher resolution ET maps provide additional insights and decision support. Here, we asked: how reliably and under what flight and field conditions can thermal and visible imaging drones estimate crop ET?
We evaluated this question in a potato (Solanum tuberosum) field in the Central Sands region of Wisconsin, USA. The goals for this study were to (1) develop a method for highresolution mapping of ET using drone imagery for potatoes over the course of a field season using the High-Resolution Mapping of Evapotranspiration (HRMET) model; (2)  (3) evaluate the uncertainty of spatial ET estimates and the impact of capture elevation on uncertainty to determine the most useful mission design in drone applications for crop ET.

Site Description
Wisconsin is the third largest producer of potatoes in the United States, with 28,000 ha harvested and over 300 million USD in sales in 2019 [36,37]. Potatoes are one of the main agricultural products of the Wisconsin Central Sands. The Wisconsin Central Sands region (8860 km 2 ) is known for its high sand content soil, humid climate, and high agricultural presence [38]. Temperatures in this region have an annual mean value of 6.6 • C with a mean January minimum of −19.2 • C and a mean August maximum of 27.1 • C [38]. On average, the area receives 830 mm of rain and 1140 mm of snow annually [38]. The storage capacity of the highly conductive sand aquifer and annual precipitation makes the area prime agricultural land. However, widespread agriculture and strong surface-ground water connection have led to water quality and groundwater depletion issues [39]. Namely, there are challenges with increasing groundwater nitrate concentrations, shallow private wells not being able to access water, and groundwater dependent surface waterbodies going dry [39]. Recent studies have shown that a reduction in frequency and magnitude of groundwater pumping for irrigation could reduce the impacts on nearby surface water in the Central Sands by reducing the amount of nitrate leached and runoff from the soil [28,40]. Precision irrigation has the potential to maintain agricultural productivity while reducing the magnitude and frequency of groundwater pumping in this region.
Drone imagery and ground reference data were collected to map ET from a 0.65-ha subsection of a 7.28-ha potato field in the Wisconsin Central Sands during the 2019 growing season ( Figure 1). The subsection was centered around an Eddy Covariance system to capture the flux tower's footprint (Section 2.7). Meteorological data were collected 3.9 km away from the field site (Section 2.5). Drone imagery and ground references were taken in partnership with an 11,000-ha commercial potato and vegetable farm. The field was planted 25 April, uniformly irrigated using a center pivot system, vine killed 18 August, and harvested 7 September. It rained almost weekly during the growing season and had a daily average reference ET of 3.62 mm, with a max of 5.30 mm on 10 July and a minimum of 0.96 mm on 26 August (Figures 2 and 3). Little spatial variability was expected in the domain due to its flatness and relatively homogenous loamy sand soil. season using the High-Resolution Mapping of Evapotranspiration (HRMET) model; (2) evaluate HRMET in potato fields through comparison with eddy covariance (EC) tower ET estimates; (3) evaluate the uncertainty of spatial ET estimates and the impact of capture elevation on uncertainty to determine the most useful mission design in drone applications for crop ET.

Site Description
Wisconsin is the third largest producer of potatoes in the United States, with 28,000 ha harvested and over 300 million USD in sales in 2019 [36,37]. Potatoes are one of the main agricultural products of the Wisconsin Central Sands. The Wisconsin Central Sands region (8860 km 2 ) is known for its high sand content soil, humid climate, and high agricultural presence [38]. Temperatures in this region have an annual mean value of 6.6 °C with a mean January minimum of −19.2 °C and a mean August maximum of 27.1 °C [38]. On average, the area receives 830 mm of rain and 1140 mm of snow annually [38]. The storage capacity of the highly conductive sand aquifer and annual precipitation makes the area prime agricultural land. However, widespread agriculture and strong surfaceground water connection have led to water quality and groundwater depletion issues [39]. Namely, there are challenges with increasing groundwater nitrate concentrations, shallow private wells not being able to access water, and groundwater dependent surface waterbodies going dry [39]. Recent studies have shown that a reduction in frequency and magnitude of groundwater pumping for irrigation could reduce the impacts on nearby surface water in the Central Sands by reducing the amount of nitrate leached and runoff from the soil [28,40]. Precision irrigation has the potential to maintain agricultural productivity while reducing the magnitude and frequency of groundwater pumping in this region.
Drone imagery and ground reference data were collected to map ET from a 0.65-ha subsection of a 7.28-ha potato field in the Wisconsin Central Sands during the 2019 growing season (Figure 1). The subsection was centered around an Eddy Covariance system to capture the flux tower's footprint (Section 2.7). Meteorological data were collected 3.9 km away from the field site (Section 2.5). Drone imagery and ground references were taken in partnership with an 11,000-ha commercial potato and vegetable farm. The field was planted 25 April, uniformly irrigated using a center pivot system, vine killed 18 August, and harvested 7 September. It rained almost weekly during the growing season and had a daily average reference ET of 3.62 mm, with a max of 5.30 mm on 10 July and a minimum of 0.96 mm on 26 August (Figures 2 and 3). Little spatial variability was expected in the domain due to its flatness and relatively homogenous loamy sand soil.         Figure 3. Growing season timeseries of irrigation (mm) and precipitation (mm) for 2019. Precipitation is in hourly increments and the timing of irrigation events are associated with the irrigation initiation (duration is not depicted). Mission sets are highlighted with black vertical lines. Imagery was not collected while the center pivot system was running.

Drone Missions and Data Processing
Eight mission sets were flown over the 2019 growing season with three different altitude (90, 60, and 30 m) flights per mission set (Table 1). Flight mission sets were scheduled to occur once a week during the least cloudy and windy day. We collected imagery over a 0.65-ha subsection of the field surrounding the Eddy Covariance tower (Figure 1). Thermal and multispectral data were collected with an 80% front and side overlap using an Altum combined multispectral-thermal camera (Micasense Inc., Seattle, WA, USA), mounted on a Matrice 210 RTK quadcopter (Shenzhen DJI Sciences and Technologies Ltd., Shenzen, Guangdong, China). The Altum camera captures five multispectral bands centered at 475, 560, 668, 717, and 842 nm. The multispectral data were captured at a resolution of 3.88, 2.59, and 1.29 cm for 90, 60, and 30 m, respectively. The Altum camera's thermal sensor is centered at 11,000 nm. The thermal imagery was captured at a resolution of 60.9, 40.6, and 20.3 cm for 90, 60, 30 m, respectively. Atmospheric corrections were not conducted on the thermal imagery due to the low elevation of the flights. Rather, the differences between the thermal imagery and ground truth data points were used in the Monte Carlo uncertainty estimation procedure to help understand whether corrections are needed for low elevation remote sensing. We processed images via the Micasense Altum processing workflow recommended for use with Agisoft Metashape Professional (Agisoft LLC, St. Petersburg, Russia). The multispectral data were calibrated in Metashape using images of a reflectance calibration panel (Micasense Inc., Seattle, WA, USA) before and after every flight. Images were then aligned geographically with high accuracy settings. A dense cloud was created using medium quality and aggressive depth filtering. A DEM was created from the dense cloud, and then used the surface for the ortho-mosaics. The ortho-mosaics were exported as geotiff files. The changes in viewing angle with elevation were considered negligible due to the low elevation of the flights and the short height of the potato canopy [41].

Thermal Ground Reference
Thermal ground reference data were collected for six of the eight mission sets. Samples were taken in four locations (three vegetation, one bare soil) immediately following the last flight in the mission set. Surface temperature values were collected using a Thermo Fisher handheld infrared thermometer (Thermo Fisher Scientific, Waltham, MA, USA). Five samples were taken at the top of the canopy (orthogonal) and averaged at each location for the six mission sets. We extracted remotely sensed temperature values using the average temperature of the 1.5 m diameter circle surrounding the ground reference sampling location. The ground reference temperatures and corresponding remotely sensed temperature values were used to calculate root mean square error (RMSE) for each flight elevation. RMSE values were 3.78, 4.42, and 3.44 • C for 90, 60, and 30 m respectively (Appendix A, Figure A1). July 8th had the highest RMSE values, while August 6th and 8th had the lowest RMSE values (Appendix A, Table A1). Atmospheric calibration was bypassed to assess uncertainty using ground reference comparisons for Monte Carlo Analyses.

Crop Phenology
We estimated leaf area index (LAI) and height values per pixel using established empirical models with the Enhanced Vegetative Index (EVI) developed for potato in the Wisconsin Central Sands region [2]. The multispectral images were aggregated to 5-m resolution to represent both plants and soils in each pixel [33]. EVI was computed from the aggregated imagery. The residual error from the empirical model was used to create ensemble data for Monte Carlo analyses. The LAI and height layers were disaggregated to the size of the thermal values and used to generate ET maps to restore resolution to the model while still accounting for the lowest resolution input (thermal).

Meteorological Data Collection
Meteorological data were collected at a nearby weather station operated by the University of Wisconsin-Madison and University of Michigan. The meteorological data were used from this station rather than the EC tower due to crucial differences in temporal resolution and a more continuous dataset. Instruments measuring air temperature (HMP155A; Campbell Scientific, Logan, UT, USA), relative humidity (HMP155A; Campbell Scientific, Logan, UT, USA), and solar radiation (LI200R Pyranometer; LI-COR Biosciences, Lincoln, NE, USA) were mounted 1.5 m above a short grass surface. The instrument measuring wind speed (03002-L Wind Sentry Set; Campbell Scientific, Logan, UT, USA) was mounted 3 m above a short grass surface. Measurements for air temperature, relative humidity, and wind speed were taken every 5 min. Measurements for solar radiation were recorded hourly. Vapor pressure was calculated using relative humidity and air temperature values. Daily reference ET was reported from the station [42].
Mean values of air temperature, vapor pressure, and wind speed used in HRMET model simulations were extracted 30 min before and after the midpoint of the respective mission set. The values for solar radiation used in HRMET model simulations were linearly regressed and calculated using the measurements recorded (1) closest to the middle of the set, (2) proceeding recording, and (3) following recording. Uncertainty in air temperature, vapor pressure, and wind speed was estimated using the standard deviation of the measurements before and after takeoff. Uncertainty in solar radiation was estimated as the root mean square error (RMSE) of the linear regression.

Remotely Sensed Maps of ET and Uncertainty Estimation
Many remotely sensed ET models rely on the designation of 'hot' and 'cold' pixels that scale the minimum and maximum limits of ET within the domain [21,23]. While this scaling requires less input data [33] and reduces uncertainty from meteorological inputs [29], it can create challenges in homogenous agricultural fields [43]. Complicating things further, irrigated fields often aim to minimize water stress, making it harder to find 'hot' pixels for scaling, especially in humid regions. This lack of a 'hot' pixel can be particularly troublesome for precision irrigation with small imagery domains. We use the High-Resolution Mapping of Evapotranspiration (HRMET) model to mitigate some of these issues. HRMET is a singular, process-based energy model that uses canopy Remote Sens. 2022, 14, 1660 7 of 21 characteristics, surface temperature, and local meteorological data to estimate latent heat flux as the residual of the energy balance, and then converted to ET. Available energy for vegetation and soil is calculated using a two-source energy algorithm while sensible heat flux is iteratively solved as a single source for each pixel. This iteration removes the need to find 'hot' and 'cold' pixels, making it useful for developing ET maps for small fields. HRMET was developed for airplane imagery and has already been used in potatoes [9], peas (Pisum sativum), sweet corn (Zea mays), and field corn (Zea mays) [43].
ET maps were generated for each flight using a version of the HRMET model calibrated for potatoes in R Studio (RStudio, Boston, MA, USA). For information on the HRMET model, refer to Zipper and Lohiede, 2014 [43]. HRMET is an open-source model and the code for R and MATLAB are available at https://github.com/samzipper/HRMET (Last accessed on 20 January 2020). HRMET is run on a pixel-by-pixel basis. The three flight elevations had the same flight domain but had different pixel counts depending on resolution. Pixel count averaged 2,120,000, 4,890,000, and 20,740,000 for the 90, 60, and 30 m flights, respectively. HRMET ET values were assessed by mean, median, minimum, and maximum across the spatial footprint of the EC tower.
We estimated the influence of uncertainty from the model inputs using a Monte Carlo approach [44]. The mean and standard deviation of all spatial, meteorological, and empirical inputs were used to create an ensemble of input data based on the standard normal distribution. HRMET was run using randomly generated values from the normal distribution ensemble until 95% of the pixels had reach stability for at least 10 runs. A running average and standard deviation were calculated after every iteration. Stability was defined when the average value for a pixel was within 0.1% its value in the regular HRMET run. It took on average 300 iterations to reach scene stability. We used the Monte Carlo analyses to generate maps of the standard deviation of ET per pixel. We ran the analyses using an Intel i9 processor (Intel, Santa Clara, CA, USA), an Nvidia Quadro P620 graphics card (Nvidia, Santa Clara, CA, USA), and 64GB Corsair Vengeance LPX RAM (Corsair, Fremont, CA, USA). Because of the long runtime for these calculations (taking approximately 1.5 weeks for 90 m missions, one month for test 60 m missions, and 3 months for test 30 m missions) and no further insights from test 60 and 30 m missions, we limited Monte Carlo analyses to only the 90 m missions. Each HRMET input and its data sources are described in Table 2.  We calculated a Relative ET index [43] to compare intra-plot variability between different flights. Relative ET was calculated by linearly normalizing the mean to the 5th and the 95th percentile for each mission. Values closer to 1 indicate the highest ET rates within the domain and the values closer to 0 indicate the lowest ET rates within the domain.

Eddy Covariance
Fluxes of sensible and latent heat were collected using an Eddy Covariance (EC) tower at a half-hourly rate [45]. Water vapor, heat, and turbulent winds were measured with a Campbell Scientific, Inc. IRGASON, along with standard meteorological, surface radiation, and soil measurements. The instruments were located at 2.7 m height and were approximated to have a footprint radius ten times the instrument height (r = 27.4 m) [46]. The time series was aggregated to an hourly scale to match the time scale of the HRMET ET algorithm. The EC ET data were then linearly regressed using the recording closest to each flight's midpoint and the recording before and after, to find the ET value nearest the flight. Ground heat flux was not collected during this period and the fluxes were assumed to have a closure of 10% from the random turbulent uncertainty [47].

Spatial Features
Notable features are present in the imagery domain (see any raster figure, i.e., Figure 4). The thick, vertical line that bisects the domain is the field road providing access to the center pivot irrigation system. The thin, horizontal line that bisects the imagery domain is a track line left behind by the center pivot system. The squared shape bare soil location in the bottom left corner of the domain is an area that the potato planter missed; this area was later planted and the canopy closed (as can be seen in the later mission sets). These features may or may not be related to changes in ET. We observed and attributed some miscellaneous dark spots that appear on some days but not at all the flight elevations as overhead cloud interference. flight. Ground heat flux was not collected during this period and the fluxes were assumed to have a closure of 10% from the random turbulent uncertainty [47].

Spatial Features
Notable features are present in the imagery domain (see any raster figure, i.e., Figure  4). The thick, vertical line that bisects the domain is the field road providing access to the center pivot irrigation system. The thin, horizontal line that bisects the imagery domain is a track line left behind by the center pivot system. The squared shape bare soil location in the bottom left corner of the domain is an area that the potato planter missed; this area was later planted and the canopy closed (as can be seen in the later mission sets). These features may or may not be related to changes in ET. We observed and attributed some miscellaneous dark spots that appear on some days but not at all the flight elevations as overhead cloud interference.

EVI, Height and LAI
Values for EVI peaked on 8 July, decreasing steadily with the senescence of the potato canopy (Appendix B, Figure A2. There was little variability in EVI between flight elevations on days that were not already associated with scattered clouds. As flight elevations

EVI, Height and LAI
Values for EVI peaked on 8 July, decreasing steadily with the senescence of the potato canopy (Appendix B, Figure A2. There was little variability in EVI between flight elevations on days that were not already associated with scattered clouds. As flight elevations decreased, there was an increase in the ability to visually discern individual rows, but no new unique features appeared in the EVI imagery with increased resolution. Values for height peaked on 8 July at 0.28 m, decreasing steadily with the senescence of the potato canopy (Appendix B, Figure A3). There was little variability between flight elevations on days that were not already associated with scattered clouds. As flight elevation decreased there was an increase in visually discernable individual rows based on their height, but no new unique features appear. Across cropped areas, the values visually appear to be more homogenous than EVI.
Values for LAI peaked on the 8 July at 4.75 m 2 m −2 , only slightly decreasing with the senescence of the potato canopy ( Figure 4). There was little variability between flight elevations on days that were not already associated with scattered clouds. As flight elevation decreased there was an increase visually discernable individual rows based on the higher resolution of the bands prior to aggregation, but no new unique features appear. Across cropped areas, the values visually appear to be more homogenous than EVI, especially around their peak days.

HRMET ET
HRMET ET was highest on 8 July at 0.9 mm h −1 and lowest on 21 June at 0.3 mm h −1 (Figure 6). ET value ranges had no observable differences across the flight elevations, with the exception of the mission sets that had cloud interference ( Table 3). As flight elevation decreased, there was an increase in visually discernable individual rows and their differences in ET. Missions that had either irrigation earlier on the day of flights (21 June and 8 July) or heavy rainfall 24 h prior to the flights (6 August), had the least variability of ET throughout the scene.

HRMET ET
HRMET ET was highest on 8 July at 0.9 mm h −1 and lowest on 21 June at 0.3 mm h −1 (Figure 6). ET value ranges had no observable differences across the flight elevations, with the exception of the mission sets that had cloud interference ( Table 3). As flight elevation decreased, there was an increase in visually discernable individual rows and their differences in ET. Missions that had either irrigation earlier on the day of flights (21 June and 8 July) or heavy rainfall 24 h prior to the flights (6 August), had the least variability of ET throughout the scene.
HRMET ET was highest on 8 July at 0.9 mm h and lowest on 21 June at 0.3 mm h ( Figure 6). ET value ranges had no observable differences across the flight elevations, with the exception of the mission sets that had cloud interference ( Table 3). As flight elevation decreased, there was an increase in visually discernable individual rows and their differences in ET. Missions that had either irrigation earlier on the day of flights (21 June and 8 July) or heavy rainfall 24 h prior to the flights (6 August), had the least variability of ET throughout the scene.    The relative ET plots accentuate the variability in the cropped areas of the plot (Figure 7). Individual rows become more visible, but at the cost of the cloud interference also becoming more visible (6 August, 60 m). There are variations in relative ET across the scene between elevations for the same mission sets, but these variations are not seen in the HRMET ET plots (17 June, 21 June, 25 June, and 2 August). Lower elevation flights provide greater delineation between higher and lower values. This replaces averaged pixels with more 'extreme' values. While across the scene, the values remain the same, new patterns emerge with higher resolution. The relative ET plots accentuate the variability in the cropped areas of the plot (Figure 7). Individual rows become more visible, but at the cost of the cloud interference also becoming more visible (6 August, 60 m). There are variations in relative ET across the scene between elevations for the same mission sets, but these variations are not seen in the HRMET ET plots (17 June, 21 June, 25 June, and 2 August). Lower elevation flights provide greater delineation between higher and lower values. This replaces averaged pixels with more 'extreme' values. While across the scene, the values remain the same, new patterns emerge with higher resolution.

Monte Carlo Uncertainty Assessment
The number of Monte Carlo simulations required to reach 95% match stability varied, though all 90-m missions converged on a solution (Figure 8). The days with greater variability in meteorological inputs tended to require more iterations to reach stability ( Figure  8). Days with variable solar radiation, such as 22 July, took almost 450 iterations. Days with low solar radiation and low wind speeds, such as 17 June and 6 August, took only 225 iterations. Using the Monte Carlo analyses, we found a per pixel standard deviation of 0.02-0.11 mm h −1 across all 90-m ET maps (Figure 9). Most missions were visually comparable with the notable exception of 22 July which had observably higher standard deviation. The spatial patterns in uncertainty across the 90-m missions signify that the dominant sources of uncertainty in the imagery are the spatial data. This would be the uncertainty of the spatial inputs (LAI, height, surface temperature) rather than the cameras.

Monte Carlo Uncertainty Assessment
The number of Monte Carlo simulations required to reach 95% match stability varied, though all 90-m missions converged on a solution (Figure 8). The days with greater variability in meteorological inputs tended to require more iterations to reach stability ( Figure 8). Days with variable solar radiation, such as 22 July, took almost 450 iterations. Days with low solar radiation and low wind speeds, such as 17 June and 6 August, took only 225 iterations. Using the Monte Carlo analyses, we found a per pixel standard deviation of 0.02-0.11 mm h −1 across all 90-m ET maps (Figure 9). Most missions were visually comparable with the notable exception of 22 July which had observably higher standard deviation. The spatial patterns in uncertainty across the 90-m missions signify that the dominant sources of uncertainty in the imagery are the spatial data. This would be the uncertainty of the spatial inputs (LAI, height, surface temperature) rather than the cameras. Uncertainty originating from the camera would have resulted in a more uniform distribution. The bare soil portions of the domain tended to show higher uncertainty. The higher uncertainty and meteorological variability on 22 July indicate that the differences in ET may be more related to temporal difference than plant stomatal control. The days with lower uncertainty and meteorological variability (21 June and 8 July) indicate that the differences in ET may be based in plant stomatal control.

Comparison of HRMET with Eddy Covariance
Across the season, the HRMET ET values fall within the range of the EC ET values ( Figure 10). The mean and median of the HRMET values performed closest to the EC ET when compared to those minimum and maximum HRMET values from the scene. Mean HRMET values tended to have higher ET values than the EC tower by a mean of 0.27, a median of 0.25, a maximum of 0.72, and a minimum of <0.00 mm h −1 across the season (Figures 11 and 12). The 30 m flight on 21 June had the closest HRMET and EC values with a difference of <0.01 mm h −1 . The 30 m flight on 8 July had the highest difference in HRMET and EC values with a difference of 0.71 mm h −1 . There were no observable differences in how ET maps from the different mission elevations compared to the EC data. The mean and median HRMET values performed similarly with an RMSE of 0.21 mm h −1 (Figure 12, Table 4). The minimum HRMET values had the weakest relationship (R 2 of 0.05) and an RMSE of 0.29 mm h −1 (Figure 12, Table 3). The maximum HRMET values had the strongest relationship (R 2 of 0.48) and an RMSE 0.32 mm h −1 (Figure 12, Table 4).    (Figure 12, Table 4). The minimum HRMET values had the weakest relationship (R 2 of 0.05) and an RMSE of 0.29 mm h −1 (Figure 12, Table 3). The maximum HRMET values had the strongest relationship (R 2 of 0.48) and an RMSE 0.32 mm h −1 (Figure 12, Table 4).   The 21 June and 8 July missions followed irrigation events (~10 and~13 mm, respectively, Figure 3). The 6 August mission had heavy rain the day before the flights (~25 mm).   The 21 June and 8 July missions followed irrigation events (~10 and ~13 mm, respectively, Figure 3). The 6 August mission had heavy rain the day before the flights (~25 mm). One-to-one line added as reference to the plots. Slope, y-intercept, R 2 and RMSE were calculated for the mean, median, minimum, and maximum HRMET ET values.

Comparison of HRMET and Eddy Covariance ET Estimates
The mean and median HRMET ET values in this study were generally higher than the EC ET values, which is consistent with past work that has found EC ET values to be lower than remotely sensed ET values from the pySEBAL, SEBs, and TSEB models [12,19]. However, our findings differ from other findings where the EC values were higher than the METRIC, MODIS, and Gleam ET values [19,30,34]. The pySEBAL, and SEBS models tend to underestimate actual ET as a result of not having true 'cold' pixels in the model domain [19]. The lack of 'cold' pixels would result in an underestimation of ET because the 'cold' pixel sets the upper limit. In contrast, HRMET does not use a 'hot' and 'cold' pixel approach, so the reason for overestimation is not related to the absence of true well-watered pixels in the scene. Alternatively, this study's assumption that the EC system's area of influence (footprint) is a static circle has also been identified as a source of inconsistency in studies comparing remotely sensed ET maps to EC measurements [48]. The EC system's area of influence changes with wind direction and atmospheric stability, which would change the appropriate HRMET domain to compare with EC measurements for a given day [48]. Future work could compute individual EC footprints for each mission set to improve scaling for comparisons. Sometimes EC measurements display a lack of closure of the energy budget which could indicate an underestimation bias of 10% that contributed to difference between HRMET and EC values [49].
The input data for HRMET and the EC tower have a difference in ideal climatic conditions. Namely, EC systems perform better with greater wind speeds while the surface temperature imagery collected for HRMET increases in variability with increasing wind and gust speeds. For example, 8 July had a lower wind speed (0.8 m s −1 ) and the EC ET was 0.71 mm h −1 lower than the HRMET values. For the EC time series on 8 July, there is an observable reduction in ET that coincides with a reduction in wind speed.
Though we assumed that there is at least~10% uncertainty in EC estimates because of turbulent sampling or energy budget closure issues [47], we did not quantify the uncertainty associated with EC estimates used for HRMET comparison, nor the uncertainties inherent in the HRMET model assumptions and equations. Additionally, the scales at which the two are collected could have an impact on how outliers impact the overall value. The EC tower gives a single value while HRMET multiple values for the same area. We compared mean, median, max, and min HRMET values to the EC ET values, but there may be other spatially weighted approaches, such as a footprint cutoff approach [50], that could be more appropriate for these comparisons. Finally, despite our attempts to match HRMET and EC estimates in time, the difference in capture period (30-min for EC and variable for dronebased image capture) makes it difficult to associate the ET values to the right moment in time. This can be particularly challenging for comparison because of dynamic atmospheric conditions. Nevertheless, our overall results show consistent temporal patterns across the two approaches, confirming the suitability of drone for ET mapping in the Midwest United States.

Adapting the HRMET Model for Very High-Resolution Imagery
HRMET was initially developed for use with aerial imagery (1-30 m 2 scale), and because of this, there are issues of scale when applying it to drone imagery (0.2-0.6 m 2 ). The resolution for the input multispectral imagery is too fine to calculate a vegetation index (VI) needed to use the regression to estimate LAI and crop height, because VIs are canopy scale measurements while the drone imagery is on the leaf scale. The imagery can be aggregated to calculate the VI and then disaggregated to maintain resolution, but the disconnect in principle must be noted. A possible solution to this would be directly measure LAI and crop height by using relationships with RGB reflectance values or using point cloud estimations, respectively [51,52].
Another issue with the input VI raster was the false presence of LAI and crop height in bare soil spots (as can be seen in the field road of Figures 3 and 4). This stems from empirical model extrapolation at very low values. This results in an overestimation of ET in bare soil areas by estimating transpiration for pixels where there is no vegetation present. A possible solution to this issue would be to use two source energy balance model that would separately estimate soil evaporation, though this may increase computational needs [53].

How High to Fly?
Overall, there were no observable differences or trends between the EVI, height, temperature, and ET estimates from the different mission elevations. There is a marginal increase in resolution and detectable features with the lower elevation flights. The influence of these additional detectable features is present at all elevations. At higher elevations, the additional features become incorporated in larger mixed pixels. Further study is needed to determine the potential for using the additional detectable features as a tool for creating management zones [54]. However, low elevation flights can create challenges with photo alignment when using multi-lens cameras [55]. Lower elevations also take longer to fly, reducing the coverage area as well as introducing the opportunity for solar radiative changes. Low elevation images also have smaller ground area to capture unique tie-points which can cause issues with aligning the imagery. In this case, there was not a need to have imagery resolution beyond what can be managed computationally. Additionally, precision irrigation can only be done as efficiently as the irrigation infrastructure in place. Based on these challenges associated with lower-elevation flights, we suggest that the flight elevation be no lower than the elevation necessary to capture the finest characteristics of the field upon which management is actionable.
Higher resolution ET maps could be more useful for perennial crops where there is high individual crop value (e.g., orchards, vineyards). The higher resolution would allow for better monitoring of water stress, health, and pests in individual trees, bushes, or vines [55][56][57]. There may be more considerations for tall orchard crops where the aspect may have a greater impact on photogrammetry and would emphasize the need to study specific mission planning to optimize image capture.

When to Fly?
When collecting drone imagery, it is best to follow current best practices and collect imagery around solar noon, fly the same time for every flight, use at least 70% front and side overlap, fly under 12 m s −1 winds, avoid clouds, and avoid collecting imagery while the field is being actively irrigated [58,59]. In practice, it is hard to meet these conditions while attempting to collect imagery on a frequent schedule. To accommodate for these demands, we collected imagery with the best conditions that were present at the time. As such, there were some conditions that were not ideal. Three missions of note were 21 June, 8 July, and 6 August where the mission domain had either just been irrigated or experienced a heavy rain event. This would cause an increase in the evaporation portion of ET and higher ET values [60,61]. However, these three days did not have the same response; 8 July and 6 August had relatively high ET values (range mm h −1 ) while 21 June had relatively low ET values (range mm h −1 ). This could be because the 21 June mission set was conducted early in the day (9:30 am) where the evaporative demand was not yet at a maximum value. Alternatively, the 8 July and 6 August mission sets were collected near solar noon at peak evaporation demand.

Limitations
Some of the conditions in which the imagery was collected resulted in the introduction of uncertainty. Discrepancies between the thermal imagery and the ground reference data (Appendix A) could indicate atmospheric interference [62]. In particular, the wind speed and variation could have temporarily cooled the surface at a given location during the capture period. The RMSE values for this imagery were higher but comparable to Simpson et al., 2021 (~3 • C); however, their data were collected over a large, controlled grass area that avoided the 'spot size effect' which led to more stable and precise readings.
Thermal ground reference data were collected after the final flight in the mission set. This led to a 30-min latency between the start of the first flight and the first ground reference point taken. Depending on how stable the surface temperature was during the sampling time, this could have introduced uncertainty. Thermal ground references of maximum and minimum temperatures were not collected, which could have set the scale for the surface temperature range [28]. The thermal ground references were also collected as spot measurements but were compared to the mean of a 1.5 m diameter circle from the thermal imagery. For future work, we recommend a larger thermal ground reference area (>3 m 2 ) that would appear in several pixels of thermal imagery collected at 90-m. The June 21st and 25th flights were flown earlier than recommended, potentially introducing more variability in solar radiation [29]. Clouds were not always avoidable (22 July) causing variability in solar radiation across the elevations and in individual plots. This variability can cause 'dark' low reflectance spots on the imagery under cloudy conditions [30].
There are more causes of spatial variability of ET to address. The variability we observed cannot be fully attributed to stomatal control. For future work, we recommend using spatially-explicit measurements of stomatal conductance and photosynthesis to complement tower-based comparisons [63].

Conclusions
Drone-based remote sensing offers the ability to rapidly generate ET maps within a season that can be used to make in-season decisions. The use of ETmaps with precision irrigation has the potential to reduce water use and nitrogen loss from agricultural fields. The High Resolution Mapping of Evapotranspiration (HRMET) model has been previously used for mapping ET from Midwestern potatoes using aerial imagery. We found that HRMET may also be useful for ET mapping with drone imagery if the flights are collected under optimal conditions. We evaluated whether there was an increase in insightful information from ET maps created using imagery collected at 90, 60, or 30 m. Though there were some additional detectable features observed from low elevation flights, the tradeoff in resources and computation does not currently justify low elevation flights for annual vegetable crops in the Midwest USA. Spatial data were the predominant source of uncertainty in the HRMET model. High solar radiation variability or low wind speeds were the biggest driver in the difference between HRMET and Eddy Covariance (EC) ET estimates. Uncertainty maps are an underutilized tool that have the potential to evaluate the strengths and weaknesses of high-resolution ET maps. Going forward, further application to the HRMET model, or models like it, could prove beneficial, especially for precision irrigation crops that may not have the ideal 'hot' and 'cold' pixels. Further study is needed to understand how HRMET performs in other climates outside the Midwest United States.

Appendix A
tower. We would also like to acknowledge K. Lawton, J. Oestreich, and A. Doherty for computational support at UW-Stevens Point. We would like to thank A. Mandel and D. Young for photogrammetry support at UC-Davis.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. Appendix A Figure A1. Thermal ground reference data vs. surface temperature values remotely collected with Micasense Altum. Ground reference data were collected within half an hour of mission launch. Table A1. Root Mean Square Error (RMSE) of thermal ground reference data vs. surface temperature (°C) values collected with Micasense Altum. Ground reference data were collected within half an hour of mission launch. Figure A1. Thermal ground reference data vs. surface temperature values remotely collected with Micasense Altum. Ground reference data were collected within half an hour of mission launch.