Scaling Flux Tower Observations of Sensible Heat Flux Using Weighted Area-to-area Regression Kriging

Sensible heat flux (H) plays an important role in characterizations of land surface water and heat balance. There are various types of H measurement methods that depend on observation scale, from local-area-scale eddy covariance (EC) to regional-scale large aperture scintillometer (LAS) and remote sensing (RS) products. However, methods of converting one H scale to another to validate RS products are still open for question. A previous area-to-area regression kriging-based scaling method performed well in converting EC-scale H to LAS-scale H. However, the method does not consider the path-weighting function in the EC-to LAS-scale kriging with the regression residue, which inevitably brought about a bias estimation. In this study, a weighted area-to-area regression kriging (WATA RK) model is proposed to convert EC-scale H to LAS-scale H. It involves path-weighting functions of EC and LAS source areas in both regression and area kriging stages. Results show that WATA RK outperforms traditional methods in most cases, improving estimation accuracy. The method is considered to provide an efficient validation of RS H flux products.


Introduction
Sensible heat flux (H) is a significant component of land surface water and energy balance determinations.It not only affects the regional heat budget, but also has a profound impact on regional water circulation [1,2].There are several methods for measuring surface water and heat fluxes, such as eddy covariance (EC) and Bowen ratio systems.However, they usually measure H at small local-area scales, which can only represent an area of a few square meters to several hundred square meters [3].For large-scale H, remote sensing (RS)-based observations are useful [4,5].H products can be determined from a range of global-coverage RS sensors including Landsat 7 (30 m spatial resolution and 16 days temporal resolution), Moderate Resolution Imaging Spectroradiometer (MODIS, 1 km; daily), and Advanced Very High Resolution Radiometer (AVHRR, 1.1 km; 14 times/day) [6][7][8].RS instruments themselves cannot directly measure H. Therefore, methods for validating RS inversion products are critical for using such products.The RS estimations are generally validated by EC measurements directly or by combing their footprint models, which involves many technical problems [9].
Since the 1990s, large aperture scintillometer (LAS) has been widely applied to measure sensible heat flux around the world, and it has also been used to validate RS estimations owing to the large scale of its measurements [10].So LAS operational principles and data processing steps have been the subject of attention to ensure data quality [11,12].Another method used to acquire area surface flux is the establishment of a multi-site observation system on a heterogeneous surface, which is done by installing flux instruments (e.g., EC or Bowen ratio systems) on typical surfaces of the experimental area and then aggregating the small local-area scale observations to estimate larger area fluxes (e.g., model grid-scale) using various upscaling methods.Experiments that have followed this method include HAPEX-MOBILHY (Hydrologic Atmospheric Pilot Experiment-Modelisation du Bilan Hydrique) [13], FIFE (The First ISLSCP Field Experiment) [5], NOPEX (Northern Hemisphere Climate Processes Land Surface Experiment) [14], EVA-GRIPS (Evaporation at Grid/Pixel Scale) [15], and BOREAS (The Boreal Ecosystem-Atmosphere Study) [16].To validate RS pixel scale H-flux with LAS systems, LAS-scale H-flux must first be validated with EC systems.Thus, in this paper, we use the EC-to LAS-scale H-flux conversion as an example to investigate the scaling method.
There are many scaling methods that can be used to incorporate observed H from EC-to LAS-scale, including area-weighted averaging, footprint-weighted averaging, numerical modeling, and kriging [17][18][19][20].Weight average methods require ECs to be placed in homogeneous regions where H is very similar everywhere in each region.Numerical modeling considers the physical mechanisms of H, often under conditions where it is not easy to determine model parameters.In a previous study [20], we proposed an area-to-area regression kriging method to convert observations from an EC observation network to the LAS scale.The method was found to improve estimations of sensible heat flux.The daily spatial trend of H was extracted by a multiple regression.Then, scaling using area-to-area kriging was applied on the residue to interpolate the LAS-scale residue.A point scale variogram was derived from an area scale variogram using a deconvolution process [21][22][23].However, path-weight functions of EC and LAS source areas were not considered in the area-to-area kriging stage.All discretized points were assumed to have an equal weight, which was inconsistent with the footprint model.
In this paper, we introduce an improved H scaling method: weighted area-to-area regression kriging (WATA RK).An example is presented that upscales the observations from an EC network to a LAS scale; however, the method itself should also be applicable for other scales.The method was constructed with a geostatistical framework.

