Evaluation of optical remote sensing to estimate actual evapotranspiration and canopy conductance

Remote Sensing of Environment 129 (2013) 250–261

Contents lists available at SciVerse ScienceDirect

Remote Sensing of Environment journal homepage: www.elsevier.com/locate/rse

Evaluation of optical remote sensing to estimate actual evapotranspiration and canopy conductance Marta Yebra a,⁎, Albert Van Dijk a, b, Ray Leuning c, Alfredo Huete d, Juan Pablo Guerschman a a

CSIRO Water for a Healthy Country Flagship, CSIRO Land and Water, Canberra, ACT, 2601, Australia Fenner School for Environment and Society, The Australian National University, Canberra, 2601, Australia CSIRO Marine and Atmospheric Research, Canberra, ACT, 2601, Australia d University of Technology Sydney, Broadway, NSW, 2007, Australia b c

a r t i c l e

i n f o

Article history: Received 14 July 2012 Received in revised form 30 October 2012 Accepted 6 November 2012 Available online 11 December 2012 Keywords: Surface conductance Penman Monteith MODIS ET EF Vegetation indices LAI fPAR

a b s t r a c t We compared estimates of actual evapotranspiration (ET) produced with six different vegetation measures derived from the MODerate resolution Imaging Spectroradiometer (MODIS) and three contrasting estimation approaches using measurements from eddy covariance flux towers at 16 FLUXNET sites located over six different land cover types. The aim was to assess optimal approaches in using optical remote sensing to estimate ET. The first two approaches directly regressed various MODIS vegetation indices (VIs) and products such as leaf area index (LAI) and fraction of photosynthetically active radiation (fPAR) with ET and evaporative fraction (EF). In the third approach, the Penman–Monteith (PM) equation was inverted to obtain surface conductance (Gs), for dry plant canopies. The Gs values were then regressed against the MODIS data products and used to parameterize the PM equation for retrievals of ET. Jack-Knife cross-validation was used to evaluate the various regression models against observed ET. The PM-Gs approach provided the lowest root mean square error (RMSE), and highest determination coefficients (R2) across all sites, with an average RMSE = 38 W m −2 and R 2 = 0.72. Direct regressions of observed ET against the VIs resulted in an average RMSE = 60 W m −2 and R 2 = 0.22, while the EF regressions an average RMSE = 42 W m−2 and R2 = 0.64. The MODIS LAI and fPAR product produced the poorest estimates of ET (RMSE > 44 W m −2 and R2 b 0.6); while the VIs each performed best for some of the land cover types. The enhanced vegetation index (EVI) produced the best ET estimates for evergreen needleleaf forest (RMSE = 28.4 W m −2, R2 = 0.66). The normalized difference vegetation index (NDVI) best estimated ET in grassland (RMSE = 23.8 W m −2 and R2 = 0.68), cropland (RMSE = 29.2 W m −2 and R 2 = 0.86) and woody savannas (RMSE = 25.4 W m−2 and R2 = 0.82), while the VI-based crop coefficient (Kc) yielded the best estimates for evergreen and deciduous broadleaf forests (RMSE = 27 W m−2 and R2 = 0.7 in both cases). Using the ensemble-average of ET as estimated using NDVI, EVI and Kc we computed global grids of dry canopy conductance (Gc) from which annual statistics were extracted to characterise different functional types. The resulting Gc values can be used to parameterize land surface models. © 2012 Elsevier Inc. All rights reserved.

1. Introduction Remote sensing is the only feasible means of spatially estimating actual evapotranspiration (ET) over large regions or continents. Various approaches developed to derive ET from remote sensing data can be broadly grouped into: (i) those incorporating satellite land surface temperature into a surface energy balance (SEB) model (Kalma et al., 2008), (ii) those using vegetation indices (VIs) (Glenn et al., 2010, 2011b) and (iii) hybrid methods that combine the surface temperature and vegetation index data (Carlson, 2007; Tang et al., 2010).

⁎ Corresponding author at: GPO Box 1666, Canberra ACT 2601, Australia. Tel.: +61 2 6246 5742; fax: +61 2 6246 5988. E-mail address: [email protected] (M. Yebra). 0034-4257/$ – see front matter © 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.rse.2012.11.004

VI based approaches are increasingly being explored, partly because SEB methods have been difficult to implement over large areas and because the number of satellite sensors that have thermal infrared bands is still limited (Glenn et al., 2010). Moreover the surface and near surface meteorological variables at the specific time that remotely sensed data is acquired required to solve a SEB model (e.g. air temperature, relative humidity and solar radiation) are difficult to obtain from daily meteorological data which highly complicates data processing (McVicar & Jupp, 1999). VI methods depend on an estimate of the density of green vegetation over the landscape (Glenn et al., 2010). Although VIs cannot detect soil evaporation nor vegetation stress except on a long time basis (Kalma et al., 2008) several studies have found that they provide better estimations of ET than thermal bands. For example, Cleugh et al. (2007) compared MODerate resolution Imaging Spectroradiometer

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

(MODIS)-based SEB and VI methods against ground measurements of ET in Australia. The SEB methods failed because small errors in land surface temperatures translated into large errors in estimates of sensible heat in the SEB equation, and hence in ET. By contrast, the VI model adequately estimated ET. Similarly, a recent study by King et al. (2011) and summarized by Glenn et al. (2011a) compared different remote sensing-based ET methods, including those based on thermal imagery (McVicar & Jupp, 1999), VIs (Guerschman et al., 2009), MODIS derived LAI (Zhang et al., 2010) and multiple remote sensing data sources (Mu et al., 2007) and concluded that the best performing method was that one based on VIs (RMSE of 0.65 against 0.87 mm d −1for the thermal method). Commonly used VIs include the normalized difference vegetation index (NDVI) (Fisher et al., 2008; Zhang et al., 2009), the enhanced vegetation index (EVI) (Leuning et al., 2005; Mu et al., 2007; Yuan et al., 2010), the normalized difference water index (NDWI) (Lu & Zhuang, 2010) and modelled satellite products such as leaf area index (LAI) (Cleugh et al., 2007; Leuning et al., 2008; Mu et al., 2007) and the fraction absorbed photosynthetically active radiation (fPAR) (Van Dijk, 2010). The VIs are typically used in one of two ways: (i) directly, to retrieve ET through an empirical relationship between ground measurements

of ET (typically from flux towers) or evaporative fraction (EF) (Nishida et al., 2003), or (ii) to parameterize the conductance term of the Penman Monteith (PM) equation (Leuning et al., 2008) (see (Glenn et al., 2010; Glenn et al., 2011a) for a comprehensive review of approaches). Despite the success of various VI-based techniques, there is no consensus on the most appropriate way to use optical remote sensing to estimate ET. The main objective of this study was to compare and evaluate the performance of three contrasting approaches and six different MODIS-derived vegetation measures to retrieve ET and thus determine the best use of optical remote sensing to estimate ET across and within land cover types.

2. Methods and data sources The general scheme of the method developed in this paper is presented in Fig. 1. Meteorological and flux data derived from eddy covariance flux towers as well as reflectances derived from MODIS were used to test three contrasting approaches to estimate ET; (i) direct regression, (ii) potential evapotranspiration (PET) scaling and (iii) PM conductance approach (PM-Gs). Each approach was tested using six different MODIS-derived vegetation measures. The estimates

Eddy covariance Flux Tower

Fluxes and Meteorological data

Meteorological variables

ET

ET/A

Evaporative fraction (EF) Approaches Tested

251

Direct regression ET = a + bVI

Inverted PM-Equation

Gs

PET scaling

PM-Gs

EF = a + bVI

Gs = a exp[b(VI − VImin )]

Remote Sensing data NDVI, EVI, Kc, NDWI

Data source

MOD15A2 (LAI, fPAR)

MCD43A4 (reflectance)

Input Process Output

MODIS

Fig. 1. Methodological flowchart. MODIS Nadir BRDF-adjusted reflectance (MCD43A4) and LAI/fPAR products (MOD15A2) were combined with flux and meteorological data average over 16-days to derive ET (direct regression), EF (PET scaling) and surface conductance (PM-Gs).

252

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

were evaluated against observed ET from the same eddy covariance flux towers where the models were calibrated. 2.1. Fluxes, meteorological data and sites

Table 2 Spectral indices calculated from MODIS including their shortened acronym, mathematical formulation and reference. ρx is the reflectance in MODIS band × (1 to 7). Index

