A Dynamic Zero-Plane Displacement Height Approach to Improve Remote Sensing-Based Maize Actual Evapotranspiration

1. Introduction

Agriculture is the largest consumer of freshwater worldwide, accounting for approximately 70% of global withdrawals [1]. In the western U.S. Maize Belt, irrigation supports roughly 58% of total maize production, placing substantial pressure on limited surface and groundwater resources [2]. Improving irrigation efficiency under these conditions requires accurate, spatially explicit estimates of crop water use or crop actual evapotranspiration (ET_a).

Measuring or estimating ET_a is fundamental for appropriate irrigation water management. Remote sensing (RS)–based surface energy balance (EB) models estimate ET_a by determining net radiation (R_n), sensible heat flux (H), soil heat flux (G), and solving the simplified EB for latent heat flux (LE) [3,4]. These EB sub-models are well suited to arid and semi-arid regions, but comparisons with ground-based eddy covariance (EC) heat fluxes data show that RS-based estimates of H and LE exhibit root mean square errors (RMSE) ranging from 28–50+ W m⁻² [5,6,7]. These expected uncertainties emphasize the importance of accurate surface aerodynamic parameters estimation, including the zero-plane displacement height (d_o) for turbulent heat exchange [8].

Recent advances in RS offer new opportunities to better represent agricultural land surface heterogeneity. Multispectral satellite imagery combined with ground-based meteorological data has been widely used to estimate EB components [9,10], while unmanned aerial systems (UAS) provide high-resolution canopy height information that captures within-field structural variability more effectively than point measurements [11,12]. Although these developments improve the characterization of canopy structure and surface conditions, aerodynamic parameters such as d_o cannot be directly observed and must still be inferred from physical representations of wind–canopy interactions.

Friction velocity (u_*, m s^-1) and aerodynamic resistance to heat transfer (r_ah, s m^-1) strongly influence the estimation of H. Both u_* and r_ah depend on d_o, which is the height (m) above the surface where the horizontal wind speed (u, m s^-1) approaches zero due to surface roughness. Estimating d_o under field conditions is challenging because wind flow interacts dynamically with vegetation structure, producing strong spatial variability, particularly in row-crop systems [13,14,15,16,17]. Most existing d_o formulations assume homogeneous canopies and relate d_o to crop height (H_c, m) or leaf area index (LAI, m² m^-2) [18,19,20], limiting their ability to represent actual wind–canopy momentum exchange in heterogeneous agricultural fields [21,22].

This study tests a contrasting hypothesis: in a row-crop system, that d_o responds dynamically to wind–canopy interactions and within-field structural heterogeneity. Accordingly, the main aim of this work is to develop and evaluate dynamic d_o models that integrate UAS-based high-resolution crop canopy structure and surface thermal information with wind characteristics.

To achieve this aim, a distributed canopy height model (CHM) was derived from a UAS-based elevation model, resolving effective and ineffective canopy areas. As part of this effort, a pseudo surface temperature model (PSTM) was developed using PlanetDove mini-satellite multispectral PlanetScope imagery and ground-based (point) thermal infrared (TIR) surface temperature measurements. These datasets were combined with u, wind direction relative to crop rows direction (${\mathsf{\theta }}_{\mathrm{rel}}$), canopy porosity (${\mathsf{\Phi }}_{\mathrm{dp}}$), vegetation fraction(f_veg), and surface–air temperature difference (ΔT), to estimate d_o dynamically within the EC upwind footprint area.

2. Materials and Methods

2.1. Study Area

The field research study (called Bowen Ratio) was conducted at the Limited Irrigation Research Farm (LIRF), operated by the United States Department of Agriculture, Agricultural Research Service (USDA-ARS), near Greeley, Colorado, USA (40.4463° N, 104.6370° W; Figure 1). The experimental field dimension was 184 m long (N-S) by 220 m wide (E-W), at an elevation of approximately 1425 m above mean sea level.

The maize field was divided into two equal sections based on irrigation management, with the western half managed under full irrigation (FI) and the eastern half under deficit irrigation (DI). Both treatments were irrigated using a subsurface drip (SSD) irrigation system, with drip lines installed 0.23 m below the soil surface and emitters spaced at 0.30 m along the irrigation lateral lines.

The site is located in a cold semi-arid climatic region (BSk) according to the Köppen climate classification [23]. Soil properties across the field were spatially variable, with bulk density ranging from 0.95 to 1.60 g cm⁻³, corresponding to textures from loamy sand to silty clay loam based on the USDA soil texture–bulk density relationships [24]. The mean soil bulk density within the effective root zone was 1.16 g cm⁻³, and the dominant soil texture was classified as sandy clay loam.

Maize (Zea mays L.), drought tolerant Pioneer p9998Q Aquamax variety, was planted on 3 May 2024, with seedling emergence observed on 16 May 2024. The final plant population density was 86,487 plants per hectare. Data collection began in mid-June 2024 and continued through October 2024.

2.2. Instrumentation and Data Collection

A spatially distributed sensor network was deployed across the experimental field to collect micrometeorological, radiative, and canopy structural data (Figure 2). The instrumentation consisted of meteorological stations, such as EC towers and EB stations, along with ground-based surface reflectance measurements, satellite observations, and UAS surveys.

2.2.1. Eddy Covariance (EC) Measurements

