Regional Estimates of Evapotranspiration over Northern China Using a Remote-sensing-based Triangle Interpolation Method

Arid and semiarid areas represent more than a third of China's mainland (Liu et al., 2008), and much of them are located in northern China. The hydrological processes in this region are important because of the increasing water resources shortage caused by aridification over the past decades (Fu and An, 2002). Evapotranspiration (ET) comprises two subprocesses: evaporation from the Earth's land surface, and transpiration from the plant canopy. As a key component of the water cycle, ET is a critical process in land-atmosphere interaction studies (Sobrino et al., 2007), and plays an important role in understanding climate dynamics and photosynthesis of the terrestrial ecosystem (Nemani et al., 2002; Nishida et al., 2003). Continuous and accurate estimation of spatial and temporal variations of ET at the regional scale is critical for improving the ability of water and land resource management, drought detection, and climate simulations (Cleugh et al., 2007; Wang and Liang, 2008; Jiang et al., 2009).

Eddy covariance (EC) flux measurement is accepted as an important method for observing latent heat flux, heat flux and other fluxes between the land surface and atmosphere at the site scale, and is widely applied in the network of near-surface flux observation (Foken and Wichura, 1996; Baldocchi et al., 2001). However, the footprint of EC measurements normally does not go beyond several square kilometers over a heterogeneous land surface (Schmid and Lloyd, 1999; Chen et al., 2008), and changes from time to time. Therefore, current field-based point measurements cannot capture reliable signals of hydrological dynamics over large areas.

Satellite remote sensing is recognized as the most viable means to retrieve large-scale distribution ET data (Li et al., 2009; Mallick et al., 2009), and is considered as a promising tool for the study of global energy balance, land surface processes, and the hydrological cycle. Remote sensing models integrate satellite and field observations to obtain a reasonable spatial distribution of land surface parameters from the landscape to regional scale. Although there have been large numbers of applications of remote sensing of ET at regional scales, at least two major problems still exist: (1) The conflicting requirements of simple data over a wide region and accuracy of those data (Cleugh et al., 2007). Regional remote sensing methods need highly accurate ground-measured near-surface data, such as wind speed, air temperature and water vapor pressure (Norman et al., 1995; Bastiaanssen et al., 1998; Bastiaanssen, 2000; French et al., 2003; Jia et al., 2003; Nishida et al., 2003; Liu et al., 2006), but assumptions about the extrapolation of atmospheric parameters and air resistance may not work over large heterogeneous areas (Venturini et al., 2004). It remains a challenge to obtain such effective regional near-surface data, and air resistance, being a key parameter in these models, is also hard to obtain over large areas.

(2) There are few regional products of ET that cover the entire northern China region. Some studies have concentrated on the local or watershed scale (Shu et al., 2011; Luo et al., 2012), and most existing global ET products are calibrated with North American and/or European data; few are validated using data from arid and semiarid China (Mu et al., 2007, 2011; Fisher and Baldocchi, 2008; Vinukollu et al., 2011), and thus may not be suitable to studies in this area. Regional-scale ET products that are well validated are badly needed for northern China.

Recently, Jiang and Islam (1999, 2001) proposed a simplified Priestley-Taylor method (Priestley and Taylor, 1972) based on the interpolation of the evaporative fraction (EF; the ratio of latent heat to available energy) from the triangular distribution of land surface temperature (LST) and vegetation index (VI). Later, (Wang et al., 2006) combined this EF interpolation method with thermal inertia to propose a new interpolation method based on the distribution of day/night LST difference (ΔT ) and VI. This interpolation method can obtain EF directly by using satellite remote sensing observations and auxiliary meteorological data. Besides, as the retrieved EF is mainly determined by the overall feature space of ΔT-VI scatter plots, this method does not rely on the high precision of surface data (Wang et al., 2006). As ET can be calculated from EF and available total energy, EF interpolation from the triangular-shaped scatter space of ΔT-VI might also be a potential method for the estimation of ET over large, heterogeneous land surface areas, such as northern China.

