Biogeosciences Novel applications of carbon isotopes in atmospheric CO 2 : what can atmospheric measurements teach us about processes in the biosphere ?

Conventionally, measurements of carbon isotopes in atmospheric CO2 (δCO2) have been used to partition fluxes between terrestrial and ocean carbon pools. However, novel analytical approaches combined with an increase in the spatial extent and frequency of δCO2 measurements allow us to conduct a global analysis of δCO2 variability to infer the isotopic composition of source CO 2 to the atmosphere (δs). This global analysis yields coherent seasonal patterns of isotopic enrichment. Our results indicate that seasonal values of δs are more highly correlated with vapor pressure deficit (r = 0.404) than relative humidity ( r = 0.149). We then evaluate two widely used stomatal conductance models and determine that the Leuning Model, which is primarily driven by vapor pressure deficit is more effective globally at predicting δs (RMSE= 1.6 ‰) than the Ball-Woodrow-Berry model, which is driven by relative humidity (RMSE = 2.7 ‰). Thus stomatal conductance on a global scale may be more sensitive to changes in vapor pressure deficit than relative humidity. This approach highlights a new application of using δCO2 measurements to validate global models.


Introduction
The isotopic composition of atmospheric carbon dioxide (δ 13 CO 2 ) is a very powerful tool for inferring sources of CO 2 to the atmosphere, as well as processes affecting the global carbon cycle.Conventionally, δ 13 CO 2 has been used for partitioning net global CO 2 uptake between the land and ocean.
Correspondence to: A. P. Ballantyne (ashley.ballantyne@colorado.edu)This application of δ 13 CO 2 is based on the greater uptake of the lighter 12 C isotope by terrestrial plants, relative to the oceans.During photosynthesis plants discriminate against the heavier isotope 13 C in atmospheric CO 2 , yielding a globally flux-weighted estimated fractionation between the atmosphere and the terrestrial biosphere (ε al ) between −14.8 ‰ and −16.5 ‰ (Fung et al., 1997;Suits et al., 2005).In contrast, isotopic fractionation during air-sea gas exchange discriminates only slightly against 13 C (ε ao = −2.0‰) and thus is an order of magnitude less than isotopic discrimination by the terrestrial biosphere.Thus this differential isotopic fractionation has been used to partition carbon fluxes between the marine and terrestrial biosphere on global scales (Ciais et al., 1995;Battle et al., 2000).With the expansion of the global δ 13 CO 2 observation network (NOAA/ESRL), we are now capable of partitioning these fluxes at regional scales to assess the spatially heterogeneous response of the biosphere to climate variations (Bousquet et al., 2000).However, δ 13 CO 2 has largely been used in the inverse mode to solve for the partitioning of fluxes between the land and ocean.Here we outline a novel application for the use of δ 13 CO 2 in an entirely different mode of model testing.
The network of δ 13 CO 2 observations is continuously expanding and thus integrating much more detailed information about carbon cycle processes.The number of sites where regular flask measurements of atmospheric δ 13 CO 2 are being made has increased from 9 sites in 1990 to over 90 sites in 2010, with 4 sites currently using non-dispersive infrared analyzers to make continuous hourly in situ measurements of CO 2 .There has also been an increase in sampling from tall towers and aircraft, which allow us to better resolve vertical gradients of CO 2 and its isotopic composition in the atmosphere (Stephens et al., 2007).This A. P. Ballantyne et al.: Novel applications of carbon isotopes in atmospheric CO 2 increase in the frequency and density of δ 13 CO 2 observations has the potential to provide new insights into the global carbon cycle at finer spatial and temporal scales.Although we only include sites from the NOAA/ESRL global flask network (http://www.esrl.noaa.gov/gmd/ccgg/) in our analysis, these flask samples could be combined with other regional sampling networks or even eddy flux measurements where δ 13 CO 2 measurements are being made.
At the local scale δ 13 CO 2 observations are a useful diagnostic for inferring ecophysiological responses to environmental changes.Using the Keeling plot approach, researchers have been able to infer changes in the isotopic signature of respired CO 2 at the scale of individual foreststands.This approach has revealed insights into the stomatal response of forests to changes in their environment over broad regional scales (Pataki et al., 2003).Using the δ 13 CO 2 composition of recently respired CO 2 , researchers have been able to infer stomatal response to atmospheric water vapor, but that this isotopic signal may take weeks to be transmitted as respired CO 2 (Ekblad and Hogberg, 2001;Bowling et al., 2002).This approach has also been used at regional scales to determine that ε al is highly sensitive to total annual precipitation amount and to a lesser extent mean annual temperature (Pataki et al., 2003).Collectively, these observations from the Keeling plot approach confirm that stomatal conductance responds to water availability, as expected.However, these observations are restricted to the regional scale (i.e.primarily North America) and there are fundamental assumptions underlying the Keeling plot approach (see 2.0 Isotopic Theory) that make it challenging to apply to the global scale.
Recent advances in isotopic theory have also taken place that allow us to gain insights into processes in the terrestrial biosphere from atmospheric observations at the global scale.An extension of the Keeling plot approach has been used by Miller and Tans (2003) to account for changes in the background concentrations of atmospheric CO 2 and δ 13 CO 2 to infer source δ 13 CO 2 to the atmosphere.This approach makes use of the increased sampling effort of atmospheric δ 13 CO 2 .Based on this approach, it has been demonstrated that latitudinal gradients in ε al are not nearly as steep as previously predicted from models and that ε al of the North American biosphere has shown a steady decline (Miller et al., 2003).Furthermore, it has been demonstrated that seasonal climatologies of source δ 13 CO 2 inferred from the atmosphere show similar patterns to seasonal climatologies δ 13 C in cellulose of trees, suggesting that patterns of ε al can be inferred from atmospheric observations (Ballantyne et al., 2010).Although these insights derived from atmospheric observations provide independent empirical evidence of spatial and temporal changes in isotopic fractionation, they do not provide detailed information about processes occurring on seasonal scales at a global scale.
Here we propose a new application of δ 13 CO 2 by analyzing the wealth of δ 13 CO 2 observations using the analytical framework proposed by Miller and Tans (2003).Using this approach, we are able to derive seasonal distributions of source δ 13 CO 2 to the atmosphere for a range of sites included in the NOAA/ESRL flask network.Lastly, we hypothesize that these seasonal distributions of δ 13 CO 2 are driven by stomatal conductance and we use two empirical models of stomatal conductance to test this hypothesis.
2 Inferring CO 2 source using isotopes

Isotopic theory
Our ability to infer the isotopic signature source CO 2 to the atmosphere is based on the conservation of mass, such that the concentration of atmospheric CO 2 (c a ): is equal to the sum of the background CO 2 concentration (c bg ) and the CO 2 contribution from recent sources, positive or negative (c s ).Because δ 13 CO 2 is also effectively conserved in the atmosphere Eq. ( 1) can be expanded to include the product of CO 2 and its isotopic composition, where δ a represents the δ 13 CO 2 composition of atmospheric CO 2 , δ bg represents the δ 13 CO 2 composition of background CO 2 , and δ s represents the δ 13 CO 2 composition of source CO 2 .These equations were first combined by Keeling (1958), to derive the familiar Keeling plot: whereby the y-intercept corresponds with source δ 13 CO 2 .However, the most common application of this relationship is to infer the isotopic composition of respired CO 2 in relatively pristine forest environments, where δ bg is assumed to be constant.This assumption is usually satisfied by sampling at night in the absence of photosynthesis when the canopy air space is stratified.However, with the greater abundance of δ 13 CO 2 observations we need not assume that δ bg is constant, but rather we can rearrange Eqs. ( 1) and (2) to formulate the following expression: based on this approach we can specify a slowly varying background concentration of CO 2 (i.e.c bg ) and its isotopic composition (i.e.δ bg ) to solve for the slope-term δ s , which corresponds to the isotopic composition of source CO 2 to the atmosphere.δ s can then be used as a diagnostic tool for processes regulating the transfer of carbon between the biosphere and the atmosphere.In theory, we can use δ s to approximate the isotopic signature of the terrestrial biosphere (δ l ), such that: where ε al is primarily regulated by stomatal conductance on shorter time-scales, according to (Farquhar et al., 1989): where a is a constant representing fractionation during diffusion (−4.4 ‰) and b is a constant representing fractionation during carboxylation by Rubisco (−27 ‰).
Although it has been shown that the isotopic signature of respired CO 2 is dependent upon the isotopic composition of the carbohydrate consumed during decarboxylation on diurnal timescales (Tcherkez et al., 2003), on seasonal timescales, at the forest stand level, the isotopic composition of respired CO 2 is thought to co-vary with the isotopic composition of recently assimilated carbon in the form of sucrose (Scartazza et al., 2004).Moreover, respired CO 2 is dominated by recently assimilated carbon at weekly to monthly timescales (Högberg et al., 2001).There is no evidence of isotopic fractionation of recently assimilated carbon during autotrophic respiration at the cellular level (Lin and Ehleringer, 1997).However, diel variations in the isotopic composition of plant respired CO 2 are often observed and cannot be explained solely by the isotopic composition of the photosynthetic substrate respired.Thus other mechanisms, such as temperature and light availability, are likely factors leading to post-photosynthetic isotopic fractionation of respired CO 2 (Werner and Gessler, 2011).Our previous analysis using a global model revealed that the coherent patterns in δ 13 C in cellulose of the biosphere and δ s inferred from the atmosphere during the growing season were driven primarily by stomatal conductance (Ballantyne et al., 2010).Therefore our estimates of δ s should be a suitable proxy for the isotopic signature of recently assimilated CO 2 (i.e.δ l ).

Analytical approach
Our approach for inferring the isotopic signature of source CO 2 to the atmosphere is contingent upon selecting a suitable background reference curve for calculating residuals.It has been demonstrated previously that the free-troposphere (2000 to 5000 m a.s.l.) represents the best background reference as it introduces the fewest artifacts when inferring a seasonal cycle in δ s (Ballantyne et al., 2010).Ideally, we would use free troposphere observations immediately above each atmospheric sampling site as our background reference; however, this is only possible for a limited number of sites.Instead Niwot Ridge, CO (NWR, 3437 m a.s.l.) has been identified as a suitable mid-continental background site for observations in North America.
To demonstrate how this analysis is performed, let us consider the isotopic variability of CO 2 observed at Wendover, UT (UTA).We see that there is considerably more variability in atmospheric CO 2 and δ 13 CO 2 observations at UTA than the more attenuated tropospheric reference curve at NWR (Fig. 1a).Atmospheric CO 2 levels at UTA are greater during winter months and lesser during summer months than those at NWR (Fig. 1a).Although the curves differ in seasonal amplitude, they tend to be in phase on seasonal timescales.Atmospheric δ 13 CO 2 values are a reflection of the CO 2 curves, with values more depleted at UTA than NWR during winter months and values more enriched at UTA than NWR during summer months (Fig. 1b).Essentially, surface observations show an amplification of the seasonal cycle in both CO 2 concentration and isotopic composition relative to the free troposphere.
If we take the residuals between CO 2 in the boundary layer and CO 2 in the free troposphere and we plot them against the residuals between the product of CO 2 and δ 13 CO 2 in the boundary layer and the product in the free troposphere, we can then determine the slope (Fig. 1c).This slope value corresponds to δ s and it is evident that δ s varies according to season.If we solve for δ s for each month, using our moving window approach, a distinct seasonal pattern emerges whereby winter months are characterized by more depleted δ s values and summer months are characterized by more enriched δ s values (Fig. 1d).The magnitude of the seasonal pattern (i.e. the signal) greatly exceeds the uncertainty in any given month (i.e. the noise), suggesting a robust seasonal isotopic signal that be an effective tool for detecting the source of fluxes from the biosphere to the atmosphere.

Site selection
For our analysis we only included terrestrial sites from the NOAA/ESRL network of atmospheric sampling sites with at least 5 yr of data.These terrestrial sites are generally located away from large fossil fuel point-sources of emissions, but some downwind sites may be affected by cumulative emissions.These criteria yielded 18 sites, 5 of which were tall tower sites, all located within the Northern Hemisphere (Table 1).Tall tower sites were sampled daily and surface sites were sampled weekly.On average these sites had 13 yr of data, which is more than sufficient for determining seasonal distributions of δ s .Although a rigorous footprint analysis of all these sites has not been performed, the footprint of the tall tower measurements tends to be much larger (∼10 4 km 2 ) than the footprint for surface measurements (<1 km 2 ) (Helliker et al., 2004).There are several potential time series to specify as background isotopic concentrations for inferring δ s , but previous analysis has indicated that observations at Niwot Ridge, CO, USA (NWR) are statistically indistinguishable from independent observations made from the free troposphere (2000 to 5000 m a.s.l.) making it an excellent background reference for the Northern Hemisphere (Ballantyne et al., 2010).

Data analysis
For each site a seasonal distribution of δ s was calculated.A moving three month window was used for calculating δ s values, such that January δ s values were calculated from observations made during December, January, and February (DJF) and February δ s values were calculated from JFM observations, etc.This approach has been shown to yield the most robust seasonal patterns with the smallest standard error estimates for any given month (Ballantyne et al., 2010).Values of δ s were calculated by first subtracting the background reference curve from time series at each site and then using a linear least-squares regression that incorporates error terms in both the x and y axes (Miller and Tans, 2003) to calculate the slope (δ s in Eq. 4).Although such a regression analysis with errors in both the x ,and y terms may lead to values of δ s that are biased (Zobitz et al., 2006), this is at ranges of CO 2 concentration below 10 ppm, which are much lower than those in our 3 month moving window of atmospheric observations.All atmospheric observations were screened for anomalous values that might contribute disproportionately to our regression.All atmospheric CO 2 and δ 13 CO 2 observations that exceeded 2σ from the mean of our 3 month moving window were excluded from our analysis.This screening of data only removed between 0 and 2 % of observations across all sites, but greatly reduced anomalous values that may have caused an over-amplification of the seasonal cycle in δ s .

Evaluation of models
To test models designed to simulate the isotopic fractionation occurring during stomatal conductance, we used the simple biosphere model SiB biosphere model (Sellers et al., 1996).The model was driven by National Centers for Environmental Prediction Reanalysis Data (Kanamitsu et al., 2002) interpolated to the model timestep for the years 1983-2006.Maps of plant functional types were derived from remote sensing products (DeFries and Townshend, 1994).Mean monthly values of assimilation-weighted leaf surface temperatures (T ) and relative humidity (RH) were calculated for each grid cell (1 • × 1 • ) encompassing a network site using the most recent version SiB3 (Baker et al., 2010).For our regression analysis we only included months when net exchange between the atmosphere and biosphere was negative, resulting in an annual cycle that was truncated to the growing season.This was done to isolate the isotopic signal attributable to carbon that had recently been assimilated by the biosphere (Miller et al., 2003).Leaf T and RH were then used to calculate saturation vapor pressure and ultimately Table 1.Seasonal correlations between δ s and metrics of atmospheric water vapor for terrestrial sites from the NOAA/ESRL network included in this analysis.For each site the corresponding station code from the NOAA/ESRL flask network and the biome according to SiB3 are reported in addition to the duration and growing season length of the dataset included in the analysis.Correlation coefficients (r) are reported between seasonal distributions of δ s and vapor pressure deficit (VPD) and relatively humidity (RH), respectively (see Fig. 2).The sign of the correlation coefficient indicates the sign of the relationship.Significant correlation coefficients (p < .05)appear in italics and highly significant correlation coefficients (p < .01)appear in bold italics.Root mean squared error (RMSE) are reported for growth season predictions of δ s from the Leuning and BWB stomatal conductance models compared with δ s values inferred from the atmosphere.Mean RMSE values across all sites are also reported.vapor pressure deficit (VPD).The assimilation-weighted values of RH and VPD were then used as the primary variables driving 2 commonly used stomatal conductance models-the Ball-Woodrow-Berry (BWB) Model (Ball, 1988): and the Leuning Model (Leuning, 1995): The effects of these models on isotopic fractionation by the biosphere were evaluated based on the framework outlined by Katul et al. (2000), whereby the equation for assimilation (Farquhar and Sharkey, 1982): was substituted into the assimilation term (A) for both models which were then solved for ratio of intercellular to atmospheric CO 2 (c i /c a ).All model parameters and the values used for the various biomes considered in this study are reported in Table 2.

Globally coherent patterns
The distribution of δ s inferred from a network of atmospheric sampling sites (Fig. 2a) reveals globally coherent seasonal patterns.The array of Northern Hemisphere sites shows a consistent pattern of enriched δ s values during summer months and more depleted values during winter months (Fig. 2b).Most of the sites included in our analysis showed maximum δ s values between July and August (see Supplement Fig. S1); however, lower latitude sites appear to reach maximum δ s values (July-August) prior to higher latitude sites (August-September).The seasonal amplitudes in δ s values vary between 2 and 7.6 ‰, with arid low latitude desert sites, such as KZD, Kazakhstan and ASK, Algeria, experiencing much greater amplitudes in seasonal variations of δ s , than mid-latitude mixed forest sites, such as BAL near the Baltic Sea and TAP in Korea (Table 1).Values of δ s inferred at most of our sites were more highly correlated with assimilation-weighted values of VPD than RH calculated at the leaf's surface (Table 1).Correlation coefficients for all sites ranged from 0.03 to 0.99, but The intercept term in the original BWB model "b" (Eq.7) can be neglected when solving for isotopic discrimination (Katul et al., 2000).b Similarly the intercept term in the Leuning model "g o " can be neglected when solving for isotopic discrimination (Katul et al., 2000) c * represents the CO 2 compensation point and varies as function of leaf surface temperature (T • C), according to Brooks and Farquhar (1985).mid-latitude and low-latitude sites with greater seasonal amplitude in δ s tended to have higher correlations with both VPD and RH.Of the 19 sites included in our analysis 16 showed significant correlations with vapor pressure deficit and only 6 showed significant correlations with relative humidity.
Although there was generally a strong relationship between seasonal δ s values and metrics of atmospheric water  vapor across sites, some sites showed a stronger response to VPD and some sites showed a stronger response to RH.Of the 16 sites that showed a significant response to VPD, 10 were more highly correlated with the natural log of VPD and of the 6 sites that showed a significant response to RH, only 2 were more highly correlated with the natural log of RH.These results indicate a non-linear response of stomatal conductance to atmospheric water vapor, especially at the sites that were more responsive to VPD.
Although some sites do not correlate well with metrics of atmospheric water vapor, if we look at the global dataset of seasonal δ s values for all sites there is a much stronger correlation with VPD than RH (Fig. 3).We found the optimal global fit to be between δ s and the natural log of VPD (r = −0.404,p-value = 2.2 × 10 −6 , DF = 152), indicating an incrementally smaller amount of isotopic fractionation at higher VPD values.There was only a slight increase in this correlation when individual non-significant sites (see PAL, OXK, and TAP in Table 1) were removed from the analysis (r = −0.408,p-value = 2.2 × 10 −6 , DF = 148).Regional differences between sites also emerged from this analysis, such that high-latitude sites tended to have a much greater response in δ s over a smaller range of VPD values, whereas δ s was not nearly as sensitive to changes in VPD values at lower latitude desert sites (Fig. 3a).Although the global relationship between δ s and RH was also significant, the relationship was not nearly as strong (r = 0.149, p-value = 0.023, DF = 152) and the optimal relationship was with a natural log transformation of RH.These global patterns suggest that if the isotopic signature of source CO 2 to the atmosphere is indeed due to stomatal conductance, then the physical mechanism responsible for these patterns of isotopic discrimination is probably VPD on a global scale.

Evaluation of models
Values of δ s inferred from atmospheric measurements can also be used to test biosphere models created to simulate the exchange of mass at the biosphere-atmosphere interface.Because stomatal conductance is the primary mechanism causing isotopic fractionation during photosynthesis (Farquhar et al., 1989) and there is no net fractionation associated with autotrophic respiration (Lin and Ehleringer, 1997), we can then use values of δ s to gain insight about factors controlling stomatal conductance.In most cases, the Leuning model driven by changes in VPD performed better at predicting δ s than the Ball-Berry model driven by changes in RH (Table 1).At almost all sites, root mean squared error (RMSE) estimates were less for the Leuning model than the BWB model.The exceptions were UTA, and WKT.Globally, RMSE values were significantly lower for the Leuning model, 1.6 ‰, than for the BWB model, 2.7 ‰ (two tailed ttest, p-value = 0.00106, DF = 27), suggesting that VPD may be more important in governing stomatal conductance than RH at global scales.
Although the Leuning model tends to outperform the BWB model, there are instances where the Leuning model deviates from observed δ s values (Fig. 3).For instance, the Leuning model tends to be fairly accurate in its predictions of δ s towards the mean of the distribution in observed δ s values; however, it appears to be biased towards more enriched values at lower and higher observed δ s values (Fig. 3c).In contrast, the BWB model tends to over-predict depleted δ s values from low-latitude sites and under-predict enriched δ s values from high-latitude sites (Fig. 3d).

Analytical approach
The approach that we have presented here, relying solely on atmospheric observations is effective for extracting seasonal information regarding biosphere-atmosphere interactions on a global scale.By specifying background concentrations of CO 2 and δ 13 CO 2 we are able to generate seasonally coherent patterns of δ s for most Northern Hemisphere sites included in this analysis.Here we have refined the original approach presented by Miller and Tans (2003) by resolving seasonal patterns on a global scale.Previous analyses of atmospheric observations and model simulations have identified the free troposphere as the most effective background reference curve for inferring δ s values at regional scales and that the high elevation mid-continental site at NWR is an effective background reference, at least for North America (Ballantyne et al., 2010).Here we have extended this approach to a wider array of Northern Hemisphere sites and our analysis reveals coherent seasonal cycles of δ s among these sites.
Although our more extensive analysis has revealed coherent seasonal patterns, there are subtle differences in these seasonal patterns that may be artifacts introduced during our analysis.For instance there appears to be a characteristic "stair-step" pattern in δ s at higher northern latitudes.This is evidenced by the dark blue curves in Fig. 2b representing seasonal patterns in δ s for Barrow, Alaska, USA (BRW); Baltic Sea, Poland (BAL), and Pallas-Sammaltunturi, Finland (PAL).Unlike lower-latitude sites that show a pronounced seasonal cycle of more enriched δ s values during summer months and more depleted values during winter months, these high-latitude sites show depleted δ s values during winter months, more enriched values during summer months, but then δ s values remain high into the fall months (see Supplement Fig. S1).This "stair-step" pattern may also be contributing to the reduced correlation coefficients between δ s and VPD and RH, as well as the less optimal fit of stomatal conductance models at higher latitudes (Table 1).Indeed, we would expect for lower correlations at higher latitudes as the length of the growing season considered in our analysis becomes greatly reduced at higher-latitudes, effectively reducing our sample size.It is difficult to ascertain whether this is a circumarctic phenomenom based on our analysis from only three sites, but it does deserve further investigation.The anomalous seasonal pattern observed at these sites may also be due to the mid-latitude background site (i.e.NWR) selected for our analysis.There may be lags introduced into the analysis as a result of the transport time of mid-latitude background air by Ferrel cells advecting air poleward.However, with increased aircraft sampling, background reference curves from higher latitudes are becoming increasingly available and more highly resolved.Some low-latitude sites also showed deviations from the seasonal patterns characteristic of mid-latitude sites.www.biogeosciences.net/8/3093/2011/Biogeosciences, 8, 3093-3106, 2011 Generally, lower-latitude sites had larger seasonal amplitudes in δ s than higher-latitude sites.However, there is a clear anomaly to the seasonal pattern in Assekrem, Algeria (ASK) where the primary peak in δ s values is observed in July, but a secondary peak in δ s is observed in November (see Supplement Fig. S1).Although this anamolous pattern may be due to our specification of the very distant NWR background reference curve for this site, it may also be due to this sites proximity to the equator.In our original analysis we included two equatorial sites (Bukit Kototaband, Indonesia: BKT and Mt.Kenya: MKN) and two southern hemisphere sites (Cape Grim, Tasmania: CGO and Gobadeb, Namibia: NMB).However, seasonal patterns were not as clear and errors were much greater at these sites, leading to their ultimate removal from our analysis.Although the equivocal results in the Southern Hemisphere are most likely due to our Northern Hemisphere background reference curve, this approach may yield ambiguous results at equatorial sites due to the inter-hemispheric transport of air masses.According to SiB3, ASK is classified as a "low-latitude C4 dessert" of Northern Africa.Recently revised global maps of C4 vegetation, however, show most C4 vegetation to be restricted to the savannahs of sub-Saharan Africa (Still et al., 2003).Furthermore, there is very little vegetation in this mountainous interior region of Algeria, based on satellite-derived estimates of vegetation cover (Foley et al., 2003), suggesting that the seasonal cycle of δ s at ASK may be an admixture of isotopic signals from Mediterranean biomes to the north and sub-Saharan biomes to the south.Therefore it is possible that tropical sites, especially those in Africa, may be subjected to both C3 and C4 vegetation signals that differ in their seasonal timing and δ s signatures.The choice of background reference curves is not trivial.In fact, our previous analysis (Ballantyne et al., 2010) comparing model simulations with observations indicated that using a reference curve from the "free-troposphere" above 3000 m a.s.l. was the optimal reference curve for inferring isotopic source signature of fluxes from the biosphere to the atmosphere.Unfortunately, such high elevation "free troposphere" reference sites are not available for comparison with all the surface sites on all continents.The alternative would be to select the appropriate marine boundary layer reference curve for each site based on latitude.Although this would address latitudinal differences in the isotopic signature of background air-masses (i.e.δ bg ),it would introduce other artifacts, such as possible seasonal patterns in isotopic signatures from air-sea gas exchange.We acknowledge that there are assumptions in using NWR as the reference curve for all sites in the Northern Hemisphere; however, using regional reference curves from surface sites would also require assumptions.Our previous analysis indicated that the choice of background reference curves did not affect the seasonal amplitude of δ s so much as the timing of peak δ s values (Ballantyne et al., 2010).Therefore our choice of NWR as the background reference curve for all Northen Hemisphere sites may have introduced temporal biases into our analysis; however, these potential temporal biases do not appear substantial enough to have impacted the significance of our results (Table 1).Moreover, these temporal biases in δ s would not differentially affect our seasonal correlations with leaf level estimates of VPD and RH.Thus as more free troposphere reference curves become available for other regions of the world; fewer assumptions will be necessary in specifying background reference curves and ultimately calculating δ s at regional scales.

Globally coherent patterns
Our analysis of a network of sites in the Northern Hemisphere has revealed seasonally coherent patterns in the isotopic composition of source CO 2 to the atmosphere.These broad global patterns are consistent with previous analyses done on individual sites.For example, Bakwin et al. (1998), employing the Keeling plot approach at ITN and HUN identified a slight enrichment of δ s during summer months compared to winter months.Bakwin et al. attributed the changes in δ s at these two sites to isotopic sources that were dominated by fossil fuel sources during winter months and dominated by terrestrial uptake during the summer months.Our results are consistent at ITN and HUN; however, we are able to take advantage of additional data to resolve a clear seasonal cycle in δ s at these two sites.Although there are numerous confounding factors that may influence δ s on seasonal time scales, the globally coherent seasonal cycle in δ s reported here indicates a common underlying physical mechanism rather than disparate local factors influencing δ s .
At local to regional scales fossil fuel emissions may have a large impact on δ s values inferred from atmospheric observations, especially near large urban areas.In fact, two of the sites that show no correlation between δ s and VPD or RH are located in heavily industrial regions of South Korea (TAP) and Eastern Germany (OXK).The strong isotopic source from fossil fuels has been clearly demonstrated in urban environments that experience strong inversions during winter months.For example in Salt Lake City, midday δ 13 CO 2 values may become depleted during winter months by as much as 5 ‰, which has been attributed primarily to increased emissions from natural gas used for home heating (Pataki et al., 2006).Although fossil fuel emissions should be taken into account when evaluating δ s values at individual sites, fossil fuel consumption and fossil fuel type show broad spatial and temporal variability.For instance, the Eastern US relies primarily on coal and the Western US relies primarily on natural gas for power generation (Pétron et al., 2008) and emissions from these sources have very different isotopic signatures.Furthermore, emissions tend to be higher during winter months in northern US states due to home heating, whereas emissions tend to be higher during summer months in southern US states due to home cooling (Gregg et al., 2010).Despite this regional variability in the timing and type of fossil fuel consumption our network of sites shows a clear and consistent seasonal cycle in δ s , especially across the continental US.Furthermore emissions from respiration (∼120 GtC yr −1 ) greatly exceed emissions from fossil fuels (∼8 GtC yr −1 ); thus a very strong seasonal cycle in fossil fuel emissions or an extremely variable isotopic signature of fossil fuel emissions would have to be invoked to account for the 6 to 8 ‰ range in the observed globally coherent δ s values.Therefore fossil fuels may help to explain the variability in seasonal δ s between sites, but fossil fuel emissions cannot explain the global patterns of δ s revealed by our analysis.
Changes in climate may also affect the seasonal cycle in δ s .In fact, a covariance between isotopic discrimination by the terrestrial biosphere and precipitation amount as mediated by El Niño events has been previously identified (Randerson et al., 2001).Although climate variability probably affects stomatal conductance and thus isotopic discrimination on inter-annual scales, this contributes to the error in our estimates of the seasonal climatology of δ s .The fact that a strong seasonal signal in δ s emerges despite the error associated with inter-annual climate variability suggests that temporal variability in δ s is actually dominated by the seasonal cycle.However, our seasonal climatologies of δ s could be used to estimate anomalies in any given year due to climate variability.For instance, if prolonged drought results in diminished stomatal conductance this could be diagnosed as an enriched departure from our seasonal climatology of δ s due to decreased isotopic discrimination.Long term variability in isotopic discrimination may also arise from changes in atmospheric CO 2 concentration and/or climate change.Model simulations of the global carbon isotopic budget indicate a decrease in isotopic discrimination by the terrestrial biosphere of approximately 0.4 ‰, primarily due to an increase in water stress (Scholze et al., 2003).Such long term trends in isotopic discrimination should probably be considered, especially at some of the sites included in our analysis with more data (∼20 yr).However, this potential change in isotopic discrimination of 0.4 ‰ yr −1 , is dwarfed by the amplitude of the seasonal signal for sites included in our analysis (between 4 and 8 ‰).
Changes in land use may also impact values of δ s over time.It has been suggested that the intensification of agriculture and the proliferation of C4 crops has offset a fraction of CO 2 emissions (Burney et al., 2010).This has been demonstrated at the regional scale in the Amazon where the widespread conversion of tropical C3 forest to C4 pasture should lead to a decrease in isotopic discrimination and such a reduction in regional isotopic discrimination could lead to the spurious conclustion that terrestrial uptake in the tropics has decreased (Townsend et al., 2002).Unfortunately, there are very few tropical sites with continuous observations of δ 13 CO 2 making it impossible to extend this analytical approach into the tropics.However, there has been widespread land-use change at temperate latitudes as well, especially the proliferation of C4 corn to meet biofuel demands.Recent analyses of δ 13 CO 2 and CO 2 observations made from a tall tower in the Midwestern US have identified a strong isotopic signal from increased corn cultivation (Griffis et al., 2010).Griffis et al. (2010) suggest that increased corn cultivation for biofuels lead to an apparent decrease in isotopic discrimination from ∼15 ‰ to ∼12 ‰.Although most of this isotopic effect due to corn production was observed during the summer months the isotopic effect and its impact on δ 13 CO 2 seem to extend into the fall as well.These independent observations from Griffis et al. (2010) are a mere 300 km. from our LEF site and may actually help to explain some of the seasonal anomalies observed at LEF.Although LEF exhibits the same seasonal pattern in δ s characteristic of most our northern hemisphere sites, there is a secondary peak in δ s that occurs in the fall and may be due to decreased isotopic discrimination as a result of increased corn cultivation (Fig. 2).Furthermore, the leaf temperature estimates from SiB3 used in our model validation exercise consider the dominant plant functional type within this region to be "mixed C3 forest".Although this is consistent for the most part with land use maps from the region (Griffis et al., 2010), the proportion and composition of crops in this region may change from year to year to meet market demands.
It is also possible that air-sea gas exchange may contribute to the observed seasonal cycle in δ s at some of our Northern Hemisphere terrestrial sites.Although isotopic fractionation associated with air-sea gas exchange (ε ao = −2.0‰) is an order of magnitude less than isotopic fractionation associated with terrestrial biosphere-atmosphere exchange (ε ao ∼ −16.0 ‰), certain terrestrial sites, especially those near the ocean, may be impacted by marine derived air masses.This marine effect would also be the strongest at mid-latitudes that are subjected to strong onshore breezes as the continents heat faster than the ocean during summer months.Such a marine effect would lead to more enriched values (i.e. less isotopic discrimination) during summer months than winter months, which is indeed what is observed in the seasonal cycle of δ s .A marine effect may also help explain the complete lack of correlation between δ s and atmospheric water vapor at some sites, such as TAP, and also the reduced amplitude in seasonal δ s at some sites, such as PAL and TAP (Table 1).However, the fact that all terrestrial sites show coherent seasonal cycles, even sites located in continental interiors (see Supplement Fig. S1), suggests a common mechanism regulating biosphere-atmosphere interactions and imparting a strong isotopic signal.
Although factors such as fossil fuels, climate variability, land use, and air-sea gas exchange may have an impact on the seasonal cycles in δ s inferred from the atmosphere, the fact that these seasonal cycles in δ s are coherent across most sites in the northern hemisphere and that they are reproducible from year to year strongly suggests that they are driven by carbon exchange between the biosphere and atmosphere.

Evaluation of models
Our results suggest that most of the seasonal variability in δ s can be explained by changes in atmospheric water vapor over the growth season.Of the 19 sites investigated, 16 showed significant correlations with VPD and 6 showed significant correlations with RH.Globally, δ s was significantly correlated with both VPD and RH, but was more highly correlated with VPD.Conceptually, we would expect increases RH to lead to increases in stomatal conductance (i.e.increased c i /c a ), resulting in greater isotopic discrimination against atmospheric δ 13 CO 2 .In contrast, we would expect increases in VPD to lead to a decrease in stomatal conductance (i.e.decreased c i /c a ), resulting in reduced isotopic discrimination against atmospheric δ 13 CO 2 .Thus we would predict a negative relationship between δ s and RH, in contrast to a positive relationship between δ s and VPD.Where the sign of these predicted relationships is consistent with two stomatal conductance models evaluated here (Eqs.7 and 8).
At every terrestrial site in the Northern Hemisphere that we evaluated VPD was positively correlated with δ s , even at sites where the correlation is not significant (Table 1).Furthermore, the distribution of these correlation coefficients does not deviate significantly from normal (W = 0.92, p-value = 0.097), indicating that stomatal conductance responds to VPD according to our conceptual model (Fig. 4a).In contrast, the correlation coefficients between RH and δ s are both positive and negative (Table 1).Moreover, the distribution of correlation coefficients between RH and δ s deviate significantly from normal (W = 0.82, p-value = 0.0015), indicating that at some sites stomatal conductance is responding to RH as expected based on our conceptual model (i.e.negative correlations), but at other sites stomatal conductance contradicts our conceptual expectations (i.e.positive correlations).Lastly, many of the optimal correlation coefficients were based on the natural log of VPD or RH, which indicates a non-linear response of δ s and thus stomatal conductance.In the Leuning model stomatal conductance responds nonlinearly to VPD, whereas in the BWB model the response of stomatal conductance to RH is strictly linear.The more significant correlations observed between δ s and VPD than δ s and RH are consistent with our model evaluations indicating that the Leuning model is more accurate at predicting δ s than the BWB model.RMSE values were significantly lower for the Leuning model than the BWB model (p-value = 0.001, DF = 27).Thus, if the global seasonal cycles observed in δ s are in fact due to stomatal response to the environment, than the Leuning Model appears to be more effective at capturing this stomatal response on seasonal timescales.Although both models showed a significant relationship with δ s globally (Table 1), both models failed to predict the extremely depleted values of δ s observed at some sites (Fig. 3).This mis-match between observed and predicted δ s was slightly rectified by removing highly depleted values of δ s during winter months when net assimilation goes below zero, but  this problem seems to persist during the shoulder seasons of fall and winter when δ s values still remain fairly depleted (see Supplement Fig. S1).However, there was not enough parameter space allowed by either model to account for these highly depleted values during winter months.In fact, in order to obtain reasonable predictions for δ s we had to set the slope term in both models (i.e.m and m L ) to 25.0, which greatly exceeds any values from the literature (Sellers et al., 1996;Leuning, 1995).This mis-match between observations and predictions could be due to extreme physiological conditions that are not well simulated by the models, or else a depletion bias in inferred values of δ s at some of our sites.The site at which this mis-match was the most pronounced was at ASK in Algeria.Although seasonal values of VPD range from 0.5 to 3.0 KPa and values of RH range from 35 to 50 %, this site does not seem to experience climatic conditions that are any more extreme than other sites evaluated such as UTA or WIS.It is also possible that the extremely depleted values of δ s observed during the winter months at ASK are an artifact of this site being very distant from the tropospheric background site specified (i.e.NWR) when calculating δ s .Therefore, is more likely that the extremely depleted values of δ s are due to possible biases introduced during this analysis and not necessarily a deficiency of the models.
There is considerable debate within the ecophysiology literature as to whether stomates respond primarily to VPD or RH.Although there is empirical evidence at the forest stand scale that stomatal conductance in some instances responds more to VPD (Bowling et al., 2002) and in other instances responds more to RH (Wang et al., 2009), a consensus has yet to emerge as to what is the primary metric of atmospheric vapor to which plants are responding.Part of this lack of consensus may be due to the fact that VPD and RH are not independent variables and thus strong empirical relationships may emerge between stomatal conductance and both of these variables.Recent efforts have turned towards combining the empirical stomatal conductance models evaluated here with optimization models to gain greater insight into stomatal sensitivity to both CO 2 and H 2 O (Medlyn et al., 2011;Katul et al., 2010).In contrast, biosphere models seem to be converging on VPD as the physical mechanism driving stomatal conductance (Medvigy et al., 2009;Cramer et al., 2001) and coupled global carbon-climate models seem to be converging on RH as the physical mechanism driving stomatal conductance (Friedlingstein et al., 2006), indicating a disconnect in how stomatal conductance is formulated at different spatial scales.Here we have analyzed and presented a global dataset of atmospheric observations that may provide new insight as to how plants respond to the soil-atmosphere water continuum.
Both of these stomatal conductance models have been derived from empirical observations at the laboratory scale and have subsequently been validated using field observations.However, they were not explicitly designed for global applications and yet they are now being used to evaluate how the Earth's biosphere will respond to future changes in atmospheric CO 2 and concomitant climate change (Friedlingstein et al., 2006).Our results suggest that the Leuning model of stomatal conductance may be more suitable for simulating stomatal conductance over a wider variety of biomes under a wider variety of climatic conditions.However, this is not so surprising as the Leuning model includes more parameters and thus has more degrees of freedom for simulating actual observations.In contrast, the BWB model is elegant in its simplicity and actually performs reasonably well for some biomes and climatic conditions.Because of its simplicity it is not surprising that it has become the default stomatal conductance parameterization for many of the next generation Earth System Models.Furthermore, in order to implement both of these conductance models the continuous surface of the Earth's terrestrial biosphere, must be transformed into individual grid cells, and the wide array of biome's must be discretized into one of several plant functional types.Although the terms driving these respective models (VPD and RH) are used almost interchangeably in the ecophysiological literature, they are very different metrics of atmospheric water vapor and thus may respond very differently to future warming scenarios.
The degree to which surface RH and VPD change in response to atmospheric warming remains uncertain.Based on the Clausius-Clapeyron relationship a 1 • C increase in temperature should increase the atmosphere's capacity to hold water by ∼7 %.This relationship seems to hold true at the global scale, where a significant increase in specific humidity (kg H 2 O vapor/kg dry air) has been attributed to anthropogenic warming over the latter half of the 20th century (Willett et al., 2007).Although it is clear that as the atmosphere warms it contains more water vapor, the response of RH is much less clear.At the global scale there does not appear to be significant trends in RH (Willett et al., 2007) and in fact spatially and temporally invariant RH seems to be an emergent property of global climate models (Held and Soden, 2000).If in fact, specific humidity is increasing (i.e. the amount of water vapor in the atmosphere) and RH (i.e. the ratio of the amount of water vapor in the atmosphere to the amount of water that the atmosphere could potentially hold) then VPD must be increasing at the global scale.However, more recent data suggests that surface RH may in fact be declining over land, possibly due to limited ocean moisture as the Earth's ocean surface warms slower than the land surface (O'Gorman and Muller, 2010).However, there is considerable regional variability in changes in RH with very limited data from tropical regions (Simmons et al., 2010).Therefore, determining whether stomatal conductance responds toVPD that is increasing globally, or to RH, that may be decreasing regionally, is critical to predicting future carbon assimilation by the biosphere and thus realistic future climate scenarios.

Conclusions
Isotopic measurements in atmospheric CO 2 have greatly enhanced our understanding of the global carbon cycle.However, the increased sampling frequency of δ 13 CO 2 (both in space and time) combined with new analytical techniques, now make it possible to address research questions that were previously intractable using isotopes.Here we have presented a novel application of δ 13 CO 2 , by which we are able to infer broad seasonal patterns of source CO 2 to the atmosphere.We have used this approach to yield a seasonal cycle of δ s that is coherent across an array of N. Hemisphere atmospheric sampling sites.The resulting pattern can be thought of as a seasonal "Keeling Plot" for the Earth.The broad seasonal coherence observed across all sites suggests a single underlying physical mechanism driving this variability.To explain this variability we test two stomatal conductance models, which both suggest that changes in atmospheric water vapor drive changes in δ s on seasonal timescales.The analytical approach here can be improved by specifying more regional background reference curves essential for calculating δ s values.This analytical approach could also be extended to the tropics and the Southern Hemisphere as atmospheric observations of δ 13 CO 2 become increasingly available.

Figure 1 .
Figure 1.Analytical approach to inferring δ 13 CO source to the atmosphere.Panel A shows the seasonal cycle of Fig. 1.Analytical approach to inferring δ 13 CO 2 source to the atmosphere.Panel (A) shows the seasonal cycle of atmospheric CO 2 observations at Wendover, Utah (UTA) compared with free troposphere observations of CO 2 at Niwot Ridge, Colorado (NWR).Panel (B) shows atmospheric δ 13 CO 2 observations at UTA compared to free troposphere observations of δ 13 CO 2 at NWR. Panel (C) shows the residual technique used for inferring changes in δ 13 CO 2 source (δ s ) to the atmosphere, where δ s is estimated as the slope.Lastly, Panel (D) shows the resultant seasonal pattern in δ s generated from monthly estimates and their associated uncertainties (1σ ).

Fig. 2 .
Fig. 2.Global coherence of δ 13 CO 2 source (δ s ) anomalies to the atmosphere.The network of sites included in our analysis (A), where the size of points corresponds with length of the dataset (see Table1) and the color of the point corresponds with latitude and also the seasonal distribution of points in panel (B).The seasonal distribution of δ s for all sites included in our analysis (B).

Fig. 3 .
Fig. 3. Correlation analysis δ 13 CO 2 source (δ s ) and metrics of atmospheric water vapor.The color of points corresponds with the color code for individual sites denoted in Fig. 1.Global correlations for growing season vapor pressure deficit and δ s (A).Global correlations for growing season relative humidity and δ s (B).Comparisons between observed δ s and predicted δ s values for the Leuning Model driven by vapor pressure deficit (C) and for the BWB Model driven by relative humidity (D).See Table 1 for global and local statistics and see Table 2 for model parameters.

Fig. 4 .
Fig. 4. Density functions of Pearson correlation coefficients.Correlation coefficients are reported for relationships between δ s and vapor pressure deficit (A.) as well as relative humidity (B.).Probability density functions (black lines) are superimposed on histograms (boxes) where positive relationships are indicated by positive correlation coefficients (grey bars) and negative relationships are indicated by negative correlation coefficients (white bars).See Table1for actual correlation coefficients.

.
See Table 1 for global and local statistics and see Table 2 for model parameters.
Table 1 for actual correlation coefficients.