Estimation of Vegetation Latent Heat Flux over Three Forest Sites in ChinaFLUX using Satellite Microwave Vegetation Water Content Index

Latent heat flux (LE) and the corresponding water vapor lost from the Earth’s surface to the atmosphere, which is called Evapotranspiration (ET), is one of the key processes in the water cycle and energy balance of the global climate system. Satellite remote sensing is the only feasible technique to estimate LE over a large-scale region. While most of the previous satellite LE methods are based on the optical vegetation index (VI), here we propose a microwave-VI (EDVI) based LE algorithm which can work for both day and night time, and under clear or non-raining conditions. This algorithm is totally driven by multiple-sensor satellite products of vegetation water content index, solar radiation, and cloud properties, with some aid from a reanalysis dataset. The satellite inputs and the performance of this algorithm are validated with in situ measurements at three ChinaFLUX forest sites. Our results show that the selected satellite observations can indeed serve as the inputs for the purpose of estimating ET. The instantaneous estimations of LE (LEcal) from this algorithm show strong positive temporal correlations with the in situ measured LE (LEobs) with the correlation coefficients (R) of 0.56–0.88 in the study years. The mean bias is kept within 16.0% (23.0 W/m2) across the three sites. At the monthly scale, the correlations between the retrieval and the in situ measurements are further improved to an R of 0.84–0.95 and the bias is less than 14.3%. The validation results also indicate that EDVI-based LE method can produce stable LEcal under different cloudy skies with good accuracy. Being independent of any in situ measurements as inputs, this algorithm shows great potential for estimating ET under both clear and cloudy skies on a global scale for climate study.


Introduction
The process of latent heat flux (LE) and evapotranspiration (ET) is associated with the exchange of water and energy between the land surface and the atmosphere [1,2].As a major component of the land surface energy budget, LE consumes about three-fifths of net radiation income, ranging from 48% to 88% based on different models [3,4].Over 80% of terrestrial ET is from plant transpiration [5].Globally, the ET from forests account for about 45% of the total ET [6].Therefore, the accurate and real-time estimation of forest LE is critical to improve our understanding of water resource management, the hydrological cycle and the associated energy balance between the terrestrial ecosystem and the atmosphere.
The conventional ground-based flux tower methods (Bowen ratio and eddy covariance) can provide relatively accurate estimates of LE within a small footprint area (200-500 m 2 radius depending on the height of tower) around flux towers at local scales with an uncertainty of about 10%-30% [5].However, the spatial representativeness of point LE estimates is limited to the regional scale.Large uncertainties will appear if point measurements are extended to regional scales because of the heterogeneity of land surfaces.
Satellite remote sensing has become the only feasible technique for scaling the point LE estimates from flux towers to mean LE over a large-scale area.Most of the current remote sensing approaches to estimate LE are based on optical vegetation indices (VIs) such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) [7].Researchers utilize these optical VIs in order to estimate foliage density and to predict ET over a vegetated surface [8].A number of satellite optical VI methods have been successfully proposed to estimate ET for multiple landscape types from the site to regional scale [9][10][11][12][13][14][15][16].Among these methods, the resistance-based Penman-Monteith (PM) model is considered as the preferred method due to its precise physical mechanism.VIs are typically used for the estimation of canopy resistance (or conductance) in the PM model [8].Previous studies have shown that optical VI-based resistance and LE estimation show good performances at the 8-day or 16-day time scale and under clear or less cloudy sky conditions [11,[17][18][19].However, since the optical VIs are retrieved based on the visible, near-infrared and shortwave-infrared bands [19] which are very sensitive to clouds and aerosols, their applicability under cloudy and overcast skies is limited.
When compared with land surface reflection at visible and infrared channels, the microwave radiative properties of land surface are less affected by atmospheric conditions such as cloud and aerosol, and the passive microwave signal is independent of solar radiation and is thus available at both the daytime and the nighttime.Microwave-based VIs thus have the advantage of being used under all sky conditions.Several microwave VIs have been proposed for monitoring vegetation, such as the frequency index (FI) [20], microwave polarization difference index (MPDI) [21] and microwave vegetation index (MVI) [22].Most of those microwave vegetation indexes are derived under clear sky only and are designed for short vegetation.It is currently still a challenge to use these microwave VIs to quantify LE with a physical scheme due to their complexities.Barraza et al. [23,24] used the combined empirical method of microwave and optical VIs to estimate conductance and LE over forest and savanna ecosystems.They illustrated the superiority of such combinations from multiple sensors for LE estimation.However, such an empirical method between canopy conductance and VIs [17][18][19]23,24] is insufficient to reflect the fast dynamics of conductance and LE affected by local environments during short-term periods (e.g., diurnal and daily scale).Furthermore, most of the validations of estimated LE are conducted at the 8-day (or 16-day) time scale and under clear or less cloudy sky since optical VIs are easily contaminated by cloud cover.In practice, the LE under cloudy sky differs significantly from that under clear sky.Thus, LE estimation under cloudy and overcast sky is equally important, particularly for the study of short-term interaction of vegetation-cloud.However, few studies have conducted ET or LE estimations under such sky conditions [25].
Min and Lin [26,27] proposed a new satellite remote-sensed microwave emissivity difference vegetation index (EDVI) at a mid-latitude forest (Harvard forest).The EDVI is defined as (MLSE 19 -MLSE 37 )/(MLSE 19 +MLSE 37 ), where MLSE 19 and MLSE 37 are the microwave land surface emissivities at 19 and 37 GHz.EDVI was found closely related to vegetation water content based on in-situ measured leave amount and their analysis on microwave radiative transfer [26,27].In addition, a stronger relationship between the EDVI and evaporative fraction (EF) than that between NDVI and EF was found [26,27].EDVI retrievals became available at the regional scale over a long time period [28].Based on the above work, a quantitative algorithm was developed on the basis of satellite EDVI and local measurements in order to estimate EF and ET fluxes at the Harvard forest [29].By using observations from multiple satellites, they even observed the diurnal variations of ET at the site as well.However, in their study, several in situ measurements were required as inputs, namely: net radiation, air temperature, photosynthetically active radiation, etc [29].Consequently, the method of Li et al. [29] was limited to sites where the above inputs are measured.
For most applications in modeling and climate study, ET and LE estimation at large scale and over a long term is required [30][31][32].A satellite-based ET retrieval algorithm independent of in situ measurements is thus required.
With the rapid progression of remote sensing technology, an increasing number of geophysical parameters are available from satellites.It is promising to develop an LE algorithm independent of in situ measurements that is driven only by satellite observations and reanalysis datasets.In this study, we tested this hypothesis and applicability of the LE method over dense vegetation in combination with the EDVI from advanced microwave scanning radiometer for EOS (AMSR-E), the net radiation flux from the clouds and earth's radiation energy system (CERES), the vegetation fraction information from moderate-resolution imaging spectroradiometer (MODIS) and the associated meteorology parameters from the reanalysis dataset of European Centre for Medium-Range Weather Forecasts (ECMWF).The results of this method were validated against the in situ measurements at selected ChinaFLUX forests sites.
In ChinaFLUX, part of a "network of regional networks" (FLUXNET), the eddy covariance (EC) technique was used to measure the H 2 O, CO 2 , and heat fluxes between the atmosphere and the ecosystems in China [33].Because the concept of EDVI was originally developed in forest, we selected three typical forest ecosystems for the validation study.These sites are the Dinghushan (DHS) subtropical evergreen broad-leaved forest, Qianyanzhou (QYZ) subtropical plantation forest, and Changbaishan (CBS) temperate deciduous mixed forest.More information is provided in Section 2.1.

