Improving the Evapotranspiration Estimation under Cloudy Condition by Extending the Ts-VI Triangle Model

Abstract

Evapotranspiration (ET) of soil-vegetation system is the main process of the water and energy exchange between the atmosphere and the land surface. Spatio-temporal continuous ET is vitally important to agriculture and ecological applications. Surface temperature and vegetation index (Ts-VI) triangle ET model based on remote sensing land surface temperature (LST) is widely used to monitor the land surface ET. However, a large number of missing data caused by the presence of clouds always reduces the availability of the main parameter LST, thus making the remote sensing-based ET estimation unavailable. In this paper, a method to improve the availability of ET estimates from Ts-VI model is proposed. Firstly, continuous LST product of the time series is obtained using a reconstruction algorithm, and then, the reconstructed LST is applied to the estimate ET using the Ts-VI model. The validation in the Heihe River Basin from 2009 to 2011 showed that the availability of ET estimates is improved from 25 days per year (d/yr) to 141 d/yr. Compared with the in situ data, a very good performance of the estimated ET is found with RMSE 1.23 mm/day and R² 0.6257 at point scale and RMSE 0.32 mm/day and R² 0.8556 at regional scale. This will improve the understanding of the water and energy exchange between the atmosphere and the land surface, especially under cloudy conditions.

1. Introduction

Evapotranspiration (ET) is the main process of water and energy exchange between the atmosphere and the land surface. At global scale, land surface ET usually accounts for 48–88% of net radiation [1] and about 60% of the total precipitation [2]. ET at regional scale plays an important role in hydrology, meteorology, and agricultural science and also in drought monitoring, water resource management, and crop yield estimation [3,4,5]. Hence, accurate estimates of ET with high spatial and temporal resolution (such as 0.01° spatial resolution and daily temporal resolution) have important theoretical and practical significance.

Remote sensing data for estimating ET mainly depends on visible-near infrared products, including surface albedo, biophysical parameters (leaf area index, fractional vegetation coverage, etc.), land cover, land surface temperature (LST), and emissivity. These products are usually used as the important input of the ET models. There have been a number of classic remote sensing ET retrieval algorithms over these years, including simplified empirical method [4], the surface temperature-vegetation index (Ts-VI) triangle method [6], the surface energy balance system (SEBS) [7], and data assimilation techniques [8]. In more recent years, new ET models have been developed, including the Atmosphere-Land Exchange Inverse model (ALEXI/DisALEXI) [9], Operational Simplified Surface Energy Balance model (SSEBOP) [10], and Satellite Irrigation Management Support (SIMS) [11]. They have been widely used to estimate regional or global ET derived by the remote sensing data [12]. However, as these ET models are derived by LST, the continuity of the ET estimates relies heavily on the availability of the LST data, which may be heavily influenced by the cloudiness [13,14].

The surface temperature and vegetation index (Ts-VI) triangle method [6,15] has been widely applied to derive ET using satellite data with promising accuracy over clear-sky conditions. Ts-VI model only needs satellite-derived land surface variables, especially suitable for the application scenario of irrigation area which may have dry and wet edges. However, the ET estimates of the Ts-VI triangle model is usually discontinuous either temporally or spatially, due to lacking of valid LST data causing by the cloudiness. Obtaining all-sky ET at the regional scale has become a great challenge for the Ts-VI triangle model.

There are several algorithms that have been proposed to fill the LST gaps, thus making all-sky ET estimation possible [16]. In theory, the filled LST can be implemented to obtain spatially complete ET distribution over large regions. However, the reliability of using the filled LST data to retrieve ET still requires further investigation. A major reason is that the LST filled by spatio-temporal reconstruction methods can only provide hypothetical clear-sky LST values but not the real LST under cloudy conditions [17,18]. In this case, it is very important to distinguish whether the satellite passes through a sunny day covered by clouds or a really cloudy day (or even a rainy day). To improve the availability of ET estimates, sunny days covered by clouds is the only application scenario for retrieving the cloud-free LST by the LST reconstruction algorithm.

The objective of this paper is twofold: (1) to propose a new method to derive spatial and temporal continuous land surface ET based on the fusion of the Ts-VI triangle model and the LST reconstruction method and (2) to verify the proposed method against the in situ ET measurements.

2. Materials and Methods

2.1. Study Area

