The MODIS (Collection V005) BRDF/albedo product: Assessment of spatial representativeness over forested landscapes

of satellite retrievals from global land surface albedo and re ﬂ ectance anisotropy products is presented. This method brings together knowledge of the intrinsic biophysical properties of a measurement site, and the surrounding landscape to produce a number of geostatistical attributes that describe the overall variability, spatial extent, strength of the spatial correlation, and spatial structure of surface albedo patterns at separate seasonal periods throughout the year. Variogram functions extracted from Enhanced Thematic Mapper Plus (ETM+) retrievals of surface albedo using multiple spatial and temporal thresholds were used to assess the degree to which a given point (tower) measurement is able to capture the intrinsic variability of the immediate landscape extending to a satellite pixel. A validation scheme was implemented over a wide range of forested landscapes, looking at both deciduous and coniferous sites, from tropical to boreal ecosystems. The experiment focused on comparisons between tower measurements of surface albedo acquired at local solar noon and matching retrievals from the MODerate Resolution Imaging Spectroradiometer (MODIS) (Collection V005) Bidirectional Re ﬂ ectance Distribution Function (BRDF)/albedo algorithm. Assessments over a select group of ﬁ eld stations with comparable landscape features and daily retrieval scenarios further demonstrate the ability of this technique to identify measurement sites that contain the intrinsic spatial and seasonal features of surface albedo over suf ﬁ ciently large enough footprints for use in modeling and remote sensing studies. This approach, therefore, improves our understanding of product uncertainty both in terms of the representativeness of the ﬁ eld data and its relationship to the larger satellite pixel.


Introduction
Land surface albedo is defined as the fraction of incident solar irradiance reflected by Earth's surface over the whole solar spectrum (Dickinson, 1983). It has long been recognized as a key land surface radiative property since it determines the amount of solar radiation absorbed by the land surface. Measuring changes in surface albedo through time provides a quantitative means for tracking land surface change and measuring that change in terms that can be directly applied by regional and global modeling efforts, in particular those that evaluate and monitor energy transfer at the surface. Variations Remote Sensing of Environment 113 (2009)   in the extent of snow cover and flooding, the phenology of natural vegetation and agricultural crops, as well as other signatures from rapidly changing surface covers (e.g. burning, clearing, and tilling) are all accompanied by significant changes in surface albedo (Schaaf et al., 2008). As such, surface albedos with an absolute accuracy of 0.02-0.05 units (Henderson-Sellers & Wilson, 1983;Sellers et al., 1995) for snow-free and snow-covered land are required by climate, biogeochemical, hydrological, and weather forecast models at a range of spatial (from 10 s of meters to 5-30 km) and temporal (from daily to monthly) scales.
Surface albedo varies strongly in space and time depending on seasonal trends, soil-vegetation contrasts, canopy chemistry and structure (Ollinger et al., 2008). If the albedo of the soil surface is significantly different from above-canopy albedo, for example, any changes in vegetation cover that may expose the soil can be important (Berbet & Costa, 2003). Rainfall events can also reduce the magnitude of above-canopy albedo by as much as 50% of its dry-peak, and are more persistent if the downpour occurs in the morning than in the afternoon, because the soil surface has time to drain and dry during the night (Samain et al., 2008). Surface albedo is also controlled by the intrinsic structural and optical properties of vegetated canopies, e.g. mutual shadowing by tree crowns and specular reflection by leaves (Lucht et al., 2000). These elements are highly solar zenith angle dependent and modulated by the Bidirectional Reflectance Distribution Function (BRDF), which describes the anisotropic reflectance of natural surfaces (Martonchik et al., 2000;Nicodemus, 1977). Daily measurements of surface albedo are also influenced by changes in atmospheric watervapor content, which absorbs strongly in the NIR (0.7-5.0 µm) region, and by other changing atmospheric conditions (e.g. clouds and haze) that influence the partitioning of direct and diffuse solar irradiance. Over bright-colored soils, vegetation cover reduces surface albedo, whereas plants on dark or litter-covered soils increases surface albedo (Rechid & Jacob, 2006). These considerations underscore the importance of understanding the natural and anthropogenic influences on surface albedo for determining changes in global climate.
A number of recent studies have evaluated the consistency of various global land surface albedo products against in-situ data at different spatial and temporal scales (Chen et al., 2008;Hautecoeur & Roujean, 2007;Jin et al., 2003a,b;Liang et al., 2002;Liu et al., 2009;Rutan et al., 2009;Salomon et al., 2006). These assessments have been carried out under the general assumption that both the satellite and tower albedo measurements are being captured from the exact same spatial domain. However, the vast differences in spatial resolution create a number of challenges for direct "point-to-pixel" comparisons. Local changes in surface albedo, within and between different ecosystems, may differ as a function of spatial resolution and, consequently, introduce scaling errors. These patterns will also change seasonally and are particularly hard to establish throughout periods of rapidly changing surface conditions. Unless the underlying surface is large and perfectly homogeneous, or a sufficient number of distributed point measurements (i.e. within a single pixel are) are made during the satellite overpass, "point" measurements alone are not sufficient to validate satellite-derived albedo retrievals (Liang et al., 2002).
Understanding the intrinsic properties and anisotropic behavior of the landscape requires more thorough assessments at a diverse range of spatial and seasonal scales. This paper presents a new methodology (illustrated in Fig. 1) for validating the MODIS BRDF/albedo product that focuses on quantifying the spatial representativeness of a given measurement site. The term spatial representativeness is defined here as the degree to which a ground-based retrieval of surface albedo is able to resolve the surrounding landscape extending to the satellite footprint. This validation technique assesses the spatial and temporal characteristics of the landscape by using multispectral high spatial resolution imagery within a geostatistical framework. This approach is related to previously introduced methods which have also employed fine spatial resolution datasets to form the intermediary basis between field measurements and MODIS data (Baccini et al., 2004;Baret et al., 2006;Tian et al., 2004).
The ability to characterize the spatial representativeness of a measurement site has significant implications for the efficiency of direct "point-to-pixel" assessments of surface biophysical properties derived from satellite data. This would enable a better understanding of product uncertainty, both in terms of the quality of the model inversions (e.g. given a limited number of cloud-free satellite observations), and their ability to capture the underlying spatial variability of the surrounding landscape. The goal of this paper is to provide indicators of spatial representativeness, as opposed to spatial heterogeneity. One could argue that the observed landscape can be categorized as both spatially heterogeneous and spatially representative depending on the structural consistency and spatial distribution of the landscape elements resolved by both the ground instrument field-of-view (be it from a small 10 m flux tower or a 500 m tall radio tower) and the surrounding region extending to a satellite pixel.