Site Descriptions
The three selected forest towers in this study are shown in Figure 1.Three years of in-situ measurements (2003)(2004)(2005) at these sites are available in this study.And at each site, we split the measurements into different time periods for calibration and validation study.These sites are in different latitude zones and cover a wide range of temperature and precipitation.The annual mean precipitation is 1956 mm, 1485 mm and 695 mm, while the annual mean temperature is 21 • C, 17.9 • C and 3.6 • C at DHS, QYZ and CBS, respectively.The terrain is also complex among the forest sites: DHS is on a steep 30 • slope, QYZ is on slightly choppy terrain, and CBS is on flat terrain [34].Previous studies have indicated that forest sites of ChinaFLUX network achieve 57%-73% energy balance closure during the daytime [33,35].More information about these sites are provided in Table 1.ChinaFLUX data in this study is available from the website: http://159.226.110.139/pingtai/LoginRe/opendata.jsp.The 30-minute in situ measurements of LE (LEobs), solar radiation, photosynthetic active radiation (PAR), air temperature (Ta), and precipitation from 2003 to 2005 are used for validation in this study.To compare with satellite observations and retrievals obtained around 13:30 local time, the averaged in-situ measurements from 12:30-14:30 are used.Temperature is a key factor affecting vegetation activity at higher latitudes.Low Ta will induce weak metabolism and suppress the LE of vegetation.Therefore, the days with midday Ta lower than 0 degrees are excluded from the validation study at CBS.This makes no impact on the samples in DHS and QYZ.

Satellite inputs of EDVI-based LE method
Table 2 provides the basic information of the employed satellite products and reanalysis datasets and variables in this study.EDVI plays a fundamental role in the current LE algorithm as it determines the spatial and temporal variations of vegetation resistance to evapotranspiration process.The satellite EDVI retrievals during 2003-2005 over China were derived from a similar method to Min et al. [28] using multiple channel microwave measurements from the AMSR-E on the Aqua satellite with a spatial resolution of ~20 km at local time 13:30 each day.It should note that the instantaneous retrievals of EDVI are conducted under no rain conditions which are identified by AMSR-E rain rate/type product [28].Since snow on the ground could significantly affect the value of  Carbon and water fluxes of the three forests are measured using open-path eddy covariance (EC) system.The system consists of an open-path infrared gas analyser (Li-7500, Licor Inc., Lincoln, NB, USA) and a 3-D sonic anemometer (CSAT3, Campbell Scientific Ltd., USA).The signals of fluxes with 10Hz sampling frequency are recorded by system and 30-min averaged fluxes are derived from black averaging within 30-min step.Other meteorological variables are measured simultaneously with 30-min temporal resolution.Solar radiation is measured using radiometers (Model CM11 and Model CNR-1, Kipp & Zonnen, Delft, Netherlands).Photosynthetic active radiation (PAR) is measured using a quantum sensor (Model LI190SB, LICOR Inc.), air temperature (Ta) is measured using shielded and aspirated probes (Model HMP45C, Campbell Scientific Inc.).Precipitation is recorded using a rain gauge above the canopy (Model 52203, Rm Young, Traverse City, MI, USA).More detailed information can be found in previous studies [33,34,36] and the references therein.
The 30-min in situ measurements of LE (LE obs ), solar radiation, photosynthetic active radiation (PAR), air temperature (Ta), and precipitation from 2003 to 2005 are used for validation in this study.To compare with satellite observations and retrievals obtained around 13:30 local time, the averaged in-situ measurements from 12:30-14:30 are used.Temperature is a key factor affecting vegetation activity at higher latitudes.Low Ta will induce weak metabolism and suppress the LE of vegetation.
Therefore, the days with midday Ta lower than 0 degrees are excluded from the validation study at CBS.This makes no impact on the samples in DHS and QYZ.

