On the spatial heterogeneity of net ecosystem productivity in complex landscapes

Micrometeorological flux towers provide spatially integrated estimates of net ecosystem production (NEP) of carbon over areas ranging from several hectares to several square kilometers, but they do so at the expense of spatially explicit information within the footprint of the tower. This finer-scale information is crucial for understanding how physical and biological factors interact and give rise to towermeasured fluxes in complex landscapes. We present a simple approach for quantifying and evaluating the spatial heterogeneity of cumulative growing season NEP for complex landscapes. Our method is based on spatially distributed information about physical and biological landscape variables and knowledge of functional relationships between constituent fluxes and these variables. We present a case study from a complex landscape in the Rocky Mountains of Montana (US) to demonstrate that the spatial distribution of cumulative growing season NEP is rather large and bears the imprint of the topographic and vegetation variables that characterize this complex landscape. Net carbon sources and net carbon sinks were distributed across the landscape in manner predictable by the intersection of these landscape variables. We simulated year-to-year climate variability and found that some portions of the landscape were consistently either carbon sinks or carbon sources, but other portions transitioned between sink and source. Our findings reveal that this emergent behavior is a unique characteristic of complex landscapes derived from the interaction of topography and vegetation. These findings offer new insight for interpreting spatially integrated carbon fluxes measured over complex landscapes.