AmeriFlux network
Tower-based albedo measurements were acquired over a broad range of forested landscapes that are part of the AmeriFlux network of sites and the global network, FLUXNET (Table 1). The AmeriFlux network provides continuous observations of ecosystem level exchanges of CO 2 , water and energy, spanning diurnal, seasonal, and interannual time scales and is currently composed of sites from North America, Central America, and South America (Law et al., 2002;Running et al., 1999). The field stations in this study cover a broad range in climate and growth form, including ponderosa pine in the Southwestern United States (US), broadleaf forest in the Midwest, and subtropical semidecidous broadleaf forests in the Amazon region (WWW1). Tower albedo data available from 16 AmeriFlux sites, ranging from 1 to 5 year measurement periods or 42 site-years (Table 1), include a wide range of settingsfrom recently harvested plantation stands to oldgrowth forests. Measurement sites were instrumented with Kipp and Zonen (CNR1, CM-3, or CM-6b), or Eppley-PSP tower albedometers. These sensors have been outfitted with clear domes to collect broadband albedo and radiation fluxes in the shortwave domain (SW) (0.3-2.8 µm). Data received from individual sites are reviewed and incorporated into the network-wide AmeriFlux database. The review process, as established by the Carbon Dioxide Information Analysis Center (CDIAC) (Cook et al., 2001;WWW2), includes checks for consistent units, naming conventions, reporting intervals, and reformatting to maintain consistency within the larger network-wide database. For this experiment, daily tower albedo values were derived by calculating the mean and standard deviation of all albedometer retrievals available within a 2-hour window centered at Local Solar Noon (LSN).