Satellite Inputs of EDVI-based LE Method
Table 2 provides the basic information of the employed satellite products and reanalysis datasets and variables in this study.EDVI plays a fundamental role in the current LE algorithm as it determines the spatial and temporal variations of vegetation resistance to evapotranspiration process.The satellite EDVI retrievals during 2003-2005 over China were derived from a similar method to Min et al. [28] using multiple channel microwave measurements from the AMSR-E on the Aqua satellite with a spatial resolution of ~20 km at local time 13:30 each day.It should note that the instantaneous retrievals of EDVI are conducted under no rain conditions which are identified by AMSR-E rain rate/type product [28].Since snow on the ground could significantly affect the value of EDVI, at the CBS site we only retrieve LE during the growing season from April to October, based on the study of Li et al. [37].The product of single scanner footprint toa/surface fluxes and clouds (SSF) of CERES provides the estimates of instantaneous daytime net and downward shortwave (nSW and dSW) and longwave (nLW and dLW) radiation fluxes at surface at ~20 km resolution (Table 2) [38].Those fluxes are estimated based on the NASA Langley parameterized shortwave/longwave algorithm (LPSA/LPLA) methods [39,40].The CERES SSF datasets are available at https://eosweb.larc.nasa.gov/project/ceres/ssf_aqua-fm3_ed4a_table. NDVI, derived from the observations of MODIS on the Aqua satellite, is also an optical VI which is highly correlated to the green foliage of vegetation.We use NDVI to estimate the vegetation fractional coverage.To do this, 16-day MODIS NDVI product (MYD13C1) with 0.05 • resolution was used.To estimate the temporal variation of vegetation fraction, the 16-day NDVI was further linearly interpolated to achieve daily NDVI.MYD13C1 is available from the Land Processes Distributed Active Archive Center (LPDAAC; https://lpdaac.usgs.gov/).
ERA-20C reanalysis dataset of ECMWF provides estimations of air temperature at 2 m (t2m) and wind at 10 and 100 m (U 10m , U 10m ) required by the ET algorithm.ERA-20C has a temporal resolution of 3 h and 0.125 • × 0.125 • spatial resolution (Table 2), and the values at 6:00 UTC (14:00 local time in China) were used in this study.The dataset is available on the ECMWF public datasets website (http://apps.ecmwf.int/datasets/).

Description of EDVI-based LE Algorithm
The primary driving force of LE is the available energy at the surface: the received net radiation Rn minus the energy transport to ground (G).An evaporative fraction (EF) is used to define the ratio of ET to the total available energy (Rn-G).This idea is often used to estimate instantaneous ET at the time of a satellite overpass [41][42][43][44][45].
In forest ecosystems, transpiration generally dominates evapotranspiration and accounts for 80-90% of the ET amount [46,47].We therefore scale the available energy in given satellite footprint to vegetation part (with the subscript veg) by the vegetation fraction coverage (VFC) and only take into account the ET from vegetation using the following formula.

Estimation of Available Energy for Vegetation
Radiation fluxes under all sky conditions are required in our algorithm.CERES SSF directly provides the all-sky surface net shortwave radiation (nSW) and surface long-wave radiation (nLW) (Table 2).Therefore, we use the sum of nSW and nLW as the input of Rn in our ET algorithm.Validations of CERES surface radiation flux with in-situ measurements are conducted in results.
Surface ground heat flux can be estimated in combination with vegetation cover and net radiation.We use the model of Su [48] to calculate G, such that: where f veg (set as 0.05) and f soil (set as 0.315) are the ratios of G to Rn over full vegetation and bare soil, respectively [48].VFC is the vegetation fractional coverage which can be obtained based on NDVI according to the method of Gutman and Ignatov [49]: where NDVI ∞ and NDVI 0 are the NDVI of full vegetation (i.e., when VFC = 1) and bare soil (i.e., when VFC = 0), respectively.The values of NDVI ∞ depend on different vegetation types, but they are relatively stable over forest types [50].Values of NDVI 0 have very small variations when VFC = 0 [51].Following Li and Zhang [50] and Zeng et al. [51], we adopt 0.90 and 0.1 as the values of NDVI ∞ and NDVI 0 , respectively.When the calculated VFC is less than 0, VFC is set to 0, and when the calculated VFC is larger than 1, VFC is set to 1.In this equation, the interpolated daily NDVI from 16-day values is used to calculate daily VFC under the assumption that vegetation foliage or greenness changes slowly and linearly within each 16-day period.We provide the variations of 16-day NDVI in Figure 2 in Section 3.1.
Remote Sens. 2019, 11, x FOR PEER REVIEW 8 of 21 that EDVI can be potentially used to retrieve LE over the three forests.As an indicator of vegetation foliage, 16-day NDVI are found to have the small and relatively slow changes within each 16-day interval, suggesting that the use of interpolated daily NDVI from 16-day NDVI for the calculation of VFC may introduce small errors.

Validation of Satellite Radiation Flux and Air Temperature Inputs
Studies have been conducted to validate the CERES surface fluxes over different surface types (i.e., island, coastal, polar, continental, and desert) and have found good accuracies of the data [38,60].In this study, since there is no in-situ measurements of net radiation (Rn) available, we therefore conduct the validation of downward shortwave radiation fluxes (dSW) between satellite observations and in-situ measurements instead.
Figure 3 shows the comparison of all-sky CERES dSW and in situ measured downward solar radiation flux at midday from 2003 to 2005.Over all forest sites, CERES dSW agrees well with the insitu measured solar radiation, with high positive correlation coefficients (R) of 0.81, 0.95, and 0.90 at DHS, QYZ, and CBS, respectively.The relative biases of all sites are less than 9%.Because dSW is the dominant term of Rn, we therefore conclude that CERES SSF estimated Rn is suitable for estimating the real net radiation for calculating ET.

Estimation of Evaporative Fraction of Forest
The estimation of EF veg is from the same method of Nishida et al. [52] and Li et al. [29]: where α is the Priestley-Taylor's parameter (1.26), ∆ is the slope of saturated vapor pressure at air temperature (hPa/K), γ is the psychometric constant (Pa/K) (66.5 Pa/K), r a is the aerodynamic resistance (s/m), and r c is the canopy resistance (s/m) which is highly related to EDVI in this study.We calculate ∆ based on the formula in study of Murray [53]: where Ta is the air temperature at 2 m (t2m) from ERA-20C.The aerodynamic resistance r a is determined by wind.Kondo [54,55] proposed two empirical formulas for forest canopy and grasslands.The formulas have been used for the estimation of ET in other studies [52,56].We thus follow their study and use equation (6) to calculate the r a of forests: grass and crops (6) where U 50m is the wind speed at 50 m (ms −1 ) which is calculated by averaging the mean values of 10-m wind (U 10m ) and 100-m wind (U 100m ) from ERA-20C.
The canopy resistance r c can be parameterized by the use of Jarvis-type equation.
where r c has two components: the stomatal resistance (r stomatal ) and the cuticle resistance (r cuticle ).r cuticle was set to be 10 5 ms −1 based on other studies [11].
The stomatal resistance r stomatal consists of five stress functions and the minimum canopy resistance (r cmin ).The five stress functions are used to quantify the impacts on r c imposed by environmental conditions.The functions of air temperature (f 1 (Ta)) and photosynthetically active radiation (f 2 (PAR)) were adopted from Jarvis [57] and Kosugi [58], respectively.The day-to-day variation of EDVI (dEDVI) represents the response of canopy leaves to environmental factors such as water vapor pressure deficit (VPD), water potential (Ψ) and ambient carbon dioxide concentration (CO 2 ).In this study, we parameterize their stress functions (i.e., f 3 f 4 f 5 ) as a whole by using their fair correlation with dEDVI based on Li et al. [29]: where a and b are the coefficients.It is worth noting that the accurate determination of a and b require the field measurements of stomatal and canopy resistance.Li et al. [29] used the field measurements of resistances and developed the relationship at forest.Since such field measurements are unavailable in this study, we thus follow their study and adopt 1.186 and 105.755 as the values of a and b for these forests.
To describe the mean seasonal variation of EDVI associated with the different stages in the growing season of the forest, we define the normalized EDVI ( N EDVI) as (EDVI-EDVI min )/(EDVI max -EDVI min ), where EDVI min and EDVI max are respectively the minimum and maximum EDVI values during growing seasons at each site [27,29].The seasonal variation of EDVI (i.e., N EDVI) represents the slow change of VWC which is highly correlated to the vegetation resistance of LE.N EDVI is thus used to determine r cmin which is the key parameter for the estimation of r stomatal .Li et al. [29] found that there was a general anti-correlation between N EDVI and canopy resistance.The minimum canopy resistance can be parameterized as: The canopy resistance will decrease when the plant leaves develop.In the case of Li et al. [29], the measurements of canopy resistance at Harvard forest were available to determine the r cmin0 .Based on the field studies of Kelliher et al. [59], the maximum vegetation conductance ranges from 13.0 mm/s (tropical rainforest) to 32.5 mm/s (cereals), while the corresponding r cmin (G −1 smax ) ranges from 76.9 s/m to 30.8 s/m.On the basis of this, r cmin0 was set to be a reference value, 50 s/m, for all three forest types in our study.

Results
Since our objective of this study is to improve the EDVI-based LE method driven by satellite and reanalysis data under all sky conditions, the inputs of the method are crucial and were first investigated in this section.Section 3.1 shows the consistency between LE obs and two vegetation indices.Section 3.2 provides the validation of accuracy of satellite and reanalysis data.The estimated instantaneous LE (LE cal ) under all skies are validate in Section 3.3.Then the accuracy of LE cal under different cloudy sky conditions are further investigated in Section 3.4.

Time Series of EDVI, NDVI and In-Situ Measured LE
The time series of midday in situ measured LE obs , satellite-based EDVI and NDVI at the three forest sites in this study are shown in Figure 2. The EDVI and NDVI show the consistent seasonal cycles with the LE obs at all forests (maximum in summer and minimum in winter).EDVI vary from -0.015 to 0.02 during growing seasons when NDVI vary from 0.3 to 0.9.The LE obs of forest ecosystems can have significant variations during the short-term period due to environmental and weather effects (e.g., precipitation and solar radiation).Daily EDVI can indicate such fast changes of vegetation status and show the similar day-to-day variations to the LE obs at all sites (Figure 2), which suggests that EDVI can be potentially used to retrieve LE over the three forests.As an indicator of vegetation foliage, 16-day NDVI are found to have the small and relatively slow changes within each 16-day interval, suggesting that the use of interpolated daily NDVI from 16-day NDVI for the calculation of VFC may introduce small errors.

Validation of Satellite Radiation Flux and Air Temperature Inputs
Studies have been conducted to validate the CERES surface fluxes over different surface types (i.e., island, coastal, polar, continental, and desert) and have found good accuracies of the data [38,60].In this study, since there is no in-situ measurements of net radiation (Rn) available, we therefore conduct the validation of downward shortwave radiation fluxes (dSW) between satellite observations and in-situ measurements instead.
Figure 3 shows the comparison of all-sky CERES dSW and in situ measured downward solar radiation flux at midday from 2003 to 2005.Over all forest sites, CERES dSW agrees well with the in-situ measured solar radiation, with high positive correlation coefficients (R) of 0.81, 0.95, and 0.90 at DHS, QYZ, and CBS, respectively.The relative biases of all sites are less than 9%.Because dSW is the dominant term of Rn, we therefore conclude that CERES SSF estimated Rn is suitable for estimating the real net radiation for calculating ET.Photosynthetically active radiation (PAR) is required in the estimation of EDVI-based canopy resistance.Studies have reported that PAR is linearly related to incident solar radiation [61][62][63].In our study, all-sky CERES dSW is used to estimate PAR by multiplying it by a constant of 1.70 μmol/W.The results in Figure 4 indicate these estimations agree well with in-situ observations across three sites.Correlation coefficients strongly vary from 0.81 to 0.95 with a small bias of no more than 11%.Photosynthetically active radiation (PAR) is required in the estimation of EDVI-based canopy resistance.Studies have reported that PAR is linearly related to incident solar radiation [61][62][63].In our study, all-sky CERES dSW is used to estimate PAR by multiplying it by a constant of 1.70 µmol/W.The results in Figure 4 indicate these estimations agree well with in-situ observations across three sites.Correlation coefficients strongly vary from 0.81 to 0.95 with a small bias of no more than 11%.Photosynthetically active radiation (PAR) is required in the estimation of EDVI-based canopy resistance.Studies have reported that PAR is linearly related to incident solar radiation [61][62][63].In our study, all-sky CERES dSW is used to estimate PAR by multiplying it by a constant of 1.70 μmol/W.The results in Figure 4 indicate these estimations agree well with in-situ observations across three sites.Correlation coefficients strongly vary from 0.81 to 0.95 with a small bias of no more than 11%.Air temperature is an important parameter for estimating forest ET.We compare the 2-meter air temperature (t2m) from ERA-20C (t2m-ECMWF) and from NCEP FNL (t2m-NCEP) with in-situ measurements.As shown in Figure 5, both t2m-ECMWF and t2m-NCEP have a high correlation with in-situ measured temperature.However, t2m-NCEP estimations are systematically lower by 3-5 K than in-situ measurements.In contrast, t2m-ECMWF are in better agreement with in situ measurements with a smaller bias (0.1-2.7K).The overall biases of t2m-ECMWF for all sites are within 3K.Therefore, we determine to use ERA-20C rather than NCEP FNL to provide complementary parameters to the satellite observations.Air temperature is an important parameter for estimating forest ET.We compare the 2-m air temperature (t2m) from ERA-20C (t2m-ECMWF) and from NCEP FNL (t2m-NCEP) with in-situ measurements.As shown in Figure 5, both t2m-ECMWF and t2m-NCEP have a high correlation with in-situ measured temperature.However, t2m-NCEP estimations are systematically lower by 3-5 K than in-situ measurements.In contrast, t2m-ECMWF are in better agreement with in situ measurements with a smaller bias (0.1-2.7 K).The overall biases of t2m-ECMWF for all sites are within 3K.Therefore, we determine to use ERA-20C rather than NCEP FNL to provide complementary parameters to the satellite observations.

The EDVI-Based LE Estimation
The time series of EDVI-based LEcal, in situ LEobs, and their differences at the three forest sites are shown in Figure 6.Statistic results are shown in Table 3.Generally, LEcal has the capability to capture the seasonality of forest LEobs from 2003 to 2005 correctly.Both of them reach maximums in the mature

The EDVI-Based LE Estimation
The time series of EDVI-based LE cal , in situ LE obs , and their differences at the three forest sites are shown in Figure 6.Statistic results are shown in Table 3.Generally, LE cal has the capability to capture the seasonality of forest LE obs from 2003 to 2005 correctly.Both of them reach maximums in the mature stage of growing seasons in mid-summer due to the ample water and solar radiation for evapotranspiration (Figure 6).There is a strong correlation between LE cal and LE obs with the overall R being 0.56-0.88.In terms of magnitude, LE cal, ranging from 0 to 500 Wm −2 at DHS and QYZ in southern China, and from 0 to 400 Wm −2 at CBS in northeastern China, matches LE obs well.Large discrepancies occur at QYZ in 2005 due to the serious imbalance of energy, which leads to a significant underestimation of in situ LE obs [64,65] and should be responsible for the large bias in this year (Table 3).Because of this, we exclude the result of 2005 in QYZ in the following discussion.The samples with the differences between LEcal and LEobs less than 150 Wm -2 account for 94%, 93%, and 92% of total samples at DHS, QYZ, and CBS, respectively (Figure 6).Our algorithm tends to underestimate the LE during transient periods with about 30-100 Wm -2 (particularly in non-snowy wintertime and early spring).In spite of this, results in Table 3 show that the method can produce the small bias varying from -30.7 to 32.8 Wm -2 with the RMSE from 56.0 to 90.9 Wm -2 , respectively.The relative bias at three forests was kept within 23% for most of the years and range from -18.6% to 22.8%.The regression lines between LEcal and LEobs are well close to a 1:1 line at all sites with slopes of 0.70-1.29 (Figure 7, Table 3).The samples with the differences between LE cal and LE obs less than 150 Wm −2 account for 94%, 93%, and 92% of total samples at DHS, QYZ, and CBS, respectively (Figure 6).Our algorithm tends to underestimate the LE during transient periods with about 30-100 Wm −2 (particularly in non-snowy wintertime and early spring).In spite of this, results in Table 3 show that the method can produce the small bias varying from −30.7 to 32.8 Wm −2 with the RMSE from 56.0 to 90.9 Wm −2 , respectively.The relative bias at three forests was kept within 23% for most of the years and range from −18.6% to 22.8%.The regression lines between LE cal and LE obs are well close to a 1:1 line at all sites with slopes of 0.70-1.29 (Figure 7, Table 3).The algorithm performs better at QYZ and CBS in terms of R (0.80-0.88) and RMSE (56-83.5 Wm - 2 ).CBS has the smallest mean RMSE (67.9 Wm -2 ), highest R (0.83) and lowest relative bias (-9.6%) for all study years (Table 3).The best LE estimations at DHS, QYZ, and CBS occur in 2005 (R=0.73,bias= 17.4 Wm -2 , RMSE=73.1Wm -2 ), 2003 (R=0.82,bias=3.0Wm - , RMSE=74.2Wm -2 ), and 2004 (R=0.88,bias= -22.7 Wm -2 , RMSE=56.0Wm -2 ), respectively.However, the estimation performed relatively poorly at DHS in 2003 (R=0.58,bias= 32.8 Wm -2 , RMSE=90.9Wm -2 ) (Figure 7, Table 3).Our algorithm could be affected by heavily rainy events occurring before or after when Aqua satellite overpasses.For these rainy days, our retrieval algorithm is expected to underestimate the LE because of the reduced microwave EDVI over wet surfaces and the omission of interception evaporation.Some rainy samples severely contaminated by such precipitation are marked by solid circles in Figure 7.If we excluded these samples, the performances would be generally improved at all sites, particularly for R and RMSE.This comparison results are shown in Figure 8.The R would The algorithm performs better at QYZ and CBS in terms of R (0.80-0.88) and RMSE (56-83.5 Wm −2 ).CBS has the smallest mean RMSE (67.9 Wm −2 ), highest R (0.83) and lowest relative bias (−9.6%) for all study years (Table 3).3).
Our algorithm could be affected by heavily rainy events occurring before or after when Aqua satellite overpasses.For these rainy days, our retrieval algorithm is expected to underestimate the LE because of the reduced microwave EDVI over wet surfaces and the omission of interception evaporation.Some rainy samples severely contaminated by such precipitation are marked by solid circles in Figure 7.If we excluded these samples, the performances would be generally improved at all sites, particularly for R and RMSE.This comparison results are shown in Figure 8.The R would be improved to 0.66-0.91 for all years and the RMSE would be reduced to 48.Additionally, a few samples with much larger LEcal than LEobs occurred during growing seasons, such as 2003 at DHS (mainly in summer with less precipitation).A possible explanation is that under high temperature and less precipitation conditions, the plant water deficit will induce leaf stomatal closure to prevent excessive water deficits in the plants during summer time [66].As a result, the real plant physiological activities, such as leaf transpiration, carbon gain, and growth, are remarkably suppressed [66].
The comparison results between LEcal and LEobs at monthly scale show the better performance (Figure 9).R can reach 0.84, 0.88, and 0.95 for DHS, QYZ, and CBS with the bias of 14.3%, 0.3%, and -12.9 %, respectively.The standard deviations of monthly mean LEcal and LEobs are comparable.As discussed, the results can be improved after removing the heavily rain-contaminated days (Figure 9).Additionally, a few samples with much larger LE cal than LE obs occurred during growing seasons, such as 2003 at DHS (mainly in summer with less precipitation).A possible explanation is that under high temperature and less precipitation conditions, the plant water deficit will induce leaf stomatal closure to prevent excessive water deficits in the plants during summer time [66].As a result, the real plant physiological activities, such as leaf transpiration, carbon gain, and growth, are remarkably suppressed [66].
The comparison results between LE cal and LE obs at monthly scale show the better performance (Figure 9).R can reach 0.84, 0.88, and 0.95 for DHS, QYZ, and CBS with the bias of 14.3%, 0.3%, and −12.9%, respectively.The standard deviations of monthly mean LE cal and LE obs are comparable.As discussed, the results can be improved after removing the heavily rain-contaminated days (Figure 9).closure to prevent excessive water deficits in the plants during summer time [66].As a result, the real plant physiological activities, such as leaf transpiration, carbon gain, and growth, are remarkably suppressed [66].
The comparison results between LEcal and LEobs at monthly scale show the better performance (Figure 9).R can reach 0.84, 0.88, and 0.95 for DHS, QYZ, and CBS with the bias of 14.3%, 0.3%, and

