Development of a 10 year ( 2001 – 2010 ) 0 . 1 ◦ dataset of land-surface energy balance for mainland China

Development of a 10 year (2001–2010) 0.1◦ dataset of land-surface energy balance for mainland China X. Chen, Z. Su, Y. Ma, S. Liu, Q. Yu, and Z. Xu Faculty of Geo-Information Science and Earth Observation, University of Twente, Enschede, the Netherlands Key Laboratory of Tibetan Environment Changes and Land Surface Processes, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing, China State Key Laboratory of Remote Sensing Science, School of Geography, Beijing Normal University, Beijing, China Plant Functional Biology and Climate Change Cluster, University of Technology, Sydney, P.O. Box 123, Broadway, NSW 2007, Australia


Introduction
As China is one of the fastest growing and urbanizing economies in the world, changes in land cover and land use can significantly influence the environment by altering landatmosphere energy and water exchanges (Suh and Lee, 2004;Lin et al., 2009).For instance, rapid urban expansion has substantially changed land-surface heat fluxes in the Pearl River delta (PRD) (Lin et al., 2009), and has increased sensible heat fluxes in the Beijing metropolitan area (C.Zhang et al., 2009).The variability of surface energy balance and its partitioning may also have an important impact on climate variability in China (Sun and Wu, 2001).Similarly, changes in surface energy fluxes have been shown to alter the intensity of the East Asian monsoon (Zhou and Huang, 2008;Qiu, 2013;Hsu and Liu, 2003).In short, understanding variation in energy fluxes is important for the study of climate change in China (Brauman et al., 2007).Nevertheless, the spatial and temporal variability of China's land-surface energy balance, and the magnitude of each, are still unknown.
While it is of critical importance to understand the partitioning of water and energy distribution across China's terrestrial surface, accurate monitoring of its spatial and temporal variation is notoriously difficult (Ma et al., 2011).Several Published by Copernicus Publications on behalf of the European Geosciences Union.
X. Chen et al.: Development of 10 year land surface fluxes for China field experiments are being carried out to monitor turbulent fluxes over selected land cover in China by using groundbased eddy covariance devices (Wang et al., 2010;Yu et al., 2006;Y. Ma et al., 2008;Li et al., 2009).However, these measurements are only representative of small areas around the locations where the measurements are being made.For this reason, the establishment of an eddy-covariance flux network cannot provide a complete land-surface heat flux picture for the entire Chinese landmass.
A number of methods can be used to derive land-surface energy fluxes.Jung et al. (2009), for example, generated global spatial flux fields by using a network upscaling method.However, their flux network included only a limited number of flux stations in China.The Global Soil Wetness Project 2 (GSWP-2) (Dirmeyer et al., 2006) produced a global land-surface product on a 1 × 1 • grid for the period 1986 to 1995.The Global Land Data Assimilation System (GLDAS) (Rodell et al., 2004) can provide a global coverage in the form of 3-hourly, 0.25 • data.Furthermore, products from the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim) (Dee et al., 2011), the National Centers for Environmental Prediction (NCEP) (Kalnay et al., 1996), Modern-Era Retrospective Analysis for Research and Applications (MERRA) (Rienecker et al., 2011) and other reanalysis data can also provide temporally continuous -but coarse -spatial resolution data sets of land-surface fluxes.Jiménez et al. (2011) made an inter-comparison of different land-surface heat flux products.When these products were applied on continental scales, the different approaches resulted in large differences (Vinukollu et al., 2011a;Jiménez et al., 2011;Mueller et al., 2011).
The problems met by using currently available flux data in climate studies of China have been reported by Zhou and Huang (2010).Zhu et al. (2012) have also reported that summer sensible heat fluxes derived from eight data sets (including NCEP, ERA, and GLDAS) of China's Tibetan Plateau region differ from each other in their spatial distribution.In addition, all the flux data sets mentioned above are based on model simulations, which have deficiencies in studying changes in water-cycle and land-air interactions in China (Y.Chen et al., 2013;Su et al., 2013;Wang and Zeng, 2012;L. Ma et al., 2008).
A spatially and temporally explicit estimate of surface energy fluxes is of considerable interest for hydrological assessments and meteorological and climatological investigations (Norman et al., 2003).Satellite-sensed data of surface variables can be used to produce maps of heat and water fluxes on different scales (Wang and Liang, 2008;X. Li et al., 2012;Liu et al., 2010;Vinukollu et al., 2011b).Remote sensing approaches to estimating surface heat and water fluxes have been largely used on regional scales (Fan et al., 2007;Ma et al., 2011;Jia et al., 2012;X. Zhang et al., 2009;Z. Li et al., 2012;Shu et al., 2011), but there is no analysis of satellite-derived data currently underway to produce a complete, physically consistent, decadal land-surface heat flux data set (Jiménez et al., 2009) for the Chinese landmass.The use of remotely sensed data offers the potential for acquiring observations of variables such as albedo, landsurface temperature, and normalized difference vegetation index (NDVI) on a continental scale for China.Is it possible to use all available satellite-observed land-surface variables directly to calculate high-resolution land-surface fluxes for the Chinese landmass, due to the reanalysis data having a coarse spatial resolution and containing large uncertainty?Since surface fluxes cannot be detected directly by satellite-borne sensors, an alternative for estimating continental water and energy fluxes can be derived by applying the aerodynamic theory of turbulent flux transfer (Ma et al., 2011) or by establishing statistical relationships between related satellite observations and land-surface fluxes (Jiménez et al., 2009;Wang et al., 2007).Most remotely sensed latent heat flux or evapotranspiration products have null values in urban, water, snow, barren and desert areas (Mu et al., 2007;Wang et al., 2007;Jiménez et al., 2009).This is due to the lack of a uniform representation of turbulent exchange processes over different types of land cover in their method.Meanwhile, the aerodynamic turbulent transfer method can describe the flux exchange through changes in surface roughness length over different land covers.Statistical methods establish relationships between satellite-sensed observations (e.g., NDVI, LST, albedo) and land-surface fluxes through various fitting techniques (Wang et al., 2007).The simple relationships established cannot give a reasonable approximation for extreme conditions such as bare soil or other types of non-canopy land cover (e.g., lakes, deserts), because land covers behave significantly differently in land-surface energy flux partitioning.Fortunately, turbulent flux transfer parameterization can overcome the shortcomings of statistical methods and produce spatially continuous distributions of land-surface energy fluxes with prepared meteorological forcing data.For this reason, we chose a more physically based method -turbulent flux parameterization -to produce the data set.
The challenge in using turbulent flux parameterization lies in the transition from regional to continental and global scales, because meteorological data of high resolution (i.e., 1-10 km) are not easily obtained for a large region.Recently, Chinese scientists produced high-resolution meteorological forcing data that can be used in our study.Another issue is the complexity met by the method when combining different spatial and temporal sampling input variables.This is discussed in detail in Sect.3.1.The last difficulty that has surrounded the application of turbulent flux parameterization on continental scales is the acquisition of roughness length.To address this difficulty, we have developed a remote-sensingbased mixing technique to estimate canopy heights on a continental scale, and use the resulting canopy height data set to derive, for the very first time, the dynamic variation of surface roughness length for the Chinese landmass.
Complex topography (shown by Fig. 1) and climatic conditions in China make it very difficult to obtain a clear picture of the distribution of energy and water fluxes with a high spatial resolution over a relatively long period for such a large area.In our study, we estimate land-surface heat fluxes with energy balance and aerodynamic parameterization formulas in a revised model of the surface energy balance system (SEBS) (Chen et al., 2013b;Chen et al., 2013a;Su, 2002;Timmermans, 2011); previous tests show that the revised model delivers better performance and improvements in cases where the type of land cover in China is bare soil, short canopy or snow (Chen et al., 2013a, b).Sensible heat flux in SEBS was derived from the difference between surface temperature and air temperature by using Monin-Obukhov similarity theory and bulk atmospheric boundary layer similarity (Brutsaert, 1999), which parameterizes ground surface momentum and heat-transfer coefficient maps to take into account surface roughness, canopy height, vegetation cover, and meteorological stability (Su et al., 2001;Su, 2002;Chen et al., 2013b).The latent heat flux can then be estimated from an energy balance model, assuming that surface net radiation and ground flux are known (Ma et al., 2002;Allen et al., 2011;Vinukollu et al., 2011b).We used high-resolution reanalysis meteorological data, which merge model outputs, remote sensing observations, and in situ measurements.In addition, we also assessed the accuracy of the surface energy balance terms (net radiation, sensible heat, latent heat, and ground heat fluxes) and their climatic trends in the preceding decade (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010).
After defining the equations of the SEBS model (Sect.2), we describe (in Sect.3) the input data and ground-truth measurements used in the study.Furthermore, we assess the capacity of the remote-sensing-based product to reproduce the range and variability of measured fluxes, by comparing them with in situ flux tower measurements, followed by trend analysis of the spatial patterns of the fluxes (Sect.4).Concluding remarks are found in Sect. 5.

