Evapotranspiration Estimates Derived Using Multi-Platform Remote Sensing in a Semiarid Region

Evapotranspiration (ET) is a key component of the water balance, especially in arid and semiarid regions. The current study takes advantage of spatially-distributed, near real-time information provided by satellite remote sensing to develop a regional scale ET product derived from remotely-sensed observations. ET is calculated by scaling PET estimated from Moderate Resolution Imaging Spectroradiometer (MODIS) products with downscaled soil moisture derived using the Soil Moisture Ocean Salinity (SMOS) satellite and a second order polynomial regression formula. The MODis-Soil Moisture ET (MOD-SMET) estimates are validated using four flux tower sites in southern Arizona USA, a calibrated empirical ET model, and model output from Version 2 of the North American Land Data Assimilation System (NLDAS-2). Validation against daily eddy covariance ET indicates correlations between 0.63 and 0.83 and root mean square errors (RMSE) between 40 and 96 W/m2. MOD-SMET estimates compare well to the calibrated empirical ET model, with a −0.14 difference in correlation between sites, on average. By comparison, NLDAS-2 models underestimate daily ET compared to both flux towers and MOD-SMET estimates. Our analysis shows the MOD-SMET approach to be effective for estimating ET. Because it requires limited ancillary ground-based data and no site-specific calibration, the method is applicable to regions where ground-based measurements are not available.


Introduction
Evapotranspiration (ET) is fundamental to understanding the regional water balance, particularly in semiarid and arid regions where ecosystem processes are often limited by the availability of water and where water loss is dominated by ET [1]. However, ET remains one of the most challenging surface fluxes to estimate due to its dependence on numerous climatological parameters as well as physical soil properties and land cover [2][3][4].
Conventional ground-based ET measurement techniques (eddy covariance, Bowen ratio) are constrained to relatively small homogeneous footprints that rarely exceed 1-2 km [5,6]. As such, they are limited in their capability to estimate fluxes on larger spatial scales due to the inherent heterogeneity of the land surface and hydroclimatological forcing. Alternatively, different methods have been proposed to estimate ET through satellite remote sensing which allows acquisition of large-scale distributed data at various spatial and temporal resolutions. These methods vary in complexity, with tradeoffs between empirical and physically-based models [7]. Methods include and 109.86° to 110.90°W, respectively. The main land cover types include shrub/scrub land (84%), evergreen forest (11%), grassland/herbaceous (2%) and areas of developed open space (1.5%) [47]. Surface elevation ranges from 830 to 2890 m. The climate is classified as semiarid, with over half of the rainfall occurring between July and September during the North American monsoon (NAM) [48]. Precipitation during the NAM is characterized by local, short-duration, high-intensity convective thunderstorms, which closely correlate with ecosystem flux responses and surface soil moisture distributions. For the study year (2013), the NAM began on 5 July and ended on 30 September [49]. Four eddy covariance flux tower sites were used for initial validation ( Table 1). Two of the AmeriFlux tower sites (Lucky Hills (US-Whs) and Kendall Grassland (US-Wkg)) are located within the USDA-ARS Walnut Gulch Experimental Watershed (WGEW). Lucky Hills (US-Whs) is surrounded by a diverse stand of Chihuahuan Desert shrubland species that dominate the area [50]. Kendall Grassland (US-Wkg) consists mainly of C4 grasses with a few scattered shrubs [51]. Vegetation surrounding Santa Rita-Creosote (US-SRC), which is located within the Santa Rita Experimental Range (SRER), is primarily mature creosote bush [52]. Vegetation at the Charleston Mesquite Woodland (US-CMW) site is predominantly dense mesquite-dominated riparian woodland with a maximum vegetation height of 10 m [53]. Vegetation water use at this site is supplemented by groundwater.

In Situ Measurements
Ground-based ET measurements were collected using the eddy covariance (EC) technique from micrometeorological towers located at each site [51]. These provide estimates of ET over footprint areas of several thousand square meters at each tower site [54] (Glenn et al., 2015). The EC method consists of three-dimensional, sonic anemometers (CSAT-3; Campbell Scientific) and open-path infrared gas analyzers (LI-7500, LI-COR) used to measure the wind velocity vector, sonic temperature and concentrations of water vapor and carbon dioxide. Data were sampled at 10 Hz Four eddy covariance flux tower sites were used for initial validation ( Table 1). Two of the AmeriFlux tower sites (Lucky Hills (US-Whs) and Kendall Grassland (US-Wkg)) are located within the USDA-ARS Walnut Gulch Experimental Watershed (WGEW). Lucky Hills (US-Whs) is surrounded by a diverse stand of Chihuahuan Desert shrubland species that dominate the area [50]. Kendall Grassland (US-Wkg) consists mainly of C4 grasses with a few scattered shrubs [51]. Vegetation surrounding Santa Rita-Creosote (US-SRC), which is located within the Santa Rita Experimental Range (SRER), is primarily mature creosote bush [52]. Vegetation at the Charleston Mesquite Woodland (US-CMW) site is predominantly dense mesquite-dominated riparian woodland with a maximum vegetation height of 10 m [53]. Vegetation water use at this site is supplemented by groundwater.