The study area is located in the middle reaches of the Heihe River Basin in Gansu Province, China. It has become an ideal experimental site for the study of regional hydrology and has a very important impact on the climate change of the whole China region and even the whole world [19]. The temperate climate is arid, and the total area is about 38,000 km² (38.5° N–40.0° N, 98.5°–102° E). Figure 1 shows the Landsat TM image of the Heihe River Basin (2010) of the study area, which reflects the current situation of the Heihe River Basin in 2010. The surface elevation in most areas is about 1200 to 1600 m. A mountain about 3000 m above sea level in the southwestern part of the study area. One eddy covariance (EC) instrument named Yingke station is located on a vegetated surface with short crops, as shown in Figure 1.

2.2. Datasets and Pre-Processing

The input satellite datasets include Moderate-resolution Imaging Spectroradiometer (MODIS) products including the daily LST/emissivity products (MOD11A1) and the 16-day cloud-free NDVI data based on MODIS NDVI (MOD13A2), and Global Land Surface Satellite System (GLASS) albedo product which is also derived from MODIS albedo data. Downward shortwave and longwave radiation from China Meteorological Forcing Dataset (CMFD) are used to calculate the net radiation. The MICLCover land cover was used to reconstruct MODIS LST.

The dataset of EC observations at the Yingke station and a dataset of ground truth of land surface ET at the satellite pixel in the Heihe River Basin will be applied in the validation. Bilinearly interpolation is used to transfer temporal resolution of albedo and NDVI from 16 or 8 days to daily. For NDVI, the bilinearly interpolation method basically meets the growth law of vegetation, which does not produce too much deviation. For albedo, the fluctuation is relatively small, which will not produce large error to the net radiation.

2.2.1. MODIS LST

The MOD11A1 (V6) data used in this study is the 1 km resolution surface temperature/emissivity product L3 product in the MODIS terrestrial standard product with a time resolution of 1 day, from 1 January 2009 to 31 December 2012.

2.2.2. Cloud-Free MODIS NDVI

The NDVI data was reconstructed using the HANTS algorithm to perform time series on the original MOD13A2 (1 km, 16 days) product [20], which eliminates gaps due to clouds. After reprojection and bilinear interpolation, continuous daily NDVI with a resolution of 0.01° was obtained.

2.2.3. GLASS Albedo

This study uses the time-space continuous global land surface albedo dataset from GLASS provided by Beijing Normal University. This data uses direct inversion algorithm [21] in combination with the statistical based temporal filtering (STF) algorithm [22] on MODIS surface reflectivity product. The data set has a spatial resolution of 1 km and a temporal resolution of 8 days. After reprojection and bilinear interpolation, continuous daily albedo with resolution of 0.01° was obtained.

2.2.4. MICLCover Land Cover Map

The land cover map of the Heihe River Basin is a subset of China 0.01° land cover map (MICLCover) that incorporates multi-source local information [23], with the IGBP land cover classification system. Based on evidence theory, it fused the 1:100,000 land use data of China in 2000, the vegetation type of China vegetation atlas (1:1 million), the 1:100,000 glacier map of China, the 1:1 million swamp wet map, and MODIS 2001 land cover product (MOD12Q1). Verification results of MICLCover show that the overall consistency with China land use map is 88.84%.

2.2.5. China Meteorological Forcing Dataset

The CMF dataset is a set of near-surface meteorological and environmental element reanalysis data [24,25]. The dataset is based on the internationally available Princeton reanalysis data, GLDAS data, GEWEX-SRB radiation data, and TRMM precipitation data, combined with conventional meteorological observation data from the China Meteorological Administration. The dataset includes seven elements with a spatio-temporal resolution of 0.1° and 3 h. The elements are near-surface air temperature, near-surface air pressure, near-surface air specific humidity, near-surface full wind speed, downward shortwave radiation, downward longwave radiation, and ground precipitation rate.

The model in this paper requires daily-scale, spatially-resolution 0.01° downward shortwave radiation, and downward longwave radiation data, so it is necessary to downscale the meteorological element data from 0.1° to 0.01°. First, hourly data for each meteorological element is averaged to obtain the daily mean value, and then the daily value is downscaled using a semi-empirical statistical downscaling method [26,27,28].

2.2.6. In Situ Measurements