Model description and development
The surface energy balance system model known as SEBS (Su, 2002) uses aerodynamic resistance to create a spatially coherent estimate of land-surface heat fluxes.Some model inputs can be obtained from remote sensing data, while others can be obtained from meteorological forcing data (e.g., GLDAS, ERA and NCEP reanalysis data).The model's equations and the required forcing variables are described in the remainder of this section.
The surface energy balance equation can be expressed as where Rn is the net radiation flux; G 0 is the ground heat flux, which is parameterized by its relationship with Rn (Su et al., 2001); H is the sensible heat flux; and LE is the latent heat flux.LE is computed by using the evaporative fraction after deriving the other three variables in Eq. ( 1) and taking into consideration energy and water limits (Su, 2002).As these fluxes were produced with a monthly average temporal resolution, energy storage in vegetation is not considered.
Net radiation flux is where α is broadband albedo, SWD is downward surface short-wave radiation, and LWD and LWU are downward and upward surface long-wave radiation, respectively.
Here, satellite-observed albedo is used.LWU is derived from land-surface temperature (LST) using the Stefan-Boltzmann law.Land-surface emissivity is derived as described in Chen et al. (2013a).LWD and SWD values are obtained from meteorological forcing data.
Sensible heat flux (H ) is computed according to the Monin-Obukhov similarity theory (MOST): where k is the von Karman constant; u * is friction velocity; ρ is air density; C p is specific heat for moist air; θ 0 is the potential temperature at the ground surface; θ a is the potential air temperature at height z; d is the zero plane displacement height; h is the stability correction function for sensible heat transfer (Brutsaert, 1999); and L is the Obukhov length.
In our study, θ a was obtained from meteorological forcing data, and θ 0 was derived from Moderate Resolution Imaging Spectroradiometer (MODIS) LST data.For more detailed information about u * and the calculation of L, see Su (2002) and Chen et al. (2013b).
The roughness height for heat transfer (z 0 h ) in Eq. ( 3) is calculated as follows: Using the fractional canopy coverage kB −1 at each pixel can be derived according to the following modification of the equation described by Su et al. (2001): where f c is fractional canopy coverage and f s is the fraction of bare soil in one pixel; kB −1 c is the kB −1 of the canopy; kB −1 s is the kB −1 of bare soil; and kB −1 m is kB −1 for mixed bare soil and canopy.As kB −1 is the most important parameter in a MOST-based calculation of sensible heat flux, kB −1 has been updated by Chen et al. (2013b).The momentum roughness length used to calculate kB −1 s was given a value of 0.004 (Chen et al., 2013b), and the heat roughness length of bare soil was calculated according to Yang et al. (2002).
The new kB −1 gives a better performance than the previous version of kB −1 (Chen et al., 2013a, b).Detailed evaluations of the new parameterization of kB −1 can be found in Chen et al. (2013b).
The roughness height for momentum transfer z om in Eq. ( 4) is derived from canopy height (HC), leaf area index (LAI) and the canopy momentum transfer model (Massman, 1997): where C 1 = 0.32, C 2 = 0.26, and C 3 = 15.1 are model constants related to the bulk surface drag coefficient (Massman, 1997).The three constants have been tested for several canopies (Chen et al., 2013b;Cammalleri et al., 2010) and evaluated as one of the best solutions for canopy turbulentflux parameterization (Cammalleri et al., 2010).C d is the drag coefficient, which typically equals 0.2 (Goudriaan, 1977); d is displacement height, which is derived from HC and the wind speed extinction coefficient (Su, 2002;Su et al., 2001).
As Chen et al. (2013b) have pointed out, HC is vital for turbulent heat simulations, which makes accurate estimation of HC for the Chinese landmass important for this study.A remote-sensing-based canopy height method (Chen et al., 2013b) was developed further to estimate canopy height distribution for the whole of China in this study.Simard et al. (2011) produced a global forest canopyheight map using data from the Geoscience Laser Altimeter System (GLAS) aboard ICESat (Ice, Cloud, and land Elevation Satellite).However, short-canopy (e.g., savanna, crop, grass, and shrub) height information cannot be acquired by laser techniques.Since short-canopy height usually varies by season throughout the year -crops are planted in spring and harvested in autumn -we calculated short-canopy height using an enhancement of the NDVI-based equation from Chen et al. (2013b): where HC max (LCT) and HC min (LCT) are the maximum and minimum short-canopy heights for a specific land cover type (LCT); HC min (LCT) is set to 0.002 m (Chen et al., 2013b); and HC max is set to 5, 2.5, 0.5, 0.5, and 0.5 m for savannas (including woody savannas), cropland, grassland, shrubland, and barren and sparsely vegetated pixels, respectively.MCD12C1 land cover type 1 in the year of 2002 is used to classify the pixels into savannas, cropland, grassland, and shrubland, barren and sparsely vegetated.NDVI min and NDVI max are minimum and maximum NDVI values during our 10-year study period.Each short-canopy pixel was given an NDVI min and NDVI max value to calculate the shortcanopy height.
The NDVI-based short-canopy height method above was used to fill cropland, grassland, shrubland, and barren pixels.The forest canopy heights (greater than 10 m) were assumed to be constant, i.e., with no seasonal change.By merging the forest canopy heights greater than 10 m and the variable short-canopy data, we constructed dynamic monthly maps of canopy heights for the Chinese landmass for the period of 2001-2010.These maps were then used to calculate heat fluxes.Figure 3 gives an example of derived canopy height at 11 Chinese flux stations.