Validation of EDVI-based LE under different cloudy sky
Since EDVI is able to indicate vegetation hydrological states under both clear and cloudy sky [26,27], EDVI-based LE method combined with all-sky satellite-retrieved radiation also has the capability of estimating LE under different cloud covers (Frc).Figure 10 shows the comparison of LEcal and LEobs under a partly cloudy sky (Frc<40%), cloudy sky (Frc >=40%) and all cloudy sky (0%<=Frc<=100%).Their corresponding statistical metrics are shown in Figure 11.In general, our

Validation of EDVI-based LE under Different Cloudy Sky
Since EDVI is able to indicate vegetation hydrological states under both clear and cloudy sky [26,27], EDVI-based LE method combined with all-sky satellite-retrieved radiation also has the capability of estimating LE under different cloud covers (Frc).Figure 10 shows the comparison of LE cal and LE obs under a partly cloudy sky (Frc < 40%), cloudy sky (Frc >= 40%) and all cloudy sky (0% <= Frc <= 100%).Their corresponding statistical metrics are shown in Figure 11.In general, our method has good performance at three sites.The slopes of fit lines are all close to 1.0 with the range of 0.88 to 1.16 (Figures 10 and 11), suggesting that the method can produce small systematic bias in the estimation of instantaneous forest LE under different cloud covers.The capability is also illustrated in the results of bias and relative bias in Figure 11.For all three sites, the method produces bias less than 35 W/m 2 (26%) at instantaneous scale when compared with in-situ measurements.These biases vary from −34.2 W/m 2 to 25.1 W/m 2 with the relative values from −25.9% to 16.7% (Figure 11), respectively.The mean bias (relative bias) under Frc < 40%, Frc ≥ 40% and all cloudy sky are well kept within 11 W/m 2 with the values of 10.6 W/m 2 (6.3%), 0.32 W/m 2 (0.8%) and 7.7 W/m 2 (5%), respectively.In addition, a good correlation between LE cal and LE obs under different cloudy skies are also found at three sites with the R of 0.62-0.80,which suggests that the seasonal dynamics of forest LE under cloud cover can be well recaptured by the EDVI-based method.Because of the relatively coarse spatial resolution (see the related discussion in Section 4), the method could produce relatively large RMSE (59 to 90 W/m 2 ) illustrated by the scattering of samples in Figure 10.
Most importantly, it can be found that the EDVI-based method is able to produce stable statistic metrics under different cloud cover conditions for three typical forests.These results indicate that the developed EDVI-based method in this study, completely driven by satellite and reanalysis datasets, can be used to estimate forest LE effectively from clear sky to cloudy sky.
LE under cloud cover can be well recaptured by the EDVI-based method.Because of the relatively coarse spatial resolution (see the related discussion in section 4), the method could produce relatively large RMSE (59 to 90 W/m 2 ) illustrated by the scattering of samples in Figure 10.
Most importantly, it can be found that the EDVI-based method is able to produce stable statistic metrics under different cloud cover conditions for three typical forests.These results indicate that the developed EDVI-based method in this study, completely driven by satellite and reanalysis datasets, can be used to estimate forest LE effectively from clear sky to cloudy sky.