MODIS (Collection V005) BRDF/albedo product
The global land surface MODIS BRDF/albedo products are operationally produced at 500 m every 8 days. The algorithm relies on multi-day, cloud-free, atmospherically-corrected surface reflectances from both Terra and Aqua to sample the surface reflectance anisotropy over a 16-day period and allow for the retrieval of a reflectance anisotropy model for each pixel. These products make use of a kernel-driven linear model of the Bidirectional Reflectance Distribution Function (BRDF), known as the Ross-Thick/Li-Sparse-Reciprocal (RTLSR) kernel-driven BRDF model (Lucht et al., 1999;Roujean et al., 1992;Wanner et al., 1995Wanner et al., , 1997, which relies on the weighted sum of an isotropic parameter and two functions (or kernels) of viewing and illumination geometry (Roujean et al., 1992).
Where: θ is the solar zenith angle; ϑ is the pixel viewing zenith angle; ϕ is the relative view-sun azimuth angle; and Λ is the solar spectrum weighted center for a given MODIS spectral band. Parameter f iso (Λ) is the isotropic scattering component and equal to the bidirectional reflectance for a view zenith angle ϑ = 0 and a solar zenith angle θ = 0. Parameter f geo (Λ) is the coefficient of the Li-Sparse-Reciprocal geometric scattering kernel K geo , which is derived from surface scattering and geometric shadowcasting theory (Li & Strahler, 1992). Parameter f vol (Λ) is the coefficient for the Ross-Thick volume scattering kernel K vol , which is derived from radiative transfer models (Ross, 1981). The best fit Ross-Thick/Li-Sparse-Reciprocal (RTLSR) (Lucht et al., 2000;Wanner et al., 1995) model parameters are then retrieved for the first seven spectral bands of MODIS. This model combination has been shown to be well suited for describing the surface reflectance anisotropy of a variety of land covers that are distributed worldwide (Privette et al., 1997). The spectral acquisitions can also be combined via narrow to broadband conversion coefficients (Liang, 2001;Liang et al., 1999) to provide broadband anisotropy information and thus broadband albedos for three additional broadbands (0.3-0.7 µm, 0.7-5.0 µm, 0.3-5.0 µm).
Once an appropriate anisotropy model has been retrieved, integration over all view angles results in a directional-hemispherical reflectance (DHR) or a black-sky albedo (BSA) at any desired solar zenith angle. A further integration over all illumination angles results in a bihemispherical reflectance (BHR) under isotropic illumination, or a white-sky albedo (WSA): where: h k (θ) is the integral of the BRDF model kernels k over a given view zenith and view-sun relative azimuth angle; H k is the integral of h k over a given solar zenith angle θ; and f k are the BRDF kernel model parameters k. These albedo quantities are intrinsic to a specific location and are governed by the character and structure of its land cover. Albeit a product based on clear-sky multi-day inputs, the MODIS BRDF/albedo algorithm can characterize the surface anisotropy and capture the daily temporal and spatial variations of albedo, as long as a concurrent knowledge of the atmospheric state is available. Accordingly, by interpolating between black-sky and white-sky albedo quantities, the following equation can be used to compute MODIS blue-sky (or actual) (or instantaneous) albedo (Lewis & Barnsley, 1994): where: D(θ s, τ) is the fraction of diffuse skylight, which is a function of aerosol optical depth and solar zenith angle (SZA): The latter term involves the diffuse sky radiation, E diffuse (Λ , τ), which is the downwelling diffuse irradiance under a perfectly absorbing lower boundary (i.e. no ground interactions). Level 2 aerosol optical depths (at 550 nm) available from the MODIS 10 km (Collection V005) Aerosol  Davis et al. (2003) Products (Remer et al., 2005), were acquired to provide daily estimates of D(θ s, τ), and thus obtain coincident MODIS blue-sky albedo retrievals at local solar noon (LSN) for each measurement site. In addition to the spectral and broadband BRDF model and albedo quantities themselves, the MODIS V005 BRDF/albedo product provides extensive quality information. A full retrieval of the parameters for the RTLSR BRDF model is attempted if sufficient high quality, cloud-free well distributed directional samples are acquired during a 16-day period. A backup algorithm is used if insufficient directional samples survive the screening process or if a robust full retrieval cannot be made. Salomon et al. (2006) found that the combination of MODIS Aqua and Terra sensor data has increased the occurrence of high quality fully data-driven retrievals and thus reduced the product's reliance upon a-priori determinations of the underlying surface anisotropy used by the backup algorithm.

Experimental methods
Recent advances in spatially-explicit environmental research have been increasingly aided by the field of geostatistics. One of the most efficiently used tools for describing the spatial continuity of primary biophysical properties is the variogram model (Carroll & Cressie, 1996;Davis, 1986;Isaaks & Srivastava, 1989). Remote sensing data enables efficient monitoring of the properties of variogram models and can further reveal highly interesting patterns of spatial variability and scaling information (Dent & Grimm, 1999). Using variogram models, a number of geostatistical attributes can be derived to quantify the spatial representativeness of a given measurement site over the spatial scales of MODIS observations (i.e. 0.5 km 2 to 1.0 km 2 ). The following two sections explain the processing steps involved in deriving these measures.

Estimating variogram model parameters from ETM+ data
The methodology for estimating broadband albedos based on empirical relations between surface total shortwave albedo measurements and Enhanced Thematic Mapper Plus (ETM+) observations was introduced by Liang (2001) and further assessed in Liang et al. (2003). Using high spatial resolution scenes, this technique requires a series of processing steps: 1. Top-of-Atmosphere (TOA) spectral reflectances were derived from radiometric calibration of Level 1 Digital Numbers (DNs) by using solar irradiance tables and determining the view and illumination conditions at the ETM+ overpass. 2. Directional surface reflectances were then produced by atmospherically correcting the TOA spectral reflectances using the 6S radiative transfer code (Kotchenova et al., 2006;Vermote et al., 1997), as well as ground or satellite retrievals of aerosol optical properties measured at the ETM+ overpass. 3. Narrowband-to-broadband (NTB) albedo conversions (Liang, 2001;Liang et al., 2003) were used to produce surface albedos equivalent to those measured by MODIS and field albedometers.
The methodology for deriving variogram functions to analyze surface albedos using ETM+ subsets as intermediaries between ground and MODIS footprints was recently introduced in Susaki et al. (2007). In this work, the variogram estimator, γ E (h), is used to obtain half the average-squared-difference between albedo values that are within certain distance classes or bins defined by multiples of 30 m (i.e. the nominal resolution of ETM+): where: z x is the surface albedo at pixel location x; and z x + h is the surface albedo at another pixel within a lag distance h. As a rule of thumb, the maximum lag distance used in each variogram is constrained by the half maximum distance of the prescribed subset and must be a factor of the minimum lag (i.e. 30 m). Thus, for a 1.0 km 2 subset h max = 690 m, for a 1.5 km 2 subset h max = 1050 m, and for a 2.0 km 2 subset h max =1410 m. The variogram model parametersrange, sill, and nugget effectcan then be modified to fit the isotropic spherical variogram model (Materon, 1963) to the variogram estimator: An illustration of isotropic spherical variogram models, and the relevant model parameters, is presented in Fig. 2. The range (a) defines the distance from a point beyond which there is no further correlation of a biophysical property associated with that point. It has also been described as the average patch size of the landscape (Cooper et al., 1997); i.e., a region that differs from its surroundings, but is not necessarily internally homogeneous. The sill (c) is the ordinate value of the range at which the variogram levels off to an asymptote. Thus, the sill describes the maximum semivariance, while the range describes the separation distance at which this occurs. If it appears that the variogram does not reach zero variance at h = 0, the apparent positive intercept is called the nugget effect (c 0 ). This is equivalent to the variance of lags smaller than that of the sampling distance. This in turn depends on variance associated with small scale variability, measurement errors, or a combination of these (Noreus et al., 1997). While the MODIS V005 BRDF/albedo product is reported on a 500 m grid, the retrievals used to obtain the BRDF (and thus the albedo) are based on many observations. Even in the best case scenario of a nadir retrieval centered directly on a 500 m pixel, the MODIS point spread function for just the sensor in the across track direction includes an area as wide as 1 km (Tan et al., 2005;Wolfe et al., 2002). Thus, MODIS observations are not necessarily centered on the precise location of the pixel in question. Such retrievals may actually be from a larger spatial domain, depending on sensor view geometry: where: P is the size of a MODIS pixel [m]; w is the MODIS instantaneous field-of-view for the reflective bands (0.6571 mrad); d is the altitude of the sensor (705 km); and θ v is the sensor's view zenith angle.
By examining the variogram model parameters at different spatial resolutions, the spatial characteristics of a given measurement site can be compared against the larger landscapes extending to a satellite pixel. Accordingly, the routines presented in Susaki et al. (2007) have been extended in this paper to further investigate the response of variogram functions using multiple spatial and seasonal thresholds. In particular, a new set of geostatistical attributes (described in Section 3.2) was created to analyze the change in variogram model parameters as a function of increased window-size (i.e. from 1.0 km 2 to 2.0 km 2 squared subsets). These spatial thresholds were specifically chosen to account for all the possible surface covers that may contribute to the directional signatures acquired by the MODIS instrument throughout a 16-day period. This is a conservative approach in the sense that it ensures that measurement sites are representative of relatively large areas to justify their suitability for evaluating the MODIS 500 m albedo retrievals. One reason this approach is conservative is that the MODIS BRDF/albedo algorithm has been shown to be a good outlier detector, and internally employs a quality assurance routine that discards the furthermost outliers from retrievals since they are obviously contaminated or of a different ground location than the rest (Schaaf et al., 2002).  An example of the variogram functions and relevant model parameters (range, sill, nugget effect) for three ETM+ subsets of size 1.0 km 2 , 1.5 km 2 , and 2.0 km 2 are shown in Fig. 3B for the Harvard Forest Environmental Measurement Site (HFEMS). The vegetation surrounding the HFEMS (seen if Fig. 3A) is predominantly a mosaic of transitional hardwoods including red oak and red maple; together with some hemlock stands, including white pine and red pine. The relatively sparse understory contains primarily red maple, yellow birch and hemlock saplings plus woody shrubs such as blueberry, and witch hazel. Fig. 3B shows that the variogram estimator (point values) for each of the ETM+ subsets are more closely aligned over shorter separation distances (h b 200 m) than for longer ones (h N 200 m). The spatial variability over Harvard Forest is also more pronounced over the larger squared regions (1.5 km 2 and 2.0 km 2 ), given the chances of finding significantly different land cover types at greater separation distances are much higher. Finally, for each of the subsets, the variogram estimator reaches an asymptote, or a constant variance among spatially uncorrelated samples, near the sample variance (var) (solid straight lines).

Relative coefficient of variation -R CV
A number of theoretical, empirical, and statistical approaches have been presented by previous studies with the common goal to characterize the spatial dependency of landscape patterns and processes (Cooper et al., 1997;Dent & Grimm, 1999;Isaaks & Srivastava, 1989;Journel & Huijbregts, 1978). Analogously, the spatial representativeness of tower albedo measurements can be characterized by producing a number of geostatistical attributes that describe the overall variability, spatial extent, strength of the spatial correlation, and spatial structure of surface albedo patterns for a given measurement site (Fig. 4).
The coefficient of variation (CV) is defined by the ratio of the standard deviation to the mean. It is a useful measure of the relative spread in the data and provides an estimate of overall variability that is independent of spatial scale. A good way to utilize this measure, for the purposes of this study, is to measure the change in CV as a function of increased field-of-view: The relative coefficient of variation, R CV , is the difference between two major terms: (1) CV x is the coefficient of variation obtained from a 1.0 km 2 ETM+ albedo subset that is centered on a given measurement site; and (2) CV 1.5x is the coefficient of variation obtained over a footprint that is 1.5 times the size of CV x (i.e. 1.5 km 2 ). If a measurement site is spatially representative, then the overall variability between the internal components of the measurement site and the adjacent landscape should be similar in magnitude and R CV should approach zero. Conversely, a large positive or negative value for R CV would correspond to higher variance within the 1.0 km 2 area immediately surrounding the tower than the larger 1.5 km 2 area. The sign of R CV thus provides a first-order estimate of the primary sources that drive the overall variability, i.e. negative for internal variations associated to CV x and positive for external variations associated to CV 1.5x . Because R CV is based on areal-mean values, comparisons at different spectral ranges (e.g. VIS 0.3-0.7 µm and NIR 0.7-5.0 µm) or between other vegetation parameters (e.g. NDVI, LAI, and LST), would thus produce comparable measures of overall variability.

Scale requirement index -R SE
While it is very challenging to resolve the footprint of a measurement site from satellite retrievals, it is essential to provide a measure of the range (a) of surface albedo with respect to both the ground instrument's field-of-view and the surrounding landscape. In principle, surface albedo retrievals that are sampled at distances greater than the variogram range can be treated as spatially independent and consequently be applied in direct assessments between ground measurements and MODIS retrievals (Susaki et al., 2007). Accordingly, a new geostatistical attribute, defined here as the scale requirement index Fig. 4. The geostatistical attributes of spatial representativeness can be used to describe the overall variability (R CV ), spatial extent (R SE ), strength of the spatial correlation (R ST ), and spatial structure (R ST ) of surface albedo patterns for a given measurement site. These measures improve our understanding of product uncertainty both in terms of the representativeness of the field data and its relationship to the larger MODIS pixel.
(R SE ), examines the variogram range by using two spatial thresholds with respect to the true spatial extent of a given measurement site: where: ]. This will be the case for a tall tower instrument setup which is overlooking a ground (circular) footprint of equal or larger extent than the squared areas defined by the range parameters at 1.0 km 2 and 1.5 km 2 (i.e. g ≥ a x and g ≥ a 1.5x ).

