Contributions of landscape heterogeneity within the footprint of eddy-covariance towers to flux measurements

Agricultural and Forest Meteorology 260–261 (2018) 144–153

Contents lists available at ScienceDirect

Agricultural and Forest Meteorology journal homepage: www.elsevier.com/locate/agrformet

Contributions of landscape heterogeneity within the footprint of eddycovariance towers to flux measurements

T

⁎

Vincenzo Giannicoa,b, Jiquan Chenb, , Changliang Shaoc, Zutao Ouyangb, Ranjeet Johnb, Raffaele Lafortezzad,a Department of Agricultural and Environmental Sciences, University of Bari “A. Moro”, Via Amendola 165/A 70126, Bari, Italy Center for Global Change and Earth Observations (CGCEO), Michigan State University, East Lansing, MI 48823, USA c National Hulunber Grassland Ecosystem Observation and Research Station & Institute of Agricultural Resources and Regional Planning, Chinese Academy of Agricultural Sciences, Beijing 100081, China d Department of Geography University of Hong Kong Centennial Campus, Pokfulam Road Hong Kong a

b

A R T I C LE I N FO

A B S T R A C T

Keywords: Eddy covariance Landscape heterogeneity Footprint Bayesian Mongolia

Flux measurements based on aerodynamic principles (e.g., eddy covariance method, or EC) assume the vegetation and landform within the footprint must be flat, large, and homogeneous, although very rarely do flux tower sites meet such requirements. Here, using two long term EC tower sites in Mongolian grasslands, we test a hypothesis that the magnitude and variation of EC flux measurements are partially dependent of landscape heterogeneity. We define landscape heterogeneity as the spatial composition and distribution of components within the footprint of a flux tower, which was quantified using high resolution WorldView-2 images to extract vegetation texture features (Contrast, Dissimilarity and Entropy). Bayesian analysis was performed to model the linkage of landscape heterogeneity with EC fluxes of CO2 in 8 intercardinal directions by dissecting it into 5 distances. We found that higher levels of landscape heterogeneity have an impact on flux measurements, especially under stable conditions. Specifically, a total of 24 Bayesian models based on the EVI-derived texture features passed the Gelman and Rubin convergence diagnostic statistic test. A positive relationship is shown between the percentage of footprint cover (Footprint%sector) and landscape texture features (Contrast and Dissimilarity), with the footprint cover acting as a function of heterogeneity under stable conditions. Negative effects were found when modeling CO2 flux (Fcsector) with Contrast under stable and unstable conditions, and with Dissimilarity and Entropy under stable conditions. With increases in high resolution remote sensing images and UAV technology (e.g., LiDAR), the results and approaches outlined in this study highlight new frontiers and opportunities for the FLUXNET community to integrate flux measurements and high-resolution remote sensing data, promoting a new generation of footprint models, and exploring the cohesive connections between flux measurements and the underlying processes (e.g., soil, physiological, ecosystem processes).

1. Introduction A fundamental assumption in flux measurements based on aerodynamic principles, including the eddy-covariance (EC) method, Bowen ratio method, and others, is that the vegetation and land form surrounding a tower is flat, large, and homogenous (Lee et al., 2006). This assumption is necessary for ensuring that the contributions from horizontal and vertical advections are minimized for accurate calculations of the vertical fluxes (e.g., the net exchange between the atmosphere and the vegetation). Otherwise, the non-zero mean advection induced by convergence or divergence of flow due to spatial source/sink in homogeneity needs to be corrected (Lee, 1998; Paw et al., 2000;

⁎

Corresponding author. E-mail address: [email protected] (J. Chen).

https://doi.org/10.1016/j.agrformet.2018.06.004 Received 29 January 2018; Received in revised form 21 May 2018; Accepted 5 June 2018 0168-1923/ © 2018 Elsevier B.V. All rights reserved.

Massman and Lee, 2002). Given these requirements, all parcels within the footprint of a flux tower will have the potential to contribute equally to the quantities measured at the sensor location so that the ecosystem-level fluxes are accurately represented. This assumption is conventionally met through site selection when installing a tower. Unfortunately, such an ideal site rarely exists across global terrestrial surfaces. We speculate that many tower sites within the FLUXNET community may not meet these requirements. Few scientific reports in the literature have attempted to assess the heterogeneity of flux towers using high resolution satellite images, nor has an effort been made to quantify the contributions of vegetation mosaics within the footprint to the vertical fluxes of materials (Xu et al., 2017).

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

South South-West (SSW), South-West West (SWW), West North-West (WNW), North-West North (NWN) (Fig. 2). Each section was successively divided into 5 sectors by distance from the tower (i.e., 50 m, 100 m, 200 m, 500 m, 1 km), yielding a total of 40 sectors around each EC tower.

Numerous efforts, however, have been made to correct the in situ flux measurements to overcome the conversion/diversion processes caused by complex terrains (Froelich and Schmid, 2002; Lee and Hu, 2002; Davis et al., 2003). The dominant approach has been based on rotating the mainstream fluxes to reflect the vertical fluxes, including the widely recognized methods of 2D, 3D, or planetary rotation (Goulden et al., 1996; Twine et al., 2000). This is achieved by advancing the classical footprint model (Schmid, 2002) to examine the deviations from the horizontal plain (Kormann and Meixner, 2001; Kljun et al., 2002; Froelich and Schmid, 2002). These corrections have improved our estimates taken on slopes and in areas with directional air movement (e.g., cold drainage). However, no study has considered the error terms caused by vegetation heterogeneity (i.e., spatial variation in vegetation composition and structure), likely due to the lack of highresolution data on the vegetation and landform within the footprint of a tower. With the availability of high-resolution data from recent satellites or UAV, it is now possible to integrate the footprint models and flux measurements with heterogeneity metrics. More importantly, connecting the vegetation heterogeneity with tower-based fluxes will provide us with a more direct and powerful avenue to explore the role of vegetation in regulating ecosystem fluxes. In this study, we define landscape heterogeneity as the spatial composition and distribution of landscape components within the footprint landscape (Forman, 1995). Such heterogeneity is often explained in terms of landscape features calculated from images using Haralick texture features (Haralick et al., 1973; Kayitakire et al., 2006; Ozdemir and Karnieli, 2011). For example, it is possible to calculate texture features related to the differences among neighborhood pixels (i.e., Contrast and Dissimilarity) and to the level of orderliness among pixel values (Entropy). The primary objective of this study is to test a hypothesis that the magnitude and variation of EC flux measurements are partially dependent of landscape heterogeneity as defined above. In particular, we aim to answer the following questions: (i) Is the variation in EC fluxes (e.g., CO2) partially dependent on landscape heterogeneity? (ii) If that dependency exists, are there distinctions among different texture features? And (iii) What are the differences between estimates under stable and unstable conditions? To answer these questions, we investigated the empirical relationships between texture features derived from highresolution data at two EC flux towers in Mongolia where the surrounding landforms were flat and the vegetation was seemingly homogeneous for at least 1 km (Shao et al., 2017). CO2 flux (Fc) was partitioned over space using a footprint model and normalized over wind direction by employing a clustering algorithm. The results and models developed through this study can be used to support decisions in evaluating new and established EC flux tower sites and understanding how within-ecosystem variation of vegetation/soil may influence ecosystem fluxes. Additionally, this study sets a pioneer step in design future footprint models that integrate landscape heterogeneity and EC fluxes for the FLUXNET community.