Summary of the method
The above results suggest that the new EDVI-based LE algorithm can effectively estimate instantaneous LE during midday at three different forests.Compared with the early EDVI method of

Summary of the Method
The above results suggest that the new EDVI-based LE algorithm can effectively estimate instantaneous LE during midday at three different forests.Compared with the early EDVI method of Li et al. [29] based on in situ measurements as inputs, the developed method in this study was completely driven by all-sky satellite-retrieved radiation fluxes and reanalysis-based meteorological data.Therefore, this method has the potential capability to estimate LE in unfrequented places, such as area in the deep forests and high mountains.
Which is central to our method is using microwave EDVI to quantify the stomatal and canopy resistance under both clear sky and cloudy sky.Most of previous resistance-based LE methods are based on optic VIs [17][18][19] which cannot be effectively implemented under cloudy sky.Their validations are thus limited to under clear or partly cloudy sky.Barraza et al. [23,24] compared the estimated surface conductance based on the different regression models of microwave and optical VIs in forest and savanna ecosystems.Although they concluded that the combined microwave-optical VIs method produced the best conductance and LE estimation at the 8-day mean scale, the performances under different sky conditions were not investigated at finer temporal resolution.EDVI-based LE method is designed for the instantaneous all-sky LE.The validations under different cloudy skies show that our method can produce stable forest LE from clear to cloudy sky with good accuracy, which could be important for the study of EF and LE under cloudy sky.