The publicly available ‘free fair-use’ FLUXNET dataset (a subset of the LaThuile dataset 2007, http://www.fluxdata.org/default.aspx, last accessed February 2012) was used to calibrate and validate the approaches described in Section 2.3 to retrieve ET. Half-hourly flux-meteorological data were selected. The data were assed for quality and gap-filled using techniques described elsewhere (Moffat et al., 2007; Papale & Valentini, 2003; Papale et al., 2006; Reichstein et al., 2005). From the total of data available (80 sites) we selected sites with (i) more than 5 years of data (since the 2000 MODIS-Terra launch) and (ii) homogeneous land cover within the lowest spatial resolution of the MODIS products used in this study (1 km). The homogeneity within a 500 m radius footprint around the flux tower was assessed using high spatial resolution aerial and satellite images from various sources (Google Earth™ http://earth.google.com, last accessed February 2012). A site was selected if 75% of the footprint was of the same land cover, judged by colour and texture. 16 sites fulfilled both selection criteria (Table 1). The sites are located across several continents and a wide range of land cover types including crops (CRO), deciduous broadleaf forest (DBF), evergreen needleleaf (ENF) and broadleaf forest (EBF), woody savannas (WSA) and grasslands (GRA) following the IGBP (International Geosphere–Biosphere Programme) classification scheme (Hansen, 2000). 2.2. Remote sensing data Since the launch of the MODIS on board the Terra and Aqua satellites, this sensor has become an invaluable source of earth observation data. Specifically, we used the MODIS/Terra + Aqua 16-day, 500 m, nadir reflectance product (MCD43A4) that is corrected for the bidirectional reflectance distribution function (BRDF) (Strahler et al., 1999) to compute four VIs commonly used to derive ET by various authors; NDVI, EVI, NDWI and Kc (Fig. 1, Table 2). NDVI is a normalized ratio of the near infrared (NIR) and red bands and is sensitive to chlorophyll. Despite its many successful uses, it exhibits scaling problems, asymptotic (saturated) signals over high biomass conditions, and is very sensitive to canopy background variations with brighter canopies negatively-biasing NDVI values (Huete et al., 2002b). EVI was developed to optimize the vegetation signal with improved sensitivity in

Normalized difference vegetation index Enhanced vegetation index Normalized difference water index Crop factor

Global vegetation moisture index

Formulation

Reference

ρ −ρ1 NDVI ¼ 2 ρ2 þ ρ1 2:5  ðρ2 −ρ1 Þ EVI ¼ ðρ2 þ 6  ρ1 −7:5  ρ3 þ 1Þ

Rouse et al. (1974) Huete et al. (2002b)

NDWI ¼

ρ2 −ρ5 ρ2 þ ρ5

Gao (1996)

Kc =Kc _max ×[1−exp(−a×EVIr∝ −b×RMIβ)], where Kc_max =0.68, a=14.12, α=2.482, b=7.991, ß=0.890, EVI−EVImin , KRMI =0.775 and EVIr ¼ EVI max −EVImin CRMI =−0.076 ðρ þ 0:1Þ−ðρ6 þ 0:02Þ GVMI ¼ 2 ðρ2 þ 0:1Þ þ ðρ6 þ 0:02Þ

Guerschman et al. (2009), Model 2b

Ceccato et al. (2002)

high biomass regions and improved vegetation monitoring through a de-coupling of the canopy background signal and a reduction in atmosphere influences (Huete et al., 2002b). Consequently, EVI is more responsive to canopy structural variations, including LAI, canopy type, plant physiognomy, and canopy architecture (Gao et al., 2000). NDWI combines two NIR channels centred approximately at 0.86 μm and 1.24 μm. Both channels are located in the high reflectance plateau of vegetation canopies and sense similar depths through vegetation canopies. However, absorption by vegetation liquid water near 0.86 μm is negligible while moderate at 1.24 μm. As a result, NDWI is sensitive to changes in liquid water content of vegetation canopies. Similarly to NDVI, NDWI does not completely remove the effects of background soil reflectance (Gao, 1996). Finally, Kc is a VI-based crop coefficient developed to estimate ET in mixed biome types at the regional scale (Guerschman et al., 2009). A crop coefficient relates the actual ET at a given stage of crop development to reference ET (ETo) calculated from meteorological data by one of a number of different equations, depending upon data availability (Glenn et al., 2011b). This crop coefficient makes use of MODIS-derived global vegetation moisture index (GVMI) and EVI. The authors included GVMI as a proxy of surface water since it combines water absorption in the short wave infrared (SWIR) wavelengths with NIR that is insensitive to leaf water. Both EVI and GVMI were combined with the Priestley–Taylor estimates of potential ET (PET) to determine ET over Australia (Guerschman et al., 2009).

Table 1 Description of eddy covariance flux tower sites used in this study. h: Canopy height, z: Measurement height. CRO: Crops; DBF: deciduous broadleaf forest; ENF and EBF: evergreen needleleaf and broadleaf forest respectively; WSA: woody savannas; GRA: grasslands. Code

Country

Lat

Lon

h (m)

z (m)

IGBP

Years

Reference

US-Bo1 US-Ne1 US-Ne2 US-Ne3 IT-Ro1 IT-Ro2 US-Ha1 US-MMS US-WCr AU-Tum FI-Hyy NL-Loo US-Ho1 US-Var AU-How US-Ton

United States United States United States United States Italy Italy United States United States United States Australia Finland Netherlands United States United States Australia United States

40.00 41.16 41.16 41.18 42.41 42.39 42.54 39.32 45.81 −35.66 61.85 52.17 45.20 38.41 −12.49 38.43

−88.29 −96.48 −96.47 −96.44 11.93 11.92 −72.17 −86.41 −90.08 148.15 24.29 5.74 −68.74 −120.95 131.15 −120.97

1 3 3 3 15 15 22 25 24 40 14 22 20 1 14 10

8 6 6 6 20 18 31 48 30 70 46 28 29 2 23 23.4

CRO CRO CRO CRO DBF DBF DBF DBF DBF EBF ENF ENF ENF GRA WSA WSA

2000–07 2001–05 2001–05 2001–05 2000–06 2002–06 2000–06 2000–05 2000–06 2001–06 2000–06 2000–06 2000–04a 2001–06 2001–06 2001–06

Meyers and Hollinger(2004) Verma et al. (2005)

a

Rey et al. (2002) Tedeschi et al. (2006) Urbanski et al. (2007) Schmid et al. (2000) Cook et al. (2004) Leuning et al. (2005) Suni et al. (2003) Dolman et al. (2002) Hollinger et al. (2004) Ma et al. (2007) Beringer et al. (2003) Ma et al. (2007)

Data from 2004 was discarded due to a problem with the units in the water flux data that reduces the values by a factor of 1000 (Bill Munger, personal communication).

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

In addition to the four VIs derived from MCD43, LAI and fPAR derived from the MODIS/Terra 8-days, 1000 m, LAI-fPAR product (MOD15A2) (Knyazikhin et al., 1999) were also explored. LAI can be defined as half the total developed area of green elements per unit of horizontal round area (Chen & Black, 1992) and it is widely used in ET algorithms to scale stomatal conductance to the canopy level. fPAR is defined as the fraction of radiation absorbed by the green vegetation elements in the 400–700 nm spectral domain under specified illumination conditions and is one of the main inputs in light use efficiency models (McCallum et al., 2009) and in ET algorithms as a proxy of fractional canopy cover (Los et al., 2000; Mu et al., 2011). Data from both MCD43 and MOD15 were extracted from the MODIS Web service (http://daac.ornl.gov/MODIS/MODIS-menu/modis_ webservice.html, last accessed February 2012). The quality control and state flags were used to remove pixels with partial or complete cloud cover or low pixel quality in the study areas. 2.3. Approaches tested to retrieve ET Three contrasting approaches representative of the main optical remote sensing ET algorithms were tested to estimate ET using the above mentioned VIs (Fig. 1): (i) direct regression, (ii) PET scaling and (iii) PM conductance (PM-Gs). The direct regression approach combined flux tower measurements of ET with time-series of VIs to scale ET over larger landscape units. This method is easy to calibrate and has been successful in a wide variety of land cover types (see examples in Table 1 of Glenn et al. (2010) and the references therein). One difficulty of this approach in local and regional scales is the scaling of instantaneous ET to daily averages (Farah et al., 2004), an issue addressed successfully by Ryu et al. (2012). The evaporative fraction (EF) represents the fraction of available energy used for evaporating water and varies between 0 and 1 under conditions of minimal advection (Macquarrie & Nkemdirim, 1991). Daytime values of EF has been found to have little variations during daytime under a wide range of environmental conditions, so studies since the early 1990s have suggested that remotely sensed data may be useful in estimating EF, which in turn would allow simple modelling of daily ET without a great deal of ancillary data (Kustas et al., 1993). Consequently, the PET scaling approach combined EF with VIs to generate spatial distributions of ET. Finally, the PM-Gs uses an empirical relationship between Gs and different VIs to parameterize the conductance term of the PM equation (Monteith, 1964) (Eq. 1). The remote sensed-derived Gs is subsequently used as input in the PM equation together with other key meteorological drivers to estimate ET. The PM-Gs approach is the most mechanistic VI–ET model and provides the best conceptual basis for ET estimation. The PM equation states that the flux of latent heat associated with ET is given by:

λE ¼


 εA þ ρC ρ =γ DGa ε þ 1 þ Ga =Gs

ð1Þ

where E is the evaporation rate (kg m −2 s −1), λ the latent heat of evaporation (MJ kg −1) Ga is the aerodynamic conductance (m s −1); ε = s/γ, in which s (Pa K −1) is the slope of the saturation vapour pressure versus temperature curve and γ (Pa K −1) is the psychometric constant; A (W m −2) is the available energy absorbed by the surface (net absorbed radiation minus soil heat flux), ρ (kg m −3) air density; Cρ (J kg −1 K −1) the specific heat at constant pressure of air; and D (Pa) the vapour pressure deficit of the air. All meteorological and flux data needed for the calibration of the three approaches were obtained from the eddy covariance towers listed in Table 1. EF was computed as the ratio between ET and A; assumed to be equal to the sum of the measured sensible (H) and

253

latent (λE) heat. This way, the lack of energy balance closure is not assumed to be caused by errors in H alone (Leuning et al., 2008). Effective daytime daily Gs (mm s −1) was calculated by rearranging the PM formulation (Eq. 1) as Gs ¼

λE Ga εA−ðε þ 1ÞλE þ ρcp D=γ

ð2Þ

Ga (m s −1) was calculated from flux tower wind speed measurements and vegetation height as follows (Montheith & Unsworth, 2007) Ga ¼

    1 z−d z−d ln ln z0 z0H k U 2

ð3Þ

where k is the von Karman constant (0.40), U wind speed (m s −1) at the measurement height (z) (m), d the zero-plane displacement (m) and z0 and z0H the roughness lengths for momentum and heat (m), respectively. The quantities d, z0 and z0H were estimated as d = 2 h/3, 0.123 h and 0.0123 h, respectively, where h is canopy height (m). Both z and h were site specific (Table 1). Meteorological and flux data were averaged over the daylight hours (defined as intervals with incoming shortwave radiation > 5 W m −2) before being used as input to Eq. (2). In this way, extreme values of conductance associated with low light, dew or high relative humidity were avoided. Additionally, a day was accepted in the analysis only if most of the noon data was not missing and had not received precipitation during the two preceding days. The latter criteria avoided the influence of evaporation of water on the canopy or soil surface. As a result the Gs values were dominated by dry surface conductance. The resulting daily Gs was further averaged for the MODIS 16-day compositing periods. The pooled dataset contained 573 16-day dry periods that were used to calibrate the relationships between the VIs and ET, EF and Gs (Table 4). ET and EF were fitted using linear regression analysis whereas Gs was better described by an exponential function of VIs. The VI value representative of bare soil (obtained from the literature) was subtracted from the observed values to force the model to estimate Gs ≈ 0 for bare soil. Model parameters were determined using the least squares minimization function MPFITFUN. This function is an implementation of the Levenberg–Marquardt nonlinear optimization algorithm adapted to IDL (Interactive Data Language, ITT Visual Information Solution, Inc.) by Markwardt (2008). The code is available from http://purl.com/net/mpfit (last, accessed May 2012). 2.4. Algorithm evaluation Predicted ET values were compared with measurements for each study site. The coefficient of determination (R 2) and the root mean square error (RMSE) were chosen as accuracy metrics. The total RMSE was decomposed into its systematic (RMSEs) and unsystematic (RMSEu) portions according to Willmott, (1982) (Table 3). A ‘good model’ is considered to have a high R 2 value, a low total RMSE and a RMSEs ≈ 0. The resulting ET estimates were first evaluated with data pooled across sites. Subsequently, the best performing approach and MODISderived vegetation measures were tested by site. The ideal model would accurately predict ET for any site. Subsequently, a Jack-Knife cross-validation was performed by excluding each respective land cover class from the calibration data and calculating the ET estimated for that class with the parameters obtained after calibrating the model against the remaining five classes. This ensures an independent validation of the trained models. Finally, our results were compared to other published studies in the discussion. These studies provide RMSE estimates against mean flux tower ET estimates in various units. In order to make the RMSE comparable, we converted all RMSE estimates to RMSE in mm mo−1;

254

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

Table 3 Measures used for model performance evaluation. A ‘good model’ is considered to have high R2, low RMSE, RMSEu ≈ RMSE and RMSEs ≈ 0. Accuracy assessment

Formulation rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ∑Ni¼1 ðP i −Oi Þ2

Root mean square error (RMSE)

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 ∑N P^ i −Oi

Systematic root mean square error (RMSEs)

i¼1

N

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
2 ∑N P i −P^ i

Unsystematic root mean square error (RMSEu)

3.2. Model performance evaluation

Acronyms Oi: observed value Pi: model-predicted value N: number of observations P^ i : estimated value based on the ordinary least-squares regression P^ i ¼ a þ bOi

N

i¼1

N

where temporal up-scaling was necessary we made two extreme assumptions to calculate a range, viz. errors are (i) entirely random (implying a linear relation between error and the square-root of sample number n) or (ii) entirely systematic (implying error increases proportionally to n). 3. Results 3.1. Alternative ET models Eighteen models (3 approaches times 6 vegetation measures) were calibrated with 573 site-periods from 16 towers over 8 years (Table 4). Direct regression and PET scaling approaches resulted in low-performing models in terms of variance in ET and EF explained. The models derived from the PM-Gs approach had consistently a better fit, with R 2 > 0.32 for all vegetation measures except LAI and fPAR which explained about 10% of the variance in Gs, ET and EF. EVI behaved differently from the other VIs (Table 4), as it explained similar variance in ET and Gs (around 45%, R 2 ≈ 0.4,) while explaining 61% (R 2 = 0.61) of the variance in EF. No unique ET–VI or EF–VI relationship was found across land cover types while a more global pattern was found for the relationships between Gs and all MOD43A4 derived VIs. Fig. 2 shows scatter plots between 16-day averages of observed ET, EF and Gs versus 16-day NDVI and EVI. ET is positively correlated to NDVI and EVI for all land cover types except EBF, where the correlation with NDVI is strongly negative and for ENF where there is no correlation at all. Positive correlations between EF and NDVI or EVI were obtained for WSA, CRO, GRA and DBF (Fig. 2b). The data for these land cover types were less scattered than in the case of ET–NDVI or ET–EVI (Fig. 2a) and the fitted equations for each land cover presented similar slopes. However, as in the previous case, the EF–NDVI and EF–EVI correlations were weak. Finally, all studied land cover types presented similar pattern when looking at the association between Gs and NDVI or EVI (Fig. 2c) with the exception of EBF for EVI. The relationship between ET, EF or Gs and the 6 VIs was ambiguous and uninterpretable when data for wet canopies were included in the Table 4 Parameter estimates as well as coefficients of determination (R2) for the 18 models (3 approaches × 6 VIs) calibrated from different approaches and VIs models (data pooled across sites). VI

NDVI EVI Kc NDWI LAI fPAR a

Direct regression ET = a + bVI

PET scaling EF = a + bVI

PM-Gs Gs = aexp[b(VI − VImin)]

a

b

R2

a

b

R2

a

b

VImina

R2

37.39 5.73 −17.46 −18.36 58.09 45.92

242.30 347.77 210.56 460.88 20.02 87.56

0.21 0.45 0.22 0.27 0.10 0.05

0.20 0.087 0.03 −0.02 0.36 0.28

1.03 1.40 0.78 1.9 0.04 0.28

0.31 0.61 0.25 0.38 0.04 0.04

2.0 2.5 0.3 1.2 4.3 2.7

4.11 3.15 5.14 6.13 0.18 1.46

0.4 0.1 0 −0.3 0.1 0.07

0.39 0.44 0.32 0.36 0.11 0.10

VImin is the VI observed for bare soil.

analysis. The VIs explained b13% of the variance in ET, EF or Gs and validation of the calibrated models against ET presented RMSE>100 W m−2.

Model performance generally improved from direct regression to PET scaling and the PM-Gs approaches in terms of ET estimations (Fig. 3). The R 2 values were lower than 0.5 and the RMSE and RMSEs greater than 55 and 37 W m −2, respectively, for the direct regression approach for all VIs. The PET approach provided better results than direct regression, with R 2 > 0.5, RMSE > 32.5 W m −2 and RMSEs 7–13 W m −2 b RMSE for all VIs. EVI offered better results than the other VIs in this approach (RMSE = 32.5 and RMSEs = 19.6 W m −2). The PM-Gs approach gave the best overall results (lowest RMSE, RMSEs and highest R 2) for all VIs, satisfying the requirement of a ‘good model’. The best performing VIs were those derived from MCD43 (NDVI, Kc, NDWI and EVI) with overall RMSE b 35 W m −2 and R 2 > 0.8, while LAI and fPAR, derived from MOD15, resulted in RMSE > 40 W m −2 and R 2 around 0.6. Based on the above comparison, the PM-Gs approach in combination with VIs derived from MCD43 was selected for further analysis. Contrary to the validation across sites where none of the MCD43 VIs outperformed the others, using the criteria of a ‘good model’, the validation per site showed NDVI best estimated ET at four sites whereas EVI and Kc were most successful for 3 and 5 sites, respectively (Fig. 4 and Table 5). Finally, NDWI gave similar results as the other VIs at most sites. Consequently, there was no clear advantage in using this VI over EVI, NDVI or Kc. Cross-validation analysis revealed that each VI performed best for some land cover types (Fig. 5 and Table 6). NDVI performed best for CRO, GRA and WSA, explaining 86%, 68% and 82% of the variance in ET, respectively. The RMSE was b 30 W m −2 and the RMSEs was much lower than the RMSEu. EVI performed the best for ENF while Kc was best for EBF and DBF. The validation per site analysis also showed that EVI was the best VI for estimating ET in ENF (Table 5). As none of the best performing VIs were uniformly superior across land cover types, we calculated an ensemble-average of ET as estimated using NDVI, EVI and Kc. By averaging the results of the various models, we follow best practice in environmental modelling (e.g. weather forecasting, climate modelling, water resources modelling, (Collins, 2007; Hagedorn et al., 2005). The underlying argument is that each of the three ensemble members performed best for some of the land cover types, suggesting that part of the error is independent among the methods used and thus averaging improves the overall estimate. Fig. 6 presents the validation results of averaging the ET estimations obtained from NDVI, EVI and Kc. The ensemble of the models accounted for 82% of the variance in ET, with RMSE of 29.8 W m −2, while the average variance accounted for any of these VIs individually was 78% with an average RMSE of 33 W m −2 (31.7, 35.1, and 32.1 for NDVI, EVI and Kc, respectively). 4. Discussion 4.1. Model performance analysis We found that the use of VIs to parameterize the conductance term of the PM equation (PM-Gs approach) provides ET estimates that are better than estimates derived from direct regressions between VIs and measured ET or EF. A unique ET–VI relationship across land cover types was not found, whereas a more global pattern was found for the Gs–VI relationships. Consequently, a global ET model across land cover types cannot be fit based on the relationships between observed ET or EF and VIs, while the PM-Gs approach presented high potential, as reflected by the higher R2 of the calibrated models (Table 4). The VIs by themselves are not sufficient to estimate ET or EF. Arguably, the PM-Gs approach is conceptually the most preferable

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

a)

300

R2 = 0.21 y = 242.324x + 37.39

R2 = 0.45 y = 347.84x + 5.7

WSA ENF EBF DBF GRA CRO

ET (Wm-2)

200

255

100

0

b)

1.5

R2 = 0.31 y = 1.03x + 0.2

R2 = 0.61 y = 1.41x + 0.09

R2 = 0.35 y = 2*exp(4.11x)

R2 = 0.43 y = 2.5*exp(3.15x)

EF

1.0

0.5

0.0

Gs (mms-1)

c)