Two EC towers were installed, with EC1 located in the DI zone and EC2 in the FI zone. Each EC system -based collected data was used to produce average H and LE heat fluxes over 30-min intervals. Three-dimensional wind components were measured at 10 Hz using a Gill WindMaster PRO ultrasonic anemometer mounted at 3.5 m above ground level. An LI-7500DS open-path gas analyzer (LI-COR Biosciences, Lincoln, NE, USA) measured high-frequency fluctuations of water vapor and CO₂. Canopy surface temperature was monitored using an infrared thermometer (IRT) mounted at each EC tower and EB station.

2.2.2. Ground-Based Canopy and Radiation Measurements

Spatial variability of maize canopy structure and spectral properties was monitored using a handheld multispectral radiometer (MSR5, CropScan Inc., Rochester, MN, USA). The instrument was mounted on a telescopic pole for nadir-view observations, at about 2.6 m, with an approximate 1-m ground footprint. Measurements were collected weekly under clear-sky conditions (R_s/R_so > 0.8) around local noon, where R_s is measured incoming shortwave radiation (W m_-2) and R_so is calculated clear sky incoming shortwave radiation (W m^-2), following the procedure in [25]. MSR data were obtained at 22 fixed stations per field, referred to as MSR stations and spaced approximately 15 m apart. Crop height (H_c) was measured with a tape measure at the same MSR5 stations and sampling days.

Four EB fixed stations were installed to monitor net radiation (R_n). Two stations were located in the FI zone (ST 1.1 at 40 m and ST 1.2 at 80 m from the northern boundary), and two in the DI zone (ST 2.1 at 40 m and ST 2.2 at 80 m). Each station was equipped with a net radiometer (CNR-01, Kipp & Zonen, Delft, The Netherlands), and two IRTs (SI-111, Apogee Instruments, Logan, UT, USA) oriented in nadir and oblique angles.

Surface Aerodynamic Temperature (SAT) towers were installed near the northwest corner of each irrigation treatment. Each tower was equipped with air temperature and relative humidity sensors (HMP45C, Vaisala, Finland) housed in twelve-plate radiation shields, and cup anemometers (03101-L, Wind Sentry, Campbell Scientific Inc., Logan, UT, USA) to measure horizontal wind speed. Air temperature (T_a, °C) and relative humidity (RH, %) were sampled every 5 seconds, and 15-min average values were recorded. An additional four-component net radiometer (CNR4, Kipp & Zonen, Delft, The Netherlands) was installed in the DI zone near ST 2.2, approximately 10 m south.

All sensor locations were surveyed using a handheld GPS to record latitude, longitude, and elevation. Meteorological data from the on-site Colorado Agricultural Meteorological Network (CoAgMet) station Greeley04 (GLY04) included solar radiation, wind speed and direction, air temperature, precipitation, and relative humidity.

2.2.3. Satellite and UAS-Based Data Observations

Multispectral satellite imagery were obtained from the Planet Dove CubeSat constellation [26], which provides near-daily surface reflectance at 3 m nominal spatial resolution in the red, green, red, blue (RGB), and NIR bands.

High-resolution aerial imagery were acquired using a UAS operated by the USDA-ARS Water Management and Systems Research Unit in Fort Collins, CO, USA. Flights were conducted at approximately 120 m above ground level, with nine cloud-free missions completed during the growing season. The UAS carried a MicaSense RedEdge-MX multispectral camera (MicaSense Inc., Seattle, WA, USA) and a thermal infrared camera. Multispectral imagery were collected with approximately 80% along-track and 70% cross-track overlap, yielding a ground sampling distance of approximately 5 cm.

2.3. Calibration and Interpolation of PlanetScope Imagery

PlanetScope multispectral imagery were processed using QGIS 3.40 ‘Bratislava’ and Python to generate spatially and temporally consistent surface reflectance datasets (layers) for the study area. Images were clipped to the area of interest (AOI), cloud-masked, and scaled to surface reflectance units. Reflectance values for the red, green, blue, and NIR bands were extracted using 5 m radius polygons centered on ground measurement locations to ensure spatial consistency between satellite, UAS, and field observations.

To reduce temporal noise in the satellite reflectance time series, band-wise surface reflectance values were smoothed as a function of day of year (DOY) using fifth-degree polynomial functions. Then, PlanetScope imagery reflectance values were calibrated using ground-based MSR5 surface reflectance measurements collected at 44 locations. For each spectral band, linear regression relationships were developed between satellite-derived and MSR-measured reflectance and were applied on a pixel-by-pixel basis to adjust the imagery values.

Because PlanetScope imagery were not available for every DOY, due to cloud presence, linear interpolation was used to fill temporal gaps between consecutive calibrated images. Interpolated pixel values were constrained to non-negative pixel values to maintain physical consistency. The resulting calibrated and temporally complete multispectral datasets were used to compute vegetation indices (VI) for subsequent analyses, including the normalized difference vegetation index (NDVI), optimized soil-adjusted vegetation index (OSAVI), and leaf area index (LAI). Their corresponding formulations are provided in Appendix A (Eqs. A1–A3).

2.4. Surface Temperature (T_s) and Net Radiation (R_n) Estimation

