Mapping Maize Evapotranspiration with Two-Source Land Surface Energy Balance Approaches and Multiscale Remote Sensing Imagery Pixel Sizes: Accuracy Determination toward a Sustainable Irrigated Agriculture

Abstract

This study evaluated the performance of remote sensing (RS) algorithms for the estimation of actual maize evapotranspiration (ET_a) using different spaceborne, airborne, and proximal multispectral data in a semi-arid climate region to identify the optimal platform that provides the best ET_a estimates to improve irrigation water management and help make irrigated agriculture sustainable. The RS platforms used in the study included Landsat-8 (30 m pixel spatial resolution), Sentinel-2 (10 m), Planet CubeSat (3 m), multispectral radiometer or MSR (1 m), and a small uncrewed aerial system or sUAS (0.03 m). Two-source surface energy balance (TSEB) models, implementing the series and parallel surface resistance approaches, were used in this study to estimate hourly maize ET_a. The data used in this study were obtained from two maize research sites in Greeley and Fort Collins, CO, USA, in 2020 and 2021. Each research site had different irrigation systems. The Greeley site had a subsurface drip system, while the Fort Collins site had surface irrigation (furrow). Maize ET_a predictions were compared to observed maize ET_a data from an eddy covariance system installed at each research site. Results indicated that the MSR5 proximal platform (1 m) provided optimal RS data for the TSEB algorithms. The MSR5 “point-based” nadir-looking surface reflectance data and surface radiometric temperature combination resulted in the smallest error when predicting hourly (mm/h) maize ET_a. The mean bias and root mean square errors (MBE and RMSE, respectively), when predicting maize hourly ET_a using the MSR5 sensor data, were equal to −0.02 (−3%) ± 0.07 (11%) mm/h MBE ± RMSE and −0.02 (−3%) ± 0.09 (14%) mm/h for the TSEB parallel and series approaches, respectively. The poorest performance, when predicting hourly TSEB maize ET_a, was from Landsat-8 (30 m) multispectral data combined with its original thermal data, since the errors were −0.03 (−5%) ± 0.16 (29%) mm/h and −0.07 (−13%) ± 0.15 (29%) mm/h for the TSEB parallel and series approaches, respectively. These results indicate the need to develop methods to improve the quality of the RS data from sub-optimal platforms/sensors/scales/calibration to further advance sustainable irrigation water management.

1. Introduction

To improve irrigation water management in agricultural fields and attain sustainability, it is critical to define the optimal irrigation time and irrigation amounts to replenish the soil vadose layers, where crop roots develop, to conserve water and soil resources. Irrigation water management practices are often based on the soil water balance (SWB) approach for irrigation scheduling development [1,2]. The SWB approach for irrigation provides a soil water volume balance that accounts for the inflow and outflow of water fluxes in the crop root zone to define the temporal changes in the soil volumetric water content [3]. The simplified daily SWB approach is given by Equation (1):

$$D_{r,i}=D_{r,i-1 }-{\left(P-RO\right)}_i-I_i-C{R}_i+E{T}_{c,i}+D{P}_i$$

(1)

where D_r,i is the water depleted in the root zone at the end of day i^th; D_r,i−1 is the water in the root zone in the previous day (i − 1)^th; P_i is the rainfall water depth; RO_i is the surface water runoff; I_i is the net irrigation water depth; CR_i is the capillary rise from shallow water table (groundwater); ET_c,i is the daily crop evapotranspiration; DP_i is the deep percolation (vertical water loss beyond the root zone). All variables in Equation (1) are given as water depth units (e.g., mm or in).

Better irrigation strategies are often related to an SWB approach that accounts for an accurate ET_c estimate through the actual crop evapotranspiration (ET_a) rates determination. This is the correct amount of water depleted in the soil throughout the plant root zone that would be replenished with irrigation when ET_c (through ET_a) has been properly determined. In general, throughout this process, water and nutrient savings or conservation are achieved because common irrigation practices tend to over-irrigate, promoting water, soil, and agro-chemical losses through land surface runoff and deep percolation, potentially contaminating groundwater and/or surface water bodies. At the local farm scale, accurate crop ET_a estimation is critical to support quasi-real-time decision-making approaches for water allocation and optimization of irrigation water management [4,5].

Modeling advancements in remote sensing (RS) of the environment have facilitated the quasi-real-time mapping of crop water requirements or ET_a for irrigation on a spatio-temporal basis, using multispectral and thermal imagery from different sensor types [6,7] since the early 1970s. Remote sensing involves the scientific measurement of emitted and reflected light across various spectral ranges, including visible, invisible, and longwave infrared (LWIR), without direct contact with the target area [8]. Optical devices mounted on aerial platforms (e.g., small aircraft or automated aerial vehicles), spaceborne systems (e.g., satellites), and proximal instruments (e.g., handheld roaming or stationary radiometers) have generated data at different temporal, spectral, and spatial resolutions, benefiting applications like irrigation water management, soil nutrient monitoring, crop growth assessment, and yield mapping [9,10,11]. The use of RS techniques to support sustainability of irrigation scheduling practices has been investigated for more than 50 years [12,13].

Remote sensing of crop ET_a approaches that use multispectral and thermal data to map crop ET_a are fundamentally based on the land surface energy balance (SEB) concept. The SEB approach for estimating ET_a calculates the energy required for evaporating water (latent heat flux, LE) as the residual term of the simplified SEB (Equation (2)).