INTRODUCTION
Landscapes are heterogeneous assemblages of biotic and abiotic structures and processes (Turner 1989). These structures and processes interact with one another and with the external environment and give rise to emergent behavior and phenomena. Such emergent behavior can be considered a key feature of aggregate complexity (Manson 2001). Micrometeorological stations, flumes and other observational tools are useful for monitoring long-term physical and biological processes over large areas (hectares to square kilometers), but in doing so they blend the finescale heterogeneity of landscape processes and underlying structures into spatially integrated values that may or may not reveal the underlying complexity of the landscape (see recent discussion by Currie 2011). Spatially integrated measurements are valuable for assessing ecological processes at the landscape scale, but such measurements can mask fundamental information regarding underlying processes and structures of complex landscapes. For example, eddy covariance flux towers are powerful tools for measuring the temporal dynamics of net ecosystem production (NEP) of carbon and for making comparisons among ecosystems, regions and biomes because they integrate NEP over areas ranging from several hectares to several square kilometers (Law et al. 2002). However, flux towers are limited in their abilities to reveal the effects of spatial heterogeneity because they assume homogeneity of land-atmosphere (i.e., vertical) fluxes within the footprint of the tower. Ancillary data collected on or beneath flux towers are critical for interpreting temporal controls on NEP (Fig. 1). Additionally, measurements of variables such as soil CO 2 efflux can provide useful information for partitioning NEP into constituent vertical fluxes. Even with these complementary data, flux towers are generally unable to deconvolve the spatial heterogeneity of land-atmosphere carbon dynamics within their individual footprints without additional information about the spatial variability of the landscape (Goulden et al. 1996). This information may or may not be necessary for simple landscapes-those with uniform, homogeneous vegetation and essentially flat terrain, for example-where the inability of flux towers to resolve sub-footprint heterogeneity of NEP and other fluxes may be of little concern. However, this information is critical for many natural landscapes that are neither flat nor covered uniformly by vegetation (Baldocchi 2003).
Many natural landscapes are characterized by complex structures (e.g., vegetation, soil, topography, geology) that affect carbon cycle processes directly or indirectly and at different spatial and temporal scales. When flux towers integrate spatially over complex landscapes, they track average characteristics of the ecosystem or landscape continuously through time at subhourly resolution, but vital information about the v www.esajournals.org behavior of different parts of the landscape is lost. To illustrate, Riveros-Iregui and McGlynn (2009) measured soil CO 2 efflux at 62 locations within a tower-measured landscape for an entire growing season and found an order of magnitude variation in the cumulative flux among the locations. Tower-based estimates of cumulative ecosystem respiration were close to the median value of cumulative soil CO 2 efflux for the 62 locations, but the tower itself provided no information about this observed range of variability in the flux.
Spatially distributed models offer one solution to the problem of lost spatial information, but models typically use flux tower data for tuning or verification, seeking by design to reconstruct rather than leverage directly the rich temporal information that flux towers provide about the landscape as a whole (e.g., Aubinet et al. 2005, Detto et al. 2006, Emanuel et al. 2010. Finer-scale carbon flux observations within a flux tower footprint offer another solution to this problem, but care must be taken to collect data from points that span the full range of each variable that contributes to landscape complexity (Goulden et al. 1996). In this paper, we describe how the latter approach can be combined with high-resolution topography and vegetation data from Lidar to map tower-based NEP across a complex landscape. Given a tower-based estimate of NEP, the framework relies on landscape variables with known spatial distributions and known mechanistic links to the components of NEP to quantify the spatial variability of NEP within a complex landscape. In other words, we seek to characterize the heretofore-undefined frequency distribution of NEP for tower-measured landscapes.
Rather than a new mechanistic model, we are proposing a simple mathematical approach that incorporates very basic knowledge about terrestrial carbon cycling to explain the spatial distribution of NEP within a complex landscape. The framework identifies landscape positions most likely to function as net carbon sources or net carbon sinks across a range of biological responses to climate conditions and the magnitude of the source or sink. It also identifies positions on the landscape most likely to transition from net carbon sources to net carbon sinks across the same range of biological responses. In doing so, this approach allows researchers to ask questions including: How much spatial variability is integrated in a flux tower estimate of NEP? To what degree is a constituent flux responsible for this spatial variability? How does landscape complexity give rise to this spatial variability? Insight from this new method can also help determine where carbon sources or sinks are likely to be found on a specific landscape and how climate variability can drive changes in the distribution of sources and sinks and their magnitudes. Here we use the approach to gain a basic understanding of the nature and magnitude of NEP heterogeneity in complex landscapes.

PARTITIONING AND DISTRIBUTING NEP
We partitioned tower-based NEP into two constituent fluxes linked functionally to landscape variables that can be derived from remote sensing data (e.g., Lidar). Tower-based NEP can be expressed as where GPP is gross primary production (i.e., total photosynthesis of the ecosystem) and R E is ecosystem respiration (i.e., the sum of all respiration components of the ecosystem), both integrated over some time period (e.g., day, month, season, year). Eq. 1 is frequently applied to one-dimensional (vertical) fluxes, which are typically orders of magnitude larger than lateral fluxes at the ecosystem scale. Positive NEP indicates that the ecosystem is a net carbon sink during the integration period and negative NEP indicates that the ecosystem is a net carbon source. Both GPP and R E carry a positive sign in this equation; however, GPP and R E are carbon fluxes of opposite direction, typically with much larger magnitudes than NEP (Moncrieff et al. 1996). The components of R E can be partitioned into aboveground and belowground respiration allowing Eq. 1 to be rewritten as where R B is belowground respiration and R A is aboveground respiration. Functionally R B comprises root respiration and microbial decomposition in soils, whereas R A comprises respiration from aboveground plant tissue (Chapin et al. v www.esajournals.org 2002). Consequently, NEP can be partitioned among aboveground and belowground components by restating Eq. 2 as where GPP À R A is the aboveground component of the land-atmosphere carbon balance and R B is its belowground component. The aboveground land-atmosphere carbon flux (GPP À R A ) cannot be measured directly at the ecosystem level, but it can be expressed mathematically as since NPP (net primary production of carbon) is the difference between GPP and all autotrophic respiration (including aboveground and belowground respiration). For some plant species, both GPP and NPP increase monotonically with vegetation structural variables such as leaf area or tree height. For example in lodgepole pine, GPP is proportional to leaf area (Lindroth et al. 2008), and leaf area is correlated with tree height (Long and Smith 1988). The ratio of NPP to GPP is often constrained within a narrow range (Makela and Valentine 2001) and there is a strong correlation between these two fluxes for many forest ecosystems (DeLucia et al. 2007), and especially for lodgepole pine (Buckley and Roberts 2006), meaning that NPP can also be assumed proportional to leaf area and also to tree height. Since NPP and GPP both increase monotonically with leaf area and tree height, Eq. 4 suggests that GPP À R A can also be represented by a monotonically increasing function of leaf area or tree height for certain species. In other words, where Z is a variable describing vegetation structure (e.g., tree height), k is a constant of proportionality, and b accounts for the nonlinearity of this relationship. Given k and b, GPP À R A can be estimated from tree height data or other structural information that can be derived easily from airborne Lidar (e.g., Lefsky et al. 2002, Zimble et al. 2003. One notable exception to a general relationship between height and GPP À R A is that in older lodgepole pines (.50 years old) leaf area declines even though trees continue growing taller, causing corresponding decreases in GPP and NPP with loss of photosynthetic tissue (Ryan and Waring 1992). In such cases, Z would represent leaf area index instead of height. The belowground land-atmosphere carbon flux R B is measured directly using soil CO 2 efflux chambers. These chambers are commonly positioned at plots distributed among experimental treatments or across environmental gradients on the landscape (Scott-Denton et al. 2003, Riveros-Iregui et al. 2008. Recent work using chamber measurements and process-based modeling has established a strong relationship between cumulative growing season R B and topography (e.g., uphill drainage area and aspect) due in large part to the role of topography in determining the radiation balance, soil temperature fluctuations and the redistribution of soil water across the landscape (Sponseller and Fisher 2008, Riveros-Iregui and McGlynn 2009, Riveros-Iregui et al. 2011. For a subalpine lodgepole forest, Riveros-Iregui and McGlynn (2009) measured R B at 62 soil plots across topographic gradients of uphill drainage area and aspect and found that where U is uphill drainage area, A is aspect (southward, a or northward, b), and C 1 . . . C 4 are empirical coefficients unique to the field site and year. They determined C 1 . . . C 4 and mapped growing season R B at 10 m resolution for a small (3.9 km 2 ) watershed, the boundary of which approximates the seasonal flux footprint of a 40 m flux tower collocated at the site. This flux tower and experimental watershed have been described previously in the literature McGlynn 2009, Emanuel et al. 2010).
Having determined cumulative growing season R B empirically as a function of U and A for every landscape position (e.g., Riveros-Iregui and McGlynn, 2009), the mean value of R B for the landscape can be determined as well. When combined with cumulative growing season NEP from a flux tower integrating the same area, Eq. 3 yields the mean value of growing season GPP À R A for the entire landscape: where brackets denote spatial means. Thus hGPP À R A i is determined, but the spatial distribution of GPP À R A remains unknown.
v www.esajournals.org Finding the spatial distribution of GPP À R A and combining it with the empirically derived spatial distribution of R B is the last step in determining the spatial distribution of NEP for a complex landscape.
Since GPP À R A is assumed to be a function of a vegetation height described by Eq. 5 and the frequency distribution of Z is provided by Lidar, the derived distribution approach (e.g., Hebson and Wood 1982) provides a basis for estimating the values of k and b and thereby determine GPP À R A as a function of Z for every landscape position. The derived distribution approach states that an unknown frequency distribution Q(y) can be derived from a known frequency distribution P(x) if y is a monotonically increasing function of x. In this case, the unknown frequency distribution function Q(GPP À R A ) can be derived from the known distribution P(Z ) given knowledge of k and b. Coefficients k and b can be determined from Eq. 5 by substituting hGPP À R A i calculated from Eq. 7 and hZi calculated from P(Z ). Because hGPP À R A i depends on both tower-measured hNEPi and chamber-measured hR B i during a given growing season, valid combinations of variables k and b include a family of power law curves that describe a range of possible vegetation behaviors for a landscape during the same growing season. A priori knowledge of the shape of this functional relationship may be used to select appropriate values for k and b for a given species. Once these values have been selected, GPP À R A can be calculated for every landscape position where tree height is known. Thus, the derived distribution approach can be used to determine GPP À R A from Z even if the exact functional relationship is unknown.
Having determined GPP À R A and R B for every position on the landscape, Eq. 3 can now be used to determine NEP for every landscape position. The mean value of the newly derived NEP distribution is exactly equal to NEP measured from the flux tower. In the next section, we apply this framework to the forested watershed described by Riveros-Iregui and McGlynn (2009)

DISTRIBUTING NEP ACROSS A FORESTED WATERSHED
Upper Stringer Creek is a 3.9 km 2 watershed located at TCEF in the Rocky Mountains of westcentral Montana (US). The watershed ranges in elevation from 2090 m to 2430 m and is dominated by lodgepole pine (Pinus contorta) forest (Mincemoyer and Birdsall 2006). On average, TCEF receives more than 70% of its 880 mm of mean annual precipitation as snow (McCaughey 1996), which melts between April and June, yielding a relatively steady dry-down of soils throughout the growing season (Riveros-Iregui et al. 2007). Soils throughout the watershed are relatively shallow and uniform (Jencso et al. 2009), and each year they transition from relatively wet at the beginning of the growing season to relatively dry at the end, although not all landscape positions dry at the same rate or to the same extent (Emanuel et al. 2010). A 40 m flux tower has measured NEP by eddy covariance at Upper Stringer Creek since October 2005. In 2005, airborne Lidar measurements were collected over approximately 50 km 2 of the TCEF, including Upper Stringer Creek by the National Center for Airborne Laser Mapping (NCALM, Berkeley, CA). First return (canopy top) and last return (ground) elevation point data were interpolated to 1-m 2 digital elevation models (DEMs). We used the ground DEM to compute slope aspect (A) and uphill drainage area (U) at 10 m resolution following the methods of Seibert and McGlynn (2007), and we used the difference between the canopy top DEM and the ground DEM to compute tree height (Z ) at 5 m resolution following Emanuel et al. (2010). We subdivided A and U into 5 m 3 5 m grid cells to match Z. Fig. 2 shows U and Z for Upper Stringer Creek as well as the location of the flux tower.
In addition to Lidar-derived spatial variables, measurements of cumulative, tower-derived hNEPi and cumulative, chamber-derived R B are v www.esajournals.org available for the 2006 growing season (9 June-30 August, 83 days total). We combined these data using the approach described in the previous section, determining hR B i from Eq. 6 and combining it with hNEPi to determine hGPP À R A i. We assumed a linear relationship between lodgepole pine tree height and GPP À R A (i.e., b ¼ 1), and used known values of hGPP À R A i and hZi to solve Eq. 5 for k. Given k, b and the spatial distribution of Z, we used Eq. 5 calculate GPP À R A at 5 m resolution for the entire 3.9 km 2 Upper Stringer Creek watershed. Given GPP À R A and R B , we used Eq. 3 to determine NEP at 5 m resolution for the entire watershed (Fig. 3).
Having mapped cumulative growing season NEP for Upper Stringer Creek and determined its frequency distribution, three key results are evident. First, NEP is highly variable across the watershed. The frequency distribution of NEP is nonparametric and contains both net sources and net sinks of carbon that are very large relative to tower-derived hNEPi, which was approximately 270 g C m À2 for the 83-day growing season. The frequency distribution of NEP bears the imprints of both the distributions of R B and GPP À R A , illustrating in a spatially explicit fashion the nature of NEP as the difference of two opposing fluxes.
Second, the variability of NEP is highly organized and is a function of landscape variables U and Z (Fig. 4). Note that although A is incorporated into the analysis through Eq. 6, the effect of this variable on soil processes is secondary to U (Riveros-Iregui and McGlynn 2009) and secondary to Z (Emanuel et al. 2010) for Upper Stringer Creek. That a relationship exists between NEP and landscape variables is somewhat unsurprising since our approach assumes and is conditioned upon a priori knowledge of the relationship between landscape variables and carbon fluxes. However, nonintuitive patterns emerge from the intersection of the two dominant landscape variables, U and Z. These patterns can be used to predict regions of the landscape likely to be strong sinks or sources (i.e., combinations of U and Z that result in relatively large positive or negative numbers) or regions of the landscape likely to shift between source and sink (i.e., combinations of U and Z relatively close to 0) as factors such as climate Finally, this approach, though simple in its representation of carbon dynamics, suggests that estimates of landscape NEP and its constituent fluxes based on small-scale (e.g., 1 m 2 ) measurements can be greatly biased by site selection. That landscape heterogeneity affects small-scale measurements has been widely recognized across the flux tower community (Goulden et al. 1996, Scott-Denton et al. 2003, Aubinet et al. 2005; however, our approach provides a simple protocol for explaining how variability in smallscale measurement varies systematically across a landscape and can be predicted based on landscape structure. For example, every location along the dashed line of Fig. 4 is the same value as hNEPi measured by the flux tower. Positions above the line have higher-than-average NEP because of increased GPP À R A , decreased R B or both. Positions below the line have lower-thanaverage NEP because of decreased GPP À R A , increased R B or both. Note that the mean landscape position [hUi,hZi] produces NEP less than hNEPi, whereas the position defined by median U and median Z produces NEP very  v www.esajournals.org close to hNEPi. There are relatively few landscape positions where NEP equals hNEPi measured by the flux tower, and the assumption that hNEPi can be disaggregated into hR B i and hGPP À R A i uniformly across the landscape is clearly invalid as demonstrated by this example.

UNCERTAINTY IN THE PRODUCTIVITY-STRUCTURE RELATIONSHIP
The relationship between tree structure (e.g., height, leaf area) and relative productivity of a landscape represented in Eq. 5 can be subject to considerable uncertainty depending on tree age, tree density and species, and for a given landscape may change from year to year as a function of climate variability. We evaluated the effects of the productivity-structure relationship on the spatial heterogeneity of NEP by examining a range of values for b for the lodgepoledominated Upper Stringer Creek watershed as an example landscape.
Equation 5 can accommodate three types of functional behaviors: slow mathematical growth (b , 1), linear growth (b ¼ 1) and fast mathematical growth (b . 1). The example in the previous section was based on linear growth (i.e., productivity is directly proportional to tree height for lodgepole pines). Since the value of b determines the sensitivity of productivity (i.e., GPP À R A ) to vegetation structure (i.e., Z ) at different structural scales, we can use b as a coarse indicator of species effects, climate effects or interactions between species and climate (e.g., vegetation water stress). By varying b systematically for a single vegetation type (here, lodgepole pine), we can account for uncertainty inherent in this simple approach, and more importantly we can begin to understand how interannual variability in the productivity-structure relationship may affect the distribution of carbon sources and sinks in a complex landscape.
We simulated GPP À R A and NEP for b ¼ 0.5, b ¼ 1 and b ¼ 2 to represent the three behaviors of Eq. 5 and assess the effects of these behaviors on the spatial distribution of GPP À R A and NEP (Fig. 5). In each case cumulative growing season hNEPi was specified by the flux tower as it was in the previous section, meaning that for a given b there was only one possible value for k. As expected, in the slow mathematical growth case lodgepole productivity is most sensitive to tree height for smaller trees, and the sensitivity diminishes with height. This case produces the least spatial variability in NEP. In the fast mathematical growth case, fewer and taller trees are responsible for an increasing amount of GPP À R A . Thus, in the b ¼ 2 case, vegetation is responsible for more of the spatial variability of NEP than in the other two cases. Although hNEPi remains the same in each case, the distributions of carbon sinks and NEP vary as the relationship between cumulative growing season productivity and vegetation structure changes. The extreme values of spatially distributed growing season NEP (61000 g C m 2 ) are not unreasonable given the range of growing season soil CO 2 effluxes observed at the site by Riveros-Iregui and McGlynn (2009), and the range of annual GPP values modeled for different forest ecosystems, including another subalpine lodgepole forest (Wang et al. 2009). v www.esajournals.org Generalizing over a wider range of scaling behavior (b ¼ 0 to 4), the influence of vegetation structure on the heterogeneity of NEP becomes more apparent (Fig. 6). When b ¼ 0, every tree in the watershed sequesters the same amount of carbon during the growing season (i.e., productivity is uniform across the landscape). In this case, the spatial variability of NEP is determined entirely by R B and its primary dependence on U. Because the frequency distribution of U is skewed (Fig. 2), approximately 60% of the watershed has above average NEP (here ''average'' is defined by the range 270 6 20 g C m À2 ). As cumulative growing season productivity becomes more sensitive to tree size (i.e., as b increases in Fig. 6), fewer trees are required to sequester the carbon necessary to offset hR B i and yield tower-measured hNEPi. Hence, for a given tower-based measurement of hNEPi, the fraction of the watershed comprising net sinks decreases as the magnitude of each remaining sink grows larger. A real-life scenario that corresponds to this case could be an extremely dry year, in which taller trees with more root biomass and deeper root penetration with subsequently more access to water than shorter trees sequester a majority of the carbon and contribute disproportionally to hNEPi.
Simulating the effects of interannual climate variability by varying b also reveals three distinct regimes of hNEPi behavior (Fig. 7). First, some parts of the landscape remain carbon sinks regardless of how sensitive productivity is to tree height. For Upper Stringer Creek, these regions comprise relatively tall trees on relatively dry hillslopes. Modeling by Emanuel et al. (2010) suggests that for the relatively dry 2006 growing  Color gradient represents likelihood of carbon sinks for transitional landscape positions (i.e., those that move from net sinks to net sources as b increases).
v www.esajournals.org season, trees in these landscape positions experienced water stress, but not to the extent that stomatal conductance was drastically reduced. We assume this to indicate that tree height is still a valid first-order control on GPP À R A , even when feedbacks from drier topographic positions are considered. Second, other parts of the landscape always remain carbon sources. They comprise relatively short trees on lower hillslope positions where abundant soil moisture persists later into the growing season (subsidized by drainage from upslope areas), facilitating R B . Trees at these landscape positions have access to relatively abundant soil moisture, but their smaller structures limit their capacity to sequester carbon compared to larger trees with more photosynthetic tissue. Third, some areas of the landscape transition from source to sink (or vice versa) from year to year depending on the relationship between productivity and structure. These areas tend to be characterized by intermediate values of landscape variables U and Z, and the constituent fluxes GPP À R A and R B are of similar magnitude. Changing the sensitivity of hNEPi to one of the landscape variables through climate mediated processes can tip the balance and cause these landscape positions to function as either net sources or net sinks of carbon.
That these distinct regimes (i.e., perpetual sinks, perpetual sources and transitional areas) exist may be intuitive, but their pattern on the landscape is non-intuitive. The spatial pattern bears the unique imprint of our field site, Upper Stringer Creek, but it also provides general insight about the ability of the NEP distribution to change from year to year. Specifically, the existence of a relatively large transitional zone characterized by intermediate topography and vegetation suggests that in these areas year to year climate variability may be the primary factor determining whether a landscape position functions as a carbon source or a sink. In contrast, areas that are perpetual carbon sources or sinks owe their source/sink status primarily to the intersection of topographic and vegetation structure rather than to interannual climate variability, although climate variability may strongly influence the magnitude of the source or sink from year to year. In the absence of vegetation disturbance or successional change (which effectively alter the spatial distribution of Z ), the source/sink status of some areas (perpetual sources and sinks) is solely a function of spatial pattern (topography and vegetation at a given position), but the source/sink status of transitional areas is a function of both spatial pattern and temporal dynamics. For the landscape as a whole, hNEPi is thus defined by the coorganization of topographic and vegetation variables on the landscape and also by climatic influence on the relationship between spatial variables and constituent fluxes of hNEPi. The existence of different NEP regimes in the same complex landscape may help to explain why it is sometimes difficult to interpret the response of tower-integrated NEP to environmental variables (e.g., Fig. 1).
Landscapes lacking complex topography, vegetation or both may contain only one or two of the regimes because, spatially, they have fewer degrees of freedom than represented in Fig. 7. For example, uniform vegetation on flat terrain could be represented as a single point in the U À Z phase space. Depending on the specific functional responses of its constituent carbon fluxes, it could be situated in any of the three regimes and exhibit source, sink or transitional behavior. As long as the vegetation remains uniform, it will remain a single point in the phase space, but it will migrate in the Z dimension each year as vegetation structure changes. This example and the preceding examples only consider year-to-year variability of the GPP À R A distribution; additional work is needed to understand (1) how the distribution of R B changes from year to year in response to climate and (2) how vegetation disturbance and successional dynamics affect the climate-mediated dependence of GPP À R A on Z.
These results suggest that the degree of landscape complexity determines, in part, whether or not carbon sinks, carbon sources and transitional regimes exist or co-exist in a given landscape. Landscapes with complex topographic and vegetation structures (e.g., Upper Stringer Creek) provide sufficient heterogeneity to support all three regimes. For Upper Stringer Creek, we define landscape complexity using two main variables, uphill drainage area and tree height. Other landscapes may derive complexity from different topographic variables (e.g., large elevation gradients, high degree of slope variability), v www.esajournals.org different vegetation variables (e.g., differences in species or functional types, density of vegetation), or additional sources (e.g., soil physical or chemical characteristics, natural or human disturbance). Landscapes that lack recognizable complexity in topography, vegetation and/or other variables have fewer spatial degrees of freedom and may be less likely (if at all) to contain multiple regimes. Thus, interactions of heterogeneous landscape variables and the biophysical processes with which these variables are associated are responsible not only for the spatial distribution of NEP, but also for the emergence of multiple NEP regimes within a landscape. Upper Stringer Creek illustrates the potential for multiple regimes to coexist within a single complex landscape, but we posit that this emergent behavior is a general phenomenon and can be considered a defining characteristic of all complex landscapes.

CONCLUSION
We have presented a simple approach for distributing cumulative growing season NEP across a complex landscape. The approach requires knowledge of the spatially integrated NEP (e.g., from a flux tower), and it requires spatially distributed information about vegetation and topographic structure for the landscape (e.g., from Lidar). The approach also requires that a relationship between soil CO 2 efflux and topographic structure exists and is known. A relationship between productivity and vegetation structure is assumed to exist as a monotonically increasing function, but its exact form need not be known to determine the spatial distribution of productivity and NEP.
We demonstrated that the spatial distribution of cumulative growing season NEP bears the imprint of the spatial variables that influence biophysical processes within a complex landscape. We also demonstrated that considerable spatial heterogeneity may underlie tower-integrated measurements of NEP, but there is structure and organization to this heterogeneity that facilitates prediction of net carbon sources and sinks within the landscape. These results suggest that uniform partitioning of NEP into constituent fluxes is not appropriate across complex landscapes.
When we simulated year-to-year climate variability by varying the relationship between vegetation structure and productivity, we found that complex landscapes fostered three distinct behavioral regimes: some areas of the landscape were perpetual carbon sinks, some areas were perpetual carbon sources, and some areas transitioned between sink and source depending on the vegetation response to climate. This can be considered a defining characteristic of complex landscapes that emerges through the interaction of topography and vegetation. The behavior of spatially distributed NEP in simpler landscapes could be described as a special case of complex landscapes with fewer spatial degrees of freedom. Analyses of spatially integrated fluxes over complex landscapes should consider the spatial heterogeneity of fluxes and whether their distributions across the landscape could have an impact on interpretations and conclusions drawn from these data.

ACKNOWLEDGMENTS
This work was funded in part by the Department of Forestry and Environmental Resources at North Carolina State University and by NSF awards EAR-0404130, EAR-0837937 and DEB-0807272. Additional support was provided by the AGU Horton Research Grant awarded to Dr. Riveros-Iregui. Airborne Lidar was provided by the NSF-supported National Center for Airborne Laser Mapping (NCALM) at the University of CA, Berkeley. We are grateful to Robert E. Keane and Elaine Sutherland of the USDA Forest Service Rocky Mountain Research Station, Tenderfoot Creek Experimental Forest for site access and logistical support. Paolo D'Odorico provided valuable feedback on the mathematical approach and two anonymous reviewers provided insightful suggestions to strengthen the paper.