Study Area and Data Description
The study area is in the middle reaches of the Heihe River Basin in northwestern China.The experiment is an important part of the Heihe Watershed Allied Telemetry Experimental Research (HiWATER) program (2012-2015), which was designed as a comprehensive ecohydrological program to capitalize on diverse interdisciplinary studies using the existing observation infrastructures in the Heihe River Basin [24].The experimental area covers 30 km × 30 km, within which our kernel study area covers a 5.5 km × 5.5 km region.Four groups of LAS (two sets in each group) were installed in 3 × 3 and 2 × 1 MODIS pixels.In the kernel area, 18 EC system sets were installed for the experiment, 17 of which were used in this study (Figure 1) [25].There are also 17 automatic weather station sets (including a super-station) installed in the field [24,25].Nine observations (06/15, 06/24, 07/10, 08/02, 08/11, 08/18, 08/27, 09/03, and 09/12) made between June and September 2012 were selected for this study.On these days, the sky was free of clouds and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) products were available.Observations and footprints of the 17 EC and four LAS groups between 12:00 and 12:30 BST (Beijing Standard Time) were averaged to obtain a mean value when the atmosphere stability was under unstable conditions with high turbulent fluxes.Wind speed data were also collected at the same time.By considering the wind direction, the wind speeds were decomposed into north and east components (WS), and then weighted averages were calculated with the corresponding EC and LAS footprint weights.In addition to these observations, synchronous ASTER RS products at 12:15 BST were downloaded from the US-National Aeronautics and Space Administration and processed according to Zhou et al. [26] and Bastiaanssen et al. [27], including the Land Surface Temperature (LST) [26], Normalized Difference Vegetation Index (NDVI), and Fractional Vegetation Cover (FVC) [27].A detailed description of the study area is provided in Ge et al. [20].

Weighted Area-to-Area Regression Kriging
The complicated non-homogeneous ground surface within the research area can be categorized into six types of land cover.Additionally, there are several roads and shelterbelts in the region, with an irrigation canal (and accompanying road) running from northwest to southeast.To undertake a geostatistical interpolation, assumptions of stationarity in the mean and in the second-order momentum are required [28,29].The stationarity in the mean assumes that the mathematical expectation of samples is constant and independent of location, whereas the stationarity in the second-order momentum assumes that covariance is invariant between any two points separated by the same distance and direction.The surface H flux with complicated ground obviously does not meet the assumptions.Based on the research of Ge et al. [20], a modified two-stage, weighted, area-to-area, regression kriging modeling method was designed to improve the upscaling estimation (Figure 2).The first stage involves spatial trend extraction using a multiple linear regression technique.By removing the spatial trend, the residual is assumed to be similar to a homogeneous surface.Then, in a second stage, a WATA RK method is used to estimate the residual of LAS H. Finally, by combining the regression and WATA results, we can estimate the LAS H.

Spatial Trend Extraction
H is the surface turbulent heat flux depending on the air flow and the surface properties.It is a function of both the air flow and the surface.According to Ge et al. [20], two surface variables (FVC and NDVI) and two ambient meteorological variables (LST, WS) can be specially selected by a stepwise regression method to extract the H spatial trend.Variables from the 17 EC systems including all nine periods are pooled together to obtain a stable regression relationship, where i is the index of the EC systems ( = 1, … ,17), and t is the index of the periods ( = 1, … ,9)., , , and are footprint model weighted average variables of the ith EC's FVC, NDVI, LST, and WS in the tsh time, respectively.and ̅ are the weighted average H and regression residue.β , β , …, β are regression coefficients.The footprint model shows that the H of an EC (or a LAS) has unequal contributions from different parts of its source region.Both the shape of the source region and the weight of the contribution are obviously affected by meteorological conditions.Variables that enter into the regression model are abstracted by the footprint models as follows.
where is the footprint weight of the jth discrete point, and n is the number of discrete points in the ith EC.Thus, the residual ̅ can be expressed as the following equation.
where is the residue of the ith EC's jth discrete point in the tth period.From Equation (3), it can be found that the residue is also a version that is weighted by the footprint model.

