The CORDOVA-ET requires specific instrumentation for collecting and handling the data from the field. Figure 1 summarizes the components, which shows the sensing nodes, the pyranometer node, and the base station described below. Figure 2 shows pictures of the different type of nodes and base station.

The calculation of evapotranspiration (ET) using energy balance is based on the fundamental equation (1), where the net radiation (Rn) equals the sum of the latent heat (LE), sensible heat (H) and the soil heat flux (G).

Manipulating the equation, it is possible to isolate LE as in eq (2). Therefore, if the terms R_n , G and H are known, it is possible to estimate LE and therefore the evaporated water from the crop.

The sensible heat (H) can be calculated using equation (3), where ρ_a is the air density, c_p is the specific heat of air, T_c is the canopy temperature, T_a is the air temperature and r_aH is the aerodynamic resistance to heat transfer. The main requirement is to measure T_c , which can be obtained using infrared thermometry, as in the CORDOVA-ET. r_aH depends mainly on the canopy architecture (plant height, density, etc.) and wind speed.

The aerodynamic resistance in neutral conditions ( r' a ) can be calculated following (4), where z is the height of wind and temperature measurements (typically 2 m); d is the zero plane displacement height ( 2 3 h ), where h is the canopy height; z_om is roughness length governing momentum transfer (depends on the canopy architecture, but can be assumed as 0.123h); z_oh is roughness length governing the transfer of heat and vapor (also depends on the canopy architecture, but can be assumed as 0.1z_om ); u_z is the wind speed at z_m ; and k is von Kármán constant (0.4).

When calculating aerodynamic resistance (r_aH ), the distinction between neutral and non-neutral conditions is crucial. Neutral conditions refer to a state where the buoyancy effects of temperature differences between the canopy and the surrounding air are negligible. In neutral conditions, the temperature gradient between the two is not significant enough to induce vertical air movement or affect the aerodynamic processes involved in the exchange of heat and momentum. On the other hand, non-neutral conditions arise when there are substantial temperature differences between the canopy and the air, leading to buoyancy effects that influence vertical air movement. In non-neutral conditions, the buoyancy component becomes significant and affects the exchange of heat and momentum between the canopy and the surrounding air. Consequently, a more sophisticated approach is required to estimate r_a accurately. Different formulations and parameterizations for r_a were evaluated by (Liu et al., 2007). In this study, we used Viney (1991), that is based in a semi-empirical function for correcting r' a [equation (5)] using Richardson number Ri_B (6).

where:

and:

where g is the gravitational acceleration (0.98 ms ⁻²).

R_n can be measured with a net radiometer, which measures all the components of the incoming and outgoing net radiation as in eq (7), where Rs is the shortwave solar radiation, measured with a pyranometer; α is the albedo, which can be measured with a second pyranometer looking downwards; L ↑ and L ↓ are the outgoing and incoming longwave radiation, which usually are measured with a pyrgeometer.

Net radiometers that measure the four components for short and longwave are a standard part of the equipment of eddy-covariance systems. However, they are not commonly available in agrometeorological weather stations. Thus modeling the albedo and longwave components is required. Estimating the incoming longwave radiation (L ↓) requires air temperature and humidity, as well as a cloudiness factor that can be calculated as a fraction of actual solar radiation (Rs) and the potential irradiance (Rs_pot ). This is the standard method for FAO and ASCE formulations (Allen et al., 1998; Walter et al., 2000). However, for the sub-daily calculations, we followed a different approach based on (Crawford and Duchon, 1999) that has been tested by Berni et al. (2009). Here, the estimate of the cloudiness is not done on every time step but calculated daily, and the same cloudiness factor is applied for all the day’s observations.

The outgoing longwave radiation (L ↑) can be estimated from the actual canopy temperature, converting the temperature to radiant energy using the Stefan–Boltzmann law: L ↑ = ϵ c · σ · T c 4 , where ϵ_c is the canopy emissivity (assumed as 0.98); σ is the Stefan-Boltzmann constant; and T_c is the canopy temperature.

One factor that requires attention is the albedo α. In the absence of two pyranometers or albedometer, we are assuming α = 0.18. However, this assumption may not be correct due to the nature of the albedo since its value depends on the nature of the canopy (e.g., ground cover, senescence) and on solar elevation. This may require further investigations in the future and modeling for different crops and conditions.

G can be measured using soil heat flux plates. However, it can also be estimated using FAO56 approach where G is calculated as a fraction of R_n that is different during the daytime (R_s ≥ 0) or nigh time as expressed in (8).

With all the components of (2) known, it is possible to calculate LE. Converting LE to ET can be performed by converting energy fluxes to mass fluxes using (9):

Where λ is the latent heat of the vaporization of water (Jkg ⁻¹), which depends on the air temperature (T_a ). The units of ET_a resulting from (9) will be kgm ⁻² s ⁻¹. Since most applications require mm/day units, mass fluxes must be aggregated throughout the day and water density can be assumed as 1m ³ = 1kg.

All these processing steps have been implemented in a Python library called OpenCropLib, which is available as open-source (Jimenez-Berni, 2023).

Another Python library, called RefET, has been used to calculate reference ET (ET_o ) from the weather station data using both daily and hourly intervals.

The wheat field used for validation was very close to an official weather station on a grass field in Campus Alameda del Obispo, Cordoba (Spain): 37.86N, 4.8W. The field was sown in November 2020 and harvested in June 2021. The dimensions of the field were 220x90m, oriented in the E-W direction.

A complete CORDOVA-ET system was installed (four nodes, a pyranometer node, and the base station) in the middle of the field on February 2021, selecting for the nodes a homogeneous area with no apparent crop issues. The air temperature and humidity sensors, and infrared thermometers were installed at 2.5 m. The weather station was installed 10m next to the official weather station at 2m height.

An eddy covariance system was installed within 10m downwind of the CORDOVA-ET sensor nodes to validate the energy fluxes. The equipment consisted of (1) a three-dimensional sonic anemometer (model CSAT3; Campbell Scientific Inc., Logan, UT, USA) to measure vertical (w), horizontal (u) wind speed and sonic air temperature T_s and (2) a high-frequency open-path infra-red gas analyzer (IRGA) to measure CO2 (Ca) and water vapor (q) concentrations in air (model LI-7500, LI-COR, Lincoln, USA). The separation between the anemometer and the gas analyzer was 0.25 m during the measurement period. The footprint of eddy fluxes was calculated using the model of Schuepp et al. (1990). Then, the sensor height of the eddy covariance system was adjusted to ensure that more than 90% of the fluxes originated inside the wheat field. Canopy height was measured weekly with a ruler.

Solar panels powered all the sensors which were connected to a datalogger (model CR1000, Campbell Scientific Inc.) operating at 10 Hz. Data were processed to obtain 30-min covariances ( w ′ C a ¯ , w ′ q ¯ , w ′ T s ¯ , w ′ u ′ ¯ ). The raw measured fluxes of carbon (F_c ), water vapor (LE), sensible heat (H), and momentum (τ) were corrected using the frequency response functions derived by Moore (1986). I.e., frequency response on sensor separation, path length averaging, and signal acquisition and processing time. F_c and LE were also corrected for density fluctuations, according to Webb et al. (1980). All the calculations for the eddy covariance system were performed using the TK3 software (Mauder and Foken, 2015). To obtain the daily averages, the semi-hourly corrected fluxes were then averaged for each 24 h, including nocturnal data with friction velocity > 0.2 m/s. A problem with the solar power supply caused a power outage from 7 to 28 April, lea.

The net radiation and its components (long and short wave, each both downwards and upwards) were measured each minute with a 4-way net radiometer (model NR01; Campbell Scientific Inc., Logan, UT, USA) connected to a datalogger (model CR1000, Campbell Scientific Inc., Logan, UT, USA) which also averaged and stored the radiation fluxes every 10-min. Next to the net radiometer. An infrared temperature (IRT) sensor (model IRTS-P, Apogee, USA) was deployed close to the neet radioameter to record continuous canopy temperature for validation. This is a widely used high-quality sensor used in multiple publications (Berni et al., 2009; Krishnan et al., 2020).

The energy balance closure of the eddy covariance energy fluxes was checked by comparing H + LE and the available energy was calculated as R_n − G, where R_n is measured with the net radiometer and G is calculated using (8). To correct for potential closure error, a correction was introduced using Bowen ratio (β = H/LE), where R_n is partitioned using β using equations (10) and (11). β was calculated on daily steps using the sum of H and LE on each day. The corresponding daily β was used in the calculation of the corrected energy fluxes.

A reference weather station was located 10m from the CORDOVA-ET and the eddy covariance systems. The weather station is part of the Irrigation Agroclimatic Information System (SIAR in Spanish) (Ministerio de Agricultura, 2023), a nationwide weather station network.

The validation of the different variables has been performed on half-hourly data averaged from the raw data recorded at higher temporal resolution. The time series of selected days have been plotted for the reference station and CORODOVA-ET sensors to show the system’s capacity for tracking diurnal conditions. For each of the evaluated variables, the determination coefficient (R ²), slope, root mean square error (RMSE), and mean absolute error (MAE) have been calculated using (12) and (13), respectively. Night-time measurements (when R_s ≤ 0) have been removed for the analysis.

Where y i ^ are the predicted values and y_i are the actual values, in this case, the variables measured by the reference station or the eddy covariance.

The different components of the energy balance (namely R_n , H and LE) were calculated using the data from the CORDOVA-ET and the methodology described above. These results have been compared with the estimates from the eddy covariance before and after applying the energy balance closure correction ( Figure 8 ). Five days in March 2021 were selected to plot the daily course of these variables after the closure correction ( Figure 9 ). Finally, a summary of the average daily values shows the evolution of the fluxes for the entire validation period (except for the data gaps) ( Figure 10 ). A detailed analysis of each component follows in the next subsections.

The high cost of accurately estimating crop evapotranspiration is a significant obstacle to quantify crop ET in agronomic studies, especially in developing countries or under limited research funding. The equipment for estimating ET using eddy covariance can cost from tens to hundreds of thousands of dollars. While there have been efforts to create low-cost systems (Markwitz and Siebicke, 2019), setting up, maintaining, and processing the data from eddy covariance instrumentation requires substantial expertise. This study aims to provide an affordable solution for researchers and technicians in countries with limited access to this technology. The use of infrared thermometry for measuring plant temperature is becoming widely adopted for estimating drought stress and irrigation management using semi-empirical methods such as the crop water stress index (CWSI) or stress time accumulation, among other methods (Maes and Steppe, 2012). However, using ground-based thermography for estimating actual evapotranspiration is still rare. While recent progress has been made in the use of machine learning for estimating ET_a (Pagano et al., 2023), the use of well tested physical models with data acquired with a simpler device can still provide robust estimates of the energy fluxes and evapotranspiration, comparable to the gold standard of eddy covariance, as demonstrated in our results.

The possibility of calculating both ET_o and ET_a with the CORDOVA-ET device with results comparable to complex instruments and established methods such as the eddy covariance, opens new avenues and applications such as studying crop coefficients under different irrigation or management strategies. The possibility of using multiple nodes for measuring the canopy temperature of various fields or treatments presents a clear advantage against eddy covariance, where the prevailing winds and instrumental setup determine the observed footprint. With this instrument, the footprint of the canopy temperature measurements can be defined based on the orientation of the sensor and the sensor height, which depending on the sensor’s field-of-view establishes the area of interest. Moreover, having a scalable system with multiple measurements in the same area can also help to characterize the spatial heterogeneity of ET in non-uniform crops. This is paramount in studies focused on validating remote sensing estimates of ET, where satellite imagery’s spatial resolution aggregates pixel heterogeneity.

The CORDOVA-ET system also provides a simple solution for applications such as precise irrigation management. The real-time calculation of both ET_o and ET_a provides a direct estimate of the water use rates, which could help farmers make decisions on irrigation scheduling. It could also be applied for water stress detection, as each node provides direct measurements of canopy temperature, air temperature and humidity, allowing a straightforward calculation of crop water stress index (CWSI) in each node location. This could be single trees (Berni et al., 2009) or uniform herbaceous crops (O’Shaughnessy et al., 2012). Since the system also provides a complete set of environmental variables, applying a physical approach to CWSI is feasible instead of using empirical baselines (Agam et al., 2013).

Different studies have focused on the sensitivity analysis of input variables in applying energy balance approaches (Maes and Steppe, 2012) or estimation of reference evapotranspiration (Estévez et al., 2009). Among the critical variables, solar radiation (R_s ) is the major contributor to the available energy required for evapotranspiration. In the early versions of this instrument setup, it was clear that the measurements from a commercial weather station were not optimal due to the lack of a cosine diffuser that integrates the total radiation from the sky hemisphere. Modifying the original instrument for integrating an optical diffuser improved the sensor’s performance and provided results similar to the reference pyranometer, with a slope of 0.99 and RMSE and MAE of 42.82 and 6.29 Wm ⁻² respectively ( Figure 12 ). However, the weather station’s pyranometer, located at 2m, is hard to inspect and clean. Eventually, the adopted solution was to incorporate a commercial pyranometer with a dedicated node. This pyranometer can be installed at a standard and comfortable height for servicing, which ensures high-quality data in the longer term.

A fundamental parameter in the energy balance modeling is the albedo. The value of α=0.23 is the standard for reference ET calculation using FAO56. However, this value is relatively high for wheat, as reported by Zhang et al. (2012). If fact, it is known that albedo changes along the day as a function of solar elevation, and that it is influenced by canopy height and canopy architecture. In this study, a fixed value of 0.18 has been selected for simplicity. This fixed value can lead to inaccuracies in the estimation of the available energy, especially at times when the solar elevation is low or when the crop starts maturing, and pigment changes can lead to significant changes in reflectance and, therefore, canopy albedo. Figure 13 . A exhibits that midday albedo evolves from 0.2 to 0.1 during the validation period, showing values below 0.18 after April ( Figure 13A ). Moreover, diurnal changes in albedo ( Figure 13B ) are also very significant and asymmetric as reported in (Zhang et al., 2012). Including an albedometer as part of our system could help track changes in albedo and use measured values instead of a fixed figure. Alternatively, a model could be developed for specific crops using phenology and solar elevation as the inputs or using proxies such as spectral reflectance sensors (Cao et al., 2018).

The aerodynamic resistance (r_a ) plays a crucial role in the H equation ((3)), but direct measurement of this parameter with instruments is not feasible. Therefore, it can only be estimated using models. In our study, we employed the model proposed by Viney (1991) to estimate r_a due to its ability to incorporate canopy temperature measurements, which allows for the inclusion of buoyancy effects. It is important to note that using this semi-empirical estimation method provides a practical approach to estimating r_a without direct measurements. However, it is essential to consider the limitations of this approach. The empirical parameters used in the model were derived under specific conditions, and their applicability to different field scenarios may vary. While the canopy-air temperature difference for buoyancy correction improves the estimation of r_a , it is still necessary to validate the model’s performance under the specific conditions encountered in the field. Furthermore, other parameters influenced by canopy architecture, such as d and z_o , could also benefit from more sophisticated parameterization techniques, like the drag partition model utilizing frontal area index (Raupach, 1992), or other simplifications (Verhoef et al., 1997). These alternative approaches may provide more accurate estimations for complex canopies, including tree crops (Berni et al., 2009).

In this instrument, the canopy temperature measured with an infrared thermometer is the basis for calculating H. Here, we assume that the canopy temperature (T_c ) measured with the IR radiometer equals the aerodynamic temperature (T ₀) originally used in equation (3). However, T ₀ is not directly measurable, and the relationship between T_c and T ₀ is complex (Colaizzi et al., 2004) and depends on viewing angle, atmospheric stability, vegetation structure, and soil temperature (Maes and Steppe, 2012). While several formulations have been proposed to deal with this discrepancy (Mahrt and Vickers, 2004; Kustas et al., 2007; Kustas and Anderson, 2009; Mallick et al., 2022), the use of off-nadir T_c observations can minimize the differences between T_c and T ₀ (Huband and Monteith, 1986; Matsushima and Kondo, 1997) and it is possible to assume for close range measurements over full canopy cover. In our case, the sensor is installed at 45 degrees off-nadir, which is close to the optimal 50-70 deg reported by (Huband and Monteith, 1986; Matsushima and Kondo, 1997), while the chances of looking at the sky or soil in trees or row crops.

It is important to note that the methodology presented relies on the big-leaf assumption and requires temperature observations from pure vegetation with minimal soil contamination. This means that if the methodology is used in crops with partial cover or in situations where there are significant differences between the temperature of the canopy and soil, the results of H may be biased. Using multi-source energy balance models could enhance this situation and provide more accurate results (Anderson et al., 1997). This modular system would be able to incorporate multiple temperature measurements for canopy and soil or even multiple viewing angles.

A significant contribution of this study is that all the components required to build the sensors, acquire the data, and process the data for obtaining energy fluxes and crop ET are open-source. Offering an open-source solution for estimating ET seeks to provide a viable alternative to researchers and technicians, particularly in developing countries. This open-source option not only reduces the entry barrier in terms of cost and knowledge. It also provides a reproducible solution where the users can completely control all the processing steps and parameterization in modeling the energy fluxes. Also, this development can serve as a base for future developments for interfacing new sensors. Similarly to how we developed a specific node for interfacing a pyranometer, other researchers could implement developments for interfacing other sensors for measuring soil moisture, canopy height, or spectral reflectance, among different possibilities.