Evapotranspiration Estimation Using Surface Energy Balance System Model: A Case Study in the Nagqu River Basin

Abstract

Calculation of actual evapotranspiration (AET) is of vital importance for the study of climate change, ecosystem carbon cycling, flooding, drought, and agricultural water demand. It is one of the more important components in the hydrological cycle and surface energy balance (SEB). How to accurately estimate AET especially for the Tibetan Plateau (TP) with complex terrain remains a challenge for the scientific community. Using multi-sensor remote sensing data, meteorological forcing data, and field observations, AET was derived for the Nagqu river basin of the Northern Tibetan Plateau from a surface energy balance system (SEBS) model. As inputs for SEBS, improved algorithms and datasets for land surface albedo and a cloud-free normalized difference vegetation index (NDVI) were also constructed. The model-estimated AET were compared with results by using the combinatory method (CM). The validation indicated that the model estimates of AET agreed well with the correlation coefficient, the root mean square error, and the mean percentage error of 0.972, 0.052 mm/h, and −10.4%, respectively. The comparison between SEBS estimation and CM results also proved the feasibility of parameterization schemes for land surface parameters and AET.

1. Introduction

Actual Evapotranspiration (AET), which is the water lost to the atmosphere from the land surface [1], is a crucial component of the terrestrial water cycle and energy balance [2,3]. In addition, accurate AET information at a regional scale is critical for quantitative understanding of land-atmosphere interactions and providing valuable means to efficiently use water resources [4,5]. However, it is difficult to measure AET through point observation in remote mountain regions, such as the Tibetan Plateau (TP), which also has harsh climate conditions. The TP, with an average elevation of more than 4000 m above sea level, is often called the ‘Third Pole’ because of excessive amounts of snow and ice existing, especially for mountainous regions. The large heat sink created during the spring and summer has a profound impact on the mid-troposphere [6,7]. The thermal forcing effectively enhances the Asian monsoon and modulates its variability [8]. In addition, for the water cycle, the TP contains the headwaters of major rivers in Asia and is often referred to as the ‘Asian Water Tower’ [9].

Several eddy-flux towers have been set up in the TP [10], but they cannot provide sufficient or accurate estimations of regional AET over a heterogeneous area. With rapid development of remote sensing technology, regional and meso-scale patterns of the characteristics of land surface parameters used for AET estimation can be mapped efficiently with satellite data [11]. Therefore, several satellite-based AET approaches are applicable for determining AET at different scales, which include simplified empirical regressions and physically based surface energy balance (SEB) models [12,13,14,15,16]. The empirical methods can be established based on in situ measurements of AET and related hydrometeorological parameters for different homogeneous regions. This kind of method has some advantages due to its simplicity. The empirical models are more accurate under the in-situ conditions it is based on. If these conditions are fairly homogenous over large scales, empirical models can be accurate enough. For heterogeneous regions, empirical models have their limitations. The physically-based models use meteorological parameters and land surface parameters retrieved from satellite images as model inputs. The calculation for these models include complex solution of the turbulent heat fluxes. However, these models have the ability to estimate meso-scale, regional scale, and even global scale AET with reasonable accuracy. They have other advantages, such as low cost, continuous monitoring, and clear physical mechanisms involved. The physically-based surface energy balance models have been widely used to estimate AET with the aid of remote sensing techniques [14,15,17,18].

Although researches on the estimation of AET have become available for various applications [11,19], most of them only focus on the regional scale or the global scale. AET is still difficult to calculate accurately because of the heterogeneous landscape and the large number of controlling factors involved. Studies on the estimation of AET at a basin scale over the TP are rather scarce. Furthermore, in recent inter-comparisons of global AET estimation results, large model uncertainties were observed [20]. First, these uncertainties come from different models and approaches for calculating AET. Second, for a heterogeneous land surface like the TP, the roughness length for heat transfer in a model can vary with geometric and environmental variables by several orders of magnitude. Large uncertainties are caused by estimating evaporation for different land surface types [20]. Third, some of the model input forcing parameters need improvement. For example, most of the previous models used vegetation indices processed by the maximum value compositing (MVC) method while the effects of cloud cover were not removed completely [21].