Data and validation
Our modeling approach makes use of a variety of satellitebased sensor data and meteorological forcing data to estimate monthly energy and water fluxes across China.The forcing data can come from satellite-based or reanalysis data sets.Due to the influence of weather, satellite-sensed visible and thermal band data (e.g., NDVI, albedo, LST) often have spatial and temporal gaps in daily data.Various temporal and spatial gap-filling algorithms have been developed to produce continuous monthly data for satellite-sensed variables (Chen et al., 2004;Moody et al., 2005).In order to avoid both spatial and temporal gaps in the final product, we selected some specific satellite-sensed data sets for this study (see Table 1).Detailed information about each input variable is described in the following subsections.
The longest period covered by the forcing data set is approximately 31 years; the shortest is about 10 years.The spatial resolution of the data set varies from 0.01 to 0.25 • , and its sample frequency from 3 h to 1 month.The meteorological forcing data developed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences (hereafter referred to as ITPCAS forcing data) (He, 2010), were constructed to study meteorological variation in China.ITPCAS forcing data cover the entire landmass of China, and have the highest temporal resolution among the input data sets used.Other variables such as LST and albedo, for example, have coarser temporal resolutions (monthly) and global coverage.When combining data of different spatial and temporal resolutions, both spatial and temporal scaling issues need to be addressed.
Estimates of land-surface energy flux can be subject to large errors, due to bias in the meteorological forcing input data.The spatial distribution of meteorological variables is closely related to topography (Li et al., 2013).When interpolating meteorological input variables to finer scales, these effects have to be accounted for (Sheffield et al., 2006), which goes beyond the scope of our study.Therefore, we chose to resample the satellite product of a high spatial resolution to a lower spatial resolution that matches the resolution of the meteorological input data.Also, the meteorological data were averaged to monthly values that have the same temporal resolution as the remotely sensed input variables.ITPCAS forcing data provide us data of the highest spatial resolution among the meteorological forcing data currently available (e.g., ERA-Interim, NCEP, GLDAS, MERRA).Taking into account all these items, our aim was to produce a monthly product of 0.1 × 0.1 • resolution land-surface heat fluxes that contains neither spatial nor temporal gaps and that can be used to study seasonal and inter-annual variability in the hydrological and energy cycles of China.