$$\mathrm{L}\mathrm{E}={\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\mathrm{H}$$

(2)

where LE is the latent heat flux; R_n is the net radiation flux; G is the soil heat flux; and H is the sensible heat flux. All terms in Equation (2) are given in W/m². The SEB LE flux is then converted to instantaneous crop ET_a (e.g., mm/h) during the RS sensor overpass. There are two common methods to determine ET_a using the SEB approach: (a) the one-source SEB (henceforth, OSEB), which considers the combined contributions of soil and vegetation to ET_a rates [14,15,16,17], and (b) the two-source SEB (or TSEB) that partitions heat fluxes and the crop ET_a in a component related to the water transpired by the plants and another related to the evaporated water from the soil [18,19,20,21].

The TSEB is a robust SEB approach suitable for estimating spatial ET_a that was initially developed by [21]. The TSEB model has two different approaches for estimating the H flux in Equation (2): the parallel surface resistances TSEB (henceforth, TSEB_par) and the (in) series surface resistances TSEB (henceforth, TSEB_ser). The TSEB_par model considers the processes of heat transfer among plants, soil, and the air above the canopy as independent of each other with two surface resistances for heat transfer. The TSEB_ser method includes the concept of heat transfer interconnection in the soil–plant–atmosphere continuum through an additional surface resistance term and a parametrization of the aerodynamic surface temperature (T_o) as a weighted-average temperature among soil, plant, and air temperatures with respective resistances as weights. Typical ET_a estimation errors, when using TSEB_par or TSEB_ser, were reported to be within 7% to 25% for row crops [22,23,24]. In regard to the desired frequency of ET_a estimates, the study by [25] indicated that a four-day RS platform overpass frequency (of usable data) would be the minimum needed for current interpolation techniques to yield meaningful daily ET_a estimates between acquired RS data. However, with high RS data acquisition frequencies, more reliable and accurate daily ET_a estimations will be possible. Therefore, more timely and accurate irrigation water amounts would be delivered to surface and pressurized systems if accurate daily ET_a maps were produced.

Examining various RS platforms that offer multispectral images of cropland fields at diverse spectral and spatial resolutions is crucial for assessing the reliability of different ET_a prediction algorithms and their accuracy when predicting ET_a values in time and space [26]. Furthermore, accurate estimation of crop ET_a, when used to optimize the irrigation water amounts and timing of application, advances environmental sustainability by decreasing topsoil erosion in agricultural areas due to reduced field surface runoff and conserves water and soil nutrients within agricultural districts, protecting the environment by reducing groundwater withdrawn rates, maintaining ecological water table levels, and preserving adequate water quality of both aquifers and surface water bodies (e.g., lakes, artificial reservoirs, and rivers). However, there have been very few studies attempting to address the performance of the TSEB RS of ET_a algorithms across different spectral and spatial scales. In a recent study, Ref. [27] explored the accuracy of the TSEB model developed by [21] using different small uncrewed aerial system (sUAS or drone) imagery pixel sizes, ranging from 0.10 m to 0.60 m, in a vineyard field located in California. The drone-captured images were subsequently aggregated to produce lower-resolution imagery with pixel sizes spanning from 3.6 m to 30 m. The results from [27] demonstrated that errors in R_n and G were relatively consistent across various RS resolutions. In contrast, errors in H and LE fluxes exhibited a clear relationship with the spatial resolution of the RS data. Another study by [28] investigated the effect of pixel heterogeneity for tree–grass when predicting ET_a using hyperspectral airborne imagery (1.5 m to 1000 m spatial resolution) and Sentinel imagery products at 20 m and 1000 m using a TSEB RS algorithm in central Spain. They found that large uncertainty, when estimating ET_a, occurred for coarse spatial resolutions.

Even though these studies have contributed to science, there have not been comprehensive studies that evaluate the differences in accuracy of the TSEB RS of ET_a algorithms using multispectral images from multiscale RS platforms such as those from proximal, airborne, and spaceborne sensors. The published studies focused only on a few RS sensors or platforms, often resampling (upscaling) their images to generate different pixel spatial resolutions. Therefore, in this study, it is hypothesized that, depending on the source of a given RS image (e.g., spaceborne, airborne, proximal platforms, sensor type, and imagery post-processing corrections), the accuracy of ET_a mapping products will vary for a given RS of the ET_a algorithm. If the stated hypothesis is valid, determining the optimal RS spectral and spatial resolution becomes necessary (critical) to better sustain irrigated agriculture by improving the estimation of ET_a when sub-optimal RS platforms (data) are used with a given RS of the ET_a algorithm.

Therefore, the objectives of the study were to (a) assess the impact (errors) on hourly ET_a estimation associated with the use of different spectral and spatial resolution RS data from multispectral spaceborne, airborne, and proximal RS sensors and when using two different TSEB algorithms, and (b) identify the RS spectral and spatial data (resolution) that provides the most accurate TSEB-based maize ET_a predictions for a specific algorithm.

2. Materials and Methods

2.1. Description of the Research Sites

2.1.1. Limited Irrigation Research Farm (LIRF)

The Limited Irrigation Research Farm (LIRF) is located in Greeley, Colorado (CO), USA, and is under the management of the United States Department of Agriculture—Agricultural Research Service (USDA—ARS). The farm is geographically located at a latitude of 40.4463° N, a longitude of 104.6371° W, and an elevation of 1432 m above mean sea level (ASL). The study involved two adjacent rectangular maize fields, each measuring 190 m by 110 m (Figure 1), where field data were collected during the periods between July and September of 2020 and 2021.

For each respective field crop growth season, each field was subjected to different irrigation water management strategies. In 2020, the West Field, hereafter referred to as Field W, was fully irrigated. In this context, “fully irrigated” represents the conditions in which frequent irrigation events were scheduled to maintain soil water content in the crop (maize) soil root zone at non-water-stress levels. Conversely, the East Field, designated as Field E, was managed as a deficit-irrigated field, resulting in crop/soil water stress conditions throughout the growing season. In 2021, the irrigation water management practices were switched between these treatment plots. Field W was transformed into the deficit-irrigated field, while Field E became the fully irrigated plot. A summary of soil wetting events, which includes irrigation and rainfall, for the years 2020 and 2021, is provided in

Table 1. Cumulative soil wetting events (irrigation and rainfall) at LIRF Fields W and E in 2020 and 2021.

Irrigation Scheduling	Growing Season	Research Field	Irrigation Events	Cumulative Gross Irrigation (mm)	Cumulative Rainfall (mm)
Fully Irrigated	2020	W	5	472	36
	2021	E	4	330	98
Deficit Irrigated	2020	E	3	309	36
	2021	W	2	176	98

.

Each maize field had the same irrigation system, a subsurface drip irrigation setup with laterals (pipes) buried at a depth of 0.23 m and emitters spaced every 0.30 m. The maize rows were north–south and were spaced 0.76 m apart. The distance between adjacent maize plants was 0.17 m. The planting density for maize remained consistent at 87,500 plants per hectare during both years. In 2020, the selected maize variety was drought tolerant NK9227-5222A (Syngenta Inc., Basel, Switzerland). Planting took place on 6 May 2020, and the harvest occurred on 13 and 14 October 2020. In 2021, there was a change in maize varieties with the introduction of other drought-tolerant maize options, P9998Q and P0157AMXT (Pioneer Hi-Bred International, Inc., Johnston, IA, USA), along with CH 194-49 DG (Channel Bio Corporation, Saint Louis, MO, USA). Planting for this year was carried out on 13 May 2021, and the harvest took place on 11 and 12 October 2021. Fields W and E, each comprising approximately 83% of their respective plots, were planted with the maize variety P0157AMXT, as indicated in Figure 2. The experimental design and data collection stations at LIRF for the years 2020 and 2021 are shown in Figure 3. It is worth noting that the prevailing wind direction remained consistent: the wind was from the south (S) to southeast (SE) direction during both years of data collection.

2.1.2. Irrigation Innovation Consortium (IIC)

This study included data collected from the Colorado State University IIC site during 2020 and 2021. Two maize fields were selected as the primary locations for data collection (Figure 4). This site is located in Fort Collins, CO, USA, at a latitude of 40.5542° N, a longitude of 105.0038° W, and an elevation of 1486 m ASL. It has a local climate characterized as a subtropical steppe with cold semiarid tendencies.

The data collection took place in two surface-irrigated (furrow) maize fields during July to September in both 2020 and 2021, as illustrated in Figure 4a. These two fields, designated as Fields F and D, had respective surface areas of 64,750 m² and 74,867 m². The crop-row orientation in Field F was east–west, while Field D had rows oriented in a north–southeast direction, with rows spaced approximately 0.17 m apart. The soil in both fields had a consistent sandy loam texture throughout the entire maize root zone system, with measured volumetric water content at field capacity (VWC_FC), permanent wilting point (VWC_PWP), and saturation (VWC_SAT) of 0.189, 0.069, and 0.410 m³/m³, respectively. The surface irrigation system employed 40 mm diameter aluminum siphon tubes and delivered water to the fields from the main on-site irrigation canal located between the two research fields. The choice of maize varieties differed between 2020 and 2021. In 2020, G02K39-3120 (Golden Harvest, Minnetonka, MN, USA) was planted on May 13, at an approximate rate of 8 seeds per m². For 2021, the NK0243-3120 and NK0314-5122 varieties (Syngenta AG, Basel, Switzerland) were planted in both fields, as indicated in Figure 4b. The seeding date was 13 May 2021, and the planting rate was 8 seeds per m².

The direction of the irrigation water flow was from east to west in Field F. In the case of Field D, the irrigation water flowed in the furrows from north to southwest. The irrigation events usually took place two to three days after obtaining water from the Sand Dike Lateral Company in Fort Collins, CO, USA. Each irrigation event extended over a duration of 6 to 12 h. The irrigation scheduling was determined based on the FAO-56 methodology [2], which was an integral component of the Water Irrigation Scheduler for Efficient (WISE) Application [1]. The WISE application uses a water balance approach for determining irrigation amounts and timing based on a dual crop coefficient approach for ET_a estimation from tabularized crop coefficient (K_c) values, a water stress coefficient (K_s), and daily alfalfa reference ET (ET_r) rates.

The K_s value was set to 1 since the maize field was constantly irrigated during the data collection season, and root zone water depletion was assumed to have been kept below allowable water depletion levels.

Throughout both years of data collection, the prevailing wind direction was consistently from the southeast (SE) to the south (S). A summary of the data related to soil wetting events, encompassing both irrigation and rainfall events for the years 2020 and 2021, is given in

Table 2. Cumulative soil wetting events (irrigation and rainfall) at IIC Fields F and D in 2020 and 2021.

Maize Growing Season	Research Field	Number of Irrigation Events	Cumulative Gross Irrigation (mm)	Cumulative Rainfall (mm)
2020	F	7	1620	34
	D	5	870
2021	F	8	1081	104
	D	4	1064

.

2.2. Crop Evapotranspiration Algorithm

2.2.1. Two-Source Surface Energy Balance (TSEB)

The TSEB estimates LE as follows (Equations (3)–(5)):

$$\mathrm{L}\mathrm{E}=\mathrm{L}{\mathrm{E}}_{\mathrm{c}}+\mathrm{L}{\mathrm{E}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$$

(3)

$$\mathrm{L}{\mathrm{E}}_{\mathrm{c}}=\mathrm{R}{\mathrm{n}}_{\mathrm{c}}-{\mathrm{H}}_{\mathrm{c}}$$

(4)

$$\mathrm{L}{\mathrm{E}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}=\mathrm{R}{\mathrm{n}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}-\mathrm{G}-{\mathrm{H}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$$

(5)

where the subscripts “c” and “soil” refer to the “canopy” and “bare soil” conditions, respectively. All terms in Equations (3)–(5) have W/m² units.

The LE is converted from W/m² to hourly ET_a (mm/h) using Equation (6):

$$\mathrm{E}{\mathrm{T}}_{\mathrm{a}\mathrm{h}}={\mathrm{C}}_{\mathrm{f}}\times {\mathrm{d}}_{\mathrm{f}}\times \frac{\mathrm{L}\mathrm{E}}{{\mathsf{\lambda }}_{\mathrm{v}}{\mathsf{\rho }}_{\mathrm{a}}}$$

(6)

where ET_ah is the hourly instantaneous actual crop ET_a (mm/h); d_f is a time-scale factor in seconds (e.g., 1 h = 3600 s); λ_v is the latent heat of water vaporization (J/kg); ρ_a is the air density (kg/m³); C_f is a length scale conversion factor (e.g., 1000 mm/1 m).

The G model used was the soil heat flux (Equation (7)) approach recommended by [21], which is a fraction of the net radiation from the soil layer (Rn_soil) during midday conditions coinciding with most spaceborne RS platform overpass times.

$$\mathrm{G}=0.35\mathrm{R}{\mathrm{n}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$$

(7)

The Rn_c and Rn_soil terms (appearing later in Equations (18) and (19)) are calculated following the approaches by [19,21] and indicated by Equations (8) and (9), respectively:

$$\mathrm{R}{\mathrm{n}}_{\mathrm{c}}=\left(1-{\mathsf{\Omega }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{a}\mathrm{r}}\right)\left(1-{\mathsf{\alpha }}_{\mathrm{c}}\right){\mathrm{R}}_{\mathrm{s}}+\exp \left(-0.95\mathrm{L}\mathrm{A}\mathrm{I}\right){\mathrm{L}}_{\mathrm{s}\mathrm{k}\mathrm{y}}+\left[1-\exp \left(-0.95\mathrm{L}\mathrm{A}\mathrm{I}\right)\right]{\mathsf{ϵ}}_{\mathrm{c}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{c}\mathrm{K}}^4-{\mathsf{ϵ}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{K}}^4$$

(8)

$$\mathrm{R}{\mathrm{n}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}={\mathsf{\Omega }}_{\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{a}\mathrm{r}}\left(1-{\mathsf{\alpha }}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}\right){\mathrm{R}}_{\mathrm{s}}+\left[1-\exp \left(-0.95\mathrm{L}\mathrm{A}\mathrm{I}\right)\right][{\mathsf{ϵ}}_{\mathrm{a}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{a}\mathrm{K}}^4+{\mathsf{ϵ}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{K}}^4-2{\mathsf{ϵ}}_{\mathrm{c}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{c}\mathrm{K}}^4]$$

(9)

where Ω_solar is the solar atmospheric transmittance (dimensionless); α_c and α_soil are the canopy and bare soil albedo (dimensionless), respectively; ϵ_c is the canopy thermal emissivity (dimensionless); ϵ_soil is the bare soil thermal emissivity (dimensionless); ϵ_a is the air emissivity (dimensionless); T_cK is the canopy temperature (K); T_soilK is the soil temperature (K); T_aK is the air temperature (K); L_sky is the longwave radiation from the sky (W/m²); σ is the Stefan–Boltzmann constant (5.67 × 10⁻⁸ W/m²/K⁴); LAI is the leaf area index (m²/m²); R_s is the incoming shortwave solar radiation (W/m²).

The α_c and α_soil were calculated following the approaches by [29] and [30], respectively. The study by [29] gave an exponential model that relates maize albedo and green LAI for a semi-arid climate region. The study by [30] provided a multivariate model that relates α_soil and the visible surface reflectance bands of the light spectrum (RED, GREEN, and BLUE). The Ω_solar variable was calculated using the nonlinear model from [31], in which the Ω_solar is a function of LAI, surface absorptivity factor [32], light extinction coefficient, and the fraction of incident photosynthetically active radiation or PAR [33,34,35].

The T_cK and T_soilK are calculated through an iterative approach using Equation (10) below [21]:

$${\mathrm{T}}_{\mathrm{s}\mathrm{K}}\cong {\left[{\mathrm{f}}_{\mathrm{c}}\left(\mathsf{\phi }\right){\mathrm{T}}_{\mathrm{c}\mathrm{K}}^4+(1-{\mathrm{f}}_{\mathrm{c}}\left(\mathsf{\phi }\right){\mathrm{T}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}\mathrm{K}}^4)\right]}^{0.25}$$

(10)

where T_sK is the equivalent nadir-looking surface temperature calculated from plant and soil temperature composites (K); f_c(φ) is the fractional green vegetation cover for a radiometric field-of-view angle φ [21].

The TSEB Series Approach for H Estimation

In the TSEB series algorithm (TSEB_ser), the heat exchange between soil, plant, and air is assumed to be interconnected [21]. The expressions for calculating H using the TSEB_ser approach are given by Equations (11)–(13) below [21]:

$${\mathrm{H}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}={\mathrm{H}}_{\mathrm{c}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}+{\mathrm{H}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}$$

(11)

$${\mathrm{H}}_{\mathrm{c}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}={\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}\mathrm{a}}\left(\frac{{\mathrm{T}}_{\mathrm{c}}-{\mathrm{T}}_{\mathrm{o}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}}{{\mathrm{r}}_{\mathrm{x}}}\right)$$

(12)

$${\mathrm{H}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}={\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}\mathrm{a}}\left(\frac{{\mathrm{T}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}-{\mathrm{T}}_{\mathrm{o}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}}{{\mathrm{r}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}}\right)$$