Since basin-scale AET data are crucial to close the water cycle, more attention should be paid to local basin scale studies such as the Nagqu river basin over the TP. The motivation for the current study is that regional AET data over the TP with reasonable accuracy are still needed. This study aims to test the feasibility of the surface energy balance system (SEBS) model to derive the AET at a local basin scale. As input parameters for AET calculations, broadband albedo was parameterized with a new formula. Normalized difference vegetation index (NDVI) was reconstructed by cloud removal before being added to the model.

2. Materials

Taking the quality and continuity of in-situ measurements into account, this study mainly focuses on the year 2003. Daily meteorological data, surface meteorological forcing data, 10-day composite SPOT Vegetation (VGT) data, and daily Moderate Resolution Imaging Spectroradiometer (MODIS) data were used as an input for 10-day mean AET calculations.

2.1. Study Area

The Nagqu river basin, which is located in the southern part of the TP, China, has a size of about 8749 km² and a mean elevation of more than 4808 m above sea level (Figure 1). It is located between the Tanggula Mountains and the Nyainqentanglha Mountains [22]. The Nagqu river basin is part of the Nu River system. Gentle hills and mountains are scattered throughout the entire basin terrain. The dominant vegetation type is meadow. Strong solar radiation and low air temperature are two characteristic climate conditions for the study area. The annual precipitation is about 430.2 mm with monthly sums ranging from 3.2 mm (January) to 103.1 mm (July). A clear demarcation between the dry and wet season is identified due to the influence of the Asian monsoon. Approximately 80% of the total precipitation falls during the wet season (from June to August).

2.2. Satellite Data

The input satellite data are shown in

Table 1. Input datasets for calculation of actual evapotranspiration by the surface energy balance system model.

Data Type	Variable	Unit	Sensors	Platform	Resolution
	Names				Spatial	Temporal
Remote sensing data	LCT	-	MODIS	Terra	1km	-
	DEM	m	ASTER	Terra	30 m	-
	LST	K	MODIS	Terra	1 km	daily
	AOD	-	MODIS	Terra	$1°\times 1°$	daily
	NDVI	-	VGT	SPOT-5	1 km	10-day
	albedo	-	VGT	SPOT-5	1 km	10-day
ITPCAS forcing data	Air temperature	K	-	-	$0.1°\times 0.1°$	3-h
	Specific humidity	kg/kg	-	-	$0.1°\times 0.1°$	3-h
	Wind speed	m/s	-	-	$0.1°\times 0.1°$	3-h
	Air pressure	Pa	-	-	$0.1°\times 0.1°$	3-h

LCT: Land Cover Type; DEM: Digital Elevation Model; LST: Land Surface Temperature; AOD: Aerosol Optical Depth; NDVI: Normalized Difference Vegetation Index; ITPCAS: Institute of Tibetan Plateau Research, Chinese Academy of Sciences; MODIS: Moderate Resolution Imaging Spectroradiometer; ASTER: Advanced Spaceborne Thermal Emission and Reflection Radiometer; VGT: Vegetation; SPOT: Satellite Pour I’Observation de la Terre.

. The datasets include: (1) 1-km MODIS land surface temperature (LST) data (MOD11A1), (2) 1-km VGT NDVI, (3) 1-km VGT narrowband albedo, (4) 30-m ASTER digital elevation model (DEM), (5) $1°\times 1°$ MODIS aerosol optical depth (AOD) (MOD08D3), and (6) 1-km MODIS land cover type (LCT) data (MOD12Q1). Since clouds are always occurring in the satellite images, we only choose images taken with total cloud cover of less than 40% in this study. With ancillary information from MOD06_L2, images with cloud cover of more than 40% were excluded. To maintain consistency in spatial and temporal resolution, all satellite data were preprocessed by temporal-spatial matching. The spatial resolution was aggregated into 1 km cells using a nearest neighbor interpolation method and the daily temporal resolution was averaged into 10-day intervals.

2.3. In Situ Data and Surface Forcing Data