2.2. High-resolution vegetation distribution We searched the available, high resolution imagery for our study sites during 2014–2015 and found one WorldView-2 image on June 6, 2013 as the closest one by date for the study period to describe the vegetation heterogeneity within the footprint of the tower. This image consists of eight color bands (red, green, blue, coastal, yellow, red edge, and two near-infrared bands) at 1.84 m spatial resolution and one panchromatic band at 0.46 m spatial resolution. The image was provided by the National Geospatial-Intelligence Agency (NGA), commercial Archive Data (https://cad4nasa.gsfc.nasa.gov/). After preprocessing the image, we extracted various texture features (Contrast, Dissimilarity and Entropy). In brief, the ENVI 5.1 tool (Exelis Visual Information Solutions, Boulder, Colorado) was employed to radiometrically and atmospherically correct the image after the orthorectification procedure. We extracted two vegetation indices (VIs) from the image: the enhanced vegetation index (EVI) and the normalized difference vegetation index (NDVI). Both VIs were used to calculate texture features using the Gray-Level Co-occurrence Matrices (GLCMs)based filter implemented in ENVI. GLCM is a matrix that summarizes the relative frequency of how often two different quantized pixel values occur in a specified spatial relationship (Haralick et al., 1973). The GLCM is calculated after quantizing the VIs into 64-integer values with a window of 3 × 3 pixels. The window size was set on the basis of a previously conducted correlation test. The GLCM was therefore used to calculate the different texture features. Texture features can be broadly divided into edge and interior features (Hall-Beyer, 2017). Edge features, including Contrast, Dissimilarity and Entropy, are likely to have greater values in areas with large differences between neighborhood pixels. Areas with these features that yield high values typically represent heterogeneous areas and natural or artificial edges. Contrast and Dissimilarity have similar behaviors and are more likely to increase where linear edges are evident although Dissimilarity increases linearly while Contrast increases exponentially. On the other hand, Entropy is expected to be higher where the pattern of heterogeneity is irregular and the edges are not continuous with such complex system edges (i.e., forest-wetland boundaries). On the contrary, interior features (e.g., Homogeneity, Correlation or Mean) have a positive relationship with homogeneity and can yield high values in areas where the pixels have very similar values to their neighboring pixels. In this study, we focused on texture features describing landscape heterogeneity (Contrast, Dissimilarity and Entropy) that are calculated as the following: Ng

Contrast =

2. Materials and methods

Ng

∑ ∑ (i−j)2 P (i, j) (1)

i=1 j=1

2.1. Study area

Ng

Dissimilarity = The two flux towers considered in our study are located in the west of Ulaanbaatar, the capital city of Mongolia (Fig. 1). The region is characterized by a continental climate; the annual mean air temperature is 1.2 °C; the mean precipitation is 196 mm that follows an irregular seasonal pattern. The first flux tower is dominated by Leymus chinensis and is classified as a meadow steppe (MDW) with permafrost substratum. The second site is classified as a short-grass typical steppe (TPL) and is dominated by Stipa krylovii and Artemisia frigida (Shao et al., 2017). The area around each flux tower was divided into eight secondary-intercardinal directions, labeled as: North North-East (NNE), North-East East (NEE), East South-East (ESE), South-East South (SES),

Ng

∑ ∑ |i−j| P (i, j) i=1 j=1

Ng

Entropy =

(2)

Ng

∑∑ i=1 j=1

P (i, j ) log(P (i, j ))) (3)

where P(i,j) represents the normalized value of the GLCM at (i,j) and Ng represents the number of quantified values of the image. All texture measures were averaged by sector (Fig. 3). Frequently, texture features are correlated with each other. Such correlation may depend on various factors such as the image’s spatial resolution or the ecosystem analyzed. Therefore we selected only those features that were not significantly 145

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

Fig. 1. (a) Location of two flux towers used for thisstudy; (b) enhanced vegetation index (EVI) and a field image of the typical steppe (TPL) area; (c) EVI and a field image of the meadow steppe (MDW) area. Fig. 2. (a) Half-hour footprint 3D plot and (b) half-hour footprint 2D plot. The shape file defining the sectors covers the 2D plot. The area around each flux tower is divided into 8 secondary -intercardinal directions – North NorthEast (NNE), North-East East (NEE), East SouthEast (ESE), South-East South (SES), South South-West (SSW), South-West West (SWW), West North-West (WNW), North-West North (NWN) – and 5 radial distances from the tower (i.e., 50 m, 100 m, 200 m, 500 m, and 1 km), resulting in a total of 40 sectors (8 intercardinal directions × 5 radial distances).

Fig. 3. Texture feature extraction framework. (a) Vegetation indices calculated over the WorldView2 image; (b_1) and (b_2) are two examples of texture features (Entropy and Homogeneity); (c) division of the area surrounding the flux tower into 8 intercardinal directions and 5 radial distances from the tower. 146

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

Fcsector = Fc Footprint%sector

correlated after Pearson’s correlation analysis among the variables.

(4)

We analyzed the relationship between Footprint%sector and the texture measures for each sector and the CO2 fluxes with the same texture metrics. Additionally, for each flux measurement we examined the relationships with different types of empirical models by daytime and nighttime.