(13)

where T_c, T_soil, and T_o^[ser] are the canopy, bare soil, and TSEB_ser aerodynamic temperatures (K), respectively; r_x and r_soil are the total boundary layer resistance of canopy leaves (s/m) and the soil resistance (s/m), respectively. The superscript [ser] refers to the TSEB_ser (series) algorithm. More details on how to calculate these resistance terms and other auxiliary variable inputs for the TSEB can be found in Appendix A and Appendix B.

The air parameters ρ_a and C_pa are calculated as indicated by Equations (14) and (15), respectively (as used in [36]):

$${\rho }_a=\left(\frac{P}{R_d{T}_a}\right)\left(1-\frac{0.378e_a}{P}\right)$$

(14)

$$C_{pa}=1004.7\times \left(1+\frac{0.522e_a}{P}\right)$$

(15)

where P is the local atmospheric pressure (Pa); e_a is actual vapor pressure (Pa); R_d is the gas constant for dry air (≈287.04 J/kg/K). The air parameters ρ_a and C_pa are given in kg/m³ and J/kg/K units, respectively.

The T_o^[ser] (given in K) is calculated as indicated by Equation (16) [21]:

$${\mathrm{T}}_{\mathrm{o}}^{\left[\mathrm{s}\mathrm{e}\mathrm{r}\right]}=\frac{\raisebox{1ex}{${\mathrm{T}}_{\mathrm{a}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{a}\mathrm{h}}$}\right.+\raisebox{1ex}{${\mathrm{T}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$}\right.\raisebox{1ex}{${\mathrm{T}}_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{x}}$}\right.}{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{a}\mathrm{h}}$}\right.+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}$}\right.+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{x}}$}\right.}$$

(16)

In this study, we followed the improved TSEB_ser approach from [18] and used the modified Penman–Monteith (PM) approach instead of the Priestley–Taylor (PT) modified model (as described in [21]) to calculate an initial T_c value to derive the surface temperature composites (T_c and T_soil). The modified PM approach for an initial T_c value is indicated by Equation (17):

$${\mathrm{T}}_{\mathrm{c}\mathrm{O}}={\mathrm{T}}_{\mathrm{a}}+\frac{\mathrm{R}{\mathrm{n}}_{\mathrm{c}}{\mathrm{r}}_{\mathrm{a}\mathrm{h}}\mathsf{\gamma }(1+\raisebox{1ex}{${\mathrm{r}}_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{a}\mathrm{h}}$}\right.)}{{\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}\mathrm{a}}\left[\mathsf{\Delta }+\mathsf{\gamma }\left(1+\raisebox{1ex}{${\mathrm{r}}_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{a}\mathrm{h}}$}\right.\right)\right]}-\frac{{\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}}}{\mathsf{\Delta }+\mathsf{\gamma }(1+\raisebox{1ex}{${\mathrm{r}}_{\mathrm{c}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{r}}_{\mathrm{a}\mathrm{h}}$}\right.)}$$

(17)

where T_cO is the initial guess for canopy temperature (K); γ is the psychometric constant (kPa/°C); r_c is the bulk canopy resistance (s/m); Δ is the slope of the saturation vapor pressure curve (kPa/°C); e_s and e_a are the saturated and actual vapor pressures in kPa, respectively. The calculation of r_c is described in Appendix B.

The TSEB Parallel Approach for H Estimation

The TSEB parallel (TSEB_par) assumes that the processes that derive the heat transfer among plants, soil, and the air above are independent and can be modeled with two separate resistances for heat transfer [18,21]. The expressions for calculating H using the TSEB_par approach are given by Equations (18)–(20) below [21]:

$${\mathrm{H}}^{\left[\mathrm{p}\mathrm{a}\mathrm{r}\right]}={\mathrm{H}}_{\mathrm{c}}^{\left[\mathrm{p}\mathrm{a}\mathrm{r}\right]}+{\mathrm{H}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}^{\left[\mathrm{p}\mathrm{a}\mathrm{r}\right]}$$

(18)

$${\mathrm{H}}_{\mathrm{c}}^{\left[\mathrm{p}\mathrm{a}\mathrm{r}\right]}={\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}\mathrm{a}}\left(\frac{{\mathrm{T}}_{\mathrm{c}}-{\mathrm{T}}_{\mathrm{a}}}{{\mathrm{r}}_{\mathrm{a}\mathrm{h}}}\right)$$

(19)