In Situ Measurements
Ground-based ET measurements were collected using the eddy covariance (EC) technique from micrometeorological towers located at each site [51]. These provide estimates of ET over footprint areas of several thousand square meters at each tower site [54] (Glenn et al., 2015). The EC method consists of three-dimensional, sonic anemometers (CSAT-3; Campbell Scientific) and open-path infrared Remote Sens. 2017, 9, 184 4 of 22 gas analyzers (LI-7500, LI-COR) used to measure the wind velocity vector, sonic temperature and concentrations of water vapor and carbon dioxide. Data were sampled at 10 Hz and statistics calculated for 30-min blocks. Previous work indicates latent heat flux (LE) from EC towers may be systematically underestimated in the region due to the lack of energy balance closure [51]. Hence, we correct eddy flux measurements for closure errors using a strategy proposed Twine et al. [55]. Observed sensible and latent heat fluxes are modified so as to sum to the available energy (R n -G) yet retain the observed Bowen ratio. Ground-based volumetric soil moisture content (m 3 /m 3 ) was estimated using Campbell Scientific CS616 water content reflectometers (Campbell Scientific, Inc.) [56]. Installation depths vary slightly between sites, with surface depths between 2.5 and 5 cm and subsurface depths (defined as root zone depth) between 12.5 and 15 cm.