Taking data continuity and quality into account, the Northern Portable Automated Meso-net (NPAM) and Bu Jiao (BJ) stations in the Nagqu river basin, together with the D105 and ANNI stations around the study area were used. The observational data for 2003 from the Coordinated Enhanced Observing Period Asia–Australia Monsoon Project on the Tibetan Plateau (CAMP/Tibet) were used as inputs for the combinatory method (CM). The in situ meteorological data included hourly multiple-layer air temperature, relative humidity, wind speed, soil heat flux, air pressure, and soil moisture. Near surface air temperature, specific humidity, and wind speed of the ITPCAS forcing dataset were selected as inputs for the SEBS model. The forcing dataset has merged the observations from 740 operational stations (among which one station was located in the study watershed) of the China Meteorological Administration (CMA) with the corresponding Princeton global meteorological forcing dataset [23]. All data were preprocessed to keep the same resolution by temporal-spatial matching (1). The forcing data were spatially interpolated to a 1-km resolution using a nearest neighbor interpolation method (2). The daily LST data of the Terra/MODIS were averaged into 10-day intervals.

3. Methods

Land surface physical parameters, such as NDVI, albedo, land surface emissivity, and downward shortwave radiation play an important role in land-atmosphere interactions. In this study, land surface characteristics were derived from satellite images and meteorological forcing data. They were used as inputs for the SEBS model.

3.1. Reconstruction of Normalized Difference Vegetation Index

The MVC method has been widely used to derive vegetation products. However, some cloud and haze effects remain in NDVI data after using the MVC method. Therefore, the NDVI time series were reconstructed by the harmonic analysis time series (HANTS) algorithm, which combines cloud removal and data smoothing in a single operation [24]. The principle of the algorithm is that vegetation growth (NDVI) has a strong seasonal effect. It can be simulated by a series of low frequency sinusoidal functions. However, cloud cover can cause negative NDVI values and was previously considered to be high frequency ‘noise’ [24]. HANTS is based on a Fourier analysis and any pixels considered ‘cloudy’ were replaced with filtered NDVI time-series values. The HANTS algorithm can be written in the following form.

$$\mathrm{y}\left(t\right)=a_0+{\displaystyle \sum }_{j=1}^n{a}_j\cos \left({\omega }_{j}t+{\varnothing }_j\right)$$

(1)

where $\mathrm{y}\left(t\right)$ is the fitted curve value at time t. $a_0$ is the average value of the time series. $n$ is the number of harmonics. $a_j$ is the Fourier amplitude component. ${\omega }_j$ is the Fourier frequency. ${\varnothing }_j$ is the Fourier phase component. In HANTS, curve fitting is applied iteratively. Based on all data points for a pixel, a least squares curve is calculated. All data points are then compared to the least squares curve. Because of cloud cover, data points that are clearly below the curve are candidates for rejection. Data points that have the greatest negative deviation are removed. The absolute error in the negative direction of the remaining observations should be smaller than 0.05 with respect to a currently fitted curve. Based on the remaining data points, a new curve is computed and the process is repeated using the remaining points until no points exceed the pre-defined error threshold [25]. There is no specific rule for determining threshold or parameters in HANTS, but some experience is usually needed before getting optimal values.

3.2. Determination of Albedo and Emissivity

In this study, we present an improvement of the model that converts the narrowband to broadband albedo for the land surface based on Liang [26]. The equation (Equation (2)) was built to calculate land surface broadband albedo by using the fitting method. Since broadband albedo is a function of the narrowband albedo, the in situ broadband albedo from BJ, D105, NPAM, and narrowband albedo (${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3,$ and ${\alpha }_4$) from SPOT/VGT were used in the fitting procedure. Then the narrowband albedo was converted to broadband albedo for all land cover types except water. Since we do not have in situ albedo measurements in the study area for water, the broadband albedo of water was determined by the method (Equation (3)) proposed by Liang [26].

$${\alpha }_{short\_land}=0.1670+2.2305{\alpha }_2-0.7477{\alpha }_3-4.5082{\alpha }_1{\alpha }_4$$

(2)

$${\alpha }_{short\_water}=-0.0022+0.3512{\alpha }_1+0.1629{\alpha }_2+0.3415{\alpha }_3+0.1651{\alpha }_4$$

(3)

where ${\alpha }_{short\_land}$ and ${\alpha }_{short\_water}$ are the broadband albedo for land surface and water body, $\mathrm{and} {\alpha }_1$, ${\alpha }_2$, ${\alpha }_3,$ and ${\alpha }_4$ are the narrowband albedos of four spectral bands, which are 0.43–0.47, 0.61–0.68, 0.78–0.89, and 1.58–1.75 µm, respectively. The observational albedo at stations D105, NPAM, BJ, and ANNI was used to evaluate the conversion formula before and after the improvement. The root mean square error (RMSE) and mean percentage error (MPE) were used to quantitatively evaluate the model accuracy.

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(X_{obs,i}-X_{model,i}\right)}^2}{N}}$$