According to Li et al. [29] which described the experiment in the study area, a total of six flux towers plus other ET stations were set up in the study area. However, for the above reasons, Huazhaizi, Zhangye, and Zhangye City stations, which located in the non-vegetated area; and Arou station, which located above the elevation of 3000 m, are not appropriate for validation. However, there are no available ways to obtain the in situ data of Linze and Linze grassland stations. Therefore, the results are validated only at Yingke station. The data set contains EC flux data from the Oasis station of Yingke irrigation area from 27 December 2007 to 31 December 2011. The site is located in the farmland of Yingke Irrigation District in Zhangye City, Gansu Province. This paper uses the EC data from 1 January 2009 to 31 December 2011 as validation data. For more information, please refer to Li et al. [29].

2.2.7. Ground Truth of Land Surface Evapotranspiration at the Satellite Pixel

Another validation dataset is the ground truth of the land surface ET at the satellite pixel in the Heihe River Basin (http://data.tpdc.ac.cn/en/data/0420dfce-2c87-4d11-b42a-1a9f26414368/, accessed on 24 May 2020), which is based on the observation data of 11 stations in the flux observation matrix of the “Multiscale Observation Experiment of Non-uniform Underlying Surface ET” in the middle reaches of the Heihe River comprehensive observation network from June to September 2012. A comprehensive method of 6 scale expansion methods was used to produce this product (named midstream flux observation matrix area MODIS satellite transit transient/Relative true value of pixel scale of surface ET (1 km, daily)), considering the impact of surface heterogeneity on the ET scale expansion. For more information, please refer to [30,31].

2.3. Methods

To apply the Ts–VI triangle model under cloudy condition for reconstructed LST, several steps are needed as shown in Figure 2. For convenience, in this study, the defective LST to be reconstructed is defined as the target LST, and the other auxiliary LST is referred to the filled LST. The improved Ts-VI triangle model consists of 2 pivotal parts, which will be introduced in detail below.

2.3.1. LST Reconstruction Method

According to the LST filling method proposed by Zeng et al. [16], a simplified LST spatio-temporal reconstruction method based on surface features is introduced. To improve the usability of the MODIS LST, Zeng et al. [16] proposes a reconstruction method based on multi-temporal data. To reduce run time of the program and maintain the excellent performance, the classification step is removed, and the reconstruction algorithm operates based on the land cover map, which is the main difference to the previous works.

(1) Multi-temporal reconstruction of LST

LST reconstruction is based on the assumption of a linear change in LST, where similar changes occur on adjacent pixels with similar surface features. Let D_T and D_F denote the detected values which are absent of the target LST and the filled LST; then, the linear relationship between them can be expressed as

(1)$D_T=a\times D_F+b$

where a and b are regression coefficients. To avoid the effects of abnormal values, the linear relationship is calculated by a robust regression using iteratively reweighed least squares, as described by O’Leary [32]. Once the final coefficients are obtained, the detected values which are absent of the target LST can be calculated through the linear relationship.

In order to guarantee the linear relationship between the target and the filled LST, the acquisition time of the filled LST must be as close as possible to that of the target LST and have the supplemental information needed for reconstruction. Hence, a procedure searching for the filled LST and similar pixel pairs can be conducted.

It is reasonable to infer that a target pixel with LST value lower than 220 K can be considered as blocked or affected by the cloud which should be reconstructed by the proposed algorithm. The algorithm can be described as follows: (1) load the LST images that are 20 days ahead of the target image and detect the valid LST pixel among the pixels that correspond to the target pixel in position, with the value greater than 220 K, and with the acquisition time closest to the target pixel; (2) if the valid LST pixel is detected, the pixel is regarded as a filling pixel whose value can be used to replace the value of the target pixel, otherwise, finish the pixel reconstruction process for the target pixel. (3) repeat the previous two steps for every target pixel, and then, a filled image can be obtained.

With the filled images, a size variable window can be used to define the searching range of a target pixel, finding points with similar land cover type in the filled image should be done. To guarantee the robustness of the regression, the minimum number of similar points is set as 20. In order to look for the nearest similarity point, a 5 × 5 window centered on the target pixel is constructed and points with similar land cover are searched in the window. Similar points that are available are defined as points whose value is greater than 220 K both in the target and the filled images. If there are less than 20 similar points in the window, expand the window by 2 pixels and continue searching until the number of similar points in the window is no less than 20. If the window is expanded to the entire image and the number of available similar points is still less than 20, the pixel reconstruction fails.

(2) Post processing of robust regression

After the above steps, most of the missing values in LST images can be reconstructed, except for some abnormal pixels. Abnormal pixels are outliers whose LST values are much higher than their surrounding pixels, because of the improper selection of similar pixels for linear regression. To eliminate the abnormal pixels, an image-based histogram analysis is used as an outlier detector through calculating the quartiles of the image LST value and then setting them as two boundaries, as described by Zeng et al. [16]. The upper and lower boundary equations are

(2)$B_L=Q_1-1.5\times (Q_3-Q_1)$

(3)$B_H=Q_3+1.5\times (Q_3-Q_1)$

where B_L and B_H represent the upper and lower bounds of the data, and Q₁ and Q₃ represent the first and third quartiles of the data, respectively. Pixels with values out of bounds will be treated as outliers. Since the range of values on the target LST is large, the original LST is divided into 50 × 50 pixel blocks, and post-processing is performed in each block separately. A 3 × 3 pixel mean filtering is implemented to remove outliers. Through the above steps, invalid pixels can be reconstructed one by one.

2.3.2. Ts-VI ET Model

Based on equilibrium evaporation, Priestley and Taylor [33] established an ET calculation method for saturated underlying surface under the assumption of no advection, and it is known as the Priestley–Taylor (P-T) model:

(4)$LE=\alpha \frac{\Delta }{\Delta +\gamma }(R_n-G)$

where LE is latent heat flux (W/m²); α is the Priestley–Taylor coefficient; R_n is net radiation (W/m²); G is soil heat flux (W/m²); Δ is saturated water vapor pressure slope (Pa/K); and γ is the wet and dry table constant. Priestley and Taylor analyzed the observations of multiple ocean and wet land sites and concluded that α = 1.26 is the most reasonable under saturated underlying conditions.

The LE obtained for the whole pixel is expressed in W/m². The corresponding ET, expressed in mm/d, are converted. Latent heat of vaporization is calculated through daily average surface temperature.

The combination of P-T equation and Ts-VI triangle space method [6,12] is a widely used method for estimation of ET with remote sensing data. α is estimated using the Ts-VI triangle space, and then, the Evaporation fraction (EF) can be directly estimated from

(5)$EF=\alpha \frac{\Delta }{\Delta +\gamma }$

The EF calculated is instantaneous. To derive a daily EF, the instantaneous EF was assumed to be constant during the daylight hours, as the EF is relatively stable in the daytime [34].

Hence, the calculation of daily ET includes two parts. The first part is the estimates of EF, which can be calculated after the determination for the Ts-VI triangle space. Here, an automatic edges determination algorithm proposed by Tang’s [12] is introduced to estimate EF. The other is the estimation of daily (R_n-G), which can be calculated by the surface energy balance equation generated by meteorological forcing dataset.

(6)$ET=\frac{EF·{(R_n-G)}_{daily}}{\lambda }$

where, $\lambda$ is the latent heat of vaporization is calculated using daily average air temperature.

(1) Automatic edges determination algorithm for the Ts–VI triangle model

Considering that NDVI is only a surface greenness parameter [35], in this paper, the fractional vegetation cover (F_r) is used in Ts-VI triangle space construction instead of NDVI because it can represent the relative proportion between soil and vegetation in the pixel, but NDVI cannot [36]. F_r is therefore estimated from NDVI product pixel by pixel using the formula proposed by Carlson and Ripley [37]:

(7)$F_r=(\frac{NDVI-NDV{I}_{\min }}{NDV{I}_{\max }-NDV{I}_{\min }})$

where NDVI_min and NDVI_max represent the minimum and maximum NDVI values that usually correspond to bare soil and fully vegetated surface, respectively. They are assigned 0.2 and 0.86, respectively, in this work just as Prihodko and Goward do [38].

In general, the dry and wet edges constitute the two physical boundaries of the Ts-F_r (Ts-VI) feature space. For a given F_r, the soil moisture increases from the dry edge to the wet edge, the surface temperature decreases from the maximum to the minimum, and EF also correspondingly decreases from the maximum to the minimum value. The EF on the dry edge varies between 0 (bare soil) and 1 (full vegetation cover surface with no water deficit in the root zone). Jiang and Islam [39] assumed that the evaporation ratio varies linearly with vegetation coverage from (EF_min = 0, α = 0, F_r = 0) to (EF_max = 1, α = 1.26, F_r = 1) on the dry edge. When the soil moisture in the root zone is in a deficit state, the surface temperature of the dry edge will become higher for the surface of a certain vegetation cover (partially or completely). For the determination of wet edges, previous studies [6,40] recommended the surface temperature of water bodies or well-irrigated farmland as the wet edge temperature with potential ET. In addition, for a certain F_r, Jiang et al. [41] theoretically confirmed the linear change of EF with surface temperature by using the overall transfer equation for estimating sensible heat flux and calculated the error range of the EF within the Ts-VI triangle in a semi-empirical way.

Assuming that the dry-surface temperature is a linear function of F_r and the wet edge is a constant, theoretically, once the highest surface temperature at F_r = 0 and the surface temperature covered by water or full vegetation are known, the dry and wet edges in the Ts-Fr (Ts-VI) feature space can be determined. However, in arid and semi-arid regions, these two types of surface types (water bodies and full vegetation cover) are difficult to identify from remote sensing sensors and may not even exist. Therefore, Tang et al. [12] developed a practical method to automatically determine dry and wet edges. Here, the flow chart in Tang’s paper is used in order to make readers better understand the automatic edges determination algorithm for the Ts–VI triangle model as shown in Figure 3. For more specific information about the automatic edge determination algorithm for the Ts–VI triangle, please refer to Tang et al. [12].

(2) Calculation of daily net radiation and post processing

At daily scale, the soil heat flux is negligible. The net radiation received by the land surface can be calculated from the radiation balance equation, using surface albedo, downward shortwave radiation (W/m²), downward longwave radiation (W/m²), daily average surface temperature, and the surface emissivity, the comprehensive emissivity for the soil-vegetation system can be calculated as follows:

(8)$\epsilon =(1-F_r)\cdot {\epsilon }_s+F_r\cdot {\epsilon }_v$

where ε_s is the emissivity of bare soil (ε_s = 0.96), and ε_v is the emissivity of vegetation (ε_v = 0.98).

The implicit assumptions in the triangle feature space include (1) soil and canopy have different effects on surface temperature. Therefore, Ts-VI is performed in the vegetation covered area; (2) changes in surface temperature in the triangle space is not due to differences cased by atmospheric conditions, only caused by changes in the available water of the soil. Moreover, the result of the LST reconstruction method is credible only on cloudy weather. Thus, when calculating ET, there are some limitation factors making pixels satisfying the following conditions that should be removed from the ET map: (1) pixels without vegetation cover, (2) pixels with snow or ice covered, (3) pixels with precipitation on that day, and (4) pixels with elevation changes out of the range of 1200–1600 m. After removing these outliers, the ET estimation is finished.

2.3.3. Performance Metrics

Quadratic performance metrics are often used to assess the accuracy of geophysical retrievals (satellite measurements) with respect to true fields [42]. Several widely used error metrics are used, including the correlation coefficient (R²), root mean square error (RMSE), relative deviation (Bias), and mean absolute error (MAE) to quantitatively assess the results. These metrics are defined as follows:

(9)$R^2={[\frac{{\displaystyle \sum _{i=1}^n(O_i-\overline{O}){\displaystyle \sum _{i=1}^n(E_i-\overline{E})}}}{\sqrt{{\displaystyle \sum _{i=1}^n{(O_i-\overline{O})}^2{\displaystyle \sum _{i=1}^n{(E_i-\overline{E})}^2}}}}]}^2$

(10)$RMSE=\sqrt{{\displaystyle \sum _{i=1}^n\frac{{(O_i-E_i)}^2}{n}}}$

(11)$MAE={\displaystyle \sum _{i=1}^n\frac{\left|O_i-E_i\right|}{n}}$

(12)$Bias=\frac{{\displaystyle \sum _{i=1}^n{E}_i}}{{\displaystyle \sum _{i=1}^n{O}_i}}-1$

where O_i is the observed ET, $\overline{O}$ is the average of the observed ET; E_i is the estimated ET, $\overline{E}$ is the average of the estimated ET and n is the number of observations.

3. Results

The algorithm described in Section 3 was applied to all the data acquired over the study area from 2009–2011. For cloudy days, the instantaneous LST was contaminated by cloud, making the in situ measured instantaneous EF cloud contaminated. Only the “clear-sky” instantaneous EF are obtained, for the LST is “clear-sky”. However, these instantaneous EF cannot be compared directly. In this study, following the assumption that the LST was only contaminated instantaneously for the satellite transit time, EF is converted from instantaneous to daily in order to compare it with the in situ measured daily ET. The validation of this research includes two parts, the LST reconstruction and the ET estimation.

3.1. Performance of LST Reconstruction

In this section, the validation results of the MODIS LST reconstruction algorithm will be presented. Using the LST reconstruction algorithm described in Section 3, most of the values occluded by the cloud can be reconstructed to obtain the hypothetical “clear-sky” LST.

The reconstructed “clear-sky” LST can be considered as the point occluded by the cloud, assuming in the cloudless state. Hence, the method of selecting the occlusion of the original clear-sky pixel value, re-constructing, and comparing the reconstructed value with the original clear-sky LST value is used as the verification method of the LST reconstruction algorithm. Choosing 2009, 100 clear sky points are randomly taken out (pixel values greater than 220 K) for each scene image, set them to 0, and reconstructed them using the LST reconstruction algorithm to reconstruct. The reconstructed LST are compared with the original LST, as shown in Figure 4.

Taking the image before and after reconstruction on the 209th day of 2009 as an example, results are shown in Figure 5. It can be seen intuitively that there was a large number of gaps (Figure 5a) before reconstruction. After reconstruction, almost all the gaps had been filled, and the spatial distribution was consistent with the original image. Consistency and good reconstruction effect proved that the method was feasible and can reconstruct LST with a large amount of cloud coverage.

The temporal fractions of LST during the 3 years of the original LST and the reconstructed LST are shown in Figure 6. It can be seen that the fractions of LST are low before reconstruction, with the maximum value less than 0.7 and the minimum value less than 0.1. It means most areas in the original LST are covered by clouds over half of the period. However, after reconstruction, the fractions of LST cover improved to 0.9 for most areas, and the minimum value resulted to be bigger than 0.7. This result shows that the reconstruction method can reconstruct the LST contaminated by clouds in almost all situations.

3.2. Overall Performance of ET Model

After the reconstruction, the MODIS LST occluded by the cloud can be reconstructed with a valid value, and the reliability of the method had been verified. Hence, the reconstructed LST can be used in the ET model to obtain ET estimates under cloudy condition. To show the difference between the Ts-VI before and after reconstruction in the ET model, the 209th day in 2009 was selected for display, as shown in Figure 7.

As can be seen from Figure 7, EF estimates using reconstructed LST were greatly improved than using original LST, since the wet and dry edge had been greatly improved. It can be seen from the two scatter plots that there is a large number of abnormal points in the original data, mainly due to the low value occluded partially by the cloud, which greatly affected the wet edge in the EF estimates. The existence of the abnormal point made the space a very large trapezoid. The minimum value of the original LST is only 238.480 K, so the final calculated wet edge only depended on it, which was severely small. Hence, it will underestimate EF, while the EF of the abnormal point covered by the cloud mentioned in the previous section will be too large. After the reconstruction process, the LST data quality had been significantly improved, which was benefited from the filtering method in pre-processing and post-processing. It can be seen from Figure 7b that the reconstructed space was a relatively concentrated approximate triangle. It can be seen that spatio-temporal reconstruction has the effect of improving quality of LST, and the subsequent ET estimates will benefit it.

Since the triangle model needs EF and Rn, EF and Rn validation are carried out before ET validation to show the reliability of the proposed method. The ground measured Rn and LE at Yingke station at the same time as satellite pass (10:30 a.m.) are used to calculate EF and validate the result. Figure 8 and Figure 9 show the scatters of Rn and EF estimates by the model against in situ measurements, respectively. For Rn, a very good agreement can be found with RMSE = 41.57W/m² and R² = 0.580. For EF, also a very good agreement can be found with RMSE = 0.1997 and R² = 0.6288. As the maximum EF is 1, the in situ data with EF bigger than 1 are deleted, and there are some days on which the EF are overestimated as 1.

Figure 10 showed a time series of ET estimates from Ts-VI model derived by original and reconstructed LST against in situ measured ET during 2009–2011. It can be seen that the ET estimates of the model had a good agreement with the in situ measurements. At the same time, the ET estimates driven by reconstructed LST can capture its changes and trends over long time series, while the ET estimates driven by original LST were very scattered, and its temporal variation cannot be captured. The temporal coverage was also improved from 75 days to 422 days, indicating that the LST reconstruction process can significantly increase the ET availability.

The scatters of the ET estimates of model derived by the reconstructed LST against in situ measurements were shown in Figure 11. A very good agreement can be found with RMSE = 1.23 mm/day and R² = 0.626. However, there were still some points with large deviation; the reason would be discussed in the next section.

In addition, the ET estimates was also compared with the dataset of ground truth of land surface ET from June to September 2012, as shown in Figure 12. A well agreement can be found with RMSE = 0.317 mm/day and R² = 0.856. Almost each point has a small error, the result is much better than the result shown in Figure 12, and reason will be discussed in the next section.

4. Discussion

4.1. Advantages of This Study

The most important contribution of this study is improving the availability of ET estimates from Ts-VI triangular model based on the LST reconstruction method. The results show that, compared with original LST derived Ts-VI triangle ET model, the estimated value of our model is significantly increased from 75 days in three years to 422 days. Meanwhile, the quality is acceptable and is consistent with the previous research results in the same area, which shows that other LST based models have the same temporal variation [43,44].

In addition, the method relies on the input of remote sensing and forcing data rather than ground data. This is another difference from the crop coefficient-based Penman–Monteith method, which obtains crop coefficient through field experiment and changes during vegetation development [45]. In addition, it is also different from the ground data driven machine learning method [46].

This study makes a key step in overcoming the issue of gaps on cloudy days to obtain spatio-temporal continued daily ET using LST-based ET model. When using other LST based ET models, this method can also be used. Hence, it provides a possible way to combine the LST reconstruction method with other LST-based ET models to improve the availability of ET estimates in the future.

4.2. Limitations of This Study

There are still some shortcomings to overcome for the algorithm. For the LST reconstruction method, most of the points that deviate from the trend line are the overestimation of the original LST by the reconstruction LST, and the more concentrated ones appear within the interval where the original temperature is lower (less than 280 K). It has been observed that these points with large errors tend to be at the edge or in the middle of the cloud in the original LST data, and the surrounding points are mostly covered by clouds. This may be due to the fact that there are some abnormal points with low LST on the image, and these pixels may be partially blocked by the cloud, resulting in a low value. Although the robust regression can avoid some abnormal situations, false reconstructed values still exist in the reconstruction results. Therefore, to further improve the proposed algorithm and achieve better reconstruction results, the data preprocessing method need to be improved, such as by introducing a corrosion algorithm to remove the boundary of the cloud and the points in the cloud, but it may also result in less available points and lower reconstruction accuracy.

Problems encountered in the determination of dry and wet edges mainly include the following: (1) LST reconstruction algorithm might introduce new errors to the determination of dry and wet edges as the reconstructed LST are not true value [18], and (2) in a short period after a heavy rainfall event, the observed dry edge in the Ts–VI triangle may not indicate pixels with minimum ET for the given Fr due to the non-zero surface soil water content. If the observed dry edge is regarded to represent the true dry edge in this case, uncertainty and even error will be introduced.

Due to the limitations of the proposed model, spatial complete ET map is not always obtained, even if a seamless LST image successfully reconstructed is used as input data (see Figure 13). The main reason is that, as mentioned in Section 3, after the removal of points which do not meet model requirements, it is hard to acquire a complete ET map at a large scale using this algorithm. However, comparing with the ET map derived by the original LST, the spatial coverage of the ET map is also obviously increased. Thus, through the LST reconstruction method, the spatial coverage of ET can be increased, and the ET value can be evaluated under cloudy conditions.

Moreover, the satellite data that was used as input is instantaneous observation, while daily ET estimates need to be obtained. In addition, the reconstructed LST can be seen as a 20-day synthesis, which is not sensitive to mutations. As can be seen in Figure 11, the days with large differences between the estimation and the observation, are always the time when the EC observation has a sudden great change. These sudden changes always occur because of irrigation since the station is located in the farmland. These time scale effects will greatly affect the accuracy of ET estimation.

More importantly, another effect that can cause uncertainty in validation results is the scale effect, which is a main error source in the remote sensing field. As the results were greatly improved from Figure 12 to Figure 13, with R² increased by 0.23 and RMSE reduced by 0.37 mm/day, it is clear that the scale effect has a big influence on the accuracy of the Ts-VI ET model. When the scale effect between the validation data and the model derived data is excluded to a certain extent, the performance of the model will get a huge improvement.

To reduce the uncertainty in the ET estimates from the satellite data, further work needs to be carried out. The geostationary satellite data may be used, such as FY-4, which is available every half an hour, so that the time scale effect can be reduced. Meanwhile, a more universal model, such as the TSEB dual source model [47], may be used to get a space continuous ET map. The new MODIS LST and emissivity product (MOD21), which has been published recently, uses a physically based algorithm to dynamically retrieve both the LST and Emissivity simultaneously for the three MODIS thermal infrared bands (29, 31 and 32) at a spatial resolution of 1 km. The MOD21 product addresses the cold bias of 3–5 K in the MOD11 heritage split-window products over arid and semi-arid regions; thus, a better result may come out using the MOD21 data in the future. In addition, for the LST reconstruction method, more geographic factors will be considered, such as altitude, to improve the LST reconstruction in areas with large fluctuations in altitude. An algorithm may also be developed to reconstruct real LST with long time series as well [18].

5. Conclusions

The large number of missing data caused by the presence of clouds always reduces the availability LST, thus making the LST based ET estimation unavailable. In this study, we developed a gap-filling algorithm using LST reconstruction method for the Ts-VI triangle model to improve the availability of daily ET in the Heihe River Basin. Firstly, continuous LST is obtained using a reconstruction algorithm, and then, the reconstructed LST is applied to the estimate daily ET using the Ts-VI model. We found the following: (1) The LST reconstruction method can well present a “clear sky” LST. (2) For the ET derived by the Ts-VI triangle model using the reconstructed LST, a good performance is found at point scale and a better performance at regional scale. (3) The number of ET estimations is increased significantly with the help of LST reconstruction method. From 2009 to 2011, the availability of ET estimates has been improved significantly. In summary, this paper proposed an ET estimation algorithm combining an LST reconstruction method and the Ts-VI triangle ET model. The algorithm can obtain ET estimates under all-sky conditions, which greatly enrich effective ET estimates in the long time series. The experiments demonstrated that the algorithm is effective and feasible, and the LST spatio-temporal reconstruction is helpful for ET estimation.

Author Contributions

Funding acquisition, Y.C.; Methodology, B.L. and Y.C.; Software, B.L. and X.G.; Supervision, Y.C.; Validation, B.L.; Writing—original draft, B.L.; Writing—review & editing, Y.C. and H.L. All authors have read and agreed to the published version of the manuscript.

Funding

This study was financially supported by the National Natural Science Foundation of China (Grant No. 41901348 and 41901346) and the Key R&D Program of the Ministry of Science and Technology, China (Grant No. 2018YFC1506504).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. A quick view of the study area and location of the eddy covariance (EC) station on the Landsat TM mosaic image of the Heihe River Basin (2010) (composited by bands 5, 4, and 3. Sources: National Tibetan Plateau Data Center, 2013).

Figure 2. Flow chart of the proposed algorithm to estimate evapotranspiration, where LST: land surface temperature; Ts-Fr: surface temperature-fractional vegetation cover; NDVI: normalized difference vegetation index; MICL: multi-source information-based China landcover.

Figure 3. Flow chart of the algorithm of the proposed dry edge determination in the Ts–VI triangle space, where Fr is the fraction of vegetation; Taver is the average value, and δ is the standard deviation of the N maximum surface temperatures (referring to Tang et al. [12]).

Figure 4. Scatter plot of reconstructed land surface temperature, where LST: land surface temperature.

Figure 5. MODIS LST image on 209th, 2009: (a) before reconstruction, (b) after reconstruction, where LST: land surface temperature.

Figure 6. Coverage of LST: (a) before reconstruction, (b) after reconstruction, where LST: land surface temperature.

Figure 7. Scatter plots of 209th in 2009 for: (a) after and (b) before reconstruction; wet and dry edges are shown as black lines, where LST: land surface temperature; Fr: fractional vegetation cover.

Figure 8. Scatters of the Rn against in situ measurements, where Rn: net radiation.

Figure 9. Scatters of the EF estimates against in situ measurements, where EF: Evaporation fraction, RS: data derived from remote sensing; EC: data observed from eddy covariance.

Figure 10. Comparisons of ET estimates from reconstructed LST derived Ts-VI model (ET_RS_after), original LST derived Ts-VI model (ET_RS_before), and measured ET, from 2009–2011. Triangles represent the original MODIS LST derived ET, where ET: evapotranspiration; RS: data derived from remote sensing; EC: data observed from eddy covariance; LST: land surface temperature.

Figure 11. Scatters of the ET estimates of model derived by the reconstructed LST against in situ measurements, where ET: evapotranspiration.

Figure 12. Scatters of the ET estimates of model derived by the reconstructed LST against ground truth of land surface ET, where ET: evapotranspiration.

Figure 13. ET map on 209th, 2009: (a) before reconstruction, (b) after reconstruction, where ET: evapotranspiration.