SMOS Satellite Observations
SMOS data consist of soil moisture estimates from the Level 3 (L3) product created by the best estimation of soil moisture and dielectric constant based on a minimization of a data quality index (DQX), as well as through temporal and/or spatial resampling or processing (Centre Aval de Traitment des Donnees SMOS (CATDS); http://catds/ifremer.fr/). Spatial resampling creates a 25 km spatial resolution soil moisture product with an intended accuracy of at least 0.04 m 3 /m 3 [57,58]. SMOS provides an added benefit to additional soil moisture satellites such as AMSR-E, and its successor, AMSR2, by measuring soil moisture in the L-band (~1-2 GHz) microwave frequency. Frequencies in the L-band are preferred as they are more sensitive to changes in soil moisture [59], less susceptible to attenuation due to the atmosphere or vegetation compared to higher frequencies such as C-band (AMSR-E and AMSR2) [60,61], and penetrate to greater depths within the surface layer than shorter wavelengths [62]. SMOS has also been found to provide improved soil moisture estimates in southern Arizona when compared to AMSR-2 and regional flux towers [63]. SMOS has a sun-synchronous orbit with local equatorial crossing times of approximately 6:00 a.m. and 6:00 p.m. in ascending and descending nodes, respectively [64]. The ascending (6:00 a.m.) node is used in the current study because surface soil layer conditions are expected to be closest to thermal equilibrium at this time [58,65,66]. SMOS soil moisture data are discarded when the quality of retrieval is poor (DQX greater than 0.07), the soil moisture value is negative, or the retrieval has failed. Radio Frequency Interference (RFI), the contamination of microwave observations, has been reported as a nonissue in the region of interest over the same time period [63].

Moderate Resolution Imaging Spectroradiometer (MODIS) Observations
Combinations of seven variables obtained from the MODIS Terra platform are used in the current study to derive both PET and downscaled soil moisture. Products include daily land surface temperature (LST) (MOD11_L2), daily emissivity (MOD11_L2), height above geoid (MOD03), daily water vapor (MOD05), daily atmospheric profile (MOD07_L2), eight-day composite albedo (MCD43B3), and 16-day composite NDVI (MOD13A2/MYD13A2). Due to the optimization of enhanced vegetation index (EVI) in improving the vegetation signal and reducing soil background influence in semi-arid regions, we substitute NDVI for EVI [67]. The phasing of both Terra and Aqua EVI products generates a combined eight-day time series of vegetation indices. In the current study, we utilize a solar zenith angle equal to local solar noon and an optical depth of 0.2 as default values when calculating blue-sky albedo based on a known black-and-white sky albedo. All MODIS products are acquired at 1 km resolution from the NASA Reverb ECHO site (http://reverb.echo.nasa.gov/reverb) in the standard hierarchical data format (HDF) for the year 2013. As with the SMOS dataset, MODIS variables indicating poor quality, failed retrieval, or contamination were omitted prior to the calculation of subsequent parameters.

Empirical ET Model
In addition to flux tower observations, MOD-SMET is also compared to an empirically derived ET at all four study sites. While many empirical ET models exist [68][69][70][71], the current study utilizes a model introduced by Scott et al. [45] and later adjusted by Bunting et al. [46] due to its dependence on only satellite derived parameters. The model follows the form: where a, b, c, d, and e are calibrated coefficients, and LST and EVI are nighttime land surface temperature and enhanced vegetation index, respectively, acquired from MODIS products. Adjustments of Equation (1) by Bunting et al. [46] led to numerous formulations requiring LST, a normalized version of EVI, and/or precipitation, depending on the site used for calibration. Models with best performance under either a riparian or upland classification, as determined by Bunting et al. [46], were used in the current study (Equations (2) and (3)).
ET(upland) = 6.93(EVI) + 0.017(PPT) − 0.507, EVI* is a normalized eight-day composite EVI, LST is the land surface temperature, and PPT is an eight-day total precipitation (mm). In the current study, an eight-day composite nighttime LST is derived from MOD11A2, while the eight-day composite EVI is derived using both MOD13A2 and MYD13A2 as described in Section 3.2.2 and in Bunting et al. [46]. The eight-day total precipitation is obtained from each flux tower site. All models were chosen as they have been calibrated specifically for three of the four sites currently under investigation. Upland sites include Wkg, SRM, and Whs and will thus be estimated using Equation (3). Riparian site CMW will be estimated using Equation (2). The eight-day empirically derived ET will be referred to as E8-ET.

NLDAS-2 Models
NLDAS-2 provides distributed hydrometeorological products over the contiguous United States at~12.5 km spatial resolution. NLDAS-2 forcing data consist of downward shortwave radiation, downward longwave radiation, 2 m air temperature, 2 m air specific humidity, precipitation, surface pressure, and 10 m wind speed and is used to drive four corresponding land surface models (LSMs). Subsequently, each LSM produces state variables and water and energy fluxes, such as ET. LSMs within NLDAS-2 include the Noah model [34], the Variable Infiltration Capacity (VIC) model [72], Mosaic model [73], and the Sacramento Soil Moisture Accounting (SAC-SMA) model [74].
The SAC-SMA uses a climatologically based PET with seasonal variation (but no intra-monthly or inter-annual variation) [75,76]. Due to the poor representation of intra-annual variation, the SAC-SMA is not included in the current study. In addition, preliminary analysis also indicated considerable underestimation of ET from VIC model simulations (not shown). Hence, the LSMs included for comparison are the Mosaic and Noah models. Derivation of ET from Mosaic and Noah are expressed as the combination of direct evaporation from the soil surface, transpiration via canopy and roots, and evaporation of precipitation intercepted by the canopy. Major differences between model ET estimates originate from differences in defined root-zone depths, which affect the magnitude of transpiration via canopy and roots. Depths range from 1 m for short vegetation to 2 m for short trees and woody vegetation in the Noah model and 40 cm for all vegetation types within the Mosaic model. For further details, see Chen et al. [34] for the Noah model and Koster et al. [73] for the Mosaic model.

PET Estimation
PET is calculated using Kim and Hogue [41] and the Priestley-Taylor equation (Equation (4)), written as: where PET is the daily average potential ET (in flux units of W/m 2 ), ∆ is the slope of saturated vapor pressure versus temperature (kPa/K), R n is the daily average net radiation (W/m 2 ), G is the daily average ground heat flux (W/m 2 ), γ is the psychrometric constant (kPa/K), and α is an empirical factor known as the Priestley-Taylor constant. A value of 1.26 is used in the current study as it best describes ET from a variety of well-watered vegetated and water surfaces (i.e., PET) [77]. Despite the simplicity of the Priestley-Taylor approach (neglects influence of vapor deficit, primarily relying on radiation and temperature as proxies), it has been used extensively and has shown to provide reasonable estimates for agricultural and hydrologic studies [78][79][80]. Priestley-Taylor estimates have also been found to correlate to Penman-based PET estimates, a more advanced resistance-based model [81,82]. The Priestley-Taylor method is advantageous for this study because it satisfies our constraint that each formula variable can be derived from satellite remote sensing information. Instantaneous net radiation is estimated through the energy balance formulation (Equation (5)), written as: where R n,inst is the instantaneous net radiation, Alb is the surface albedo (MCD43B3), R sw is the instant downward shortwave radiation, and R lw,down and R lw,up are instant downward and upward longwave radiation, respectively. Estimation of downward shortwave radiation stems from Zillman [83] and modifications by Bisht and Bras [84]. This parameterization scheme uses near surface vapor pressure (e 0 ) (MOD07 and MOD05) and solar zenith angle (θ) (MOD03) to estimate downward shortwave radiation as follows: where S 0 is the solar constant at the top of the atmosphere (1367 W·m −2 ) and β is an empirical coefficient determined to be 0.2 by Niemelä et al. [85] and Bisht et al. [86]. Upward longwave radiation is expressed using the Stefan-Boltzman equation: where ε s is the surface emissivity (MOD11), σ is the Stefan-Boltzman constant (5.67 × 10 −8 W·m −2 ·K −4 ), and T s is the surface temperature (K) (MOD11). Downward longwave radiation is based on a parameterization scheme by Brutsaert [87] and is estimated as: where ε a is the air emissivity (determined by water vapor pressure (MOD05) and air temperature (MOD07)), and T a is interpolated air temperature (K) (MOD07). Equations (6)-(8) are used within Equation (5) to determine instantaneous net radiation at time of satellite overpass. Instantaneous net radiation estimates are converted to daily average net radiation (R n,daily ) through a sinusoidal function that assumes R n values become positive at sunrise and begin to decline at sunset (Equations (9)) [41,86,88].
Remote Sens. 2017, 9, 184 7 of 22 R n,daily and R n,inst are daily and instantaneous net radiation, respectively, t sunrise and t sunset are sunrise and sunset times obtained from the U.S. Naval Observatory, and t i is the satellite overpass time.
Ground heat flux (G) is estimated through an empirical formulation as a fraction of daily net radiation, utilizing radiometric surface temperature, surface albedo, and EVI as proposed by Bastiannssen [10] (Equation (10)).
In the current study, T s is a radiometric surface temperature derived from MOD11 in degrees Celsius, Alb is the surface albedo (MCD43B3), and EVI is the Enhanced Vegetation Index (MOD13A2 and MYD13A2) (replacing the originally utilized NDVI in Bastiannssen [10]).
Lastly, slope of the saturation vapor pressure versus temperature (∆) is calculated as: where T a is interpolated air temperature ( • C) (MOD07) and e s is the saturated vapor pressure (kPa) written as:

Soil Moisture Estimation
Daily downscaled soil moisture estimates at a 1 km 2 spatial resolution are derived through a second order polynomial regression formula relating soil moisture availability, land surface temperature (T s ), and vegetation indices [15,[42][43][44]. The regression relation, proposed by Carlson et al. [15], can be written as: where θ sfc is the estimated soil moisture, a i,j are the model fitting coefficients, and EVI* and T s * are the normalized enhanced vegetation index (MOD13A2/MYD13A2) and normalized land surface temperature (MOD11), respectively, defined as: where EVI min and EVI max are the minimum and maximum MODIS-derived EVI values determined over the study domain, and similarly, where T s,min and T s,max are the minimum and maximum MODIS derived T s values determined over the study domain. This approach, and similar variations, has been applied in many regions across the globe with promising results [42][43][44]89,90]. Specifically, the method above was shown to provide reasonable downscaled surface soil moisture estimates in the southern Arizona region by Knipper et al. [63]. Information on saturated and residual soil moisture content (θ sat and θ res , respectively) were spatially derived through STATSGO soil texture data using established pedo-transfer functions [91,92] to calculate topsoil effective saturation (θ sfc,eff ) at 1 km grid using the following equation first proposed by van Genuchten [93]: where θ sfc,eff , θ sfc, θ res , and θ sat represent the effective saturation, downscaled soil moisture estimates, and residual and saturated moisture content at 1 km pixel, respectively.

Actual Evapotranspiration
ET is derived using PET and a simple soil moisture function, f(θ) [91] (Equation (16)): where ET is the actual evapotranspiration and the soil moisture function, f(θ), is a dimensionless variable estimated by a simple linear model: where θ rz is the volumetric soil moisture content at root zone depth and θ fc is the field capacity derived spatially though STATSGO soil texture using established pedo-transfer functions [38,91]. The root zone soil moisture (θ rz ) is estimated using a methodology presented by Bastiaanssen et al. [38], in which subsurface saturation is derived through an empirical relationship between Leaf Area Index (LAI) and remotely-sensed surface soil moisture. In the original formulation proposed by Bastiaanssen et al. [38], LAI is calculated using the Normalized Difference Vegetation Index (NDVI) and remotely-sensed surface moisture is obtained through estimates made by the AMSR-E satellite. In the current study, we substitute MODIS derived EVI for LAI as EVI has been found to be more responsive to canopy structural variations and is optimized to improve vegetation signal and reduce soil background influence [67]. Moreover, EVI was found to better correlate with observed sub-surface soil moisture estimates taken at root zone over the current study region. We also substitute derived 1 km effective saturation (θ sfc,eff ) for AMSR-E soil moisture estimates given the improvement SMOS estimates have shown over AMSR derived estimates in the region [63]. Root zone soil moisture is then estimated as: where θ rz represents the root zone soil moisture, EVI* is a normalized EVI (Equation (14)), and θ sfc,eff is the 1 km effective saturation derived using SMOS soil moisture estimates as described in Section 4.2.

Annual Total Actual Evapotranspiration Estimates
Spatial distribution of annual MOD-SMET ET for the study region is estimated between 300 mm/year and 1500 mm/year, depending on local physical attributes such as vegetation and soil physical parameters ( Figure 2). Computation of annual ET from MOD-SMET involves linear interpolation, where ET values for unavailable dates (due to error in the remote sensor or cloud-contaminated images) are linearly interpolated between two image acquisition dates. This method is generally suitable when remotely sensed images are available at regular intervals and each image captures the overall pattern of variation in ET [94]. For more information regarding interpolation between clear sky days, please refer to Singh et al. [94]. Time between cloud/error free images is six days on average (standard deviation of five days), with the largest gaps in data occurring during the monsoon season where cloudy days increase significantly. Daily ET values are then summed to an annual total ET estimate at 1 km 2 spatial resolution ( Figure 2). each image captures the overall pattern of variation in ET [94]. For more information regarding interpolation between clear sky days, please refer to Singh et al. [94]. Time between cloud/error free images is six days on average (standard deviation of five days), with the largest gaps in data occurring during the monsoon season where cloudy days increase significantly. Daily ET values are then summed to an annual total ET estimate at 1 km 2 spatial resolution ( Figure 2). MOD-SMET shows varying ET patterns within the region of interest ( Figure 2). As expected, regions of bare soil or minimal vegetation have the lowest mean annual total ET (<400 mm/year), with regions defined by more expansive vegetation and increased soil moisture availability reporting higher mean annual total ET (>1100 mm/year) ( Figure 2; Table 2). Large standard deviations (wide spatial variations) are also reported for each land cover class (Table 2), with values of 248 mm for evergreen forest, 224 mm for shrub/scrubland, and 210 mm for grassland. These three MOD-SMET shows varying ET patterns within the region of interest ( Figure 2). As expected, regions of bare soil or minimal vegetation have the lowest mean annual total ET (<400 mm/year), with regions defined by more expansive vegetation and increased soil moisture availability reporting higher mean annual total ET (>1100 mm/year) ( Figure 2; Table 2). Large standard deviations (wide spatial variations) are also reported for each land cover class (Table 2), with values of 248 mm for evergreen forest, 224 mm for shrub/scrubland, and 210 mm for grassland. These three specific land use/land cover types were chosen as they make up 97% of the total land use/land cover of the study area. As noted in Singh et al. [94], this large variability in mean annual ET for each specific land use/land cover classification may be attributed to the thematic accuracy of the NLCD map, where land cover classes may be misclassified, resulting in a wide range of annual ET. An additional source of spatial variability in ET between land cover classes may be due to the convective and episodic nature of summer precipitation in the region and the subsequent role moisture pulses have in the ET flux over the region [50,95,96]. Regions such as southern Arizona exhibit significant spatial trends and variability in precipitation at the annual scale [97], thus contributing to large variability in ET. A histogram of annual ET for the three specified land use/land cover classes shows large spatial variability (Figure 3). Roughly a quarter of the evergreen forests report an annual ET between 900 and 1000 mm/year, with another quarter attributing to estimates greater than 1200 mm/year (Figure 3). Shrub/scrubland, which accounts for an overwhelming majority of the region (84%), shows the largest spatial variability, with roughly 20% of pixels having annual total ET rates of either 400-500 mm/year or 700-900 mm/year. About half of the grassland (53%) has an annual ET between 400 and 600 mm/year, with 30% of pixels reporting estimates of 800-900 mm/year. variability (Figure 3). Roughly a quarter of the evergreen forests report an annual ET between 900 and 1000 mm/year, with another quarter attributing to estimates greater than 1200 mm/year ( Figure  3). Shrub/scrubland, which accounts for an overwhelming majority of the region (84%), shows the largest spatial variability, with roughly 20% of pixels having annual total ET rates of either 400-500 mm/year or 700-900 mm/year. About half of the grassland (53%) has an annual ET between 400 and 600 mm/year, with 30% of pixels reporting estimates of 800-900 mm/year. Values reported here are larger than those reported by Singh et al. [94], where a similar analysis was performed for the Colorado River Basin using the operational Simplified Surface Energy Balance (SSEBop) model [14]. Singh et al. [94] used a larger study domain, including regions in the more northern latitudes. These northern regions contain portions of Colorado and the Colorado Rockies. In comparison with southern Arizona, these regions experience cooler temperatures, a winter (non-growing) season, and receive reduced radiation. These differences may explain the lower annual total ET estimates from Singh et al. [94]. It is important to note, however, that the region of interest in the current study is characterized as a water limited system. Therefore, Values reported here are larger than those reported by Singh et al. [94], where a similar analysis was performed for the Colorado River Basin using the operational Simplified Surface Energy Balance (SSEBop) model [14]. Singh et al. [94] used a larger study domain, including regions in the more northern latitudes. These northern regions contain portions of Colorado and the Colorado Rockies. In comparison with southern Arizona, these regions experience cooler temperatures, a winter (non-growing) season, and receive reduced radiation. These differences may explain the lower annual total ET estimates from Singh et al. [94]. It is important to note, however, that the region of interest in the current study is characterized as a water limited system. Therefore, amplified radiation or surface temperatures over the study domain in comparison to more northern latitude regions do not directly imply an increase in ET.
High annual total ET over a majority of the study region are attributed to the influence of land surface temperature and vegetation greenness in the MOD-SMET derivation. These remotely sensed parameters remain elevated (especially MODIS-derived EVI due to its eight-day revisit time) during portions of the year where observed ET rates are much more variable, likely propagating into larger than expected annual totals. For example, an increase in duration between clear-sky days during the monsoonal period (increase in cloud cover) is likely to cause elevated estimates of ET for interpolated days that may not be representative. High evaporative demand and precipitation characterized by local, short-duration convective thunderstorms can cause rapid changes in soil moisture status and ET rates between satellite overpass times. Due to the increase in duration between clear-sky days during this time period, MOD-SMET is unable to capture the variability of ET occurring at the surface, including precipitation free days that likely exhibit drier conditions and subsequently lower ET rates. Interpolation related issues will be further discussed when comparing MOD-SMET to flux tower estimates at the point scale.

Soil Moisture Analysis
Analysis is also undertaken to better understand the role of subsurface soil moisture in MOD-SMET estimates. Observed subsurface estimates are reported at 10 cm for SRM and 15 cm for sites Whs, Wkg, and CMW. Most soil moisture dynamics and roughly 70% of the total root mass in the soil in the region is confined to depths shallower than 15 cm [50,98]. As such, observed subsurface depths between 10 and 15 cm are defined as representative of root zone soil moisture in the current study. One exception to this, however, is the trees at CMW that have been shown to access substantial quantities of groundwater at 11 m depth [53].
Examination of the influence of soil water availability on ET indicates a positive correlation between observed root zone soil moisture and observed ET at sites Whs, SRM, and Wkg ( Figure 4, top row). However, the riparian site CMW indicates an insignificant negative correlation (−0.08) when comparing in situ subsurface soil moisture to observed ET. This is likely attributed to the vegetation at the site and its ability to access groundwater at greater depths, as previously mentioned. Comparison between modeled subsurface soil moisture and its influence on MOD-SMET estimates indicate similar trends to observed, with positive correlations present at all sites, including CMW (Figure 4, bottom row). The positive correlation at CMW between modeled subsurface soil moisture and MOD-SMET is likely due to pixel heterogeneity and inclusion of non-riparian surface estimates within the MOD-SMET pixel. This concept of pixel inclusion at CMW will be further elaborated on in subsequent sections.   Lower correlation over regions of denser or healthier vegetation may explain the overestimation of ET in these regions. MOD-SMET is predicated on soil water availability being the limiting factor in the ET estimation process (see Section 4.3). As such, the soil moisture function places a large amount of weight on soil water availability in scaling PET (Equation (16)). However, in higher elevation evergreen forest where we are seeing lower correlation, soil water availability is no longer the limiting factor. Modeled subsurface soil moisture estimates in these regions are consistently larger than modeled subsurface soil moisture estimates made in neighboring lower A spatial analysis of the correlation between derived subsurface soil moisture and MOD-SMET is also performed to understand where subsurface soil moisture has the most influence on derived ET estimates ( Figure 5). Higher correlations are found in regions of minimal or less healthy vegetation (lower values of vegetation indices) where soil water availability is likely more limited ( Figure 5). Lower correlations are present over regions with denser or healthier vegetation (higher values of vegetation indices) and where soil water availability is likely less limited ( Figure 5). These results are consistent with prior studies, where positive correlations between soil moisture and ET are found when soil water availability is insufficient [99][100][101].  Lower correlation over regions of denser or healthier vegetation may explain the overestimation of ET in these regions. MOD-SMET is predicated on soil water availability being the limiting factor in the ET estimation process (see Section 4.3). As such, the soil moisture function places a large amount of weight on soil water availability in scaling PET (Equation (16)). However, in higher elevation evergreen forest where we are seeing lower correlation, soil water availability is no longer the limiting factor. Modeled subsurface soil moisture estimates in these regions are consistently larger than modeled subsurface soil moisture estimates made in neighboring lower elevation shrub/scrubland pixels (not shown). Therefore, given the formulation of MOD-SMET, ET is occurring closer to its potential rate in these higher elevation, greener regions. It is also important to note the large variability in correlation values for land cover defined as shrub/scrubland and grassland. Specifically, areas of shrub land in the east central portion (directly on and east of the riparian corridor) or southwestern portion of the study domain report lower correlations than areas of shrub land in the central and west central portions of the region ( Figure 5). Higher correlations reported over the central region of shrub/scrub land are associated with more reasonable estimates Lower correlation over regions of denser or healthier vegetation may explain the overestimation of ET in these regions. MOD-SMET is predicated on soil water availability being the limiting factor in the ET estimation process (see Section 4.3). As such, the soil moisture function places a large amount of weight on soil water availability in scaling PET (Equation (16)). However, in higher elevation evergreen forest where we are seeing lower correlation, soil water availability is no longer the limiting factor. Modeled subsurface soil moisture estimates in these regions are consistently larger than modeled subsurface soil moisture estimates made in neighboring lower elevation shrub/scrubland pixels (not shown). Therefore, given the formulation of MOD-SMET, ET is occurring closer to its potential rate in these higher elevation, greener regions. It is also important to note the large variability in correlation values for land cover defined as shrub/scrubland and grassland. Specifically, areas of shrub land in the east central portion (directly on and east of the riparian corridor) or southwestern portion of the study domain report lower correlations than areas of shrub land in the central and west central portions of the region ( Figure 5). Higher correlations reported over the central region of shrub/scrub land are associated with more reasonable estimates of annual total ET for this type of land cover (~400 mm/year) (Figure 2). Areas of decreased correlation show elevated estimates of annual total ET (~700-900 mm/year) (Figure 2). Reported discrepancies or gaps in annual total ET for the same land cover (shrub/scrub or grassland) is evident in the histogram provided in Figure 3 and may ultimately lead to skewed high estimates of annual total ET for both land cover types. Differences in correlation for defined land cover types are attributed to variations in soil physical parameters derived through STATSGO soil texture data. Estimated soil moisture stress functions used to scale PET estimates are consistently lower in regions associated with lower correlations, generating larger estimates of ET in these regions.

Ground-Based Observations
Seasonal evolution of ET shows pre-monsoon conditions in May and June are characterized by lower ET (Figure 6) with the exception of CMW where the vegetation has access to groundwater. The onset of the monsoon season near the beginning of July causes greening and an observed increase in ET. Evapotranspiration remains elevated during the monsoon season, with peak values occurring around August before steadily declining following the monsoon. Although decreasing, post-monsoon ET rates remain larger at CMW ( Figure 6D) where the deep roots and near surface water table allow transpiration to continue at higher rates after the monsoon season compared to upland sites Whs, SRM, and Wkg ( Figure 6A-C, respectively).
The onset of the monsoon season near the beginning of July causes greening and an observed increase in ET. Evapotranspiration remains elevated during the monsoon season, with peak values occurring around August before steadily declining following the monsoon. Although decreasing, post-monsoon ET rates remain larger at CMW ( Figure 6D) where the deep roots and near surface water table allow transpiration to continue at higher rates after the monsoon season compared to upland sites Whs, SRM, and Wkg ( Figure 6A-C, respectively).  Overall, modeled ET estimates compare well against observed values, with correlations ranging from 0.63 to 0.83 and RMSE errors between 40 and 96 W/m 2 (Table 3). Sites Whs and Wkg show slight positive bias over the entire time period (10 W/m 2 and 0.6 W/m 2 , respectively) while SRM and CMW report considerable negative biases (−20 and −56 W/m 2 , respectively) over the same period. MOD-SMET estimates compare much better to observed values during the drier portions of the year (non-monsoon). RMSE errors decrease to between 8 and 26 W/m 2 , with the considerable negative biases reported at SRM and CMW decreasing to −1.5 and −9.0 W/m 2 , respectively, during the non-monsoon time period. Extremely low observed and modeled ET estimates at upland sites Whs, SRM, and Wkg during the driest portion of the year (May and June) are likely due to the limited vegetation cover and shallow rooted grasses and shrubs that dominate the areas surrounding each site, subsequently resulting in extremely low observed and modeled ET estimates.
Both net radiation and air temperature are key variables in the ET (PET) scheme. We found that each are well correlated to in situ observations, with correlation coefficients ranging from 0.67 to 0.79 for net radiation and 0.92-0.93 for air temperature. Positive bias in ET reported for Whs and Wkg can be attributed to the positive bias in net radiation reported at each site (13 W/m 2 and 15 W/m 2 , respectively). Similarly, negative bias in ET reported for SRM and CMW coincides with negative bias in derived net radiation (−9.0 W/m 2 and −46 W/m 2 ). RMSE values range between 48 W/m 2 and 70 W/m 2 between sites, with the largest occurring at CMW. Due to sparse spatial information surrounding each individual flux tower site, it is difficult to explicitly identify sources of error and uncertainties. However, a clear source of uncertainty arises from the scale differences between remotely-sensed satellite estimates and ground-based in situ measurements. In addition, topographical effects and uncertainties stemming from instrumental measurements or satellite retrieval are also a factor [102]. Relative to results from Kim and Hogue [4], which report RMSE and R values of 93 W/m 2 and 0.67, on average, when comparing similarly derived net radiation estimates at the same study sites, results of the current study show relatively lower error and stronger correlation. Table 3. Correlation coefficient (R), root mean square error (RMSE) (W/m 2 ), bias (W/m 2 ), slope (-) and intercept (W/m 2 ) between observed ground-based ET estimates and MODIS ET, the eight-day empirical model, Noah, and Mosaic models. For perfect agreement, slope = 1, intercept = 0. As previously stated, riparian site CMW shows a negative bias (−56.0 W/m 2 ) over the entire study period. However, it can be seen that bias does not remain consistent, with lower bias during the beginning and end of the year and a higher negative bias found prior to, during, and shortly after the monsoon season ( Figure 6D). Negative bias reported at CMW is attributed to an observed underestimation of net radiation (−46 W/m 2 ) from pixel heterogeneity. Site CMW features predominant green leaves and dark woody stems with ample branching, leading to low albedo and subsequently high observed net radiation values. However, the limited extent of the mesquite woodland is considerably smaller than the 1 km 2 MODIS pixel size, allowing inclusion of higher albedo land surfaces to skew net radiation towards lower values that are not representative of the riparian region. Misrepresentation of net radiation during this time period leads to the underestimation of ET reported at CMW ( Figure 6D). In general, ET rates at CMW are consistently higher than those reported at the other sites, reflecting the influence of a more consistent source of soil water. The regression line between MOD-SMET and observed ET shows negative intercepts at upland sites Whs, SRM, and Wkg (positive intercept reported at riparian site CMW) with slopes ranging between 0.46 and 1.34 (the smallest of the two, 0.46, occurring at CMW) ( Table 3).

Eight-Day Empirical and NLDAS-2 Model Comparison
MOD-SMET and eight-day average E8-ET are compared against observed daily and eight-day average observed ET, respectively (Figure 7). E8-ET shows stronger correlations between sites (0.87 on average), while also reporting improved RMSE errors (16.6 W/m 2 on average) ( Table 3). MOD-SMET bias estimates are larger than those reported by E8-ET. Lower bias values reported by E8-ET may be due to the use of precipitation within the empirical models specified for upland sites. ET is strongly coupled to moisture pulses at these sites [50,95,96]. However, rapid changes in soil moisture due to frequent dry antecedent conditions and/or high evaporative demand, which is characteristic of the region, can cause rapid changes in soil moisture status between satellite overpass times. On average, slopes are better for MOD-SMET (0.93 on average) when compared to E8-ET (0.72 on average). Intercept values are positive for E8-ET at upland sites Whs, Wkg, and SRM, while being slightly negative at riparian site CMW. The opposite trend is true for MOD-SMET, which reports negative intercept values at upland sites and a positive slope at riparian site CMW. Despite larger RMSE errors and slightly lower correlations, MOD-SMET provides valuable spatially distributed ET information on an improved temporal scale (available daily when no clouds are present). It is also important to note that while MOD-SMET performs slightly worse than E8-ET, it does not require site specific calibration and can therefore be more confidently applied for spatial analysis compared to E8-ET.
E8-ET (0.72 on average). Intercept values are positive for E8-ET at upland sites Whs, Wkg, and SRM, while being slightly negative at riparian site CMW. The opposite trend is true for MOD-SMET, which reports negative intercept values at upland sites and a positive slope at riparian site CMW. Despite larger RMSE errors and slightly lower correlations, MOD-SMET provides valuable spatially distributed ET information on an improved temporal scale (available daily when no clouds are present). It is also important to note that while MOD-SMET performs slightly worse than E8-ET, it does not require site specific calibration and can therefore be more confidently applied for spatial analysis compared to E8-ET. Mosaic consistently simulates higher ET values, followed by Noah modeled ET. This is likely attributed to a shallower root zone soil moisture depth defined by Mosaic, which better matches the true root zone soil moisture in the region. These findings are consistent with those reported by Xia et al. [75], who performed an analysis comparing ET climatologies for both LSMs over six different regions in the United States.  (Table 3). Noah shows the best correlations (0.83 on average) while reporting a similar RMSE to Mosaic (32.8 W/m 2 on average for Noah and 33.3 W/m 2 for Mosaic). Slopes for both are much lower than those reported by MOD-SMET (0.40 and 0.48 for Noah and Mosaic on average, respectively), with small positive intercepts indicative of their negative bias and overall underestimation of ET. Mosaic consistently simulates higher ET values, followed by Noah modeled ET. This is likely attributed to a shallower root zone soil moisture depth defined by Mosaic, which better matches the true root zone soil moisture in the region. These findings are consistent with those reported by Xia et al. [75], who performed an analysis comparing ET climatologies for both LSMs over six different regions in the United States.
A comparison between annual total ET for each study site indicates MOD-SMET best represents annual ET at sites Whs (2% difference) and CMW (17% difference) ( Table 4). MOD-SMET largely overestimates annual total ET at Whs, where all three additional models perform relatively well (19% difference on average between Noah, Mosaic, and E8-ET modeled estimates). Both Noah and Mosaic severely underestimate annual totals at riparian site CMW (Table 4). By contrast, E8-ET shows a slightly closer estimate, which is likely attributed to the site specific calibration performed in E8-ET. MOD-SMET estimates show the lowest percent difference to observed annual total ET when averaged over all sites (37% difference on average between sites), followed by E8-ET (49% difference on average between sites). However, the distinction between MOD-SMET estimates and the additional three models is the positive bias in MOD-SMET (as noted in Section 5.2). This is most evident at sites Whs and Wkg, where estimated annual total ET is roughly double that of observations. Conversely, SRM does not show a large positive bias in MOD-SMET annual total ET. Discrepancies between Whs/Wkg and SRM originate from differences in both MODIS-derived albedo values and estimated soil physical properties. MODIS-derived albedo estimates at and surrounding sites Whs and Wkg are, on average, 20% lower than values reported for SRM. Lower albedo estimates contribute to an increase in net radiation, and by association, PET and subsequently ET. Moreover, estimated soil properties derived through STATSGO soil texture data differ between SRM and Whs/Wkg. Specifically, the derived soil moisture stress function used to scale PET reports consistently larger values at SRM when compared to the same soil moisture stress function derived for Whs and Wkg, generating larger estimates of ET at sites Whs and Wkg. A comparison between annual total ET for each study site indicates MOD-SMET best represents annual ET at sites Whs (2% difference) and CMW (17% difference) ( Table 4). MOD-SMET largely overestimates annual total ET at Whs, where all three additional models perform relatively well (19% difference on average between Noah, Mosaic, and E8-ET modeled estimates). Both Noah and Mosaic severely underestimate annual totals at riparian site CMW (Table 4). By contrast, E8-ET shows a slightly closer estimate, which is likely attributed to the site specific calibration performed in E8-ET. MOD-SMET estimates show the lowest percent difference to observed annual total ET when averaged over all sites (37% difference on average between sites), followed by E8-ET (49% difference on average between sites). However, the distinction between MOD-SMET estimates and the additional three models is the positive bias in MOD-SMET (as noted in Section 5.2). This is most evident at sites Whs and Wkg, where estimated annual total ET is roughly double that of observations. Conversely, SRM does not show a large positive bias in MOD-SMET annual total ET. Discrepancies between Whs/Wkg and SRM originate from differences in both MODIS-derived albedo values and estimated soil physical properties. MODIS-derived albedo estimates at and surrounding sites Whs and Wkg are, on average, 20% lower than values reported for SRM. Lower albedo estimates contribute to an increase in net radiation, and by association, PET and subsequently ET. Moreover, estimated soil properties derived through STATSGO soil texture data differ between SRM and Whs/Wkg. Specifically, the derived soil moisture stress function used to scale PET reports consistently larger values at SRM when compared to the same soil moisture stress function derived for Whs and Wkg, generating larger estimates of ET at sites Whs and Wkg.
Again, we note the impact of increases in duration between clear-sky days during the monsoonal period. Despite relative agreement between observed ET and MOD-SMET ET estimates during the drier portions of the year (Figure 6), MOD-SMET tends to overestimate during the monsoon period (for available days) ( Figure 6). Larger data gaps in cloud-free images during the monsoonal period is likely to cause elevated estimates for interpolated days, especially at sites Whs and Wkg. Given the shallowness of infiltration, high evaporative demand, and precipitation that is  Again, we note the impact of increases in duration between clear-sky days during the monsoonal period. Despite relative agreement between observed ET and MOD-SMET ET estimates during the drier portions of the year (Figure 6), MOD-SMET tends to overestimate during the monsoon period (for available days) ( Figure 6). Larger data gaps in cloud-free images during the monsoonal period is likely to cause elevated estimates for interpolated days, especially at sites Whs and Wkg. Given the shallowness of infiltration, high evaporative demand, and precipitation that is characterized by convective thunderstorms, it is likely that soil moisture and subsequently ET undergo rapid changes between satellite overpass times. Due to an increase in duration between clear-sky days during the monsoon season (as mentioned in Section 5.1), MOD-SMET is unable to capture the variability of ET occurring at the surface, including precipitation free days that likely exhibit drier conditions and subsequently lower ET rates.

Conclusions
The current study develops a method to determine actual ET (MOD-SMET) using multi-platform remote sensing products over a semiarid region, evaluating the product in southeastern Arizona for 2013. Large spatial variations in total annual ET were observed, corresponding to vegetative characteristics in the region. Despite well correlated spatial representation, MOD-SMET overestimates annual total ET when compared to both observed estimates and previously published studies. This may be attributed to both discrepancies in spatial soil physical properties for defined land cover classification and an increase in duration between clear-sky days during the monsoonal period, which causes elevated estimates of daily ET for interpolated days during this time period. As such, optimal application of MOD-SMET occurs when the time between cloud-free images is minimal. Spatial analysis between modeled subsurface soil moisture and MOD-SMET ET estimates indicate increased correlations over regions of insufficient soil water availability. In contrast, lower correlations are present over regions of sufficient soil water availability. Soil water availability is considered the limiting factor in the MOD-SMET estimation process. However, in regions with higher soil moisture relative to lower elevation shrubland, such as higher elevation forests, using soil water availability as the limiting factor is less valid. Therefore, MOD-SMET assumes unrealistic ET rates closer to the potential rate, leading to an overestimation in these regions. Regions defined as shrub/scrub land report large variability in correlation between modeled subsurface soil moisture and MOD-SMET, with areas reporting higher correlations coinciding with more reasonable estimates of annual total ET and those areas reporting lower correlations coinciding with elevated estimates.
Point scale comparison indicates derived MOD-SMET values are well correlated with observed values, while showing slight positive bias at Whs and Wkg and more pronounced negative bias at SRM and CMW. Discrepancies between the sites are likely attributed to site specific vegetation and soil properties, access to groundwater, and inclusion of pixel heterogeneity. This is most evident at CMW, where the extent of the riparian vegetation is considerably smaller than the 1 km 2 MODIS pixel. MOD-SMET estimates compare slightly worse than those made by the eight-day empirical ET model developed for this region. However, by comparison, MOD-SMET does not require site-specific calibration and can be more confidently applied when requiring a spatial interpretation of ET rates at the surface. NLDAS-2 LSM simulated ET shows negative bias at all sites, with Mosaic estimates slightly less bias than those reported by Noah. This is likely attributed to a shallower root zone soil moisture depth parameterized within Mosaic, which better matches root zone soil moisture in the region. Similar to spatial annual total ET rates, MOD-SMET overestimates annual total ET at each of the four sites. In contrast, Noah, Mosaic, and E8-ET all underestimate annual totals. MOD-SMET reports the lowest percent difference between modeled total ET and observed on average between the sites.
The developed approach shows that a simple ET model based on PET with a downscaled soil moisture product is an effective alternative to more complex surface-atmosphere models for estimating actual ET. Despite the limitation of only being available under clear sky conditions, the method described here can still prove beneficial in an operational setting. Specifically, the method has the ability to aid regional drought monitoring, improving water allocation and decision-making in a region where water conservation is crucial. The proposed methodology also requires limited ancillary ground-based data (only soil texture data from STATSGO is used), site specific calibration, or subjective specifications, allowing it to be transferable to ungauged basins located in water-limited regions. As such, the proposed method can help guide model validation and assimilation in data sparse regions, allowing improved predictions of land-atmosphere interactions in semiarid regions.