Area-averaged evapotranspiration over a heterogeneous land surface : Aggregation of multi-point EC flux measurements with high-resolution landcover map and footprint analysis

The determination of area-averaged evapotranspiration (ET) at the satellite pixel scale/model grid scale over a heterogeneous land surface plays a significant role in developing and improving the parameterization schemes of the remote sensing based ET estimation models and general hydro-meteorological models. The Heihe Watershed Allied Telemetry Experimental Research (HiWATER) flux matrix provided a unique opportunity to build an aggregation scheme for area-averaged fluxes. On the basis of HiWATER flux matrix dataset and high-resolution land-cover map, this study focused 15 on estimating the area-averaged ET over a heterogeneous landscape with footprint analysis and multivariate regression. The procedure is as follows: Firstly, quality-control and uncertainty-estimation for the data of the flux matrix, including 17 eddy-covariance (EC) sites and 4 groups of large aperture scintillometer (LAS), were carefully done. Secondly, the representativeness of each EC site was quantitatively evaluated; footprint analysis was also performed for each LAS path. Thirdly, based on the high-resolution land-cover map derived from aircraft remote sensing, a flux aggregation method was 20 established combining footprint analysis and multiple-linear regression. Then, the area-averaged sensible heat fluxes obtained from the EC flux matrix were validated by the LAS measurements. Finally, the area-averaged ET of the kernel experimental area of HiWATER was estimated. Compared with the formerly used and rather simple approaches, such as the arithmetic average and area-weighted methods, present scheme is not only with a much better database but also has a solid grounding in physics and mathematics in the integration of area-averaged fluxes over a heterogeneous surface. Results from 25 2 this study, both instantaneous and daily ET at the satellite pixel scale, can be used for the validation of relevant remote sensing models and land surface process models. Furthermore, this work will be extended to the water balance study of the whole Heihe River basin.


