**1. 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<sup>2</sup> 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 di fference 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–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–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 a ffected 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 di fference 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–19,23,24] is insu fficient to reflect the fast dynamics of conductance and LE a ffected 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 di ffers 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 di fference vegetation index (EDVI) at a mid-latitude forest (Harvard forest). The EDVI is defined as (MLSE19-MLSE37)/(MLSE19+MLSE37), where MLSE19 and MLSE37 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–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 H2O, CO2, 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.

#### **2. Data and Method**

## *2.1. Site Descriptions*

The three selected forest towers in this study are shown in Figure 1. Three years of in-situ measurements (2003–2005) at these sites are available in this study. And at each site, we split the measurements into di fferent time periods for calibration and validation study. These sites are in di fferent 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.

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



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 (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.

#### *2.2. 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 a ffect 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].


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

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 (U10m, U10m) 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/).

#### *2.3. 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–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.

$$\rm ET\_{\rm veg} = EF\_{\rm veg} \cdot (R\mathbf{n} - G)\_{\rm veg} = EF\_{\rm veg} \cdot (R\mathbf{n} - G) \cdot \rm VFC \tag{1}$$

#### 2.3.1. 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:

$$\mathbf{G} = \text{Rn} \left[ \mathbf{f}\_{\text{veg}} + (\mathbf{1} - \text{VFC}) \left( \mathbf{f}\_{\text{soil}} - \mathbf{f}\_{\text{veg}} \right) \right] \tag{2}$$

where fveg (set as 0.05) and fsoil (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]:

$$\text{VFC} = \frac{\text{NDVI} - \text{NDVI}\_0}{\text{NDVI}\_\text{os} - \text{NDVI}\_0} \tag{3}$$

where NDVI ∞ and NDVI0 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 di fferent vegetation types, but they are relatively stable over forest types [50]. Values of NDVI0 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 NDVI0, 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.

**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.

#### 2.3.2. Estimation of Evaporative Fraction of Forest

The estimation of EFveg is from the same method of Nishida et al. [52] and Li et al. [29]:

$$\text{EF}\_{\text{veg}} = \frac{\alpha \Delta}{\Delta + \gamma (1 + \mathbf{r}\_{\text{c}}/2\mathbf{r}\_{\text{a}})} \tag{4}$$

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), ra is the aerodynamic resistance (s/m), and rc 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]:

$$\Delta = \frac{26,297.76}{\left(\text{Ta} - 29.65\right)^2} \exp\left(\frac{17.67(\text{Ta} - 273.15)}{\text{Ta} - 29.65}\right) \tag{5}$$

where Ta is the air temperature at 2 m (t2m) from ERA-20C.

The aerodynamic resistance ra 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 ra of forests:

$$\mathbf{r\_a} = \begin{cases} \frac{1}{0.008 \mathbf{U\_{50m}}} & \text{forest types} \\ \frac{1}{0.003 \mathbf{U\_{1m}}} & \text{gass and corps} \end{cases} \tag{6}$$

where U50m is the wind speed at 50 m (ms−1) which is calculated by averaging the mean values of 10-m wind (U10m) and 100-m wind (U100m) from ERA-20C.

The canopy resistance rc can be parameterized by the use of Jarvis-type equation.

$$\mathbf{r\_c} = \left[ \frac{\mathbf{f\_1(Ta)f\_2(PAR)f\_3(VPD)f\_4(\Psi)f\_5(CO\_2)}}{\mathbf{r\_{cmin}}} + \frac{1}{\mathbf{r\_{cuitcle}}} \right]^{-1} \tag{7}$$

where rc has two components: the stomatal resistance (rstomatal) and the cuticle resistance (rcuticle). rcuticle was set to be 10<sup>5</sup> ms<sup>−</sup><sup>1</sup> based on other studies [11].

The stomatal resistance rstomatal consists of five stress functions and the minimum canopy resistance (rcmin). The five stress functions are used to quantify the impacts on rc imposed by environmental conditions. The functions of air temperature (f1(Ta)) and photosynthetically active radiation (f2(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 (CO2). In this study, we parameterize their stress functions (i.e., f3f4f5) as a whole by using their fair correlation with dEDVI based on Li et al. [29]:

$$\mathbf{f\_3(VPD)f\_4(\Psi)f\_5(CO\_2) = [\mathbf{a} - \mathbf{b} \times \text{dEDVI}]^{-1}} \tag{8}$$

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 (NEDVI) as (EDVI-EDVImin)/(EDVImax-EDVImin), where EDVImin and EDVImax are respectively the minimum and maximum EDVI values during growing seasons at each site [27,29]. The seasonal variation of EDVI (i.e., NEDVI) represents the slow change of VWC which is highly correlated to the vegetation resistance of LE. NEDVI is thus used to determine rcmin which is the key parameter for the estimation of rstomatal. Li et al. [29] found that there was a general anti-correlation between NEDVI and canopy resistance. The minimum canopy resistance can be parameterized as:

$$\mathbf{r}\_{\rm cmin} = \mathbf{r}\_{\rm cmin0} / ^{\rm N\_E}\text{EDVI} \tag{9}$$

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 rcmin0. 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 rcmin (G−<sup>1</sup> smax) ranges from 76.9 s/m to 30.8 s/m. On the basis of this, rcmin0 was set to be a reference value, 50 s/m, for all three forest types in our study.