(4)

$$\mathrm{MPE}=\frac{100\%}{N}{\displaystyle \sum }_{i=1}^N\frac{X_{obs,i}-X_{model,i}}{X_{obs,i}}$$

(5)

where $X_{obs,i}$ stands for in-situ measurements. $X_{model,i}$ stands for model derived values and N is the sample number.

Based on derived NDVI and ${\alpha }_{short}$, emissivity was derived from the following equation.

$$\mathsf{\epsilon }=P_V·{\epsilon }_v+\left(P_C-P_V\right)·{\epsilon }_{water}+\left(1-P_C\right)·{\epsilon }_s$$

(6)

where $\mathsf{\epsilon }$ is the emissivity of mixed pixels, which is considered as a mixture of vegetation, water, and bare soil. ${\epsilon }_v$, ${\epsilon }_{water},$ and ${\epsilon }_s$ are the emissivity values for vegetation, water, and bare soil, respectively. The ${\epsilon }_v$, ${\epsilon }_{water}$, and ${\epsilon }_s$ are 0.98831, 0.994685, and 0.972785 [27,28]. $P_C$ is the ratio of mixed-pixel albedo to bare soil albedo. $P_V$ is the fractional vegetation cover. It can be calculated using the equation below.

$$P_V={\left[\frac{NDVI-NDV{I}_{min}}{NDV{I}_{max}-NDV{I}_{min}}\right]}^2$$

(7)

where $NDV{I}_{max}$ and $NDV{I}_{min}$ are the $NDVI$ values for full vegetation cover and bare soil, respectively [29].

3.3. Parameterization of Land Surface Heat Fluxes

The surface energy balance system (SEBS) is a single source model used to calculate the land surface heat fluxes. The detailed framework of SEBS can be found in Su [30]. A dynamic model for thermal roughness [31], the Monin-Obukhov similarity theory (MOST) and the bulk atmospheric similarity (BAS) theory [32] were integrated in the SEBS model. SEBS can be used at local and regional scales for all atmospheric stability conditions. AET is estimated as the residual of the SEB equation (Equation (8)).

$$R_n=H+\mathsf{\lambda }\mathrm{E}+\mathrm{G}$$

(8)

where $R_n$ ($W/m^2$)is the net radiation flux, H ($W/m^2$) and G ($W/m^2$) are the turbulent sensible heat flux and soil heat flux, respectively, $\mathsf{\lambda }\mathrm{E}$ ($W/m^2$) is the turbulent latent heat flux, and $\mathsf{\lambda }$ ($\mathrm{J}/\mathrm{kg}$) is the latent heat of vaporization.

The net radiation can be calculated using the equation below.

$$R_n=\left(1-\alpha \right)·R_{swd}+\epsilon ·R_{lwd}-\sigma ·\epsilon ·T_S^4$$

(9)

where $\alpha$ is the albedo, $T_S$ is the land surface temperature, $\epsilon$ is the surface emissivity, $\sigma$ is the Stefan-Bolzmann constant, $R_{swd}$ is the downward shortwave radiation, and $R_{lwd}$ is the downward longwave radiation.

The soil heat flux is calculated by the formula below.

$$G=R_n·\left[0.05+0.265·\left(1-P_V\right)\right]$$

(10)

The SEBS model constrains the surface heat flux estimation by introducing wet limit cases and dry limit cases. In wet limit cases, the evapotranspiration is only limited by the available energy. The SEB equation can be written by using the equation below.

$$\lambda E_{wet}=R_n-G-H_{wet}$$

(11)

In dry limit cases, the latent heat flux becomes zero because of soil water shortage. The SEB equation can be given by the equation below.

$$H_{dry}=R_n-G$$

(12)

The relative evaporation is shown below.