Uncertainties in the Results
In spite of these, some uncertainties in our EDVI-based LE method should be noted.The retrieved EDVI (EDVI) over a landscape is the integration of EDVI from vegetation (EDVIveg) and bare soil (EDVIsoil).During growing seasons, EDVI is dominated by EDVIveg and can be used as a good indicator of vegetation states.Thus, the LE estimated from EDVI approximates the in-situ measured LE.However, for some vegetation regions such as deciduous forests, EDVI can be affected by soil signals to different extent during transient periods.
Since the EDVI-based LE method is implemented over forests in this study, the evaporation of intercepted water and bare soil are simply omitted in the forest LE estimation.This simplification could result in some underestimation, particularly during the non-growing seasons of deciduous forests when soil evaporation may dominate the total LE.Further, as we discussed in Section 3.3, the algorithm could be affected by heavily rainy events occurring near the time when Aqua satellite overpasses.The wet canopy surface can reduce EDVI and thus result in the underestimation of LE in our method.In spite of these, results in this study (Table 3) indicate that our current method for instantaneous LE estimation at forests with the overall mean bias of 16% (DHS), 11.4% (QYZ) and −9.6% (CBS) (Table 3) is well within the typical error range of 30% for satellite VI-based LE method [8].
As the surrogates for in situ measurements, the quality of the satellite remotely sensed and reanalyzed data, which are the inputs of our algorithm, particularly net radiation, which is the direct driving force for evapotranspiration, would highly affect the accuracy of the LE estimate, and the effect of this would vary at different sites (Figure 3).Yan et al. [60] compared CERES surface radiation fluxes with in situ measurements over Loess Plateau and found CERES surface downward radiation fluxes have higher accuracies in clear sky than those in cloudy sky.The standard deviations of the dSW differences rise from ~30 W m −2 to ~130 W m −2 when cloud coverage increases from 5% to 80%.The errors of CERES estimated net radiation will be directly transferred to the results of EDVI-based LE cal.
Due to the inhomogeneity of vegetation coverage on the ground and the inconsistency of spatial resolution among the input datasets (Table 2), the matching of these data over the selected area (0.15º × 0.15º) around each forest flux tower may introduce errors.In addition, spatial mismatch of the geolocation between satellite field of view to the forest sites also can result in some poor performances in validations.For example, the flux tower at DHS is on a steep 30 • slope and is close to a city, and there is a river located in the southeast.However, in our algorithm, the input data averaged over the 0.15º × 0.15º area surrounding the flux tower are used to estimate the regional mean LE.This certainly introduced additional discrepancies between the satellite EDVI-based LE cal and the in situ measured LE obs .
In addition, it should be noted that the validation results are also affected by the accuracies of the in situ measurements.The EC method may produce potential errors or uncertainty of 10%-30% [5].Also, the EC flux towers frequently suffer energy closure errors.In this study, the three ChinaFLUX forest sites achieve approximately only 70% energy balance closure during daytime, with about 30% of variation unexplained [35].
Although affected by the above uncertainties, the accuracies of instantaneous LE estimations in this study are comparable with those in Li et al. [29] which utilized both local in-situ measurements and satellite retrieved EDVI as data input.The results of this study demonstrate that it is feasible to operate the EDVI-based LE algorithm under both clear and cloudy skies with all satellite observations as data input.