30

20

10

0 0.0

0.2

0.4

0.6

NDVI-NDVImin

0.8 .0

0.2

0.4

0.6

0.8

EVI-EVImin

Fig. 2. Scatter plots of 16-days observed averages of a) ET, b) EF and c) Gs versus 16-days NDVI (left) and EVI (right) derived from MODIS.

since VIs are used to derive Gs, which is then combined with meteorological data within the PM framework. An exception was found with EVI. EVI was the only VI where EF estimations provided ET values as accurate as the PM-Gs approach. These findings disagree with Sjöström et al. (2011) who found that EVI only accounted for a negligible to moderate part of the variability in EF during the growing season (R 2 ≤ 0.50) across African ecosystems. The only individual site where the authors observed R 2 similar to those reported in our study was a grassland site (R 2 = 0.69). However, Nagler et al. (2007) derived significant (P b 0.01, R 2 not reported) linear regression equations between EF and EVI over sparsely vegetated grasslands and shrubland sites in a semiarid watershed in south eastern Arizona. This implies that ET can still be retrieved with a reasonable precision directly from remote sensing data when meteorological data are not available.

Wang et al. (2006) observed a general increase in EF with the NDVI for 11 sites in the Southern Great Plains of the United States, but they found that the regression slope varies from sites to site, in agreement with our findings. Several studies have shown that MODIS EVI is better correlated with ground-based ET measurements than MODIS NDVI across a variety of land cover types. For example, Nagler et al. (2005) reported correlation coefficients of 0.76 and 0.68 for EVI and NDVI with ET measured at four flux tower riparian sites on the Middle Rio Grande (New Mexico). Nagler et al. (2007) also found stronger correlations of ET with EVI than NDVI over sparsely vegetated grasslands and shrublands in a semiarid watershed in south-eastern Arizona, while Nagler et al. (2009) found that the correlation between ET and EVI was significant (P = 0.036) but of low predictive power (R 2 = 0.55), similar to the R 2 reported in this study.