Relative strength of the spatial correlation -R ST
Another good indicator of spatial variability is the strength of spatial dependence (or autocorrelation) over the range (SD) (Cooper et al., 1997): This quantity measures the slope of the ascending limb of the variogram estimator after standardizing the semivariance by dividing through by its maximum value. By standardizing the semivariance before calculating the slope, a measure that depends on the range (a) and the nugget variance (c 0 ) (Cooper et al., 1997;Dent & Grimm, 1999) can be obtained. In using SD as an indicator of spatial representativeness, this measure can be used to assess changes in surface albedo over short separation distances: where: Unlike SD, the relative strength of the spatial correlation, R ST , provides a spatially-explicit representation of where the most different source (s) of albedo variability (e.g. lakes, small ponds, clear-cuts, and bare areas) are likely to be situated with respect to the measurement site. Thus, a negative value for R ST can be attributed to internal outliers associated to ST x , whereas positive values are attributed to external outliers associated to ST 1.5x . Li and Reynolds (1995) introduced the proportion of structural variation:

Relative proportion of structural variation -R SV
SH is a measure that describes the amount of landscape variability that is attributable to spatial (as opposed to random) effects. This measure can be obtained by subtracting the variogram nugget (c 0 ) from the sill (c) and then dividing by the sill. Since the sill represents the maximum (overall) variation, and the nugget represents pure random variation, subtracting both terms results in a measure of spatiallycorrelated variation. In order to provide a relative measure of structural variability, the relative change in SH can be measured as a function of increased field-of-view: where: The relative proportion of structural variation (R SV ) is a function of the range (a), nugget (c 0 ), sill (c), and the variogram estimator (γ E (h)). R SV is also dependent on separation distance (h), which allows for a full depiction of spatial patterns that may emerge at distances smaller than the range.