Pros and Cons of this method
In comparison to Li et al. [29], the biggest improvement of this method is that its inputs do not depend on any in-situ measurements.This feature makes it applicable to any places with available satellite observations.Further, with proper calibration and validation, it can be developed to be a global algorithm in the future.
Compared to other optical VIs-based methods, the most significant advantage of this method is its capability for estimating ET related LE under cloudy sky.As shown in Section 3.4, no matter in what cloud fraction condition, the retrieved LE keeps small bias and good consistency with the in situ measurements.The most serious disadvantage of this method could be its poor spatial resolution, i.e., ~20 km.It is certain that the vegetation states and the evapotranspiration rate can vary significantly over this scale since the heterogeneous surfaces.This shortcoming may offset its merit of being useful under all weather conditions.A study of downscaling this method to the finer spatial resolution of optical VI is undergoing in our lab.
It should be noted that some parameterization schemes in our EDVI-based ET method need to be further improved in the future, especially for those resistance estimations based on EDVI and dEDVI.Since the selected three forests are similar to the Harvard forest where the initial method was developed [29], EDVI-based resistance schemes (e.g., a and b in Equation ( 8)) are thus assumed to be the same in this study.In addition, this method is recommended in flat or moderate terrain due to the associated errors with large scale satellite imagery of steep slope and LE estimation.
The real LE in the forest can be significantly affected by heavy rainfall due to the enhanced evaporation from interception water on the leaves.However, such an effect cannot be captured by the current method because it only takes into account one evapotranspiration source from inner vegetation water.In contrast, the plant water deficit will induce leaf stomatal closure to prevent excessive water deficits.However, this response is not described in the current model.A further study to consider more ET sources and the phenology response of vegetation to drought should be conducted to improve the retrieving performance.
Besides EDVI, there are several other vegetation water content related indexes that have been developed.For example, the normalized difference water index (NDWI) [67] and land surface water index (LSWI) [68].Due to the differences in their physical connections to vegetation water and their distinct spatial and temporal resolution, they may be used to estimate ET and LE with particular advantages and disadvantages, respectively.A comprehensive comparison study among them will be valuable for improving the satellite based global estimation of evapotranspiration and latent heat flux.