Meteorological forcing data
In studies previous to ours, reanalysis data have been applied in many different ways, for example to construct land-surface forcing data (Sheffield et al., 2006), to detect climate trends (Taniguchi and Koike, 2008), and to investigate water and energy cycles on regional and continental scales (Roads and Betts, 2000).Reanalysis data have also been applied by the remote sensing community to derive estimates of global terrestrial evapotranspiration and gross primary production (Mu et al., 2007;Yuan et al., 2010).Few studies, however, have used reanalysis data together with remotely sensed ground data to derive global land-energy fluxes (sensible heat flux, latent heat flux, net radiation, etc.).
Researchers have developed several kinds of reanalysis data.Comparisons and evaluations of these reanalysis products with in situ observations have been performed for individual sites, specific regions, and the entire globe (Wang and Zeng, 2012;Decker et al., 2011).It is well known that inaccuracies existing in reanalysis forcing data may have substantial impacts on the simulation of land-surface energy partitioning.It is difficult to choose which reanalysis data are better for use as forcing data.Additionally, the spatial resolution of all of the above reanalysis/forcing data sets is not as high as that of remote sensing data.The ITPCAS forcing data set was produced by merging a variety of data sources.This data set benefits in particular from the merging of information from 740 weather stations operated by the China Meteorological Administration that have not been used in other forcing data.The data set has already been used to run land-surface models, and has been shown to be more accurate than other forcing data sets (Chen et al., 2011;Liu and Xie, 2013).ITPCAS meteorological forcing data include variables such as instantaneous near-surface air temperature (T a ), near-surface air pressure (P ), near-surface air-specific humidity (Q), near-surface wind speed (W s ) at a temporal resolution of 3 h, and 3-hourly mean downward surface short-wave (SWD) and downward surface long-wave (LWD) radiation.The time period covered is from 1979 to 2010; the spatial resolution has a grid size of 0.1 × 0.1 • .After spatially interpolating the monthly mean LST from 0.05 × 0.05 • to 0.1 × 0.1 • resolution, we picked out LST values of pixels that included the 11 flux tower stations from which in situ measurements were gathered.The time series comparisons of LST with the ground measurements were shown by Fig. 2. It shows that the processed monthly LST can present the seasonal variations in LST over different land covers very well.The pixel values were validated against the in situ LST measurements.Detailed information about each station is given in Sect.3.2.The linear correlation (R = 1.0),RMSE (= 1.9 K) and MB (mean value of the satellite data minus in situ observation equals 0.5 K) indicate that the quality of the merged remotely sensed monthly LST data in China is high.They also show that the monthly LST captures the in situ LST variability of different elevations and land surfaces, which described in Sect.4.2.

Albedo
Land-surface albedo determines the fraction of short-wave radiation absorbed by the ground, thus influencing the surface energy budget.Studies of land-surface energy balance require temporal and spatial albedo input data without gaps.Several research projects have been devoted to pro-ducing long-term time series of surface albedo from various satellite-borne sensors (Riihel et al., 2013;Muller et al., 2012;N. F. Liu et al., 2013).However, most of the albedo products do not provide gap-filled time-series albedo maps.Taking the MODIS MCD43B albedo product as an example, 20 to 40 % of the pixels of global landmass miss valid albedo values every year (N.F. Liu et al., 2013).Twenty percent invalid values in albedo input data will result in the same number of empty values in output, an issue that limits the albedo data that can be used in our study.After checking several albedo products (including GlobAlbedo (Muller et al., 2012), CMSAF cLouds, Albedo and RAdiation Surface Albedo (CLARA-SAL albedo) (Riihel et al., 2013), and MCD43B), we decided to use GlobAlbedo, as its data do not contain spatial or temporal gaps.This albedo data set is based on a monthly sample, and has a spatial resolution of 0.05 • , which we interpolated to a 0.1 • resolution for our study.

