3.1.

In Situ Measurements

Study flux towers consist of a net radiometer instrument and several soil heat flow transducers with soil temperature/moisture sensors for estimating the available energy. Sensible and latent heat fluxes are measured by eddy covariance (consisting of a sonic anemometer and a fast-response open-path $\text{water vapor}/{\mathrm{CO}}_2$ sensor) at all study sites. Soil heat fluxes are measured with soil heat flux plates and are corrected for heat storage above the plate. Thus, relevant components of the surface energy balance, including measurements of net radiation, soil heat flux, sensible and latent heat fluxes, are provided at all study flux towers. In addition, ancillary meteorological observations (air temperature, wind speed, vapor, atmospheric pressure, precipitation, etc.) are available at the four study sites. Measurements at the grassland and shrubland upland sites (i.e., Kendall and Santa Rita) are made at 2 to 3 m above the local terrain. Measurements at the riparian sites (i.e., Charleston and Lewis Spring) are made at 1.5 to 2.5 m above mean canopy height. All sites have data recorded at 30-min intervals. When the daily average is calculated, days with less than 75% of data available are excluded.

3.2.

Satellite Observations (MODIS)

The Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, with 36 spectral bands (20 reflective solar and 16 thermal emissive bands), provides unprecedented information regarding vegetation and surface energy, which are critical for developing a remotely sensed ET model.³⁴ The proposed MODIS-based ET model requires a total of 15 variables obtained from MODIS atmospheric and land surface products for the calculation of net radiation and ground heat flux, as well as vegetation index (VI) and land surface temperature (LST) products for the development of the EF (Table 2). The required Aqua MODIS data are extracted from the NASA EOS Gateway ( https://wist.echo.nasa.gov/api/) and the Land Processes Distributed Active Archive Center (LP DAAC) ( http://modis.gsfc.nasa.gov) with a standard hierarchical data format (HDF). The albedo data used in the current study is a Terra and Aqua combined product (MCD43A1). For simplicity, we utilize a solar zenith angle equal to local solar noon and an optical depth of 0.2 as default values in calculating the actual (blue-sky) albedo based on a known black-and-white sky albedo ( http://daac.ornl.gov/MODIS/MODIS-menu/MCD43.html). Daily (MYD11A1) and 8-day composite (MYD11A2) LST products were obtained at a 1-km spatial resolution. Both daytime and nighttime products from MYD11A2 are acquired in order to calculate the difference of daytime and nighttime of LST. The VI was obtained using 16-day composite datasets with a 250-m resolution from both Terra (MOD13Q1) and Aqua (MYD13Q1). The MODIS composite technique (constrained-view angle-maximum value composite) includes determination of the two greatest VI values for each pixel per 16-day composite interval. The VI value observed with the nearest to nadir view is then selected for inclusion in the composite product.³⁵ An 8-day phasing in the production of the 16-day composites between Terra and Aqua allows the generation of a combined 8-day time series of VI data. We note that the compositing process may result in phenological changes going undetected, particularly where localized short-duration, high-intensity convective thunderstorms are characterized.³⁶

Table 2

Summary of MODIS products used in this study with relevant spatial and temporal information.

Since enhanced vegetation index (EVI) is optimized to improve the vegetation signal and reduce soil background influence, we utilize EVI instead of normalized difference vegetation index (NDVI) in the current study.³⁷ We note that the selection of VI may affect the shape of the final EVI-LST triangle; however, the EVI-calculated triangle method has been shown to perform as well as, or better than, than NDVI-calculated triangle method in most cases.³⁸

4.1.

Estimation of Available Energy ($Rn-G$)

To estimate the daily available energy ($Rn-G$), which in combination with EF returns the daily actual ET, we propose a satellite-based stand-alone methodology.

The net radiation for clear sky conditions is estimated according to

Upward longwave radiation is expressed using the Stefan–Boltzmann equation:

A parameterization scheme derived by Brutsaert⁴⁰ is chosen to estimate downwelling longwave radiation. This methodology is chosen based on the comparison among eight longwave parameterization schemes.^40–47

In the Brutsaert scheme, air emissivity (${\epsilon }_{\mathrm{a}}$) is determined by water vapor pressure ($e_0$) and air temperature ($T_{\mathrm{a}}$). Water vapor pressure is obtained from MYD05, and air temperature is retrieved through Eqs. (6) and (7).

Under cloudy conditions, we use a formulation of downward shortwave radiation proposed by Bisht and Bras.⁴⁸

The downward longwave radiation for cloudy pixels is estimated as proposed by Bisht and Bras.⁴⁸

To estimate the upward longwave radiation for cloudy-sky pixels, we obtain the surface temperature ($T_{\mathrm{s}\_06}$) from MYD06 to substitute for the MYD11 data product, and surface emissivity (${\epsilon }_{\mathrm{s}\_\mathrm{A}2}$) is obtained from MYD11A2, which is an 8-day composite product.

Using Eqs. (1)–(8), we can now calculate net radiation estimates for all sky conditions at the satellite overpass (instantaneous) time. Finally, we convert the net radiation estimates from the instantaneous time to a daily average net radiation estimate based on a sinusoidal function, assuming that net radiation begins to rise (become positive) at sunrise and declines at sunset.^23,49,50

Soil heat flux ($G$) is estimated from the net radiation estimates and the MODIS EVI, as proposed in the Kim and Hogue²³ methodology. However, due to the relatively poor performance of the ground heat flux scheme noted in Kim and Hogue²³ at the semi-arid study site, we incorporate a modified ground heat flux scheme into our previous algorithm for application in the Arizona region. Given that empirical schemes have strong sensor dependency, we updated the previously formulated Moran and Jackson² relationship used in our prior work. We also evaluate the performance of various MODIS-based VI (NDVI, VI, and EVI) in approximating flux tower observations of ground heat flux for the study sites. Huete et al. demonstrated that EVI is optimized to improve the vegetation signal and reduce the soil background influences in the semi-arid regions.³⁷ EVI also provided the best performance in our study, with optimized parameters of 0.22 and $-1.4$, among the vegetation indices tested. The revised ground heat flux scheme is developed as follows:

4.2.

Estimation of Evaporative Fraction

The EF is defined as the ratio of LE to available energy.¹⁵

We utilize a model developed by Wang et al., which derives EF through interpretation of the Priestley–Taylor parameter $\alpha$ using the trapezoidal distribution of a $\mathrm{\Delta }T$ (aqua day-night temperature difference)—NDVI spatial variation.¹⁷ Again, in the current study, we incorporate EVI instead of NDVI in our model. We also chose to utilize the temporal variation of $T_{\mathrm{s}}$, rather than $T_{\mathrm{s}}$ itself, to avoid a significant absolute bias from the MODIS LST product as noted in previous studies.^49,51 EF is formulated as

Due to the instability of air temperature retrieval from the MOD07 product and no significant difference (less than 5%) between the use of air and surface temperatures reported by Wang et al., surface temperature is used to estimate $\mathrm{\Delta }$ instead of air temperature.¹⁷

At each pixel, ET consists of surface evaporation from a bare soil and transpiration from vegetation. For a given EVI, variations of surface temperature occur primarily as a result of stressed surface soil water content. At $\mathrm{\Delta }T=0$, on the cold edge, bare soil and vegetation will evaporate at the potential rate, implying that ET is not limited by soil moisture. As $\mathrm{\Delta }T$ increases, evaporation from bare soil decreases and become negligible on the warm edge.

Based on the above interpretation of the ($\mathrm{\Delta }T$ and EVI) space, we can derive

From Jiang and Islam,¹² ${\alpha }^{\max }(=1.26)$ and ${\alpha }^{\min }$ can be calculated as¹²

By substituting ${\alpha }_{\mathrm{i}}$ from Eq. (15) into Eq. (12), we obtain EF for any point within the boundary of the triangle domain as

Equation (17) can also be rewritten as

Since a range of LST and vegetation indices are required to correctly represent the triangle boundary, we employ an 8-day composite LST (MYD11A2; the average LST of all cloud-free data in the compositing window) which minimizes cloud obstacles in determining $\mathrm{\Delta }{T}^{\max }$ and $\mathrm{\Delta }{T}^{\min }$. However, $\mathrm{\Delta }{T}_{\mathrm{i}}$ is derived from the daily LST (MYD11A1) to represent the LST at the satellite overpass. Thus, Eq. (18) can ultimately be expressed as

Finally, daily LE is derived from Eq. (20), given estimated net radiation ($Rn$) [Eq. (9)], ground heat flux ($G$) [Eq. (10)] at a daily basis and evaporative fraction [Eq. (19)] with the assumption of constant EF during a day.⁵²

5.1.

Validation of Net Radiation

The initial net radiation model developed in the Kim and Hogue²³ methodology incorporates the cloud product from MOD08 ($1\times 1\mathrm{deg}$) using regional fitting coefficients for downward shortwave radiation under cloudy conditions; however, this study uses the cloud product from MOD06 ($1\times 1\mathrm{km}$ and $5\times 5\mathrm{km}$) and does not require regional calibration. Bisht and Bras⁴⁸ discussed this method in detail and validated their algorithm against a range of flux towers throughout the United States. Although a semi-arid site (Desert Rock in Nevada) is included in the Bisht and Bras⁴⁸ validation, the single semi-arid/dry site was not sufficient to provide strong conclusions on the performance of the net radiation model for semi-arid regions. Thus, before application of the proposed ET model, net radiation estimates derived in the current study are further evaluated against observed net radiation from the flux tower sites.

First, we present results of the comparison between estimated instantaneous downward shortwave radiation for all sky conditions against ground-measured instantaneous solar radiation. For clear sky conditions, the shortwave radiation is estimated using the PS model from Kim and Hogue,²³ but for shortwave radiation under cloudy sky conditions, we calculate shortwave radiation using the Eq. (4) suggested by Bisht and Bras,⁴⁸ along with a regional parameterized scheme from the Kim and Hogue²³ methodology. In the comparison during clear sky days, shortwave radiation estimates have bias errors that range from 1 to 5%, and correlations are on the range of 0.94 to 0.98 across all study sites (Fig. 2). Although the Bisht and Bras⁴⁸ shortwave radiation scheme has significant advantages in that it is independent from ground-based measurements, the scatter plot between the estimated and observed quantities (Fig. 2) indicates that the result from the Bisht and Bras⁴⁸ model tends to be more scattered than the shortwave radiation estimates derived from Kim and Hogue²³ across all study sites. Bias errors in instantaneous shortwave radiation obtained from the Kim and Hogue²³ method range between 12 and 17% of the mean value, while errors from the Bisht and Bras⁴⁸ method were slightly higher, between 17 and 22% across the study sites. Root mean square errors (RMSEs) reported in Kim and Hogue²³ were also lower than those reported in Bisht and Bras.⁴⁸ Our results are also similar to a previous study, in which GOES-based solar insolation estimates have errors ranging from 5 to 10% and 15 to 30% for clear and all sky conditions, respectively.⁵³

Fig. 2

Comparison of downward shortwave radiation observations and the MODIS-derived shortwave radiation using Bisht and Bras,⁴⁸ referred as BB10, and Kim and Hogue,²³ referred as KH08, at four flux tower sites in Arizona. Both the correlation coefficient ($R$) and root mean square error (RMSE) are shown.

The downward longwave radiation scheme,⁴¹ as mentioned in Kim and Hogue,²³ showed larger uncertainty in semi-arid regions.²³ Therefore, we compared several existing algorithms that use more readily available meteorological observations such as air temperature and humidity.⁵⁴ Although these simpler algorithms may have more uncertainty relative to the more complex methods, they may be more useful for a variety of disciplinary satellite applications. Eight longwave parameterization schemes^40–47 were assessed at the Santa Rita site, which is the only site that provides longwave radiation observations.^40–47 The Brutsaert,⁴⁰ Garratt,⁴⁵ and Keding⁴⁶ algorithms show the best overall performance in terms of correlation and RMSE (Fig. 3). Among them, the Brutsaert⁴⁰ algorithm showed the best performance with a correlation of 0.94 and an RMSE of $30\mathrm{W}/{\mathrm{m}}^2$. Due to the limited observed longwave radiation data available, we are not able to perform a more comprehensive assessment of downscaling algorithms across diverse biomes including deserts, semi-arid rangelands, and riparian corridors. However, based on the results noted above, we incorporate the Brutsaert⁴⁰ longwave radiation scheme into our net radiation model for our study area.

Fig. 3

Comparison of tested downward longwave radiation schemes at the Santa Rita site.^40–47 Both the correlation coefficient ($R$) and RMSE are shown.

Comparison between the observed net radiation and the modeled net radiation at satellite overpass time [Fig. 4(a)] tends to be scattered due to the cumulative effect of bias in both the shortwave and longwave radiation schemes. However, this bias was mostly related to the underestimation in the shortwave radiation. Although not shown here, the MODIS-derived upward longwave radiation is slightly overestimated (up to $52\mathrm{W}/{\mathrm{m}}^2$ for the Santa Rita site) due to the positive bias of the LST 1-km product MODIS (MYD11). Previous studies note that the generalized split-window LST algorithm used in the 1-km MODIS LST retrieval overestimates LST in semi-arid and arid regions.^51,55,56 The overall RMSE of net radiation derived from the MODIS products ranges from 153 to $234\mathrm{W}/{\mathrm{m}}^2$. The largest bias is observed at the Charleston site, where the modeled net radiation is, on average, $197\mathrm{W}/{\mathrm{m}}^2$ underestimated compared with the observations during the MODIS overpass time. Derived daily scatter plots for each site from the sinusoidal model are shown in Fig. 4(b). Results from the daily net radiation comparison are almost similar or slightly improved relative to the instantaneous net radiation comparisons in terms of %Bias with 32 and 20% at the Lewis Spring and Kendall sites, respectively, in the daily comparison. The RMSE error at Kendall is $69\mathrm{W}/{\mathrm{m}}^2$ and 97 to $116\mathrm{W}/{\mathrm{m}}^2$ for the remaining flux tower sites. Correlation coefficients from 0.65 to 0.69 are found at the daily scale. Due to the lack of detailed spatial information around the flux towers, it is difficult to explicitly identify sources of errors and uncertainties. However, we note that significant uncertainty may result from the obvious scale differences between the MODIS-based value and in situ data. Several model inputs (i.e., MOD06 cloud fraction and surface temperature and MOD11 surface temperature) are relatively coarse (1 to 5 km), which ultimately adds uncertainty and reduces accuracy when evaluating at the flux tower scale. In addition, topographic effects and uncertainties from instrumental measurements and remote-sensing retrieval are also a factor.⁵⁷ Relative to results from the semi-arid site (Desert Rock, Nevada) in the Bisht and Bras⁴⁸ study, which reported an RMSE of $41\mathrm{W}/{\mathrm{m}}^2$ and a correlation ($R$) of 0.88, results at our study sites show relatively higher errors. However, a recent study by Tang et al. noted an RMSE of 57 and $84\mathrm{W}/{\mathrm{m}}^2$ and a correlation of 0.50 and 0.34 for MODIS-based net radiation model over the Kendall and Audubon Research Ranch flux tower sites during satellite overpass time in clear days.⁵⁷

Fig. 4

Comparison of net radiation observations and the MODIS-derived net radiation (a) at satellite overpassing time and (b) daily estimate. Correlation coefficient ($R$), RMSE, bias, and mean of observation (mean) are also presented.

5.2.

Sensitivity of Evaporative Fraction Model

Before validation of the final ET product, we evaluate the sensitivity of the triangle relationship (i.e., LST and EVI responses). Since the trapezoidal domain may be different during the growing (mid-May to early September) or senescent seasons, we separate results and explore triangle development during the two distinct periods. Figure 5(a) shows the relationship between the EF and the difference between daytime and nighttime LST from the Aqua platform, and Fig. 5(b) shows the correlation between the EF and EVI at individual sites. The observed EF is calculated from observed (flux tower) net radiation, ground heat flux, and latent heat flux through the relationship $\mathrm{EF}=\mathrm{LE}/(Rn-G)$. First, correlations between EF versus LST and EF versus EVI are very different during the growing and nongrowing seasons. For all sites, the correlations are significantly higher during the growing season than during the nongrowing season, due to the fact that EF is closely related to root-zone soil moisture during conditions of extreme soil water deficit (nongrowing season), which is difficult to represent only through either LST or EVI variables. But, in general, LST appears better correlated with EF than EVI during the nongrowing season, because the influence of soil moisture is more dominant when vegetation is small. However, during the growing season, the correlation of both LST and EVI with EF are reasonably high, ranging in magnitude from 0.63 to 0.86 and 0.50 to 0.86, respectively, indicating that LST and EVI can be used to parameterize EF. More importantly, it is clear that the slopes in the scatter plots (Fig. 5) are different for the upland grass (Kendall), savanna (Santa Rita), and riparian sites (i.e., Charleston and Lewis Spring), indicating that the influence of LST and EVI relative to EF varies from site to site. Given these results, we note that a single-triangular LST-EVI relationship may not be satisfactory for all of our study sites. Since riparian vegetation is distinctly different from the upland vegetation (density and species and better access to water compared with upland vegetation), the dominant factors affecting ET patterns are different. Ultimately, precipitation may be a more appropriate variable for predicting ET at these dry upland sites, as noted in Nagler et al.⁵⁸

Fig. 5

Regression equations for (a) evaporative fraction (EF) versus land surface temperature (LST) and (b) evaporative fraction (EF) versus EVI at the four flux tower sites. Data points are MODIS-derived 8-day composite LST (MYD11A2) and 8-day composite EVI (MOD13Q1 and MYD13Q1). Correlation coefficient ($R$) and regression line (gray solid and dashed line) are presented.

We evaluate the EF derived from the triangle method against observed (flux tower) EF [Fig. 6(a)]. The EF comparison demonstrates that the performance of the proposed triangle model in riparian areas is reasonably good, but the comparison is not as good at the semi-arid rangelands or upland sites. Particularly, the lowest correlation occurs for the Kendall site, which has the smallest vegetation cover and limited water supply, which ultimately results in extremely low EF even during growing season. Comparison at the Savanna mesquite (Santa Rita) site shows a relatively good correlation (0.65) during the growing season. At the riparian sites, the correlation varies from 0.62 to 0.80 and mean absolute error (MAE) varies from 0.10 to 0.15 (unitless) during the growing season. However, the triangle-derived EF during the nongrowing season is generally overestimated at all sites. Similarly, results reported by Tang et al. showed that the triangle method significantly overestimates EF when compared with observed EF from in situ measurements from May to September during 2004 to 2007 at two flux tower sites in Arizona.⁵⁷ Another study by Wang et al. showed that the correlation coefficients of EF ranged from $-0.081$ to 0.89 at a range of sites over the southern Great Plains with data collected from May to October in 2004.¹⁷ They also presented that EF is nonuniform regardless of the vegetation uniformity having low correlations when the mean soil moisture content is low.

Fig. 6

(a) Scatter plot of EF for the 8-day composite and (b) the latent heat (LE) at the daily scale at four flux sites. The correlation coefficient ($R$), RMSE, mean absolute error (MAE), and mean of observation (mean) are also presented for growing and nongrowing seasons.

5.3.

Validation of Evapotranspiration

Based on the net radiation estimates (Sec. 5.1) and the triangle-derived EF (Sec. 5.2), we calculate the actual ET (LE) and validate against the ground-based LE at the four study sites [Fig. 6(b)]. Evaluating our MODIS-derived net radiation results against observed net radiation from flux towers for all sites, we hypothesize that the MODIS-derived net radiation contributes significant uncertainty or bias in the ET model. Compared with the result from the EF [Fig. 6(a)], the estimated daily LE is slightly improved for all sites, likely due to the increased number of days used in the comparison. Particularly, Lewis Spring has the highest correlation (0.82) and low RMSE ($36.67\mathrm{W}/{\mathrm{m}}^2$). Results from Charleston also show similar values (0.81, $51\mathrm{W}/{\mathrm{m}}^2$ and $55\mathrm{W}/{\mathrm{m}}^2$ for correlation, RMSE, and MAE, respectively). Overall, the LE rates for the riparian areas are consistently higher and show lower variability, reflecting the typical features of riparian LE and the influence of a consistent supply of soil water. However, we note that underestimation of LE in riparian sites during the growing season may be related to the inclusion of LE from dry, nonriparian portions within the MODIS pixels. Unlike the riparian sites, the LE at the upland sites tends to be more strongly tied to precipitation events. During the summer monsoon season, the variability of LE is especially dependent on precipitation. We also observe that the LE is moderately correlated with precipitation (correlations of 0.58 and 0.63 for Kendall and Santa Rita sites, respectively). In addition, both of the upland sites (Santa Rita and Kendall) are sparsely vegetated and contain more bare soil exposure. Hence, the patchy distribution of vegetation at the upland sites likely causes a higher level of noise due to the scale difference between the satellite-derived VI ($\sim 250\mathrm{m}$) and the flux tower (footprint).

It is important to note the additional degree of uncertainty in this analysis due to the accuracy of eddy covariance systems and our comparison to these estimates. From previous work in this region,⁵⁹ LE measurements from eddy covariance tower are noted to be systematically underestimated due to the lack of energy balance closure.³⁰ In the current study, we applied a Bowen ratio correction method (BR method), in which the available energy is repartitioned into H and LE by conserving their ratio to enforce closure in the energy components.⁵⁷ We observed that the correct LE is consistently greater than the observed LE (MAE: 6 to $7\mathrm{W}/{\mathrm{m}}^2$ for Kendall and Santa Rita). These differences were much larger at both riparian sites showing about RMSE of $30\mathrm{W}/{\mathrm{m}}^2$ and MAE of 19 to $22\mathrm{W}/{\mathrm{m}}^2$, which is about 18 to 26% of the mean. Our results are similar to previous studies,³⁰ noting that LE is systematically underestimated by 17 to 27% from the lack of energy balance closure at these sites (i.e., Kendall and Santa Rita in Arizona).⁵⁹ Also, due to the long, narrow strips of vegetation along the riparian corridor, flux measurements at the riparian sites may not meet the required fetch (i.e., upwind area), which can be a significant source of error.⁶⁰ Performance of the proposed triangle model for an 8-day time series is highlighted in Fig. 7. The results are significantly improved over the daily results. The correlation is 0.72 on average for all four sites, and the RMSE is also decreased. The multiple relationships between ET, LST, and EVI are complex in semi-arid regions. However, the various relationships may be categorized depending on the land surface type (upland versus riparian). Water use efficiency is also variable depending on plant species, as seen at the grass-dominated Kendall site and the Mesquite Santa Rita site.

Fig. 7

Eight-day time series (a) and scatter plots (b) of observed LE and the MODIS-derived latent heat (LE) at four flux tower sites along with precipitation ($\mathrm{mm}/\text{day}$). The correlation ($R$), RMSE, MAE, and mean of observation (mean) are shown.