Results and discussion
4.1. Spatial representativeness of the AmeriFlux forest sites Fig. 5 shows a group of eight 18.0 km ETM+ subsets, with close-up images showing the circular (tower) footprints of each measurement site, as well as the 1.0 km 2 and 1.5 km 2 spatial boundaries. Results for R CV for each of these cases, and all other study locations, are summarized within Table 2. Fig. 6 also shows the corresponding variogram functions, for each measurement site presented in Fig. 5. This paper specifically focuses on the use of variogram functions over the native scales of MODIS observations. Thus, geostatistical attributes like R CV , are defined by information extracted from variogram functions in the range of 1.0 km 2 -1.5 km 2 . Nonetheless, the variogram responses using 2.0 km 2 subsets have been kept to further demonstrate the utility of this technique for analyses at other spatial scales. Fig. 5A-D shows a selection of University of California-Irvine (UCI) measurement sites that are a chronological series of field stations representative of secondary succession growth stages after large stand replacement fires. Black spruce trees undergo a slow growth process enabling the accurate determination of the chronosequence of stand age disturbance (Goulden et al., 2006). Additionally, boreal forests make up approximately 25% of forest ecosystems on Earth. With both of these in mind, the UCI sites provide an excellent location to study the impacts of surface albedo as a function of sequential wildfires. The UCI-1850 and UCI-1930 sites (Fig. 5A-B) are surrounded by closed canopies of tall black spruce, nearly 100% feather moss cover, and an open understory with a few alders, Labrador tea, and willow (McMillan & Goulden, 2008). When the size of the ETM+ footprint was varied from 1.0 to 1.5 km, R CV increased by +3.79% over UCI-1850 and decreased slightly by − 1.06% over UCI-1930. These are reasonably stable changes associated with smaller landscape components in the surrounding regions. A small road circumvents the northeastern section of UCI-1930 and just about crosses the 1.0 km 2 boundary. Since both the internal (CV 1 km ) and external (CV 1.5 km ) regions include this road, its effect on the overall variability is negligible, hence the low value for R CV . A similar degree of overall stationarity between regions was also found at the Morgan Monroe State Forest (MMSF) (Fig. 5F), where R CV increased slightly by 0.03%. This station is located over a secondary successional broadleaf forest within the maple-beech to oak hickory transition zone of the eastern deciduous forest, and is dominated by 60-80 year-old deciduous trees (e.g. sugar maple, tulip poplar, sassafras, white oak, and black oak) that have survived a selective logging period that ended over the past 10 years .  R CV . The same issue applies to UCI-1981, where a lake on the west side of the measurement site remains external to the imposed 1 km 2 site boundary.
At a height of 396.0 m, the WLEF-ChEAS Park Falls tall tower (Fig. 5G) has the largest instrument footprint of all the study sites. In the vicinity of the station, the forest is predominantly a mix comprising of 70% deciduous trees, including aspen, birch, maple, basswood, alder, and 30% conifers, including balsam fir, jack pine, black spruce, and white cedar. Upland areas are occupied by broadleaf deciduous and coniferous tree species, while the lowlands areas are occupied by tree species, shrubs, and a variety of grass and sedges. A further examination of the scale requirement index (R SE ) using five ETM+ scenes over the ChEAS site, with overpass periods during 17 March, 2001, 8 April, 2003, 26 May, 2003, 12 September, 2002, and 12 November, 2001, resulted in all R SE values b 0.1%. These results suggest that a direct assessment between ground-based estimates and MODIS retrievals over the ChEAS site will be suitable throughout all seasons. A similar trend was also found at the Tapajos National Forest (TNF) Km83-Logged station (Fig. 5E), where R SE b 0.1%. This station is categorized as a closed tropical forest with canopy emergents (canopy height = 35-40 m) from flat, upland terrain (Hernandez-Filho et al., 1993). The forest extended 5 km to the east, 8 km to the south, and 40 km to the north before reaching pasture (Miller et al., 2004), hence the low value for R SE .
The ground footprints of UCI-1850 and UCI-1930 are almost double the size of the younger stations (UCI-1964 andUCI-1981). It has been a common guideline among field programs to install tower albedometers at reasonable heights above the canopy (Loescher & Munger, 2006). Since the tallest trees at UCI-1964 andUCI-1981 extend to no more than 7.0 m, the stations were consequently instrumented at much lower heights (10-12 m). This arrangement constrains measurement of the surrounding landscape extending to a MODIS pixel. A spatial footprint analysis is even more challenging at UCI-1981, where the range of the measurement site varied from 437.48 to 703.36 m.
At the University of Michigan Biological Station (UMBS) (Fig. 5H), results for R ST show a very positive trend (70.94%). This field station lies on lake-border plains in northern Lower Michigan right in the transition zone between mixed hardwoods and boreal forests. In the vicinity of the tower (within a 1 km 2 radius), the forest canopy is predominantly deciduous broadleaf bigtooth aspen and trembling aspen, with significant numbers of red oak, beech, sugar maple, white pine, and hemlock. The understory is dominated by bracken fern and saplings of red maple, red oak, beech, and white pine (Curtis et al., 2002). A close assessment of the ETM+ retrievals over UMBS reveals that slope of the ascending limb at the 1.5 km 2 threshold is slightly sharper than at 1.0 km 2 . The primary source region that drives the strength of the variability can be attributed to early greening of a group of deciduous broadleaf trees, in addition to the conifer stands that are located at about 1.0 km southeast of the UMBS tower. For this particular ETM+ overpass (26 April, 2000), this sector appeared to have much darker albedos than the relatively dormant regions surrounding the station. An even more pronounced change in the strength of the variability can be observed at the 1.5 km 2 to 2.0 km 2 spatial threshold. Unlike the 1.0 km 2 vs. 1.5 km 2 case, these external differences of albedo are linked to the presence of lake water in the northeastern corner of the 2.0 km 2 region.
Results for R SV for each of the cases illustrated on Fig. 5 and all other study locations are summarized within Table 2. Large positive values for R SV are present over measurement sites where more patches, or clumps of trees in conjunction with bare areas, are observed when measuring the surrounding landscape extending to a The ground footprint, g, is a function of is the height of the field albedometer and its field-of-view (see Eq. (11)).
1.5 km 2 footprint. For example, by increasing the spatial footprint at the UCI-1850 site, the southwestern region merged with a green patch consisting of a combination of wetlands and/or an aspen grove. Both wetlands and aspen are deciduous covers which, in this particular region, reflects strongly in the NIR (0.7-5.0 µm). At UMBS, a cluster of early greened-up trees, just southeast of the station, also increased in size. These results suggest that R SV can detect spatial patterns that go against the overall trends captured by R CV . In using R SV as an indicator of spatial representativeness, this geostatistical attribute can also be used to examine the structural variability of a measurement site to resolve the relative magnitude of change patterns at different seasonal periods throughout the year. Fig. 7 shows ETM+ subsets over the Bartlett Experimental Forest for three time periods illustrating seasonal variations in greenness (A -26 August, 2000), to early senescence (B -27 September, 2000), and dormancy (C -20 October, 2000). The Bartlett Experimental Forest is located within the White Mountain National Forest in northcentral New Hampshire, USA. At low-to mid-elevations, vegetation is dominated by northern hardwoods (American beech, sugar maple, and yellow birch), with some red maple and paper birch also present. Fig. 7. Top-of-Atmosphere shortwave reflectance composites (ETM+ Bands 7-4-2), and the corresponding variogram plots, centered over Bartlett Experimental Forest for three time periods, illustrating conditions of greenness (A -26 August, 2000), early senescence (B -27 September, 2000), and dormancy (C -20 October, 2000). Trees and bushes are in shades of green (both light and dark tones) and bare areas are seen in light-lavender, magenta, and pale-pink. Results for each geostatistical attribute (i.e. R CV , A, R ST , and R SV ) are presented on Table 2. Conifers (eastern hemlock and eastern white pine) are occasionally found intermixed with the more abundant deciduous species but are generally confined to the lowest elevations and across the edges of the Albany Brook tributary, which crosses just west of the measurement site. In the vicinity of the tower, the forest is predominantly red maple, sugar maple, and American beech forest types (Jenkins et al., 2007).
The three ETM+ subsets in Fig. 7 illustrate how the effects of the fall foliage can significantly affect the structural variability of the Fig. 9. Time series and scatter plots of surface albedo at Local Solar Noon (LSN) for the Flagstaff-Managed (above) and Unmanaged (below) stations from late 2005 to late 2007. Setup is the same as Fig. 8. observed landscape. Both the internal and external footprints of the Bartlett measurement site are covered by bare areas, caused by clearcutting (shown by the light-pink patches). As the landscape enters into dormancy, these areas gradually merge with the fully dormant tree regions. Consequently, both R CV and R SV decreased during this progression as the landscape entered into dormancy. Fig. 10. Time series and scatter plots of surface albedo at Local Solar Noon (LSN) for the UCI-1850 (above) and UCI-1930 (below) stations from late 2001 to late 2004. Setup is the same as Fig. 8. In addition, the blue and red asterisks indicate snow-covered retrievals using the main and backup algorithms, respectively. However, the magnitude of the change in R SV (179.75%) was much larger than that of R CV , (11.19%). These results suggest that R SV was able to resolve the reduction in within-biome variability between the already existing bare areas and the fully dormant trees. A closer look at the variogram functions on Fig. 7 also show that the range resolved at the 1.0 km 2 boundary was more stable (160 ± 33 m), throughout Fig. 11. Time series and scatter plots of surface albedo at Local Solar Noon (LSN) for the UCI-1964 (above) and UCI-1981 (below) stations from late 2001 to late 2004. Setup is the same as Fig. 8. In addition, the blue and red asterisks indicate snow-covered retrievals using the main and backup algorithms, respectively. this 55-day period, than the 1.5 km 2 boundary (384 ± 115 m). Thus, the degree of "patchiness"as indicated by the presence of distinct patches between trees and bare areaswas much stronger in the immediate region surrounding the measurement site.