256

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

0.9

70

0.8 0.7

50

0.6

40

0.5

30

0.4

R2

RMSE (W m-2)

60

0.3

20

0.2 10

0.1

0

Direct regression

PET scaling

RMSE

fPAR

LAI

NDWI

Kc

EVI

NDVI

fPAR

LAI

NDWI

Kc

EVI

NDVI

fPAR

LAI

NDWI

Kc

EVI

NDVI

0

PM-Gs

R2

RMSEs

Fig. 3. Evaluation of the calibrated ET models against observed ET across sites. Shown are the total and systematic root mean square error (RMSE, RMSEs, respectively in W m−2, left axis), as well as the R2 values between estimated and observed ET (right-hand axis). The best models are those with the lowest RMSE (squares), RMSEs ≈ 0 (triangles) and the highest R2 values (asterisks). Model performance improves from direct regression to PET scaling and the PM-Gs approaches. The best performing VIs were those derived from MCD43 with RMSE b 35 W m−2 and R2 > 0.8.

Within all the VI–Gs models, those using VIs from MCD43 reflectance product (NDVI, Kc, NDWI and EVI) estimated ET more accurately than the models derived from MOD15 vegetation product (LAI and fPAR). This suggests that it is not only leaf area or cover fraction but also greenness per unit of area that contributes to ET predictions. Alternatively, the overall poor performance of LAI and fPAR to estimate ET using any of the tested approaches can be explained by

erroneous estimates of LAI, especially over evergreen needleleaf and mixed forests (Garrigues et al., 2008; Shabanov et al., 2005) or errors in the retrievals due to the short and noisier 8-day composite period of MOD15. Finally, the different spatial resolution of LAI/fPAR (1 km) and reflectance (500 m) MODIS products also enhances the performance of the latter. Our findings have implications for the large number of algorithms that use MODIS derived LAI to scale up leaf level

60

RMSE (Wm-2)

50 40 30 20 10

EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI

0

US-Bo1

US-Ne1

US-Ne2

US-Ne3

IT-Ro1

IT-Ro2

US-Ha1

US-MMS

RMSE (Wm-2)

50 40 30 20 10

EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI EVI Kc NDVI NDWI

0

US-WCr

AU-Tum

FI-Hyy

NL-Loo

RMSE

US-Ho1

US-Var

AU-How

US-Ton

RMSEs

Fig. 4. Evaluation of PM-Gs approach against observed ET and pooling data by site. Only VIs derived from MCD43 were used. Shown are RMSE and RMSEs, in W m−2 (left-hand side) as well as the R2 between estimated and observed ET (right-hand side). Dotted line at 35 W m−2 is the average accuracy of the methods resulting from the analysis using data pooled across sites. For some sites one VI clearly outperformed the others. The best VI for each site is shaded lowest RMSE (filled squares) and RMSE > RMSEs (triangles).

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261 Table 5 Accuracy assessment for the model which best estimates ET for each studied site. Shown are RMSE, RMSEs and RMSEu (W m−2), as well as the R2 between estimated and observed ET. N is the number of 16-day periods with valid tower and MODIS derived measurements.

257

Table 6 Accuracy assessment for the model which best estimates ET for each vegetation type. Shown are RMSE, RMSEs and RMSEu (W m−2), as well as the R2 between estimated and observed ET. N is the number of 16-day periods with valid tower and MODIS derived measurements.

IGBP

Site

Best index

N

RMSE

RMSEu

RMSEs

R2

IGBP

Best index

N

RMSE

RMSEu

RMSEs

R2

CRO

