Regional Actual Evapotranspiration Estimation with Land and Meteorological Variables Derived from Multi-Source Satellite Data

Abstract

Evapotranspiration (ET) is one of the components in the water cycle and the surface energy balance systems. It is fundamental information for agriculture, water resource management, and climate change research. This study presents a scheme for regional actual evapotranspiration estimation using multi-source satellite data to compute key land and meteorological variables characterizing land surface, soil, vegetation, and the atmospheric boundary layer. The algorithms are validated using ground observations from the Heihe River Basin of northwest China. Monthly data estimates at a resolution of 1 km from the proposed algorithms compared well with ground observation data, with a root mean square error (RMSE) of 0.80 mm and a mean relative error (MRE) of −7.11%. The overall deviation between the average yearly ET derived from the proposed algorithms and ground-based water balance measurements was 9.44% for a small watershed and 1% for the entire basin. This study demonstrates that both accuracy and spatial depiction of actual evapotranspiration estimation can be significantly improved by using multi-source satellite data to measure the required land surface and meteorological variables. This reduces dependence on spatial interpolation of ground-derived meteorological variables which can be problematic, especially in data-sparse regions, and allows the production of region-wide ET datasets.

1. Introduction

Actual evapotranspiration (ET) is a major component of the water cycle and plays a key role in surface energy balance by linking the hydrological cycle, carbon budget, and energy transfer [1,2,3,4]. ET is the transfer of water from land surfaces to the atmosphere primarily through turbulence which depends on the aerodynamic characteristics of the surface as determined by the land surface cover type [5], as well as climatic conditions. Consequently, a number of remote sensing-based models for regional ET estimation incorporate both land [6,7,8,9] and meteorological variables [10,11,12,13]. Validation of existing ET datasets estimated from those models suggests varying performances over different landcovers, climate, and elevation zones [14,15,16,17], with the errors in large regional ET estimates often caused by land surface heterogeneity derived from terrain roughness, landcover, vegetation community, soil context, etc. [18].

Land variables such as landcover, soil moisture content, albedo, vegetation indices, and land surface temperature (LST) are usually derived from single sensor or satellite data sources such as Landsat Thematic Mapper (TM) [6,7,8] or the Moderate Resolution Imaging Spectroradiometer (MODIS) [10,11,19,20], or similar satellites [21]. Vegetation cover [22], topographical terrain variation [2], and surface roughness derived from radar backscatter [23] control the aerodynamic roughness length [4,24,25,26] which influences the exchange from land to the atmosphere and, consequently, the movement of heat, momentum, and atmospheric moisture [24,25]. But temporal changes in the vegetation canopy structure also need to be addressed [27].

Albedo, vegetation indices, and surface temperature can be used to quantify soil heat flux [10,11,28,29,30]; however, it is often assumed to be zero or ignored by many existing satellite-based models over an annual [6,20,31,32] or daily [10,33] period, or is considered to be a constant proportion of instantaneous net radiation (Rn) or sensible heat flux (H) for instantaneous surface-soil heat flux [4,6,7,34,35]. However, instantaneous soil heat flux is strongly sensitive to soil texture [28,29,30]. There are also seasonal variations of daily soil heat flux that should not be ignored [28] and need to be considered to generate accurate ET data [11].

Meteorological variables such as air temperature, pressure, humidity, sunshine duration, and wind speed are generally acquired from daily meteorological data [12,13], or the reanalysis of these datasets [10,11,36]. However, meteorological variables measured as point samples do not effectively reflect spatial variation [37]. To enhance these meteorological inputs, spatial interpolation methods are applied in order to convert the point sample data to the satellite data pixel level [38]. Under some circumstances, these interpolation methods do not provide an accurate spatial distribution of variables that are important for evapotranspiration estimation at the regional scale [39,40]. Vapor pressure deficit (VPD), as determined from air temperature and humidity [41], is an important meteorological variable since stomatal conductance decreases as the vapor pressure deficit (VPD) increases for many plant species [42,43]. The VPD from a reanalysis dataset has been incorporated in the generation of MODIS ET products [10]. However, the VPD derived from spatial interpolation of meteorological measurements [4,44,45,46] can result in a large relative error [47].

The height, temperature, humidity and pressure of the atmospheric boundary layer (planetary boundary layer, PBL) are other important and sensitive meteorological variables in the calculation of ET [48]. Atmospheric boundary layer height is variable but most existing satellite-based ET models use a fixed value at 850 hpa [12,20], and therefore they cannot capture the actual boundary layer variability [4,10]. Air friction velocity, humidity, and temperature measured by weather balloons have been used as alternative solutions [4]. However, the limited number of ground balloon sites and differences in observation times between balloon release and satellite passage seriously impair the usefulness of this approach for the acquisition of regional boundary layer measures.

Daily net radiation combines land and meteorological variables, affecting energy balance closure and, particularly, the accuracy of evapotranspiration models [49,50]. Albedo is a land variable for net radiation [4,6,10,20]. The Ångström–Prescott model has been adopted to estimate global solar irradiation [4,51,52,53,54] with sunshine duration interpolated from ground measurements and the empirical Ångström coefficients incorporated as constants for annual and daily calculations [34,52,55]. However, the sunshine duration measurements could markedly limit the spatial distribution of global solar irradiation, consequently, influencing the estimation of net radiation [4,56,57,58,59]. In addition, there are overestimations of downward shortwave radiation on rough surfaces from some meteorological reanalysis datasets [42].