In this study, we estimated ET over northern China using the modeling method proposed by (Wang et al., 2006). With auxiliary near-surface meteorological data, we first retrieved EF from interpolation based on the contextual space of the ΔT-VI distribution. Then, we simulated regional ET for the year 2008 at a spatial resolution of 0.5 km and temporal resolution of eight days, validating the results with in-situ data from EC flux measurements. Finally, we analyzed the spatial patterns of EF and ET over the study area. Dynamic changes in ET across main land cover types in the study area were also analyzed. The aim of the work was to find a simplified method for estimating ET over large and heterogeneous areas, and to some extent improve our understanding of the regional hydrological cycle by integrating remote sensing and auxiliary meteorological data.

Northern China extends latitudinally from 31.39 °N to 53.55°N and longitudinally from 73.45°E to 135.08°E, covering about 4.3 million square kilometers. Much of the area is covered by arid and semiarid regions, and a west-east humidity gradient is shaped by precipitation.

2.2.1 Flux data for validation

Figure 1.  Geographical location of the study area and the EC flux sites. The map of the study area was generated from the land cover type product (MCD12Q1) of the MODIS Land Team. Biome class key: water=water body; ENF=evergreen needleleaf forest; DNF=deciduous needleleaf forest; DBF=deciduous broadleaf forest; MF=mixed forest; Oshrub=open shrubland; Grass=grassland; Crop=cropland; UBU=urban and built-up; SOI=snow and ice; BSV=barren or sparsely vegetated.

Field-measured EF and ET were calculated from 14 EC flux sites operating under the Enhanced Coordinated Observation of Land-Atmospheric Interactions project in northern China (Wang et al., 2010). This network was established with the aim to enhance understanding of the spatial pattern and seasonality of water, energy and carbon flux in arid and semiarid areas of China since 2008. The EC flux equipment consists of a 3D sonic anemometer and an open-path fast response infrared gas analyzer. The 14 sites are widely distributed and represent very well the main land cover types and dominant biomes over the region (Fig. 1).

A uniform and intensive calibration of the infrared gas analyzers was performed weeks before the observation period (July-September 2008) to ensure the quality of data and the comparability between them. The raw flux data were recorded at a high rate of 10 Hz by a fast response data logger (Model CR5000, Campbell Scientific Inc., Logan, Utah, USA). Data processing of all the flux records from the 14 sites were performed with EdiRe software (Robert Clement, University of Edinburgh, UK), which included coordinate rotation, the WPL term (Webb et al., 1980), and sonic air temperature correction.

2.2.2 Meteorological data