Weighted Area-to-Area Kriging
Residues are assumed white noise in a classic linear regression model.However, when they are distributed on a geographical space in some pattern, there is spatial autocorrelation between them.Then, near points have similar residues.It can be used to estimate an unknown location's regression residue.Exploration data analysis of the H regression model showed that there are rather strong spatial autocorrelations in the regression residues.It would be helpful to improve the estimation accuracy of H by adding back the interpolated residue.After removing the spatial trend, the residue can be treated as a random spatial field and assumed to meet the stationarity in the second-order momentum.For a selected LAS, the true H residue (denoted as ) over the entire area at a specific period can be given by: where m is the number of LAS L0 for the discretized points; and are the ith point's footprint weight and residue, respectively.When nearby EC systems are available and there is a strong spatial autocorrelation in the random field, the can be estimated by nearby EC systems with proper weights: where ̂ is an estimation of , N is the number of ECs, and λ and are the weight and residue of the ith EC, respectively.
A good geostatistical estimation should be unbiased and with a minimum estimation variance, i.e., where E(•) is the mathematical expectation.Assuming a first-order stationary state, i.e., E[ ] = E = , from E( ̂ ) = we can get ∑ λ = 1.This means that the weights of all ECs should be summed to one.The left part of Equation ( 6) can then be written as where C(•,•) is the covariance operator.All three items in the equation are area-to-area covariances.They cannot be calculated directly, but should be viable point scale covariances.The three items can then be further expanded separately as follows.
where C , is the covariance between the ith EC's kth discretized point and jth EC's lth discretized point at some period, while and are footprint weights of the two points, respectively.
is the footprint weight of the jth discretized point of target LAS L0.C , is the covariance between the ith EC's kth discretized point and the jth discretized point of L0.C , is the covariance between the ith and jth points of L0.In these equations, λ ( = 1, … , ) are unknowns to be solved.Point-scale covariance can be calculated only after a point-scale variogram is deconvolved from the area-scale variogram, as shown in Section 2.2.3.Similar to ordinary kriging, Equation ( 6) is a typical linear optimization problem that can be solved by the Lagrange multiplier method.Then, we can define the following linear equation, where μ is the introduced Lagrange multiplier.To minimize Equation ( 9), we compute the partial derivatives of with respect to the weights λ ( = 1, … , ) and the Lagrange multiplier μ, and equate each of them to 0, i.e., = 0 and = 0.By rearranging the derived N+1 terms, we obtained the following linear WATA kriging system.
where , are discretized points corresponding to the ith and jth EC, respectively.This is a linear system with N+1 unknowns and N+1 equations.The unknowns can be easily solved by a basic matrix inversion or Gaussian elimination method.Then, the LAS L0 residue can be estimated from Equation (5), and the estimation variance can be calculated from Equations ( 7) and (8).
There is a simplified form of estimation variance that can avoid the × terms in the quadruple summation.By multiplying each of the first equations in Equation ( 10) by λ , then summing these equations together and rearranging, we obtain the following result [29]: Substituting this into Equations ( 7) and ( 8) allows us to express the minimized variance as Note that the estimation variance is mainly determined by point-scale variance within LAS L0 and point scale variance between L0 and nearby ECs.

Point-to-Point Variogram Reconstruction
Because the source regions of all 17 EC systems are at the area scale, a point-scale variogram of H residue cannot be directly calculated from them.According to ATA kriging, we have derived point-scale covariance values for the H residue by a deconvolution method (Figure 3).First, an empirical area scale variogram is calculated by pooling the residue data from all nine periods.At this step, each EC is treated as a single point by its centroid.A theoretical area variogram is fitted with a Matérn model, which is the "true" area variogram model referred to by later regularized variograms.Then, a point-scale variogram model is initialized with empirical parameters.Based on this point variogram model, a regularized area variogram can be obtained as follows: The regularized area variogram can be compared with the theoretical area variogram.The point variogram model is then modified to decrease the difference.The process is iterated to minimize the difference between a regularized area variogram model and the corresponding theoretical area variogram model.The model parameters can also be solved by maximum likelihood estimation [30].