2.3. Flux measurements and the footprint model An open-path EC system was installed at 2.0 m above the ground on the towers at the MDW and TPL. The system is comprised of an infrared gas analyzer (IRGA, LI-7500, LI-COR, Lincoln, NE) and a CSAT3 threedimensional (3D) sonic anemometer (Campbell Scientific Inc., CSI, Logan, UT) to calculate net ecosystem exchange of CO2, sensible (H) and latent (LE) with the atmosphere. The 3D wind velocities, sonic temperature, and CO2 and H2O concentrations were sampled at a 10-Hz frequency, and the raw time series (TS) data were recorded. The IRGA was calibrated before field setup and at the beginning of the growing season of each year (see more details in Shao et al., 2017). For the present study, we used the EC data during June-August in 2014 and 2015 when the spatial variation of the vegetation is best represented by the Worldview-2 images. Although there exists a time mismatch between EC data and satellite imagery, our sites were not disturbed by any major natural or human activities. Therefore, we assumed that the landscape structure in 2013 remained the same during 2014–2015. Half-hourly NEE fluxes were calculated using the EdiRe workflow, following the protocols of Chu et al. (2014). In brief, the raw TS data quality was checked, with spikes removed. Time lags between measured scalars and vertical velocity were also removed. The planar fit method was applied to rotate the three velocity components into the mean streamline coordinate system (Wilczak et al., 2001). The raw sonic temperature was corrected, a 30-min blocking average without detrending was used (Moncrieff et al., 2004), and a Webb-PearmanLeuning (WPL) correction was applied to correct for air density influences (Webb et al., 1980). The periods with poor turbulent development were calculated to filter the stationarity, integral turbulence characteristics, and friction velocity (u*) of each 30-min flux. All the successive processing operations were conducted using R statistical programming language (Version 3.3.2) (R Core Team, 2016). To reduce the effect of covariance between landscape flux and wind direction, we selected only those observations that were distributed equally across all wind directions. The observations were clustered through an unsupervised method using the variables that contribute to the footprint model calculation as inputs (i.e., u*, Obukhov length, standard deviation of cross wind speed and horizontal wind speed). Hence, we used a clustering algorithm based on finite normal mixture modeling, implemented in the R package “mclust” (Fraley et al., 2012). The procedure selects the optimal Gaussian finite mixture model fitted by the expectation maximization (EM) algorithm according to the Bayesian information criterion (BIC) and assigns the same class to observations that are similar according to a given set of input variables. After the data were clustered, we selected those observations of the same class that were equally distributed across all wind directions. The footprint for each half-hourly flux was calculated using the Kormann and Meixner (2001) model by rotating the correspondence to the wind direction (Fig. 2). The field-measured grass height was similar for both of the study areas. Therefore, the same measurement height, calculated as the geometric height above the ground minus the displacement height, was used as the input parameter for the footprint calculation. The roughness length was calculated from measured wind speed and friction velocity assuming a logarithmic wind profile. For each of the 40 sectors, the footprint percentage was calculated using a graphical integration approach. Specifically, we used the tool developed by Agroscope ART (Swiss Federal Research Station, Air Pollution/ Climate Group, CH-8046 Zürich, Switzerland) (Neftel et al., 2008) that integrates the footprint density function over defined surface areas. To obtain the footprint of a flux measurement (i.e., CO2) for each sector, we multiplied the flux measurement by the corresponding percentage of the footprint (Footprint%sector) in each of the half-hourly measurements (Chen et al., 2008, 2012). The variable Fcsector was then calculated as:

2.4. Modeling the effects of landscape heterogeneity on fluxes Bayesian analysis was performed to model the linkage of landscape heterogeneity with fluxes. To accomplish this, we developed various Bayesian models that predicted the changes in Fcsector, and Footprint %sector using texture metrics (e.g., Entropy, Homogeneity). The advantage of using Bayesian analysis is to utilize the probability distribution of a model parameter rather than a point estimation of the mean. All the processes were run in the R environment (R Core Team, 2016). The models were constructed using JAGS version 4 (Plummer, 2016), with the CODA packages (Plummer et al., 2006) used for convergence diagnostics, jagsUI (Kellner, 2016) for posterior predictive checks (PPCs), and functions of Kruschke (2014) to create plots that compare actual and simulated data from Bayesian model in order to visualize the predictive capability of the model. We used a t-distribution to model the distributions of fluxes. The tdistribution was preferred over a normal distribution because it: 1) has heavier tails and allows us to accurately describe the data with outliers; 2) is conventionally used as a sampling distribution for hypothesis testing and has also been widely used to model data (Lange et al., 1989; Tsionas, 2002; Jones and Faddy, 2003; Kruschke, 2014). In JAGS, the tdistribution is defined by three parameters: (1) a location parameter “μ” that is equivalent to the normal distribution mean; (2) a precision parameter that can be expressed as a function of the standard deviation (1/σ2), which in turn is comparable with the normal distribution standard deviation; and (3) a shape parameter “ν” that controls the heaviness of the distribution tails. The distribution is heavily tailed when ν is low and becomes normal as ν approaches infinity. The shape parameter is commonly known as ‘degrees of freedom’ in classic null hypothesis tests, but in this context it can be expressed as a normality parameter (Kruschke, 2014). The models are described as:

1 y ∼ t ⎛μ, 2 , ν ⎞ σ ⎝ ⎠

(5)

μ = β0 + β1 x

(6)

β0 ∼ N(0,10)

(7)

β1 ∼ N(0,10)

(8)

σ ∼ unif(0.001, 1000)

ν=a+1 a ∼ exp(

(9) (10)

1 ) 29

(11)

where t, N, unif and exp is the student t alue, normal, uniform and exponential distributions, respectively, y is the flux measurement, and x is a quantitative variable of vegetation texture. All dependent and independent variables were standardized before running the model (i.e., mean 0, and standard deviation 1) and examined for non-informative prior distributions. The prior distribution for the normality parameter v equals (a + 1), where a has an exponential distribution with a parameter equaling to 1/29. This number is chosen due to the fact that the t-distribution becomes normal when the parameter approaches 30 and, conversely, cannot be considered normal if the parameter is < 30. The exponential distribution with a parameter of 1/29 indicates that the mean of the distribution is 29 and therefore the prior distribution has v of (a + 1) and falls between 1 (i.e., extremely different from a normal distribution) and 30 (i.e., resembling 147

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

only EVI models are presented in this paper. The posterior predictive checks (PPC) for the Posterior probability density functions (PPD) yielded an acceptable result (∼0.5) for the majority of the MDW models (Table 3). The lower and higher PPCs were found for Fcsector flux measurement with Entropy (i.e., 0.42 for daytime observations) as the predictor. In contrast, the PPCs for the TPL models were much higher than 0.5 (Fig. 4). Our final models included nighttime data that had PPCs of ∼ 0.5, with other models posting PPC value of > 0.7 or < 0.15 (Table 4). The PPDs of the normality parameter (ν) indicated a median close to 30 for the majority of the MDW models. Instead, for the TPL models, the ν value was low, indicating that the distribution was far from normal. The PPDs of the slope parameter (β1) for the MDW models appeared significantly different by direction and distance for nighttime and daytime models (Fig. 7). All models for nighttime showed high relationships although only in three models were the HDIlow and HDIhigh (i.e., highest density interval for 95% probability) not overlapping zero (Table 3). In particular, greater positive effects were found for the flux measurement models with Contrast to predict Footprint%sector, whereas more significant negative effects (HDIlow and HDIhigh not overlapping zero) were found in the models that predicted Fcsector with Contrast, Dissimilarity and Entropy as predictors. For daytime models, there seems a united conclusion: greater negative effects were found in the models that predicted Footprint%sector. Interestingly, only one TPL model showed quartile intervals not overlapping with 0 for β1 of PPDs with Dissimilarity as the predictive variable (Fig. 8).

Table 1 Texture features selected after correlation analysis for the meadow steppe (MDW) and typical steppe (TPL) towers. Each row represents the mean (standard deviation) of one of the texture feature (Contrast, Dissimilarity and Entropy) calculated over a circular area (1 km radius) centered on the correspondent EC flux tower. Texture

MDW

TPL

Contrast Dissimilarity Entropy

2.45 (1.59) 1.00 (0.56) 58.57 (4.14)

1.95 (0.25) 2.04 (0.04) 34.49 (2.23)

a normal distribution). Posterior probability density (PPD) functions were sampled from four Markov Chain Monte Carlo (MCMC) chains with 12,500 iterations using 1000 burn-in and 500 adaption iterations. Convergence was assessed by applying the Gelman and Rubin diagnostic statistic (Gelman and Rubin, 1992), while the model PPDs were evaluated using posterior predictive checks (PPC). We then used the PPD functions to estimate 1st and 3rd quantiles and posterior highest density interval for 95% probability (HDIlow and HDIhigh). 3. Results 3.1. Texture features and flux measurement selection All the texture features were selected after correlation analysis (i.e., Contrast, Dissimilarity, and Entropy) (Table 1), with their Pearson coefficients of < 0.80. The optimal number of components (i.e., clusters) that resulted from the model-based clustering algorithm was 9 for meadow steppe (MDW) and 13 for short-grass typical steppe (TPL), while the BICs were −2324.3 and −1174.1, respectively. After visual examinations of clusters, we selected measurements that distributed equally by wind direction (Table 2). These included a mean ( ± standard deviation) Fc of −3.37 ( ± 4.84) μmol m−2 s−1 for MDW and −2.89 ( ± 3.18) μmol m−2 s−1 for TPL. The area surrounding the MDW site presented a much higher variation in texture then those at the TPL site.