Comparison between ground-measured and MODIS-derived albedos
Figs. 8-11 show a number of plots that evaluate the daily performance of the MODIS V005 BRDF/albedo product (MCD43A3) against coincident ground data from the AmeriFlux sites. Each figure has been divided into pairs of measurement sites that are in close proximity to each other and share the same multi-year retrieval periods. For each time series plot, the dark points with error bars are daily retrievals of surface albedo. These values were derived by calculating the mean and standard deviation of all ground retrievals available within a 2-hour window centered at Local Solar Noon (LSN). Ground retrievals with large error bars are linked to intra-daily variations in downwelling irradiance, usually resulting from periods of increased cloudiness. Since the MODIS V005 BRDF/albedo product is routinely used as a clear-sky product, daily ground albedo records with estimates of cloud fraction N0.75 have been removed.
The blue circles and red triangles, on the time series plots, indicate snow-free blue-sky albedo retrievals at LSN derived from MCD43A3 using the main (full inversion) and backup algorithms, respectively. The blue and red asterisks, on Figs. 10 and 11, indicate snow-covered retrievals, as detected by the snow flag in the MODIS V005 albedo product embedded QA (MCD43A2), using the main and backup algorithms, respectively.
A series of scatter plots, located at the center of each figure, evaluate the daily blue-sky albedos from MCD43A3 against field measurements using four seasons: January-February-March (JFM) in blue, April-May-June (AMJ) in green, July-August-September (JAS) in red, and October-November-December (OND) in yellow. MODIS data from full inversion retrievals (in open shapes) were also separated from the "Full + Backup" inversion cases (in asterisks) to investigate the inversion quality differences when both methods are applied. The solid line is the one-to-one line and the dashed lines are ±0.02 and ±0.05 units. The statistical results for the absolute error (i.e. Ground mean-MODIS mean) and RMSE values calculated for the entire multi-year retrieval period are also displayed for each scenario (i.e. Full/Full + Backup). Table 3 provides summary statistics (RMSE values of Full/Full + Backup retrievals) for each seasonal class, plus a final assessment combining all data records from all measurement sites.
The daily retrievals from MCD43A3 agreed closely with the daily albedos over broadleaf forests, such as the Oak Ridge sites, Chestnut Ridge and Walker Branch (Fig. 8) as well as evergreen needleleaf forest sites, such as the Flagstaff stations (Fig. 9) and the UCI sites ( Figs. 10 and 11). At the Oak Ridge forest sites, surface albedo was stable throughout the seasons, with consistently small peaks reached throughout the dry "snow-free" winter periods of 2006-2007. MODIS albedos were usually within ± 16.5% of the ground-measured seasonal average albedo during the spring and summer periods, and within ±17.0% during the fall and winter. Both the Walker Branch and Chestnut Ridge towers are located on one of East Tennessee's long parallel ridges about 5 km apart from each other. The Flagstaff sites, located within Northern Arizona University's Centennial Forest, consist of a densely forested unmanaged site and a similar managed site that was partially restored to historical stand density conditions using mechanical thinning in September 2006 (Finkral & Evans, 2008;Sullivan et al., 2008). The Flagstaff-Unmanaged (or control) site represents a typical dense ponderosa pine stand in northern Arizona, and has not experienced forest management in over a century. The Flagstaff-Managed site is located about 6.7 km north from the control site. The former site was manually thinned in September 2006, reducing tree density by 67% (Finkral & Evans, 2008;Sullivan et al., 2008). MODIS albedos at the Flagstaff sites were usually within ±22.5% of the ground-measured seasonal average albedo during the spring and summer periods, and ±28.0% during the fall and winter. Daily retrievals over the Flagstaff-Managed site did not show a change in surface albedo associated with the thinning in September 2006. A comparison of the MODIS retrievals at the control site also denoted a consistent divergence from the ground measurements on both stations. The underling variations were most likely caused by periods of increased cloudiness during the Arizona Monsoon. MODIS albedos were often greater than ground-measured albedos during the cloudy and rainy monsoon season (July through September). Since the MODIS V005 BRDF/albedo product is routinely used as a clear-sky product, long-term periods of increased cloudiness usually cause an Table 3 RMSE (Full/Full + Mag) between field-measured albedos at LSN and MODIS V005 albedos (MCD43A3) sites, using the SWIR (0.3-5.0 µm) broadband, for four seasonal classes: January-February-March (JFM), April-May-June (AMJ), July-August-September (JAS), and October-November-December (OND).  Table 4 Correlation matrixes, for leaf-on (top) and leaf-off (bottom) conditions, comparing the geostatistical attributes of spatial representativeness (R CV , R SE , R ST , and R SV ) and the absolute RMSEs obtained for MODIS retrievals using all data records from each measurement site. overestimation of the albedo. This happens because only a few number of retrievals obtained prior to the rain season were used to describe underlying reflectance anisotropy, and thus the albedo, of subsequent periods. The MODIS algorithm performed much better during the same period in 2007, particularly over the Flagstaff-Managed Forest station, where a series of full and backup retrievals followed the sudden dips in surface-measured albedo. The MODIS algorithm also performed well throughout the seasons at the UCI sites. MODIS albedos were usually within ±16.5% of the ground-measured seasonal average albedo during the spring and summer periods and within ±21.2% of the ground-measured seasonal average albedo during the spring and summer periods. Over these older stations (UCI-1850 and UCI-1930), the percent tree cover was at least ∼15% higher than the younger UCI-1964 andUCI-1981 sites. These marked differences in tree density had a significant effect in the magnitude of the albedo throughout the three recorded periods of Processing and data flow diagram illustrating the use of the geostatistical attributes of spatial representativeness to identify, rate, and rank a group of measurement sites based on their ability to appropriately represent a large enough footprint to validate satellite-derived retrievals of surface albedo. A ranking procedure selects either a standard (ST) or RAW score, depending on variogram model accuracy. Under the standard scenario (i.e. Case 1) an "optimal" tower height can also be determined by minimizing for the ground footprint parameter, g, (Eq. (11)). snow cover (i.e. 2002-2004). For instance, MODIS albedos over the older UCI sites remained under 0.25 throughout snow periods, while the younger UCI sites reached a peak albedo of 0.45. Snow cover periods were also consistently shorter, by 20-30 days, at the older UCI sites. Thus, it appears that snow had less of an effect on albedo in the more mature forests. Table 4 shows two correlation matrixes, for leaf-on and leaf-off conditions, comparing the geostatistical attributes of spatial representativeness (R CV , R SE , R ST , and R SV ) and the absolute RMSEs obtained for MODIS retrievals using all data records from each measurement site. The correlation matrixes provide an indication of the relationship of the geostatistical attributes to one another. By further including the RMSEs, the geostatistical attributes that are playing a major role in explaining the differences in quality between the MODIS albedos and tower measurements can be identified. A first look at the results shows a lot of colinearity between R CV , R SE , and R ST . These correlations become even more pronounced during leaf-off conditions, where they also show a higher correspondence with the RMSE results. This appears to be the case when comparing the RMSE results between pairs of measurement sites. For instance, the scale requirement index (R SE ), was on average, 26.47% higher at the younger UCI sites during leaf-on conditions, and 36.67% higher during leaf-off conditions. Conversely, the RMSE levels the UCI sites were consistently lower (by ∼0.0133) at the older UCI-1850 and UCI-1930 sites. This degree of correspondence with the quality of MODIS albedo retrievals was also apparent over the Oak Ridge sites, where the geostatistical attributes at Chestnut Ridge were similar in magnitude for those at Walker Branch during leaf-on conditions (by ∼2.0%), but were much lower (by 20.32%) during leaf-off conditions. The RMSEs at Chestnut Ridge were also consistently lower throughout the year (by 0.0018 when compared to Walker Branch), and were even lower during leaf-off conditions (by ∼ 0.007). A closer look at the ETM+ subsets over the Oak Ridge forest sites ( Fig. 12A-D) shows a very similar spatial distribution of landscapes across Oak Ridge Forest. Despite the fact that the ground footprint over Chestnut Ridge was 50% larger (by ∼250 m) than the footprint over at Walker Branch, the average patch size, as defined by the range of the variogram at the station (see Fig. 13), was consistently higher by 82.35%. This resulted in similar magnitudes between the geostatistical attributes of both stations during the leaf-on periods, but more favorable results for Chestnut Ridge during leaf-off periods. A closer comparison of the remaining geostatistical attributes over the Flagstaff stations (Table 2) suggests that while the magnitudes of the scale requirement index (R SE ) were almost the same throughout the year, the differences in R CV , R ST , and R SV were as equally consistent in predicting which site (i.e. Flagstaff-Managed) was the most spatially representative of its underlying landscape.
The Oak Ridge sites show the expected seasonal pattern in a deciduous-dominated forest; with more vegetation during leaf-on conditions ( Fig. 12A-B) and less vegetation during leaf-off conditions ( Fig. 12C-D). In contrast, the Flagstaff sites do not show the expected pattern for the two time periods. In particular the leaf-on stage in June 2000 ( Fig. 12E-F) shows less vegetation and more bare ground than the leaf-off stage in November 1999 ( Fig. 12G-H). This unexpected pattern at Flagstaff is likely due to a pronounced drought in the winter and spring of 2000 that greatly limited green-up of herbaceous vegetation in the spring. Thus, what is labeled as "Leaf-On" in Fig. 12E-F for Flagstaff in June 2000 actually was a time of little herbaceous growth. These seasonal variations are confirmed in Adams and Kolb (2004), which documented a severe drought for 2000 based on time series analyses using the Palmer Drought Severity Index (PDSI) (Palmer, 1965).