Figure 1 .
Figure 1.Locations of three ChinaFLUX forests sites.More information regarding these pictures is available at: http://www.chinaflux.org/.

Figure 1 .
Figure 1.Locations of three ChinaFLUX forests sites.More information regarding these pictures is available at: http://www.chinaflux.org/.

Figure 2 .
Figure 2. Time series of EDVI, in-situ measured LE (averaged over 2 hours around the satellite overpass), and 16-day NDVI from 2003 to 2005 at three forest sites of ChinaFLUX.QYZ forest site had a serious imbalance of energy in 2005.

Figure 2 .
Figure 2. Time series of EDVI, in-situ measured LE (averaged over 2 h around the satellite overpass), and 16-day NDVI from 2003 to 2005 at three forest sites of ChinaFLUX.QYZ forest site had a serious imbalance of energy in 2005.

Figure 3 .
Figure 3.Time series and scatter plots of in-situ measured solar radiation (averaged over 2 hours around the satellite overpass) and matched surface downward shortwave radiative fluxes (dSW) from CERES SSF under all-sky conditions.

Figure 3 .
Figure 3.Time series and scatter plots of in-situ measured solar radiation (averaged over 2 h around the satellite overpass) and matched surface downward shortwave radiative fluxes (dSW) from CERES SSF under all-sky conditions.

Figure 3 .
Figure 3.Time series and scatter plots of in-situ measured solar radiation (averaged over 2 hours around the satellite overpass) and matched surface downward shortwave radiative fluxes (dSW) from CERES SSF under all-sky conditions.

Figure 4 .
Figure 4. Comparisons between in-situ measured PAR and estimated PAR (dSW*1.70).dSW is surface downward shortwave radiative fluxes from CERES SSF in all sky conditions.

Figure 4 .
Figure 4. Comparisons between in-situ measured PAR and estimated PAR (dSW*1.70).dSW is surface downward shortwave radiative fluxes from CERES SSF in all sky conditions.

Figure 5 .
Figure 5.Time series and scatter plots of in-situ measured air temperature (Ta), 2-meter temperature from ECMWF reanalysis (t2m-ECMWF) and 2-meter temperature from NCEP reanalysis (t2m-NCEP).In-situ measurements at DHS, QYZ, and CBS are at 20 m, 23 m, and 26 m, respectively.In scatterplots, t2m-ECMWF and t2m-NCEP are marked by black and blue color, respectively.R1 represents the correlation coefficient (R) of in-situ Ta and t2m-ECMWF.R2 represents the R of in-situ Ta and t2m-NCEP.BIAS1 represents the difference between mean t2m-ECMWF and in-situ Ta.BIAS2 represents the difference between mean t2m-NCEP and in-situ Ta.

Figure 5 .
Figure 5.Time series and scatter plots of in-situ measured air temperature (Ta), 2-meter temperature from ECMWF reanalysis (t2m-ECMWF) and 2-meter temperature from NCEP reanalysis (t2m-NCEP).In-situ measurements at DHS, QYZ, and CBS are at 20 m, 23 m, and 26 m, respectively.In scatterplots, t2m-ECMWF and t2m-NCEP are marked by black and blue color, respectively.R1 represents the correlation coefficient (R) of in-situ Ta and t2m-ECMWF.R2 represents the R of in-situ Ta and t2m-NCEP.BIAS1 represents the difference between mean t2m-ECMWF and in-situ Ta.BIAS2 represents the difference between mean t2m-NCEP and in-situ Ta.

Figure 6 .
Figure 6.Time series of in-situ LE (LE obs ), estimated LE (LE cal ), and their differences at three forest sites from 2003 to 2005.LE obs was averaged over 2 h around the satellite overpass.The dashed lines are ±150 Wm −2 .

Figure 6 .
Figure 6.Time series of in-situ LE (LEobs), estimated LE (LEcal), and their differences at three forest sites from 2003 to 2005.LEobs was averaged over 2 hours around the satellite overpass.The dashed lines are ±150 Wm -2 .

Figure 7 .
Figure 7. Comparisons between daily LEobs and LEcal at three sites from 2003 to 2005.The two dashed lines are the 1:1 line±RMSE.Solid circles are samples severely contaminated by precipitation.

Figure 7 .
Figure 7. Comparisons between daily LE obs and LE cal at three sites from 2003 to 2005.The two dashed lines are the 1:1 line±RMSE.Solid circles are samples severely contaminated by precipitation.

21 Figure 8 .
Figure 8.Comparison results of all LE estimations (white bars) and the LE estimations after removing the heavily rain-contaminated days (black bars).

Figure 8 .
Figure 8.Comparison results of all LE estimations (white bars) and the LE estimations after removing the heavily rain-contaminated days (black bars).

Figure 9 .
Figure 9. Comparisons between monthly mean LEobs and LEcal at DHS, QYZ, and CBS.Horizontal and vertical error bars stand for the standard deviations of LEobs and LEcal.Numbers in the brackets are the statistic results after removing the heavily rain-contaminated samples.

Figure 9 .
Figure 9. Comparisons between monthly mean LE obs and LE cal at DHS, QYZ, and CBS.Horizontal and vertical error bars stand for the standard deviations of LE obs and LE cal .Numbers in the brackets are the statistic results after removing the heavily rain-contaminated samples.

Figure 11 .
Figure 11.Statistical metrics at three sites under different cloud cover.Slope is the slope of fit line between LE cal and LE obs .

Table 2 .
Multiple data sources used in EDVI-based LE method.

Table 3 .
Summary of comprehensive metrics for the results of three forest sites.LE cal is the mean value of all estimations.LE obs is the mean value of all in-situ measurements.BIAS = LE cal − LE obs .