$${ʌ}_r=1-\frac{\lambda E_{wet}-\lambda E}{\lambda E_{wet}}$$

(13)

By the use of the above limit conditions, SEBS avoids the tedious process of selecting hot and cold pixels, which is different from other SEB models. The evaporative fraction is given by the formula below.

$$ʌ=\frac{\lambda E}{H+\lambda E}=\frac{\lambda E}{R_n-G}=\frac{{ʌ}_r·\lambda E_{wet}}{R_n-G}$$

(14)

where ${ʌ}_r$ is relative evaporation and $\lambda E_{wet}$ is the latent heat flux at the wet limit. The detailed computational schemes for $\lambda E_{wet}$ can be found in Reference [30].

The surface solar radiation is an important input in SEBS. For a mountainous region like the TP, $R_{swd}$ varies with different geometric relationships between the solar position and the tilt of land surface, which includes factors of altitude, surface slope, and aspect. $R_{swd}$ is divided into direct ($I_b$), diffuse ($I_d$), and reflected ($I_r$) radiation, which can be computed by Equations (15)–(17).

$$I_b=I_0·{\tau }_c·\cos \theta$$

(15)

$$I_d=I_0·{\tau }_d·{(\cos s)}^2/(2\sin a)$$

(16)

$$I_r=r·I_0·{\tau }_r·{(\sin s)}^2/(2\sin a)$$

(17)

where $I_0$ is the solar irradiance at the top of the atmosphere, $a$ is the solar altitude angle, $\theta$ is solar incidence angle, $s$ is the slope of the ground surface, and $r$ is the ground albedo. ${\tau }_c$, ${\tau }_d,$ and ${\tau }_r$ are the transmittance for solar beam radiation, the solar diffuse radiation, and the reflected radiation, respectively. The detailed solutions for the parameters in Equations (15)–(17) were presented in Reference [33].

3.4. The Combinatory Method

In the absence of measured latent and sensible heat flux data by eddy-covariance instruments, AET, which is determined by the combinatory method, were used for validation with estimation results obtained by SEBS. Since CM needs gradient observation of wind speed, air temperature, and humidity as inputs, only in situ measurements from four stations (NAPM, BJ, D105, and ANNI) were used to calculate AET at point scale. The combinatory method combines the aerodynamic method with the energy balance method and determines the turbulent fluxes in the surface layer based on MOST [5,34]. $\mathsf{\lambda }\mathrm{E}$ can be derived by using Equations (18)–(22).

$$H_0=\rho C_p{k}^2{Z}_A^2\frac{\partial U}{\partial Z}\frac{\partial \theta }{\partial Z}$$

(18)

$${\mathsf{\lambda }\mathrm{E}}_0=\rho \mathsf{\lambda }{k}^2{Z}_A^2\frac{\partial U}{\partial Z}\frac{\partial q}{\partial Z}$$

(19)

$$G=G_Z+{{\displaystyle \int }}_0^Z{C}_S\frac{\partial T}{\partial t}\mathrm{dz}$$

(20)

$$\mathrm{F}=\frac{R_n-G}{H_0+{\mathsf{\lambda }\mathrm{E}}_0}$$

(21)

$$\mathsf{\lambda }\mathrm{E}={\mathsf{\lambda }\mathrm{E}}_0\cdot \mathrm{F}$$

(22)

In this case, $H_0$ and ${\mathsf{\lambda }\mathrm{E}}_0$ $G\left(\mathrm{W}/m^2\right)$ are the turbulent sensible heat flux and turbulent latent heat flux before correcting. $G_Z\left(W/m^2\right)$ is the soil heat flux measured at a specific soil depth, $\mathrm{F}$ is the stratification influence function, $\mathsf{\lambda }\mathrm{E}$ ($W/m^2$) is the latent heat flux after correction, and $\mathrm{E} \left(\mathrm{mm}\right)$ is evapotranspiration. $C_p$($\mathrm{J}/\mathrm{kg}/\mathrm{K}$) is the specific heat capacity of air, $\rho$ ($kg/m^3$) is air density, and k = 0.4 (Von Karman’s constant). $Z_A=\sqrt{Z_i\cdot Z_{i+1}}$, $Z_i\left(m\right)$ and $Z_{i+1}\left(m\right)$ are the observation heights of two levels. $C_S$ is specific heat capacity of soil. T ($\mathrm{K}$) and U$\left(m/s\right)$ are the measured air temperatures and wind speed while $q$ and $\theta$ are specific humidity and potential temperature, respectively [5].