Results and Discussion
After abstracting FVC, NDVI, LST, and WS with the footprint models, nine regression models were built for each of the periods based on the proposed multiple linear regression model.The four independent variables were also abstracted for the four LAS footprints.By comparing regression results with the observed LAS H, spatial trends in LAS H could be captured well, especially for LAS1 (Figure 3a).The black line in Figure 3 is the 1:1 line.The red dashed line is the regression line between the observations and estimations, with no intercept, which means that it is the regression line in force that passes the origin (0,0).H is over-estimated for LAS2 and LAS4, whereas it is under-estimated for LAS3.The residual part of LAS H is then estimated from the 17 EC H residuals by the above WATA kriging model.The final estimated result is obtained by adding the regression and interpolation parts.A comparison with observed values is shown in Figure 3b.There are obvious improvements compared with the regression.For both LAS1 and LAS2, the red dashed lines are almost coincident with the black lines.The slope value of LAS4 is 1.056 (i.e., very close to 1).From Figure 3a, it can be found that the overestimation is mainly caused by the regression part in the two-step modeling.After applying the WATA kriging to the regression residue and adding the result to the regression result, the overestimation is decreased at both LAS2 and LAS4.However, the under-estimation phenomenon still exists in LAS3, even though the slope is increased from 0.776 to 0.805.The most important reason for this may be a systematic bias between the LAS and EC values owing to the large number of buildings in the LAS source area [20].The high H resulting from buildings can be captured by the LAS, but it is not captured by EC because almost all EC systems are located within corn fields.
Area-and footprint-weighted models are two other H upscaling methods often used for aggregating EC-to LAS-scale H [17,18].Our H estimations have also been compared with the results of the two models undertaken by Liang [31].Unlike a WATA regression kriging model, for which the EC weight is determined by solving a linear system to guarantee a best unbiased estimation, for area-and footprint-weighted models, the research area has to be partitioned into small homogeneous strata where the H are considered to be the same everywhere.Additionally, there is at least one EC system in each stratum.The LAS source areas are partitioned into strata in the same way.Then, the LAS H is the weighted sum of H for all strata.For area-weighted models, the weights are the area ratios of each stratum in the LAS source area.For footprint-weighted models, the weights are the footprints of each stratum.Sometimes, it can be difficult for a complex research area to be partitioned into homogeneous H strata.However, the proposed method does not require such strata, which makes it more appropriate for use in complex areas.Liang [31] compared area-and footprint-weighted models to ATA RK using the same input data (Table 1).The slopes between observations and estimations are all less than 1, which suggests that they are more likely to underestimate the LAS H in this research area.Multiple linear regression and ATA RK perform better than the other two methods.However, in general, WATA RK outperforms all of them, except in the LAS4 because its slope is slightly greater than that of ATA RK.The root-mean-square-error (RMSE) and mean-bias-error (MBE) are also calculated and compared in Table 1.Generally, area-weighted and footprint-weighted models have larger RMSE than the other three models.When four LAS's results are considered together, there is no obvious different RMSE for the multiple linear regression, ATA RK, and WATA RK models.If two models of the three are selected, each model has two RMSE values smaller than the other one, while at the same time it also has two RMSE values larger than the other one.However, WATA RK has a slightly smaller MBE than multiple linear regression and ATA RK.
EC measurement is thought to be accurate for heat flux.It is usually used to validate larger scale measurements and products.Zeweldi et al. [32] installed two EC systems and a set of LAS systems in a homogeneous dry semi-arid region.The linear correlation coefficients between the measured sensible heat flux from LAS and the measured values from the two ECs were 0.79 and 0.89, respectively.Xu et al. [25] compared both LAS and EC measurements.Measurements from various ECs agreed very well, along with measurements from various LASs.On a relatively homogenous surface, they found that the regression slope between LAS and EC had average measurements of 0.92 with R 2 of 0.89.The discrepancy is considered to be caused mainly by the heterogeneity of the surface.An improved aggregation method might enable higher correlation and R 2 coefficients.The discrepancies in measurements of LAS and EC are more obvious for heterogeneous surfaces.This might be the result of an energy imbalance phenomenon in response to the different source areas of the LAS and EC measurements [3].The energy imbalance of EC was found in almost all the experiments around the world, and the reasons leading to the imbalance were still under-debated [33].Previous studies usually considered H from EC was correct in the energy budget when compared with models [34].The energy balance closures at ASTER passing time were around 0.8 in our study area.The footprint model can quantitatively define the contribution of different source areas.The usefulness of the models might not be obvious when surfaces are homogeneous.However, this research shows that on a heterogeneous surface, because different types of sources have different water and heat conditions, weighted averages become more important in estimating upscaled heat flux.Compared with the above upscaling methods, the proposed regression kriging considered three important aspects of important information at the same time: (1) footprint weights of EC and LAS are considered in the aggregation of both H-fluxes and related independent variables; (2) a multiple linear model is used to capture daily large-scale spatial trends and patterns; (3) spatial autocorrelation is found between the local scale residues, and it is used to interpolate the corresponding LAS-scale residue.

Conclusions
Sensible heat flux has attracted great attention in RS because of its importance in quantifications of land surface water and heat balance.Different scales of observation are possible when measuring H. EC and LAS are two of the most popular measurements.RS provides numerous products [4,5,35]; however, the validation of large-scale observations and products is a necessary first step in their application to different domains.Although many scaling methods are available for aggregating EC-scale H to LAS-scale H (e.g., area-weighted, footprint-weighted, and ATA kriging), we propose a WATA RK model for H aggregation, which was found to improve the accuracy of H estimations.The object of the model is to upscale EC-scale H to LAS-scale H but not for mapping the spatial distribution of H flux.The method includes three main steps: first, the spatial trend of H should be extracted using a weighted regression model; then, the residual of LAS H can be estimated with a WATA kriging; finally, the LAS H can be estimated by combining the two parts.We suggest that this method is applicable for other scales of H aggregations and validations of RS products.In particular, it is suitable for complicated and irregular surfaces where the component source areas have different water and heat conditions.We will further validate and test different scale RS products based on the proposed method.

Figure 2 .
Figure 2. Process flowchart of WATA RK for sensible heat fluxes.

Figure 3 .
Figure 3.Comparison of estimations and observations for LAS H. (a) Estimated by regression; (b) estimated by WATA regression kriging.

Table 1 .
Estimated results comparison between different models.