The development and evaluation of a long-term high-resolution dataset of potential and actual evapotranspiration for Mexico based on remote sensing data are described. Evapotranspiration is calculated using a modified version of the Penman–Monteith algorithm, with input radiation and meteorological data from the International Satellite Cloud Climatology Project (ISCCP) and vegetation distribution derived from Advanced Very High Resolution Radiometer (AVHRR) products. The ISCCP data are downscaled to ⅛° resolution using statistical relationships with data from the North American Regional Reanalysis (NARR). The final product is available at ⅛°, daily, for 1984–2006 for all Mexico. Comparisons are made with the NARR offline land surface model and measurements from approximately 1800 pan stations. The remote sensing estimate follows well the seasonal cycle and spatial pattern of the comparison datasets, with a peak in late summer at the height of the North American monsoon and highest values in low-lying and coastal regions. The spatial average over Mexico is biased low by about 0.3 mm day−1, with a monthly rmse of about 0.5 mm day−1. The underestimation may be related to the lack of a model for canopy evaporation, which is estimated to be up to 30% of total evapotranspiration. Uncertainties in both the remote sensing–based estimates (because of input data uncertainties) and the true value of evapotranspiration (represented by the spread in the comparison datasets) are up to 0.5 and 1.2 mm day−1, respectively. This study is a first step in quantifying the long-term variation in global land evapotranspiration from remote sensing data.
Evapotranspiration (ET) is a key variable in the terrestrial water, energy, and carbon cycles and provides the link between the land surface and the atmospheric boundary layer (Guo et al. 2006; Betts 2007; Lawrence et al. 2007). It is also potentially the largest component of the land water budget after precipitation, especially in semiarid regions, and is therefore central to quantifying water cycle variability and managing water resources. However, at the spatial scales necessary for making informed decisions about, for example, water allocation, reservoir management, and irrigation scheduling (Bastiaanssen et al. 2005), ground observations of ET are essentially nonexistent. At regional-to-continental scales, this lack of data also hampers the evaluation of large-scale variability and change in the terrestrial water cycle, and hinders improvement of land parameterizations in climate forecast and projection models. An emerging potential source of data is from remote sensing that has the advantage of large-scale and near-temporally continuous coverage (Entekhabi et al. 1999). This is finding real-world applications in the agricultural sector (Droogers and Bastiaanssen 2002; Bandara 2003; Consoli et al. 2006; Tasumi and Allen 2007; Santos et al. 2008), where high-resolution data are required to quantify water use efficiency, and for climate research, where long-term, large-scale datasets are needed (Choudhury et al. 1998; Mu et al. 2007; Fisher et al. 2008).
While evapotranspiration itself is a conceptually simple process, it has proven difficult to accurately retrieve it over natural surfaces using satellite observations. The difficulties include accurately representing the physical mechanisms controlling the process; heterogeneities at the land surface; the interacting influences and feedbacks in the soil, vegetation, and meteorological continuum; and imprecise knowledge of variables controlling evapotranspiration. Typical approaches for estimating evapotranspiration from remote sensing are based on direct empirical regressions that relate evapotranspiration to surface temperature and vegetation indices (e.g., Nemani and Running 1989; Nishida et al. 2003a); predictive equations that approximate the main evaporative processes (e.g., Monteith 1965; Priestley and Taylor 1972) or surface energy balance considerations that use radiometric surface temperatures to solve the coupled equations of sensible, latent, and available heat energy directly (Bastiaanssen et al. 1998; Su 2002; Norman et al. 1995; Anderson et al. 1997); or through assimilation (Caparrini et al. 2004). A comprehensive overview of methods for estimating surface heat fluxes using satellite radiometric data is given by Kalma et al. (2008).
Significant progress has been made in recent years, with experiments showing systematic agreement between estimates from ground and aircraft measurements over several test areas (French et al. 2005; Su et al. 2005; McCabe and Wood 2006; Kustas et al. 2005). However, much of the remote sensing results (e.g., Norman et al. 2003; Anderson et al. 2004; McCabe et al. 2008) are at small scales using high-resolution thermal imaging with low temporal repeat and a focus on agricultural areas. There are few results of regional estimates that only utilize remote sensing measurements (e.g., Nishida et al. 2003b; Houborg and Soegaard 2004; Anderson et al. 2007a,b; Pan et al. 2007; McCabe et al. 2008), and these are generally limited to subseasonal periods. More recently however, these methods have been applied to estimate evapotranspiration at larger scales, in particular because of interest in Earth Observing System (EOS)-era data to provide satellite remote sensing inputs at continental-to-global scales of the required vegetation, radiation, and near-surface meteorological data. Cleugh et al. (2007) compared the Penman–Monteith (PM) model (Monteith 1965) and an energy balance model based on surface temperature taken from the Moderate Resolution Imaging Spectroradiometer (MODIS) at two tower sites in Australia. They found the PM model to perform better mainly because of the high sensitivity of the energy balance approach to small errors in the surface temperature and no balance constraint on sensible heat. Mu et al. (2007) built upon the work of Cleugh et al. (2007) with adjustments to the stomatal conductance based on environmental controls and applied the algorithm globally. Leuning et al. (2008) replaced the earlier empirical models of surface conductance of Cleugh et al. (2007) and Mu et al. (2007) with a biophysical two-parameter model and calibrated the PM model for 15 flux towers globally. In contrast, Zhang et al. (2008) used estimates of mean annual evapotranspiration estimated from long-term observations of precipitation minus runoff for 120 catchments in southeast Australia to calibrate the PM model. Fisher et al. (2008) developed a model based on the reduction of potential evapotranspiration (PE) from a Priestley–Taylor model (Priestley and Taylor 1972) to ET using biophysical and soil moisture controls derived from remote sensing and developed global estimates of ET for 1986–93 based on Advanced Very High Resolution Radiometer (AVHRR) and International Satellite Land Surface Climatology Project (ISLSCP) initiative data.
Similarly, the authors have retrieved continental-scale land evapotranspiration using EOS data with both the PM method (Ferguson et al. 2008) and a surface energy balance method (Sheffield et al. 2008) using the Surface Energy Balance System (SEBS; Su 2002). These approaches show considerable potential for monitoring ET and surface fluxes based on current EOS data. With the development of legacy datasets, such as the ISLSCP Initiative II (ISLSCP-II) (Hall et al. 2006) and the International Satellite Cloud Climatology Project (ISCCP) (Zhang et al. 2004), the opportunity exists to extend these estimates back in time to form multidecadal, albeit low spatial resolution, but high temporal resolution (subdaily) time series. Such long-term datasets can help improve estimates of the large-scale water cycle, its variability, and feedbacks with climate change, and provide useful benchmarks for evaluating and improving land parameterizations in climate forecast and projection models.
In this paper, we test the viability of generating a long-term dataset of ET from mostly remote sensing data using the ISCCP database and of obtaining high-resolution estimates through downscaling of the input data. The downscaling provides estimates at spatial scales that reflect the heteorogeneity of the land surface (McCabe and Wood 2006), which are relevant to water resources planning (Tang et al. 2009). Our domain of interest is Mexico, which is a generally moisture-limited region but encompasses a wide range of biogeographical, topographic and climatic variation, thus providing a good test of the approach. There also exists a range of long-term datasets of in situ, reanalysis, and modeled data for the region that can be used for evaluation. We use the algorithm of Mu et al. (2007) that is based on the PM method (Monteith, 1965). The PM method is generally considered the most accurate method for estimating PE and hence ET over large vegetated regions, when all necessary input data are available (e.g., Garatuza et al. 1998). The input data and calculated PE and ET are evaluated against other large-scale estimates from reanalysis and offline hydrologic modeling. We also calculate a reference crop evapotranspiration ETrc and compare these to local estimates from a large database of evaporation pans. We use these evaluations and the uncertainties in the input data to estimate the uncertainty in the derived ET.
We use the modified Penman–Monteith algorithm of Mu et al. [2007; remote sensing PM (RS-PM)], which is designed specifically for use with remote sensing data. To produce estimates at the finescale, we used a downscaling scheme described later to bring the ISCCP input data to ⅛° resolution. Using the downscaled data, ET was calculated using the RS-PM model at daily, ⅛° resolution for 1984–2006. Here PE and ETrc are calculated using standard versions of the PM algorithm (Shuttleworth 1993) and are hereafter also referred to with the RS-PM prefix when forced by remote sensing data. Next we describe the details of the RS-PM model, the method for downscaling the input data, and the approach used for the evaluation and uncertainty analysis.
a. Remote sensing Penman–Monteith algorithm
The RS-PM formulation presented in Mu et al. (2007) is based on the PM equation (Monteith 1965) and is a revised version of the method presented by Cleugh et al. (2007). The PM equation given in Eq. (1) models the diffusion of energy from plants or soil against stomatal and aerodynamic resistance, given inputs of net radiation, temperature, and humidity:
where Δ (Pa K−1) is the slope of the vapor pressure deficit versus air temperature relationship, Rnet is the net radiation (W m−2), ρ is the density of air (kg m−3), cp is the specific heat of air (J kg−1 K−1), VPD is the vapor pressure deficit (Pa), and γ is the psychometric constant (Pa K−1). Here ra and rs are the aerodynamic and canopy (or surface) resistances (s m−1), respectively. ET (W m−2) is converted into water equivalent units of millimeters per day via the latent heat of vaporization and density of water, and these units are used hereafter. Canopy resistance is derived by upscaling stomatal resistance via leaf area index (LAI), and the RS-PM invokes additional environmental adjustments to resistance as a function of VPD and minimum temperature. As recognized by Leuning et al. (2008), stomatal resistance is a key parameter that varies by vegetation type. In contrast to Mu et al. (2007), who used a constant minimum stomatal resistance for all vegetation types, we implement vegetation-dependent values as used in the VIC model and developed by Maurer et al. (2002). Evaporation directly from the soil surface is also calculated as a function of relative humidity and VPD. Note that evaporation from canopy-intercepted water is not modeled and that the model parameters, such as minimum stomatal resistances, are not calibrated. Potential evaporation and reference crop evapotranspiration are estimated following Shuttleworth (1993):
Here u is the wind speed (m s−1) and λ is the latent heat of vaporization (J kg−1). Equation (2) represents the potential evaporation from an extensive free water surface and is equivalent to Eq. (1), but with rs = 0 s m−1 (i.e., no resistance to transpiration) and standard measurement heights for the meteorological data. Reference crop evapotranspiration [Eq. (3)] is an idealized estimate that we compare to the pan observations. It is defined as the rate of evapotranspiration from an idealized grass crop of height 0.12 m, an albedo of 0.23, and a surface resistance of 69 s m−1. It is often related to pan evaporation with a pan coefficient (ETrc/pan) that varies by the type of pan, surrounding ground cover (e.g., short grass), and mean humidity and wind speed/fetch conditions. The coefficient is generally of the order of 0.4–0.85 for U.S. class “A” pans and up to 1.1 for sunken Colorado-type pans under low wind conditions (Shuttleworth 1993).
Required inputs for calculating ET, PE, and ETrc are net radiation Rnet, surface shortwave (SWdown, SWup) and longwave (LWdown, LWup) radiation, surface meteorological data (humidity q, air temperature Ta, surface temperature Ts, pressure Ps, and wind speed u), and surface characteristics (vegetation distribution, LAI, emissivity ɛ, and albedo α). Values of emissivity and albedo are only required if upward SW and LW fluxes are not available. The radiation data, meteorological data, emissivity, and albedo were obtained from the ISCCP database at daily resolution. Wind speed is not available as part of the International Satellite Cloud Climatology Project (ISCCP) database and so was taken from the latest version of the dataset of Sheffield et al. (2006). Vegetation distribution and LAI values were taken from AVHRR products as used by Zhu and Lettenmaier (2007). Ground heat flux is estimated as a function of Rnet and depends on the fractional vegetation coverage, varying linearly between 5% (fully vegetated) and 35% (bare ground) of net radiation (Su 2002).
b. Spatial downscaling
To produce high spatial resolution estimates, we developed a simple downscaling scheme to bring the ISCCP input data to ⅛° resolution, with the North American Regional Reanalysis (NARR) providing the underlying spatial distribution. The downscaling scheme works as follows for a particular variable, such a shortwave radiation, and is repeated for each ISCCP variable:
At the monthly time scale, the NARR data (32 km) are aggregated to 2.5°.
The ratio of the NARR to the ISCCP data at 2.5° are bilinearly interpolated to ⅛° resolution and smoothed in space using a 9-point weighted filter.
These ratios are then used to scale the NARR monthly data at ⅛° to produce a high-resolution dataset that matches the ISCCP data when aggregated to 2.5°.
The ISCCP daily data are bilinearly interpolated to ⅛° and then scaled so that their monthly mean matches the scaled ⅛° monthly NARR data calculated in the previous step.
c. Evaluation and uncertainty analysis
The RS-PM estimates were evaluated by comparison with estimates taken from offline land surface modeling, reanalysis, and station measurements. The offline ET data are taken from a long-term simulation using the Variable Infiltration Capacity (VIC) hydrologic land surface model. The VIC model was calibrated to measured streamflow at a number of unmanaged basins across Mexico (Zhu and Lettenmaier 2007; Munoz-Arriola et al. 2009) and is therefore expected to give reasonable estimates of ET (at least at annual scales) as it was forced by observed precipitation and temperature. Zhu and Lettenmaier (2007) reported Nash–Sutcliffe efficiency values of greater than 0.5 for most basins and generally greater than 0.7 for more humid basins, with biases less than 25%. To test how well the VIC data replicate observation-based ET, we compared the RS-PM and VIC estimates of ET over these basins against ET calculated as a residual of the long-term land water balance (ET = P − Q), where P is precipitation and Q is streamflow, and storage change is assumed to be negligible. Observations of P and Q are taken from Zhu and Lettenmaier (2007). This comparison is a test of the RS-PM against an observation-based estimate of ET and a test of the VIC data in providing a best estimate of ET over Mexico.
The reanalysis ET data are taken from the NARR, which is a state-of-the-art high-resolution reanalysis of the atmospheric–land water and energy cycles for North America. The NARR is also used to infer ET as a residual of the atmospheric water balance (ET = P − MC − dw/dt), where MC is atmospheric total column water vapor convergence and dw/dt is the change in column water vapor. The NARR assimilates gauge precipitation and radiosonde data of atmospheric temperature and humidity and is therefore likely to do a reasonable job of replicating the variation in atmospheric water vapor and thus ET as a residual, at least over large time and space scales (Sheffield et al. 2009). We also evaluated the data against observations from a large database of station measurements, including pan evaporation. Uncertainty in the RS-PM estimates is derived from uncertainties in the forcing data (radiation, meteorological data, and surface characteristics) as well as variations in the parameterizations used in the algorithm (e.g., the calculation of canopy resistance). The effects of these are explored in section 4e. We also discuss the likely effect of evaporation of intercepted rainfall (canopy evaporation), which is not modeled in the RS-PM algorithm, on the ET estimates.
The primary data sources for the evapotranspiration estimates are the ISCCP satellite–based product that provides radiation and other meteorological variables and surface parameters, and the AVHRR vegetation products. We also describe the evaluation datasets (VIC, NARR, and the station data) and two remote sensing datasets that are used as alternative input sources to assess the uncertainties in the ET estimates [the National Aeronautics and Space Administration (NASA) Langley’s Surface Radiation Budget (SRB) dataset and MODIS-based vegetation products]. Details of the various datasets used are given and are summarized in Table 1.
a. ISCCP remotely sensed radiation, meteorological data, and surface characteristics
The ISCCP (Schiffer and Rossow 1983) was established in 1982 to collect and analyze satellite radiance measurements to infer the global distribution of clouds and their role in radiation and the hydrologic cycle (Rossow and Dueñas 2004). We use the Climatological summary product [flux data-surface (FD-SRF)] Radiative Flux (RadFlux) version of the dataset, which contains surface radiative fluxes and a summary of the physical quantities used to calculate them. Shortwave and longwave radiative flux profiles are calculated using remotely sensed and modeled datasets to specify the properties of earth’s atmosphere and surface. We used the daily averaged data and resampled to 2.5° equal-angle resolution.
b. SRB remotely sensed radiation
NASA Langley Research Center’s version 2.5 (V2.5) product is based on satellite data from the ISCCP C1 data product and the Earth Radiation Budget Experiment (the recently released V3 extends to 2007). Two versions are available for shortwave [SRB-SW, SRB-Quality Check (QC)SW] and longwave (SRB-LW, SRB-QCLW) radiation, depending on the retrieval algorithm used. Comparison of these products with surface measurements has indicated that no one product is globally superior, although regional comparisons have indicated better performance for the SRB-SW and SRB-QCLW products, as discussed in Sheffield et al. (2006). We therefore use the V2.5 SRB-SW and SRB-QCLW in this paper. The SRB data are resampled from their native nested grid to a 1.0° equal-angle grid and at daily time steps for 1984–2006.
c. AVHRR and MODIS vegetation distribution and LAI
The distribution of vegetation fractional cover (Fc) is taken from the AVHRR-based dataset of Hansen et al. (2000). Values of LAI are specified for each vegetation type that exists in each grid cell by resampling the dataset of Myneni et al. (1997), which is based on AVHRR normalized difference vegetation index (NDVI) values. The LAI values are specified for each month but do not vary from year to year. This is consistent with the climatological inputs to the VIC and NARR (in terms of greenness fraction) and because the AVHRR data do not go beyond mid-2001. We ran a RS-PM simulation using annually varying AVHRR LAI data for 1984–2000 that showed daily differences in ET averaged over Mexico of less than 0.2 mm day−1, and negligible differences when averaged over all years. Related work (Ferguson et al. 2010) has indicated that differences in vegetation datasets (including LAI) have a much larger effect on the ET estimates. To assess the effect of uncertainties in vegetation, we used land cover and LAI data developed for the North America Carbon Program (NACP) as derived from MODIS. The NACP LAI data have been gap filled and smoothed by Gao et al. (2008) to allow for cloud contamination and other factors that affect the quality of the retrievals. These data were sampled from their native 1-km resolution to ⅛° in a similar manner to the AVHRR data.
d. North American Regional Reanalysis
Mesinger et al. (2006) describe the development of the NARR, including an overview of the improvements relative to the first National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis (Kalnay et al. 1996) and initial evaluations. The largest improvement is deemed to be the assimilated observed precipitation, and the knock-on effect to land–atmosphere coupling and land surface hydrology, which also benefited from improvements to the land surface model, Noah (Mitchell et al. 2004). The precipitation is assimilated in terms of latent heating profiles that then force precipitation in the model, producing a close match to observed precipitation (Mesinger et al. 2006). The NARR surface radiation and near-surface meteorological data are used as a basis for the spatial downscaling of inputs to the RS-PM. ET taken directly from the NARR, and calculated as a residual of the NARR atmospheric water, is used to evaluate the RS-PM ET.
e. VIC hydrologic simulation
The VIC land surface model (Liang et al. 1994, 1996; Cherkauer et al. 2002) simulates the terrestrial water and energy balances and distinguishes itself from other land surface schemes through the representation of subgrid variability in soil storage capacity as a spatial probability distribution. Horizontally, VIC represents the land surface by a number of tiled land cover classes, which are specified by their grid fractional area, LAI, canopy resistance, and vertical root distribution. ET is calculated using a PM formulation with adjustments to canopy resistance to account for environmental factors, including available soil moisture. Canopy evaporation is also calculated. We ran the VIC model using the data of Zhu and Lettenmaier (2007), who previously ran the VIC model to assess the long-term (1925–2004) hydrologic cycle over Mexico. The simulation was forced by their gauge-based dataset of precipitation and temperature, with wind speed taken from the NCEP–NCAR reanalysis. Other required meteorological data (humidity, pressure, radiation) were calculated using predictive relationships (Thornton and Running 1999) based on the precipitation and temperature data.
f. Mexican National Meteorological Service station data
Station data were obtained from the Servicio Meteorológico Nacional (SMN) of the Mexican Ministry of Environment and Natural Resources [Secretaría de Medio Ambiente y Recursos Naturales (SEMARNAT)]. The data consist of daily precipitation, maximum and minimum temperature, and pan evaporation for an average of 1834 stations with periods ranging from 1961–2007. Class A pans were initially set up in 1922, reaching a relatively constant and widespread distribution after the 1970s. Between 1984 and 2006, the number of operating pans decreased from 2931 in 1984 to 1181 in 2006.
4. Results and discussion
We first compare the main inputs to the RS-PM (ISCCP radiation and temperature) with other available datasets (SRB, NARR, VIC) at the ISCCP resolution of 2.5°. This is followed by an evaluation of the downscaling method in terms of how well it replicates the spatial distribution of the NARR data and how well the NARR represents the spatial distribution of measured station data. We then show the high-resolution PE, ETrc, and ET fields calculated using the RS-PM and evaluate these against data from the small basin water balance estimates, pan measurements, VIC, and NARR. Lastly, we discuss the effect of the uncertainties and simplifications of the method and the choice of input data.
a. Comparison of input data
The main forcing of ET is net radiation, as derived from the residual of the balance of downward and upward shortwave (SWup, SWdown) and longwave (LWup, LWdown) radiation fluxes. We compared the ISCCP, SRB, and NARR radiation data to assess uncertainties in the ET estimates through uncertainties in radiation. The ISCCP and SRB fluxes share some common inputs but use different retrieval algorithms. Zhang et al. (2004) showed that the errors in ISCCP monthly-mean surface fluxes relative to station measurements from the Global Energy Balance Archive database (Ohmura and Gilgen 1991) and the Baseline Surface Radiation Network (Ohmura et al. 1998) are of the order of 2 to 9 and −15 to 2 W m−2 for SW and LW, respectively, which are within estimated measurement accuracy. For northern Eurasia, Troy and Wood (2009) found the ISCCP to compare well with station measurements from the World Radiation Data Centre, with errors within the ranges reported by Zhang et al. (2004). The bias for stations within 15°–35°N (latitude band of Mexico) was reported to be higher by Zhang et al. (2004), at 15 and 5 W m−2 for SW and LW, respectively. Nevertheless, these results show that the ISCCP data are a reasonable estimate of large-scale surface radiative fluxes. By inference, validation of shortwave fluxes indicates that the ISCCP albedo is also reasonable (Wang et al. 2006).
Comparisons of the ISCCP, SRB, and NARR surface radiation components, albedo, and surface temperature data over Mexico are shown in Figs. 1 and 2 and summarized in Table 2. The ISCCP SWdown is biased high compared with the SRB by about 10 W m−2 on average, with the bias largely concentrated in the summer and reducing to near zero in the winter. The NARR is biased high by more than 30 W m−2 on average, and this is fairly uniform over the year. Interestingly, the bias in the NARR SWup is comparable to that for SWdown because of the high bias in its albedo [mean high bias of 0.09 (90%) and 0.12 (58%) compared with ISCCP and SRB, respectively], and thus the effect of the higher SWdown on Rnet is negligible. Differences in the LWdown flux are small between the three datasets (∼5 and 10 W m−2 high bias relative to SRB and NARR, respectively). For LWup, the ISCCP is generally lower than the other datasets by less than about 15 W m−2 over the year, although during the summer it is biased low by more than 30 W m−2, which is related to the low bias in surface temperature relative to the NARR (−4.4 K for June–August; surface temperature is not available from the SRB dataset). In terms of Rnet, the differences in the individual components translate into the ISCCP being biased high by about 30 W m−2. In the summer, this bias is about 48–50 W m−2. This high bias is mainly due to the low bias in the ISCCP LWup (annual bias of −10 to −13 W m−2; summer bias of −21 to −27 W m−2). An example of the spatial distribution of the differences between the ISCCP and NARR (the differences between ISCCP and SRB are similar) is shown in Fig. 2 for Rnet, which also illustrates the coarseness of the ISCCP 2.5° grid over Mexico and the need for downscaling.
We also compared the near-surface air temperature (Ta) fields between the ISCCP, NARR, and VIC datasets (Fig. 3; Table 2). Here Ta data are not included in the SRB database. The VIC data are based on spatially interpolated station observations and so are likely the best estimates of the regional distribution. Over all Mexico, the mean monthly ISCCP data are biased low by 0.3 K (rmse = 1.47 K), with peak gridcell biases up to 19 K during the summer. The bias and rmse relative to the NARR are 0.2 and 1.5 K, respectively.
b. Evaluation of downscaled input data
Figure 4 shows examples of the downscaled radiation and meteorological fields for July and compares these to the NARR, which provided the finescale spatial distribution for the downscaling. In general, the downscaling method captures the spatial variability given by the NARR, especially for LWdown, Ta, q, Ps, and u. For other variables, the differences are overwhelmed by the biases between the original ISCCP data and the NARR, such as for SWdown, α, Ts, and thus Rnet. For Ps and q, there is a clear delineation of the sign of the differences between high elevations and coastal regions, implying that some further correction may be necessary, perhaps based on intervariable consistency and elevation effects (Sheffield et al. 2006).
The accuracy of the downscaling method relies on how well the method can reproduce the spatial distribution of the NARR data and the accuracy of the spatial distribution of the NARR data itself. As well as visual inspection of the downscaled data, we also evaluated the method by comparing three sets of data. First, we compared the NARR with a downscaled version of itself to validate the method. In this case, the NARR data was aggregated to the 2.5° ISCCP grid and then downscaled to ⅛°. Second, we compared the downscaled ISCCP data with that obtained from simple interpolation methods (bilinear interpolation and box averaging). Third, we compared the NARR data and the downscaled ISCCP data with gridded station data analyses. This is a direct comparison of the NARR data with observational data and is thus an evaluation of whether the spatial distribution of the NARR data is reasonable and can serve as the basis for the downscaling.
Figure 5 shows the mean seasonal cycle of the spatial correlation between the NARR data (interpolated from 32 km to ⅛°) for the four downscaled datasets (the downscaled NARR data, the downscaled ISCCP data, the interpolated ISCCP data, and the gridbox-averaged ISCCP data). Comparison of the downscaled NARR data with the original NARR data (solid line) shows a near-perfect match as expected. The gridbox-averaging method (dotted line) generally gives the lowest correlation, closely followed by the interpolation method (dotted–dashed line), indicating the high spatial variability in the NARR data that cannot be captured by simple averaging or interpolation. The downscaled ISCCP data (dashed line) show a much higher correlation that is rarely less than 0.7. The correlations for albedo are generally the least (annual average correlation = 0.53), and this is reflected in the biases in the actual values between the ISCCP and NARR (Table 2). This is expected given the high spatial variability in albedo as driven by spatial variability in land cover type and also the within-type variability, such as for grassland (Gao et al. 2005) and bare ground, as dependent on soil type (Tsvetsinskaya et al. 2006). Although the downscaling method will force the ISCCP data to match the spatial distribution of the NARR data at the monthly time scale, they will still differ because of the differences between the NARR and ISCCP at the coarse (2.5°) resolution.
We evaluated the spatial distribution of the NARR data to see how well it replicates the spatial distribution of available measurements and thus whether it is appropriate to use as a basis for the downscaling (Fig. 6). This was done against gridded station observations of air temperature used to force the VIC simulation, which is the only meteorological variable (in addition to precipitation) for which relatively dense station observations are available over all Mexico. The station data are also used as an independent evaluation of the downscaled ISCCP data. The mean seasonal cycle and spatial distribution of the NARR and VIC data are similar with a mean bias and spatial correlation of −0.3 K (NARR lower) and 0.9, respectively, indicating that the NARR data give a reasonable representation of the finescale variability in air temperature. The downscaled ISCCP data compare well with both the NARR and VIC data in the winter and autumn (r > 0.8), but they drop (to r ≈ 0.6) in the early summer and are biased low by about −2 K relative to both datasets in July. At the gridcell scale, the scatterplots (Fig. 6, bottom panel) indicate the greater spread for the VIC comparison. In the absence of station data for the other variables (radiation, pressure, humidity, among others.), we assume that the NARR is our best estimate for these variables.
c. High-resolution estimates of potential evaporation
Using the downscaled data, PE was calculated using Eq. (2) and is shown in Figs. 7 and 8 compared with the NARR (PE is not available as a standard output from the VIC model). Figure 8 also compares the RS-PM ETrc and measured pan evaporation, and this is discussed later. The seasonal cycle of PE generally follows the cycle of net radiation, which peaks in the summer, but it is tempered by high humidity (low VPD) at this time of year and thus peaks in the late spring. Peak values in the spring tend to be clustered along the coasts and then shift to the north of the country in the summer following the North American monsoon (Higgins et al. 2006). The spatial patterns and annual cycle are similar between the RS-PM and NARR, with the RS-PM lower in the spring and higher in the autumn. These biases tend to cancel each other out over the year to give an average bias of 0.44 mm day−1. This is despite the ISCCP net radiation being 28% higher on average than the NARR (Table 2).
We carried out an experiment that forced the RS-PM with the downscaled version of the NARR net radiation (all other forcings unchanged), which showed an expected reduction in PE (not shown). The difference between the two datasets is therefore the result of other factors, including differences in the other meteorological variables. A second experiment that forced the RS-PM with only NARR downscaled radiation and meteorological data (not shown) showed that even with the same forcings (albeit monthly), the RS-PM underestimates PE compared with NARR by 1.5 mm day−1 (the effect of using monthly data versus daily data was explored by comparing the RS-PM output using monthly average and daily ISCCP data, which showed very little difference). There are a number of reasons that may cause the difference in NARR PE and NARR-forced RS-PM PE, including foremost the differences in the specification of the aerodynamic resistance. In the NARR, PE is calculated based on the Mahrt and Ek (1984) modified version of the PM equation that includes the influence atmospheric stability on turbulent water vapor transfer. Mahrt and Ek (1984) found this to increase the aerodynamic term by up to 2 times for afternoons of modest instability for their test location in Australia. The differences may also be the result of ignoring the diurnal covariability of atmospheric variables that can lead to differences between subdaily integrations of the PM equation (as for the NARR) versus calculations using daily averaged values (as for the RS-PM approach used here). Mahrt and Ek (1984) estimate that for the data they looked at, this was a small effect because of effective calibration in the PM, but it may be an influence here. Furthermore, the available energy in the NARR may differ because of the online calculation of ground heat flux and account for the latent heat released from precipitation.
Figure 8 also shows the RS-PM ETrc and measured pan evaporation. The daily pan measurements are averaged to monthly values and then sampled on the same ⅛° grid as the other datasets. Where more than one station lies within a grid cell, their average value is used. The areally averaged seasonal cycle of the RS-PM ETrc follows the RS-PM PE (Fig. 8) but is uniformly lower because of the canopy resistance term used in its calculation, which models the reduction of PE to ET. The RS-PM ETrc and pan measurements are similar, with the RS-PM being slightly lower during February–April and higher from June–November. Figure 9 shows seasonal maps of the RS-PM ETrc and the pan data. There is a tendency for the RS-PM data to be higher in coastal regions, especially along the Caribbean coast, and within about 1 mm day−1 in the central highlands. We investigate the spatial relationship between the ETrc and pan data further in Fig. 10, which shows the spatial distribution of their annual values and their ratio (pan coefficient). There is a delineation in the ratio values between the drier highland/inland and the more humid coastal regions, with the latter being within the generally accepted range of pan coefficients for U.S. class A pans (0.4–0.85; Shuttleworth 1993). In the drier inland regions, the ratios can reach one or above, approaching the upper limit of the accepted range for Colorado-type pans (<1.1). Note that the downscaled ISCCP humidity is lower than the NARR data along coastal regions (see Fig. 4) and therefore would tend to overestimate PE and therefore ETrc.
d. High-resolution estimates of actual evapotranspiration
We first compared the RS-PM and VIC estimates with those from the small basin water balance as derived from measured P − Q (Fig. 11). The VIC data are generally close to the P − Q values, with a high bias of 0.1 mm day−1 (4.9%) on average across all basins (area weighted). The largest differences are in the more humid basins in the south (e.g., basins 11 and 14). For basin 11 (Poza Rica), the VIC model simulates very little baseflow as compared to measured streamflow, which accounts for much of the bias in ET. For basin 14 (Las Perlas), the VIC model tends to underestimate peak monthly streamflow in some years, which may be due to errors in the precipitation. The RS-PM estimates are biased low by about −0.25 mm day−1 (−14.7%) across all basins.
Next we compared the calculated RS-PM ET to estimates from the VIC model, the NARR, and inferred values from the NARR atmospheric water budget. Their seasonal cycles over Mexico are compared in Fig. 12. The comparisons are limited to 1984–2003 because of the period of the VIC simulation. Overall, there is general agreement in the amplitude and phase of the seasonal cycle. There is a high bias in the spring compared to the VIC and NARR datasets. During the summer, the RS-PM falls somewhere in between the two other datasets, and there is a low bias in the autumn and winter. The overall bias and rmse are given in Table 3. Also shown in Fig. 12 is the mean seasonal cycle derived from the atmospheric water balance using NARR precipitation, moisture convergence, and change in atmospheric moisture. This comparison shows good agreement during the peak months and end of year, but the inferred values are much lower during the spring. Inspection of the full time series of NARR atmospheric variables that are used to infer evapotranspiration reveals a switch to lower mean precipitation from the 1994 drought onward coupled with an increase in moisture convergence, which reduces the inferred evapotranspiration considerably. The lower precipitation is not apparent in the VIC station–based series. Replacing the NARR precipitation with the VIC precipitation increases the range of the inferred evapotranspiration (not shown).
Figure 13 gives maps of the RS-PM- and VIC-simulated ET, which shows that their spatial distributions are similar. The tendency for the RS-PM data to be lower in the summer and autumn, as seen in the time series in Fig. 12, is generally confined to the warmer and more humid coastal regions. The maps in Fig. 14 show that the differences of the RS-PM with the NARR are more randomly distributed, although the RS-PM is again generally lower in coastal regions. A summary of the error statistics for the RS-PM ET is given in Table 3. The overall bias for 1984–2006 over Mexico is about 0.1–0.2 mm day−1 (RS-PM low) with rmse values of about 0.5 and 0.3 mm day−1 at monthly and annual time scales, respectively.
e. Uncertainties in ET
Uncertainties in the RS-PM ET are attributable to uncertainties in the input forcings and vegetation parameters. It should be noted that this is different in concept to errors in the ET estimates, which are relative to the true value of ET (which itself is unknown). The errors are derived from errors in the forcings, but also from errors in the RS-PM formulation and how well it represents the true physical processes. In this section, we are concerned only with the uncertainties in the data and not with the errors relative to other estimates (e.g., VIC and NARR). We therefore take the RS-PM as our best representation of the evapotranspiration process when using remote sensing data; however, we explore the uncertainties by using different sources of inputs, specifically net radiation and LAI, to which the RS-PM is most sensitive. As we are interested in using remote sensing data, we limit this interest to using data from the SRB dataset for radiation, which is the only comparable long-term remote sensing-based dataset, and the NACP/MODIS vegetation parameters as an alternative to the AVHRR data. Using the SRB net radiation decreases ET by −0.13 mm day−1 (−9.0%) annually over Mexico, with the smallest absolute decrease in February (−0.04 mm day−1, −5.5%) and largest in August (−0.29 mm day−1, −10.0%). This is expected given the high bias in the ISCCP Rnet relative to the SRB (and NARR; section 4a). Using the NACP/MODIS vegetation parameters, ET increased by 0.27 mm day−1 (18.6%) annually and ranged from an increase of 0.52 mm day−1 (39.4%) in May to a slight decrease of −0.03 mm day−1 in August (−1.0%). These changes are mostly a reflection of the earlier and smoother seasonal cycle in the MODIS LAI compared with the AVHRR.
Figure 15 shows an estimate of the uncertainty in the RS-PM based on the uncertainties in the input Rnet and vegetation data and compares this with the uncertainty in our best estimate of large-scale ET from NARR and VIC. We assume that the uncertainties in the RS-PM estimates are normally distributed with standard deviation σ derived from the spread of the RS-PM estimates, and differ by month. Similarly, the uncertainty in the true value of ET is based on the difference between the VIC and NARR data. Monthly values of σ vary between 0.06 and 0.11 mm day−1 for the RS-PM estimates and between 0.15 and 0.27 mm day−1 for the VIC and NARR data, with the highest values in the summer. Figure 15 shows the mean of the RS-PM estimates and the 95% confidence limits calculated as ±2σ and applied as a relative error. Overall, the RS-PM and NARR–VIC estimates overlap in all months but less so in the spring and autumn, when the RS-PM is biased high and low, respectively. The uncertainty in the RS-PM values range from 0.19 to 0.49 mm day−1, whereas the uncertainty in the NARR–VIC values range from 0.31 to 1.24 mm day−1.
f. Effect of canopy evaporation
One potential drawback of the RS-PM method is that it does not calculate direct evaporation of water lying on the vegetation canopy (interception), which may explain part of the low bias in the RS-PM relative to the NARR and VIC. In general, canopy evaporation may account for 25%–75% of total evapotranspiration from forests and up to 60% of total rainfall (David et al. 2005). For the VIC simulation, canopy evaporation accounts for about 30% of total evapotranspiration on average at the annual time scale, indicating that it can be significant, although this may be somewhat overestimated because the simulation was forced by daily precipitation. Figure 16 shows scatterplots of gridcell annual mean ET from RS-PM versus VIC (black circles) and RS-PM versus VIC minus canopy evaporation (gray circles). The latter is perhaps a fairer comparison with the RS-PM data, because it compares transpiration plus soil evaporation only. Notice that the RS-PM and VIC data (black circles) are well correlated (r = 0.81), which reflects the spatial agreement in Fig. 13, but there is some spread in the values (rmse = 0.56 mm day−1). The slight low bias in the RS-PM is −0.21 mm day−1. For the comparison with the VIC data minus canopy evaporation (gray circles), there is a large reduction in the VIC data. The correlation does not change much (r = 0.75), but the RS-PM is now larger on average than the VIC data (bias = 0.36 mm day−1) with essentially the same spread (rmse = 0.57 mm day−1). This suggests that the RS-PM method is actually overestimating transpiration plus soil evaporation, but the fact that it lies somewhere in between the two VIC datasets is somewhat encouraging.
This paper describes the development of a method for estimating high-resolution ET from mostly remote sensing data sources and applies this to generate a long-term dataset for Mexico that can help answer questions about water variability and change. The method is based on the modified Penman–Monteith algorithm of Mu et al. (2007) that is focused on remote sensing applications. Key to estimating ET is determining the value of canopy resistance, which is scaled from stomatal resistance using LAI and adjusted to reflect the daily variation of humidity and temperature. Daily inputs of surface radiation and meteorological data are taken from the ISCCP database, which are derived from satellite radiances and modeled data. Surface characteristics, including albedo and emissivity, are taken from the ISCCP, and vegetation cover and LAI are resampled from AVHRR products. The ISCCP data are downscaled from 280 km to ⅛° using a simple weighting scheme based on the NARR data. The downscaled data are a significant improvement over that obtained from simple interpolation methods, especially for highly spatially variable data such as albedo, and show spatial detail that reflects the heterogeneity of the landscape.
Using the RS-PM method and downscaled ISCCP input data, daily time series of PE, ETrc, and ET were generated at ⅛° resolution for 1984–2006 for the whole of Mexico. The data were evaluated against other large-scale estimates from offline land surface modeling (VIC), atmospheric reanalysis (NARR), and pan evaporation measurements. These evaluations indicate good agreement in the seasonal cycle and spatial distribution of both PE and ET with a slight underestimation by the RS-PM. For PE, the RS-PM is lower on average by less than 0.5 mm day−1 compared with the NARR, despite the ISCCP net radiation being higher by 28%. This may be in part attributed to the differences between the RS-PM and the NARR in the formulation for aerodynamic resistance and accounting for atmospheric stability in the NARR. For ET, the RS-PM has a low bias in the summer/autumn relative to VIC and NARR and a slight high bias in the spring. Overall, there is general agreement in the amplitude and phase of the seasonal cycle, with the largest differences restricted to warmer and more humid regions. The monthly bias and rmse relative to VIC (which may be our best estimate of evapotranspiration, at least annually) over the whole of Mexico are −0.32 and 0.48 mm day−1, respectively. Comparison against small basin water balance measurements indicates that the VIC ET is reasonable and confirms the underestimation by the RS-PM. The NARR results are in contrast to comparisons over United States, which indicate that the NARR evapotranspiration appears to be biased high relative to other estimates from offline modeling and inferred values (Nigam and Ruiz-Barradas 2006; Sheffield et al. 2009).
Despite the reasonable match to our best estimates of regional ET, uncertainties in the input data as well as uncertainties in the true value of ET contribute to uncertainty in the RS-PM estimate and its evaluation. An analysis of the uncertainty as a result of using alternative sources of input net radiation (SRB) and vegetation parameters (MODIS) shows that 95% confidence limits on the RS-PM estimates can vary by as much as 0.49 mm day−1. Uncertainty in the true value of ET as estimated by the spread between the NARR and VIC is much larger, and a comparison of the RS-PM and NARR/VIC uncertainties shows considerable overlap in all months but still a tendency for overestimation in the spring and underestimation in the late summer and autumn. These remaining differences may be in part due to the physical relationship between evapotranspiration and soil moisture, which is not directly accounted for in the RS-PM approach, but they may also be the result of biases in the NARR and VIC data. Further analysis (including assessments of variability at interannual and submonthly time scales, and error propagation and enhancement in the downscaling) and comparisons with other data sources (such as high-resolution retrievals based on Landsat data) are required to fully address the differences. One potential drawback of the RS-PM method is that it does not calculate direct evaporation of water on the vegetation canopy, which can be a significant part of the total evapotranspiration and thus may further explain part of these differences. There is potential to include a parameterization for canopy evaporation based on remote sensing data, either using current inputs [e.g., as a function of relative humidity as in Fisher et al. (2008)], or a more direct approach based on the canopy storage capacity (e.g., Liang et al. 1994) using remotely sensed precipitation.
To the authors’ knowledge, this is the first regional, high-resolution, long-term dataset of evapotranspiration derived from remote sensing data. Given the importance of understanding the variability of ET over large scales, this presents an advance in demonstrating the utility of using remote sensing to quantify the hydrologic cycle and to help answer questions about water cycle variability and change. The 23 years of the dataset allows assessment of the interannual variability in evapotranspiration over Mexico and the analysis of trends in evapotranspiration under the constraints of uncertainty in the temporal consistency of the input data and assumptions about land use change and anthropogenic influences. This will be the subject of a future study. The availability of ISCCP data for the globe naturally leads to the potential of developing estimates in other regions and globally, and this is the subject of continuing work. In the absence of high-resolution data (e.g., NARR) outside of North America, the question arises as to how the input data can be downscaled. An alternative is to use physiographic data as the basis for downscaling, such as elevation, to provide lapse-rate-based adjustments of meteorological data (Sheffield et al. 2006).
This work was supported by the National Aeronautics and Space Administration (Grants NNX08AN40A and NNH07ZDA001N-NEWS) and the World Bank through the Mexican Ministry of Environment and Natural Resources. We thank the Servicio Meteorológico Nacional of Mexico for providing the station data. The ISCCP data were obtained online (available at http://isccp.giss.nasa.gov) and the SRB data were obtained from the NASA Langley Research Center’s Atmospheric Sciences Data Center. We also thank the three anonymous reviewers for their helpful comments.
Corresponding author address: Justin Sheffield, Department of Civil and Environmental Engineering, Princeton University, Princeton, NJ 08544. Email: email@example.com