4. Results

4.1. Reconstructed Normalized Difference Vegetation Index

Two different pixels (located at 31.79° N, 91.94° E and 31.09° N, 91.71° E) were chosen randomly for comparison between the original NDVI processed by MVC and reconstructed NDVI by HANTS (Figure 2). The MVC process was used to recalibrate for the removal of cloud, but some abnormal values caused by residual cloud cover still remain visible in Figure 2a,d. As can be seen on DOY 233 (Figure 2a), the NDVI value of the pixel marked by a red point suddenly decreased from 0.45 to 0.1. This signal was recorded as a low NDVI value in blue in Figure 2b,e. The change took place in a very short time period of 10 days, which was not caused by a vegetation growth process but by clouds. After reconstruction (Figure 2c,f), some dark green areas diminished and were replaced by normal NDVI values. As shown by Figure 2a,d, an improved smoothing of time series can be achieved. Negative outliers of NDVI observations that are clearly below the curve due to cloud cover or atmospheric haze, were removed while the characteristics of the original trend was maintained.

4.2. Improvement of Broadband Albedo

As shown in Figure 3a,b, the improved method (Equation (2)) has higher R (0.655), lower RMSE (0.078), and MPE (−7.21%) than the original method (Equation (3)) (with R, RMSE, and MPE values of 0.485%, 0.1%, and 10.629%, respectively), which means the new conversion model could provide a better estimation of broadband albedo than the original one. Since snow-covered situations have not been taken into account, some errors might occur when albedo is relatively high. According to global land cover 2000 (GLC2000), the dominant land cover in the study area is alpine and subalpine plain grass. The error of estimated albedo could be large only when snowfall happens, which only takes up a small proportion. This can also be identified from Figure 3.

To further validate the improved model, independent observational data from the ANNI station were selected for cross validation (Figure 4). As shown by Figure 4a,b, R, MPE, and RMSE of the original method are 0.511, 27.5%, and 0.108 while the statistical indexes from the improved method are 0.606, 9.8%, and 0.076, respectively. Therefore, the improved model could provide a more accurate estimation of broadband albedo in the Nagqu river basin.

4.3. Estimation of Actual Evapotranspiration

Land surface 10-day AET derived from SEBS were validated with the 10-day AET determined using in situ observations taken in the Nagqu river basin (Figure 5). Since the validation sites are built on relatively flat and homogenous underlying surfaces, the observations are representative of a 1-km pixel. For the time being, how to remove cloud effects on LST and downward radiation flux remains a challenging issue. As shown in Figure 5, the model-estimated AET values are very close to CM results with R, RMSE, and MPE of 0.972, 0.052 mm/h, and −10.4%, respectively.

Maps of hourly AET in the Nagqu river basin in January, April, May, and December are shown in Figure 6. AET in January (0 to 0.29 mm/h, Figure 6a) and December (0 to 0.43 mm/h, Figure 6d) were significantly lower than those in May (varies from 0 to 0.77 mm/h, Figure 6c). Because of the increasing soil moisture during the period of frozen soil melt, AET increases suddenly from April to May, with average AET increasing from 0.184 mm/h to 0.479 mm/h. In the SEBS model, the uncertainty in AET is greatly reduced by taking the energy balance at the limiting cases into account. In comparison with previous studies [12,35], which derived the latent heat flux only through LST and surface layer meteorological conditions and are not constrained by limiting cases, the associated uncertainty in the derived ET is small. Moreover, instead of using fixed values as before, the dynamic roughness height for heat transfer was incorporated in SEBS.

5. Discussion