4. Discussions High-resolution remote sensing satellite imagery is a cutting-edge resource that sheds light on the contributions of landscape heterogeneity within the footprint of EC flux tower measurements. By means of a modeling analysis using high-resolution satellite images, we tested the hypothesis that EC flux measurements are partially dependent on landscape heterogeneity. The results obtained through our acquired data are in accordance with our hypotheses and provide useful information about the potentials of high-resolution remote sensing data on landscape heterogeneity estimation and subsequently on EC fluxes measurements. A stronger relationship with texture features was observed for the MDW flux tower compared to the TPL tower. Although the two flux towers are only approximately 15 km apart, they differ in certain respects. The MDW canopy structure is more complex and commonly provides a greater percentage of cover, and features a higher canopy and therefore a higher quantity of biomass. The melting of the permafrost substratum that supplies additional amounts of water induces such characteristics, whereas the TPL canopy structure is purely rainfed. Furthermore, both areas surrounding the MDW and TPL towers are grazed. However, TPL is much more affected by grazing, which results in a lower percentage of canopy cover and height and possibly homogenizes the vegetation across the landscape. This is a consequence of permafrost and accessibility, as the MDW tower is located in a slight depression while the area surrounding TPL could be considered flat (John et al., 2015). Texture feature extraction through satellite imaging confirmed such differences (Fig. 1). Overall, higher values of mean and standard deviation were found in the MDW study area for each of the texture features extracted, excluding only the Dissimilarity feature (Table 1). The latter was, on average, higher in the TPL study area, although there was still a sensibly higher standard deviation in the MDW area. Such correlations are comparable with those found in previous studies. For example, Sibanda et al. (2017) used texture features extracted from WorldView-2 to estimate the biomass of grasslands in southern Africa with an RMSE of 0.55 kg m−2, while Kayitakire et al. (2006) found a correlation coefficient of determination of 0.76 using the Correlation texture feature to estimate canopy height of a forested area located in eastern Belgium. In our study, the selected texture

3.2. Landscape heterogeneity effects on flux measurements A comprehensive representation of the relationship between the variables used in this study (i.e., landscape heterogeneity metrics and EC flux measurements) for the MDW and TPL towers are shown in Figs. 5 and 6, respectively. Bayesian models developed in this study performed well in describing the relationships between texture and fluxes (Figs. 7 and 8 and Tables 3 and 4). A total of 24 Bayesian models (2 towers * 3 texture features * 2 flux-derived measurements: Footprint%secto and FCsector * 2 day/night measurements) based on the EVI-derived texture features passed the Gelman and Rubin convergence diagnostic statistic test. Despite the unpronounced differences between the models from EVI and NDVI, the EVI textures performed slightly better overall. Therefore, Table 2 The mean (std) of the flux measurements selected observations through the clustering algorithm and number of selected observations in each one of the 8 sectors: North North-East (NNE), North-East East (NEE), East South-East (ESE), South-East South (SES), South South-West (SSW), South-West West (SWW), West North-West (WNW), North-West North (NWN). Flux tower

MDW TPL

Fc (μmol m−2 s−1)

Number of selected observations NNE

NEE

ESE

SES

SSW

SWW

WNW

NWN

−3.37 (4.84) −2.89 (3.18)

66 54

40 89

64 85

14 13

39 49

114 31

81 24

119 247

148

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

Table 3 Results of the Bayesian modeling for the tower MDW. For each model the Posterior predictive check (PPC), 1st, 3rd quartiles, median, low and posterior highest density interval for 95% probability (HDIlow and HDIhigh) for the parameter β1 and median value for the parameter ν, are presented. y (Flux measurement)

Footprint%sector Footprint%sector Footprint%sector Footprint%sector Footprint%sector Footprint%sector Fcsector Fcsector Fcsector Fcsector Fcsector Fcsector

x (Texture feature)

Contrast Contrast Dissimilarity Dissimilarity Entropy Entropy Contrast Contrast Dissimilarity Dissimilarity Entropy Entropy

Time

Day Night Day Night Day Night Day Night Day Night Day Night

β1

ν

Median

1st quartile

3rd quartile

HDIlow

HDIhigh

Median

0.062 0.415 −0.132 0.278 −0.323 −0.075 −0.268 −0.422 −0.060 −0.537 0.189 −0.596

−0.106 0.275 −0.288 0.127 −0.462 −0.217 −0.406 −0.563 −0.217 −0.660 0.044 −0.713

0.218 0.554 0.012 0.424 −0.181 0.080 −0.125 −0.289 0.112 −0.410 0.330 −0.476

−0.419 −0.011 −0.603 −0.194 −0.764 −0.527 −0.697 −0.837 −0.518 −0.927 −0.258 −0.952

0.504 0.829 0.327 0.696 0.112 0.390 0.192 −0.030 0.420 −0.161 0.617 −0.215

25.392 28.819 25.704 27.872 25.798 27.536 23.500 20.981 22.401 22.719 19.296 22.070

PPC

0.550 0.530 0.540 0.540 0.530 0.550 0.440 0.570 0.440 0.550 0.420 0.550