Ground-based IRT measurements collected at MSR stations were corrected for surface thermal emissivity (Ɛ_s) and background radiation effects following [27] procedures. To ensure data quality, only observations acquired under clear-sky conditions were retained. Corrected IRT-based surface temperature (T_s) was computed as in [28,29]:

$\begin{array}{c}\begin{array}{cccc}& {\mathrm{T}}_{\mathrm{s}}={\left[{\displaystyle \frac{({\mathrm{T}}_{\mathrm{I}\mathrm{R}\mathrm{T}}+273.15{)}^4-(1-{\mathsf{\epsilon }}_{\mathrm{s}}\left)\right({\mathrm{T}}_{\mathrm{b}\mathrm{g}}+273.15{)}^4}{{\mathsf{\epsilon }}_{\mathrm{s}}}}\right]}^{{\displaystyle \frac{1}{4}}}-273.15& & \end{array} (1)\end{array}$

where T_IRT is the measured surface temperature (°C), T_bg is the background (sky) temperature (−15 °C), and Ɛ_s was derived from VI.

Because PlanetScope imagery were acquired near local noon (T_rmt,12), MSR-based measurements collected at different times (~11 am to 1 pm MST) were converted to a noon time using the scaling approach developed by [30] as:

$\begin{array}{c}\begin{array}{cccc}& {\mathrm{T}}_{\mathrm{r}\mathrm{m}\mathrm{t},12}={\mathrm{T}}_{\mathrm{e}}+\left({\mathrm{T}}_{\mathrm{r}\mathrm{m}\mathrm{t},\mathrm{t}}-{\mathrm{T}}_{\mathrm{e}}\right){\displaystyle \frac{\left({\mathrm{T}}_{\mathrm{r}\mathrm{e}\mathrm{f},12}−{\mathrm{T}}_{\mathrm{e}}\right)}{\left({\mathrm{T}}_{\mathrm{r}\mathrm{e}\mathrm{f},\mathrm{t}}−{\mathrm{T}}_{\mathrm{e}}\right)}} & & \end{array}\left(2\right)\end{array}$

where T_rmt,t is the measured canopy temperature at time t, T_ref,12 and T_ref,t are reference-station temperatures at noon and time t, and T_e is the early-morning temperature.

The corrected T_s and vegetation condition relationship was modeled using OSAVI. Surface temperature was then estimated for the entire field using an exponential relationship between T_s and OSAVI as:

$\begin{array}{cccc}& {\mathrm{T}}_{\mathrm{s}}=\mathrm{a}\cdot \mathrm{O}\mathrm{S}\mathrm{A}\mathrm{V}{\mathrm{I}}^{\mathrm{b}}& & \end{array}\left(3\right)$

where a and b are empirically derived coefficients. Separate relationships were developed for the crop growth phases based on seasonal variation in OSAVI. These relationships were applied to PlanetScope OSAVI raster imagery to generate a spatially continuous pseudo-surface temperature model (PSTM) for the growing season. Noon PSTM maps were temporally extended to half-hourly intervals between 11:00 and 13:00 using [30].

Net radiation was estimated on a pixel basis using the land surface net radiation energy balance equation [31]:

$\begin{array}{cccc}& {\mathrm{R}}_{\mathrm{n}}=\left(1-\mathsf{\alpha }\right){\mathrm{R}}_{\mathrm{s}}+{\mathsf{\epsilon }}_{\mathrm{a}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{a}}^4-{\mathsf{\epsilon }}_{\mathrm{s}}\mathsf{\sigma }{\mathrm{T}}_{\mathrm{s}}^4& & \end{array}\left(5\right)$

where R_s (W m⁻²) is incoming shortwave radiation, α is surface albedo, Ɛ_a is air thermal emissivity, T_a and T_s are air and surface temperatures (K), and σ is the Stefan–Boltzmann constant equal to $5.670\text{×}{10}^{-8} {\mathrm{W} \mathrm{m}}^{-2} {\mathrm{K}}^{-4}$. The accuracy of distributed RS-based R_n estimates were evaluated with field net radiometer R_n observations.

2.5. Soil Heat Flux (G) Estimation

Soil heat flux (G, W m⁻²) was estimated using an empirical RS–based model developed by [32], which expresses G as a function of a VI and R_n as:

$\begin{array}{cccc}& \mathrm{G}=\left[0.3811\exp \left(-2.3187\cdot \mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\right)\right]{\mathrm{R}}_{\mathrm{n}}& & \end{array}\left(6\right)$

where R_n (W m⁻²) was calculated using PSTM-based raster imagery. This G model was originally developed and validated in maize fields under SSD irrigation comparable to the conditions of the present study [33,34].

2.6. Crop Height Estimation from UAS-Derived Elevation Data

Crop canopy height was estimated using UAS–derived digital surface models (DSMs) acquired at multiple stages of the growing season. A reference bare-soil digital elevation model (DEM) was generated from the DOY 155 raster image in QGIS, when vegetation cover was minimal. Bare-soil pixel elevation values were identified using NDVI = 0.25 and were interpolated to produce a bare-soil elevation surface by producing point clouds through the inverse distance weighting (IDW) method. Crop height for each flight date was calculated as the difference between the DSM and the reference DEM as:

