Temporal Upscaling and Reconstruction of Thermal Remotely Sensed Instantaneous Evapotranspiration

Abstract

Currently, thermal remote sensing-based evapotranspiration (ET) models can only calculate instantaneous ET at the time of satellite overpass. Five temporal upscaling methods, namely, constant evaporative fraction (ConEF), corrected ConEF (CorEF), diurnal evaporative fraction (DiEF), constant solar radiation ratio (SolRad), and constant reference evaporative fraction (ConET_rF), were selected to upscale the instantaneous ET to daily values. Moreover, five temporal reconstruction approaches, namely, data assimilation (ET_EnKF and ET_SCE_UA), surface resistance (ET_SR), reference evapotranspiration (ET_ET_rF), and harmonic analysis of time series (ET_HANTS), were used to produce continuous daily ET with discrete clear-sky daily ET values. For clear-sky daily ET generation, SolRad and ConET_rF produced the best estimates. In contrast, ConEF usually underestimated the daily ET. The optimum method, however, was found by combining SolRad and ConET_rF, which produced the lowest root-mean-square error (RMSE) values. For continuous daily ET production, ET_ET_rF and ET_SCE_UA performed the best, whereas the ET_SR and ET_HANTS methods had large errors. The annual ET distributions over the Beijing area were calculated with these methods. The spatial ET distributions from ET_ET_rF and ET_SCE_UA had the same trend as ETWatch products, and had a smaller RMSE when compared with ET observations derived from the water balance method.

1. Introduction

Evapotranspiration (ET) is one of the most important components of surface energy budgets and hydrologic cycles. Flux observation networks (e.g., FluxNet) and enhanced experiments (e.g., HiWATER) have been set up to monitor long-term surface heat fluxes over different kinds of land cover [1,2,3]. However, in situ measurements are expensive and difficult to extend to continental scales. Thus, models have been developed to map ET over continental scales based on remotely sensed information [4,5,6,7,8,9,10,11]. Unfortunately, thermal remotely sensed models can only detect a snapshot of ET spatial distributions at the time of satellite overpass. The snapshot of the ET map cannot satisfy practical applications, such as, irrigation judgment, and water resource planning and management. To satisfy such applications, the instantaneous ET needs to be upscaled to daily values for clear-sky days. Moreover, continuous daily ET values should be generated with temporal reconstruction approaches to fulfill long-term hydrologic requirements.

There are many methods to upscale instantaneous ET to daily values. The most popular approach is the constant evaporative fraction (EF) method, which assumes that EF (EF is defined as the ratio of latent heat flux and available energy) is constant during the daytime [12]. Thus, the daily ET can be obtained with daily available energy (difference of net radiation, Rn, and soil heat flux, G₀) and remotely sensed instantaneous EF. However, the daytime ET is usually underestimated by 5%–10% with the EF constant method [13,14]. The diurnal variations of EF have been found to be constant during the daytime (9:00–16:00 local time), with a slight increase during the afternoon [15]. The underestimation can be corrected by adding 10% ET directly [16]. The estimation errors using the constant EF method are mainly caused by: (1) the difference between instantaneous EF and daytime EF; (2) the difference between daytime and daily available energy; and (3) the difference between daytime and daily latent heat flux. Daily ET estimates can be corrected by eliminating these differences [17]. Another ET upscaling method is to assume the ratio between latent heat flux and other variables (e.g., solar radiation, reference evapotranspiration) is constant during the daytime. Latent heat flux and solar radiation have similar diurnal variations (e.g., sine functions), and daily ET can be estimated with the relationship between latent heat flux and solar radiation [18,19]. The reference ET (ET_r) can depict the diurnal variations of net radiation, air temperature, wind speed, and relative humidity. The ratio between instantaneous ET and ET_r is assumed to be constant during the daytime, and daily ET can be estimated with daily ET_r [20,21]. Most of these methods have been compared over different land cover types, and the feasibility of each approach has been noted [22,23,24].

The upscaling methods can only be used to estimate daily ET with the instantaneous ET values over clear-sky days. This can be an issue, because it is difficult to obtain cloud-free thermal remote sensing data for every day of the year. It is, therefore, necessary to develop algorithms to reconstruct continuous daily ET variations with discrete clear-sky daily ET estimates. The Penman-Monteith equation has been used to calculate global ET with the revised surface resistance algorithm [25,26]. For clear-sky days, the surface resistance of Penman-Monteith can be retrieved by thermal remotely sensed daily ET. Then, continuous daily surface resistance can be gap filled with the two neighboring clear-sky surface resistances. Finally, continuous daily ET values can be estimated with the reconstructed surface resistance, leaf area index (LAI), and meteorological elements. This ET temporal reconstruction method, namely, the surface resistance approach, has been used for long-term regional and global ET production [25,26,27].

The reference evapotranspiration approach is a commonly used temporal reconstruction method [28]. This approach assumes that the ratio (ET_rF) between clear-sky daily ET and ET_r has liner variations between two nearby clear-sky days. The cloudy-day ET_rF can be linearly interpolated with the clear-sky day ET_rF, which is like generating the seasonal crop coefficient (Kc) with discrete clear-sky values. Finally, continuous daily ET can be obtained with ET_r and the interpolated ET_rF. The reference evapotranspiration method has been used in the METRIC (mapping evapotranspiration at high resolution with internalized calibration) and SEBAL (surface energy balance algorithm for land) models to produce regional total ET that is vital to irrigation management [28,29].

The harmonic analysis of time series (HANTS) method was first developed to reconstruct cloud contaminated NDVI observations with clear-sky data at prescribed times [30,31]. The HANTS algorithm considers only the most significant frequencies of the time series data with the least squares curve fitting method [30]. With time iteration, the large positive and negative values are rejected and reconstructed to a smoother data curve. The seasonal ET curve can also be reconstructed with the HANTS algorithm by inputting clear-sky ET, and the reconstructed monthly ET has been found to be more reliable than the daily results [32].

Data assimilation is a method used to produce a more reliable land surface state by assimilating spatially sparse in situ observations or temporal instantaneous remote sensing data [33,34]. Various data assimilation algorithms (e.g., ensemble Kalman filter (EnKF), shuffled complex evolution method developed at The University of Arizona (SCE_UA), particle filter, variational method, etc.) have been developed during the last two decades [35,36,37,38]. Data assimilation methods have been used to generate temporal continuous soil moisture, soil temperature, surface fluxes, NDVI and LAI, by assimilating remote sensing data into land surface models [39,40,41,42,43,44,45,46]. For surface flux gap filling, a data assimilation scheme has been developed to reconstruct surface flux data from eddy covariance measurements [47].