NDVI
The normalized difference vegetation index (NDVI) is regarded as a reliable indicator of vegetation parameters.NDVI has been used widely to explore vegetation dynamics and their relationships with environmental factors (Piao et al., 2006).NDVI data from the Systeme Pour l'Observation de la Terre (SPOT) VEGETATION sensor, distributed by Vito, have a spatial resolution of 1 km × 1 km and a temporal resolution of 10 days (synthesized on days 1, 11 and 21 of each month).In order to reduce noise resulting from clouds, the maximum NDVI value in a month for each pixel is selected to represent the canopy status of that month.

Canopy fraction
Canopy fraction (f c ) is defined as the fraction of ground surface covered by the vegetation canopy (varying from 0 to 1).f c in SEBS is used to distinguish the contributions of vegetation and soil to the roughness parameterization.Here, f c was derived from NDVI data using the following equation: (9)

Validation data
The product generated by our model needed to be validated by comparing it with an independent observational data set.
The energy balance measurement system (eddy covariance, four-component radiation and ground heat flux) at flux sites is widely accepted as a method for direct measurement of energy and fluxes, and is widely applied for assessing global evapotranspiration products (Zhang et al., 2010;Jung et al., 2011;Yan et al., 2012;Fisher et al., 2008).
Half-hourly fluxes were processed using standardized quality control procedures, which are described in the litera-ture references for each station.The half-hourly H , LE, and four-component radiation was then averaged to monthly values.Monthly average values derived from less than 70 % of the flux data in each month were not used in the validations.Gap filling was not used for the flux measurement data.

Canopy height assessment
We checked the canopy height variations at the 10 flux stations produced by Eq. ( 8), and GLAS forest height (Fig. 3).The derived canopy height for AL is about 0, which is reasonable for the local land cover.YC and WS stations, located in northern China, represent typical agricultural land, where crops mature twice per year.The highest canopy height is around 2.2 m, a similar magnitude to the height of maize in summer.The step decrease in canopy height in June at these two stations is due to the fact that wheat/maize is harvested and new seeds are sown during this period.This step variation in the canopy height also causes similar step changes in sensible and latent heat flux (shown by Fig. 5).These canopy height assessments at the observation sites enable us to consider that the developed method in this work is an appropriate one for solving scarcity of canopy height information in a continental area.

Validation against flux tower measurements
The accuracy of remote-sensing-based land-surface heat fluxes is questionable without validation against ground-   based measurements (Meir and Woodward, 2010).This subsection describes the validation of the SEBS model against heat flux measurements from a diverse range of climates.
In order to analyze the source of flux calculation errors, variables related to surface radiation fluxes were all validated against flux station observations.Table 3 shows that H and LE have RMSE values slightly less than 22 W m −2 , which is lower than the RMSE values of products of other statistical methods (see Table 7 in Wang et al., 2007, andTable 5 in Jiménez et al., 2009).Indeed, Kalma et al. (2008) assessed 30 published LE validation results obtained by using ground flux measurements, and reported an average RMSE value of about 50 W m −2 and relative errors of 15-30 %.The RMSE of our LE data set is significantly lower than their averaged RMSE value.
We also compared our validation results with those of other, similar products produced by a previous version of SEBS.Vinukollu et al. (2011b), for instance, produced global land-surface fluxes with RMSE values of 40.5 W m −2 (sensible flux) and 26.1 W m −2 (latent flux) (calculated from Table 4 in Vinukollu et al., 2011b), which are larger than those in our study.The difference could be due to the model improvement and more accurate meteorological forcing data set used in our study.Table 3 lists the values of the statistical parameters for the validation of a data product produced by GLDAS (which has the highest spatial resolution compared with other available terrestrial energy-flux data sets) against the same measurements from the Chinese flux stations as used in our study.According to the mean values of the statistical variables, the quality of our flux data set is comparable to the GLDAS model and data assimilation results.These comparisons of accuracy demonstrate that our revised model is efficient for producing a high-resolution data set of land-surface energy fluxes for China.
Net radiation has relatively higher RMSE and MB values than H , LE and G 0 in the data set, because its accuracy is dependent on the accuracy of the other variable estimates (albedo, LST, SWD, LWD, LWU, etc.).Any errors in these variables can cause bias in net radiation.LWD, for example, has a linear-fitting slope value of 0.9, with most points located around the fitting line (Fig. 4) instead not 1 : 1 line.The correlation coefficient is as high as 0.98, thus demonstrating that there is still room for improvement of the LWD algorithms.LWD in ITPCAS was calculated with algorithms developed from measurements from across the Tibetan Plateau.The LWD algorithms may not, therefore, be Table 3.Comparison of the accuracy of our flux data product and GLDAS against in situ measurements from 11 Chinese flux towers.MB is the mean of the observation minus the model simulation.