$\begin{array}{cccc}& {\mathrm{H}}_{\mathrm{c}\text{\_}\mathrm{DOY}} = \mathrm{DS}{\mathrm{M}}_{\mathrm{DOY}} - \mathrm{DE}{\mathrm{M}}_{\mathrm{bare}}& & \end{array}\left(7\right)$

The resulting crop height models (CHMs) provided spatially distributed estimates of canopy height across the field.

To improve the agreement between UAS-derived and field-measured H_c values (H_c,field), CHM values extracted at ground observation locations were calibrated with field measured H_c (at 44 sites) using a power-law relationship as:

$\begin{array}{cccc}& {\mathrm{H}}_{\mathrm{c},\mathrm{f}\mathrm{i}\mathrm{e}\mathrm{l}\mathrm{d}}=\mathrm{a}\cdot {\mathrm{H}}_{\mathrm{c},\mathrm{U}\mathrm{A}\mathrm{S}}^{\mathrm{b}}& & \end{array}\left(8\right)$

where a and b are empirical coefficients determined from regression analysis.

To ensure that canopy height estimates represented only effective vegetation, NDVI-based binary vegetation masks were applied to the corrected CHMs using Python. Pixels classified as non-vegetated were assigned zero height, reducing the influence of bare soil, residues, and non-crop surfaces. For dates with incomplete UAS coverage, missing canopy height areas were reconstructed using nearby CHMs with complete coverage while preserving spatial patterns. A continuous raster time series of canopy height for the crop growing season was generated using pixel-wise linear interpolation between adjacent UAS observation dates as:

${\mathrm{H}}_{{\mathrm{c}}_{\mathrm{interpolated}}} = {\mathrm{H}}_{{\mathrm{c}}_{\mathrm{prev}}}+ {\displaystyle \frac{({\mathrm{H}}_{{\mathrm{c}}_{\mathrm{next}}}-{\mathrm{H}}_{{\mathrm{c}}_{\mathrm{prev}}})}{({\mathrm{DOY}}_{\mathrm{next}}-{\mathrm{DOY}}_{\mathrm{prev}})}} \text{×} ({\mathrm{DOY}}_{\mathrm{interpolated}} - {\mathrm{DOY}}_{\mathrm{prev}}) \left(9\right)$

The resulting crop height and LAI values derived from OSAVI [35] were then used to assess three widely applied d_o formulations [18,19,20] under heterogeneous canopy conditions (from the FI and DI treatments).

2.7. Sensible Heat Flux and Surface Aerodynamic Terms Calculation

2.7.1. EC Footprint Area (Pixel) Weighting

To ensure consistency between EC tower-based heat flux observations and spatially distributed raster inputs, the EC footprint upwind area was used to compute footprint-weighted means of key surface variables (e.g., T_s, ΔT, H_c, f_veg, and ${\mathsf{\Phi }}_{\mathrm{dp}}$) for each 30-min interval. These footprint-weighted variables were used in the aerodynamic terms calculations and in the dynamic d_o modeling workflow.

2.7.2. Vegetation and Surface Inputs

Vegetation and surface-state variables required for heat flux modeling were obtained from RS–derived raster products (e.g., VIs and T_s products) and summarized over the EC footprint for each timestep, ensuring that modeled aerodynamic terms reflect the effective heat flux source area contributing to the EC measurements.

2.7.3. Aerodynamic Temperature (T_aero) Calculation

An empirical T_aero model developed by [36] was used, which is needed in the calculation of H. This is, this step links surface thermal heterogeneity (raster PSTM) to the bulk aerodynamic calculations required for H, along with the inclusion of u_* and r_ah.

2.7.4. Psychrometric and Aerodynamic Variables Calculations

Supporting atmospheric variables (e.g., psychrometric terms and atmospheric stability) were derived from the surface aerodynamic tower observations and were used to compute the atmosphere surface-layer exchange quantities.

Friction velocity, for neutral atmospheric conditions, was computed as:

$\begin{array}{cccc}& {\mathrm{u}}_{\mathrm{*}}={\displaystyle \frac{\mathrm{k}\mathrm{u}}{\ln \left(\frac{{\mathrm{z}}_{\mathrm{u}}-{\mathrm{d}}_{\mathrm{o}}}{{\mathrm{z}}_{\mathrm{o}\mathrm{m}}}\right)}}& & \end{array}\left(10\right)$

where, k is the von Karman constant (0.41), Z_u is the height (m) at which the wind speed was measured, and Z_om is the roughness length for momentum transfer (m).

Aerodynamic resistance, for neutral atmospheric conditions, to heat transfer was computed as:

$\begin{array}{cccc}& {\mathrm{r}}_{\mathrm{a}\mathrm{h}}={\displaystyle \frac{\ln \left(\frac{{\mathrm{z}}_{\mathrm{u}}-{\mathrm{d}}_{\mathrm{o}}}{{\mathrm{z}}_{\mathrm{o}\mathrm{h}}}\right)}{\mathrm{k}{\mathrm{u}}^{\mathrm{*}}}}& & \end{array}\left(11\right)$

where, Zoh is the roughness length for heat transfer (m).

Under neutral atmospheric conditions, sensible heat flux was estimated as:

$\begin{array}{cccc}& \mathrm{H}=\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}}{\displaystyle \frac{({\mathrm{T}}_{\mathrm{a}\mathrm{e}\mathrm{r}\mathrm{o}}-{\mathrm{T}}_{\mathrm{a}})}{{\mathrm{r}}_{\mathrm{a}\mathrm{h}}}}& & \end{array}\left(12\right)$