Overall, poor spatial depiction of some land and meteorological variables fails to adequately describe energy and moisture transport between the land surface and the adjacent atmosphere [60,61,62,63]. This could be one of the main sources of error in estimating ET regionwide, particularly over rough surfaces which extract more sensible heat from passing air than relatively smooth surfaces [5]. Additional uncertainty derives from using limited satellite data to estimate meteorological variables in areas with fluctuating daily climatic conditions [64], which may not provide sufficiently precise information regarding water and energy fluxes at river basin scales [16,39,65,66].

Currently, there are many sensors or satellites providing land and atmospheric products such as BRDF from MODIS [39,67,68], daily multi-pass surface temperature [29,39], cloud products from geo-stationary meteorological satellite [69,70,71], as well as atmospheric column total precipitable water and atmospheric profile products from polar satellites [22,72,73]. Long-sequence atmospheric profile retrieval products (such as MYD07 and AIRS) provide new solutions for obtaining atmospheric boundary layer parameters without reliance on ground measurements [37,74]. Use of these data should better represent the spatial characteristics of both the land and meteorological variables controlling surface energy balance and atmospheric turbulence, and thus better quantify fluxes between soil/vegetation and the atmospheric boundary layer.

This paper addresses the above-mentioned challenges and limitations by presenting algorithms to generate continuous ET estimates using multi-source satellite data for key land and meteorological variables. Section 2 describes the study area, the data used, and its processing. Section 3 explains the proposed algorithms. Results and validation are discussed in Section 4 using data for the Heihe River Basin in northern China where a water crisis has emerged over the past two decades. Discussion of the data analysis is provided in Section 5 and conclusions in Section 6.

2. Study Area and Data

2.1. Study Area

The Heihe River (from 97° 24′ to 102° 10′ E and 37° 41′ to 42° 42′ N) of northwest China forms a closed basin originating in the Qilian mountains and flowing through the Zhangye oasis in its middle reaches, the Gobi Desert, the Ejina oasis downstream, and finally terminating in Juyanhai Lake (Figure 1). The basin is divided into upstream, middle, and downstream reaches by the Yingluo and Zhengyi gorges. The maximum altitude is approximately 5000 m in the upper reaches where the main types of landcover are alpine meadow and forest, the annual precipitation exceeds 350 mm, and the minimum air temperature approaches −40 °C in winter. The elevation ranges from 2000 to 3000 m in the middle reach where the main types of landcover are the Gobi Desert and oases and the annual precipitation is 100 to 250 mm. It also forms the center of the Hexi corridor of Gansu Province. The elevation is about 1000 m in the downstream area with mostly desert landcover with some oases irrigated from groundwater. The maximum annual precipitation does not exceed 50 mm, and air temperatures reach highs of 40 °C in summer. Over the past few decades, many scientific experiments and projects have been carried out focusing on ecological and hydrological aspects of the Heihe River Basin, including water conservation in irrigation, water rights allocation and other integrated water resources management issues [75,76,77,78]. Therefore, abundant ground-based measurement datasets are available including hydrological, meteorological, water and heat flux, and various eco-environmental surface parameters.

2.2. Field Observation Data

Ground-based measurements from six stations were collected from 2008 to 2014 by the WATER [77] and HiWATER [78] projects in the Heihe River Basin. Eddy covariance (EC) towers and automatic weather stations were standard equipment at each station. All of the data is freely available from the National Earth System Science Data Platform of the National Natural Science Foundation of China, including turbulent fluxes, air temperature, air humidity, air pressure, wind speed, net radiation, four radiation components, soil temperature profiles, soil moisture profiles, and soil heat flux plates. The six observation stations are described in

Table 1. List of fields station in the Heihe River Basin used in this study.

Station	Long	Lat	H	Ecosystem	Location
Arou	100.46	38.04	3033	Alpine meadow	Upstream
Yingke	100.41	38.86	1519	Agriculture, Maize	Midstream
Daman	100.37	38.85	1556	Agriculture, Maize	Midstream
Wetland	100.45	38.98	1460	Wetland	Midstream
Huazhaizi	100.32	38.77	1731	Desert	Midstream
Shenshawo	100.49	38.79	1594	Desert	Midstream
Sidaoqiao	101.14	42.00	873.00	Sparse tamarix	Downstream

Note: Long, longitude (°E); Lat, latitude (°N); and H, elevation (m).

and detailed information for each station can be found in [77,78].

EdiRe software was used to process the raw EC data, including spike detection, lag correction of H₂O/CO₂ relative to the vertical wind component, sonic virtual temperature correction, coordinating rotation using the planar fit method, correction for density fluctuation (WPL-correction), and frequency response correction [79,80,81]. Details of the data processing steps for the raw EC data and the automatic weather stations can be found in [79,80,81,82].

2.3. Remote Sensing Data

The daily MODIS 1B clear-sky data, MYD07 and AIRS, were downloaded from the NASA Goddard Space Flight Center (GSFC) Distributed Active Archive Center (GDAAC), from 2000 to 2014 [4]. The built-in ground control point (GCP) calibration coefficients were used for geometric correction, and the reflectance or radiance calculated. Atmospheric correction was carried out based on 6S atmospheric radiative transfer models [29,57]. The threshold values method with multiple properties was then used to detected cloud pixels [83]. Daily water-vapor profile data with different layers can be extracted directly from MYD07 and AIRS products. A gap-filling method was used to calculate missing pixel values from the atmospheric profiles data [74].

The daily AMSR-E brightness temperature data at 25 km resolution, from 2000 to 2014, coupled with surface temperature and NDVI data at 1 km resolution from MODIS [4,83,84,85], were used to estimate daily soil moisture contents.

