Assessment of Spatial Representativeness of Eddy Covariance Flux Data from Flux Tower to Regional Grid

Combining flux tower measurements with remote sensing or land surface models is generally regarded as an efficient method to scale up flux data from site to region. However, due to the heterogeneous nature of the vegetated land surface, the changing flux source areas and the mismatching between ground source areas and remote sensing grids, direct use of in-situ flux measurements can lead to major scaling bias if their spatial representativeness is unknown. Here, we calculate and assess the spatial representativeness of 15 flux sites across northern China in two aspects: first, examine how well a tower represents fluxes from the specific targeted vegetation type, which is called vegetation-type level; and, second, examine how representative is the flux tower footprint of the broader landscape or regional extents, which is called spatial-scale level. We select fraction of target vegetation type (FTVT) and Normalized Difference Vegetation Index (NDVI) as key indicators to calculate the spatial representativeness of 15 EC sites. Then, these sites were ranked into four grades based on FTVT or cluster analysis from high to low in order: (1) homogeneous; (2) representative; (3) acceptable; and (4) disturbed measurements. The results indicate that: (1) Footprint climatology for each site was mainly distributed in an irregular shape, had similar spatial pattern as spatial distribution of prevailing wind direction; (2) At vegetation-type level, the number of homogeneous, representative, acceptable and disturbed measurements is 8, 4, 1 and 2, respectively. The average FTVT was 0.83, grass and crop sites had greater representativeness than forest sites; (3) At spatial-scale level, flux sites with zonal vegetation had greater representativeness than non-zonal vegetation sites, and the scales were further divided into three sub-scales: (a) in flux site scale, the average of absolute NDVI bias was 4.34%, the number of the above four grades is 9, 4, 1 and 1, respectively; (b) in remote sensing pixel scale, the average of absolute NDVI bias was 8.27%, the number is 7, 2, 2 and 4, respectively; (c) in land model grid scale, the average of absolute NDVI bias was 12.13%, the number is 5, 4, 3 and 3. These results demonstrate the variation of spatial representativeness of flux measurements among different application levels and scales and highlighted the importance of proper interpretation of EC flux measurements. These results also suggest that source area of EC flux should be involved in model validation and/or calibration with EC flux measurements.