$${\mathrm{H}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}^{\left[\mathrm{p}\mathrm{a}\mathrm{r}\right]}={\mathsf{\rho }}_{\mathrm{a}}{\mathrm{C}}_{\mathrm{p}\mathrm{a}}\left(\frac{{\mathrm{T}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}-{\mathrm{T}}_{\mathrm{a}}}{{{\mathrm{r}}_{\mathrm{a}\mathrm{h}}+\mathrm{r}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}}\right)$$

(20)

where the superscript [par] refers to the TSEB_par (parallel) algorithm.

For the TSEB_par algorithm, the initial assumption regarding the initial value for LE_c was based on the original work from [21] and indicated by Equation (21), as follows:

$$\mathrm{L}{\mathrm{E}}_{\mathrm{c}\mathrm{i}}=\left[1.3{\mathrm{f}}_{\mathrm{g}}\left(\frac{\mathsf{\Delta }}{\mathsf{\Delta }+\mathsf{\gamma }}\right)\right]\left[\mathrm{R}{\mathrm{n}}_{\mathrm{c}}+\mathrm{R}{\mathrm{n}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}-\left(\mathrm{R}{\mathrm{n}}_{\mathrm{c}}+\mathrm{R}{\mathrm{n}}_{\mathrm{s}\mathrm{o}\mathrm{i}\mathrm{l}}\right)\times \mathrm{e}\mathrm{x}\mathrm{p}(0.90\ln (1-{\mathrm{f}}_{\mathrm{c}}\left)\right)\right]$$

(21)

where LE_ci is the initial guess value for iterating T_c and T_soil (W/m²), and f_g is the green fraction of LAI (dimensionless), calculated as indicated by [21].

2.3. Vegetation Indices Calculation

The normalized difference vegetation index (NDVI) and optimized soil-adjusted vegetation index (OSAVI) [37] were calculated by Equations (22) and (23), respectively:

$$\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}\mathrm{D}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}\mathrm{D}}$$

(22)

$$\mathrm{O}\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}=\frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}\mathrm{E}\mathrm{D}}{\mathrm{N}\mathrm{I}\mathrm{R}+\mathrm{R}\mathrm{E}\mathrm{D}+0.16}\times 1.16$$

(23)

where NIR and RED are the surface reflectance values (decimals), for the near-infrared and red bands, provided by a given RS platform (dimensionless).

The f_c model used in the TSEB_ser and TSEB_par approaches is presented by Equations (24)–(28), as follows [21]:

$${\mathrm{f}}_{\mathrm{c},\mathrm{o}}=1-\mathrm{e}\mathrm{x}\mathrm{p}(-0.50\mathrm{L}\mathrm{A}\mathrm{I})$$

(24)

$$\mathrm{L}\mathrm{A}{\mathrm{I}}_{\mathrm{L}}=\mathrm{L}\mathrm{A}\mathrm{I}/{\mathrm{f}}_{\mathrm{c},\mathrm{o}}$$

(25)

$${\mathrm{f}}_{\mathrm{s}}=1+{\mathrm{f}}_{\mathrm{c},\mathrm{o}}\times \exp \left(-0.50\mathrm{L}\mathrm{A}{\mathrm{I}}_{\mathrm{L}}\right)-{\mathrm{f}}_{\mathrm{c},\mathrm{o}}$$

(26)

$$\mathrm{C}\mathrm{F}=-\mathrm{l}\mathrm{n}\left(\frac{{\mathrm{f}}_{\mathrm{s}}}{0.50\mathrm{L}\mathrm{A}\mathrm{I}}\right)$$

(27)

$${\mathrm{f}}_{\mathrm{c}}=1-\mathrm{e}\mathrm{x}\mathrm{p}(-0.50\times \mathrm{C}\mathrm{F}\times \mathrm{L}\mathrm{A}\mathrm{I})$$

(28)

where f_c,o is the initial f_c value before adjustments (dimensionless); LAI_L is the local LAI (m²/m²); f_s is the soil fractional cover (dimensionless); CF is the vegetation clumping factor (dimensionless).

The LAI is calculated using the model from [38], an exponential model calibrated for maize and sorghum and indicated by Equation (29) below:

$$\mathrm{L}\mathrm{A}\mathrm{I}=0.263\times \mathrm{e}\mathrm{x}\mathrm{p}(3.813\times \mathrm{O}\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I})$$

(29)

Maize h_c was estimated through an exponential model for maize and soybeans and indicated by Equation (30) below [39]:

$${\mathrm{h}}_{\mathrm{c}}=\left(1.86\mathrm{O}\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}-0.20\right)\times \left[1+4.82\times 10^{-7}\mathrm{e}\mathrm{x}\mathrm{p}\left(17.69\mathrm{O}\mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}\right)\right]$$

(30)

2.4. Remote Sensing Platforms

2.4.1. Spaceborne

Landsat-8

Landsat-8 is a satellite-based RS platform jointly managed by the United States Geological Service (USGS) and the National Aeronautics and Space Administration (NASA). A Landsat-8 satellite is equipped with an operational land imager (OLI) and a thermal infrared sensor (TIRS). These instruments capture images of Earth’s surface, with the OLI providing data at a spatial pixel resolution of 30 m and the TIRS at 100 m, resampled to 30 m [40]. These images are acquired every 16 days. The OLI sensor captures shortwave multispectral data, while the TIRS camera records longwave infrared (LWIR) thermal radiation images. The study sites, LIRF and IIC, are strategically located in the overlapping region of Landsat-8 scenes with path/row designations of 33/32 and 34/32, respectively. As a result, the temporal resolution for data acquisition for Landsat-8, in this study, was once every eight days depending on sky cloudiness conditions near noon time.

Landsat-8 follows a sun-synchronous orbit around Earth, orbiting at an altitude of 705 km, with equator crossings occurring at approximately 11:30 a.m. local time. The original radiometric resolution of Landsat-8 imagery is 12 bits, but this is enhanced to 16 bits after post-processing by USGS/NASA. To convert digital numbers (DN) into ground surface reflectance and nadir-looking temperature (Landsat-8 Level-2 imagery), linear calibration coefficients are provided in the metadata of the imagery file. Level-2 images undergo rigorous calibration procedures, eliminating the need for additional post-processing once the final surface reflectance and temperature images are accurately derived from the original DN values, following the methods outlined by Roy et al. [40] in 2014. Further details regarding atmospheric corrections applied to Landsat-8 imagery can be found in [41].

Table 3. Landsat-8 multispectral bands used in this study at both LIRF and IIC sites.

Bands	Central Wavelength (nm)	Bandwidth (nm)	Spatial Resolution (m)
BLUE	480	60	30
GREEN	560	60	30
RED	655	30	30
NIR	870	30	30
LWIR	1090	60	100 (Resampled to 30)

presents the spectral characteristics of the Landsat-8 bands considered in this research. The pixels that overlapped with the measurement stations at each research site were considered to be representative for the estimation of maize ET_a using the TSEB algorithm.

Sentinel-2

The Sentinel-2 satellites are under the care and operation of the European Space Agency (ESA), an intergovernmental organization representing 22 European countries. The management of the Sentinel satellite missions falls within the purview of the Copernicus Programme. This satellite constellation comprises two units, Sentinel-2A (S2A) and Sentinel-2B (S2B), taking Earth’s landscape images at noon. Each of these satellites takes multispectral image scenes, covering an area of 290 km by 290 km. They orbit the Earth, providing imagery every 10 days for a single satellite device around the equator. When both satellites are combined, this interval shortens to 5 days. For areas located at mid-latitudes, such as the LIRF and IIC research facilities, the revisiting time is as frequent as every 2 to 3 days near local noon time. Additional details concerning the satellite’s design, operation, and components can be found in the work of [42].

Both the S2A and S2B satellites follow sun-synchronous orbits, maintaining an average altitude of 786 km, with equator crossings occurring around noon local time (12 p.m.). It is important to note that Sentinel-2 satellites currently do not provide thermal imagery. The spatial resolution of Sentinel-2 images varies and depends on the specific multispectral bands in use, ranging from 10 m to 60 m. For this study, we have considered only Sentinel-2 bands 2, 3, 4, and 8, which offer a spatial resolution of 10 m, as these bands are provided at their original spatial resolution, as documented in

Table 4. Sentinel-2 multispectral bands used in this study at LIRF and IIC research sites.

Bands	Central Wavelength (nm)	Bandwidth (nm)	Spatial Resolution (m)
BLUE	492	66	10
GREEN	560	36	10
RED	665	31	10
NIR	833	106	10

.

The original radiometric resolution of Sentinel-2 images is 12 bits. However, ESA enhances these images to a 16-bit radiometric resolution through post-processing. Sentinel-2 Level-2 images undergo calibration and pre-processing, aiming to provide ground-based surface reflectance. The surface reflectance images provided (downloaded) include a scaling factor of 10,000. To obtain surface reflectance decimal values, the Sentinel downloaded images are divided by 10,000. The atmospheric corrections are carried out using a radiative transfer algorithm developed by ESA, known as Sen2Cor. For details regarding the use of Sen2Cor in atmospherically correcting S2A and S2B satellite images, refer to the work of [43]. The pixels that contained the stations of measurements for each research site were considered representative for the calculations of maize ET_a using the TSEB approaches.