Fengyun geostationary meteorological satellite (FY-2C/D/E) hourly cloud products, from 2005 to 2014, at 5 km resolution were downloaded from the National Satellite Meteorological Center (NSMC), Fengyun Satellite Data Center.

2.4. Meteorological and Radiation Data

Daily data for five national radiation stations (Gangcha, Minqin, Jiuquan, Dunhuang, and Ejin Banner) and twenty-one meteorological stations were collected from the China National Meteorological Bureau from 2000 to 2014 [4]. Radiation data included sunshine hours, solar radiation, net radiation, and reflected radiation. Meteorological data included maximum air temperature, minimum air temperature, average air pressure, average air humidity, average wind speed, and sunshine hours. Two different spatial interpolation methods were used to interpolate these point data to raster data at 1 km resolution for correspondence with the satellite data. Inverse distance squared was used for air temperatures and air pressure based on digital elevation model (DEM) data, while thin plate splines were used for air humidity, wind speed, and sunshine hours [4,86].

Raster rainfall data from 2000 to 2014 were also collected from the China National Meteorological Bureau. This daily precipitation analysis at a resolution of 0.25° x 0.25° was produced by the National Meteorological Information Center (NMIC) based on daily precipitation observation from over 2419 gauges across China [87,88,89]. Error analysis showed that approximately 91% of the absolute error was less than 1 mm/d [87].

3. Algorithms

In order to maximize information obtained from the multi-source satellite data, an integrated set of algorithms was used characterize the land surface, soil, vegetation and atmospheric boundary layer, to quantify fluxes between soil/vegetation and the boundary layer, and to determine the individual components of the surface energy balance as well as some key land and meteorological variables. Aerodynamic roughness length and surface soil heat flux are land variables, vapor pressure deficit (VPD) and the height of the atmospheric boundary layer are meteorological variables, and daily net radiation combines land and meteorological variables. Net radiation, surface-soil heat flux, and sensible heat flux are individual components of the surface energy balance [4,6,9,10,20,39]. Algorithms for the individual components of the surface energy balance and key land and meteorological variables are integrated into the surface energy balance as:

$$R_n=\lambda ET+H+G_0$$

(1)

λLE is computed as the residue of the energy balance. ET can be converted from λLE based on the temperature dependent latent heat of vaporization and the density of water. Please refer to Appendix A for variable definition for all formulas in this paper.

3.1. Daily Net Radiation (R_n,daily)

Daily net radiation was estimated from sunshine duration derived from the FY-2C/D/E cloud products from 2005 to 2014 using monthly Ångström coefficients (a_s and b_s) calibrated from ground-based meteorological and radiation measurements from 2000 to 2014 derived through spatial interpolation [56,57,58].

$$R_{n,daily}=(1-\alpha )\cdot \left(a_s+b_s\frac{n}{N}\right)R_{\mathrm{a}}-R_{nl}$$

(2)

$$\begin{array}{l}{R}_{nl,daily}=\sigma \left(\frac{T_{\max }^4+T_{\min }^4}{2}\right)\left(0.33+0.01LAI-0.15\sqrt{e_a}\right)\left(0.84\frac{R_s}{R_{so}}+0.15\right)forLAI<3\\ R_{nl,daily}=\sigma \left(\frac{T_{\max }^4+T_{\min }^4}{2}\right)\left(0.36-0.15\sqrt{e_a}\right)\left(0.84\frac{R_s}{R_{so}}+0.15\right)forLAI\ge 3\end{array}$$

(3)

$$n={\displaystyle \sum _{i=t_r+0.25}^{i=t_s-0.25}{F}_i}*t_g$$

(4)

R_nl can be calibrated with ground-based measurements [57], n is sunshine hours and can be obtained from the FY-2C/D/E hourly cloud product from sunrise to sunset at 1 km resolution based on the index (F) for cloud-type factors for FY-2C/D/E cloud type using the following equation (4) [58].

3.2. Surface Soil Heat Flux (G₀)

The instantaneous surface soil heat flux is estimated using the following equation, with the shortwave infrared band used to better reflect the influence of soil texture [29].