Guidelines for use in validation and modeling activities
In developing the geostatistical attributes of spatial representativeness, one of the main objectives of this study is to develop the means to identify, rate, and rank a given number of measurement sites based on their ability to appropriately represent a large enough footprint to validate satellite-derived retrievals of surface albedo. As a secondary objective, this method should also provide ancillary information to improve site representativeness. However, the accuracy of these assessments depends on the quality of the spherical variogram model parameters and, in particular, their ability to provide a good fit to the variogram estimator. Only then will the confidence in the geostatistical attributes be strong enough to provide an accurate assessment. To address these conditions, two fitting scenarios with varying degrees of variogram model accuracy have been identified (Fig. 14). The differences between the "current" tower height and the "optimal" tower height were used to identify measurement sites that are appropriately representative of a large enough footprint to validate satellite-derived retrievals of albedo. A recommended increase in tower height over field stations that resulted in negative differences should improve site representativeness.
The first case focuses on the specific retrieval scenario when the spherical variogram model parameters, at both internal (i.e. 1.0 km 2 ) and external (i.e. 1.5 km 2 ) scales, are correctly fitting the variogram estimator. Using a combination of geostatistical attributes, with R SE as the primary marker of spatial representativeness and applying the other geostatistical measures as secondary weights, the ST score (or standard score) can be used to evaluate a group of measurement sites: Accurate estimates of R SE provide a true measure of the spatial extent of surface albedo with respect to both the ground instrument's fieldof-view and its surrounding landscape. Accordingly, a scenario that fits well, enables the use of R SE to determine the "optimal" tower height (i.e. the minimum height at which the requirements of spatial representativeness are satisfied) by minimizing for the ground footprint parameter, g, in Eq. (11). The second case focuses on the particular retrieval scenario when the spherical variogram model parameters, at both internal and external scales, are not providing a good fit to the variogram estimator. In such cases, the variograms are usually fitted against a linear function, which overestimates the values of the range, nugget variance, and sill. This occurs when the true value of the range parameter is larger than the half maximum distance of the measured subsets. The fact that the semivariance increases dramatically with range also implies the presence of an underlying trend. Thus, accurate estimates for ST score and the optimal tower height cannot be retrieved from these retrieval scenarios. However, because the relative coefficient of variation (R CV ) is only based on areal-mean estimates that are independent of variogram model parameters, a first-order (or RAW) score can be used as a replacement for ST score : Both ST score and RAW score are directly proportional to site representativeness. Thus, when selecting and ranking a number of measurement sites, the ones with the highest scores should also be the most suitable for use in direct "point-to-pixel" comparisons to the co-located satellite data. Since it is often difficult to determine the true shape of the variogram estimator (i.e. linear vs. spherical), a proper ranking of measurement sites (based on either the ST score or RAW score ) should be applied on a case-by-case basis. If for instance a measurement site returns a set of spherical variogram model parameters that provide a good fit to the variogram estimator, then the ST score should be used in the final rankings. Otherwise, if the parameters are fitted against a linear function (thus providing an improper fit to the variogram estimator), then the RAW score should be used instead. It is therefore recommended to calculate both ST score and RAW score and then select the most appropriate value by determining whether or not the spherical variogram model parameters are providing a good fit to the experimental data. In the ideal case, the ST score should be retrieved over all measurement sites. This allows for more accurate assessments of spatial representativeness, which in turn improves the estimates of optimal tower height. Preference should therefore be given to high resolution datasets that provide a sufficient number of samples both within the ground footprint of each measurement site and the surrounding landscape extending to the satellite pixel. Table 5 provides the final rankings of the AmeriFlux Forest Sites, based on ST score , during leaf-on and leaf-off conditions. The RAW score was also calculated, but was not used in determining the final rankings, since all the selected retrievals provided good spherical model fits. To assess the status of each measurement site, the difference between the "current" tower height and the "optimal" tower height was also determined. This measure can be used to identify measurement sites that are appropriately representative of a large enough footprint to be used for direct validation of satellite-derived retrievals of albedo. Conversely, for the field station heights that result in negative differences, the recommended increase in tower height can be used to improve their representativeness. In particular, a slight increase in the instrument's height at UCI-1964, from 12 to 17 m, should be sufficient to capture the true spatial extent over that measurement site. This would be a more difficult task over UCI-1981, where a minimum increase in tower height by 34 m is needed to fulfill the same requirements as UCI-1964. These results demonstrate the significance of this methodology; since it enables users to allocate their time, resources, and instrument efforts as efficiently as possible. Results also show that three of the topfive stations (i.e. Morgan Monroe State Forest, Santarem-Km83, and WLEF-ChEAS) were ranked as the most spatially representative during both leaf-on and leaf-off conditions. The same applies to the least representative sites (UMBS, UCI-1964, UCI-1981, andHowland Forest West Tower), which were ranked at the bottom five throughout the year.