Planet CubeSat

Planet CubeSat is a cost-effective commercial constellation of micro (Dove) satellites, managed by Planet Labs, Inc. in San Francisco, California, USA. Comprising over 130 CubeSat units, these satellites observe Earth’s landscapes, providing high temporal (daily) and spatial resolutions (3 m). Planet CubeSat microsatellites capture multispectral imagery in the visible and near-infrared (NIR) portions of the light spectrum. The radiometric resolution of Planet CubeSat imagery starts at 12 bits during image acquisition but is enhanced to a 16-bit resolution through post-processing before it is made accessible for download. Planet CubeSat satellites are notably smaller and lighter (0.10 m × 0.10 m × 0.30 m and 4 kg, respectively) compared to their counterparts, like Landsat-8 and Sentinel-2. Operating in a sun-synchronous orbit with an altitude ranging from 450 km to 580 km, these satellites pass over the equator between 9:30 and 11:30 a.m. local time. We used imagery files from 11:30 a.m. since it was the closest to the other time of RS data acquisition regarding the other RS sensors in this study. Planet CubeSat’s capabilities are exclusively limited to multispectral imagery, as detailed in

Table 5. Planet CubeSat multispectral bands used in this study at LIRF and IIC research sites.

Bands	Central Wavelength (nm)	Bandwidth (nm)	Spatial Resolution (m)
BLUE	491	60	3
GREEN	566	90	3
RED	666	80	3
NIR	867	80	3

. Unfortunately, thermal imagery is still not provided by these constellations of microsatellites.

The surface reflectance images obtained by Planet CubeSat undergo pre-processing and calibration, also incorporating a scaling factor of 10,000. The imagery pre-processing and calibration process encompass adjustments for radiation scattering due to atmospheric gases, aerosol concentration, and their altitude-dependent variations between Earth’s surface and the satellite’s at-sensor camera in space. To enhance their calibration, data from MODIS (Moderate Resolution Imaging Spectroradiometer), such as water vapor, ozone, and aerosol quality control products, are utilized, along with the 6SV2.1 radiative transfer model. However, the atmospheric correction process is still evolving, primarily because Planet’s imagery calibration does not currently address effects like stray light, haze, and the influence of thin cirrus clouds. The approach by Planet Labs assumes that Earth’s landscapes behave as Lambertian surfaces, scattering light uniformly in all directions, and that all scenes are effectively at sea level. Following image acquisition, geometric corrections are meticulously executed using sensor telemetry, ground control points (GCP), and finely detailed digital elevation models (DEM). Furthermore, the Planet Team released a harmonized version of Planet imagery, including CubeSat data, aligning the quality of multispectral data with (calibration to) Sentinel-2 standards. For further details on the image harmonization processes, refer to the works of [44,45]. In this study, the primary data source used was the CubeSat harmonized images.

In our study using the data from the LIRF and IIC using all listed spaceborne RS sensors, we specifically focused on clear-sky images. This deliberate choice was made to ensure that cloud cover conditions did not interfere with the accuracy of ground-based surface reflectance and temperature values in our research fields. For Landsat-8, the evaluation of evapotranspiration (ET_a) included two different datasets. The first dataset involved the use of the original Landsat-8 platform, encompassing imagery data from OLIS (multispectral) and TIRS (thermal). The second Landsat-8 dataset combined OLIS data with ground-based infrared temperature (IRT) data. The utilization of TIRS and ground-based T_s data was specific to the Landsat-8 spaceborne platform since Planet CubeSat and Sentinel-2 relied on ground-based nadir-looking IRT T_s data for input in the TSEB RS of the ET_a algorithms. This choice, in the case of the Landsat-8 imagery, was made due to the pivotal role of its TIRS data, albeit its original pixel size, in providing original spaceborne T_s information for maize ET_a estimation. The Planet Cubesat pixel data that had the measurement stations at LIRF and IIC were used for the estimation of maize ET_a using the TSEB algorithms.

2.4.2. Proximal

At the LIRF and IIC research sites, proximal surface reflectance and nadir-looking T_s data was obtained using a handheld multispectral radiometer (MSR5, CropScan Inc., Rochester, MN, USA). The MSR5 radiometer is a compact device consisting of a quasi-cubic radiometer measuring 0.80 × 0.80 × 0.10 m. It features an integrated IRT sensor from Exergen Corporation in Watertown, Massachusetts. The radiometer has a field of view (FOV) of 28 degrees and captures readings at an altitude of 2.2 m above ground level (AGS). This measurement setup covers an area on the ground equivalent to a 1-meter-diameter circle, with a 2V:1H aspect ratio (

Table 6. The MSR5 multispectral bands used in this study at the LIRF and IIC research sites.

Bands	Central Wavelength (nm)	Bandwidth (nm)	Spatial Resolution (m)
BLUE	485	70	1
GREEN	560	80	1
RED	660	60	1
NIR	830	140	1

). The attached Exergen IRT, with a FOV ratio of 3V:1H, results in a spatial footprint of 0.80 m in diameter. The MSR5 is a passive sensor, relying on natural sunlight for data collection of surface temperature. It replicates the spectral characteristics of Landsat-5, obtaining data in the visible, NIR, and MIR (mid-infrared) light spectrum. For this study, we deployed MSR5 units at both the LIRF and IIC sites. The multispectral data from these MSR5 devices were sampled about once a week, with four readings at each measurement location across all research sites: more specifically, two readings within the crop rows and two readings in the maize inter-row spaces per site. These readings around solar noon (11:30 a.m. to 1:30 p.m.) were then averaged to obtain the final measurements per sampling location in each research site per field visit. The MSR5 measurements’ locations were the reference for the selection of pixels from the different images/platforms or scales, considering areas with homogenous soil texture and canopy cover.

Airborne

The USDA-ARS Water Management and Systems Research Unit and the CSU Drone Center scheduled UAS missions for the research sites around local solar noon (11:30 a.m. to 1:30 p.m. MST). At the LIRF, the USDA-ARS team was responsible for conducting UAS missions, while the CSU Drone Center took charge of the UAS missions at the IIC site. These UAS images were acquired using a MicaSense RedEdge-MX multispectral camera (MicaSense Inc., Seattle, WA, USA), seamlessly integrated into the airborne platform. The RedEdge-MX detector captures data across the visible and invisible light spectrum, with bands including BLUE (475 nm, 32 nm bandwidth), GREEN (560 nm, 27 nm bandwidth), RED (668 nm, 14 nm bandwidth), and NIR (842 nm, 57 nm bandwidth).

The UAS’s surface reflectance imagery data serve as another RS sensor for this study, and they were complemented by nadir-looking T_s data obtained from point-based measurements conducted at each of the research sites. This combined dataset is utilized as input for estimating hourly maize ET_a using the TSEB approaches. For an overview of the UAS missions carried out in 2020 and 2021 at LIRF and IIC, refer to

Table 7. The UAS mission summary for the USDA-ARS and CSU Drone Center at all sites.

	USD—ARS	CSU Drone Center
UAS Unit	DJI S900	DJI M600
Flight Altitude (m)	120	100
UAS Speed (m/s)	5	5
Temporal Resolution	Weekly	Weekly
Imagery Pixel Size (m)	0.03	0.08
Overlap/Sidelap Percentage (%)	88/70	80/70
Calibrated Reflectance Panel	Yes	Yes
Orthorectified Coordinate System	WGS84 UTM	WGS84 UTM
Post-processing Imagery Software	Agisoft Metashape	Pix4D v4.5.6

.

2.5. Field Data Collection

The experiment was replicated in both the LIRF and IIC sites to obtain similar datasets for the evaluation of the airborne, spaceborne, and proximal platforms’ derived RS data when used in the prediction of ET_a using TSEB RS algorithms. A total of three field measurement stations provided the ground-based input data to estimate and evaluate maize ET_a. The following data were measured at each station: R_n, G, nadir-looking surface radiometric temperature (T_s), incoming shortwave solar radiation (R_s), shallow soil temperature, and volumetric water content. A flux tower provided measurements of H and LE in each research site (Figure 3 and Figure 4a).

2.5.1. Surface Heat Fluxes

At the IIC site, a two-way NR-Lite and two four-way CNR1 net radiometers (Kipp and Zonen, Delft, The Netherlands) measured R_n at the height of 3.3 m AGS. At the LIRF site, each of the three net radiometers was a two-way NR-Lite in 2020 and 2021. The CNR1 radiometer was installed on Field F (west station). NR-Lite radiometers measure net shortwave and longwave radiation within a spectral range from 0.2 to 100 μm, temperature dependency of 0.12%/°C, and a directional error of less than 30 W/m² at, at least, 1000 W/m² [46]. The CNR1 radiometer provides data regarding all four terms of the net radiation budget, and it has a measurement uncertainty that is within 10 to 35 μV/W/m² and a directional error of 25 W/m² at 1000 W/m² [47].