where, Cp is the specific heat of dry air (J kg^-1 K^-1), ρ is the air density (kg m^-3), T_aero is the aerodynamic temperature (K), and Ta is the air temperature (K), measured at 3.5 m.

Atmospheric stability was characterized using the Monin–Obukhov stability length theory (L, m) [37]:

$\begin{array}{cccc}& \mathrm{L}=-{\displaystyle \frac{\mathsf{\rho }{\mathrm{C}}_{\mathrm{p}}{\mathrm{T}}_{\mathrm{a}}{{\mathrm{u}}_{\mathrm{*}}}^3}{\mathrm{k}\mathrm{g}\mathrm{H}}}& & \end{array}\left(13\right)$

Atmospheric stability corrections functions followed procedures by [38] for unstable conditions and [39] for stable conditions.

2.8. EC Energy-Balance Closure Correction and Evaluation

2.8.1. EC Energy-Balance Closure (EBC) Correction

EC measurements commonly exhibit incomplete surface energy balance closure (EBC), such that the measured turbulent fluxes (H+LE) are less than the available energy (AE=R_n-G). According to [40], the EB closure mismatch is more reasonably attributed to underestimation of EC turbulent heat fluxes than to errors in AE. Because the probable error in AE was only about 5-10%, whereas the observed EC lack of EBC was typically 10–30%. In this study, the EBC ratio [(H+LE)/AE] at both stations ranged from 0.60 to 0.95, indicating that the EC system underestimated turbulent heat fluxes by 5–40%. Therefore, the EC fluxes were adjusted using the following procedure:

The simplified EB of land surfaces equation is:

$\begin{array}{cccc}& R_n=H+LE+G+\mathsf{\Delta }& & \end{array}\left(14\right)$

where ∆ is a residual. Next, ∆ was partitioned between H and LE by adopting the Bowen ratio (β) method [41] (Equation (15)), following the approach described [40], yielding the residual components H and LE (Equations (15)–(18)).

$\begin{array}{cccc}& \mathsf{\beta }={\displaystyle \frac{\mathrm{H}}{\mathrm{L}\mathrm{E}}}={\displaystyle \frac{\mathsf{\Delta }\mathrm{H}}{\mathsf{\Delta }\mathrm{L}\mathrm{E}}}& & \end{array}\left(15\right)$

$\begin{array}{cccc}& \mathsf{\Delta }=\mathsf{\Delta }\mathrm{H}+\mathsf{\Delta }\mathrm{L}\mathrm{E}& & \end{array}\left(16\right)$

$\begin{array}{cccc}& \mathsf{\Delta }\mathrm{H}=\mathsf{\beta }\times \mathrm{L}\mathrm{E}& & \end{array}\left(17\right)$

$\begin{array}{cccc}& \mathsf{\Delta }\mathrm{L}\mathrm{E}={\displaystyle \frac{{\mathrm{R}}_{\mathrm{n}}-\mathrm{G}-\left(1+\mathsf{\beta }\right)\mathrm{L}\mathrm{E}}{1+\mathsf{\beta }}}& & \end{array}\left(18\right)$

This procedure yields complete (corrected) turbulent heat fluxes that satisfy EBC while retaining the observed partitioning implied by β.

2.9. Latent Heat Flux (LE) Estimation

Remote sensing-based LE (ETa) was derived as the residual of the simplified surface EB approach (Equation (14)), after estimating R_n, G, and H.

2.10. Dynamic d_o Modeling

2.10.1. Measured Calculation

Measured d_o was inferred by using the EC-derived u_* values and solving the u_* equation for d_o. In the process, the EC footprint-weighted canopy structural and wind–canopy interaction predictors, defined in Section 2.10.2, were considered and the resulting d_o values were bounded using canopy-height-based limits as 0 ≤ d_o ≤ H_c. These bounds were used to keep dynamic d_o estimates physically plausible across growth stages and spatial heterogeneity, because d_o represents an effective height associated with wind–canopy interactions and therefore cannot be negative or exceed the canopy height.

2.10.2. Crop Dynamic Porosity and Fractional Vegetation Cover

Surface aerodynamic terms required for the dynamic d_o formulations were derived from EC footprint-weighted canopy structural and thermal variables. Canopy structure within each 30-min EC heat fluxes footprint was characterized using UAS-derived CHM products. A weighted dynamic canopy porosity was defined as:

$\begin{array}{cccc}& {\mathsf{\Phi }}_{dp}=1-{\displaystyle \frac{V_c}{V_f}}& & \end{array}\left(19\right)$

where V_c (m³) is the summed canopy volume computed from pixel-wise canopy height contributions within the EC footprint area, and V_f (m³) is the total footprint volume defined as the footprint area multiplied by the mean canopy height within that area.

The weighted fractional vegetation cover within each footprint was computed as:

$\begin{array}{cccc}& f_{veg}={\displaystyle \frac{N_{veg}}{N_{pix}}}& & \end{array}\left(20\right)$

where N_veg is the total number of vegetated pixels, and N_pix is the total number of pixels in the footprint area.

The temperature gradient between the surface and air was expressed as:

$\begin{array}{cccc}& \mathsf{\Delta }T=T_s-T_a& & \end{array}\left(21\right)$

where T_s (K) is the footprint-weighted surface temperature extracted from PSTM products, and T_a (K) is the air temperature obtained from the EC system, at 3.5 m, for the corresponding 30-min intervals.

Wind direction effects relative to crop row orientation were quantified using the angular deviation between the wind azimuth angle and the row direction, defined as:

$\begin{array}{cccc}& {\theta }_{rel}=\theta -{\theta }_{row}& & \end{array}\left(22\right)$

where θ (°) is the wind azimuth angle, and θ_row is the plant row orientation. Being θ_row = 0 when the crop is aligned with the north–south orientation. The trigonometric components of θ_rel were subsequently used to characterize wind flow alignment, with sin (θ_rel) representing cross-row flow.

These footprint-weighted canopy structural, thermal, and directional predictors were paired with concurrent wind measurements.

2.10.3. Model Formulation, Training, Testing, and Validation

Dynamic d_o models were formulated by relating EC-derived d_o (measured). The combined DI and FI dataset spanning all crop growth stages was randomly partitioned, with 80% of the observations used for model calibration and the remaining 20% reserved for independent validation. The calibrated models were subsequently tested independently for the DI and FI fields to assess robustness under varying canopy structure and aerodynamic conditions.

2.10.4. Statistical Evaluation

Model performance was assessed using standard goodness-of-fit and error metrics, including the coefficient of determination (R²), mean absolute error (MAE), mean absolute deviation (MAD), root mean square error (RMSE), normalized RMSE (NRMSE), Nash–Sutcliffe Efficiency (NSE), normalized NSE (NNSE), and the index of agreement (d_r) [42,43]. These statistical parameters were applied consistently across timesteps and footprint-weighted areas of modeled vs. reference d_o quantities.

3. Results

Irrigation was the primary water source input during the growing season, with uniform (same) applications for the DI and FI fields until before DOY 160; thereafter, DI and FI treatments were irrigated differently. Rainfall amounts were very small events (< 5 mm). Local noon conditions showed strong seasonality, with T_a of 18–37 °C and R_s of 900–1000 W m⁻² during mid-season.

3.1. Seasonal PlanetScope Reflectance Patterns and MSR-Based Calibration

PlanetScope surface reflectance images exhibited pronounced seasonal pixel value variation across all spectral bands. The FI treatment maize field consistently showed higher NIR reflectance and lower Red reflectance than the DI treatment, indicating greater canopy density and vigor. However, the observed reflectance bands contained short-term variability unrelated to crop phenology, primarily due to atmospheric effects, illumination differences, and sensor noise. Thus, a fifth-degree polynomial curve was fitted to reflectance values, for each band, extracted within 5 m radius buffers around fixed and MSR stations. The obtained downloaded Planet reflectance products yielded strong agreement with observed MSR reflectance (R² > 0.80 for FI and ≥ 0.65 for DI). Then, surface reflectance calibration curves were fitted to enable estimation of reflectance for missing observation dates, resulting in temporally continuous daily surface reflectance imagery. Calibration with ground-based MSR5 reflectance measurements improved the agreement between PlanetScope and MSR5 reflectances (R² = 0.75–0.88 across bands), yielding a spatially and temporally consistent dataset.

3.2. Performance Evaluation of PSTM and Net Radiation Models

Corrected IRT-based T_s values for Ɛ_s and background temperature effects, showed clear contrasts between irrigation treatments. In the FI field, T_s decreased from early to mid-season as spatial variability declined, and as f_veg increased, with a low interquartile range (IQR), variability of mean values, narrowing from approximately 5 to 2 °C, as the crop growth season progressed, and remaining low through the late season. In contrast, the DI field showed a higher T_s value and substantially greater spatial variability during mid-season (IQR ≈ 6–10 °C), followed by a reduction in both temperature and variability toward the late season.

A strong inverse relationship was observed between OSAVI and IRT-derived T_s during the growing season. The relationship was described using an exponential decay function as:

$T_s=a.e^{-b\cdot OSAVI} \left(23\right)$

where T_s represents the corrected surface temperature, and a and b are fitted coefficients. Stage-specific parameter values varied across crop growth phases, with the strongest relationship observed during early canopy development with R² = 0.97, followed by mid-season R² = 0.82 and late season R² = 0.80.

Using fixed stations-based sensor data, OSAVI vs. IRT-based T_s relationships, daily raster imagery for surface temperature were generated for the entire growing season. The PSTM-based images were able to capture the spatial contrasts in T_s between the FI and DI fields and the seasonal decline during canopy development and senescence. The resulting T_s products for both irrigation treatments/fields were generated at 30-min intervals between 11:00 and 13:00 local time (MST) to match the data interval of EC observations.

Modeled R_n computed using PSTM-derived T_s, showed a strong agreement with field measurements across all fixed stations. The ratio of measured to modeled Rₙ approached one, with most values ranging between 0.8 and 1.3. Net radiation was generated as raster imagery for the entire field, capturing spatial and temporal variability. Larger deviations occurred during the late season (after DOY 260), when ratios occasionally exceeded 1.5 or were below 0.7. These deviations likely reflect increased canopy senescence and partial soil exposure during the late growth stage, which enhanced surface thermal heterogeneity and reduced the reliability of vegetation-based surface temperature estimation in the PSTM approach. Despite these deviations, the late-season observations were retained to preserve the full seasonal variability in the analysis. Overall, the agreement between measured and modeled Rₙ remained high across stations, with R² ≥ 0.90.