Summary
A new validation framework for the estimation of the spatial representativeness can successfully indicate whether tower measurements are properly capturing the albedo retrievals over a large enough area to be suitable for use in direct "point-to-pixel" comparisons to the co-located satellite data. The approach combines knowledge of biophysical, spatial, and seasonal signatures of measurement sites and their surrounding landscape by utilizing multispectral high-resolution imagery as an intermediate step between ground and satellite retrievals. By first establishing variogram model parameters of the locale at multiple fields-of-view from high spatial resolution imagery, the geostatistical attributes of spatial representativeness can be generated. If a measurement site is spatially representative, then the overall variability between the internal components of the measurement site and the adjacent landscape should be similar in magnitude and relative coefficient of variation (R CV ) should approach zero. Similarly, the scale requirement index (R SE ) provides a measure of the spatial extent (or range) of surface albedo with respect to both the ground instrument's field-of-view and its surrounding landscape. A measure of the relative strength of the spatial correlation (R ST ) captures the maximum variance in the surrounding landscape and identifies the primary sources that drive the spatial variability over a given measurement site. Finally, the relative proportion of structural variation (R SV ) describes the amount of landscape variability that is attributable to spatial (as opposed to random) effects, and can further identify how seasonal changes, within and across the landscape, are affecting the structural variability of a given measurement site.
Results suggest that tower-based albedo measurements acquired over a broad range of forested landscapes are generally capturing the footprint of MODIS observations. However, measurement sites that were fitted with tower albedometers at very low heights (b20 m) above ground level (e.g. UCI-1964 andUCI-1981) or in close proximity to significantly different surface conditions (e.g. UMBS) will confront difficulties in capturing the albedo retrievals from a satellite footprint N1.0 km.
The geostatistical attributes of spatial representativeness have broad utility, but are particularly useful for: (1) identifying measurement sites that are appropriately representative of a large enough footprint to validate satellite-derived retrievals of albedo; and, if needed, recommend a series of instrumental upgrades to further improve their utility for evaluating satellite-derived retrievals of albedo; (2) producing a pixel-specific measure of product uncertainty both in terms of the quality of the algorithm inversions (e.g. given a limited number of satellite retrievals), and their ability to capture the underlying spatial and seasonal variability at local (b1.0 km 2 ) scales; and (3) identifying measurement sites that are appropriately representative of the broader regional ecosystems. For the latter, efforts are currently under way to further close the scaling gaps between the point-based and region-based assessments of spatial representativeness through the development of a similar set of upscaling routines to assess the spatial representativeness of other key terrestrial essential climate variables (e.g. surface temperature, snow cover, and leaf area index) at different spatial and seasonal scales. These approaches will ensure that the intrinsic spatial and seasonal features of flux towers are truly representative of their regional landscape, thus providing the necessary quality controls required by rigorous modeling efforts.