Surface G data were determined using the soil heat flux plate method. At LIRF, two HFT3-L soil heat flux plates (Radiation and Energy Balance Inc., Bellevue, WA, USA) were buried at 0.08 m between and below maize rows at each measurement station. At the IIC site, the HFT3-L plates were placed between two consecutive maize rows due to the flooded furrow during irrigation events. One 5TE soil water content sensor (Decagon Devices Inc., Pullman, WA, USA) was buried at 0.04 m. Two T107 temperature probes (Campbell Scientific Inc., Logan, UT, USA) were installed at 0.02 and 0.06 m below ground surface (BGS) to determine average soil temperature and for the calculation of soil heat storage above the 0.08 m soil layer from the HFT3-L plates. The HFT3-L sensors have thicknesses and diameters equal to 3.91 and 38.2 mm, respectively. The measurement uncertainty of soil heat flux from the plates is 5% [48].

Measured LE and H data were produced from acquired high-frequency wind speed, air temperature, and water vapor measurements using an eddy covariance (EC) system installed at each research site. At the IIC (Figure 5a), the EC system consisted of an LI-7500A open-path CO₂/H₂O gas analyzer (LI-COR Biosciences, Lincoln, NE, USA) and a CSAT three-dimensional (3D) sonic anemometer (Campbell Scientific Inc., Logan, UT, USA). At LIRF (Figure 5b), an LI-7500DS open-path gas analyzer (LI-COR Biosciences, Lincoln, NE, USA) and a Gill WindMaster three-dimensional (3D) sonic anemometer (Gill Instruments, Lymington, Hampshire, UK) provided measurements of LE and H, respectively. Both EC systems at LIRF and IIC were installed at 3.5 m AGS, positioned facing the prevailing wind direction at each site (135° azimuth angle), and set to a sampling frequency equal to 10 Hz. The EC turbulent fluxes and respective ancillary data were recorded as 15-min and half-hour averages at the IIC and LIRF sites, respectively. The EC system often provides imbalanced turbulent fluxes regarding SEB closure [49,50,51], with the closure SEB ratio “(H + LE)/(R_n − G)” ranging from 70 to 90% [52,53]. To improve the representativeness of H, LE, and ET_a measurements from the EC system, the residual-LE closure approach was chosen in this study to ensure the closure of the surface heat fluxes. The work of [54] indicated that the residual-LE method calculates measured LE as the difference among measured R_n, G, and H (from the EC system) and that most of the unresolved EC system closure issues are due to LE rather than H.

Table 8. Correction methods applied to the EC data at LIRF and IIC.

Correction Method	Source	Research Site
Wind coordinate or tilt correction	[55,56]	LIRF and IIC
Air density fluctuation—the Webb–Pearman–Leuning (WPL) correction	[57]	LIRF and IIC
Humidity correction of sonic temperature	[58,59]	LIRF and IIC
Statistical analysis of data screening	[60]	LIRF
The angle of attack correction for 3D wind components	[61]	LIRF

shows the corrections performed in the high-frequency EC data at LIRF and IIC in 2020 and 2021.

A two-dimensional (2D) EC heat flux (source) footprint analysis was performed to filter the EC-derived heat flux data to consider only flux source areas contributing to H and LE fluxes coming strictly from the maize fields at both the LIRF and IIC sites (as described in [62]), which is an analytical heat-flux-source approach that provides 2D footprint extents based on turbulence characteristics of the air flow and surface, such as Monin–Obukhov atmospheric stability length (L_MO, m), friction velocity (u_*, m/s), the standard deviation of lateral velocity (σ_v), Z_u, Z_om, and the atmospheric boundary layer height (H_L). To compare the predictions of ET_a at each station of measurement with the hourly and daily EC ET_a data, it was assumed that the fixed measurement instrumentation stations that were within the 2D EC footprint were representative of observed ET_a data from the EC system during the RS platform overpass date and time (Figure 5). At the IIC site, the EC flux tower was located at the northwest corner of the field in 2020. Since the west and east measurement stations (Field F on Figure 4a) were farther from the footprint area for H and LE fluxes, the data from the west station was assumed to represent a maize ET_a comparison between the EC data and the remote sensing of ET_a predictions since it was the closest station to the flux tower. Figure 6 and Figure 7 show the 2D EC footprints that served as a reference to filter the EC data at the LIRF and IIC sites, respectively.

2.5.2. Micrometeorological Data

Micrometeorological data were obtained at the EC heat flux or weather tower in each field. An HMP45C sensor probe (Vaisala, Helsinki, Finland) measured air temperature (T_a, °C) and relative humidity (RH, %) at 3.5 m AGS. The HMP45C probe has a 1000-ohm thermometer with measurement uncertainty within 0.20 to 0.30 °C at ambient temperatures varying from 20 to 40 °C. Within the HMP45C sensor, a HUMICAP H-chip measures RH with an uncertainty of approximately ±1% when T_a equals 20 °C [63]. Wind speed and direction were measured using a three-dimensional (3D) sonic anemometer installed at each research site, 3.5 m AGS. At IIC, the 3D sonic anemometer was a CSAT (Campbell Scientific Inc., Logan, UT, USA). In contrast, at LIRF, the 3D sonic anemometer was a Gill WindMaster 3D sonic anemometer (Gill Instruments, Lymington, Hampshire, UK). Incoming shortwave solar radiation was measured using an LI-200X pyranometer (LI-COR, Lincoln, NE, USA) at LIRF (3.3 m AGS) and the CNR01 net radiometer at the IIC site. The on-site weather data were recorded every minute and averaged every 15 min in CR1000 and CR3000 dataloggers (Campbell Scientific Inc., Logan, UT, USA) at the IIC and LIRF sites, respectively.

2.5.3. Surface Temperature Data

Nadir-looking T_s data were measured using SI-111 IRT sensors (Apogee Instruments, Logan, UT, USA). Each measurement station had one SI-111 IRT sensor installed 1 m above the canopy, nadir-looking, during the data collection campaign at IIC and LIRF. The IRT sensors were placed in a 4-meter-tall vertical post and raised to higher heights until the canopy reached maximum h_c, always keeping the 1-meter distance between the sensor and the top of the canopy. The FOV of SI-111 IRT sensors is a 22° half-angle. The sensors have a fast response time (<1 s) and a small uncertainty (±0.20 °C) when the target temperatures are between −20 °C and 65 °C. The ground-based area where most of the thermal radiation is sensed by the IRTs was equivalent to a 3-meter-diameter circumference. If we consider the upper part of the canopy, the footprint of the area sensed by the IRTs was about 1 m in diameter. The T_s data, alongside measured R_n and buried sensors to calculate surface G, were recorded every minute and averaged every 15 min in either a CR1000 or CR3000 datalogger (Campbell Scientific Inc., Logan, UT, USA) in each measurement station at LIRF and IIC. The SI-111 IRT sensors provided the input data to run the TSEB approaches in this study considering all RS sensors but the MSR.

2.6. Statistical Data Analysis

The following statistical variables have been considered to compare the performance of the different ET_a models across the spaceborne and airborne RS platforms: mean bias error (MBE), root mean square error (RMSE), normalized MBE (NMBE), normalized RMSE (NRMSE), and the coefficient of determination (R²). Equations (31)–(34) indicate MBE, NMBE, RMSE, and NRMSE, respectively:

$$MBE=\left(\frac{1}{n}\right)\sum _{i=1}^n\left(E_i-O_i\right)$$

(31)

$$NMBE=\left(\frac{MBE}{\overline{O}}\right)\times 100\%$$

(32)

$$RMSE=\sqrt{\left(\frac{1}{n}\right)\sum _{i=1}^n{\left(E_i-O_i\right)}^2}$$

(33)

$$NRMSE=\left(\frac{RMSE}{\overline{O}}\right)\times 100\%$$

(34)

where $\overline{O}$ is the mean of the observed data; n is the sample size; E_i and O_i are the estimated and observed values, respectively. NMBE and NRMSE are given in percentages, while Equations (32) and (34) provide statistical indicators with the same units of the primary variables. Based on the guidelines in [64], the performances of the ET_a models have been classified into one of the following categories: excellent (NRMSE ≤ 10%), good (10% < NRMSE ≤ 20%), fair (20% < NRMSE ≤ 30%), and poor (NRMSE > 30%).

The R², in the context of model performance assessment, informs about the degree of variability in the observed data explained by the modeling approach. Equation (35) gives the mathematical expression for R²:

$${\mathrm{R}}^2=\frac{\sum ({\mathrm{E}}_{\mathrm{i}}-\overline{\mathrm{E}})({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}})}{\sqrt{\left[\sum {\left({\mathrm{E}}_{\mathrm{i}}-\overline{\mathrm{E}}\right)}^2\right]\left[\sum {\left({\mathrm{O}}_{\mathrm{i}}-\overline{\mathrm{O}}\right)}^2\right]}}$$

(35)

where $\overline{E}$ is the mean value of the predictions. This study defines the optimal remote sensing platform as the source of multispectral data with the smallest NRMSE. In case two or more platforms have identical NRMSE, the highest d_r index between the two platforms is considered the optimal data for a given remote sensing of the ET_a algorithm. The modified index of agreement (d_r) was calculated to assess model performance by comparing the sum of the residuals to the total difference between observed values and the respective mean of the measured data. Higher d_r values indicate that the predicted values have more statistical agreement with the observed data, showing better model performance [65].