Introduction
Land surface evapotranspiration (ET) is not only a key component in the regional water circulation, but is also essential in the surface energy balance and land surface process.Under the condition of increasing shortage of water resources, high-precision estimation of ET at regional scale is essential for such applications, as the management of river basin water resources, regional planning, and the sustainable development of agriculture (Wang et al., 2003).Currently, the commonly used methods for acquisition of regional ET are ground-based observation, remote sensing based estimation, and model simulation, respectively.The Earth's surface is always characterized by spatial heterogeneity.Large land surface heterogeneity affects greatly the exchanges of momentum, heat, and water between the Published by Copernicus Publications on behalf of the European Geosciences Union.
F. Xu et al.: Aggregation of multi-point EC flux measurements land surface and atmosphere (Mengelkamp et al., 2006).Indeed, the surface heterogeneity caused either by the contrast in soil moisture or vegetation type generates a large spatial variability of fluxes, which limit the use of the eddycovariance (EC) system, unless one deploys a network of EC devices (Ezzahar et al., 2009b).The flux tower group can quantify the turbulent exchange of energy and mass between the atmosphere and a variety of surface types (Sellers et al., 1995), and these local point measurements need to be aggregated to provide meaningful area-averaged fluxes (André et al., 1986).If special aggregation rules for local flux measurements are applied, measurements can provide averaged fluxes at model grid scale (Beyrich et al., 2006;Mahrt et al., 2001).But given the EC network's high price and the requirement for their continuous maintenance, the largeaperture scintillometer (LAS) is a useful alternative method for direct measurements of area-averaged sensible heat fluxes on the scale of 1-5 km (Ezzahar et al., 2009b;Ezzahar and Chehbouni, 2009).
Satellite has been considered a promising data source for deriving regional ET with the development of remote sensing techniques (Ezzahar et al., 2009a).In response to increasing demand for spatially distributed hydrologic information, many satellite-based approaches have been developed for routine monitoring of ET at a regional scale (Anderson et al., 2012).Nevertheless, the effectiveness of the remote sensing based methods for estimating ET must be fully assessed by ground-based area-averaged flux measurements, mainly due to the uncertainties of model inputs and parameterization schemes (Wang et al., 2003).Furthermore, there may be a bias in directly comparing a remote sensing based ET estimation with in situ measurements, because of their spatial-scale mismatch and spatial heterogeneity at the subpixel scale (Jia et al., 2012).
General atmospheric-hydrological models (e.g., numerical weather prediction) can adequately describe the interaction between the atmosphere and the underlying surface using complex parameterization schemes.The development and validation of these models are usually based on measurements performed over homogeneous land surfaces.While the assumption of homogeneity might be justified at the local scale (10-10 3 m), it is often violated at the scale of the grid resolution of current regional atmospheric models (about 10 4 m) (Beyrich et al., 2006;Beyrich and Mengelkamp, 2006).Therefore, it is significantly important to determine the area-averaged surface fluxes at the satellite pixel scale/model grid scale (10 3 -10 4 m) for the evaluation of general hydro-meteorological models and relevant remote sensing models.The estimation of area-averaged fluxes is usually not straightforward, especially for heterogeneous land surfaces.
A number of international field experiments have been performed over heterogeneous land surfaces in different geographical and climate regions of the Earth in recent decades (Mengelkamp et al., 2006;Beyrich et al., 2006;Wang, 1999), such as HAPEX-MOBILHY (André et al., 1986), FIFE (Sellers et al., 1988), HAPEX-SAHEL (Goutorbe et al., 1994), BOREAS (Sellers et al., 1995), NOPEX (Halldin et al., 1998), and LITFASS-2003(Mengelkamp et al., 2006).In these experiments, surface fluxes at the model grid scale, estimated from multi-point flux observations using various flux aggregation techniques, were compared with those obtained from LAS systems and remote sensing estimation methods.In the former studies, the most commonly used and rather simple flux aggregation methods mainly include the arithmetic average method, the area-weighted method, and the footprint-weighted method (Liu et al., 2016).These studies revealed, under careful data processing and quality control (Charuchittipan et al., 2014) as well as analysis of the energy balance closure for flux data (Foken et al., 2006(Foken et al., , 2010)), that the integration of the multi-site EC flux measurements and area-averaged fluxes from scintillometers and aircraft observations can provide reasonable estimates over a heterogeneous landscape (Mahrt et al., 2001;Beyrich et al., 2006;Liu et al., 2016).
However, the integration schemes of the aforementioned methods are applicable for relatively uniform sites, of which the local flux measurements are representative of the individual surface types.For the interpretation of tower flux measurements over a heterogeneous land surface, operational footprint analysis is an essential approach (Schmid, 2002).The development of footprint models provides diagnostic tools to quantify the representativeness of tower flux measurements for selected sites (Horst and Weil, 1992;Kim et al., 2006).Besides, it had been demonstrated that the footprint climatology can be combined with information on the spatial variability of vegetation types provided by satellite images (Kim et al., 2006;Chen et al., 2008).Land cover reflects the combined effects of vegetation, climate, soil, and topography; some relationship could be found between land cover and measured surface fluxes (Ogunjemiyo et al., 2003).Ran et al. (2016) proposed four indicators with footprint analysis and a land-cover map to improve the representativity of EC towers and correct the EC flux measurements.But this method did not obtain the surface fluxes of individual land-cover types, but just corrected the EC observations with some prior coefficients.Some previous studies have successfully related the aircraft observed fluxes to surface cover types with the integration of footprint analysis and a satellitebased land-cover map (Ogunjemiyo et al., 2003;Kirby et al., 2008;Hutjes et al., 2010).Among these works, a flux disaggregation method (Hutjes et al., 2010), developed from a former study presented by Ogunjemiyo et al. (2003), would be a promising method for integrating multiple tower-based flux measurements to satellite pixel or grid scale on account of its theoretical framework.The application of this method in attributing heterogeneous EC flux measurements to separate land-cover classes will hopefully be a way to have insight into the component fluxes from various land-cover types.It also provides a chance to develop a flux aggregation scheme for exploring the extension of multiple EC flux observations to satellite pixel/gird scale.
A multi-scale observation experiment on evapotranspiration over a heterogeneous land surface was conducted in the middle reaches of the Heihe River basin during the project of HiWATER (Heihe Watershed Allied Telemetry Experimental Research) in 2012 (Li et al., 2013;Liu et al., 2016).A comprehensive flux matrix, consisting of 17 EC sites and four groups of LAS systems within a 5 × 5 km 2 area, was specifically designed to capture the multi-scale characteristics of ET over a heterogeneous landscape during the experiment.The HiWATER flux matrix, with abundant multi-scale flux data, provided a unique opportunity to build an aggregation scheme for area-averaged fluxes over a heterogeneous land surface.The objective of this study is to integrate multi-point EC flux measurements to area-average with high-resolution land-cover map and footprint analysis.The main issues were as follows: (1) the representativeness of the EC flux matrix was quantitatively evaluated; (2) a flux aggregation scheme was established and used for estimating area-averaged sensible heat fluxes from the EC flux matrix, taking LAS measurements as a reference to check the integration algorithm; and (3) the developed scheme was applied to determine the area-averaged evapotranspiration over a heterogeneous land surface.
2 Study sites and data