Thermal remotely sensed instantaneous ET should be extended to daily scale over cloud-free days and reconstructed to continuous daily values, to fulfill the requirements of water planning and management, and other applications. Although temporal upscaling and reconstruction methods have been used broadly, they have not been compared extensively, especially the temporal reconstruction approaches. In this study, five ET upscaling methods (from instantaneous ET to daily over a clear-sky day) were compared with in situ observations and an optimum method is proposed. To reconstruct clear-sky daily ET to continuous daily, data assimilation methods based on the Penman-Monteith equation were developed with two data assimilation techniques (EnKF and SCE_UA). The developed data assimilation methods were also compared with three commonly used temporal reconstruction approaches (surface resistance method, reference evapotranspiration approach, and HANTS package). Finally, the annual remotely sensed ET over Beijing area was estimated with these methods and compared with the water balance ET data.

2. Methodology

2.1. Remotely Sensed ET Model

In this study, a remotely sensed ET model was used to estimate the instantaneous ET during satellite overpass time [5]. The model’s parameterization scheme is based on the SEBS model [4] and revised according to the heterogeneous land surface of the Beijing area. Using this model, the instantaneous net radiation (Rn), ground heat flux (G) and sensible heat flux (H) can be calculated with atmospheric interpolated forcing data (wind speed, air temperature/humidity, etc.) and remote sensing data (land surface temperature, albedo, emissivity, NDVI, etc.). With Rn, G, and H, the latent heat flux (LE) can be obtained with the residual method via LE = Rn − G − H.

2.2. Instantaneous ET Upscaling Methods

Five methods were selected to upscale instantaneous ET to daily values and are summarized below.

2.2.1. Constant Evaporative Fraction Method (ConEF)

The ConEF method treats the evaporative fraction (EF, EF = LE/Rn − G) as a constant during the daytime and that daily ET can be calculated using daily available energy via [12],

λET_d = [EF_i(Rn − G)_d]

(1)

where λ is latent heat of vaporization (J∙kg⁻¹), ET_d is the daily ET (mm/day), EF_i represents the evaporative fraction during satellite overpass time i (-), (Rn − G)_d represents daily available energy (W∙m⁻²).

2.2.2. Corrected Constant Evaporative Fraction Method (CorEF)

The ConEF method usually leads to a 5%–10% underestimation of daily ET using EF around noon instead of daily EF [13,14]. This method adds 10% of EF to correct the error caused by ConEF via [16],

λET_d = [1.1 × EF_i (Rn − G)_d]

(2)

2.2.3. Diurnal Evaporative Fraction Method (DiEF)

The DiEF method introduced three parameters to correct ConEF method, including (1) the difference between instantaneous and daytime EF, β_EF; (2) the difference between daytime-averaged and daily-averaged available energy, β_A; (3) the difference between daytime and daily ET, β_LE. This method can be expressed as [17],

λET_d = EF_i (Rn − G)_d × (β_LE × β_A × β_EF)

(3)

$${\beta }_{EF}=1.26\frac{a_A\text{sin}\left(\text{π}\frac{t_{*}}{\text{∆}{t}_d}\right)+c_A}{\text{sin}\left(\text{π}\frac{t_{*}}{\text{∆}{t}_d}-0.1\right)}$$

(4)

$${\beta }_A=1.26\frac{0.5\text{∆}{t}_d}{\text{∆}{t}_n\left[c_A-\left|c_A\right|{\text{cos}}^5\left(\frac{\pi Y_d}{365}\right)\right]+0.5\text{∆}{t}_d}$$

(5)

where t_* is the length of time from sunrise (9:00), ∆t_d is the length of time from sunrise to sunset (9 h), ∆t_n is the length of time from sunset to sunrise (∆t_n = 24 − ∆t_d), Y_d is the day of year, a_A and c_A are empirical parameters that are set to 0.83 and −0.03, respectively, and β_LE is set to 1.1.

2.2.4. Constant Solar Radiation Ratio Method (SolRad)

The SolRad method assumes the ratio of latent heat flux and solar radiation is constant at the daily time scale, which is expressed as [18,19],

ES = LE/Rs

(6)

λET_d = ES_i(Rs)_d

(7)

where ES_i is the ratio between latent heat flux and solar radiation at time i (satellite overpass time), and Rs is the solar radiation.

2.2.5. Constant Reference Evaporative Fraction Method (ConET_rF)

The ConET_rF method assumes the ratio of ET and ET_r is constant on the daily time scale. The daily ET can be calculated using [20,21],

ET_rF = ET/ET_r

(8)

ET_d = (ET_rF)_i(ET_r)_d

(9)

where ET_rF is the ratio between ET and ET_r, and ET_r can be computed following the FAO-56 paper [48].

2.3. Continuously Daily ET Temporal Reconstruction Approaches

In this section, five approaches were introduced to reconstruct continuous daily ET with clear-sky daily ET estimates.

2.3.1. Data Assimilation Method

In the developed data assimilation scheme, the Penman-Monteith equation was used as the model operator and the clear-sky daily ET estimates were assimilated into the Penman-Monteith equation to generate the continuous daily ET. The Penman-Monteith equation can be expressed as [25],

$$\lambda E{T}_d=\frac{\text{∆}\left(R_n-G\right)+{\rho }_a{c}_P\frac{(e_s-e_a)}{r_a}}{\text{∆}+\gamma (1+\frac{r_s}{r_a})}$$

(10)

where ∆ is the slope of the curve relating saturated water vapor pressure to temperature (Pa∙K⁻¹), ρ_a is the air density (Kg∙m⁻³), c_p represents the specific heat capacity of air (J∙Kg⁻¹∙K⁻¹), r_s is the surface resistance (s∙m⁻¹), r_a is aerodynamic resistance (s∙m⁻¹), e_a is the actual water vapor pressure (Pa), e_s is the saturated water vapor pressure (Pa), and γ is the psychrometric constant (Pa∙K⁻¹). The parameterization scheme of the surface resistance of this study follows Mu et al. [25]. The ET over the open water surface can be calculated as,

$$\lambda E{T}_d=\frac{\text{∆}}{\text{∆}+\gamma }\left(R_n-G\right)+\frac{\gamma }{\text{∆}+\gamma }{E}_a$$

(11)

where E_a is the evapotranspiration rate (mm/day).

Before constructing the data assimilation scheme, the extended Fourier amplitude sensitivity test (EFAST) [49] was used to rank the model variables that are important for ET estimates with the Penman-Monteith equation. According to the sensitivity test, the available energy (Rn − G) and surface resistance (r_s) were the most important variables in the Penman-Monteith equation. Thus, two tunable parameters were added along with the available energy (Rn − G) and surface resistance (r_s), and the data assimilation scheme was constructed as,

$$\lambda E{T}_d=\frac{\text{∆}\alpha \left(R_n-G\right)+{\rho }_a{c}_P\frac{(e_s-e_a)}{r_a}}{\text{∆}+\gamma (1+\frac{\beta r_s}{r_a})}$$

(12)

Over the open water surface, the data assimilation scheme was constructed as,

$$\lambda E{T}_d=\frac{\alpha \text{∆}}{\text{∆}+\gamma }\left(R_n-G\right)+\frac{\beta \gamma }{\text{∆}+\gamma }{E}_a$$

(13)

where α and β are two tunable parameters in the data assimilation scheme.