Outliers have been excluded from the analysis based on the median absolute deviation approach (MADA). The MADA method for filtering extreme values in a dataset uses the median instead of the mean as a central tendency measure. The median allows for flagging points that do not conform with the sampled data’s trends [66]. The MADA index is defined by Equation (36) when a Gaussian distribution assumption is considered for the data without the influence of extreme values [67].

$$\mathrm{M}\mathrm{A}\mathrm{D}\mathrm{A}=1.4826\times \mathrm{M}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left[|x_i-Median\left(x\right)|\right]$$

(36)

where x_i is the value of a given variable at a specified timestep; Median(x) is the median of the variable’s sample size. In this study, the criteria for filtering the data for potential outliers was the recommendation by [66]. The median ± 2.5 times the MADA index is the cutoff value expected in each sampled dataset.

3. Results

3.1. Error Analysis from the TSEB Algorithms

3.1.1. TSEB Data Analysis Combined (LIRF and IIC 2020–2021)

The results from TSEB_par and TSEB_ser of the ET_a algorithms indicated that data used from the MSR5 (1 m) RS sensor resulted in the best performance when estimating hourly maize ET_a compared to the use of other RS platforms’ data, since the smallest NRMSE was 11% and 14%, respectively, for both TSEB algorithms used. The overall errors in predicting hourly maize ET_a using MSR5 (1 m) data were −0.02 (−3%) ± 0.07 (11%) mm/h and −0.02 (−4%) ± 0.09 (14%) mm/h for the TSEB_par and TSEB_ser, respectively. The errors associated with Landsat-8 (30 m) TIRS data used in the RS of ET_a estimation were −0.03 (−5%) ± 0.16 (29%) mm/h and −0.07 (−13%) ± 0.15 (29%) mm/h, while reductions of 31% and 34% in NRMSE were observed when the TIRS was replaced by the ground-based IRT data for the TSEB_par and TSEB_ser, respectively (Figure 8 and Figure 9).

The NRMSE errors in estimating H using Landsat-8 multispectral visible and NIR and TIRS data were greater than 30% for the TSEB_par and TSEB_ser, while the NRMSE from the remaining spaceborne, airborne, and proximal platforms were within 20 to 30% (Figure 10 and Figure 11). Using proximal IRT data to estimate maize ET_a with Landsat-8 (30 m), surface reflectance improved the estimation of hourly maize ET_a, with reduced model performance errors of 0.02 (4%) ± 0.10 (19%) mm/h and −0.04 (−8%) ± 0.11 (20%) mm/h for the TSEB_par and TSEB_ser, respectively. Most of the improvements in hourly maize ET_a using TSEB_par and TSEB_ser were due to better H predictions than the accuracy on the estimation of R_n and G fluxes. For the TSEB_par, the performance of the R_n model did not significantly change when using TIRS and IRT data since the NRMSE was 13% and 14%, respectively. The TSEB G model NRMSE ranged from 20 to 45% across all the RS platforms in this study, which had a fair to poor performance since NRMSE > 20%. However, since the magnitude of G is much smaller than R_n and H, the error propagation in LE estimates due to G is also smaller than the other two SEB input fluxes. Furthermore, the original TSEB G model assumes that G varies linearly with the R_n flux associated with bare soil. However, the work of [18] indicated that G and bare soil R_n have a time-phased difference between the two fluxes that requires a specific local calibration.

A consistent trend of maize ET_a underestimation was present for all spaceborne platforms/sensors (data) when considering the TSEB_ser RS of the ET_a algorithm. The TSEB_ser hourly maize ET_a underestimation range was within −3% (CubeSat) to −13% (Landsat-8 TIRS). When considering the TSEB_par algorithm, there was an underestimation of hourly maize ET_a when using UAS (−3%), MSR5 (−3%), Sentinel-2 (−3%), and Landsat-8 TIRS (−5%). The slight overestimations of hourly maize ET_a for CubeSat (3 m) and Landsat-8 IRT (30 m) were 2% and 4%, respectively. Sentinel-2 (10 m) and Planet CubeSat (3 m) had similar performances since the NRMSE for TSEB_par and TSEB_ser had a 6% difference between the two high-resolution spaceborne platforms.

The underestimation trend shows a concurrent bias in the TSEB_ser algorithm’s ability to predict actual ET_a rates across various platforms. The largest underestimation by Landsat-8 TIRS could be attributed to its thermal infrared sensor’s spatial resolution and spectral characteristics, which may not fully capture the heterogeneous surface temperature and moisture conditions within agricultural fields. This result suggests that there are limitations, inherent to this sensor, in accurately detecting surface temperatures for relatively small maize fields displaying some surface heterogeneity. The slight overestimations noted with CubeSat and Landsat-8 IRT suggest that these sensors, at their respective resolutions, might capture more of the variance in surface temperatures and moisture levels than the TIRS sensor. These discrepancies underscore the importance of the spectral and spatial resolution’s role in capturing field heterogeneity.

The comparative performance of Sentinel-2 and Planet CubeSat, with a minimal difference in NRMSE for both TSEB algorithms, indicates that high-resolution spaceborne platforms are capable of providing reliable ET_a estimates, as compared to Landsat-8. The similar performance of these two sensors, despite their different operational designs and data characteristics, suggests that high-spatial-resolution spaceborne sensors have better results for small agricultural fields, compared to coarse satellite pixel resolutions. It also highlights the critical role of spatial resolution in capturing the detailed variability of agricultural landscapes, which is essential for promoting conditions to enhance sustainable irrigation water management in cropland fields.

In addition, the trend of underestimation observed with the TSEB_par algorithm across UAS, MSR5, Sentinel-2, and Landsat-8 TIRS (with relatively minor variations) suggests the need for algorithmic adjustments or enhancements to further improve ET_a estimates. The TSEB RS of the ET_a models were originally developed and validated using airborne data at a 5 m spatial resolution. Given the large number of input parameters in the TSEB RS of ET_a, some of the uncertainty in the final maize ET_a estimation might be due to sub-models that do not have locally calibrated parameters for the semi-arid region across different climate conditions. Nevertheless, the performance of the ET_a algorithm is acceptable.

In our study, the range of accuracy of both TSEB approaches for maize ET_a (10–30%) was consistent with previously reported values in the current existing literature for other row crops such as vineyards when using sUAS imagery aggregated to different local spatial scales (e.g., 4 m, 7 m, 14 m, and 30 m) [27] and soybeans and maize when using Landsat-7 (30 m) and UAS (6 m) imagery spatial scales [68]. The similar results obtained between maize ET_a observed and estimated values from previous studies using different spatial and spectral resolution RS sensors reinforces the reliability of selecting finer RS spatial resolution sensors to use as input for predicting crop ET_a with a TSEB RS of the ET_a algorithm to advance sustainable agricultural water management. Both Ref. [27] and Ref. [68] indicated that the smallest errors observed when predicting LE and crop ET_a occurred at spatial resolutions less than 15 m.

The comparative analysis of multiple RS sensors presented a relevant understanding of their data quality as input for predicting maize ET_a using the TSEB RS of the ET_a algorithms evaluated in the study. Despite the overall agreement in ET_a accuracy observed in the existing literature, it is important to consider factors such as spectral resolution, revisit frequency, and sensor-specific limitations when deciding which RS sensor (data) is better, given the site-specific conditions and sustainable irrigation needs. Since the MSR5 (proximal RS sensor) has limitations regarding the spatial coverage of large fields in a timely manner, operational costs, and data processing, its use could be hindered throughout different stages of the growing season. Given that the high-spatial-resolution spaceborne platforms (e.g., Sentinel-2 and Planet CubeSat) had similar accuracy performances (NMBE and NRMSE) compared to the MSR5 (Figure 10 and Figure 11), their use could be justified when proximal RS devices are not ideal for collecting data over large agriculture fields with significant areas of canopy heterogeneity that make on-site data acquisition challenging and unrealistic.

3.1.2. TSEB Data Analysis Separately for LIRF and IIC 2020–2021 Data

When evaluating the LIRF 2020–2021 data alone, the results were also consistent with the previous combined data from LIRF and IIC regarding both the TSEB_ser and TSEB_par algorithms. The proximal platform MSR5 (1 m) outperformed all other platforms since it had the smallest NRMSE equal to 10% (TSEB_par) and 14% (TSEB_ser) see

Table 9. Error analysis from the TSEB_par hourly maize ET_a evaluation for the LIRF and IIC 2020–2021 data.

Site	Platform	Spatial Resolution (m)	n	MBE (mm/d)	NMBE (%)	RMSE (mm/d)	NRMSE (%)	R2
LIRF	Landsat-8 TIRS	30	14	−0.04	−7%	0.13	25%	0.51
	Landsat-8 IRT	30	14	0.02	3%	0.09	18%	0.73
	Sentinel-2	10	35	−0.02	−4%	0.09	16%	0.85
	CubeSat	3	63	0.02	4%	0.10	18%	0.69
	MSR5	1	17	−0.01	−1%	0.06	10%	0.91
	UAS	0.03	11	−0.03	−5%	0.09	15%	0.67