$$G_0{}_{,i}=R_{n,i}\cdot [T_s\cdot (a_1-a_{2}RVI+a_{3}RV{I}^2)\cdot (b{}_1+b_2{\alpha }_3+b_3{\alpha }_4)\cdot e^{-0.25/\cos (soz)}$$

(5)

The daily soil heat flux is obtained from the daily integral of the instantaneous surface soil heat flux, with diurnal variation of the surface temperature simulated with land surface temperatures from multiple temporal satellite passages. A half-order derivative solution of the heat flow equation is then used to estimate the diurnal variation of surface soil heat flux, with the general value of β set at 0.5 [29,30].

$$G\left(t\right)=\frac{\mathsf{\Gamma }}{\sqrt{\pi }}{\displaystyle {\int }_0^t\frac{dT\left(s\right)}{\sqrt{t-s}}{e}^{\left[-\beta \left(LAI/\cos \tau \right)\right]}}$$

(6)

3.3. Sensible Heat Flux (H)

The sensible heat flux is estimated with Equations (7) and (8) using aerodynamic roughness length, vapor pressure deficit [37], and atmospheric boundary layer [74] without involving the parameter KB⁻¹ [90].

$$H=\frac{\rho C_p\left(T_s-T_a\right)\left(e_s^{*}-e_s\right)f\left(LAI\right)}{r_{c\min }\left(e_s-e_a\right)F_1^{-1}{F}_2^{-1}{F}_3^{-1}{F}_4^{-4}}$$

(7)

$$f\left(LAI\right)=1.2+0.3*LAI$$

(8)

Four factors are currently considered for the aerodynamic roughness length z_0m. The calibration of the wind gradient reflects the contribution of various land surface components.

$$z_{0m}=(z_{0m}^v+z_{om}^{nir})\cdot z_{om}^T+z_{0m}^r$$

(9)

The roughness component of the vegetation condition z^v_0m was computed as:

$$z_{0m}^v={\left(b_1+b_2((NDVI)>0/NDV{I}_{\max })\right)}^{b_3}$$

(10)

The factors b₁, b₂, and b₃ are empirically determined coefficients with values of 0.001, 0.5, and 2.5 and NDVI_max is the maximum NDVI during the growing season for each year.

The roughness of the terrain is calculated as:

$$z_{0m}^T=\frac{slope-a}{b}$$

(11)

The slope was calculated based on the 90 m digital elevation model data from Shuttle Radar Topography Mission (SRTM) aggregated to 1 km spatial resolution. Factors a and b are 4 and 15, respectively, determined empirically.

The aerodynamic roughness component deriving from non-vegetated roughness features such as cities and micro-terrain was incorporated to reflect the influence of the roughness length on momentum movement:

$$\log \left(z_{0m}^r\right)=-1.221+0.0906\left({\sigma }_0\right)$$

(12)

where σ₀ represents the radar backscattering coefficient (dB) data acquired from Sentinel-1 data aggregated to 1 km spatial resolution using the nearest neighbor method.

The bidirectional reflectance distribution function (BRDF) is linked to changes in the vegetation canopy structure over time [67] and is a linear fitting of the Ross-thick and LiSparseR kernels [68]. The geometric roughness caused by the change of the vegetation canopy structure was incorporated in the aerodynamic roughness length calculation by:

$$z_{0m}^{nir}=e^{c_1\frac{k_2}{k_0}+c_2}$$

(13)

where c₁ and c₂ are empirically determined coefficients which equal 10 and 0 here; k₀ is the weighting factor of the isotropic scattering component and can be calculated as the bidirectional reflectance with both the zenith view and solar zenith angles set to 0; and k₂ is the weighting factor from the Li-sparse-reciprocal geometric scattering kernel, obtained based on surface scattering and the shadow casting theory [68].

The variable F₃ in Equation (7) is the environmental constraint factor function using vapor pressure deficit (VPD). The VPD is estimated based on the atmospheric profile product data, NDVI and surface temperature [37]. The detailed equations are:

$$VPD=e^{*}(T_a)-e_a$$

(14)

$$e^{*}(T_a)=0.6108\cdot \exp \left[17.27\cdot T_a/\left(237.3+T_a\right)\right]$$

(15)

$$T_a=T_s-k_{12}\cdot W-k_{11}\cdot NDVI-k_{10}$$

(16)

$$e_a=0.126\times {10}^6\cdot a\cdot (1+T_a/273.15)$$

(17)

where k₁₂ is a constant determined by statistical analysis, and k₁₀, k₁₁ are regression coefficients.

The variable F₄ in Equation (7) is the atmospheric boundary layer constraint factor function. A sharp decrease in the vertical profile of the water vapor mixing ratio based on atmospheric profiles and the minimum vertical gradient of the water vapor mixing ratio is used to estimate the boundary layer mixing height and retrieving temperature, humidity, and pressure at the height [74].

3.4. Land Surface Conditions

The following algorithms were adopted to address specific surface conditions based on recommendations from the World Meteorological Organization (WMO) and other sources:

• for the sublimation of snow and ice estimation [91]

  $$E{T}_{ice}=(0.18+0.098U)(e_s-e_a)$$

  (18)
• for water surface evaporation [92]

  $$E{T}_{\begin{array}{l}water\\ surface\end{array}}=\frac{\mathsf{\Delta }(R_{\begin{array}{l}n,water\\ surface\end{array}}-G_{\begin{array}{l}0,water\\ surface\end{array}})+6.43\cdot \gamma (1+0.536u)(e_s-e_a)}{\mathsf{\Delta }(\mathsf{\Delta }+\gamma )}$$

  (19)

  $$G_{\begin{array}{l}0,water\\ surface\end{array}}=0.225R_{\begin{array}{l}n,water\\ surface\end{array}}$$

  (20)
• for bare land [49], which was identified using an empirically determined NDVI threshold

  $$L{E}_{\mathrm{b}are}=\frac{\mathsf{\Delta }}{(\mathsf{\Delta }+\gamma )}(R_n-G)-\epsilon \sigma (T_s^4-T_a^4)+G\ln (T_s/T_a)$$

  (21)

3.5. Daily Surface Resistance (r_s,daily)

The Penman–Monteith model is commonly used as the gap-filling algorithm for inversing surface resistances and the ET calculation [4,10,93]. Environmental and biophysical parameters have proven adequate for the calculation of daily ET in this algorithm, including minimum air temperature, m(T_min), vapor pressure deficit, m(VPD), smoothed leaf area index, (LAI), soil moisture content, and net radiation, [9,10,59,94,95,96,97]. Of these, net radiation, soil moisture, and wind speed are more relevant with higher R² values [59]. The daily surface resistance (r_s,daily) is extended from neighboring clear days r_s:

$$r_{s,daily}=\frac{r_{\min ,clear}\times LA{I}_{clear}\times R{n}_{clear}\times S{M}_{clear}\times U_{clear}}{LA{I}_{daily}\times S{M}_{daily}\times R{n}_{daily}\times U_{daily}}$$

(22)

SM_daily can be obtained from AMSR-E data.

4. Results and Validation

The proposed algorithms were integrated into the ETWatch [4], and then daily evapotranspiration dataset from 2000 to 2014 at 1 km resolution across the Heihe River Basin was produced (Figure 2 and

Table 2. Annual rainfall and estimated evapotranspiration (ET) by landcover type from 2000 to 2014 for the Heihe River Basin (in mm).

Year	Rainfall	ET by Landcover Type						Average ET	Bias
		Forest	Grassland	Wetland	Cultivated Land	Artificial Surface	Desert
2000	106.3	402.9	313.7	371.4	436.6	347.3	89.3	132.9	−20%
2001	101.7	349.5	263.6	323.7	348.6	291.1	78.3	114.1	−11%
2002	132.2	360.3	282.2	348.1	445.3	349.4	97.0	136.0	−3%
2003	143.1	401.2	342.8	411.0	465.9	361.1	88.3	135.8	5%
2004	110.3	356.2	292.6	384.1	429.9	333.0	92.8	132.6	−17%
2005	136.4	306.7	249.4	346.3	393.9	304.5	94.8	128.4	6%
2006	118.4	322.6	265.4	370.2	415.8	312.7	84.9	122.9	−4%
2007	163.7	355.0	302.3	343.3	419.3	323.6	85.5	127.0	29%
2008	145.1	313.2	291.5	395.0	474.4	357.1	95.4	136.1	7%
2009	127.7	315.1	278.0	378.8	401.4	313.0	91.9	129.0	−1%
2010	166.0	351.3	296.8	426.1	389.7	293.6	103.1	140.4	18%
2011	130.2	397.3	326.1	398.2	384.4	298.0	94.7	136.4	−5%
2012	145.5	385.5	305.4	399.6	435.0	328.7	99.1	140.5	4%
2013	140.6	356.7	296.2	399.0	465.0	346.3	102.9	143.4	−2%
2014	141.6	353.1	298.7	404.4	425.3	319.0	97.9	137.9	3%
Average	133.9	355.1	293.6	379.9	422.0	325.2	93.1	132.9	1%

).

4.1. Spatial and Temporal Variations of ET Across the Heihe River Basin

Heihe River Basin is a typical arid inland river basin with an annual precipitation of 101 to 166 mm, averaging 133.9 mm, over the period 2000 to 2014. The average annual ET is 132.9 mm/year over the same period, showing a 99% closure in the overall balance of precipitation and consumption in this closed basin. Annually, this balance varies from −20% to +29%, with a surplus of water in the wet years (2007 and 2010), and insufficient in dry years (2000 and 2004).

The annual ET for irrigated land is more than 422 mm/year, with 93.1 mm/year for the desert. The annual ET is approximately 300 to 600 mm/year over the Ejin Oasis area. The maximum ET is observed at East Juyanhai Lake with an annual ET of more than 1000 mm/year. There is also an inter-annual pattern variability for annual ET (Figure 2). The East Juyanhai Lake is visible every year, but the West Juyanhai Lake was only visible in 2003 and 2004, and again in 2008 and 2010, when rainfall was higher than in other years. Rainfall was highest in 2010 (Table 2), and therefore the West Juyanhai Lake was larger than in other years. There was almost no water available in downstream portions of the river in 2001, since rainfall was only 101.7 mm which was the lowest of the 15 years.

The spatial ET pattern reflects both vegetation cover and climate. In the eastern part of the upstream area, the annual ET at 1 km resolution ranges from 150 to 600 mm/year; the lowest ET is located at the western part of the upstream area which has arid conditions at less than 100 mm/year (Figure 3). The areas in the eastern region with higher annual ET values have more vegetation cover and more humid conditions. Those central and eastern parts of the upstream area which have lower values of ET are mainly covered with seasonal snow, whereas the western part of the upper reaches (the black areas) are mainly covered with glaciers and permanent snow. The seasonal snow cover completely melts in summer and the ground becomes bare soil with a small ET, and therefore annual sublimation from glaciers and permanent snow in the western part of the upstream area is much higher than that of the seasonal snow in the eastern region. The ET values with large spatial variation found in the western part of the upper reaches are partly due to heterogeneity of the surface landcover.

4.2. Validation with In Situ Ground Measurements

Table 3. Validation results of daily evapotranspiration estimated by the proposed algorithms at six eddy covariance (EC) sites.

Site	Time	R	MAE (mm/day)	RMSE (mm/day)	NSE
Arou	2008.06–2014.12	0.890	0.489	0.670	0.786
Yingke	2008.01–2010.12	0.915	0.540	0.736	0.832
Daman	2012.06–2014.12	0.919	0.582	0.768	0.827
Sidaoqiao	2013.07–2013.12	0.843	0.468	0.468	0.700
Huazhaizi	2012.08–2014.12	0.880	0.278	0.278	0.761
Shenshawo	2012.08–2014.12	0.842	0.229	0.375	0.705

shows the validation results for daily evapotranspiration over the 2008 to 2014 period estimated by the proposed algorithms at 1 km resolution and the ground measurements at six EC stations, each with different land surface types, including alpine meadow in the upstream areas; irrigated oases land, and Gobi and desert in the middle reaches; and sparse vegetation (tamarix) in the downstream areas (see Table 1). The correlation coefficients (R) between observed and estimated daily evapotranspiration are 0.84 or greater, indicating a high correlation between observed and estimated daily evapotranspiration. The Nash–Sutcliffe efficiency coefficient (NSE) values are also high (above 0.7). The mean absolute errors (MAE) vary from 0.229 to 0.582 mm/day and the root mean square errors (RMSE) from 0.278 to 0.768 mm/day. These small error values between observed and estimated daily evapotranspiration indicate that the proposed algorithms can accurately estimate daily evapotranspiration.

The Yingke site in the Zhangye oasis in the middle reaches of Heihe River Basin was selected to validate the daily ET estimates. Located in an irrigated area, six irrigation events occurred during the experiment. Figure 4 (left) and Figure 5 show scatter plots of estimated and measured daily ET in 2008. The coefficient of determination (R²) between measured and estimated ET is 0.86, the RMSE is 0.69 mm, and the MAE is 0.49. Figure 5 shows that evapotranspiration decreased gradually before irrigation events as the soil dried up and, then, increased following irrigation, confirming that the proposed algorithms can capture the impact of irrigation on daily evapotranspiration ranging from less than 1 mm/day to as high as >5 mm/day.

Validation results for monthly ET against EC observations at the six EC sites are presented in Figure 6 and

Table 4. Performance for monthly evapotranspiration estimation by the proposed algorithm scheme at each of the six EC stations.

Site	R2	MRE (%)	RMSE (mm/Month)	MAPE (%)
Arou	0.96	−8.85	1.09	13.00
Yingke	0.99	−9.88	1.83	10.29
Daman	0.95	−9.88	2.13	13.52
Sidaoqiao	0.81	−3.96	5.47	28.31
Huazhaizi	0.78	−16.97	1.02	22.86
Shenshawo	0.85	15.97	1.19	31.04

. The estimated monthly ET is slightly lower than the observed EC at four sites, but more different at the two desert sites, i.e., considerably lower at Huazhaizi and higher at Shenshawo. For Arou, in the highest reach of the basin, the coefficient of determination (R²) is 0.96, the RMSE is 1.09 mm, the mean relative error (MRE) is −8.85%, and the mean absolute percentage error (MAPE) is 13%, indicating a strong relationship between estimated and measured monthly ET. For the Daman and Yingke sites in the midstream oasis area of the Heihe River Basin, and Sidaoqiao in the downstream oasis area, the coefficients of determination between measured and estimated monthly results are 0.95, 0.99, and 0.81, respectively. Desert is the major landcover type in the Shenshawo and Huazhaizi stations and the coefficients of determination are 0.85 and 0.78, respectively. Overall, the coefficients of determination are all relatively high suggesting that the proposed algorithms are able to accurately estimate evapotranspiration. Although there is some variation in this predictive ability, it does not appear to be consistently related to landcover type suggesting that the algorithms themselves are able to adequately account for landcover and related differences.

Figure 7 compares the overall monthly ET estimates from the proposed algorithms and the observed values for all six stations combined, which have different climates, soil environments, and rainfall conditions. Again, the predictive ability is good between the proposed algorithms and the EC measurements. The root mean square error is 0.80 mm/month, the average relative error is −7.11%, the mean absolute percentage error is 16.41%, and the coefficient of determination (R²) is 0.9425.

4.3. Comparison with Water Balance Estimates

The Heihe River Basin is a closed river basin without outflows, which implies mean annual evapotranspiration can be assumed to equal mean annual precipitation. As noted previously (Section 4.1), there is only a 1% difference between average annual evapotranspiration and rainfall, indicating there is good balance over the 15-year study period and implying there is about 1.0 mm of rainfall available per year on average to recharge groundwater (Table 2). The maximum surplus was 29% in 2007 when rainfall was relatively high (163.7 mm), and the greatest deficit was 20% in 2007 when rainfall was only 106.3 mm.

Palugou, a small closed watershed located in the upstream portion of the Heihe River Basin, was used to further investigate water balance. If there is no change in catchment water storage, the difference between rainfall and runoff can be taken as the actual evapotranspiration over the entire catchment. Therefore, water balance data between 2012 and 2013 were obtained for the Palugou watershed (

Table 5. Water balance measurements over the Palugou watershed (2012 to 2013).

Name	Area (km2)	Rainfall (mm/year)	Runoff (mm/year)	Water Balance ET (mm/year)	Proposed Algorithms (mm/year)	Deviation (%)
Palugou	2.89	458.59	62.83	395.76	437	9.44

). The overall deviation between estimated ET and the water balance estimates of ET from 2012 to 2013 was 9.44%, which is reasonably small, indicating that the proposed algorithms produced an accurate estimate of catchment-scale evapotranspiration.

5. Discussion

This paper presents a set of algorithms for estimating regional actual evapotranspiration using key land and meteorological variables derived from multiple satellite data sources. These data characterize the land surface, soil, vegetation and atmospheric boundary layer, and quantify fluxes between soil/vegetation and the boundary layer. This algorithm scheme takes advantage of multi-sourced satellite datasets to obtain more accurate information for key land and meteorological variables, particularly those that best represent on a region-wide basis the land surface and climate characteristics that are most relevant for evapotranspiration.

5.1. Validation of ET Estimates

This algorithm scheme was tested in a typical inland river basin, the Heihe River Basin of northwest China, where a high density of flux towers was available to provide validation data. The flux stations used in this study represent different climatic and landcover conditions, and good agreement was obtained between modeled and measured evapotranspiration. The results (Figure 6) suggest that, collectively, the algorithms accurately represent the surface energy balance components and atmospheric characteristics that are directly relevant to evapotranspiration. Good agreement between 1 km resolution ET estimates and the flux data was obtained, with RMSE and MRE values of 0.80 mm and −7.11%, respectively.

Previous studies have validated the individual components of energy balance, as well as some key land and meteorological variables, using independent ground observations. The validation of net radiation had a R² of 0.85 and the averaged Nash–Sutcliffe equation was 0.84 in research by Wu et al. [57]. The R² of surface soil heat flux was 0.96 or better and the index of agreement was more than 0.98 under different underlying surface types in research by Zhu et al. [29]. Sensible heat flux, validated using data from a large aperture scintillometer (LAS), had RMSE and MAE values less than 12.08 w/m² and 9.90 w/m², respectively, again under different underlying surface types, in research by Zhuang et al. [90].

For this algorithm scheme, only a few parameters need to be recalibrated if new data are used or different land surface types or climate systems are involved. The net radiation algorithm needs to be recalibrated when using cloud products from other geostationary satellites such as GOES or MSG for the western hemisphere. The boundary layer and VPD also need to be recalibrated for study areas with different climatic systems. This is also the case for soil heat flux if soil textures differ.

The estimated ET slightly underestimated measured ET (Figure 4, Figure 5 and Figure 7). Some of this can result from the differing formats of the estimated and observed ET values. Estimated ET is a pixel value averaged for a 1 km x 1 km area [4] which comprises a mix of many different landcovers such as bare soil, vegetation, water, roads, and buildings. Observed EC comprises a point value at a meteorological station which likely has a range of only ten to a few hundred meters. This mismatch in scales between the remote sensing pixel and flux station measurements is especially problematic for areas with heterogeneous land use [4,57,58].

The algorithms scheme was also validated in several other arid and semi-arid regions in China using the parameters derived from the Heihe River Basin. An ET dataset for the Ebinur Lake basin, derived from the proposed algorithms scheme, was used to assess ecosystem restoration under various scenarios [98]. Monthly ET, validated from July 2014 to August 2015 with the EC data, had an adjusted R² of 0.95 and a RMSE of 3.59 mm/month, indicating a strong correlation between the measured and estimated monthly ET. An ET dataset for the Turpan Depression, also derived from the proposed algorithms scheme, was used to assess a water-saving project [99]. Validation based on EC flux station observation data from 2012 to 2016 had an R² value of 0.85 and an RMSE of 20.8 mm/month. An ET dataset from the North China Plain was used to crosscheck the estimated ET against a triangle model [100]. It showed that the triangle ET and ET estimated from the proposed algorithms can provide accurate daily and monthly estimates of ET at a pixel by pixel scale. The same ET dataset was also used to validate simulated evapotranspiration from a coupled surface water-groundwater model using the MIKE SHE code for the alluvial North China Plain [101]. An R² value of 0.93 for the entire plain again demonstrated good performance from the proposed algorithms.

Further validation is provided by comparing the ET estimates derived here with the MOD16 and SSEBop satellite-derived ET products. Figure 8 relates overall ET from the algorithm scheme proposed in this paper, the MOD16 ET product [10,11], and the SSEBop ET product [102,103], all at 1 km resolution, to the measured ET at five stations variously located in the upstream, middle, and downstream reaches of the Heihe River Basin. The observations are monthly from 2008 to 2014. The proposed algorithms have the smallest RMSE and MRE values and the highest coefficient of determination as compared with those for the MOD16 and SSEBop data.

5.2. Water Balance and Related Evaluations

The performance of the algorithm scheme at a regional scale can also be assessed by a comparison with mean annual evapotranspiration from the water balance. Since the water balance does not depend upon the energy balance and is independent of remotely sensed ET, this comparison is more persuasive on a regional scale. The results in Table 5 show a deviation of 9.44% between the water balance estimate and the prediction from the algorithm scheme. This is encouraging but the data covers a relatively short time span. Additional time is required to verify further the remotely sensed ET data.

Previously published research also helps to verify the algorithm scheme proposed here. Gao and Yang [104] and Qin and Yang [105] developed a geomorphology-based ecohydrological model (GBEHM) to simulate the hydrological processes associated with vegetation dynamics in the upper Heihe River Basin. They showed that annual average ET between the model simulation results and ET from the proposed algorithms are very similar with 311 mm for the model simulation and 307 mm for the proposed algorithms. Tian and Zheng [106] proposed an integrated surface water-groundwater model (GSFLOW) for the middle and lower reaches of the Heihe River Basin. Compared to the monthly basin-wide ET estimated by the GSFLOW model, the ET estimated from the algorithms proposed here matched fairly well, However the GSFLOW model provided higher ET estimates for farmlands than the proposed algorithm scheme.

The 1% water balance deviation reported in Table 2 for the Heihe River Basin overall from 2000 to 2014 is encouraging in that there appears to be no significant ground water overdraft. The most important water issue in the basin is the unevenly temporal and spatial water distribution and consumption between basin reaches [107,108,109]. The upper reach is the runoff producing area, the middle and lower reaches are the major water consumption areas. The middle reaches of the basin have experienced an increase in population and a large expansion of irrigated agriculture leading to a higher consumption of water than allowed by the water division agreement between the middle and lower reaches [108]. This has resulted in water availability problems in the lower reaches of the basin [81], and less inflow into its terminal lake. In fact, in some dry years Lake Juyanhai has disappeared (see Figure 2) [107,108,109].

5.3. Use of Multi-Source Satellite Data in the Algorithm Scheme

Multi-source satellite data are used here to estimate the key land and meteorological variables and energy fluxes. This reduces the dependence on spatial interpolation of meteorological variables that can be problematic, especially in data-sparse regions. As a consequence, evapotranspiration estimation accuracy can be greatly improved. Figure 4 indicates that the daily ET at Yinke is changed from overestimated to slightly underestimated by using satellite rather than interpolated ground data, and the coefficient of determination is increased by nearly eight percentage points. Similarly, in Figure 5 estimated daily ET derived solely from satellite data is closer to observed ET (measured by EC) than using atmospheric variables derived from interpolation of meteorological measurements between 1 January 2018 and 31 December 2018.

Although each satellite-derived variable was separately incorporated when developing the algorithm scheme, no systematic inclusion and exclusion process was undertaken to test the independent contribution of each to improve estimation accuracy. However, there are good reasons to believe that several variables are particularly helpful. For example, to estimate sunshine radiation for daily net radiation, the hourly cloud product from a low-resolution geostationary satellite was used to capture cloud changes over time, which more appropriately reflects the spatiotemporal change of sunshine duration. This method better depicts the spatial and temporal distribution of surface shortwave radiation, as well as the net radiation [57,58].

Atmospheric boundary variables (such as height, temperature, humidity, and pressure at boundary layer height) and near surface VPD are also important and sensitive variables in ET calculation. Up to a certain point, the higher the boundary layer, the larger the difference between surface temperature and boundary layer temperature, resulting in a smaller ET due to the inverse relationship between the ET and this difference. Therefore, an accurate estimation of atmospheric boundary layer variables is important to improve ET estimation performance. However, most existing models assume that the height corresponding to 850 hPa of atmospheric pressure is the height of the atmospheric boundary layer [4]. This value is usually obtained from weather balloon data from sounding stations. Generally, there is a limited number of sounding stations in a region and, consequently, it is hard to reflect the spatial heterogeneity of the atmospheric boundary layer. Furthermore, the weather balloons only observe the boundary layer at 24:00 and 12:00. A sine function should be used to adjust the air temperature in the boundary layer from 12:00 to 13:30 (the MODIS/Aqua overpass time). However, the variation of atmospheric layer temperature does not completely follow a sine function, and therefore this conversion method can generate additional errors for the atmospheric boundary layer variables. Our algorithm scheme uses daily MODIS moisture profile products (MOD07 and AIRS) which have more than a hundred layers of atmospheric temperature and water-vapor profile data. These are used to directly obtain the spatial distribution of the atmospheric boundary layer variables [74], which also solves the problem of time inconsistency. This provides a solution for regional atmospheric boundary layer variable estimation without reliance on ground measurements.

Aerodynamic roughness length is another variable that has a strong impact on the estimation of ET. Aerodynamic roughness length is strongly influenced by the geometric characteristics of the roughness elements (buildings, roads, vegetation canopy, hills, etc.) of urban and rural surfaces. This can lead to large uncertainties in ET estimates. Currently, there are major limitations to a purely empirical aerodynamic roughness length calculation based on either landcover data or just roughness height. For vegetated areas, the height of vegetation and differences in vegetation structure are not well reflected by NDVI. The BRDF data from the MODIS satellite provides a potentially superior solution by making full use of optical, infrared, topographic, and micro-topographic data for calculation of the aerodynamic roughness length.

In addition, to better reflect the influence of soil texture on surface soil heat flux variation, shortwave infrared reflectance data were obtained for the more than 200 days each year with an average cloud coverage less than 30% and were used to estimate the surface soil heat flux.

NDVI data were used to estimate aerodynamic roughness length, vapor pressure deficit, and daily soil moisture content. Since NDVI is sensitive to soil background reflectance at low fractional cover and saturates at moderate to high fractional cover, fractional vegetation cover can provide substantially more information for surface parameters and ET estimation [110]. VPD is an important parameter of sensible and latent heat estimation at the soil-vegetation-atmosphere interface. However, recent research shows that climatic water deficit (CWD) better reflects vegetation water stress and produces good results for predicting plant stress and mortality, especially in arid and semi-arid regions [111]. Both the fractional vegetation cover and CWD should be explored to further develop the ET algorithms in order to reduce the influence of parameter uncertainty on ET estimation.

6. Conclusions

In this study, we present a scheme for regional evapotranspiration estimation using land and meteorological variables primarily derived from multi-source satellite datasets, and validate these estimates using field measurements. The multi-source satellite data used in the proposed algorithm scheme better depict regional spatial patterns and temporal characteristics than more traditional ground-based observations since they provide more accurate estimates of the surface parameters which are relevant to evapotranspiration processes.

The proposed algorithms, integrated into the ETWatch, accurately estimate daily and monthly evapotranspiration from areas with different climatic and surface characteristics which is encouraging because these areas represent alpine meadow, agriculture, and Gobi Desert with very different evapotranspiration characteristics. The methodology was validated with data collected from 2008 to 2014 across a typical inland river basin, the Heihe River Basin in the northwest of China. The estimated ET from 2008 to 2014 in the Heihe River Basin at a spatial resolution of 1 km were evaluated against different eddy covariance (EC) flux observations at the local scale and compared with the annual precipitation at the regional scale. The overall deviations for the monthly ET data as compared with ground observation data were small with RMSE and MRE of 0.80 mm and −7.11%, respectively.

It can be concluded that the proposed algorithms are reliable for estimating regional evapotranspiration under complex surface conditions. Indeed, the evapotranspiration datasets developed in this study have already demonstrated their utility. They have been used by more than 50 research groups to carry out water consumption structure analysis and irrigation management in oases in the middle reaches of the Heihe River Basin, as well as water resources estimation and evaluation, eco-hydrological process research, and studies of sustainable water resource utilization. Future studies should focus on data covering more climatic systems and land surface conditions, as well as the comparison of remote sensing estimates of regional evapotranspiration or global evapotranspiration with simulation results from distributed hydrological models.