Site description
This study was based on ground-based observation dataset, collected from the multi-scale flux matrix of HiWATER from May to September 2012.The kernel experimental area (5 × 5 km 2 ) of the multi-scale observation experiment was located in the Yingke and Daman irrigation district within Zhangye oasis.The land-cover types in the kernel experimental area were dominated by maize (72 %), vegetables (5 %), orchard and shelterbelt (8 %) as well as residential area and roads (15 %).The small square maize fields were staggered with windbreak trees, roads and irrigation ditches; the surface status of this oasis was actually very heterogeneous.As shown by the numbers 1-17 in Fig. 1, 17 sites were installed according to the distribution of crop planting structure and land cover.Each of them was equipped with an eddycovariance system (with two layers in site 15) and an auto- matic weather station (AWS) to capture the exchange process of surface water and energy budget at the local scale.Spatial distribution of EC/AWS systems is shown in Fig. 1, with site 1 of vegetable (pepper) field, site 4 of residential area, site 17 of apple orchard, and the other 14 sites are spatially distributed on the dominated maize fields.Key micrometeorological observations at each AWS included four-component radiation, one or two levels wind/temperature/relative humidity, soil temperature/moisture and soil heat flux.Among these sites, site 15 was a superstation equipped with two levels of EC system, seven-level wind speed/direction and air temperature/humidity. Four groups of large-aperture scintillometers were installed across the experimental district (see Fig. 1).Details of the EC and LAS systems in HiWATER flux matrix are given in Tables 1 and 2, respectively.
2.2 Data collection, processing and quality control

Flux data processing and quality control
Data on the typical clear days of 29-30 June 2012 were selected for the following analysis, including EC data from 16 towers (except site 3 and the highest level (34 m) of site 15) and four groups of LAS data as well as the multi-point micrometeorological data listed above.The last round of irrigation in each plot was done before 26 June; during the 2 days, there was almost no irrigation in the flux matrix.Firstly, AWS data sampled at 10 min were averaged to a 30 min period.
Careful data processing and quality control for EC and LAS raw data were then performed so as to ensure a high-quality flux dataset.The EddyPro software developed by LI-COR (Lincoln, Nebraska USA, www.licor.com/eddypro)was used to process the 10 Hz raw EC data into a half-hourly averaged flux data, by procedures including spike removal, angle of attack correction (for Gill), time lag correction, coordinate rotation (2-D rotation), frequency response correction, sonic virtual temperature correction, and corrections for density fluctuation (Webb-Pearman-Leuning, WPL).Data quality assessment was performed for the turbulent flux every 30 min using the flagging system with three different levels (0, 1, and 2) (Mauder and Foken, 2015).Detailed information on the processing steps can be found in Wang et al. (2015) and Xu et al. (2013).For this study, only the flux data of flag 0 (the best) were used.Flux data of flag 2, as well as the data at night when the friction velocity was below 0.1 m s −1 , were discarded (Blanken, 1998;Liu et al., 2011).To obtain daily ET, at first, a gap-filling method, based on the nonlinear regression (establishing the relationship between the latent heat flux and net radiation), for the 30 min latent heat fluxes (LE) was used.Then, the daily ET was calculated by summing the half-hourly gap-filled ET to 24 h totals.
For the EC systems used in the data analysis, we have tried to reduce the systematic errors to a minimum with a pre-observation inter-comparison, and careful maintenances during the observation period (Xu et al., 2013).The random errors were also analyzed by a separate research, which can be minimized in an ensemble average (Wang et al., 2015).The energy balance closure ratio (EBR) for the EC data of the flux matrix was also carefully assessed.Generally, the EBR during the 3 and half months was good.For the 17 EC stations in the intensive observation area, the average EBR was about 0.92.Except the lowest EBR (0.78) in orchard site, values in other sites were scattered without clear relation to the surface status.For site 15 with two heights of EC system, the relevant EBR were 0.89 (at 4.5 m) and 1.03 (for 34 m), respectively (Xu et al., 2017).For daytime conditions (global radiation Rg > 20 W m −2 ), the systematic error of the scalar fluxes (δ sys F ) can be quantified indirectly through the surface energy balance closure (EBR), which is defined as (Mauder et al., 2013).So the systematic error of the turbulent flux for HiWATER flux matrix was generally about 8 %.
The LAS system provided a measurement of the structure parameter for the refractive index of air (C 2 n ) with an output period of 1 min.The raw data were firstly checked, mainly rejecting the saturated cases when C 2 n < 0.193R −8/3 λ 1/3 D 5/3 (where R is the path length, D the optical aperture, and λ the wavelength) (Ochs and Wilson, 1993).Then, the data were averaged to 30 min, and the path-average sensible heat fluxes were iteratively calculated based on Monin-Obukhov similarity theory (MOST) (Andreas, 1988).The parameters used in this calculation, like the roughness height and zero-plane displacement, were obtained following Martano (2000); other parameters, including wind speed, air pressure and temperature, Obukhov length, and Bowen ratio, were directly obtained from relevant EC measurements.For this study, the sensible heat fluxes from the LAS systems at daytime (8:30-15:30, Beijing Standard Time, BST; the time difference between local time and BST is approximately +1 h 18 min) were selected.
As for the eddy-covariance systems, quality control for the flux data from the four groups of LAS was also done.The systematic errors from data processing, e.g., the larger effects of Bowen-ratio correction in this oasis area, were carefully minimized.We checked the sensible heat fluxes (H) from the four groups of LAS with that from the nearer ECs.Except for LAS3 (under its path there were clearly some village buildings, so the H_LAS is higher), others agreed very well with that of ECs.These are consistent with the other two studies in the Heihe River basin: Liu et al. (2011) have reported that the LAS measured sensible heat flux over an alpine meadow was consistent with that of EC, and Xu et al. (2013) have concluded that both the EC and LAS system measurements were reliable and comparable during the intensive observation period of HiWATER.