US-Bo1 US-Ne1 US-Ne2 US-Ne3 IT-Ro1 IT-Ro2 US-Ha1 US-MMS US-WCr AU-Tum FI-Hyy NL-Loo US-Ho1 US-Var AU-How US-Ton

NDVI NDWI EVI Kc NDWIa Kc NDWI Kc NDVI Kc EVI EVI EVIb NDVI Kc NDVI

28 22 25 29 42 20 37 30 23 60 23 36 27 47 64 60

23.9 28.2 24.5 28.2 24.0 17.8 31.8 27.4 10.0 27.2 23.9 29.7 30.0 23.8 29.8 20.9

20.1 24.8 18.2 26.7 17.4 13.7 25.7 24.2 9.7 22.8 16.5 28.4 14.5 23.2 29.6 18.8

12.9 13.3 16.4 9.1 16.6 11.3 18.7 12.8 2.3 14.7 17.3 8.7 26.3 5.2 4.0 9.2

0.76 0.85 0.91 0.85 0.72 0.87 0.56 0.78 0.95 0.68 0.56 0.57 0.90 0.74 0.79 0.80

CRO DBF EBF ENF GRA WSA

NDVI Kc Kc EVI NDVI NDVI

104 152 60 86 47 124

29.2 27.7 27.2 28.4 23.8 25.4

22.0 25.7 22.8 23.4 23.2 20.3

19.2 10.3 14.7 16.0 5.2 15.3

0.86 0.66 0.70 0.66 0.68 0.82

DBF

EBF ENF

GRA WSA

(Chuvieco et al., 2009; Yebra & Chuvieco, 2009; Yebra et al., 2008). Low leaf chlorophyll concentration is often coordinated with lower stomatal conductance (Matsumoto et al., 2005) and explains the good performance of NDVI to estimate Gs and ET of GRA and WSA since these land cover types present strong seasonality in LAI, chlorophyll content and red reflectance. Such seasonal variations are less evident in high biomass regions such as evergreen forest, where NDVI saturates asymptotically (Huete et al., 2002b). Consequently, the availability of NDVI to resolve changes in Gs and ET in this land functional type is limited. By contrast, EVI saturates less and is more responsive to canopy structural variations and changes in NIR reflectance (Gao et al., 2000; Huete et al., 2002a). This could explain why EVI had better ability to estimate ET in ENF than NDVI. Leaf reflectance in the NIR is determined primarily by properties of tissue and cell structure (Grant, 1987) and changes in NIR reflectance in ENF are likely caused by changes in the internal structure of leaf cells, turgidity of the guard cells and stomatal conductance. Castro-Esau et al. (2006) found strong positive correlations (R2 > 0.7) between NIR reflectance and percentage of intercellular space in nine tropical dry forest tree species. Bowman (1989) showed a strong negative relationship (R2 = 0.72) between NIR and leaf turgor pressure when vegetation is under stress conditions in cotton leaves. Kc combines the EVI and GVMI, an index which like NDWI is responsive to leaf water and soil moisture content to estimate the actual ET from PET (Guerschman et al., 2009). The simultaneous use of both indices makes Kc sensitive to changes in Gs in dense canopies (see Figure 2 in Guerschman et al., 2009). The ensemble-average estimates of ET derived from NDVI, EVI and Kc accounted for more of the variance in ET (82%) and lowest in the RMSE (29.8 W m −2) (Fig. 6) compared to the VIs individually (Fig. 3).

a At these sites RMSEs is very close to RMSEu so the errors caused by the predictors are similar to those caused by uncontrolled factors. b RMSEs >RMSEu but is the only index which reaches RMSE b 35 W m−2 and high R2.

stomatal conductance to the canopy and surface level. Overestimates of LAI in some biomes result in overestimates of ET even if other inputs such as the meteorological data are relatively accurate (Mu et al., 2011). Furthermore, the quality flags (MODIS QA) are insufficient to reduce the weight of unreliable satellite products, especially over tropical forests (Demarty et al., 2007). Some authors have reported that temporal filling of unreliable MODIS data (including LAI, fPAR and albedo) as well as smoothing filters can greatly improve the accuracy of inputs (Zhang et al., 2008). However the filled values are artificial and also contain uncertainties (Mu et al., 2011). Forcing agreement between water and PM-Gs estimates of ET from catchments at annual time scales also minimizes errors introduced by those in LAI and meteorological data (Zhang et al., 2008). None of the VIs derived from MCD43 performed consistently best across all land cover types. EVI produced considerably better ET estimates for ENF, Kc for EBF and DBF, and NDVI for GRA and WSA (Fig. 5). NDVI responds strongly to variations in red reflectance and consequently, is very sensitive to chlorophyll and LAI variations that follow leaf drying in many species, particularly in grasslands

60

1 0.9 0.8 0.7

40

0.6 30

0.5

R2

RMSE (Wm-2)

50

0.4 20

0.3 0.2

10

0.1 0

NDVI

RMSEs

WSA

GRA

ENF

EBF

Kc

EVI

RMSE

DBF

CRO

WSA

GRA

ENF

EBF

DBF

CRO

WSA

GRA

ENF

EBF

DBF

CRO

0

R2

Fig. 5. Evaluation of MCD43-derived VIs to estimate ET over different land cover types. Estimations are derived from the PM-Gs equations obtained by excluding each respective land cover class from the calibration data and calculating the ET estimated for that class with the parameters obtained after calibration the model against the remaining five classes. Dotted line at 35 W m−2 is the average accuracy of the methods resulting from the analysis using data pooled across sites. For some sites one VI clearly outperformed the others. The best VI for each site is shaded lowest RMSE (filled squares) and RMSE >RMSEs (triangles).

258

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

300

systematic and total RMSE that can be inferred from our results (13 vs. 30 W m −2, respectively).

R2= 0.82 y=−8.58+1.01x

ET Observed (Wm−2)

4.2. Implications for large scale ET estimation and land surface modelling 200

WSA ENF EBF DBF GRA CRO

100

0

0

100

200

300

ET Predicted (Wm−2) Fig. 6. Predicted versus measured 16-day ET (W m−2) derived from the ensemble of Gc- NDVI,-EVI,-kc models combined with meteorological data within the PM framework.

RMSE values found in this study are commensurate with those reported previously (Table 7), although a conclusive comparison between studies is not possible due to the many differences between them. These include: i) the calibration and validation procedures (RMSE sometimes reflects site-based or multi-site parameter fitting, sometimes independent cross-validation); ii) in the input data used (the MODIS product and collection, flux tower or large scale gridded meteorological data); iii) in the number and location of field sites considered; and iv) in flux data selection and processing (e.g., in the present study we selected those days expected to be dominated by transpiration only). Systematic differences between model predictions and observations were lower (RMSEs = 13.1 W m −2) than unsystematic errors (RMSEu = 26.7 W m −2). Guerschman et al. (2009) showed that the relative error in ET estimation using their method decreased from around 20% per month to 10% per year if the RMSE was entirely random, compared to a range of 6–20% estimated when the errors are systematic. From these numbers, it can be inferred that about 29% of the total monthly RMSE was systematic and the remainder random. These numbers agree with the ratio of 33% between

Quantifying ET from the global land surface is critical for further understanding of Earth's climate system (Jiménez et al., 2011). Satellite VIs provide spatially comprehensive and temporally frequent observations of surface conditions. We have shown that several VIs correlate well with Gs, allowing the use of remote sensing observations to estimate this parameter and ET under the PM equation framework. One important advantage of the estimation approach developed here is that identical coefficients can be applied globally to all land cover types. This avoids artefacts and errors induced by categorical land cover mapping. A global map of mean dry canopy conductance Gc (mm s −1) (2001–2011) is shown in Fig. 7. We refer here to Gc instead to Gs because the relationships between Gs and NDVI, EVI and Kc derived in this study were for dry canopies only. Gc is highest for tropical rainforests (>9 mm s−1), lowest for arid and semi-arid regions (b4 mm s −1), and intermediate for temperate and boreal forests. In large-scale applications of these Gc estimates to derive ET large-scale grids rather than tower measurements of meteorology are needed. These are readily available, however, and other studies have shown that the degradation in accuracy is typically relatively small (Cleugh et al., 2007; Leuning et al., 2008; Mu et al., 2007). To derive estimates of total ET rather than ET from dry canopies alone, all additional evaporative fluxes need to be considered. This can be achieved by using the method developed here to parameterise Gs in a land surface model and/or by using additional remote sensing observations to estimate soil moisture, open water and wet canopy evaporation rates (see Mueller et al., 2011; van Dijk & Renzullo, 2011, for various examples). Time series of remotely-sensed Gc can be integrated with hydrological models that use the PM equation. For the MODIS period (2000 to present), Gc estimates based on remote sensing time series are the most appropriate to use but for applications outside this period, a monthly or seasonal average pattern is arguably more suitable. To this end, 2001–2011 mean values of global monthly Gc were calculated based on EVI, NDVI and Kc (see Supplementary material; data available via http://eos.csiro.au/). Alternatively, the estimates produced here may be used in parameter look up tables. For the land cover types considered here, some of the Gc distribution statistics are listed in Table 8 as well as indicators of spatial variation in these metrics. Land cover types with strong seasonality such as WSA, DBF and CRO present large differences in maximum and minimum Gs values while evergreen forest does not