Because of complicated weather conditions over the TP, only several satellite images with total cloud cover less than 40% were used in this study. Although the derived results were generally in agreement with field measurements, some discrepancies still exist. Several reasons may account for this. First, there are errors in the satellite-derived parameters, such as LST, albedo, emissivity, and NDVI. The surface meteorological forcing data of air temperature, humidity, and wind speed also have errors, which has been proved by several studies [31,36]. Second, scale issues always exist between satellite-derived results and in-situ measurements [37]. Since there is no absolute homogenous underlying surface, the error becomes much larger with surface heterogeneity. Third, the SEBS model itself has assumptions and limitations. For example, the model is based on the Monin-Obukhov similarity theory, which produces larger error under stable atmospheric stratification [17,18,30].

The MVC method is widely used to remove cloud cover on satellite images [38,39,40]. It can be used to produce time series of land surface characteristic parameters, such as NDVI, LST, and the leaf area index (LAI) [21,25,41]. However, evidence has shown that artefacts caused by residual cloud cover and atmospheric haze remain visible in MVC-processed time series, which has also been proved in our study. The HANTS algorithm was used to reconstruct one-year cloud-free time series based on the original MVC time series data. The robustness of the HANTS algorithm was confirmed by validating at two different vegetation sites. The HANTS algorithm is superior to other time series processing methods but improvements related to its numerical stability in case of long gaps with no data are still desirable. Taking superiority of the HANTS algorithm into account, we would like to suggest using the above method to further process MVC time series to avoid some uncertainties especially for the estimation of AET and vegetation trend analyses.

The improved land surface characteristic parameter retrieval algorithms, which were incorporated in SEBS in this study, have the potential to be used to derive time series of AET and other land surface heat fluxes (such as sensible heat flux, radiation flux, and soil heat flux) as well. It can also be applied to regions at a larger spatial scale, such as the Tibetan Plateau, China, and even the globe. Additionally, the SEBS model itself can be applied to geostationary satellites to estimate AET with a fine temporal resolution [18]. The derived AET, together with other components of the water cycle, will help close the water cycle at a basin scale and promote the understanding of energy and the water cycle over the TP. Another good example is that the derived AET and in situ hydrometeorological parameters were used to investigate the inner relationship among AET and energy parameters, hydrological parameters, and dynamical parameters [5]. The results show that the energy-related factors are most important for the AET in areas with high elevation like the Nagqu river basin. Both the precipitation and wind speed have positive and negative feedback on the AET.

On the application side, SEBS has great potential. It can be suitable for an operational monitoring system if the input forcing data can be provided in a timely manner. SEBS results may be used for initialization or validation of hydrological, meteorological, and ecological models that usually require SEB components at different spatial and temporal scales. The second valuable application is soil water deficit monitoring [42,43]. The model can provide actual sensible heat flux as well as sensible heat flux under limiting cases. Although cloud contaminations have been removed successfully in NDVI, albedo, and emissivity data, the application of optical remote sensing data is still limited by cloud contamination because the cloud-contaminated LST and the cloud effects on radiation have not been taken into account yet. In the future, broadband albedo estimation would be improved by taking snow-covered situations into account. Moreover, the estimation of the cloud effects on LST and radiation should be focused on as well. Lastly, more extensive observation networks equipped with eddy covariance instruments need to be constructed over the TP, especially for the vast western region, to provide essential validation information for different satellite estimations and land surface models.

6. Conclusions

In this study, the land surface parameters in the Nagqu river basin, such as NDVI and surface albedo, were retrieved from satellite images. Based on meteorological forcing data and derived land surface parameters, the SEBS model were used to estimate AET with a fine spatial resolution of 1 km. The following conclusions were drawn.

(1) A cloud-free time series of NDVI was successfully reconstructed using the HANTS algorithm.

(2) A land surface broadband albedo parameterization scheme was improved. The albedos derived from the improved method were in good agreement with in situ measurements with a correlation coefficient, a root mean square error, and a mean percentage error of 0.606%, 0.076%, and −9.8%, respectively.

(3) AET with a fine spatial resolution of 1 km was derived in the Nagqu river basin based on the SEBS model with TERRA/MODIS, SPOT/VGT, ASTER/DEM data, and meteorological forcing data as model inputs.

(4) The feasibility and potentiality of using SEBS model to derive mesoscale AET has been tested by comparing it with results from CM in this paper. The validation results indicated that the model-estimated AET agreed well with CM AET with a correlation coefficient, a root mean square error, and a mean percentage error of 0.972, 0.052 mm/h, and −10.4%, respectively.