Energy flux
Radiation flux accurate for other parts of China (K.Yang, personal communication, 2013).This underlines the need for more accurate LWD radiation fluxes in order to improve the accuracy of turbulent fluxes and evapotranspiration.
In addition to the statistical evaluation of model results against observations, seasonal and inter-annual changes in the model results also need to be checked.Yucheng station, which is an agricultural experimental station with winter wheat and summer maize as dominant crops, was taken as an example (Fig. 5).Crops at Yucheng station mature twice per year, which is representative of warm temperate farming cropland, typical for the North China Plain.A 2-year flux data set was used to compare against values extracted from our model-derived product.The inter-annual and sea-sonal LST and LWU data closely match the in situ observations.The SWD term also successfully captures seasonal variations.LWD is systematically lower than observations.The LE produced at Yucheng station not only captures seasonal variation, but also responds at step stages, which occur when the wheat is harvested or maize seeds have just been sown (from June to August).The increased sensible heat and decreased latent heat flux observed in July 2003 were caused by the wheat harvest; however, this signal change is not captured by the model result.The simulated sensible and latent heat produced by SEBS has a 1-month lag when compared to reality.This phenomenon is caused by adopting a maximum monthly NDVI value, resulting in faulty representation of canopy status changes in the month of June.The Semi-Arid Climate and Environment Observatory of Lanzhou University (SC station) is situated on China's Loess Plateau, at 1965.8 m a.s.l.Annual mean precipitation there is 381.1 mm, and annual evapotranspiration is 1528.5 mm (Huang et al., 2008).Being typical of stations operating under arid conditions, its flux measurements were compared with the grid point values extracted from the model product (Fig. 6).In 2008, the land surface around the station was covered by snow from 19 January to 20 February.Consequently, the GlobAlbedo value was high for February.Unexpectedly, albedo was relatively low for January, which could be caused by the coarse temporal sampling of the station pixel by the satellite sensor.The calculated monthly sensible heat and latent heat in January 2008 have biases of −11.7 (with an observed monthly mean sensible heat of 15 W m −2 ) and −7.6 W m −2 (with an observed monthly mean latent heat of 4.8 W m −2 ), respectively.The relatively large bias for SC station when covered with snow may be caused by the mixed pixel around the station.
The results of other stations have been included in the Supplement submitted with this paper.Comparison with the results of these other stations shows that model estimates of surface energy balance variables match the magnitude and seasonal variation observed at stations in several contrasting ecosystems.Comparisons between the flux-tower-measured and modeled fluxes show that latent fluxes were more accurate than sensible fluxes.Comparisons with other studies, which are presented in Table 4, show that the accuracy of our data set is one of the best among high-resolution data sets of land-surface fluxes.

Spatial distribution of land-surface energy fluxes.
Using maps of average annual land-surface radiation and energy fluxes, we analyzed the spatial patterns of radiation and energy fluxes for the Chinese landmass and compared them with other products, such as GLDAS.The highest values of downward surface solar radiation (Fig. 7a) are located in the southwest of the Tibetan Plateau, while the lowest values occur in the Sichuan Basin (SB).The highest levels of upward short-wave radiation (Fig. 7c) occur around the snowcovered peaks of the Himalaya (HM), Karakorum (KRM), Kunlun (KLM), Qilian (QLM) and Nyainqentanglha (NQM) mountain ranges.The strongest net solar radiation (SWD minus SWU) on the Chinese landmass occurs in the southern part of the Tibetan Plateau (see supplementary materials).The downward and upward long-wave radiation (Fig. 7b and  c) on the Tibetan Plateau are the lowest for the entire Chinese landmass.Southern China has the highest levels of upward and downward long-wave radiation.The highest values of net long-wave radiation (LWU minus LWD) occur in the southern and western parts of the Tibetan Plateau (see Supplement).
Figure 8 shows that northwestern China (NWC), the western Tibetan Plateau (TP), the inner Mongolian Plateau (MP) and the Loess Plateau (LP) have the highest yearly average values for surface-sensible heat flux.Croplands of the northern China Plain (NCP, including the lowlands of Shandong, Henan, and Hebei provinces) and the northeastern China Plain (NEP, including the lowlands of Liaoning, Jilin, and Heilongjiang provinces) have low average yearly values for sensible heat flux.The Pearl River delta (PRD) and Tarim (TRB) and Sichuan (SCB) basins also have low levels of sensible heat flux, as do the Yinchuan (YCB) and the inner Mongolian basins (IMB) along the Yellow River.This spatial distribution is consistent with GLDAS results (see Supplement).
Simulated annual latent heat fluxes (Fig. 8b) exhibit a southeast-to-northwest decreasing gradient, which is consistent with other studies (Y.Liu et al., 2013).The southeastern Tibetan Plateau has high levels of annual latent heat flux.The Gobi Desert, in the northwest of China (NWC), has the lowest annual latent heat flux, followed by the western Tibetan Plateau and the inner Mongolian Plateau (MP).Lake regions along the Yangtze River and the region of basins along the Yellow River have relatively high levels of latent heat flux.
The highest levels of annual average surface net radiation (Fig. 8c) can be found in southwestern China and the Lhasa Basin (LB); the lowest levels occur in the Sichuan Basin (SCB) and the Junggar Basin (JB).The highest levels of annual average ground-heat flux (Fig. 8c) are to be found in western China, due to large amounts of incoming solar radiation that occur under dry conditions.The monthly average of G 0 is negligible when compared with other fluxes.
The role of plateau heating in Asia's monsoons is being discussed vigorously (Qiu, 2013;Wu et al., 2012;Boos and Kuang, 2010).Figure 9 shows seasonal comparisons of H between boreal winter (DJF), spring (MAM), summer (JJA) and autumn (SON).The largest area of positive sensible heating occurs in spring.Lee et al. (2011) have shown that contrasting sensible heat fluxes between the Chinese landmass and the seas surrounding it during the pre-monsoon period (April-May) affect monsoon development in East Asia.Figure 9a shows that sources of sensible heating in spring occur over the Tibetan and several other plateaus in China.During summer, the highest sensible heat fluxes are to be found on the western Tibetan Plateau, the eastern Loess Plateau (LP) and in northwestern China (NWC).
LE in summer has the largest area of high latent heating, followed by that in spring, autumn and winter (Fig. 10).Latent heat in summer is highest in southeastern and southern China, as a result of abundant rainfall in these regions.Similarly on irrigated land, such as that found in Yinchuan (YB), the inner Mongolian basin (IMB) and the downstream basins of the Tianshan (TM) and Kunlun (KLM) mountains, latent heat and evapotranspiration are high, due to the ample supply of water in summer.Latent heat fluxes in autumn and winter are significantly lower than those of the other two seasons.The magnitudes and spatial patterns of LE in China of our product are generally consistent with other reports (Yao et al., 2013;Mu et al., 2007;Jung et al., 2010).Net radiation in summer has the highest values of the four seasons.Most of the Chinese landmass acts as a source of surface energy for the atmosphere (Fig. 11).