High-resolution land-cover map from remote sensing
Based on the airborne hyper-spectral images acquired by the Compact Airborne Spectrographic Imager (CASI) on 29 June 2012 and the Canopy Height Model (CHM) data from the lidar data collected on 9 July 2012, a land-cover classification map with 1 m spatial resolution was derived using an object-based classification method.This was done mainly for the kernel experimental area.The overall accuracy of the land-cover classification map is up to 90 %, and the Kappa coefficient is approximately 0.9.The detailed classification process of the map can be found in Liu and Bo (2015).Land-cover misclassification still occurred on this map because of spectral similarity, especially on the edges of different surface-cover types.To obtain a more accurate landcover map, the misclassified patches of the land-cover map were visually and manually revised, according to the highresolution Charge-Coupled Device (CCD) images (acquired on 26 July) and the Google Earth imagery (on 3 September 2012).Finally, for the aim of this study, the refined 12 kinds of land classification types, of which most were different vegetables of small areas, were merged into four kinds (maize, woods, vegetables, and non-vegetation types) in accordance with crop species and surface types, as shown in Fig. 1.Among the four land-cover types, the non-vegetation type mainly contains the buildings, roads, as well as plastic film covered greenhouses, while the woods type consists of orchard and shelterbelts.land-cover classes (Hutjes et al., 2010).
where F is the total flux of any scalar (e.g., the sensible and latent heat flux in our case) for a specified area, A k is the fractional coverage of an individual land-cover class k within that area, and F k is the flux from the individual land-cover class k; n is the number of land-cover classes that is distinguished in the specified area.
The observed flux (F obs ) at height z m can be closely related to the true surface flux from upwind measurement point through the footprint function, in continuous form (Leclerc and Foken, 2014): Here x obs , y obs are the site coordinates, z m is the effective observation height, defined as z m = z − d(where z is the sensor height, d the zero-plane displacement).The footprint function w (x, y, z m ) describes the flux portion seen at the sensor position (x obs , y obs , z m ).Equation ( 2) can be discretized for a uniform grid over a landscape, as in a land-cover classification map based on a remote sensing image.Leaving out the height dependence for simplification, it becomes where each pixel x y of the map is assumed to be homogeneous, which is uniquely classified as belonging to class k.
The fraction of the kth land-cover type in the footprint (fp) is then defined as Combining Eqs. ( 3) and ( 4), the multi-linear model for the flux becomes A critical assumption under the flux aggregation method is that each land cover k (area A k ) has a constant source strength (F k ).Thus, as in Eq. ( 1), the flux (F ) for a specific area is a weighted aggregation of its various landcover classes.Based on multi-point tower flux measurements (F obs ), multiple-linear regression equations can be formulated by overlaying the flux footprint with a land-cover map (X fp,k ), as follows in Eq. ( 5).In this study, the multiple-linear regression method (using the "regress" and "robustfit" algorithms from the Matlab ® statistical toolbox) determined the regression coefficients (estimates of the specific flux for each land-cover class in the case F k ) by minimizing the squared residuals, and the standard error of each coefficient estimate was also quantified.So the specific flux (F k ) not only contains estimates, but also its uncertainties in regression.For each LAS path, the measured flux (i.e., sensible heat flux) can also be disaggregated into the component flux by the relevant footprint function as in Eq. ( 5).This can be taken as a validation of the former step.The accuracy of this method is highly dependent on four most important aspects: (1) better flux data for all EC sites; (2) a better land-cover classification map; (3) more precise flux footprint analysis; and (4) good flux and footprint data for LAS.So properly processed flux data, an accurate highresolution land-cover map, and appropriate footprint models are the foundation of formulating a better multiple-linear regression.Sometimes the multi-linear regression equations may not be solvable.When this problem happens, the classification accuracy of the land-cover map should be carefully checked, and the selected footprint model should also be verified with an alternative method.