Table 7 Accuracy assessment of different method used to estimate ET in the literature. All the reported errors were converted to a common unit (mm mo−1) for the comparison. Where temporal up-scaling was necessary we made two extreme assumptions to calculate a range, viz. (i) errors are entirely random (implying a linear relation between error and the square-root of sample number n) or (ii) entirely systematic (implying error increases proportionally to n). Method

Reference

Time interval (days)

Sites

Reported RMSE

RMSE (mm mo−1)

Universal parameters PM EVI DR EVI DR EVI PM LAI PM LAI EF NDVI, SAVI EF Kc PM LAI

This study Yang et al. (2006) Nagler et al. (2005) Mu et al. (2007) Mu et al. (2011) Fisher et al. (2008) Guerschman et al. (2009) Cleugh et al. (2007)

16 8 16 8 1 30 30 8

16 19 4 19 46 15 5a 2

30 W m−2 0.62 mm d−1 1.6 mm d−1 27.3 W m−2 0.84 mm d−1 16 mm mo−1 17 mm mo−1 27 W m−2

23–31 17–21 35–49 15–29 5–25 16 17 15–29

Site-specific parameters PM LAI

Leuning et al. (2008)

1

15

0.82 mm d−1

4–25

a

VIs

Errors reported excluding Fogg &Hume Dams (a floodplain and a reservoir lake).

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

150° W

120° W

90° W

60° W

30° W

0°

30° E

60° E

259

90° E

120° E

150° E

180°

90°

60° N

60° N

30° N

30° N

0°

0°

30° S

30° S

60° S

60° S

90°

90° 150° W

0

120° W

90° W

1

2

60° W

3

30° W

4

0°

5

30° E

6

60° E

7

90° E

8

120° E

9

150° E

10

2001-2011 Average Gc (mms-1) Fig. 7. Mean global Gc for 2001–2011 derived from the ensemble of NDVI, EVI and Kc calculated from MCD43C4 data.

Table 8 Annual climatology of spatial means of Gc (mm s−1) for different land cover types (IGBP) using a base period of 2001–2011. Shown are the annual mean, 5% percentile (low values of Gc) and 95% percentile (high values of Gc). Gc was retrieved as the average of the estimates using NDVI, EVI and Kc. IGBP

Spatial Statistic

High (95%Percentile)

Mean

Low (5%Percentile)

GRA

Mean Median Q25 Q75 Mean Median Q25 Q75 Mean Median Q25 Q75 Mean Median Q25 Q75 Mean Median Q25 Q75 Mean Median Q25 Q75

3.61 3.02 2.23 4.32 7.85 8.31 6.59 9.41 10.28 10.46 9.72 10.99 10.52 11.22 8.76 12.16 7.46 7.54 6.60 8.31 7.32 7.41 5.61 8.93

2.41 1.99 1.62 2.65 5.34 5.45 4.42 6.37 9.33 9.55 8.59 10.21 6.35 6.43 5.75 6.95 5.89 5.81 5.25 6.39 4.11 3.93 3.10 4.91

1.71 1.41 1.28 1.70 2.77 2.57 2.06 3.31 8.19 8.61 6.95 9.50 2.86 2.67 2.31 3.28 4.77 4.66 4.17 5.17 2.35 2.04 1.63 2.72

WSA

EBF

DBF

ENF

CRO

and soil under optimum meteorological conditions arising from differences in plant physiology and LAI) around 20 mm s −1 for natural vegetation and 30 mm s −1 for agricultural crops. We derive lower maximum Gc values of around 12 mm s −1. This reflects Gc values derived from flux tower ET measurements that were generally lower than those reported by Kelliher et al. (1995), although higher values did exist in the unexplained data scatter (>20 mm s −1, particularly for crops, Fig. 2c). In addition, we selected ‘dry’ days only to avoid

15

Gc (mms-1)