Firstly, two auxiliary meteorological datasets—daily air temperature and daily shortwave solar radiation, from stations in and near the study area—were interpolated to the grid-data layers using kriging. The air temperature data were retrieved from 390 meteorological stations provided by the China Meteorological Data Service system (http://data.cma.gov.cn), while the radiation data were retrieved from 52 out of these 390 meteorological stations. The daily raster datasets were then averaged to the same time intervals of the MODIS land product as the 8-day composites.

2.2.3 Remote sensing data

The Digital Elevation Model (DEM) data were obtained from the U.S. Geological Survey (USGS) global DEM dataset (http://rmmcweb.cr.usgs.gov/el- evation/) in order to analyze the complex topography and calculate the atmospheric pressure in the region.

MODIS land products (C5) developed by its Land Team were used to derive surface reflection (MOD09A1), day/night land surface temperatures (MOD11A2), and land cover types (MCD12Q1). MOD09A1 was used to generate time series of NDVI (Normalized Difference Vegetation Index) and albedo over the study region for the whole of 2008. MOD11A2 was used to retrieve the difference of land surface temperature between day and night. MCD12Q1 was employed to distinguish between the features of EF and ET for different land cover categories as follows:

where I indicates NDVI and albedo (α) was derived from five bands of the MOD09A1 dataset (Liang, 2000). Ρ₁, Ρ₂, Ρ₃, Ρ₄, Ρ₅, Ρ₇ are the reflectance bands in MOD09A1: red band (620-670 nm); near infrared band (841-876 nm); blue band (459-478 nm), green band (545-565 nm); and two shortwave infrared bands (1230-1250 nm and 2105-2155 nm), respectively. These MODIS products have been improved in terms of spatial and spectral resolutions, atmosphere correction and calibration algorithm, as compared to the old generation sensors, such as the Advanced Very High Resolution Radiometer (AVHRR).

We considered two major factors—the footprint of EC flux measurements and the geolocation accuracy of satellite products—to decide the size of subsets, and the flux footprints were retrieved from an Eulerian analytical approach Kormann2001. Then, subsets of 3×3 pixels of MOD09A1 (approximately 1.5×1.5 km²) that centered on each eddy flux tower site were selected.

: A general form of ET can be expressed by the following formulation (Parlange and Albertson, 1995):

where L is latent heat; φ is the Budyko-Thomwaite-Mather parameter; A and B are dependent by models; Δis the slope of saturated vapor pressure on air temperature; γ is the psychometric constant; R_n is net radiation; G is soil heat flux; f(u) is a function of wind speed; e_a^* is the saturation vapor pressure; and e_a is the vapor pressure. Regional estimation of ET with Eq. (3) is strongly restricted by the difficulty in obtaining its representative parameters over a highly heterogeneous landscape. However, it can be simplified to a form of Priestley-Taylor formulation by setting A=1 and B=0:

where E_f indicates EF, and φ is a complex effects parameter (dimensionless) that accounts for soil moisture and wind speed. In this study, ET was estimated with an algorithm based on Eq. (4), and the value of φ varying from 0 to (Δ+λ)/ Δ. Thus, Eq. (4) was considered as an extended Priestley-Taylor equation, which allows the application of φ over large heterogeneous areas with different water conditions, while Priestley-Taylor evaporation is primarily applicable for only well-watered areas. The maximal value of φ was set to 1.26 here, because it is relatively insensitive to small changes in the absolute value of atmospheric parameters such as air temperature (Wang et al., 2006).

The interpolation scheme of φ based on the ΔT-NDVI contextual space for each pixel in an image was divided into two steps. The first step was to determine two bounds controlled by physical states of the land surface: the upper decreasing dry and lower horizontal wet edges of the φ value. These two edges were generated fromΔT-NDVI 2D scatter plots (Fig. 2). The second step was to interpolate within the contextual space. Pixels were ranked into several classes according to the sequence of their NDVI value. Then, φ was retrieved from the linear interpolation between zero and the highest temperature difference within each NDVI class. Further details of the interpolation scheme for image pixels can be found in Jiang and Islam (1999, 2001).

Surface net radiation is the sum of surface shortwave and longwave net radiation, and is estimated by the following equation:R_n=(1-α)R_a+ε_aσT⁴_a-ε_gT⁴_a (5)

Where α is surface albedo; R_a is shortwave radiation; ε_a is atmospheric emissivity; σ is the Stefan-Boltzmann constant (5.67×10^-8 W m² K^-4); ε_a is surface emissivity; and T_a and T_a are atmospheric temperature (K) and surface temperature (K), respectively. Soil heat flux (G) is the energy that is utilized in heating the soil, which can be estimated from net radiation and the vegetation fraction (Su, 2002):G=0.35R_n(1-P_v)+0.05R_nP_v (6)

where I_min and I_max are the minimal and maximal NDVI values, respectively, over the entire study region.

Using these key factors as inputs, ET was estimated by partitioning the available energy from the equations above. The whole process can be expressed as shown in Fig. 3.

Figure 2.  Triangular-shaped scatter plots of the difference of day/night land surface temperature (ΔT) and normalized difference vegetation index (NDVI) values.

We evaluated the accuracy of the estimated EF and ET by comparing them to EC flux tower measurements at the 14 sites mentioned above. The observation period was from July to September, and covered most of the growing season in the study area. The estimated EF and ET from remote sensing were retrieved from the mean value of 3×3 pixel subsets (1.5×1.5 km²) around flux sites. Liner regression and relative errors with corresponding field-measured EF and ET were performed. The estimated EF and ET agreed well with flux tower measurements (Fig. 4). For all sites, this agreement for both EF and ET was highly statistically significant (p<0.01), and the overall correlation coefficient of EF (R²=0.64) was relatively higher than ET (R²=0.57). Validation with EC-measured ET showed the mean RMSE and bias was 0.78 mm d^-1 (22.03 W m^-2) and 0.31 mm d^-1 (8.86 W m^-2), respectively. Meanwhile, the flux sites were categorized into three main land cover types: grassland, cropland, and forest. A comparison of means and the relative errors between EC-measured and model-simulated results was performed (Table 1). Crop sites had the lowest relative errors, followed by grass sites, while forest sites had the largest relative errors. Most of these relative errors were less than 20% and the layers of EF and ET were simulated with reasonable accuracy in general.

To further evaluate the performance of output layers, the annual value of satellite-estimated EF and ET over the 3_×3 pixel subsets surrounding each flux site were ordered along multi-year average precipitation (Table 2). Both simulated EF and ET showed an increasing trend along the gradient of precipitation, and annual ET of these sites corresponded well to the annual precipitation. In croplands, EF ranged from 0.26 to 0.36 and ET ranged from 357.91 mm yr^-1 to 462.56 mm yr^-1. For grasslands, both EF and ET were greatly heterogeneous, with EF ranging from 0.15 to 0.31 and ET from 204.59 mm yr^-1 to 583.76 mm yr^-1, respectively. As for forest sites, the two deciduous broadleaf forest sites (CW and MY) had very close EF (0.4 vs 0.42) and ET (565.24 mm yr^-1 vs 568.92 mm yr^-1), while the evergreen needleleaf forest site (DYK) was much lower for both parameters (0.34 and 322.48 mm yr^-1, respectively). Compared to the value of EF, sites with high elevation, such as AR and DYK, had less annual ET, probably due to stronger low-temperature-stress than other sites.

Figure 5.  Spatial pattern of annual mean EF in the year 2008 over northern China.

Annual mean EF and annual total ET in the year 2008 were generated and are shown in Figs. 5 and 6, respectively. The above results suggest that, although current estimates may still contain biases, general spatial and temporal patterns of EF and ET over the region were well captured by our model simulation. Spatial patterns of EF and ET were quite well matched to each other, except in areas with high elevation or high latitude. ET in those areas was lower than other areas with similar EF, probably due to the stress of low temperature. Generally, the annual mean layer of ET captured the spatial structure well. The spatial pattern of ET also corresponded to distributions of the major biome types. The ET over the study area increased along a clear longitudinal gradient from lowland areas in the east to the high plateaus of the west. This spatial pattern was probably controlled by the gradient of precipitation, green vegetation fractions, and the intensity of human activities. In 2008, the mean ET was 453.87 mm and the maximal ET was 1351.36 mm over the entire region. As expected, forest and cropland that dominated in the eastern part of the study area had the highest ET values, while bare soils, open shrublands and sparse grasslands from desert areas or short-growing areas in the western part had the lowest estimates. In the central part, dense grasslands in the east were gradually replaced by sparse grasses and shrubs farther inland, and ET dropped gradually from east to west. Meanwhile, the modeled ET also captured the heterogeneity of the spatial distribution well. The heavily degraded vegetation in the eastern part, as well as croplands, grasslands, meadow steppe and evergreen needle forest existing as oasis or sub-alpine regions in the western part, and big water bodies throughout the whole study area, were clearly identified by the modeled ET.

According to MODIS MCD12Q1 data in 2008, main land cover types over the region were classified into eight types: evergreen needleleaf forest (ENF), deciduous needleleaf forest (DNF), deciduous broadleaf forest (DBF), mixed forest (MF), open shrubland (Oshrub), grassland (Grass), cropland (Crop), and barren or sparsely vegetated (BSV).

Time series plots of these main land cover types were generated to further analyze the seasonal pattern of ET, which clearly identified the seasonal variability of ET (Fig. 7). For different land cover types, ET varied greatly at the time of onset, at the end of the growing season, and in terms of the magnitude of variation. On average, at the annual scale, BSV had the lowest ET, which related to its low vegetation fraction. Seasonal variation was not significant and the magnitude of variation was greater than 1 mm d^-1 only in the rainy season (May-September). Grass and Oshrub were mainly distributed in northern and western parts of the study area. Their ET was slightly greater than BSV, and may be limited by low temperature and precipitation. Crop showed similar values to forests, but the time for peak ET was one week later, probably delayed by farmland management. Forests, including DBF, ENF, DNF and MF had the greatest ET and also showed great magnitude of variation.

Figure 7.  Time series of eight-day ET for the main land cover types over northern China.

The average values of EF and ET for the major land cover types are summarized in Table 3. Forests accounted for about 10% of the study area and had the largest ET (MF: 828.27 mm; DBF: 831.03 mm; DNF: 698.42 mm; ENF: 939.33 mm). Crop was next, accounting for about 20% of the study area (658.31 mm). Grass (446.58 mm) and Oshrub (539.71 mm) were less than Crop and accounted for about 30%, while BSV had the lowest ET (223.89 mm) and accounted for about 30% of the study area. EF had the same distribution as ET among the land cover types.

Although the simplified interpolation from the triangular LST (ΔT)-VI scatter space may introduce errors to the model results, this model efficiently estimates ET with limited surface ancillary data (Venturini et al., 2004) and avoids complicated parameterizations over large areas (Jiang et al., 2009). According to the validation result mentioned above, this model is practical and suitable for large areas such as northern China. Moreover, EF was retrieved based on a con textual relationship between remotely sensed temperature and vegetation index data in this model, and the retrieved EF is insensitive to changes in absolute values of surface temperature, which can be affected by local atmospheric conditions and sensor viewing angles (Venturini et al., 2004). It can be inferred that this triangle interpolation method could reduce model uncertainties caused by the errors of precision from the input remote sensing data. (Jiang et al., 2004) evaluated the error properties of the model based on Gaussian Error Propagation (GEP) and found that uncertainties of the LST (ΔT)-VI model are not larger than 20%, exhibiting a similar or better error bound than other ET models.

In the LST (ΔT)-VI model, a major assumption of a negative relationship between latent heat flux and LST or changes in LST may not work at high latitudes and in cold areas, where temperature is a main limitation of ET. Therefore, this model is more suitable for the growing season in mid-latitude areas (Wang and Dickinson, 2012). In our study, the flux data for validation were collected during the growing season (July-September), but the seasonal variation of ET (Fig. 7) showed reasonable dynamic curves for the whole year.

As the outliers and spurious dry points may affect the determination of the dry edge (Tang et al., 2010), it is necessary to filter the scatters of the LST (ΔT)-VI triangle space. In our study, we found that outliers with highΔT were mainly distributed in low vegetation fraction areas such as the Gobi Desert. These scatters were filtered out when determining the warn edge and the φ of the scatters beyond the dry edge was set as zero.

Besides uncertainties in the model, several factors may contribute to the bias between remotely-sensed estimates of ET and field-measured ET. The first factor is the mismatch in spatial scales between the satellite footprints and the EC measurement footprints. The footprint of EC measurement is a function influenced by observation height (Schmid, 1997), canopy structure, and other local environmental conditions (Yi, 2008; Yi et al., 2010). Therefore, the footprint of EC can be seen as a changing variable over both temporal and spatial scales, and the size of the EC footprint changes mainly from 0.01 to 10 km², while the satellite footprint can be seen as a static constant with fixed pixel size. In heterogeneous areas, the differing scales of the tower and satellite ET can introduce different surface information from different land cover types. This can make footprints from the two observation platforms hard to match with each other, and can further influence the performance of ET modeling (Vinulollu et al., 2011). Issues related to mismatches between different observation platforms need to be further explored.

Another source of bias for ET modeling is the flux measurement itself. Currently, EC flux measurement suffers from two major problems. The first is the uncertainties in field measurements. Potential factors influencing measurement accuracy are systematic errors from turbulence, random errors from sample size, vertical and horizontal advection, and subjective decisions on gap filling (Loescher et al., 2006). EC flux measurements can have uncertainties of about 10%-30%, even at the same site (Mu et al., 2011). Another problem with EC flux measurements is the energy balance closure issue. That is, the phenomenon that total surface layer heat flux (the sum of latent and sensible heat) is commonly less than the available energy (the sum of net radiation and ground heat flux) in EC flux observations. The ratio between them is approximately 60%-90% (Franssen et al., 2010; Liu et al.,2011). Several factors may contribute to the enclosure issue, including different source areas of energy flux measurement, loss of large eddy turbulent fluxes due to high frequency sampling, and the absence of other energy terms (Kidston et al., 2010).

The third source may be cloud contamination and other atmospheric disturbances. The constant EF method is widely used to scale up instantaneous ET to daily or even larger timescales, such as eight-day values with the assistance of total solar radiation, and it is accepted that instantaneous EF data at satellite overpass times have good agreement with daily averaged EF, and thus the bias between them can be neglected (Jiang and Islam, 2001; Sun et al., 2012). However, cloud contamination can lead to major errors if no specific correction is performed to remove atmospheric disturbances. Currently, MODIS land product integrated satellite observation data operate over a timescale of 8 to 16 days in clear sky conditions. The reflective spectral information of cloud covered surfaces cannot be captured by optical or thermal infrared sensors. As EF can be affected by different weather conditions, ET retrieved from clear skies may not represent other weather conditions well. In this study, we adopted accumulative radiation data to better reflect the temporal changes in radiation. However, owing to possible effects of cloud contamination, the modeled EF from satellite data could only represent values in clear-sky days. Therefore, uncertainties in satellite-based EF modeling may introduce new errors by affecting energy partitioning. Future work should address the use of geostationary satellite data to integrate polar satellite data for better capturing diurnal changes in EF (Anderson et al., 2011; Shu et al., 2011; Sun et al., 2012).

We explored the application of an EF interpolation method at the regional scale to assess the validity of this remote-sensing-based algorithm. ET was estimated with auxiliary near-surface meteorological data and the EF retrieved from the contextual information of the 2D scatter space between LST and NDVI. We generated eight-day layers of ET over northern China for the year 2008. These modeled layers adequately captured the spatial and seasonal variation of ET. The estimated eight-day EF and ET were validated with 14 EC flux measurements in the study area, and these agreed well with flux tower measurements.

Compared to the traditionally used residual method, which estimates latent heat flux indirectly, with a small error in the input surface air temperature sometimes introducing large errors in ET output results, the newly proposed EF interpolation method is operationally feasible and has the potential to provide reasonable ET estimation over large heterogeneous areas with a limited number of meteorological input parameters.

We evaluated the triangle interpolation method and discussed the potential factors that contribute to the discrepancy of ET between remotely sensed modeled results and EC flux measured results. In short, these factors can be summarized as the mismatch between the source areas of satellite observations and EC flux field observations, the uncertainties that lie in EC flux measurements and post-processing methods, and cloud contamination and other atmospheric disturbances. Future work will address these issues to improve regional ET modeling.