Introduction
The eddy covariance (EC) technique is a micrometeorological tool for direct measurements of the exchanges of carbon, water, and energy between land surface and atmosphere [1].Temporal continuous flux data for a certain vegetation type provided by EC is essential for quantifying biogeochemical cycles [2,3] and for understanding key processes of land-atmosphere interactions [4,5].As the number of flux sites increases sharply, a larger regional/global coverage of EC flux sites and denser EC flux network, intergraded with remote sensing and/or land surface models, greatly benefited the research of quantifying magnitudes and dynamics of large scale carbon and water budgets over different eco-regions in the world [6][7][8].However, direct use of in-situ flux measurements can lead to scaling bias due to heterogeneous nature of the vegetated land surface, the changing source areas and the mismatching between source areas and grids modeled by remote sensing or land surface models [9,10].
A feasible scheme to resolve such problems is to assess spatial representativeness of flux data measured by EC.EC flux data were mainly applied in two aspects: one is using EC data to represent the flux of target vegetation type, which is called vegetation-type level; another is using EC data to represent the flux in different spatial scope, which is called spatial-scale level.Therefore, the flux data should also be assessed based on these two aspects.For the vegetation-type level, spatial representativeness was defined as how well the flux measurements represent the target vegetation types in the source area.It can be expressed as fraction of target vegetation type (FTVT) [11,12].For spatial-scale level, flux measurements were assessed to identify how large an area the flux measurements could be up-scaled to.It can be expressed by a bias estimation of vegetation indices such as Normalized Difference Vegetation Index (NDVI) from EC source areas and from remote sensing or land surface model simulated grids [13][14][15][16].NDVI generally correlates with vegetation activity, is an indicator of LAI, biomass and photosynthesis [17][18][19].
The source area of EC measurements can be calculated from footprint modeling.The footprint distribution is the relative contribution of different parts in the upwind source area to a measured scalar EC flux [20], which is a function of observation height [21], vegetation properties, surface roughness, and other local environmental conditions [22].The footprint probability distribution function (PDF) is the source strength of flux elements over a heterogeneous surface.Therefore, combining footprint modeling and remote sensing indices such as FTVT and NDVI is a feasible approach to assess and to understand the spatial representativeness of EC flux measurements.
Large-scale FLUXNET-driven model-observation comparisons are becoming very common and some advances have been made in assessment of surface characteristics near flux towers [11,12,[23][24][25][26].However, in the view of data users, it is still not clear that whether the spatial representativeness of flux measurement is suitable for their own application.Because the current assessments are not focused on a certain application or spatial level.Besides, although the number of flux sites has been increased remarkably, the spatial representativeness of EC flux measurements in China has not yet been fully considered so far.
In this study, based on different applications of EC flux measurements, two approaches to determine the spatial representativeness of EC flux measurements were proposed: (1) based on what fraction of the flux measurement is related to the target vegetation type, which is vegetation-type level; and (2) the spatial scope that the flux measurement can be extended to, which is spatial-scale level.We select FTVT and NDVI as key indicators to assess the spatial representativeness of 15 EC sites across northern China based on remote sensing analysis and footprint modeling.Here, the key questions are: (1) how well does the tower represent flux footprint from the targeted vegetation type; and (2) how representative is the flux tower footprint of the broader landscape and regional extents (in terms of NDVI).Furthermore, the spatial scope varies and can be divided into three sub-scales for spatial-scale level: (1) flux site scale, ranging from 0 to 1 km 2 [2,27,28]; (2) remote sensing pixel scale, which ranges from 1 km 2 to 9 km 2 , is most pixel sizes of moderate resolution satellite remote sensing images [29][30][31][32][33]; and (3) land model grid scale, which ranges from 9 km 2 to 324 km 2 , is roughly from 0.01 • × 0.01 • to 0.2 • × 0.2 • , was set to be comparable with the scale of regional or global land surface models [6,34,35].Our object is to assess the EC flux data in different aspects and scales, to understand the variation of spatial representativeness of EC flux measurements, to further reduce the uncertainties in up-scaling and improve the spatial representativeness.

Flux Sites and Data
The EC flux data in growing season (June to September 2009) were collected from 15 sites under a Coordinated Observation of Land-Atmospheric Interactions project in northern China [36].Much of the area is under arid and semi-arid climate, and is influenced by a west-east humidity gradient associated with variation in precipitation.These 15 flux sites represent the dominant vegetation/land cover types in the region: temperate grassland, cropland, deciduous broadleaf forests, and evergreen needle leaf forests (Figure 1).Characteristics of EC flux towers are summarized in Table 1.
the variation of spatial representativeness of EC flux measurements, to further reduce the uncertainties in up-scaling and improve the spatial representativeness.

Flux Sites and Data
The EC flux data in growing season (June to September 2009) were collected from 15 sites under a Coordinated Observation of Land-Atmospheric Interactions project in northern China [36].Much of the area is under arid and semi-arid climate, and is influenced by a west-east humidity gradient associated with variation in precipitation.These 15 flux sites represent the dominant vegetation/land cover types in the region: temperate grassland, cropland, deciduous broadleaf forests, and evergreen needle leaf forests (Figure 1).Characteristics of EC flux towers are summarized in Table 1.
Raw data of eddy covariance were recorded at a high rate of 10 Hz by a fast response data logger (Model CR5000, Campbell Scientific Inc., Logan, UT, USA) from all the 15 sites for this study.Data processing of all the flux records were performed with EdiRe software (Robert Clement, University of Edinburgh, UK), which include removing spikes from raw data, coordinate rotation, and frequency corrections.The average fraction of missing data due to gaps in the raw data was below 10% (minimum 2.1%, and maximum 13.3%).* "MAP" is the abbreviation of "Mean Annual Precipitation".
Raw data of eddy covariance were recorded at a high rate of 10 Hz by a fast response data logger (Model CR5000, Campbell Scientific Inc., Logan, UT, USA) from all the 15 sites for this study.Data processing of all the flux records were performed with EdiRe software (Robert Clement, University of Edinburgh, UK), which include removing spikes from raw data, coordinate rotation, and frequency corrections.The average fraction of missing data due to gaps in the raw data was below 10% (minimum 2.1%, and maximum 13.3%).

Remote Sensing Data
Cloud-free multi-spectral data of Landsat 5 TM (Thematic Mapper) at 30 m resolution were standard L1t files (georegistered, orthorectified) acquired over the similar period of coordinated observation with 18 km × 18 km area centered on each EC tower (Table 2) from the U.S. Geological Survey.The digital number (DN) data of these images was converted to top-of-atmosphere (TOA) reflectance, and further atmospherically corrected to land surface reflectance.

Representativeness of Flux Measurements
For the vegetation-type level, FTVT was calculated from land cover map and PDF.PDF was calculated as footprint climatology.Land cover classification was derived from Landsat-5 TM bands.According to field survey and local land use records, supervised classification (maximum likelihood) was performed to acquire vegetation cover classification around flux sites.The results of the classifications were validated with the areas surveyed in field or high spatial resolution images from Google Earth.In general, the accuracy of classification can be accepted (kappa index ranged from 0.79 to 0.98).Masks were then generated from land cover maps to screen out other land cover or vegetation types, and FTVT were further weighted by PDF.The higher FTVT means the better representativeness in vegetation-type level.
For spatial-scale level, based on adopting NDVI as a surrogate of land surface flux, spatial representativeness was expressed as bias of NDVI.The function of bias (δ) was proposed by Schmid and Lloyd [37] and further modified by Chen et al. [23] which was defined as to what extent the PDF-weighted NDVI in the source area matches with the un-weighted NDVI in a certain spatial scope of interest: where λ Φ is PDF-weighted NDVI of the source area.λ is un-weighted NDVI, which is calculated from the average value of NDVI over a rectangle window with an increasing size from 0.03 km to 18 km at 0.03 km step.NDVI is calculated from spectral reflectance in red and near infrared bands [38,39]: where R nir is the reflectance within the near-infrared (Band 4 in TM) and R red is the reflectance within the red band (Band 3 in TM).The lower bias of NDVI means the better representativeness in spatial-scale level.

Footprint and Footprint Climatology
EC flux measurements are averaged by footprint of upwind surface flux.Flux footprint describes the contribution of each element in the upwind source area to the measured flux.In this study, we calculated the flux footprints with an Eulerian analytical approach [40] at 30-min time steps at a pixel resolution of 30 m × 30 m (consistent with the Landsat 5).This approach assumes stand power law to obtain vertical profile of the horizontal wind velocity and the eddy diffusivity.Linkage of power-law profiles to Monin-Obukhov similarity profiles leads to an analytical solution of the flux footprint, the upwind distance, and the diffusion height.The flux footprint of EC flux measurement, f(x,y,z m ) was expressed as follows: where x is the downwind distance pointing against the average horizontal wind direction, y is the crosswind wind distance, z m is the measurement height, f y (x, z m ) is the crosswind integrated footprint, and D y (x, y) is the Gaussian crosswind distribution function of the lateral dispersion [41].
Then, the footprints for each site were accumulated to footprint climatology during the entire coordinated observation period (June-September 2009).The accumulated flux footprint was normalized by the cumulative footprint area of flux measurements to yield the footprint climatology contribution for each pixel as: where Ω Π is the cumulative footprint area [23].As the contribution of flux signals from each element within the source area is not equal, a weighted distribution map was generated from footprint climatology.Then, according to this map, FTVT and NDVI were weighted for each site.

Classification of EC Flux Measurements
Four grades of representativeness for the flux measurements were proposed from high to low in order: (1) homogeneous measurements; (2) representative measurements; (3) acceptable measurements; and (4) disturbed measurements.For vegetation-type level, grades were ranked by the PDF-weighted value of FTVT: 95% or more is homogeneous, 80% to 95% is representative, 50% to 80% is acceptable, and less than 50% is disturbed measurements [12].For spatial-scale level, grades were ranked by cluster analysis with k-means algorithm that classified the 15 sites into four groups based on changing trends and magnitudes.The ranking was taken in the 3 sub-scales: (1) flux site scale, ranged from 0 to 1 km 2 ; (2) remote sensing pixel scale, which ranges from 1 km 2 to 9 km 2 ; and (3) land model grid scale, which ranges from 9 km 2 to 324 km 2 .

Footprint Climatology
Contour lines (60%, 80%, 90%, 95% and 99%) of footprint climatology were overlaid on 3 km × 3 km color composite (Landsat TM bands 7, 4 and 2) rectangle images centered at each flux site (Figure 2).The composite images roughly illustrated the surface features such as land cover and topography around the sites.The footprints of the 15 sites differed substantially from each other in size and shape.For most sites, footprint climatology was distributed in an irregular shape, except for DS and MQ, where footprints were distributed symmetrically around the EC instruments and formed a roughly circular shape (Figure 2).Wind directions near the flux sites played an important role in shaping the footprint distribution.For each site, spatial distribution of prevailing wind direction presented similar spatial pattern as footprint climatology (Figure 3).

Footprint Climatology
Contour lines (60%, 80%, 90%, 95% and 99%) of footprint climatology were overlaid on 3 km × 3 km color composite (Landsat TM bands 7, 4 and 2) rectangle images centered at each flux site (Figure 2).The composite images roughly illustrated the surface features such as land cover and topography around the sites.The footprints of the 15 sites differed substantially from each other in size and shape.For most sites, footprint climatology was distributed in an irregular shape, except for DS and MQ, where footprints were distributed symmetrically around the EC instruments and formed a roughly circular shape (Figure 2).Wind directions near the flux sites played an important role in shaping the footprint distribution.For each site, spatial distribution of prevailing wind direction presented similar spatial pattern as footprint climatology (Figure 3).The area of footprint climatology showed exponential increase with the increasing cumulative footprint percentages (Figure 4).The area of 99% of footprint climatology varied from 0.52 km 2 to 1.65 km 2 , and 60% was limited in the range from 0.0048 km 2 to 0.23 km 2 , indicating that the signals captured by EC sensor was mainly concentrated around areas near to the sensor location.Areas of footprint climatology were mainly controlled by the effective height of EC sensor and surface roughness.DYK, MY, and GT, the three sites with higher effective height than others, also had much larger source areas than others.The areas of DYK and MY, two forest sites, were 1.21 km 2 and 1.56 km 2 , respectively, while larger footprint was also found at cropland site GT with a sensor height of 15.6 m.The high altitude of sensor and relatively low surface roughness resulted in the GT site having the largest footprint climatology with an area of 1.65 km 2 .The area of footprint climatology showed exponential increase with the increasing cumulative footprint percentages (Figure 4).The area of 99% of footprint climatology varied from 0.52 km 2 to 1.65 km 2 , and 60% was limited in the range from 0.0048 km 2 to 0.23 km 2 , indicating that the signals captured by EC sensor was mainly concentrated around areas near to the sensor location.Areas of footprint climatology were mainly controlled by the effective height of EC sensor and surface roughness.DYK, MY, and GT, the three sites with higher effective height than others, also had much larger source areas than others.The areas of DYK and MY, two forest sites, were 1.21 km 2 and 1.56 km 2 , respectively, while larger footprint was also found at cropland site GT with a sensor height of 15.6 m.The high altitude of sensor and relatively low surface roughness resulted in the GT site having the largest footprint climatology with an area of 1.65 km 2 .

Spatial Representativeness at Vegetation-Type Level
At vegetation-type level, the average PDF-weighted FTVT was 0.83 and the PDF-weighted FTVT ranged from 0.29 to 1 (Table 3).Eight flux sites were greater than 0.95 (green cells), which were ranked into homogeneous measurements and four flux sites were greater than 0.8 (yellow cells), which were ranked into representative measurements.As to the other three sites, DYK was greater than 0.5 (blue cell), which was ranked into acceptable measurement, and CW and MY (red cells) were less than 0.5, which were ranked into disturbed measurements.As a contrast, in un-weighted FTVT classification, many of these sites had lower value of FTVT than the weighted ones.The number of homogeneous and representative measurements decreased to three and four, respectively, while both the number of acceptable and disturbed measurements increased to four.Table 3. Spatial representativeness of the sites at vegetation-type level and spatial-scale level, ranked by fraction of target vegetation type (FTVT) and cluster analysis of NDVI biases (%), respectively.The color of each cell represents its grades: green cell is homogeneous measurements; yellow cell is representative measurements; blue cell is acceptable measurements; and red cell is disturbed measurements.

Spatial Representativeness at Spatial-Scale Level
At spatial-scale level, bias of NDVI for each site was repetitively generated and collected (Figure 5).The reference areas were then divided into three sub-scales: flux site scale, remote sensing pixel scale, and land model grid scale.For each sub-scale, all 15 sites were ranked into four grades based on cluster analysis.

Spatial Representativeness at Vegetation-Type Level
At vegetation-type level, the average PDF-weighted FTVT was 0.83 and the PDF-weighted FTVT ranged from 0.29 to 1 (Table 3).Eight flux sites were greater than 0.95 (green cells), which were ranked into homogeneous measurements and four flux sites were greater than 0.8 (yellow cells), which were ranked into representative measurements.As to the other three sites, DYK was greater than 0.5 (blue cell), which was ranked into acceptable measurement, and CW and MY (red cells) were less than 0.5, which were ranked into disturbed measurements.As a contrast, in un-weighted FTVT classification, many of these sites had lower value of FTVT than the weighted ones.The number of homogeneous and representative measurements decreased to three and four, respectively, while both the number of acceptable and disturbed measurements increased to four.Table 3. Spatial representativeness of the sites at vegetation-type level and spatial-scale level, ranked by fraction of target vegetation type (FTVT) and cluster analysis of NDVI biases (%), respectively.The color of each cell represents its grades: green cell is homogeneous measurements; yellow cell is representative measurements; blue cell is acceptable measurements; and red cell is disturbed measurements.At spatial-scale level, bias of NDVI for each site was repetitively generated and collected (Figure 5).The reference areas were then divided into three sub-scales: flux site scale, remote sensing pixel scale, and land model grid scale.For each sub-scale, all 15 sites were ranked into four grades based on cluster analysis.In flux site scale, the average of absolute NDVI bias was 4.34%, and nine sites were categorized as homogeneous measurements due to their lowest biases among four groups and the decreasing trends of biases with increasing reference areas.LZ, NM, TYG and YK were ranked into representative measurements.Biases of NDVI at JZ and SPT showed large amplitude of changes and tended to increase with window size.Therefore, JZ and SPT were ranked into acceptable and disturbed measurement, respectively (Figure 6a).In remote sensing pixel scale, the average of absolute NDVI bias was 8.27%, seven sites were ranked into homogeneous measurements due to their lowest biases in the four grades.GT and YZ also had lower biases, but their variations of biases trended to increase, therefore, these two sites were ranked into representative measurements.The remaining six sites had larger biases, NM and TYG were ranked into acceptable measurements and JZ, LZ, SPT and YK were ranked into disturbed measurements (Figure 6b).In land model grid scale, the average of absolute NDVI bias was 12.13%, five sites (CW, DS, MY, NM and TYG) were ranked into homogeneous measurements due to their lowest biases and steady trends.Although GT, MQ, TYC and SPT had larger variations of biases than the five sites above, their variations trended to be steady and therefore ranked into representative measurements.Based on the variations and the changing trends of biases, AR, DYK and YZ were ranked into acceptable measurements and JZ, LZ and YK were ranked into disturbed measurements (Figure 6c).

Weighted
In general, five sites (CW, DS, MQ, MY and TYC) kept high spatial representativeness in all three sub-scales.Six sites (AR, DYK, GT, NM, TYG and YZ) had high spatial representativeness in a certain spatial scale.While the remaining four sites (JZ, LZ, SPT and YK) present low spatial representativeness at all sub-scales (Table 3).With the increasing area of spatial scope, spatial representativeness of the sites decreased (Table 4).In flux site scale, the average of absolute NDVI bias was 4.34%, and nine sites were categorized as homogeneous measurements due to their lowest biases among four groups and the decreasing trends of biases with increasing reference areas.LZ, NM, TYG and YK were ranked into representative measurements.Biases of NDVI at JZ and SPT showed large amplitude of changes and tended to increase with window size.Therefore, JZ and SPT were ranked into acceptable and disturbed measurement, respectively (Figure 6a).In remote sensing pixel scale, the average of absolute NDVI bias was 8.27%, seven sites were ranked into homogeneous measurements due to their lowest biases in the four grades.GT and YZ also had lower biases, but their variations of biases trended to increase, therefore, these two sites were ranked into representative measurements.The remaining six sites had larger biases, NM and TYG were ranked into acceptable measurements and JZ, LZ, SPT and YK were ranked into disturbed measurements (Figure 6b).In land model grid scale, the average of absolute NDVI bias was 12.13%, five sites (CW, DS, MY, NM and TYG) were ranked into homogeneous measurements due to their lowest biases and steady trends.Although GT, MQ, TYC and SPT had larger variations of biases than the five sites above, their variations trended to be steady and therefore ranked into representative measurements.Based on the variations and the changing trends of biases, AR, DYK and YZ were ranked into acceptable measurements and JZ, LZ and YK were ranked into disturbed measurements (Figure 6c).
In general, five sites (CW, DS, MQ, MY and TYC) kept high spatial representativeness in all three sub-scales.Six sites (AR, DYK, GT, NM, TYG and YZ) had high spatial representativeness in a certain spatial scale.While the remaining four sites (JZ, LZ, SPT and YK) present low spatial representativeness at all sub-scales (Table 3).With the increasing area of spatial scope, spatial representativeness of the sites decreased (Table 4).

Assessment of Spatial Representativeness
In general, spatial representativeness of flux measurements from different sites varied among different application levels and scales.Most of the sites in this study had high spatial representativeness in limited spatial scale or in a certain vegetation type.Sites with high representativeness at vegetation-type level did not necessarily have high representativeness at spatial-scale level.For example, CW and MY are sites with disturbed measurements at vegetation-type level, but they are homogeneous measurements at spatial-scale level.
At vegetation-type level, grass and crop sites had higher FTVT than forest sites.For these sites, dominant vegetation types in the source areas were the specific target vegetation types and therefore the flux measurements represent the vegetation types well.Towers of forest sites were mounted on complex surfaces containing different vegetation types (Figure 7).Flux measurements from these sites were mixed by different vegetation types and therefore they were considered not suitable for flux studies on a specific target vegetation type.However, they may not be set up for observing a specific vegetation type, but for a large scale observation with mosaic land cover types [24], as their source areas had no dominant vegetation type [42].

Assessment of Spatial Representativeness
In general, spatial representativeness of flux measurements from different sites varied among different application levels and scales.Most of the sites in this study had high spatial representativeness in limited spatial scale or in a certain vegetation type.Sites with high representativeness at vegetation-type level did not necessarily have high representativeness at spatial-scale level.For example, CW and MY are sites with disturbed measurements at vegetationtype level, but they are homogeneous measurements at spatial-scale level.
At vegetation-type level, grass and crop sites had higher FTVT than forest sites.For these sites, dominant vegetation types in the source areas were the specific target vegetation types and therefore the flux measurements represent the vegetation types well.Towers of forest sites were mounted on complex surfaces containing different vegetation types (Figure 7).Flux measurements from these sites were mixed by different vegetation types and therefore they were considered not suitable for flux studies on a specific target vegetation type.However, they may not be set up for observing a specific vegetation type, but for a large scale observation with mosaic land cover types [24], as their source areas had no dominant vegetation type [42].At spatial-scale level, spatial representativeness generally decreased along with the increasing area of interest.This is consistent with the results of Chen et al. [23].An exceptional site is SPT, which At spatial-scale level, spatial representativeness generally decreased along with the increasing area of interest.This is consistent with the results of Chen et al. [23].An exceptional site is SPT, which is a steppe desert site near to the Yellow River (Figure 7).With the increasing area of interest, the fraction of water body decreased and the influence of water body on flux measurements was also weakened, therefore, the representativeness of this site trended to increase (Figures 7 and 8).Besides, zonal vegetation sites had greater representativeness than non-zonal vegetation sites.DS and MQ are zonal vegetation sites with desert steppe and sub-alpine meadow steppe respectively.JZ, LZ and YK are sites with non-zonal vegetation type.LZ and YK are crop sites surrounded by desert in oasis ecosystem, while JZ is a crop site near urban area.There is a huge contrast of vegetation greenness between cropland and bare soil or urban area.As to CW, MY and TYC, although surroundings near to the tower were mixed with different vegetation types, but these vegetation types had similar greenness and vegetation indices, therefore these sites still had high representativeness.The similarity of vegetation activity surrounds the sensor, expressed as NDVI in this study, play an important role in determining the representativeness of site (Figures 8 and 9).For each site in this study, weighted and un-weighted NDVI was compared (Table 5).NDVI weighted by footprint climatology were closer to real flux measurement than un-weighted NDVI.Areas near to towers had greater fraction of target vegetation, and therefore were more important than other areas.The difference between weighted and un-weighted NDVI were less than 10% for most sites except for JZ, MY, NM, SPT and TYC.These five sites had greater complexity of vegetation types, were more heterogeneous in NDVI distribution, and were also considered to have greater uncertainties in flux up-scaling.Caution should be taken into these sites because the uncertainty in long term flux measurements can be introduced by the footprint variation.is a steppe desert site near to the Yellow River (Figure 7).With the increasing area of interest, the fraction of water body decreased and the influence of water body on flux measurements was also weakened, therefore, the representativeness of this site trended to increase (Figures 7 and 8).Besides, zonal vegetation sites had greater representativeness than non-zonal vegetation sites.DS and MQ are zonal vegetation sites with desert steppe and sub-alpine meadow steppe respectively.JZ, LZ and YK are sites with non-zonal vegetation type.LZ and YK are crop sites surrounded by desert in oasis ecosystem, while JZ is a crop site near urban area.There is a huge contrast of vegetation greenness between cropland and bare soil or urban area.As to CW, MY and TYC, although surroundings near to the tower were mixed with different vegetation types, but these vegetation types had similar greenness and vegetation indices, therefore these sites still had high representativeness.The similarity of vegetation activity surrounds the sensor, expressed as NDVI in this study, play an important role in determining the representativeness of site (Figures 8 and 9).For each site in this study, weighted and un-weighted NDVI was compared (Table 5).NDVI weighted by footprint climatology were closer to real flux measurement than un-weighted NDVI.Areas near to towers had greater fraction of target vegetation, and therefore were more important than other areas.The difference between weighted and un-weighted NDVI were less than 10% for most sites except for JZ, MY, NM, SPT and TYC.These five sites had greater complexity of vegetation types, were more heterogeneous in NDVI distribution, and were also considered to have greater uncertainties in flux up-scaling.Caution should be taken into these sites because the uncertainty in long term flux measurements can be introduced by the footprint variation.The variation of spatial representativeness implicates the importance of understanding spatial representativeness of flux measurements before use.Flux data should be assessed first to determine whether the measurements could be used to reflect actual surface condition from a certain spatial scale or a specific vegetation type [11,43].Data users should also be aware of the spatial representativeness when applying flux measurements.

Implications for Up-Scaling
Surface flux can be up-scaled by integrating flux measurement with remote sensing pixel or land model grid.With the development of algorithms, remote sensing and land surface model are becoming main method for acquiring regional or even global scale distribution of carbon, water and energy flux.EC flux measurements play an important role in calibration and validation of these models.However, EC flux data from a single set cannot directly be compared with them due to the influences of heterogeneity on flux measurements [44,45].Mismatching between local scale EC flux observations and regional measurements can increase uncertainties in surface flux up-scaling.Hence, source area of EC flux should be involved in model validation and/or calibration with EC flux measurements [46].Fu et al. (2014) have up-scaled sixteen flux sites from the Canadian Carbon Program and AmeriFlux networks located in North America to landscape and regional scales based on an improved method which taking the PDF-weighted flux measurements of EC tower into account [47].Liu et al. (2016) applied a combined method which introducing source area distribution and land surface temperature to up-scale evapotranspiration measurements from site scale to satellite pixel scale [48].
Meanwhile, spatial representativeness of flux measurements can be improved by applying a footprint-based filter to improve their spatial representativeness over large area.This is extremely important for sites in monsoon region due to their sharply different wind directions in different seasons.Sánchez et al. (2010) improved the energy closure by weighting flux with target vegetation type in a FLUXNET boreal forest site of Finland [49].They reported the energy closure is related to wind direction when the land surface is heterogeneous.Spatial weighting by footprint climatology can improve the usability of flux measurements by introducing more target vegetation information extracted from the areas near the tower.Barcza et al. (2009) extracted crop-specific carbon dioxide fluxes from tall tower EC measurements under mixed land cover types.They found the extracted fluxes exhibited strong covariance with crop-specific NDVI time series [24].

Limitation of Remote Sensing Images
Remote sensing images play an important role in assessment of spatial representativeness by providing important surface vegetation information such as vegetation types and vegetation indices.Although different spatial resolution images ranging from 1 m to 250 m were adopt [11,46,47], Landsat images were still considered as the most suitable remote sensing data set for global coverage and free of charge [23].More remote sensing images of the study period will benefit the assessment by capturing detailed seasonality of underlying vegetation types.However, the availability of Landsat images was limited due to the 16-day revisit-cycle and cloud contamination.Our work relied on the assumption that seasonalities of surface flux for different vegetation types within the source area are similar enough.One solution to address the issue is to generate more images based on image fusion.Fu et al. [47] had generated high temporal-spatial resolution net ecosystem exchange (NEE) based on Enhanced Spatial and Temporal Adaptive Reflectance Fusion Model (ESTARFM).Images from Landsat and Moderate Resolution Imaging Spectroradiometer (MODIS) were merged to generate Landsat-like images with both high spatial (30 m) and temporal (8-day) resolution [47,50].

Conclusions
Assessment of spatial representativeness of EC flux measurements is important for understanding the effects of surface heterogeneity around flux towers on flux up-scaling.However, the current assessments were not focused on a certain application or spatial level.The users are still not clear whether the spatial representativeness of flux measurement is suitable for their own application.In this study, we integrated satellite derived vegetation information with footprint modeling to calculate and assess the spatial representativeness of 15 EC flux sites across northern China at vegetation-type level and spatial-scale level respectively.These sites were then ranked in four grades from high to low in order: (1) homogeneous; (2) representative; (3) acceptable; and (4) disturbed measurements based on FTVT or cluster analysis.
In general, spatial representativeness of flux measurements from different sites varied among different application levels and scales.Most of the sites had high spatial representativeness in limited spatial scale or in a certain vegetation type.Sites with high representativeness at vegetation-type level did not necessarily have high representativeness at spatial-scale level.At vegetation-type level, FTVT ranged from 0.29 to 1 and the average FTVT was 0.83.Grass and crop sites had higher representativeness than forest sites.At spatial-scale level, spatial representativeness generally decreased along with the increasing area of interest and the scales were further divided into three sub-scales: flux site scale, remote sensing pixel scale and land model grid scale.Five sites kept high spatial representativeness in all three sub-scales.Six sites had high spatial representativeness in a certain spatial scale, and four sites present low spatial representativeness at all sub-scales.Flux sites with zonal vegetation had greater representativeness than non-zonal vegetation sites.
These results illustrated the heterogeneity around the EC flux tower and highlighted the importance of proper interpretation of EC flux measurements.These results also suggest that source area of EC flux should be involved in model validation and/or calibration with EC flux measurements.Spatial representativeness of flux measurements can be also improved by applying a footprint-based filter.To make a more comprehensive assessment, future efforts should be paid in generating more high temporal-spatial resolution images to better reflect seasonal changes of surface flux in complex underlying surfaces around flux towers.

Figure 1 .
Figure 1.Geographical location of the study area and the 15 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.

Figure 1 .
Figure 1.Geographical location of the study area and the 15 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.

Figure 2 .
Figure 2. Color composite images (bands 7, 4 and 2) of Landsat at a 30 m resolution for each area (3 km × 3 km) centered at individual flux sites and the contours of cumulative annual footprint climatology of the observation period.The keys of the contour line for percentage of total footprint are: purple, 60%; cyan, 80%; green, 90%; blue, 95%; and red, 99%.

Figure 2 .
Figure 2. Color composite images (bands 7, 4 and 2) of Landsat at a 30 m resolution for each area (3 km × 3 km) centered at individual flux sites and the contours of cumulative annual footprint climatology of the observation period.The keys of the contour line for percentage of total footprint are: purple, 60%; cyan, 80%; green, 90%; blue, 95%; and red, 99%.

Figure 3 .
Figure 3. Rose charts of wind direction frequency of the 15 flux sites in this study.For each site, spatial distribution of prevailing wind direction presented similar spatial pattern as footprint climatology (Figure 2).

Figure 3 .
Figure 3. Rose charts of wind direction frequency of the 15 flux sites in this study.For each site, spatial distribution of prevailing wind direction presented similar spatial pattern as footprint climatology (Figure 2).

Figure 4 .
Figure 4. Cumulative footprint and their corresponding areas for the 15 eddy covariance (EC) sites.

Figure 4 .
Figure 4. Cumulative footprint and their corresponding areas for the 15 eddy covariance (EC) sites.

Figure 5 .
Figure 5. Variations of biases of Normalized Difference Vegetation Index (NDVI) along with the increasing area of reference.

Figure 5 .
Figure 5. Variations of biases of Normalized Difference Vegetation Index (NDVI) along with the increasing area of reference.

Figure 7 .
Figure 7. Land cover classification of the 15 sites at a 30 m resolution over grid area of 3 km × 3 km centered at individual flux sites.

Figure 7 .
Figure 7. Land cover classification of the 15 sites at a 30 m resolution over grid area of 3 km × 3 km centered at individual flux sites.

Figure 8 .
Figure 8. Land cover maps of the 15 sites at a 30 m resolution for each area (18 km × 18 km) centered at individual flux sites.

Figure 8 .Figure 9 .
Figure 8. Land cover maps of the 15 sites at a 30 m resolution for each area (18 km × 18 km) centered at individual flux sites.

Figure 9 .
Figure 9. NDVI maps of the 15 sites at a 30 m resolution over grid of 18 km × 18 km centered at individual flux sites.

Table 1 .
Main characteristics of the 15 flux sites in the study region.

Table 2 .
Basic information of the acquired images (Landsat-5 TM).

Table 4 .
Number of sites in different spatial-scale level and different grades.

Table 4 .
Number of sites in different spatial-scale level and different grades.

Table 5 .
Weighted by footprint and un-weighted NDVI.

Table 5 .
Weighted by footprint and un-weighted NDVI.