For example, Wang et al. (2006) used a decomposition procedure based on footprint models and ecosystem models to distinguish fluxes belonging to different stand types. On the contrary, in this study, we analyzed the spatial heterogeneity within the ecosystem. Indeed, the ecosystem analyzed can be considered homogeneous (Shao et al., 2017). Furthermore, while the approach of Wang et al. (2006) requires prior knowledge on the landscape within the footprint of a tower, our approach does not require a priori information, and heterogeneity is estimated directly from high resolution images through texture features. Overall, Dissimilarity and Entropy were found to be the most suitable in describing spatial heterogeneity. While both texture features are related to edges and therefore to rapid changes in pixel value (i.e., EVI) across a small distance, Entropy is higher when the pattern of heterogeneity is irregular (non-linear edges). Despite the results obtained here, texture feature can be highly site-specific (Hall-Beyer, 2017); therefore, we strongly suggest to calculate all the texture features related to heterogeneity (i.e., Contrast, Dissimilarity and Entropy) in future studies.

features accurately measured the different levels of heterogeneity of the landscape surrounding the flux towers. 4.1. Texture features as descriptors of landscape heterogeneity A positive relationship is shown between the percentage of footprint cover and landscape texture features, with the footprint cover acting as a function of heterogeneity during stable conditions. The Footprint %sector variable presented a higher relationship for both texture features related to variability (Contrast and Dissimilarity) within a more heterogeneous area (i.e., MDW flux tower) (Fig. 5). This finding is consistent with the well-known footprint characteristics (Schmid and Lloyd, 1999) and in accordance with a more recent study evaluating experimental flux footprint models (Heidbach et al., 2017). At night, the footprint area increased considerably and was most likely heterogeneous. Credible effects (HDIlow and HDIhigh not overlapping 0) (Table 3) were found in 3 of the 6 Fc fluxes models. Within the MDW area, the Fc flux established a negative relationship with Contrast under stable and unstable conditions and with Dissimilarity and Entropy texture features under stable conditions (Fig. 5). This can lead to a substantial overestimation of CO2 uptake across heterogeneous areas during the day. Errors in CO2 flux measurements can become magnified when propagated, if the aim is to scale up measurements using, for instance, medium-resolution remote sensing data. As a case in point, Ran et al. (2016) found that upscaling fluxes over a heterogeneous agricultural landscape in northwestern China using single towers can result in an overestimation of net ecosystem productivity ranging from 14% to 33%. The influence of landscape heterogeneity on CO2 fluxes has been studied to understand the biases due to differences between ecosystems.

4.2. Uncertainty and limitations In this study we used the analytical footprint model of Kormann and Meixner (2001). A considerable number of footprint models have been developed over time. Most notably, analytical models, Lagrangian stochastic dispersion models and large-eddy simulations are some of the most frequently used. Analytical models, such as the one developed by Schmid (1994), or the model used in our study, are based on the diffusion model of Van Ulden (1978). Although analytical models have

Table 4 Results of the Bayesian modeling for the tower TPL. For each model the Posterior predictive check (PPC), 1 st, 3rd quartiles, median, low and posterior highest density interval for 95% probability (HDIlow and HDIhigh) for the parameter β1 and median value for the parameter ν, are presented. y (Flux measurement)

Footprint%sector Footprint%sector Footprint%sector Footprint%sector Footprint%sector Footprint%sector Fcsector Fcsector Fcsector Fcsector Fcsector Fcsector

x (Texture feature)

Contrast Contrast Dissimilarity Dissimilarity Entropy Entropy Contrast Contrast Dissimilarity Dissimilarity Entropy Entropy

Time

Day Night Day Night Day Night Day Night Day Night Day Night

β1

ν st

rd

Median

1 quartile

3

−0.104 −0.106 −0.140 −0.153 −0.133 −0.115 0.069 −0.131 0.046 −0.178 −0.008 −0.145

−0.197 −0.251 −0.236 −0.299 −0.256 −0.253 −0.014 −0.275 −0.030 −0.311 −0.085 −0.270

0.004 0.020 −0.038 −0.027 −0.018 0.014 0.141 0.004 0.123 −0.029 0.081 −0.007

149

quartile

HDIlow

HDIhigh

Median

−0.441 −0.549 −0.483 −0.570 −0.504 −0.534 −0.218 −0.568 −0.236 −0.601 −0.283 −0.573

0.250 0.288 0.222 0.268 0.227 0.284 0.319 0.282 0.303 0.256 0.299 0.260

2.434 24.505 19.941 2.533 2.265 2.214 18.068 2.671 5.580 2.092 6.609 23.072

PPC

0.750 0.600 0.750 0.610 0.730 0.620 0.150 0.590 0.150 0.590 0.150 0.610

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

Fig. 4. Example of a posterior predictive check (PPC) plot for: (a) meadow steppe (MDW), and (b) typical steppe (TPL) model. The plot for both towers is generated from the posterior probability densities of the Bayesian model with the Fcsector variable as y and the Entropy variable as x. The plots compare actual and simulated data in order to visualize the predictive capability of the model. In the MDW model, the actual and simulated data exhibit a higher correspondence, resulting in a PPC of 0.42. In the TPL model, the actual and simulated data reveal a higher discrepancy, resulting in a PPC of 0.15.

Fig. 5. Scatterplots of satellite-derived texture metrics (Contrast, Dissimilarity and Entropy) as a function of the EC derived measurements (Fc and Footprint percentage) for the observations selected using cluster analysis for MDW tower during daytime (a) and nighttime (b). The dotted line represents the linear relationship while r is the Pearson correlation coefficient.

Fig. 6. Scatterplots of Satellite-derived texture metrics (Contrast, Dissimilarity, and Entropy) as a function of the EC derived measurements (Fc and Footprint percentage) for the observations selected using cluster analysis for TPL tower during daytime (a) and nighttime (b). The dotted line represents the linear relationship while r is the Pearson correlation coefficient.

150

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

Fig. 7. Boxplots of the β1 coefficient for all the Bayesian models developed for the MDW tower. Each box represents the 1st and 3rd quantiles of the coefficient distribution of the selected texture features (Contrast, Dissimilarity and Entropy). The horizontal line represents the distribution median, the whiskers indicate the limits of the 1.5*interquartile range, while the solid circles represent the outliers. Fc is CO2 flux and Footprint% is the footprint percentage in the sector.

Fig. 8. Boxplots of the β1 coefficient for all the Bayesian models developed for the TPL tower. Each box represents the 1st and 3rd quantiles of the coefficient distribution of the selected texture features (Contrast, Dissimilarity, and Entropy). The horizontal line represents the distribution median, the whiskers indicate the limits of the 1.5*interquartile range, while the solid circles represent the outliers. Fc is CO2 flux and Footprint% is the footprint percentage in the sector.