Trend analysis
The ability to capture the inter-and intra-annual variation for each land-surface energy variable is of interest to researchers of monsoon phenomena and climate change (Zhu et al., 2012).Indeed, understanding these variations is essential for studies on climate change and water-resource-related issues.We have calculated annual average values for each flux variable.The nonparametric Mann-Kendall (MK) test is one of the most widely used methods for hydro-meteorological time series analysis (Z.Liu et al., 2013;Gan, 1998).The MK method was applied to the series of annual average fluxes to  check variations during the period 2001-2010.The resulting slope indicates that downward surface short-wave radiation increased during that decade over the majority of the Tibetan Plateau (Fig. 12).
The ground solar measurements at China Meteorological Administration (CMA) stations from 2003 to 2006, as shown in Fig. 1b of Yang et al. (2012), confirm the increasing trend of downward surface short-wave radiation found in our study.The annual mean visibility measured at these sta-tions also displays an increasing trend (Fig. 2a of Yang et al., 2012), while ERA-40 reanalyzed precipitable water and station-observed specific humidity show a decreasing trend from 2000 to 2006 (Fig. 3a of Yang et al., 2012).These results indicate that the atmosphere over the plateau is becoming drier, which would explain why SWD has increased during the decade.
The upward short-wave radiation over the Himalaya (HM), Ganges (GM), Karakorum (KRM), Qilian (QLM)  and Nyainqentanglha (NQM) mountain ranges has also decreased over the last 10 years, which may be caused by the glacial retreat that has occurred in these areas (Scherler et al., 2011;Yao et al., 2004).The Lhasa Basin (LB) has the steepest rising trend in LWU, perhaps because of the relatively greater degree of anthropogenic (e.g., urbanization) activity occurring in this area.The trend analysis did not reveal any clear spatial pattern in downward long-wave radiation.Net radiation over several high mountain ranges (including the Himalaya, Ganges, Karakorum, Qilian and Nyainqentanglha mountain ranges) increased by approximately 5 W m −2 between 2001 and 2010 (Fig. 13).The strongest increase in net radiation occurred in the central part of the Tibetan Plateau.As Matthew (2010) has pointed out, soil moisture in the central Tibetan Plateau showed an increasing trend from 1987 to 2008.Wetter soil can cause the ground surface to absorb more net radiation and thus increase latent heat flux.Moreover, wetter soil can increase soil heating capac- ity (Guan et al., 2009) and so further increase ground heat flux.The increases in net radiation and soil moisture may also explain a rising trend in latent heat in the central Tibetan Plateau.Clearly, the plateau is experiencing accelerated environmental changes (Zhong et al., 2011;Salama et al., 2012).Indeed, land-surface radiation and energy trend analyses also show that the Tibetan Plateau is experiencing a relatively stronger change in land-surface radiation -verified by Tang et al. (2011) -and energy exchange than other parts of China.