3.3. Evaluation and Reconstruction of UAS-Derived CHM

Canopy height estimates from UAV-DSM exhibited a systematic underestimation relative to field measurements in both FI and DI fields, with larger deviations observed at higher values. Following the calibration approach described in Section 2.6, we obtained:

$H_{c,field}=1.507{*H}_{c, UAS}^{0.127} \left(24\right)$

This H_c calibration resulted in near one-to-one agreement between UAS-derived (H_c,UAS) and field-measured H_c (H_c,field), and a marked reduction in error (RMSE ≈ 0.25 m; mean bias error (MBE) ≈ 0.05 m). The spatial H_c calibration was considered very good.

3.4. Comparison of Existing d_o Models

Seasonal daily d_o estimates obtained from [18,19,20] formulations exhibited distinct temporal behavior over the growing season when calculated using H_c,field values. The [18] formulation produced a smooth, height-proportional d_o trend with minimal variability, while the [19] model showed modest additional variability within a narrow seasonal range. In contrast, the [20] formulation exhibited the largest temporal variability, particularly during early and late growth stages. These periods occurred during early and late growth stages, when H_c was low and inter-row gaps were exposed.

3.5. Dynamic d_o Models: Formulation and Performance Evaluation

Two dynamic d_o models were developed to represent sub-seasonal variability in canopy–atmosphere interaction. The first formulation expresses d_o as a function of Φ_dp, u, wind direction, and ΔT (Equation (25)). While the other (Equation (26)) uses f_veg instead of Φ_dp.

$d_o=1.72-0.361u{\mathsf{\Phi }}_{dp}-0.038\mathsf{\Delta }T-0.322\mathrm{u}\sin \left({\theta }_{rel}\right) \left(25\right)$

$d_o=1.16+0.617f_{veg}-0.039\mathsf{\Delta }T-0.359\mathrm{u}\sin \left({\theta }_{rel}\right) \left(26\right)$

Both d_o dynamic models showed comparable predictive performance. Error magnitudes were similar for both models, with an MAE of approximately 0.33 m and an NRMSE of approximately 43%. Additional performance metrics, including NSE (≈ 0.48) and d_r (≈ 0.68) for the test dataset, indicated moderate and consistent predictive skill.

3.5.6. Performance Evaluation of H Estimation

In terms of H estimation, both d_o formulations showed similar performance under DI and FI conditions. Use of Equation 25 in the estimation of H produced the following coefficient of determination when evaluated with EC H values, R² = 0.70 (DI) and 0.79 (FI). While Equation 26 produced H evaluation R² = 0.70 (DI) and 0.80 (FI). In both cases, modeled H values clustered tightly around the 1:1 line in comparison with EB-corrected EC measurements (Figure 3), indicating consistent correspondence and improved agreement, especially under FI conditions.

The agreement of H modeled vs. observations was further reflected by “d” _”r” values. The Φ_dp-based H model achieved dᵣ of 0.70 (DI) and 0.74 (FI), while the f_veg-based H model yielded “d” _”r” = 0.69 (DI) and 0.74 (FI). These values indicate strong alignment between modeled and observed H in both magnitude and temporal structure, with consistently higher agreement under FI conditions. In contrast, the established displacement height formulations (Pereira et al., 1999; Shaw and Pereira, 1982; Brutsaert, 1982) exhibited lower R², higher NRMSE, and reduced “d” _”r” values, particularly under DI conditions, indicating weaker agreement and greater scatter. However, absolute error metrics were comparable, with differences in MAE and RMSE remaining small (within 2–6 W m⁻²). A comparative representation through a heat map is shown in Table 1.

Relative error metrics further confirmed the similarity between the two dynamic formulations. For the Φ_dp-based d_o model, H estimation NRMSE was 25.13% (DI) and 28.71% (FI), compared with 24.74% (DI) and 28.39% (FI) for the f_veg-based model. NMAE values ranged from 21.41–23.40% for the Φ_dp-based d_o model and 21.04–23.24% for the f_veg-based d_o model, indicating nearly identical relative errors across irrigation conditions.

3.5.7. Performance Evaluation of Dynamic “d” _”o” Models for LE Estimation

In terms of LE, both “d” _”o” formulations showed strong correspondence with EC-based LE observations under DI and FI conditions. Both do models produced similar results for the LE estimation, R² = 0.88 (DI) and 0.93 (FI). Thus, LE values clustered closely around the 1:1 line in comparison with EBC-corrected EC measurements (Figure 4), indicating a high level of agreement under both irrigation conditions.

Relative error metrics further confirmed the similarity between the two dynamic do formulations. For the Φ_dp-based do model, LE estimation resulted with a NRMSE of 15.98% (DI) and 12.83% (FI), compared with 16.16% (DI) and 13.09% (FI) for the f_veg-based do model. NMAE values ranged from 10.64–13.89% across both irrigation conditions, indicating low relative error and limited sensitivity to irrigation regime.