However, such models require high computational costs compared to analytical models. On the other hand, large-eddy simulation models compute flow characteristics using direct numeric simulations without requiring known flow characteristics. These models are very accurate, but can be costly in terms of computations (Sagaut, 2006). In this study, we tested our approach on the analytical footprint model of Kormann and Meixner due to its widespread use, the very low computational cost and for the fact that, differently from more complex models (e.g.,

customarily been preferred in practical applications due to lower computational times, such models do not account for heterogeneity. The Lagrangian stochastic model calculates the trajectories of a given particle from the source to the sensor (Horst and Weil, 1992), or vice versa, in a backward time frame (Kljun et al., 2002) and then extracts samples at sensor height. Backward Lagrangian stochastic models do not require homogeneous conditions and, hence, could potentially be used over heterogeneous terrain if flow characteristics are known. 151

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al.

measurements (i.e., chamber-based measurements), our findings can be used in EC fluxes scaling-up more recent approaches (i.e., environmental response functions or ERFs). ERFs, originally developed for aircraft-based EC measurements (Metzger et al., 2013), has been recently found to be suitable in scaling-up flux tower measurements (Xu et al., 2017). ERFs require biophysical surface drivers to describe spatial variations which are usually extracted from moderate resolution data. In this context, we think that the texture features proposed here could be suitable as an additional biophysical surface driver to describe heterogeneity.

Backward Lagrangian stochastic models), it does not account for heterogeneity. For this reason, the approach used in our study has a great potential to be applied for the broader FLUXNET community. In our work, a source of uncertainty can be identified in the use of half-hourly observations with no time averaging. Numerous authors have investigated the importance of the averaging period finding that a longer period would reduce EC flux measurement random errors (e.g., Richardson et al., 2006). Time-averaging has also implications in the energy balance closure problem. For example, Finnigan et al. (2003) found that averaging longer time periods could close the energy balance. On the other hand, is important to note that our study used heterogeneity as a treatment, which should have an impact on the energy balance closure. For example, Shao et al. (2014) conducted an experiment at three grassland sites and detected the spatial variability of soil heat fluxes, which is connected with the community heterogeneity. Within this study, we have provided further evidence of and quantified this spatial community variability and its effects on EC fluxes.

5. Conclusions In order to conduct flux measurements using the EC technique, the area surrounding a tower must be flat, large, and homogeneous, although very rarely do flux tower sites meet such requirements. In this study, we tested the hypothesis that uncertainties in EC flux measurements are partially dependent on landscape heterogeneity within the footprint of an EC tower. To accomplish this objective we modeled fluxderived measurements obtained in a grassland biome, including footprint percentage and Fc as a function of landscape heterogeneity, which in turn was derived from texture features extracted from high-resolution satellite images. We found that higher levels of landscape heterogeneity have an impact on flux measurements, especially during stable conditions. Our results provide a starting point for integrating highresolution remote sensing data in the correction procedures of EC flux tower estimates. However, further research is warranted to evaluate the potential of such an approach in a larger number of tower sites. Indeed, high-resolution satellite images would be insufficient in more complex tower sites such as forest ecosystems, where essential heterogeneity features (e.g., vertical forest structure) are unlikely to be described using two-dimensional images. Therefore, additional studies should focus on the integration of three-dimensional high-resolution data (e.g., LiDAR) in the modeling procedure.

4.3. Further research We investigated the influence of heterogeneity on EC flux measurements in a grassland biome. Although very little research has been conducted to examine this issue in comparable biomes, a relatively large number of studies have been conducted in forested areas (Froelich and Schmid, 2002; Lee and Hu, 2002; Chu et al., 2013; Griebel et al., 2016). The analytical approach used in this investigation demonstrates a distinct relationship between heterogeneity and EC flux measurements. However, we realize that such an approach can be difficult to apply in the context of forests because identifying properties (e.g., forest vertical structure, canopy height) is a difficult task using optical remote sensing data. Yet, active sensors like LiDAR could be an excellent choice in estimating forest metrics (Giannico et al., 2016), and can be a valid alternative to optical data when applying an analytical approach similar to ours. Heterogeneity features derived from high-resolution remote sensing data should be included as ancillary data in EC flux tower networks (e.g., AmeriFlux, AsiaFlux, USCCC, and FLUXNET). Such networks contain auxiliary information about the site in addition to flux data. For example, FLUXNET – the largest dataset of CO2, water, and energy fluxes available (Chu et al., 2017) – provides a large number of siterelated variables coupled with flux data. Information on topography (e.g., aspect and slope), vegetation (e.g., phenology and plot conditions) and soil chemistry is usually provided. Yet landscape features derived from high-resolution remote sensing are seldom available. For example, texture analysis to determine spatial heterogeneity at the plot level is restricted, probably owing to lack of availability of high-resolution data. Landscape heterogeneity is indeed often limited to the Landsat classifications with a 30-m resolution. However, the relative ease of hyper-spatial data acquisition across flux tower sites has been made possible by increased competition among commercial high-resolution data providers, which have reduced costs. Next-generation remote sensing technologies, such as Cube-Sats, further suggest that the cost of purchase will be reduced in the future (Woellert et al., 2011). The proliferation of drone-mounted multispectral and LiDAR sensors will also facilitate the analysis of fine-scale heterogeneity studies around current and proposed flux tower sites (Zhang et al., 2016). More importantly, connecting landscape heterogeneity within the footprint of a flux tower will open many opportunities for the scientific community to examine the direct contributions and processes of biophysical and ecosystem (e.g., physiological, soil, and ecosystem processes) to the net ecosystem exchanges of materials and energy from the EC towers. Our findings demonstrated the dependence of EC flux measurements on heterogeneity. Although more research is needed to better understand this dependence using ground truth data derived from alternative