show great variation in these quantities. The annual mean value of Gc is higher in EBF and lowest in GRA and CRO with intermediate values in WSA, DBF and ENF (Fig. 8). Kelliher et al. (1995) reported mean values of maximum bulk surface conductance (Gsmax including plant canopy

5%Percentile Mean 95%Percentile

10

5

0 GRA

WSA

EBF

DBF

ENF

CRO

Fig. 8. Box plots of the spatial variability of 5% percentile (low values), mean and 95% percentile (high values) of annual Gc (mm s−1) for 2001–2011 and different land cover types. Gc was retrieved as the average of values predicted based on NDVI, EVI and Kc calculated from MCD43C4 data.

260

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261

overestimation of Gc due to non-transpiration fluxes, whereas such fluxes may have been influenced the estimates of Kelliher et al. (1995). 5. Conclusions Over recent years, several methods have been published that use VIs to estimate ET but there is as yet no consensus on the most appropriate method to use VIs in ET estimation. We analysed flux tower observations from 16 sites world-wide in conjunction with time series of several spectral VIs and vegetation products derived from the MODIS sensor to test which is the most appropriate way to use optical remote sensing to estimate ET. We conclude that the use of VIs (in particular, EVI, NDVI and Kc) to estimate Gs within the PM framework provides the best basis for ET estimation, when compared to direct regression or PET scaling approaches or the two MODIS vegetation products (LAI and fPAR). PET scaling using EVI was the only alternative approach that produced similarly good results, although it is arguably less desirable from a conceptual point of view. The poor performance of the vegetation products was surprising but may in part be due to noise and bias in the retrievals. This remains to be proven. None of the three most successful VIs (NDVI, EVI and Kc) were uniformly superior across all land cover types. We speculate that this may be because different physiological attributes dominate temporal and spatial Gs patterns in different land cover types. Based on the mean of the three best VIs, we generated gridded and tabulated conductance data dominated by the canopy (Gc) to facilitate application of our results in large scale ET estimation. Further work is needed to extend our conclusions to other land cover types not examined in this work. Additionally, model performance should be tested fully independently by applying the calibrated models to an independent observational dataset. The public FLUXNET dataset is expected to soon double the available data (http://www.fluxdata.org/default.aspx, last accessed April 2012), which will create further opportunities. Acknowledgements This work used eddy covariance data acquired by the FLUXNET community and in particular by the following networks: AmeriFlux (U.S. Department of Energy, Biological and Environmental Research, Terrestrial Carbon Program (DE‐FG02‐04ER63917 and DE‐FG02‐ 04ER63911)), CarboEuropeIP, CarboItaly and OzFlux. We acknowledge the financial support to the eddy covariance data harmonization provided by CarboEuropeIP, FAO‐GTOS‐TCO, iLEAPS, Max Planck Institute for Biogeochemistry, National Science Foundation, University of Tuscia, Université Laval and Environment Canada and the US Department of Energy and the database development and technical support from Berkeley Water Center, Lawrence Berkeley National Laboratory, Microsoft Research eScience, Oak Ridge National Laboratory, University of California — Berkeley, and University of Virginia. We thank Eva Van Gorsel for her help and useful comments in an early phase of the manuscript preparation, David Blondeau-Patissier for providing assistance with IDL programming and David Fanning for the great job he has done to make IDL a bit more friendly (http://www.idlcoyote.com/documents/programs.php). Marta Yebra was supported by a CSIRO OCE postdoctoral scholarship. Appendix A. Supplementary data Supplementary data to this article can be found online at http:// dx.doi.org/10.1016/j.rse.2012.11.004. References Beringer, J., Hutley, L. B., Tapper, N. J., Coutts, A., Kerley, A., & O'Grady, A. P. (2003). Fire impacts on surface heat, moisture and carbon fluxes from a tropical savanna in northern Australia. International Journal of Wildland Fire, 12, 333–340.

Bowman, W. D. (1989). The relationship between leaf water status, gas exchange, and spectral reflectance in cotton leaves. Remote Sensing of Environment, 30, 249–255. Carlson, T. (2007). An overview of the “triangle method” for estimating surface evapotranspiration and soil moisture from satellite imagery. Sensors, 7, 1612–1629. Castro-Esau, K. L., Sánchez-Azofeifa, G. A., Rivard, B., Wright, S. J., & Quesada, M. (2006). Variability in leaf optical properties of Mesoamerican trees and the potential for species classification. American Journal of Botany, 93, 517–530. Ceccato, P., Gobron, N., Flasse, S., Pinty, B., & Tarantola, S. (2002). Designing a spectral index to estimate vegetation water content from remote sensing data: Part 1 theoretical approach. Remote Sensing of Environment, 82, 188–197. Chen, J. M., & Black, T. A. (1992). Defining leaf-area index for non-flat leaves. Plant, Cell & Environment, 15, 421–429. Chuvieco, E., Gonzalez, I., Verdu, F., Aguado, I., & Yebra, M. (2009). Prediction of fire occurrence from live fuel moisture content measurements in a Mediterranean ecosystem. International Journal of Wildland Fire, 18, 430–441. Cleugh, H. A., Leuning, R., Mu, Q. Z., & Running, S. W. (2007). Regional evaporation estimates from flux tower and MODIS satellite data. Remote Sensing of Environment, 106, 285–304. Collins, M. (2007). Ensembles and probabilities: a new era in the prediction of climate change. Philosophical Transactions. Series A: Mathematical, Physical, and Engineering Sciences, 365, 1957–1970. Cook, B. D., Davis, K. J., Wang, W. G., Desai, A., Berger, B. W., Teclaw, R. M., et al. (2004). Carbon exchange and venting anomalies in an upland deciduous forest in northern Wisconsin, USA. Agricultural and Forest Meteorology, 126, 271–295. Demarty, J., Chevallier, F., Friend, A. D., Viovy, N., Piao, S., & Ciais, P. (2007). Assimilation of global MODIS leaf area index retrievals within a terrestrial biosphere model. Geophysical Research Letters, 34, L15402. Dolman, A. J., Moors, E. J., & Elbers, J. A. (2002). The carbon uptake of a mid latitude pine forest growing on sandy soil. Agricultural and Forest Meteorology, 111, 157–170. Farah, H. O., Bastiaanssen, W. G. M., & Feddes, R. A. (2004). Evaluation of the temporal variability of the evaporative fraction in a tropical watershed. International Journal of Applied Earth Observation and Geoinformation, 5, 129–140. Fisher, J. B., Tu, K. P., & Baldocchi, D. D. (2008). Global estimates of the land-atmosphere water flux based on monthly AVHRR and ISLSCP-II data, validated at 16 FLUXNET sites. Remote Sensing of Environment, 112, 901–919. Gao, B. C. (1996). NDWI. A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment, 58, 257–266. Gao, X., Huete, A. R., Ni, W. G., & Miura, T. (2000). Optical-biophysical relationships of vegetation spectra without background contamination. Remote Sensing of Environment, 74, 609–620. Garrigues, S., Lacaze, R., Baret, F., Morisette, J. T., Weiss, M., Nickeson, J. E., et al. (2008). Validation and intercomparison of global Leaf Area Index products derived from remote sensing data. Journal of Geophysical Research, 113, G02028. Glenn, E. P., Doody, T. M., Guerschman, J. P., Huete, A. R., King, E. A., McVicar, T. R., et al. (2011a). Actual evapotranspiration estimation by ground and remote sensing methods: The Australian experience. Hydrological Processes, 25, 4103–4116. Glenn, E., Nagler, P., & Huete, A. (2010). Vegetation index methods for estimating evapotranspiration by remote sensing. Surveys in Geophysics, 31, 531–555. Glenn, E. P., Neale, C. M. U., Hunsaker, D. J., & Nagler, P. L. (2011b). Vegetation index-based crop coefficients to estimate evapotranspiration by remote sensing in agricultural and natural ecosystems. Hydrological Processes, 25, 4050–4062. Grant, L. (1987). Diffuse and specular characteristics of leaf reflectance. Remote Sensing of Environment, 22, 309–322. Guerschman, J. P., Van Dijk, A. I. J. M., Mattersdorf, G., Beringer, J., Hutley, L. B., Leuning, R., et al. (2009). Scaling of potential evapotranspiration with MODIS data reproduces flux observations and catchment water balance observations across Australia. Journal of Hydrology, 369, 107–119. Hagedorn, R., Doblas-Reyes, F. J., & Palmer, T. N. (2005). The rationale behind the success of multi-model ensembles in seasonal forecasting — I. Basic concept. Tellus Series A-Dynamic Meteorology and Oceanography, 57, 219–233. Hansen, M. C. (2000). A comparison of the IGBP DISCover and University of Maryland 1 km global land cover products. International Journal of Remote Sensing, 21, 1365–1373. Hollinger, D. Y., Aber, J., Dail, B., Davidson, E. A., Goltz, S. M., Hughes, H., et al. (2004). Spatial and temporal variability in forest-atmosphere CO2 exchange. Global Change Biology, 10, 1689–1706. Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002a). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83, 195–213. Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X., & Ferreira, L. G. (2002b). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83, 195–213. Jiménez, C., Prigent, C., Mueller, B., Seneviratne, S. I., McCabe, M. F., Wood, E. F., et al. (2011). Global intercomparison of 12 land surface heat flux estimates. Journal of Geophysical Research, 116, D02102. Kalma, J. D., McVicar, T. R., & McCabe, M. F. (2008). Estimating land surface evaporation: A review of methods using remotely sensed surface temperature data. Surveys in Geophysics, 29, 421–469. Kelliher, F. M., Leuning, R., Raupach, M. R., & Schulze, E. D. (1995). Maximum conductances for evaporation from global vegetation types. Agricultural and Forest Meteorology, 73, 1–16. King, E., van Niel, T., van Dijk, A., Wang, Z., Paget, M., Raupach, T., et al. (2011). In C.W.f.a.H.C.N.R. Flagship (Ed.), Actual evapotranspiration estimates for Australia intercomparison and evaluation (Canberra). Knyazikhin, Y., Glassy, J., Privette, J. L., Tian, Y., Lotsch, A., Zhang, Y., et al. (1999). MODIS Leaf Area Index (LAI) And Fraction Of Photosynthetically Active Radiation Absorbed By Vegetation (FPAR) Product (MOD15). Algorithm Theoretical Basis Document. http://eospso.gsfc.nasa.gov/atbd/modistables.html

M. Yebra et al. / Remote Sensing of Environment 129 (2013) 250–261 Kustas, W. P., Schmugge, T. J., Humes, K. S., Jackson, T. J., Parry, R., Weltz, M. A., et al. (1993). Relationships between evaporative fraction and remotely-sensed vegetation index and microwave brightness temperature for semiarid rangelands. Journal of Applied Meteorology, 32, 1781–1790. Leuning, R., Cleugh, H. A., Zegelin, S. J., & Hughes, D. (2005). Carbon and water fluxes over a temperate Eucalyptus forest and a tropical wet/dry savanna in Australia: measurements and comparison with MODIS remote sensing estimates. Agricultural and Forest Meteorology, 129, 151–173. Leuning, R., Zhang, Y. Q., Rajaud, A., Cleugh, H., & Tu, K. (2008). A simple surface conductance model to estimate regional evaporation using MODIS leaf area index and the Penman–Monteith equation. Water Resources Research, 44. Los, S. O., Pollack, N. H., Parris, M. T., Collatz, G. J., Tucker, C. J., Sellers, P. J., et al. (2000). A global 9-yr biophysical land surface dataset from NOAA AVHRR data. Journal of Hydrometeorology, 1, 183–199. Lu, X., & Zhuang, Q. (2010). Evaluating evapotranspiration and water-use efficiency of terrestrial ecosystems in the conterminous United States using MODIS and AmeriFlux data. Remote Sensing of Environment, 114, 1924–1939. Ma, S., Baldocchi, D. D., Xu, L., & Hehn, T. (2007). Inter-annual variability in carbon dioxide exchange of an oak/grass savanna and open grassland in California. Agricultural and Forest Meteorology, 147, 157–171. Macquarrie, P., & Nkemdirim, L. C. (1991). Potential, actual, and equilibrium evapotranspiration in a wheat field. Water Resources Bulletin, 27, 73–82. Markwardt, C. B. (2008). Non-linear least squares fitting in IDL with MPFIT. Quebec, Canada: Astronomical Data Analysis Software and Systems XVIII. Matsumoto, K., Ohta, T., & Tanaka, T. (2005). Dependence of stomatal conductance on leaf chlorophyll concentration and meteorological variables. Agricultural and Forest Meteorology, 132, 44–57. McCallum, I., Wagner, W., Schmullius, C., Shvidenko, A., Obersteiner, M., Fritz, S., et al. (2009). Satellite-based terrestrial production efficiency modeling. Carbon Balance and Management, 4, 8. McVicar, T. R., & Jupp, D. L. B. (1999). Estimating one-time-of-day meteorological data from standard daily data as inputs to thermal remote sensing based energy balance models. Agricultural and Forest Meteorology, 96, 219–238. Meyers, T. P., & Hollinger, S. E. (2004). An assessment of storage terms in the surface energy balance of maize and soybean. Agricultural and Forest Meteorology, 125, 105–115. Moffat, A. M., Papale, D., Reichstein, M., Hollinger, D. Y., Richardson, A. D., Barr, A. G., et al. (2007). Comprehensive comparison of gap-filling techniques for eddy covariance net carbon fluxes. Agricultural and Forest Meteorology, 147, 209–232. Monteith (1964). Evaporation and environment. Cambridge: Cambridge University Press. Montheith, J., & Unsworth, M. H. (2007). Principles of environmental physics (3rd ed.). : Elsevier. Mu, Q., Heinsch, F. A., Zhao, M., & Running, S. W. (2007). Development of a global evapotranspiration algorithm based on MODIS and global meteorology data. Remote Sensing of Environment, 111, 519–536. Mu, Q., Zhao, M., & Running, S. W. (2011). Improvements to a MODIS global terrestrial evapotranspiration algorithm. Remote Sensing of Environment, 115, 1781–1800. Mueller, B., Seneviratne, S. I., Jimenez, C., Corti, T., Hirschi, M., Balsamo, G., et al. (2011). Evaluation of global observations-based evapotranspiration datasets and IPCC AR4 simulations. Geophysical Research Letters, 38. Nagler, P. L., Cleverly, J., Glenn, E., Lampkin, D., Huete, A., & Wan, Z. M. (2005). Predicting riparian evapotranspiration from MODIS vegetation indices and meteorological data. Remote Sensing of Environment, 94, 17–30. Nagler, P. L., Glenn, E. P., Kim, H., Emmerich, W., Scott, R. L., Huxman, T. E., et al. (2007). Relationship between evapotranspiration and precipitation pulses in a semiarid rangeland estimated by moisture flux towers and MODIS vegetation indices. Journal of Arid Environments, 70, 443–462. Nagler, P., Morino, K., Murray, R. S., Osterberg, J., & Glenn, E. (2009). An empirical algorithm for estimating agricultural and riparian evapotranspiration using MODIS enhanced vegetation index and ground measurements of ET. I. Description of method. Remote Sensing, 1, 1273–1297. Nishida, K., Nemani, R. R., Glassy, J. M., & Running, S. W. (2003). Development of an evapotranspiration index from aqua/MODIS for monitoring surface moisture status. IEEE Transactions on Geoscience and Remote Sensing, 41, 493–501. Papale, D., Reichstein, M., Aubinet, M., Canfora, E., Bernhofer, C., Kutsch, W., et al. (2006). Towards a standardized processing of net ecosystem exchange measured with eddy covariance technique: Algorithms and uncertainty estimation. Biogeosciences, 3, 571–583. Papale, D., & Valentini, A. (2003). A new assessment of European forests carbon exchanges by eddy fluxes and artificial neural network spatialization. Global Change Biology, 9, 525–535. Reichstein, M., Falge, E., Baldocchi, D., Papale, D., Aubinet, M., Berbigier, P., et al. (2005). On the separation of net ecosystem exchange into assimilation and ecosystem respiration: Review and improved algorithm. Global Change Biology, 11, 1424–1439.

261

Rey, A., Pegoraro, E., Tedeschi, V., De Parri, I., Jarvis, P. G., & Valentini, R. (2002). Annual variation in soil respiration and its components in a coppice oak forest in Central Italy. Global Change Biology, 8, 851–866. Rouse, J. W., Haas, R. W., Schell, J. A., Deering, D. H., & Harlan, J. C. (1974). Monitoring the vernal advancement and retrogradation (Greenwave effect) of natural vegetation. Greenbelt, MD, USA: NASA/GSFC. Ryu, Y., Baldocchi, D. D., Black, T. A., Detto, M., Law, B. E., Leuning, R., et al. (2012). On the temporal upscaling of evapotranspiration from instantaneous remote sensing measurements to 8-day mean daily-sums. Agricultural and Forest Meteorology, 152, 212–222. Schmid, H. P., Grimmond, C. S. B., Cropley, F., Offerle, B., & Su, H. -B. (2000). Measurements of CO2 and energy fluxes over a mixed hardwood forest in the mid-western United States. Agricultural and Forest Meteorology, 103, 357–374. Shabanov, N. V., Dong, H., Wenze, Y., Tan, B., Knyazikhin, Y., Myneni, R. B., et al. (2005). Analysis and optimization of the MODIS leaf area index algorithm retrievals over broadleaf forests. IEEE Transactions on Geoscience and Remote Sensing, 43, 1855–1865. Sjöström, M., Ardö, J., Arneth, A., Boulain, N., Cappelaere, B., Eklundh, L., et al. (2011). Exploring the potential of MODIS EVI for modeling gross primary production across African ecosystems. Remote Sensing of Environment, 115, 1081–1089. Strahler, A. H., Muller, J. P., & Sciences Team Members, Modis (1999). MODIS BRDF albedo product: Algorithm theoretical basis document version 5.0. (53 pp.). Suni, T., Rinne, J., Reissell, A., Altimir, N., Keronen, P., Rannik, U., et al. (2003). Long-term measurements of surface fluxes above a Scots pine forest in Hyytiala, southern Finland, 1996–2001. Boreal Environment Research, 8, 287–301. Tang, R., Li, Z. -L., & Tang, B. (2010). An application of the Ts–VI triangle method with enhanced edges determination for evapotranspiration estimation from MODIS data in arid and semi-arid regions: Implementation and validation. Remote Sensing of Environment, 114, 540–551. Tedeschi, V., Rey, A. N. A., Manca, G., Valentini, R., Jarvis, P. G., & Borghetti, M. (2006). Soil respiration in a Mediterranean oak forest at different developmental stages after coppicing. Global Change Biology, 12, 110–121. Urbanski, S., Barford, C., Wofsy, S., Kucharik, C., Pyle, E., Budney, J., et al. (2007). Factors controlling CO2 exchange on timescales from hourly to decadal at Harvard Forest. Journal of Geophysical Research, 112, G02020. Van Dijk, A. I. J. M. (2010). The Australian water resources assessment system. Technical report 3. Landscape model (version 0.5) technical description WIRADA technical report. CSIRO Water for a Healthy Country Flagship (Canberra). van Dijk, A. I. J. M., & Renzullo, L. J. (2011). Water resource monitoring systems and the role of satellite observations. Hydrology and Earth System Sciences Discussions, 15, 39–55. Verma, S. B., Dobermann, A., Cassman, K. G., Walters, D. T., Knops, J. M., Arkebauer, T. J., et al. (2005). Annual carbon dioxide exchange in irrigated and rainfed maize-based agroecosystems. Agricultural and Forest Meteorology, 131, 77–96. Wang, K. C., Li, Z. Q., & Cribb, M. (2006). Estimation of evaporative fraction from a combination of day and night land surface temperatures and NDVI: A new method to determine the Priestley–Taylor parameter. Remote Sensing of Environment, 102, 293–305. Willmott, C. J. (1982). Some comments on the evaluation of model performance. Bulletin of the American Meteorological Society, 63, 1309–1313. Yang, F. H., White, M. A., Michaelis, A. R., Ichii, K., Hashimoto, H., Votava, P., et al. (2006). Prediction of continental-scale evapotranspiration by combining MODIS and AmeriFlux data through support vector machine. IEEE Transactions on Geoscience and Remote Sensing, 44, 3452–3461. Yebra, M., & Chuvieco, E. (2009). Linking ecological information and radiative transfer models to estimate fuel moisture content in the Mediterranean region of Spain: Solving the ill-posed inverse problem. Remote Sensing of Environment, 113, 2403–2411. Yebra, M., Chuvieco, E., & Riano, D. (2008). Estimation of live fuel moisture content from MODIS images for fire risk assessment. Agricultural and Forest Meteorology, 148, 523–536. Yuan, W. P., Liu, S. G., Yu, G. R., Bonnefond, J. M., Chen, J. Q., Davis, K., et al. (2010). Global estimates of evapotranspiration and gross primary production based on MODIS and global meteorology data. Remote Sensing of Environment, 114, 1416–1431. Zhang, Y. Q., Chiew, F. H. S., Zhang, L., Leuning, R., & Cleugh, H. A. (2008). Estimating catchment evaporation and runoff using MODIS leaf area index and the Penman– Monteith equation. Water Resources Research, 44, W10420. Zhang, K., Kimball, J. S., Mu, Q., Jones, L. A., Goetz, S. J., & Running, S. W. (2009). Satellite based analysis of northern ET trends and associated changes in the regional water balance from 1983 to 2005. Journal of Hydrology, 379, 92–110. Zhang, Y. Q., Leuning, R., Hutley, L. B., Beringer, J., McHugh, I., & Walker, J. P. (2010). Using long-term water balances to parameterize surface conductances and calculate evaporation at 0.05 degrees spatial resolution. Water Resources Research, 46.