Conclusions and discussion
In view of China's highly fragmented landscape, highresolution land-surface heat flux maps are necessary for hydrological studies.As China includes arid, semi-arid, humid, and semi-humid regions, quantifying its water and energy budgets is a challenge.We have developed the surface energy balance system (SEBS) further to produce a land-surface heat flux data set on a continental scale of higher resolution than data sets derived using other methods.Generally, the global surface energy flux data sets, including reanalysis data, do not have enough spatial and temporal resolution when looking at the national-level fluxes.The surface flux data sets from reanalysis data sets still contain large uncertainty, partly due to the deficiency in their land-surface process model that simulates land-surface temperature by solving soil thermal transport equations (Chen et al., 1996), and usually result in a large error in LST simulation (Chen et al., 2011;Wang et al., 2014) if the model is not properly calibrated by measure-ments (Hogue et al., 2005).Therefore, the hypothesis tested in this paper is whether it is possible to neglect the complex process in the soil by using satellite-observed land-surface temperature directly to calculate the land-surface fluxes on a continental scale.This study has demonstrated a benchmark on how to use satellites to derive a land-surface flux data set for a continental area on a personal laptop, which is absolutely not feasible for the land-surface process modeler to do in such a time-and resource-economic way.We have overcome the shortcomings of previous remotely sensed evapotranspiration products that have null values in barren and desert areas.Usually, the surface roughness length is given a fixed value in numerical models.Here, we also found a solution on how to produce a dynamic surface roughness length due to variations in the canopy height for a continental area.This work will provide suggestions on canopy height to the numerical modelers.In summary, using remote sensing data and surface meteorological information, an independent data product of monthly resolution has been developed for landsurface heat flux analysis.We have validated our remotesensing-based approach with in situ observations from 11 flux stations in China.Taking into account the limitations of the available spatial data and computing resources, we applied the model to the entire Chinese landmass using a 0.1 • resolution meteorological data set, MODIS LST, vegetation indices and other variables to generate a climatological data set of land-surface energy balance for a 10-year period.The modeling results for both pixel-point and spatial distribution demonstrate that this approach meets our aims in terms of (a) being robust across a variety of land cover and climate types, and (b) performing well for the temporal and spatial scales of interest.The spatial distribution maps generated for each variable of surface energy balance give important background information on the terrestrial hydrology and energy cycles.This product also demonstrates the impact of topography and climatic conditions on land-air energy and moisture exchanges in China.
The applicability of remote-sensing-based estimates of land-surface fluxes is hampered by the limited temporal coverage of satellite sensors (Ryu et al., 2012).Remote sensing data are snapshots of the land-surface status at a particular point in space and time (Ryu et al., 2011).It is challenging to compare remote-sensing-based monthly flux data with ground measurements that are made on timescales ranging from half-hourly through to monthly.The accuracy of land-surface heat fluxes is largely dependent on the remotely sensed land-surface temperature.Here, we have made an assumption that the averaged Aqua and Terra sensorsensed LST in each month can represent the monthly average LST.The Terra satellite sensor passes twice a day (at about 10:30 a.m. and 22:30 p.m. local time); the Aqua satellite also passes twice a day (at about 01:30 a.m. and 13:30 p.m. local time); so, MODIS has maximally four samples each day.The samples may not be enough for calculating the monthly LST, also due to the cloud noise.Besides, the time length of MODIS data sets is not longer than 15 years, which may limit the application of our data set in climate analysis.Additionally, the sensible heat flux over forest is underestimated by the present turbulent flux parameterization method in SEBS, which does not take the roughness sublayer over high canopy (Bosveld, 1999) into consideration.The low bias in the wind speed of the ITPCAS forcing data set (not shown here) could also be one reason for the lower estimation of sensible heat flux by our method.
The energy flux product we have developed has a spatial resolution of approximately 10 km, while flux towers have a footprint of tens to hundreds of meters.The tower footprint may not be representative of the larger pixel of the product, and this mismatch will result in errors if the mean of the satellite pixel is different from that of the flux tower footprint.Remote-sensing-based studies stress that direct comparison is a challenge because scale mismatch (Norman et al., 2003) and the heterogeneity of the land surface reduce the spatial representativeness of ground-site measurements (Mi et al., 2006).Another challenge is validating the grid-box-based simulation results on the scale of the Chinese landmass, since reliable observations of flux data are only available from a few sites in the simulated region.
Potential effects of changes in land-surface heat fluxes on the monsoon over East Asia (Lee et al., 2011) as a result of China's recent urbanization can be studied further using our product.As an independent satellitebased product, it can also be used as a data source for evaluating land-surface models.We also produced an evapotranspiration product for the China land area using the data set in this paper.The land-surface fluxes and evapotranspiration product can be downloaded from the URL.The recent product will be shared when the input data set is available: https://drive.google.com/folderview?id= 0B7yGrB1U9eDec2JFbnA5eldlVHc&usp=sharing.
The Supplement related to this article is available online at doi:10.5194/acp-14-13097-2014-supplement.

Figure 2 .
Figure 2. Time-series comparison of monthly averaged LST derived from MOD11C3, MYD11C3 and in situ measurements.

Figure 3 .
Figure 3. Monthly variation of canopy height at the 10 flux stations.

Figure 5 .
Figure 5. Time-series comparison of SEBS input and output variables against measurements at Yucheng station.Black lines are SEBS results; red lines are measured values.

Figure 6 .
Figure 6.Time-series comparison of SEBS input and output variables against measurements at SC station.Black lines are SEBS results; red lines are measured values.

Table 1 .
Input data sets used for calculating land-surface fluxes for China (see Sects. 2 and 3 for an explanation of abbreviations).
(Coll et al., 2009;Wan and Li, 2008), with an average bias of less than 1 K(Coll et al., 2009;Wan and Li, 2008).The MODIS V5 monthly LST product, MOD11C3 and MYD11C3, has a 0.05 • grid size, without gaps, and covers the period March 2000 to the near-present.It provides monthly daytime and night-time LST values.In our study, we averaged the daytime and night-time values of MOD11C3 and MYD11C3 to represent monthly means.

Table 2 .
Flux tower sites supplying measurement data for product validation.

Table 4 .
Comparison of statistical values reported in similar studies.