Acknowledgements This study was supported by the “Dynamics of Coupled Natural and Human Systems (CNH)” Program of the NSF (#1313761), the LCLUC program of NASA (NNX15AD10G), and the USCCC. We appreciate the careful editing by Connor Crank on earlier versions of the manuscript. We thank an anonymous reviewer for providing very valuable suggestions to improve the quality of the manuscript. References Chen, B., Black, T.A., Coops, N.C., Hilker, T., Trofymow (Tony), J.A., Morgenstern, K., 2008. Assessing tower flux footprint climatology and scaling between remotely sensed and eddy covariance measurements. Bound.-Layer Meteorol. 130, 137–167. http://dx.doi.org/10.1007/s10546-008-9339-1. Chen, B., Coops, N.C., Fu, D., Margolis, H.A., Amiro, B.D., Black, T.A., Arain, M.A., Barr, A.G., Bourque, C.P.-A., Flanagan, L.B., Lafleur, P.M., McCaughey, J.H., Wofsy, S.C., 2012. Characterizing spatial representativeness of flux tower eddy-covariance measurements across the Canadian Carbon Program Network using remote sensing and footprint analysis. Remote Sens. Environ. 124, 742–755. http://dx.doi.org/10.1016/ j.rse.2012.06.007. Chu, H.-S., Liang, N., Lai, C.-W., Wu, C.-C., Chang, S.-C., Hsia, Y.-J., 2013. Topographic effects on CO2 flux measurements at the Chi-Lan Mountain forest site. Taiwan J. For. Sci. 28, 1–16. Chu, H., Chen, J., Gottgens, J.F., Ouyang, Z., John, R., Czajkowski, K., Becker, R., 2014. Net ecosystem methane and carbon dioxide exchanges in a Lake Erie coastal marsh and a nearby cropland. J. Geophys. Res. Biogeosci. 119http://dx.doi.org/10.1002/ 2013JG002520. 2013JG002520. Chu, H., Baldocchi, D.D., John, R., Wolf, S., Reichstein, M., 2017. Fluxes all of the time? A primer on the temporal representativeness of FLUXNET. J. Geophys. Res. Biogeosci. 122http://dx.doi.org/10.1002/2016JG003576. 2016JG003576. Davis, K.J., Bakwin, P.S., Yi, C., Berger, B.W., Zhao, C., Teclaw, R.M., Isebrands, J.G., 2003. The annual cycles of CO2 and H2O exchange over a northern mixed forest as observed from a very tall tower. Glob. Change Biol. 9, 1278–1293. http://dx.doi.org/ 10.1046/j.1365-2486.2003.00672.x. Finnigan, J.J., Clement, R., Malhi, Y., Leuning, R., Cleugh, H.A., 2003. A re-evaluation of

152

Agricultural and Forest Meteorology 260–261 (2018) 144–153

V. Giannico et al. long-term flux measurement techniques part I: averaging and coordinate rotation. Bound.-Layer Meteorol. 107, 1–48. http://dx.doi.org/10.1023/A:1021554900225. Forman, R.T.T., 1995. Some general principles of landscape and regional ecology. Landsc. Ecol. 10, 133–142. http://dx.doi.org/10.1007/BF00133027. Fraley, C., Raftery, A.E., Murphy, T.B., Scrucca, L., 2012. mclust v. 4 for R: Normal Mixture Modeling for Modelbased Clustering, Classification, and Density Estimation. Dept. Stat. Univ. Wash. Tech. Rep. 587. Froelich, N.J., Schmid, H.P., 2002. 10.8 An investigation of advection and gully flows in complex forested terrain. Twenty Fifth Conference on Agricultural and Forest Meteorology. The Society, p. 175. Gelman, A., Rubin, D.B., 1992. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 457–472. Giannico, V., Lafortezza, R., John, R., Sanesi, G., Pesola, L., Chen, J., 2016. Estimating stand volume and above-ground biomass of urban forests using LiDAR. Remote Sens. 8, 339. http://dx.doi.org/10.3390/rs8040339. Goulden, M.L., Munger, J.W., Fan, S.-M., Daube, B.C., Wofsy, S.C., 1996. Measurements of carbon sequestration by long-term eddy covariance: methods and a critical evaluation of accuracy. Glob. Change Biol. 2, 169–182. http://dx.doi.org/10.1111/j. 1365-2486.1996.tb00070.x. Griebel, A., Bennett, L.T., Metzen, D., Cleverly, J., Burba, G., Arndt, S.K., 2016. Effects of inhomogeneities within the flux footprint on the interpretation of seasonal, annual, and interannual ecosystem carbon exchange. Agric. For. Meteorol. 221, 50–60. http://dx.doi.org/10.1016/j.agrformet.2016.02.002. Hall-Beyer, M., 2017. Practical guidelines for choosing GLCM textures to use in landscape classification tasks over a range of moderate spatial scales. Int. J. Remote Sens. 38, 1312–1338. http://dx.doi.org/10.1080/01431161.2016.1278314. Haralick, R.M., Shanmugam, K., Dinstein, I., 1973. Textural features for image classification. IEEE Trans. Syst. Man. Cybern. SMC-3 610–621. http://dx.doi.org/10.1109/ TSMC.1973.4309314. Heidbach, K., Schmid, H.P., Mauder, M., 2017. Experimental evaluation of flux footprint models. Agric. For. Meteorol. 246, 142–153. http://dx.doi.org/10.1016/j.agrformet. 2017.06.008. Horst, T.W., Weil, J.C., 1992. Footprint estimation for scalar flux measurements in the atmospheric surface layer. Bound.-Layer Meteorol. 59, 279–296. http://dx.doi.org/ 10.1007/BF00119817. John, R., Chen, J., Kim, Y., Ou-yang, Z., Xiao, J., Park, H., Shao, C., Zhang, Y., Amarjargal, A., Batkhshig, O., Qi, J., 2015. Differentiating anthropogenic modification and precipitation-driven change on vegetation productivity on the Mongolian plateau. Landsc. Ecol. 31, 547–566. http://dx.doi.org/10.1007/s10980-015-0261-x. Jones, M.C., Faddy, M.J., 2003. A skew extension of the t-distribution, with applications. J. R. Stat. Soc. Ser. B Stat. Methodol. 65, 159–174. http://dx.doi.org/10.1111/14679868.00378. Kayitakire, F., Hamel, C., Defourny, P., 2006. Retrieving forest structure variables based on image texture analysis and IKONOS-2 imagery. Remote Sens. Environ. 102, 390–401. http://dx.doi.org/10.1016/j.rse.2006.02.022. Kellner, K., 2016. jagsUI: A Wrapper Around “rjags” to Streamline “JAGS” Analyses. Kljun, N., Rotach, M.W., Schmid, H.P., 2002. A three-dimensional backward Lagrangian footprint model for a wide range of boundary-layer stratifications. Bound.-Layer Meteorol. 103, 205–226. http://dx.doi.org/10.1023/A:1014556300021. Kormann, R., Meixner, F.X., 2001. An analytical footprint model for non-neutral stratification. Bound.-Layer Meteorol. 99, 207–224. http://dx.doi.org/10.1023/ A:1018991015119. Kruschke, J., 2014. Doing Bayesian Data Analysis: A Tutorial With R, JAGS, and Stan. Academic Press. Lange, K.L., Little, R.J.A., Taylor, J.M.G., 1989. Robust statistical modeling using the t distribution. J. Am. Stat. Assoc. 84, 881–896. http://dx.doi.org/10.1080/01621459. 1989.10478852. Lee, X., 1998. On micrometeorological observations of surface-air exchange over tall vegetation. Agric. For. Meteorol. 91, 39–49. http://dx.doi.org/10.1016/S01681923(98)00071-9. Lee, X., Hu, X., 2002. Forest-air fluxes of carbon, water and energy over non-flat terrain. Bound.-Layer Meteorol. 103, 277–301. http://dx.doi.org/10.1023/ A:1014508928693. Lee, X., Massman, W., Law, B., 2006. Handbook of Micrometeorology: A Guide for Surface Flux Measurement and Analysis. Springer Science & Business Media. Massman, W.J., Lee, X., 2002. Eddy covariance flux corrections and uncertainties in longterm studies of carbon and energy exchanges. Agric. For. Meteorol. FLUXNET 2000 Synth. 113, 121–144. http://dx.doi.org/10.1016/S0168-1923(02)00105-3. Metzger, S., Junkermann, W., Mauder, M., Butterbach-Bahl, K., Trancón y Widemann, B., Neidl, F., Schäfer, K., Wieneke, S., Zheng, X.H., Schmid, H.P., Foken, T., 2013. Spatially explicit regionalization of airborne flux measurements using environmental response functions. Biogeosciences 10, 2193–2217. http://dx.doi.org/10.5194/bg10-2193-2013. Moncrieff, J., Clement, R., Finnigan, J., Meyers, T., 2004. Averaging, detrending, and filtering of eddy covariance time series. In: Lee, X., Massman, W., Law, B. (Eds.),

Handbook of Micrometeorology, Atmospheric and Oceanographic Sciences Library. Springer, Netherlands. http://dx.doi.org/10.1007/1-4020-2265-4_2. pp. 7–31. Neftel, A., Spirig, C., Ammann, C., 2008. Application and test of a simple tool for operational footprint evaluations. Environ. Pollut. 152, 644–652. http://dx.doi.org/10. 1016/j.envpol.2007.06.062. Ozdemir, I., Karnieli, A., 2011. Predicting forest structural parameters using the image texture derived from WorldView-2 multispectral imagery in a dryland forest, Israel. Int. J. Appl. Earth Observ. Geoinform. 13, 701–710. http://dx.doi.org/10.1016/j.jag. 2011.05.006. Paw, U.K.T., Baldocchi, D.D., Meyers, T.P., Wilson, K.B., 2000. Correction of eddy-covariance measurements incorporating both advective effects and density fluxes. Bound.-Layer Meteorol. 97, 487–511. http://dx.doi.org/10.1023/A:1002786702909. Plummer, Martyn, 2016. rjags: Bayesian Graphical Models Using MCMC. Plummer, M., Best, N., Cowles, K., Vines, K., 2006. CODA: convergence diagnosis and output analysis for MCMC. R. News 6, 7–11. R Core Team, 2016. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Ran, Y., Li, X., Sun, R., Kljun, N., Zhang, L., Wang, X., Zhu, G., 2016. Spatial representativeness and uncertainty of eddy covariance carbon flux measurements for upscaling net ecosystem productivity to the grid scale. Agric. For. Meteorol. OasisDesert Syst. 230–231, 114–127. http://dx.doi.org/10.1016/j.agrformet.2016.05.008. Richardson, A.D., Hollinger, D.Y., Burba, G.G., Davis, K.J., Flanagan, L.B., Katul, G.G., William Munger, J., Ricciuto, D.M., Stoy, P.C., Suyker, A.E., Verma, S.B., Wofsy, S.C., 2006. A multi-site analysis of random error in tower-based measurements of carbon and energy fluxes. Agric. For. Meteorol. 136, 1–18. http://dx.doi.org/10.1016/j. agrformet.2006.01.007. Sagaut, P., 2006. Large Eddy Simulation for Incompressible Flows: An Introduction. Springer Science & Business Media. Schmid, H.P., 1994. Source areas for scalars and scalar fluxes. Bound.-Layer Meteorol. 67, 293–318. Schmid, H.P., 2002. Footprint modeling for vegetation atmosphere exchange studies: a review and perspective. Agric. For. Meteorol., FLUXNET 2000 Synth. 113, 159–183. http://dx.doi.org/10.1016/S0168-1923(02)00107-7. Schmid, H.P., Lloyd, C.R., 1999. Spatial representativeness and the location bias of flux footprints over inhomogeneous areas. Agric. For. Meteorol. 93, 195–209. http://dx. doi.org/10.1016/S0168-1923(98)00119-1. Shao, C., Li, L., Dong, G., Chen, J., 2014. Spatial variation of net radiation and its contribution to energy balance closures in grassland ecosystems. Ecol. Process. 3 (7). http://dx.doi.org/10.1186/2192-1709-3-7. Shao, C., Chen, J., Chu, H., Lafortezza, R., Dong, G., Abraha, M., Batkhishig, O., John, R., Ouyang, Z., Zhang, Y., Qi, J., 2017. Grassland productivity and carbon sequestration in Mongolian grasslands: the underlying mechanisms and nomadic implications. Environ. Res. 159, 124–134. http://dx.doi.org/10.1016/j.envres.2017.08.001. Sibanda, M., Mutanga, O., Rouget, M., Kumar, L., 2017. Estimating biomass of native grass grown under complex management treatments using WorldView-3 spectral derivatives. Remote Sens. 9, 55. http://dx.doi.org/10.3390/rs9010055. Tsionas, E.G., 2002. Bayesian inference in the noncentral student-t model. J. Comput. Graph. Stat. 11, 208–221. http://dx.doi.org/10.1198/106186002317375695. Twine, T.E., Kustas, W.P., Norman, J.M., Cook, D.R., Houser, P.R., Meyers, T.P., Prueger, J.H., Starks, P.J., Wesely, M.L., 2000. Correcting eddy-covariance flux underestimates over a grassland. Agric. For. Meteorol. 103, 279–300. http://dx.doi.org/10. 1016/S0168-1923(00)00123-4. Van Ulden, A.P., 1978. Simple estimates for vertical diffusion from sources near the ground. Atmosp. Environ. 1967 (12), 2125–2129. http://dx.doi.org/10.1016/00046981(78)90167-1. Wang, W., Davis Kenneth, J., Cook Bruce, D., Butler Martha, P., Ricciuto Daniel, M., 2006. Decomposing CO2 fluxes measured over a mixed ecosystem at a tall tower and extending to a region: a case study. J. Geophys. Res. Biogeosci. 111. http://dx.doi.org/ 10.1029/2005JG000093. Webb, E.K., Pearman, G.I., Leuning, R., 1980. Correction of flux measurements for density effects due to heat and water vapour transfer. Q. J. R. Meteorol. Soc. 106, 85–100. http://dx.doi.org/10.1002/qj.49710644707. Wilczak, J.M., Oncley, S.P., Stage, S.A., 2001. Sonic anemometer tilt correction algorithms. Bound.-Layer Meteorol. 99, 127–150. http://dx.doi.org/10.1023/ A:1018966204465. Woellert, K., Ehrenfreund, P., Ricco, A.J., Hertzfeld, H., 2011. Cubesats: cost-effective science and technology platforms for emerging and developing nations. Adv. Space Res. 47, 663–684. http://dx.doi.org/10.1016/j.asr.2010.10.009. Xu, K., Metzger, S., Desai, A.R., 2017. Upscaling tower-observed turbulent exchange at fine spatio-temporal resolution using environmental response functions. Agric. For. Meteorol. 232, 10–22. http://dx.doi.org/10.1016/j.agrformet.2016.07.019. Zhang, J., Hu, J., Lian, J., Fan, Z., Ouyang, X., Ye, W., 2016. Seeing the forest from drones: testing the potential of lightweight drones as a tool for long-term forest monitoring. Biol. Conserv. 198, 60–69. http://dx.doi.org/10.1016/j.biocon.2016.03. 027.

153