Footprint models
The Eulerian analytical footprint model, which developed by Kormann and Meixner (2001), was used for estimating the single time flux footprint of EC measurements, due to its ease of use and wide range of stability as well as its numerical stability (Leclerc and Foken, 2014).Besides, as we have checked, its footprint estimates were in good agreement with the calculations of more sophisticated backward Lagrangian footprint models, such as the Kljun scheme (Kljun et al., 2002;Kljun et al., 2015).The footprint function w(x, y, z) can be expressed in terms of a crosswind integrated flux footprint function f y (x, z) and a Gaussian crosswind distribution function D y (x, y).The analytic solution of Kormann and Meixner (2001) is depicted by Eq. ( 6).More details on the derivation of f y (x, z) and D y (x, y) as well as the relevant parameters can be seen in Kormann and Meixner (2001).
The flux contribution source area of LAS measurements was estimated by combining the footprint function w(x, y, z) for point flux measurement with the path-weighting function W (x) of LAS (Meijninger et al., 2002).For equal transmitter and receiver apertures, this path-weighting function is symmetrical bell-shaped having a center maximum and tapering to zero at the transmitter and receiver end.For the LAS footprint calculation, the approach of Korman and Meixner (2001)  estimation.The equation of the LAS footprint function is where x 1 , x 2 are the positions of LAS transmitter and receiver, respectively.x, y represent the locations of points along the path of LAS.x , y are the coordinates of each of the upwind points.z LAS is the effective height of LAS measurements.
To obtain the averaged footprint of EC flux measurements (such as daily and monthly), the flux-weighted footprint climatology method was applied for each pixel (Liu et al., 2016).The equation of the weighted footprint climatology is Here i denotes the time step (e.g., 30 min), N is the total number of 30 min periods within the selected time frame (such as daily scale), Flux(i) is the EC observed flux at the i time step (half-hourly LE in our case), and w i (x, y, z) represents a half-hourly footprint estimate calculated via Eq.( 6).
The inputs of the analytical footprint model mainly include the measurement height, wind direction, wind speed, and the Obukhov length, which can be easily derived from flux tower measurements.The daily-averaged flux footprint of the EC observations was calculated by Eq. ( 8).The flux contribution of the total source area was set to 90 % for both EC and LAS measurements.The normalized daily-averaged footprint of ECs and half-hourly footprint estimates of LASs were separately overlaid with the land-cover map to determine the footprint-weighted contribution of each land-cover class to the measured flux from the EC and LAS systems.

Data processing flow for the determination of area-averaged fluxes
The overall data processing flow for determining the areaaveraged evapotranspiration over a heterogeneous land surface mainly includes three steps (Fig. 2).Firstly, the spatial representativeness of 16 EC sites within the 5×5 km 2 experimental area was quantitatively assessed by overlaying flux footprint climatology with a highresolution land-cover map.Detailed analyses on this aspect are going to be presented in the following section.
Secondly, the sensible heat fluxes we aggregated were evaluated via the LAS measured path-weighted areal flux.Specifically speaking, based on footprint analysis and a highresolution land-cover map, the sensible heat flux for individual land-cover types was firstly disaggregated from multiple EC flux measurements by performing a multiple-linear regression analysis (Eq. 5).To obtain area-averaged fluxes representative of the LAS source area, the EC-disaggregated fluxes for all land-cover classes were aggregated again according to the fractional weight of each land-cover class in the LAS footprint (Eq.4).Finally, the EC-aggregated fluxes were compared with LAS observations.At last, the area-averaged evapotranspiration over a heterogeneous land surface was estimated with EC measurements and the flux integration scheme proposed and verified in the last step (Eq.1).