The “Φ” _”dp” -based do model resulted in a LE estimated vs. observed agreement of “d” _”r” = 0.78 (DI) and 0.83 (FI), and the “f” _”veg” -based do model produced identical values, indicating strong agreement on LE model performance. In comparison, the established/existing d_o formulations exhibited lower LE dᵣ values and slightly higher NRMSE, particularly under DI conditions, despite comparable absolute error magnitudes. A comparative representation through a heat map is shown in Table 2.

3.5.8. The New Dynamic d_o Models Improvement Evaluation

Relative to established empirical formulations, both dynamic d_o models yielded measurable improvements in the estimation of H. The Φ_dp-based d_o model reduced NRMSE by 15.6% under DI and 21.1% under FI, while the f_veg-based d_o model achieved reductions of 16.9% (DI) and 21.9% (FI). Agreement metrics showed concurrent enhancement, with d_r increasing by 0.7% (DI) and 3.08% (FI) for the Φ_dp-based d_o model and by 0.4% (DI) and 3.2% (FI) for the f_veg-based model. Improvements were also evident in variation around the 1:1 line of H estimated vs. observed values, with R² increases of 9.8% (DI) and 7.17% (FI) for the Φ_dp-based do formulation and 9.8% (DI) and 8.3% (FI) for the f_veg-based d_o formulation. Collectively, these results indicate enhanced representation of both magnitude and variability of H across irrigation regimes.

Uncertainty analysis provides context for these gains. Reported uncertainties in R_n (±6%) and G (±15%) propagate to an approximate ±10% uncertainty in AE [40,44]. As EC heat flux measurements typically underestimate turbulent fluxes by 10–30%, corresponding to 20–35% uncertainty in uncorrected H and LE. After EBC, using the Bowen ratio approach, corrected EC heat fluxes exhibited an uncertainty of approximately ±10–15% relative to the AE. Because EC-based u_* was used to calculate the ‘measured’ d_o, there is an associated uncertainty that propagates into the d_o estimation for the developed d_o models of ~17–25%.

4. Discussion

The results presented here demonstrate that the dynamic parameterization of d_o directly improves the representation of canopy–atmosphere exchange processes, and consequently the estimation of H and LE. The Φ_dp- and f_veg-based d_o models consistently outperformed static/fixed crop-height and LAI-based d_o formulations because newly modeled d_o explicitly accounted for temporal and spatial variability in canopy density, wind direction, and surface–air temperature gradients, processes that strongly control surface-atmosphere aerodynamic interactions, including r_ah.

The static crop height-based d_o model [18] produced larger errors because d_o was scaled only with H_c and therefore it could not adjust to changes in canopy porosity, spatially variability, directional wind-canopy interactions, and stress effects. LAI-based d_o formulations [19,20] partially improved the heat fluxes estimation performance by incorporating canopy density, but remained insensitive to within-season structural changes and footprint-weighted canopy heterogeneity. In contrast, the Φ_dp- and f_veg-based d_o formulations dynamically adjusted evolving canopy structures within the effective heat flux footprint area, leading to reduced relative errors for both H and LE.

The improvement in heat flux estimation was most evident under heterogeneous canopy conditions, where dynamic d_o reduced systematic bias and enhanced energy partitioning consistency. This result was reflected in improved closure of the surface EB, with modeled H+LE closely matching AE across irrigation treatments. The integration of pseudo-surface temperature and UAS-derived canopy height further strengthened this response by providing spatially resolved thermal and structural inputs, allowing the dynamic d_o models to better represent field-scale variability in aerodynamic forcing.

5. Conclusions

This study demonstrates that the dynamic parameterization of d_o, using maize canopy porosity and fractional vegetation cover, from UAS RS, improves heat flux estimation accuracy relative to commonly used d_o approaches. The dynamic d_o models captured spatial and temporal variability in row crop structure and microclimate, achieving higher agreement and lower relative error for both H and LE across contrasting irrigation conditions. Key outcomes include:

(a)

reduction of H and LE estimation errors by up to ~20% compared with fixeheight do formulations,

(b)

improved EBC, with modeled H + LE closely matching AE, and consistent d_o model performance under both deficit and full irrigation, supporting the robustness of the dynamic d_o approach for heterogeneous agricultural canopies.

It is hypothesized that further d_o model improvements can be achieved by enhancing early-season canopy characterization through high-frequency UAS multispectral (RGB) observations, to better capture bare-soil reference conditions and initial canopy development. Extending calibration and evaluation of the Φ_dp- and f_veg-based d_o models across additional crop structure types and climatic regions is recommended to assess transferability and scalability. Incorporating multi-year datasets and full heat flux-footprint area representations would further strengthen model robustness for operational evapotranspiration and energy-balance applications.

Appendix A.1

NDVI was calculated as indicated by Equation (A1):

NDVI = (NIR – RED) / (NIR + RED)DakonggeDakonggeDakonggeDakonggeDakonggeDakongge(A1)

OSAVI was developed by [45], indicated by Equation (A2):

OSAVI = 1.16 (NIR – RED) / (NIR + RED + 0.16)DakonggeDakonggeDakonggeDakongge(A2)

LAI was calculated by the relationship proposed by [35], indicated by Equation (A3):

$\begin{array}{c}LAI=0.263\times exp(3.813\bullet OSAVI) \left(\mathrm{A}3\right)\end{array}$