To update the two parameters with clear-sky ET estimates, two data assimilation algorithms were used, namely, the EnKF and SCE_UA method. The EnKF and SCE_UA algorithms are popular methods in data assimilation [34]. The detailed implementation of EnKF and SCE_UA can be found in Xu et al. [42]. The two parameters were updated with previous clear-sky ET observations when using the EnKF method. However, the SCE_UA method can use clear-sky ET observations before and after the reconstruction day. Thus, the SCE_UA approach can make full use of more observations. In this study, two parameters were updated with four neighbored clear-sky ET observations when using the SCE_UA method, and one previous ET observation when using the EnKF method. The data assimilation methods are defined as ET_EnKF and ET_SCE_UA hereafter.

2.3.2. Surface Resistance Method

The surface resistance method was built, based on the Penman-Monteith equation, to reconstruct ET estimates over cloudy conditions [25]. The surface resistance (r_s) in the Penman-Monteith equation represents wet conditions of the land surface and the stomatal conductance of vegetation. Over a clear-sky day, r_s can be retrieved with the Penman-Monteith equation based on remotely sensed ET estimates. The cloudy-sky day r_s can be computed with the r_s neighboring clear-sky days, LAI, and the daily reduction functions for stomatal conductance to minimize the air temperature m(T_min) and vapor pressure deficit m(VPD). The surface resistance on cloudy-sky days can be calculated using the following equation,

$$r_{s\_cld}=\frac{LA{I}_{clr}\times r_{s\_clr}}{LA{I}_{cld}\times m\left(T_{min}\right)\times m\left(VPD\right)}$$

(14)

where r_{s_cld} and r_{s_clr} are the surface resistance on cloudy-sky and clear-sky days (s∙m⁻¹), and LAI_cld and LAI_clr represent the leaf area index on cloudy-sky and clear-sky days (m²∙m⁻²). A table of m(T_min) and m(VPD) functions based on different land cover types was built and used in this study [25]. The surface resistance method is defined as ET_SR hereafter.

2.3.3. Reference Evapotranspiration Method

This method is the same as the ConET_rF upscaling method described in Section 2.2. The ratio of ET and ET_r on cloudy-sky days can be interpolated by ET_rF on neighboring clear-sky days. The equation for this method is shown below,

ET_cld = (ET_rF)_clr(ET_r)_cld

(15)

where ET_cld is the cloudy-sky ET estimates, (ET_rF)_clr is the interpolated clear-sky ET_rF, and (ET_r)_cld represents cloudy-sky ET_r. The reference evapotranspiration method is defined as ET_ET_rF hereafter.

2.3.4. Harmonic ANalysis of Time Series (HANTS) Method

The HANTS method is a time series analysis method based on Fourier transforms [31]. The HANTS method allows the user to choose the frequencies of periodic functions to estimate the time series of variables with discrete observations. The gaps in ET can be filled with discrete clear-sky ET values. The parameters used for the HANTS method are displayed below in Table 1.

Table 1. HANTS parameters setting for two study sites: Guantao and Arou.

Parameters	Guantao	Arou
Minimum Valid Data	0	0
Maximum Valid Data	550	320
Number of Frequencies	6	4
Period(s)	360 130 90 70 50 30	360 200 150 100
Fit Error Tolerance	2	2

Table 1. HANTS parameters setting for two study sites: Guantao and Arou.

Parameters	Guantao	Arou
Minimum Valid Data	0	0
Maximum Valid Data	550	320
Number of Frequencies	6	4
Period(s)	360 130 90 70 50 30	360 200 150 100
Fit Error Tolerance	2	2

3. Data Sets

3.1. In Situ Meteorological and Validation Data