Results and discussion
4.1 The characteristics of the surface heat and water vapor fluxes The sensible heat flux over residential area (EC04) was higher than that over the vegetated surfaces (Fig. 3a) and reached a maximum of about 150 W m −2 at afternoon, while the latent heat flux was smallest among all sites (Fig. 3b).
Over the vegetated surfaces (maize, orchard, vegetable), the sensible heat flux was normally less than 100 W m −2 because of sufficient irrigation (Fig. 3a), and it was significantly different between different vegetated surfaces (Fig. 4a).Besides, difference in latent heat flux was also clear.
The maize fields showed higher latent heat fluxes and lower sensible heat fluxes than that of other two vegetated surfaces.One of the possible reasons is that both of the orchard area and vegetable field are comparatively sparsely vegetated than the maize field during the crop growing period.Over maize field, the sensible heat flux was relatively small and even negative in the midafternoon (known as the "oasis effect") (Figs.3a, 4a); while the latent heat flux was quite large, with a maximum value up to 600 W m −2 (Figs.3b, 4b).There was also a difference in sensible and latent heat fluxes among maize sites (Fig. 4).The values of standard deviation (SD) of LE and H for 13 maize sites were about 43 and 8 W m −2 , respectively.

Analysis of the representativeness of the multi-point EC flux measurements
To further understand the variability of surface energy fluxes between different sites in a heterogeneous landscape, footprint analyses for the spatial representativeness of 16 EC sites were performed by superimposing flux footprint with a landcover map derived from aircraft images (Fig. 1).The fraction of all land-cover classes present in the daytime-averaged footprint is shown in Fig. 5. Due to the differences of observation height and some variations of wind direction/speed, the EC source areas of differ-  ent sites were distinctly different (Fig. 1).For each EC flux measurements, the relative contribution of each land cover to the total measured flux was presented in Fig. 5.For sites 1 and 17, the dominated surface types in the source area were vegetable and orchard, respectively.For site 4, the fractional weight of the non-vegetation type in the footprint was about 0.79, and that for woods was around 0.21.For other sites, the relative contribution of maize field to the EC measured flux was generally higher than 0.9.However, for site 2 and site 9, the non-vegetated surfaces covered 0.15 to 0.20 of the footprint area; while at site 10, the contribution of vegetable type to the flux measurements reached about 0.12.
The above analysis shows that the EC-matrix measurements are generally representative of the dominated surface types in the HiWATER intensive observation area.The "averaged" fluxes determined by weighting the upwind surface flux by individual land-cover classes and the relevant flux footprints will be representative of the target area.