IIC	Landsat-8 TIRS	30	12	−0.02	−3%	0.18	33%	0.24
	Landsat-8 IRT	30	12	0.03	5%	0.11	21%	0.76
	Sentinel-2	10	32	−0.01	−2%	0.09	15%	0.73
	CubeSat	3	62	0	0%	0.08	14%	0.76
	MSR5	1	28	−0.03	−4%	0.08	12%	0.82
	UAS	0.03	13	−0.01	−2%	0.09	14%	0.58

and

Table 10. Error analysis from the TSEB_ser hourly maize ET_a evaluation for the LIRF and IIC 2020–2021 data.

Site	Platform	Spatial Resolution (m)	n	MBE (mm/d)	NMBE (%)	RMSE (mm/d)	NRMSE (%)	R2
LIRF	Landsat-8 TIRS	30	14	−0.07	−13%	0.14	27%	0.48
	Landsat-8 IRT	30	14	−0.07	−13%	0.10	20%	0.81
	Sentinel-2	10	35	−0.01	−2%	0.08	15%	0.83
	CubeSat	3	63	−0.03	−6%	0.09	16%	0.79
	MSR5	1	17	0	0%	0.08	14%	0.82
	UAS	0.03	11	−0.04	−8%	0.08	14%	0.78
IIC	Landsat-8 TIRS	30	12	−0.07	−13%	0.16	30%	0.41
	Landsat-8 IRT	30	12	−0.01	−2%	0.12	21%	0.71
	Sentinel-2	10	32	−0.05	−8%	0.11	17%	0.75
	CubeSat	3	62	0.01	1%	0.08	13%	0.78
	MSR5	1	28	−0.04	−6%	0.09	14%	0.75
	UAS	0.03	13	−0.09	−13%	0.11	17%	0.72

. Underestimation (NMBE) of TSEB_par hourly maize ET_a predictions ranged from −1 to −7% for the Landsat-8 TIRS (30 m), Sentinel-2 (10 m), MSR5 (1 m), and UAS (0.03 m). Overestimation of TSEB_par maize ET_a was observed in the case of RS data used from the Planet CubeSat (3 m) and Landsat-8 IRT platforms, with respective NMBE equal to 4% and 3%. Regarding the TSEB_ser RS of the ET_a algorithm, similar trends of underestimation and overestimation of hourly maize ET_a were observed compared to the TSEB_par algorithm. There were 28% and 26% reductions in NRMSE, respectively, for the TSEB_par and TSEB_ser RS of ET_a results when Landsat-8 surface reflectance and IRT T_s data were used to estimate hourly maize ET_a. When evaluating the IIC 2020–2021 data alone, the results were again consistent with the previous analysis regarding LIRF and the combined LIRF-IIC data. Similarly, the MSR5 (1 m) outperformed the spaceborne and airborne platforms/sensors, with an overall error equal to −0.03 (−4%) ± 0.08 (12%) mm/h and −0.04 (−6%) ± 0.09 (14%) for the TSEB_par and TSEB_ser, respectively.

4. Discussion

The observed results were related to differences in RS sensor types, the assumptions of the TSEB ET_a model and inherent uncertainty, and the complex physical processes that derive the heat and water vapor transfer between the surface and atmosphere. Regarding the RS data characteristics, the spatial resolution significantly impacts the accuracy of the hourly maize ET_a.

The RS platforms with higher spatial resolution (<10 m) capture finer surface feature details within agricultural fields. These high-spatial-resolution data have the potential to better characterize spatial variability in soil and vegetation conditions, which is essential for accurate crop ET_a estimation using the TSEB RS of ET_a algorithm. The RS sensors with coarse spatial resolution, such as Landsat-8 (30 m), have limitations in providing relevant multispectral data that represents well local variations of the Ts and surface reflectance values for smaller agricultural fields. These limitations can lead to increased uncertainty in ETa estimates, particularly in row crop fields.

The integration of ground-based measurements, such as IRT nadir-looking Ts data, has been shown to significantly improve the accuracy of RS-based ET_a estimates, especially when coarse Ts spatial-resolution data are not representative of local field conditions. The observed underestimation of the ET_a estimation, when using the proximal, airborne, and spaceborne RS surface reflectance data, may be related to uncertainties in the TSEB model parameters and the simplifications of the surface energy balance equations. Another challenge regarding the use of different RS sensors is the temporal resolution. Limited revisit frequencies and local atmospheric effects (e.g., aerosols) can introduce uncertainty in data acquisition and quality, ultimately impacting ET_a estimates.

Accurate, spatio-temporal ET_a predictions are essential for increasing crop yields while mitigating water scarcity issues, a critical factor in securing water sustainability within a diverse range of water stakeholders. While remote sensing data, particularly through the TSEB model, has significantly advanced our understanding of crop ET_a, challenges persist in implementing sustainable solutions within agricultural settings.

Despite the potential of RS to advance sustainable water management in irrigated fields, limitations exist. Spaceborne data such as those from Landsat-8, Sentinel-2, and Planet CubeSat often lack consistent daily imagery acquisition for ideal ET_a modeling conditions, and this inconsistency can create challenges in implementing daily irrigation scheduling based on RS data inputs. To overcome the limitations, we propose an approach that integrates multiscale RS imagery from different platforms to generate RS data on a daily basis with similar radiometric quality to be readily available for water stakeholders (e.g., water management authorities, engineers, and agronomists). This integrated multiscale approach can provide consistent daily RS data for calculating crop ET_a and supporting water management decisions at local and large scales in agricultural settings.

Beyond the RS platforms’ spectral resolution, revisit frequency, and sensor-specific limitations and calibration, other factors merit consideration for a comprehensive understanding and effective adoption of RS technologies in agriculture. Atmospheric variables such as humidity, temperature, concentration of gases, dust, and wind speed can significantly affect the accuracy of RS data, particularly for spaceborne sensors. Nonetheless, a multi-source RS data analysis can compensate for the limitations of individual RS sensors/platforms, offering an expanded perspective of crop water requirements that can lead to more sustainable irrigation practices in cropland.

5. Conclusions

This research was conducted in a semi-arid climate area, in maize fields irrigated with sub-surface drip and furrow irrigation systems, at two research sites in northern Colorado, USA. We aim to investigate the performance of two TSEB remote sensing of ET_a algorithms when using input data from different (multiscale) remote sensing sensors/platforms. The hypothesis was that the accuracy of RS of ET_a estimation depended on both the pixel spatial and spectral/radiometric resolutions of the multispectral data used and on the specific parameters within the RS of the ET_a algorithms. The primary conclusion is that, for both TSEB approaches evaluated (TSEB_par and TSEB_ser), the best remote-sensing-based surface reflectance and temperature data for predicting maize hourly ET_a were those from the handheld MSR5 radiometer. The second-best RS data were multispectral surface reflectance images from the UAS, Planet, and Sentinel-2 RS platforms (plus surface temperature from stationary IRT sensors). However, using RS data from Landsat (optical and TIRS) resulted in larger ET_a estimation errors.

While it is possible to estimate crop ET_a using various remote sensing platforms, selecting the most suitable RS data for a given ET_a algorithm has the potential to significantly enhance irrigation water management by using more accurate ET_a estimates. In this study, it was found that the accuracy of the ET_a predictions was not the same across the different remote sensing sensors.

The use of the appropriate remote sensing data (i.e., MSR5) with the TSEB remote sensing of ET_a algorithms, to optimize maize irrigation scheduling, presents a significant contribution toward advancing sustainability in irrigated agriculture. The combination of MSR5 multispectral and thermal data to determine the contributions of soil and vegetation components to ET_a can offer a more accurate understanding of water consumption in cropland ecosystems, compared to the most common Landsat data use.

To improve the effectiveness of sustainable solutions using remote sensing data, future research in sustainability should focus on refining the TSEB algorithms, integrating diverse datasets within the same data analysis context, and addressing challenges associated with scaling from local (e.g., farms) to regional (e.g., irrigation districts and watersheds) levels.

This study highlights the need for further research aimed at improving the data quality of sub-optimal remote sensing platforms/sensors when only those data are available. It is critical to develop imagery calibration protocols to improve the quality of the remote sensing data needed for the prediction of crop ET_a under different surface and climate conditions. This would help enable the use of the most desirable remote sensing data with high accuracy for effective irrigation water management. Also, we recognize the need for more research including a wider range of commercial crops to strengthen the analysis of how the TSEB approaches perform when estimating crop ET_a for other crop types.

Additionally, the role of advanced spatial data analysis and machine learning algorithms in processing and interpreting RS data could be an alternative to explore to improve the quality of the RS data for the sensors that did not perform better than the MSR5. These technologies can provide a possible framework to address complex patterns and relationships within large imagery datasets, facilitating more applied and predictive approaches to crop water use and stress levels. By leveraging the computational capabilities of these artificial intelligence models, researchers and practitioners can refine the application of TSEB RS algorithms and determine irrigation scheduling practices to meet the water requirements of specific crops under local field conditions.