Data from two sites namely Guantao and Arou, located in different natural zones (the eastern monsoon area, and the Qinghai-Tibet alpine area) of China, were selected to validate various temporal upscaling and reconstruction methods described in Section 2. Data from another two sites, namely, Daxing and Miyun, located in Beijing, China were selected for access to regional application results. Guantao is a cropland site with winter wheat/maize and cotton rotation, and is located in the Hebei Province. Daxing is a cropland site covered with winter wheat/maize and vegetables located southeast of Beijing. Miyun is an orchard site covered with orchards and maize, and is located in a valley northeast of Beijing. Ground measurements at Guantao, Daxing, and Miyun sites consist of multiscale surface flux and meteorological measurements in the Hai River Basin [50]. Arou is a grassland site covered with alpine meadows, and is located in the Qinghai Province. The in situ measurements at the Arou site are a part of the Watershed Allied Telemetry Experimental Research (WATER) [51]. The vegetation growing time was determined according to vegetation phenology data provided by the Chinese Ecosystem Research Network (CERN) (http://www.cern.ac.cn/).

The main characteristics of the four sites are summarized in Table 2. The meteorological data (wind speed/direction, air temperature/humidity, air pressure, precipitation, shortwave/longwave radiation) were measured with automatic weather stations (AWS) mounted on flux towers. Eddy covariance (EC) instruments were also equipped on flux towers, and measured latent heat flux data were used as input and validation data. Original EC data were collected at a sampling frequency of 10 Hz, and processed using the post processing software Edire (http://www.geos.ed.ac.uk/abs/research/micromet/EdiRe). The post processing included spike removal, lag correction of H2O/CO2 relative to the vertical wind component, sonic virtual temperature correction, the performance of the planar fit coordinate rotation, corrections for density fluctuation (WPL-correction), and frequency response correction, etc. [50,52]. The ground heat flux was measured by a soil heat flux plate buried at a specific depth and was corrected to the soil surface (0 cm) with soil temperature and moisture observations [53]. A large aperture scintillometer (LAS), consisting of a transmitter and receiver, was used to measure the sensible heat flux over a larger spatial scale at the Daxing and Miyun sites. At the Daxing site, the LAS transmitter and receiver were installed on two opposing towers, positioned in the southwest and northeast. At the Miyun site, the LAS transmitter and receiver were installed on two opposing hills along the valley in the northeast and southwest. Data were recorded at a sampling frequency of 1 Hz using a Kipp & Zonen LAS at the two sites. With measured net radiation and ground heat flux, the LAS latent heat flux can be calculated as the residual of the surface energy balance (LE = Rn − G₀ − H). The detailed height/depth of surface heat flux measurements and LAS path length at the experiment sites were summarized in Table 3.

All the meteorological and flux data were processed in 30-min intervals with strict quality control procedures [50,52]. The EC-derived ET data in 2010 at the Guantao site, and in 2009 at the Arou site, were selected for use in the method comparisons. The 30-min data at noon (12:00 PM) on clear-sky days were used to upscale instantaneous ET to daily value. The accumulated daily ET on clear-sky days was used to access the accuracy of temporal upscaling methods. The continuous daily ET was reconstructed with the accumulated clear-sky daily ET. The 30-min ET data for cloudy-days were processed to daily values to validate the temporal reconstruction results. ET derived from LAS instruments in the year 2009 at the Daxing and Miyun sites was used to validate the remotely sensed ET estimates.

Table 2. Summary of the surface characteristics of the four sites.

Site	Latitude	Longitude	Location (Province)	Land Cover	Vegetation Growing Season (DOY)
Guantao	36.52°N	115.13°E	Hebei	Winter wheat/maize and cotton	100–161 & 182–283
Arou	38.04°N	100.91°E	Qinghai	Alpine meadow	105–273
Daxing	39.62°N	116.42°E	Beijing	Winter wheat/maize and vegetables	100–161 & 182–283
Miyun	40.63°N	117.32°E	Beijing	Orchard and maize	100–283

Table 2. Summary of the surface characteristics of the four sites.

Site	Latitude	Longitude	Location (Province)	Land Cover	Vegetation Growing Season (DOY)
Guantao	36.52°N	115.13°E	Hebei	Winter wheat/maize and cotton	100–161 & 182–283
Arou	38.04°N	100.91°E	Qinghai	Alpine meadow	105–273
Daxing	39.62°N	116.42°E	Beijing	Winter wheat/maize and vegetables	100–161 & 182–283
Miyun	40.63°N	117.32°E	Beijing	Orchard and maize	100–283

Table 3. Summary of surface heat flux measurement height/depth (m).

Instrument	Variable	Guantao	Arou	Daxing	Miyun
EC	Sensible and latent heat flux	15.6	3.2	3.0	26.7
LAS	Sensible heat flux			27.0 Path length 2480.0	35.9 Path length 2420.0
Soil heat flux plate	Soil heat flux	0.02	0.05	0.02	0.02

Table 3. Summary of surface heat flux measurement height/depth (m).

Instrument	Variable	Guantao	Arou	Daxing	Miyun
EC	Sensible and latent heat flux	15.6	3.2	3.0	26.7
LAS	Sensible heat flux			27.0 Path length 2480.0	35.9 Path length 2420.0
Soil heat flux plate	Soil heat flux	0.02	0.05	0.02	0.02

Figure 1. The temporal distribution of clear-sky daily evapotranspiration (ET) over the whole year.

Figure 1. The temporal distribution of clear-sky daily evapotranspiration (ET) over the whole year.

The data from clear-sky days were selected for temporal upscaling method comparisons at the Guantao and Arou sites, according to the following three conditions: (1) the data must have a double-zero MODIS quality assurance flag, indicating no cloud cover (QA = 00); (2) no missing data during a day; and (3) an energy balance ratio larger than 0.8 for EC flux data (EBR = (H + LE)/(Rn − G₀)). After this data selection procedure, there were 71 and 67 clear-sky days at the Guantao and Arou sites, respectively. In addition, data from 273 and 239 cloudy-sky days were used for ET temporal reconstruction approach comparisons at the Guantao and Arou sites, respectively. The temporal distributions of clear-sky days and continuous daily ET, measured at the Guantao and Arou sites, are shown in Figure 1.

3.2. Regional Application Data

The ET temporal upscaling and reconstruction methods were applied over the Beijing area, which is located in the northern part of the North China Plain. The land surface data of Beijing consist of distinct land cover types (cropland, garden-plot, forest, grassland, urban, road, water, and bare land), as shown below in Figure 2. The locations of four experiment sites and locations of Beijing are also shown In Figure 2.

Figure 2. Site locations and land use of Beijing in China, 2009.

Figure 2. Site locations and land use of Beijing in China, 2009.

The hourly meteorological data (wind speed, air temperature/humidity, and sunshine percentage) were collected from 20 weather stations within the Beijing area in 2009. The data required temporal and spatial interpolation for regional applications. The cubic spline interpolation method was used to interpolate the hourly air temperature/humidity observations and the satellite overpass time. The linear interpolation method was used to interpolate the hourly wind speed observations and the satellite overpass time. The Kriging and inverse distance weighted (IDW) methods were compared for regional air temperature and wind speed interpolation. Through comparison, the Krigring method had higher accuracy for air temperature/humidity interpolation than the IDW method, while the IDW method had higher accuracy for wind speed interpolation than Krigring method (not shown). Thus, the Kriging method was used to interpolate the air temperature/humidity observations over sparse data sites to the remote sensing pixel scale (1 km in this study) over Beijing, and the IDW method was used to interpolate the wind speed observations to the same scale. The lapse rates were allowed in the Kriging method for spatial interpolation by including DEM data (obtained from United States Geological Survey, USGS) [54].

The MODIS products (collection 5) were used to estimate clear-sky instantaneous ET with the remotely sensed model proposed by Liu et al. [5]. These products were land surface temperature (LST, MOD11A1), albedo (MOD09A1), NDVI (MOD13A2), and leaf area index (LAI, MOD15A2). The spatial resolution of the LST, NDVI, and LAI products was 1 km and the albedo product was aggregated to the same spatial resolution of the other products. The ET products from MODIS (MOD16A3) [26] and ETWatch [27] were also acquired for comparison over the Beijing area. The land cover map was obtained by artificial interpretation using MODIS visible and infrared data.

The total ET was calculated with hydrological data over Beijing using the water balance equation ET = P − R + ∆S, in which P represent precipitation, R means runoff (including water outflow Q₁, water inflow Q₂, and emergency water diversion Q₃, R = Q₁ − Q₂ − Q₃), and ∆S is change of terrestrial water storage (including change of water storage in reservoirs ∆S₁, and change of ground water resources ∆S₂, ∆S = ∆S₁ + ∆S₂) [55]. According to Beijing Water Resource Bulletin 2009 published by the Beijing Water Authority (http://www.bjwater.gov.cn/pub/bjwater/), P = 448 mm, R = 2.6 × 10⁸ m³ (Q₁ = 8.23 × 10⁸ m³, Q₂ = 3.03 × 10⁸ m³, and Q₃ = 2.6 × 10⁸ m³), and ∆S = −7.66 × 10⁸ m³ (∆S₁ = −1.32 × 10⁸ m³, and ∆S₂ = −6.34 × 10⁸ m³). R and ∆S were transformed to the same unit of ET by dividing the area of Beijing (1.68 × 10¹⁰ m²), which equals 15.5 mm and −45.6 mm, respectively. With these data and the water balance equation, the annual ET of Beijing in 2009 was about 387 mm.

3.3. Regional Remotely Sensed ET Validation

The pixel-scale ET estimated by the remotely sensed model were validated with LAS measurements. One transmitter and one receiver installed on a pair of opposing towers form a complete LAS system, which can measure the average value of sensible heat flux between the two towers. The LAS measurements at the Daxing and Miyun sites were used as the validation data for model estimates. The source area of the LAS measurements was calculated using a footprint model as follows [50,52],

$$f_{LAS}\left(x,y,z_{eff}\right)={{\displaystyle \int }}_{x_2}^{x_1}W\left(x^{\prime }\right)f\left(x^{\prime }-x,y^{\prime }-y,z_{eff}\right)d{x}^{\prime }$$

(16)

where W(x’) is the path-weighting function of the LAS, x₁, x₂ are the locations of the LAS transmitter and receiver, respectively, x’, y’ are the points along the optical length of the LAS, x, y are the coordinates upwind of each point (x’, y’), and z_eff is the effective measurement height of the LAS.

The calculated LAS source areas at the Miyun and Daxing sites were overlaid with remote sensing pixels as shown in Figure 3. The relative weight of the source area was determined by using the footprint model (Equation (16)). The relative weight was highest near the center of the LAS optical path, and decreased as distance from the center increased. The ET estimates at each covered pixel were spatially aggregated into one value that can be compared with LAS measurements via [55],

$$Y_{weighted}={\displaystyle \sum }_{i=1}^n\left(x_i\times P_i\right)$$

(17)

where Y_weighted is the aggregated remote sensing ET, which is spatially the same as the LAS measurements, x_i is the weight of pixel i, P_i is the ET estimate at pixel i, and n is the number of pixels covered by the LAS source area.

Figure 3. LAS source area overlaid with 1 km MODIS pixel (left is the Miyun site, right is the Daxing site; LAS_T is the LAS transmitter and LAS_R is the LAS receiver; the background pix is derived from MODIS LAI products, MOD15A2).

Figure 3. LAS source area overlaid with 1 km MODIS pixel (left is the Miyun site, right is the Daxing site; LAS_T is the LAS transmitter and LAS_R is the LAS receiver; the background pix is derived from MODIS LAI products, MOD15A2).

4. Results and Discussions

4.1. In Situ Validation

The five temporal upscaling methods and five temporal reconstruction approaches were tested and compared with EC-derived ET observations at the Guantao (for 2010) and Arou (for 2009) sites. The mean relative error (MRE), root mean square error (RMSE) and correlation (R) were used to assess the results.

The temporal upscaling results from the five methods were compared with the daily clear-sky ET observations (see Table 4). The RMSE values of the five methods were all within 1 mm/day and the R values were all larger than 0.89. The SolRad and ConET_rF methods were the two best upscaling methods, as indicated by the smallest RMSE values and largest R values. The MRE values indicated that the ConEF method significantly underestimated the daily ET by 19.1% and 14.1% at the Guantao and Arou sites, respectively. Since many studies have indicated that “self-preserving” will lead to an approximately 5%–10% underestimation of the daytime ET, 10% was added to the ET values for the CorEF method, to correct the underestimation, which led to an 8.1% and 8.6% reduction in MRE at the Guantao and Arou sites. However, the CorEF method still underestimated the daily ET, because nighttime ET was not accounted for. The DiEF method accounted for the nighttime ET by using three factors that enhance the daily ET estimates with the daytime instantaneous data. The MRE values for the DiEF method were near zero. Similarly, the SolRad and ConET_rF methods had MRE values around zero, which indicates that these three methods had little bias in ET temporal upscaling. Moreover, SolRad and ConET_rF methods had smaller RMSE than DiEF, indicating that the ratio between latent heat flux (ET) and solar radiation or reference ET was more stable than EF. According to Gentine et al. [15], EF is higher during early morning and late afternoon than it is midday. Here, Similar to the finding of previous studies [23,24], this phenomenon caused underestimation of the daytime average EF from midday EF observations with the ConEF method. However, the ConET_rF method upscales ET with daily ET_r, which is better at capturing the impacts of advection, changing wind, and changing humidity conditions during the day with the reference ET included [28].

Table 4. The statistics of the five temporal upscaling methods at the Guantao and Arou sites.

Site	Season	Statistics	ConEF	CorEF	DiEF	SolRad	ConETrF
Guantao	Vegetation, growing	MRE (%)	−13.0	−4.3	4.4	−2.1	−1.4
		RMSE (mm/day)	0.79	0.80	0.96	0.70	0.65
		R (-)	0.93	0.93	0.93	0.93	0.93
	Vegetation, dormant	MRE (%)	−29.5	−22.4	−15.4	−11.2	0.91
		RMSE (mm/day)	0.58	0.51	0.54	0.33	0.33
		R (-)	0.96	0.96	0.96	0.98	0.98
	Overall	MRE (%)	−19.1	−11.0	−2.7	−6.9	−1.3
		RMSE (mm/day)	0.76	0.71	0.79	0.58	0.65
		R (-)	0.94	0.94	0.94	0.95	0.93
Arou	Vegetation, growing	MRE (%)	−15.4	−6.9	2.2	−3.9	0.79
		RMSE (mm/day)	0.61	0.54	0.60	0.29	0.49
		R (-)	0.89	0.89	0.89	0.97	0.98
	Vegetation, dormant	MRE (%)	−45.1	−39.6	−37.6	−5.9	−18.7
		RMSE (mm/day)	0.17	0.15	0.14	0.07	0.11
		R (-)	0.98	0.98	0.98	0.98	0.97
	Overall	MRE (%)	−14.1	−5.5	3.2	−3.4	0.9
		RMSE (mm/day)	0.64	0.65	0.70	0.21	0.39
		R (-)	0.89	0.89	0.89	0.99	0.96

Table 4. The statistics of the five temporal upscaling methods at the Guantao and Arou sites.

Site	Season	Statistics	ConEF	CorEF	DiEF	SolRad	ConETrF
Guantao	Vegetation, growing	MRE (%)	−13.0	−4.3	4.4	−2.1	−1.4
		RMSE (mm/day)	0.79	0.80	0.96	0.70	0.65
		R (-)	0.93	0.93	0.93	0.93	0.93
	Vegetation, dormant	MRE (%)	−29.5	−22.4	−15.4	−11.2	0.91
		RMSE (mm/day)	0.58	0.51	0.54	0.33	0.33
		R (-)	0.96	0.96	0.96	0.98	0.98
	Overall	MRE (%)	−19.1	−11.0	−2.7	−6.9	−1.3
		RMSE (mm/day)	0.76	0.71	0.79	0.58	0.65
		R (-)	0.94	0.94	0.94	0.95	0.93
Arou	Vegetation, growing	MRE (%)	−15.4	−6.9	2.2	−3.9	0.79
		RMSE (mm/day)	0.61	0.54	0.60	0.29	0.49
		R (-)	0.89	0.89	0.89	0.97	0.98
	Vegetation, dormant	MRE (%)	−45.1	−39.6	−37.6	−5.9	−18.7
		RMSE (mm/day)	0.17	0.15	0.14	0.07	0.11
		R (-)	0.98	0.98	0.98	0.98	0.97
	Overall	MRE (%)	−14.1	−5.5	3.2	−3.4	0.9
		RMSE (mm/day)	0.64	0.65	0.70	0.21	0.39
		R (-)	0.89	0.89	0.89	0.99	0.96

Among the five upscaling methods, SolRad performed the best during the season of dormant vegetation, and ConET_rF produced the best results during the season of active vegetation growth. Similar findings were obtained by Colaizzi et al. [22], which proved that the ConET_rF method performs better during periods of active vegetation growth. Thus, the optimum method is a combination of the SolRad method when vegetation growth is dormant, and the ConET_rF method during periods of active vegetation growth. RMSE values for the optimum method were 0.58 mm/day and 0.25 mm/day for the Guantao and Arou sites, respectively. The corresponding MRE values were −3.7% and 0%. The ET estimates from the optimum method are plotted in Figure 4. Figure 4 shows that the ET estimates from the optimum method were close to observational values.

Figure 4. The comparisons of clear-sky ET estimates from the optimum method with in situ observations.

Figure 4. The comparisons of clear-sky ET estimates from the optimum method with in situ observations.

The continuous daily ET was reconstructed with discrete accumulated clear-sky daily ET values measured by EC instruments at the Guantao and Arou sites. The cloudy-day ET temporal reconstruction results from the developed data assimilation methods (ET_EnKF and ET_SCE_UA), ET_SR, ET_ET_rF, and ET_HANTS method, were compared with the ET observations (see Figure 5). As shown, the data assimilation methods and the ET_ET_rF method performed the best, as indicated by the scattering about the 1:1 line. The ET_SR method underestimated the cloudy-day ET significantly, and the ET_HANTS method overestimated.

Figure 5. The comparisons of cloudy-day ET between the model estimations and observations. (Guantao site result (top); Arou site result (bottom)).

Figure 5. The comparisons of cloudy-day ET between the model estimations and observations. (Guantao site result (top); Arou site result (bottom)).

The statistics of each temporal reconstruction method are shown in Table 5. Over the whole time period, the RMSEs of most of the temporal reconstruction approaches were within 1.0 mm/day, except for the ET_HANTS method at the Guantao site (1.06 mm/day). The ET_ET_rF method produced the lowest RMSE values, which were 0.61 and 0.37 mm/day at Guantao and Arou, respectively. MRE values of ET_SCE_UA method were around zero (−5.4% and −2.6% at the Guantao and Arou sites), indicating the lowest ET estimation bias. Of the two data assimilation methods, the ET_SCE_UA approach outperformed the ET_EnKF. As indicated in Section 2.3.1, the ET_SCE_UA approach can make full use of four neighbored clear-sky ET observations to update the model parameters, while the ET_EnKF method only uses one previous ET observations. Thus, the ET_SCE_UA approach can produce lower errors than ET_EnKF method. The performances of these methods were also analyzed during periods of vegetation growth and dormant vegetation, as shown in Table 5. The periods of vegetation growth usually had larger RMSEs than periods of vegetation dormancy, and over the entire period, because the ET value was larger in summer than the other seasons, which was caused by higher solar radiation. The MRE values of most methods were closer to zero during periods of vegetation growth than periods of dormant vegetation, except for the ET_HANTS method. The data assimilation, ET_SR, and ET_ET_rF methods were all based on Penman-Monteith equation. Since the key parameter (surface resistance) of the Penman-Monteith equation is parameterized by leaf area index, it may get more accurate results for periods of vegetation growth than periods of vegetation dormancy [25].

Table 5. The statistics of temporal reconstruction approaches at the Guantao and Arou sites.

Site	Season	Statistics	ET_EnKF	ET_SCE_UA	ET_SR	ET_ETrF	ET_HANTS
Guantao	Vegetation, growing	MRE (%)	−7.4	−0.2	−13.8	9.6	41.1
		RMSE (mm/day)	1.14	0.81	1.20	0.80	1.53
		R (-)	0.75	0.84	0.73	0.83	0.40
	Vegetation, dormant	MRE (%)	−53.8	−22.4	−56.1	14.0	−28.0
		RMSE (mm/day)	0.43	0.51	0.44	0.41	0.44
		R (-)	0.85	0.79	0.85	0.82	0.84
	Overall	MRE (%)	−18.3	−5.4	−23.6	10.7	24.9
		RMSE (mm/day)	0.81	0.66	0.85	0.61	1.06
		R (-)	0.87	0.91	0.85	0.91	0.83
Arou	Vegetation, growing	MRE (%)	12.7	0.0	−19.7	−6.5	21.5
		RMSE (mm/day)	0.60	0.48	0.80	0.45	0.98
		R (-)	0.90	0.94	0.88	0.90	0.59
	Vegetation, dormant	MRE (%)	−3.2	−31.8	−73.5	−19.1	−39.2
		RMSE (mm/day)	0.51	0.33	0.52	0.30	0.30
		R (-)	0.57	0.63	0.51	0.68	0.76
	Overall	MRE (%)	9.6	−2.6	−25.4	10.5	12.9
		RMSE (mm/day)	0.57	0.41	0.66	0.37	0.73
		R (-)	0.93	0.96	0.91	0.95	0.88

Table 5. The statistics of temporal reconstruction approaches at the Guantao and Arou sites.

Site	Season	Statistics	ET_EnKF	ET_SCE_UA	ET_SR	ET_ETrF	ET_HANTS
Guantao	Vegetation, growing	MRE (%)	−7.4	−0.2	−13.8	9.6	41.1
		RMSE (mm/day)	1.14	0.81	1.20	0.80	1.53
		R (-)	0.75	0.84	0.73	0.83	0.40
	Vegetation, dormant	MRE (%)	−53.8	−22.4	−56.1	14.0	−28.0
		RMSE (mm/day)	0.43	0.51	0.44	0.41	0.44
		R (-)	0.85	0.79	0.85	0.82	0.84
	Overall	MRE (%)	−18.3	−5.4	−23.6	10.7	24.9
		RMSE (mm/day)	0.81	0.66	0.85	0.61	1.06
		R (-)	0.87	0.91	0.85	0.91	0.83
Arou	Vegetation, growing	MRE (%)	12.7	0.0	−19.7	−6.5	21.5
		RMSE (mm/day)	0.60	0.48	0.80	0.45	0.98
		R (-)	0.90	0.94	0.88	0.90	0.59
	Vegetation, dormant	MRE (%)	−3.2	−31.8	−73.5	−19.1	−39.2
		RMSE (mm/day)	0.51	0.33	0.52	0.30	0.30
		R (-)	0.57	0.63	0.51	0.68	0.76
	Overall	MRE (%)	9.6	−2.6	−25.4	10.5	12.9
		RMSE (mm/day)	0.57	0.41	0.66	0.37	0.73
		R (-)	0.93	0.96	0.91	0.95	0.88

No soil moisture or precipitation data were used in any of the temporal reconstruction approaches. However, soil moisture may increase suddenly due to precipitation and thermal remote sensing is not able to detect this change during rainy days. The variations of soil moisture will lead to change in ET distributions. Thus, errors may be cut down by developing ET temporal reconstruction methods that can assimilate microwave soil moisture (e.g., SMAP), or precipitation observations. For the ET_HANTS method, the cloudy-sky ET can be estimated based on clear-sky ET observations. However, the clear-sky ET is usually higher than cloudy-sky ET, which is a result of larger values of solar radiation. Without an input for forcing data (e.g., solar radiation, precipitation), the ET_HANTS method cannot detect a sudden decrease in ET caused by clouds. Thus, the estimated cloudy-sky ET from ET_HANTS is usually larger than in situ measurements, as was found here and is indicated by Figure 5. Moreover, the ET_HANTS method produced smooth seasonal variation of cloudy-sky daily ET. Thus, the scattering from the ET_HANTS method was sparser than the others (see Figure 5). The surface resistance in the ET_SR method with cloudy-sky conditions was obtained with the clear-sky values and ancillary parameters (LAI, m(T_min) and m(VPD) functions). Because threshold values of m(T_min) and m(VPD) were constructed based on vegetation types globally [25], they may not suitable for all specific land cover and climatic condition and should be calibrated with local observations.

Figure 6. The monthly ET variations from five temporal reconstruction approaches.

Figure 6. The monthly ET variations from five temporal reconstruction approaches.

The monthly variation of ET at the Guantao and Arou sites is shown in Figure 6. From Figure 6, the five methods show similar trends. At the Guantao site, all the methods display dual-peak ET values, caused by different growth patterns of winter wheat and summer corn. The low ebb of the ET estimates is indicative of a crop rotation period, where the winter wheat was cut down and the summer corn was at the early stages of growth. The ET_HANTS method estimated an ET that was larger than the other methods for May and August, which are the months of peak ET for both winter wheat and summer corn. The ET_HANTS overestimated the ET by 24.9% (see Table 5), which was caused by the large error for these two months at the Guantao site. At the Arou site, all methods depict a single-peak ET trend throughout the year. The ET_HANTS overestimated the ET by 12.9%, which was caused by the large error from May to August.

4.2. Regional Application

In this study, the regional instantaneous clear-sky ET values were estimated by the remotely sensed ET model proposed by Liu et al. [5]. The instantaneous ET estimates were upscaled to daily values with the combination of the SolRad and ConET_rF methods (the “optimum method”). The ConEF method was also applied for comparison. Over land surface covered with vegetation (cropland, garden-plot, forest, and grassland in Figure 2), the SolRad and ConET_rF methods were used for ET upscaling, during periods of vegetation dormancy and growth, respectively. Over land surface without vegetation (urban, water, road, and bare land), the SolRad method was used for ET upscaling over the entire year. In Beijing, the season of vegetation growth is from day of year (DOY) 100–283, for garden-plot, forest, and grassland areas, but for cropland, the season of vegetation growth is from DOY 100–161 and 182–283, according to vegetation phenology. As discussed in Section 4.1, the ET_ET_rF and ET_SCE_UA methods performed the best out of the ET temporal reconstruction approaches. Thus, those approaches were selected to produce the continuous daily ET with upscaled clear-sky daily ET for the Beijing area. The ET_SR method was also applied in Beijing for comparison. The comparisons between model estimates and the LAS observations are shown in Figure 7 and Figure 8.

Figure 7 shows the temporal upscaling results of the instantaneous ET estimates with the optimum method. The model estimates were close to LAS daily ET observations and the scatter plots approximately form 1:1 lines. At the Daxing site, the MRE (RMSE) value was 5.4% (0.72 mm/day). The corresponding MRE (RMSE) for the Miyun site was 0.3% (0.54 mm/day). Thus, the optimum method showed minor bias and low RMSE values at the two sites, which indicates that the optimum method is a suitable thermal ET upscaling method.

Figure 7. The comparisons of clear-sky daily ET estimates from the optimum method with LAS observations.

Figure 7. The comparisons of clear-sky daily ET estimates from the optimum method with LAS observations.

The continuous daily ET was reconstructed based on the discrete clear-sky daily ET upscaling results. The comparisons between reconstructed cloudy-day estimates and LAS observations are shown in Figure 8. Compared with observations, the scatter plots from the ET_SCE_UA and ET_ET_rF methods are close to 1:1 lines. The ET_SR method overestimated the continuous daily ET compared with observational values. The statistics of the temporal reconstruction approaches are summarized in Table 6. The ET_ET_rF and ET_SCE_UA methods produced lower RMSEs and MREs, and showed higher correlations than the ET_SR method. The ET_ET_rF method performed the best at both the Daxing and Miyun sites. At the Miyun site, the ET_ET_rF method underestimated the continuous daily ET when the value was larger than 3 mm/day (Figure 8), which was caused by the underestimation of clear-sky ET, derived from the remotely sensed model and the upscaling method (Figure 7).

Figure 8. The comparisons of cloudy-day ET between model estimates and in situ observations. (Daxing site result (top row); Miyun site result (bottom row)).

Figure 8. The comparisons of cloudy-day ET between model estimates and in situ observations. (Daxing site result (top row); Miyun site result (bottom row)).

Table 6. The statistics of temporal reconstruction approaches at the Daxing and Miyun sites.

Site Name	Statistics	ET_SCE_UA	ET_ETrF	ET_SR
Daxing	MRE (%)	−14.43	−12.98	72.34
	RMSE (mm/day)	0.74	0.56	1.60
	R (-)	0.83	0.87	0.73
Miyun	MRE (%)	−23.98	−16.24	42.74
	RMSE (mm/day)	0.89	0.75	1.32
	R (-)	0.74	0.77	0.74

Table 6. The statistics of temporal reconstruction approaches at the Daxing and Miyun sites.

Site Name	Statistics	ET_SCE_UA	ET_ETrF	ET_SR
Daxing	MRE (%)	−14.43	−12.98	72.34
	RMSE (mm/day)	0.74	0.56	1.60
	R (-)	0.83	0.87	0.73
Miyun	MRE (%)	−23.98	−16.24	42.74
	RMSE (mm/day)	0.89	0.75	1.32
	R (-)	0.74	0.77	0.74

The comparisons of monthly accumulated ET estimates and LAS observations are shown in Figure 9. Results from the three methods had the same trends as the observations. The Daxing site exhibits weak dual-peak ET variation, and the Miyun site shows single-peak ET variation. At the Daxing site, ET estimates from the three methods were higher than the observations. At the Miyun site, ET estimates from the three methods were higher than the observations from Apr. to Jun., but lower than the observations in Aug. The ET_ET_rF method had the lowest RMSE values (8.2 mm/month and 12.3 mm/month), and MRE values from ET_SCE_UA were closer to zero (−7.84% and −3.56%) at the Daxing and Miyun sites, respectively.

Figure 9. The comparisons of monthly ET between model estimates and observations at the Daxing and Miyun sites.

Figure 9. The comparisons of monthly ET between model estimates and observations at the Daxing and Miyun sites.

The regional annual ET estimates from the three methods are shown in Figure 10, and the ET products from MOD16A3 and ETWatch are also shown, for comparison. Over the whole region, the estimations from the ET_ET_rF and ET_SCE_UA methods had the same trend as the ET products from ETWatch. In the downtown area of Beijing, ET values were smaller than in other areas, as expected, due to paved land surfaces. The cropland and forest distributed around the downtown area showed high annual ET values, caused by irrigation and vegetation transpiration. The largest annual ET values were found over open water in the Miyun reservoir area, north of Beijing. However, the ET_SR method performed poorly over some of the urban areas, which showed large ET values. The ET_SR method is also adopted by MOD16A3, but no values were obtained from urban and water areas over Beijing.

The annual ET values of different land cover types are shown in Table 7. The ET values over water, cropland, garden-plot, forest, and grassland were in the first class with larger ET from the ET_SCE_UA, ET_ET_rF, MOD16A3, and ETWatch methods. Over the water area, most of the available energy is evaporated to water vapor that leads to largest ET estimates. Irrigation of vegetated lands, e.g., croplands and garden-plots, made vegetation transpiration values much larger than values from other vegetated lands that were unirrigated, e.g., unirrigated grassland. Forested land is generally located in mountainous regions, which experience high levels of precipitation, but a large portion of the rainfall becomes runoff and does not add to ET. Thus, forested land tends to have less soil water, and lower ET values than cropland and garden-plot. The smallest ET values were found in the urban areas of Beijing, because most of the land is paved, not vegetated, and precipitation is routed as storm water into the city’s sewer system. Since the Penman-Monteith method is better suited for use with densely vegetated land cover, the ET_SR method did not perform well over urban areas. The corresponding values for the ET_SCE_UA, ET_ET_rF, ET_SR, and ETWatch methods were 347 mm, 356 mm, 580 mm, and 470 mm, respectively. The results from the ET_ET_rF and ET_SCE_UA methods were closer to the water balance derived value (387 mm), and the ET_SR and ETWatch methods had larger errors. Thus, annual ET from the ET_SCE_UA method was in the same order of magnitude as the water balance method, and ET estimates of different land covers were more reasonable than the ET_ET_rF method (ET from forest was larger than water).

Figure 10. Annual ET distributions for the Beijing area in 2009.

Figure 10. Annual ET distributions for the Beijing area in 2009.

Table 7. Annual ET over different land covers in 2009.

Land Cover Type	ET_SCE_UA	ET_ETrF	ET_SR	MOD16A3	ETWatch
cropland	405	446	604	295	483
garden plot	404	463	626	290	505
forest	334	488	534	428	549
grassland	362	338	499	390	529
urban	153	146	609	-	382
road	94	120	504	-	340
water	705	480	656	-	522
bare land	320	368	608	207	447

Table 7. Annual ET over different land covers in 2009.

Land Cover Type	ET_SCE_UA	ET_ETrF	ET_SR	MOD16A3	ETWatch
cropland	405	446	604	295	483
garden plot	404	463	626	290	505
forest	334	488	534	428	549
grassland	362	338	499	390	529
urban	153	146	609	-	382
road	94	120	504	-	340
water	705	480	656	-	522
bare land	320	368	608	207	447

The averaged ET values over vegetated land covers (cropland, garden-plot, forest and grassland) were 376 mm, 434 mm, 566 mm, 350 mm, and 516 mm for ET_SCE_UA, ET_ET_rF, ET_SR, MOD16A3, and ETWatch, respectively. ET estimates were the highest from the ET_SR and ETWatch, while ET from MOD16A3 was the smallest, which indicates that ET was underestimated over the Beijing area. Moreover, ET values from forest and grassland were larger than the values from cropland and garden-plots, which did not seem reasonable for MOD16A3 ET products.

5. Conclusions

Long-term regional evapotranspiration is important for water resources planning and management and water-saving agricultural. The development of remotely sensed ET models enabled the mapping of regional ET at the time of satellite overpass. However, thermal remote sensing data are often contaminated by clouds. Thus, temporal upscaling methods are needed to upscale the instantaneous ET into daily values. Moreover, continuous daily ET should be reconstructed with discrete clear-sky ET. In this study, five temporal upscaling methods were selected to extend the instantaneous ET to daily values, and five temporal reconstruction approaches were used to generate continuous daily ET values for one year. These methods were compared extensively at two sites in China, namely, Guantao and Arou, which were selected based on distinct climate conditions. For regional application, annual ET distributions over the Beijing area for 2009 were produced with these methods, and the total annual ET was compared with the water balance method. Additionally, ET products from MOD16A3 and ETWatch were obtained for regional comparisons. LAS observations from two sites within Beijing (Daxing and Miyun) were acquired for ET estimate validations.

The five temporal upscaling methods were compared with EC measurements at the Guantao and Arou sites. Comparisons showed that the SolRad and ConET_rF methods performed the best, as indicated by having produced the smallest RMSE values. The commonly used ConEF method underestimated the daily ET, due to the increasing EF in the early morning and late afternoon. The CorEF and DiEF methods attempted to correct the underestimation by adding 10% of the ET, and by using three factors that enhance daily ET estimation with instantaneous EF observation. Because reference ET can capture impacts from changes in wind and humidity, the ConET_rF method was more stable than the ConEF and its related methods (CorEF and DiEF methods). Moreover, the ConET_rF method produced the smallest RMSE values during the season of vegetation growth, where Penman-Monteith equation performs well. The SolRad method performed best during the season of vegetation dormancy, and for non-vegetated land cover. The combination of the SolRad and ConET_rF methods, each optimized with respect to land cover and time period, produced RMSE values of 0.58 mm/day and 0.25 mm/day, at the Guantao and Arou sites, respectively. The corresponding MRE values were −3.7% and 0%. Thus, the combination of the SolRad and ConET_rF methods (the “optimum method”) performed the best overall.

The five temporal reconstruction approaches were used to generate continuous daily ET values with discrete clear-sky daily ET at the Guantao and Arou sites. Through comparison, the reference evapotranspiration method (ET_ET_rF), and the data assimilation method of the SCE_UA algorithm (ET_SCE_UA), were found to perform the best at the two sites. The ET_ET_rF method produced the lowest RMSE values of 0.61 and 0.37 mm/day at Guantao and Arou, respectively. The MRE values of the ET_SCE_UA method were around zero (−5.4% and −2.6% at the Guantao and Arou sites), which turned out to be the lowest ET estimation biases. The surface resistance method (ET_SR) and the HANTS method (ET_HANTS) had large errors compared to observational values. The temporal reconstructed continuous daily ET values might have had larger errors than monthly ET values, because the random errors of continuous daily ET values cancelled each other out. Rainfall that occurred during cloudy-days increased the soil moisture, EF and ET_rF, which may have led to errors in reconstructing continuous daily ET with discrete clear-sky ET values. In future studies, soil moisture information should be included to mitigate the variations of wet land surface conditions.

For regional application, a remotely sensed ET model was selected to calculate the instantaneous regional ET based on MODIS data. The optimum method was used to upscale the instantaneous ET to daily values. The ET_SCE_UA and ET_ET_rF methods were used to reconstruct continuous daily ET based on clear-sky daily ET. The distribution of calculated annual ET over the Beijing area had a similar trend to the MOD16A3 and ETWatch products. In 2009, the annual ET values over the Beijing area for the ET_SCE_UA, ET_ET_rF, ET_SR, and ETWatch methods were 347 mm, 356 mm, 580 mm, and 470 mm, respectively. The results from the ET_ET_rF and ET_SCE_UA methods were closer to the water balance derived annual ET for Beijing (387 mm). The ET estimates from ET_SCE_UA were more reasonable than those from ET_ET_rF over different land covers.

Although ET temporal upscaling and reconstruction methods have been compared at two typical sites and applied over the Beijing area in China, they need to be validated over broader flux sites around the world to assess broader application.