Evaluation of the EC-aggregated fluxes
Firstly, the sensible heat fluxes for each land-cover type were disaggregated from EC flux observations via the scheme mentioned in Sect.3.1.Figure 6 gives the variation of the EC-disaggregated component fluxes and their uncertainties in the regression.The diurnal cycle of the sensible heat fluxes for each land-cover class was highly significant, and the difference in fluxes between different land-cover types was also distinct.The mean standard error (SE) of estimates of sensible heat flux for maize fields was the lowest, with a value of about 4 W m −2 , mainly because of its dominant land-cover type and the majority of EC sites, while the SE values for land-cover types of non-vegetation, vegetable, and woods were larger, with values of about 15, 12, and 13 W m −2 , respectively, mainly due to the heterogeneity of the underlying surfaces, especially for site 4.
Then, the EC-disaggregated fluxes for all land-cover classes were aggregated again to obtain area-averaged fluxes representative of the LAS source area.Figure 7 illustrates a scatterplot of half-hourly area-averaged sensible heat fluxes estimated from the EC flux matrix (hereafter referred to as H_ECagg) versus LAS measurements (H_LAS), as well as the linear regression parameters (equations and R 2 ).The different statistics between them are also listed in Table 3.
For LAS1 (see Fig. 7a and Table 3), a good agreement is found between EC aggregated fluxes and LAS measurements, with high correlation and a low root mean square error (RMSE) value (R 2 = 0.75, RMSE = 1 W m −2 ).The scatter points in the graph are nearly close to the 1 : 1 line.The MBE and MAPE values were almost 4 W m −2 and 10 %, respectively.For LAS2 (Fig. 7b), there was a little more scatter and a little higher RMSE (R 2 = 0.59, RMSE = 7 W m −2 ), probably owing to a slight parts of village buildings under its path center (Fig. 1).For LAS4 (Fig. 7d), the sensible heat flux aggregated from the EC matrix was a little larger than that from LAS measurements (R 2 = 0.65, RMSE = 13 W m −2 ).For LAS3 (Fig. 7c, Table 3), there was a slightly weak relationship between LAS measured sensible heat fluxes and that from multi-point EC measurements.The scatter points were overall below the 1 : 1 line, and both RMSE and MBE were larger than those of LAS1, LAS2, and LAS4.This would be from a larger area of residential buildings distributed under the central part of the LAS path (Fig. 1), and their representative sensible heat flux derived from a single site 4 would be insufficient.Besides, as the closure ratio of the surface energy balance (EBR) is comparatively lower over heterogeneous areas, such as for site 4 (Xu et al., 2017), the energy flux from large eddies or secondary circulations induced from surface heterogeneities may not be captured by a singlepoint EC, but may be measured via the LAS system.Thus, the observations might be able to close the surface energy balance better the EC method (Foken, 2008;et al., 2010).
The uncertainties of LAS, induced in its data processing by use of EC measurements (such as Obukhov length and Bowen ratio), also have considerable influence on the comparison.When using a constant Bowen ratio as in former processing of the LAS data (Xu et al., 2013), the MAPE of the comparison between LAS and EC aggregation would be increased by about 20 %, and the value of the RMSE would be more than 10 W m −2 .
The present flux aggregation method relies on the using of footprint model and high-resolution land-cover map (Eq.5).We have checked the footprint results between the Kormann and Meixner (2001) model we used with that of Kljun's scheme (Kljun et al., 2002(Kljun et al., , 2015)), the differences were apparently minor.Hutjes et al. (2010) also reported that the un-certainties about the real size of the footprint do not affect the disaggregating of flux estimates too much as long as the uncertainty does not lead to significantly different fractional covers.Besides, the error in the land-cover classification map we used is no more than 3 %.So the effect of uncertainties in X fp,k on estimates of F k is rather small.As mentioned in Sect.2.2.1, the systematic error of EC fluxes was also small.Therefore, disaggregation from fluxes over a heterogeneous area into individual land-cover types is feasible.The present scheme is also practical over a truly complex landscape or a simple "binary" land (Hutjes et al., 2010).
In the comparison with LAS measurements, both the arithmetic-averaged and area-averaged methods (Liu et al., 2016) displayed a relatively larger scatter and larger values of RMSE and MAPE than the present aggregation method, even for a relatively homogeneous land surface (e.g., near the path of LAS1).The present scheme not only has a much better database, but also has a solid grounding in physics and mathematics in the aggregation of area-averaged fluxes over a heterogeneous surface.

Estimation of area-averaged evapotranspiration (ET)
The present flux aggregation scheme, as evaluated in Sect.4.3, was adopted to determine the area-averaged ET over the study area with multi-point EC flux measurements and high-resolution land-cover map.The EC disaggregated daily ET for four land-cover classes on 2 typical clear days is shown in Fig. 8; the standard errors of the ET estimates from least squares regression are also plotted as error bars.As can be seen, the daily ET values for maize field were the highest (7-8 mm) during the crop growing season.For woods areas the value of daily ET was about 6.5 mm; and for vegetable field it was 6-7 mm.On the contrary, the daily ET for nonvegetation type varied largely, with values about 2.8 mm on 29 June and 1.6 mm on 30 June, respectively.Table 4 lists the total ET for different land-cover classes and their proportions of the total area.The total ET for the study area was approximately 169.7 × 10 3 m 3 day −1 on 29 June and 153.6 × 10 3 m 3 day −1 on 30 June.The ET of maize fields occupied more than 80 % of the total area, due to its domination in the study area, while that for woods, vegetables, and non-vegetation fields occupied about 8, 5, and 5 %, respectively.
The area-averaged ET over the kernel experimental area of HiWATER was finally estimated, with values of approximately 7 mm day −1 on 29 June and 6 mm day −1 on 30 June, respectively.

Summary and conclusions
On the basis of an accurate high-resolution land-cover classification map and multi-point ground-based flux measurements from 16 EC systems and four groups of LAS systems during the intensive observation period of HiWATER, a flux aggregation method for determining area-averaged flux was established through the combination of footprint analysis and multiple-linear regression.The method was applied to estimate the area-averaged surface fluxes over a heterogeneous surface from multi-point EC flux measurements, and its results were verified by the LAS measurements.Ultimately, the integration method was applied to estimate area-averaged ET over the study area.
Robust quality control and uncertainty estimation for the EC and LAS data, done through careful data processing and inter-comparison as well as assessment of the energy balance closure, ensure the accuracy of the flux dataset used in data Hydrol.Earth Syst.Sci., 21, 4037-4051, 2017 www.hydrol-earth-syst-sci.net/21/4037/2017/ analysis.For the deep interpretation of the surface fluxes over different land surfaces, the combination of footprint analyses for the representativeness of EC flux measurements and high-resolution land-cover map can be a practical way, and it is also the foundation for the establishment of the flux aggregation algorithm.With a high-quality flux dataset (EC & LAS), precise flux footprint estimates, and an accurate land-cover classification map, a flux aggregation method can be successfully established by multiple linear regression analysis.It achieves the goal of determining the area-averaged fluxes over heterogeneous areas from the EC flux matrix.However, the agreement between the results of the flux integration method and LAS measurements partly relates to the heterogeneity of the land surface, as it influences greatly the energy balance closure of the EC flux dataset.On the other hand, the bias may be partly attributed to the insufficient distribution of flux stations on some dominant land covers; the uncertainty of LAS measurements was also a factor that affects the result of the comparison.
In spite of the limitations mentioned above, the current flux integration scheme provides a unique opportunity to discompose the heterogeneous land surface fluxes into their single components, and the disaggregation process has the potential to scale up multiple EC measurements to an oasis landscape, even to a whole river basin.Besides, compared with the formerly used and rather simple approaches (e.g., the arithmetic average and area-weighted methods), the present scheme is not only with a much better database, but also has a solid grounding in physics and mathematics in the integration of area-averaged fluxes over a heterogeneous surface.Results from this study, such as daily ET at the satellite pixel scale, can be applied for the validation of flux estimates of meso-γ -scale (1-20 km) models.Furthermore, this work will be extended to the water balance study of the whole Heihe River basin, which is quite practical for hydrological modeling and basin water resource management.
method combining footprint analysis and multivariate regression It is generally accepted that an area-averaged flux equals the area-weighted sum of the component fluxes from individual F. Xu et al.: Aggregation of multi-point EC flux measurements

FFigure 3 .
Figure 3.Diurnal cycle of the sensible heat fluxes (a) and latent heat fluxes (b) between different sites during 29 and 30 June 2012.

Figure 3
Figure 3 depicts the diurnal cycle of the sensible and latent heat fluxes at different sites during two clear days.It not only reveals the energy exchange of different sites but also the significant differences in the magnitude of the sensible and latent heat fluxes between different surface types during the growing season.The sensible heat flux over residential area (EC04) was higher than that over the vegetated surfaces (Fig.3a) and reached a maximum of about 150 W m −2 at afternoon, while the latent heat flux was smallest among all sites (Fig.3b).Over the vegetated surfaces (maize, orchard, vegetable), the sensible heat flux was normally less than 100 W m −2 because of sufficient irrigation (Fig.3a), and it was significantly different between different vegetated surfaces (Fig.4a).Besides, difference in latent heat flux was also clear.

Figure 4 .
Figure 4. Diurnal cycle of the mean sensible (a) and latent (b) heat fluxes for 13 maize field sites and different types of vegetation; the error bar is the standard deviation.

Figure 5 .
Figure 5.The fractional weight of each land-cover class in the daytime-averaged flux footprint of each EC flux measurement for 29-30 June 2012.

Figure 6 .
Figure 6.The diurnal cycle of the sensible heat flux for each land-cover class for 29-30 June 2012; the error bars are standard errors of sensible heat flux estimated in the regression.

Figure 7 .
Figure 7.The comparison between LAS observed fluxes (x axis) and EC aggregated fluxes (y axis).

FFigure 8 .
Figure 8.The disaggregated daily ET for each land-cover class in the kernel experimental area of HiWATER on 29 and 30 June 2012; the error bar is the standard error of ET estimates in the regression.

Table 1 .
Details of eddy-covariance systems in the HiWATER flux matrix.

Table 2 .
Details of large-aperture scintillometers in the HiWATER flux matrix.

land cover class Meteorology data When F k is LE/ET H of each land cover class When F k is H
was still used for the single-point footprint Hydrol.Earth Syst.Sci., 21, 4037-4051, 2017 www.hydrol-earth-syst-sci.net/21/4037/2017/ Figure 2. Schematic illustration of data processing steps; LC: land-cover class; SA: source area; H: sensible heat flux; LE: latent heat flux; ET: evapotranspiration.

Table 3 .
Different statistics between LAS observed flux and EC aggregated flux at LAS sites.O i /n; P i is the EC aggregated value; O i is the LAS observed value; O is the mean measured value; n is the number of samples.RMSE is root mean square error, MAPE is mean absolute percentage error, and MBE is the mean bias error.

Table 4 .
ET for each land-cover class and their proportion of the kernel experimental area ET.