A possible global covariance between terrestrial gross primary production and 13C discrimination: Consequences for the atmospheric 13C budget and its response to ENSO

It is well known that terrestrial photosynthesis and 13C discrimination vary in response to a number of environmental and biological factors such as atmospheric humidity and genotypic differences in stomatal regulation. Small changes in the global balance between diffusive conductances to CO2 and photosynthesis in C3 vegetation have the potential to influence the 13C budget of the atmosphere because these changes scale with the relatively large one‐way gross primary production (GPP) flux. Over a period of days to years, this atmospheric isotopic forcing is damped by the return flux consisting mostly of respiration, Fire, and volatile organic carbon losses. Here we explore the magnitude of this class of isotopic disequilibria with an ecophysiological model (SiB2) and a double deconvolution inversion framework that includes time‐varying discrimination for the period of 1981–1994. If the net land carbon sink and plant 13C discrimination covary on interannual timescales at the global scale, consistent with El Niño‐induced drought stress causing a decline in global GPP and C3 discrimination, then less interannual variability in ocean and land net carbon exchange is required to explain atmospheric trends in δ13C and CO2 as compared with previous studies that assumed discrimination was invariant.


Introduction
[2] El Niño-induced changes in ocean circulation and global climate affect atmospheric CO 2 concentrations through several different mechanisms. In the eastern Pacific, a shutdown in equatorial upwelling caused by a relaxation of trade winds leads to a decrease in mixed layer pCO 2 levels and in the ocean-to-atmosphere CO 2 flux [Feely et al., 1997[Feely et al., , 1999Winguth et al., 1994]. In terrestrial ecosystems, decreased strength of the Asian monsoon imposes widespread drought stress in southeast Asia and eastern Australia, the northern part of South America, Central America, and the Sahel in Africa [Dai et al., 1997;Janicot et al., 1996;Kripalani and Kulkarni, 1997]. Droughts in these regions are not fully compensated by wet spells in other regions: at the global scale precipitation decreases on land during El Niño events [Dai et al., 1998]. The combination of drought stress and increased temperatures appears to increase carbon loss on land [Braswell et al., 1997;Tian et al., 1998;Potter and Klooster, 1999], although relative contribution of changes in gross primary production (GPP), ecosystem respiration, and fires is not well known.
[3] On El Niño-Southern Oscillation (ENSO) timescales (2 -7 years), changes in atmospheric d 13 C (in conjunction with the growth rate of CO 2 ) have been used frequently as a diagnostic of interannual variability in ocean and land carbon fluxes [Ciais et al., 1995;Enting et al., 1995;Francey et al., 1995;Keeling et al., 1995;Battle et al., 2000]. Known as a ''double deconvolution,'' the d 13 C-derived record of terrestrial net carbon exchange has served as target (standard) for terrestrial ecosystem models attempting to reproduce interannual and decadal flux variability [Maisongrande et al., 1995;Thompson et al., 1996;Gerard et al., 1999;Potter and Klooster, 1999]. Previous applications of the double deconvolution approach have assumed that plant discrimination against 13 C during photosynthesis remains constant from year-to-year. This approximation has been widely employed on many timescales [Joos and Bruno, 1998;Trudinger et al., 1999;Battle et al., 2000]. However, the potential of time varying discrimination to induce isotopic disequilibria has been recognized for several years [Lloyd and Farquhar, 1994;Francey et al., 1995;Ciais et al., 1995;Flanagan et al., 1997]. Studies of environmental and biological controls on 13 C discrimination during photosynthesis have revealed a number of factors that cause discrimination to vary [Farquhar et al., 1989]. These include atmospheric humidity, solar radiation, drought stress, and plant type, all of which may be expected to respond to interannual climate variability.
[4] With the assumption of constant discrimination, the governing equations describing the atmospheric 13 C budget simplify so that discrimination affects only the net land flux term and not the gross photosynthetic flux (GPP) term. Here we assess evidence that the global terrestrial biosphere, as a whole, experiences drought-stress during El Niño events. We also summarize theoretical and empirical evidence that drought-stress leads to predictable responses for GPP and 13 C discrimination in C3 plants. We then provide a general framework for assessing the magnitude of isotopic disequilibria induced by time varying terrestrial 13 C discrimination. We show that isotopic budgets for land need to include discrimination anomalies (much less than 1%) interacting with the one-way total photosynthetic flux (which at the global scale has been estimated between 100-150 Pg C/yr). We suggest that a significant fraction of the accelerated decline in atmosphere d 13 C during latter stages of El Niño events may be attributed to a global decline in discrimination in tropical C3 ecosystems. Allowing discrimination to covary with GPP, in a double deconvolution inversion (presented below) suggests that land/ocean interannual flux variability may be substantially smaller than that predicted in previous analyses.

Climate Variability and Global Ecosystem Function
[5] To assess interannual changes in drought-stress at the global scale, we generated a monthly time series of precipitation (PPT), surface air temperature, and vapor pressure deficit (VPD) that were weighted spatially by net primary production (NPP) fluxes from the terrestrial biosphere. The rational for this weighting function followed from our interest in the global budget of atmospheric 13 C: anomalies in climate over land may have greater impact on the atmospheric budget when they occur in places with high plant productivity.
[6] We used the Global Precipitation Climatology Project Version 2 monthly product for PPT [Huffman et al., 1997;Susskind et al., 1997], the Goddard Institute for Space Studies data set for surface air temperature [Hansen and Lebedeff, 1987;Hansen et al., 1999], and NCEP reanalysis products of 2-m relative humidity and air temperature for our VPD estimates [Kistler et al., 2001]. For the VPD estimates, we used surface air temperature as a proxy for internal leaf temperatures and monthly mean temperature and humidity data from NCEP (hence the relatively low estimates of VPD presented here). For the NPP-weighting function, we used the Carnegie-Ames-Stanford Approach (CASA) terrestrial biosphere model [Randerson et al., 1997]. NPP in CASA is derived from satellite-derived estimates of absorbed photosynthetically active radiation (PAR) [Bishop and Rossow, 1991;Los et al., 1994] and a light use efficiency term (units of g C/MJ PAR) that has a partial dependence on local environmental conditions. We used normalized estimates of the Southern Oscillation Index from the Australian Bureau of Meteorology (www.bom.gov. au/climate/current/soihtm1.shtml).
[7] We constructed the NPP-weighted time series of climate variables (PPT, temperature, or VPD; denoted generically by W ) using the following equation: where P(x,i) represents the climatological mean NPP value for each grid cell, x, and month, i, and C(x, i, t) represents the climate variable (PPT, temperature, or VPD) at each x, i, and year, t. We analyzed the 1980 to 2000 period, which included four El Niño events (1982 -1983, 1987, 1991-1992, and 1997 -1998), and two major volcanic eruptions (1982 and 1991) [Hansen et al., 1999].

C Discrimination
[8] Variability in 13 C discrimination, D ab , by photosynthesis results from variability in discrimination along the pathway of CO 2 flux from the atmosphere through surface boundary layer, the stomatal pore, the cellular milieu and ultimately from fixation by rubisco, the primary carboxylation enzyme in photosynthesis.
where g b , g s and g m are the conductances to CO 2 transport in the leaf boundary layer, stomatal pores and the cellular liquid phase (in the mesophyll) respectively. C a , C s , C i and C c are the CO 2 mole fraction in the atmosphere, leaf surface, intercellular spaces and chloroplasts respectively. b, s, m, and e represent discrimination factors associated with each step in the pathway (see Farquhar et al. [1989] for values). Rubisco discriminates strongly against 13 CO 2 but overall discrimination is modulated by the 13 CO 2 levels within the chloroplasts that in turn are influenced by discrimination by diffusion along the gas and liquid phase pathways leading to the chloroplasts. C4 plants exhibit lower 13 C discrimination because rubisco is isolated from diffusive interaction with the atmosphere as a result of the CO 2 concentrating mechanism found in these plants which discriminates relatively little [Farquhar et al., 1989].
[9] Studies of leaves and plants reveal that though discrimination capacity of rubisco is relatively constant among species, the 13 C/ 12 C ratio of CO 2 in the chloroplasts may vary in response to environmental conditions and species characteristics causing overall discrimination to vary. Perhaps the most well studied and significant source of variability in C3 plants occurs from the diffusion of CO 2 through the stomata [Berry, 1987]. Lower stomatal conductances cause the intercellular CO 2 concentration during steady state photosynthesis to fall and the 13 C/ 12 C ratio of that CO 2 to increase relative to the outside atmosphere. Stomatal conductance is strongly coupled to photosynthesis (i.e., rubisco activity) in a way that tends to maintain intercellular CO 2 concentrations (and thus the internal 13 C/ 12 C) at fairly constant levels [Wong et al., 1979]. However, as the vapor pressure gradient between the leaf and the atmosphere increases, stomatal conductance declines relative to photosynthetic capacity reducing intercellular CO 2 , increasing the 13 C/ 12 C of internal CO 2 , countering discrimination by rubisco and lowering overall discrimination by photosynthesis. An equation developed by Farquhar et al. [1982] provides a quantitative basis for predicting these changes in discrimination as a function of intercellular CO 2 . Other evidence suggests that water stress and high solar radiation levels can also lead to lower overall discrimination [Farquhar et al., 1989]. There is also a tendency for herbaceous species to discriminate more than woody species as a result of interspecific differences in stomatal conductance relative to photosynthetic capacity [Ehleringer, 1993].
[10] To illustrate and quantify how environmental and physiological conditions affect discrimination, we used the SiB2 land surface model that includes photosynthesis and stomatal conductance parameterizations [Sellers et al., 1996b] allowing diagnosis of component and overall discrimination during photosynthesis at hourly time steps. The previously described version of SiB2 has been updated to include a diffusive conductance for CO 2 between the intercellular spaces and the sites of rubisco catalyzed carboxylation. This conductance is in series with CO 2 diffusion across the leaf boundary layer and the stomata (N.S. Suits et al., Simulating seasonal and spatial variations in global concentrations and carbon isotopic ratios of atmospheric CO 2 , in preparation for Global Biogeochemical Cycles, 2002). The value of mesophyll conductance, g m , is a function of the maximum capacity for photosynthesis (rubisco activity or V max [see Evans and Loreto, 2000]), a water availability scaling factor, w(q) that depends on water filled pore space in the soil, and the canopy integration factor, Å [Sellers et al., 1996b]. The canopy integration factor, Å, scales the mesophyll conductance for a leaf at the top of the canopy to a bulk canopy value analogous to the approach used to calculate canopy photosynthesis and stomatal conductance [Sellers et al., 1996a]. Here n is a scaling parameter tuned so that C i À C c was about 50 ppm mole fraction for light saturated unstressed photosynthesis.
In SiB2 stomatal conductance is represented by where h s and b c are canopy bulk leaf surface humidity and minimum conductance respectively. m bb is an empirical constant that can be viewed as the reciprocal of the intrinsic water use efficiency (and in the definition of terms given here it also includes the diffusion coefficient of CO 2 relative to water vapor). The smaller the value of m bb the higher the water use efficiency, the lower C i /C a and the lower D ab becomes. Canopy photosynthesis, A c , is a function of C c , short wave solar radiation, temperature, and soil water availability (and V max ). In SiB2 A c is equivalent to GPP minus leaf respiration.
[11] We used meteorological conditions [Sellers et al., 1989] and vegetation parameters [Sellers et al., 1996a] from an Amazonian forest to test the sensitivity of discrimination to environmental conditions and model parameters. Sensitivities (S) are expressed as where V is the dependent variable [e.g., discrimination ( 13 D ab ) or canopy photosynthesis (A c )] and P is a parameter or independent variable. We also calculated the sensitivity of D ab for a given change in A c caused by a change in P.
While not strictly a sensitivity estimate, it does reflect how discrimination and photosynthesis change relative to one another which is important for evaluating the influence of vegetation on atmospheric CO 2 and d a as will be shown below. Environmental conditions (atmospheric humidity, solar radiation) and model parameters (FPAR, V max , intrinsic water use efficiency (1/m bb ), and n) were reduced individually by 20% for the calculation of respective S values. Water stress response was estimated by initializing model simulations at lower soil water content causing the water stress factor to increase by 14% and reducing A c by about 10%. Model parameters (FPAR, V max , m bb , and g m ) can be considered as surrogates for physiological/biological controls on discrimination. Dependent variables were expressed as averages of hourly values weighted by A c .
2.3. Implications for the Global Atmospheric 13 C Budget 2.3.1. Conventional Double Deconvolution Approach [12] The atmospheric mass balance of total C in CO 2 can be described by the following equation: where the atmospheric rate of change d(C a )/dt is equal to fossil fuel (F f ) release and uptake by ocean (N o ) and land (N b ) net carbon sinks. With equation 1, N o and N b are unknown and so an additional constraint (from O 2 /N 2 ratios [Keeling et al., 1996] or 13 C [Tans et al., 1993]) is required to distinguish between the two sinks. For the case of 13 C, a second equation can be written in d notation that captures the difference in isotopic fractionation associated with photosynthesis and air-sea gas exchange : where the rate of change in atmospheric 13 C is approximately equal to isotopic changes imposed by a terrestrial carbon sink, terrestrial disequilibria forcing (caused by the changing atmospheric isotope composition), fossil fuel isotope forcing, ocean uptake, and ocean disequilibria forcing (see Table 1 for a definition of the various terms in equations (6) - (9)).
[13] Using atmospheric measurements, along with characterization of fossil fuel 13 C/ 12 C composition, fractionation by photosynthesis and ocean exchange, and models of carbon turnover in terrestrial and oceanic reservoirs, it is possible to solve these two equations to estimate terrestrial and ocean carbon sinks at various spatial and temporal scales.

Proposed Framework for Variable Plant Discrimination
[14] The representation of biosphere-atmosphere 13 C fluxes described by equation (7) includes an implicit assumption that discrimination against 13 C during photosynthesis is invariant. A more general description of biosphere-atmosphere exchange allows for the possibility of a changing discrimination: where the net land flux, N b , is now replaced with the sum of GPP (G b ) and the biosphere-atmosphere return flux (R b ). While R b is largely composed of plant and heterotrophic respiration (R a and R h ), fire and volatile organic carbon losses are also included in this term. Discrimination is now defined as the sum of an invariant mean, D ab -, and a time $À125 Pg C/yr GPP gross atmosphere-terrestrial biosphere flux (GPP) G b 0 $ À32 Pg C/yr ( Figure 5) In Figure 5 we use G b 0 instead of G b for illustration purposes: We assume that all C3 plant respiration and 1/3 of heterotrophic respiration occurs during the same year as G b , and so this component of the flux does not impact the annual 13 C or CO 2 budgets illustrated in Figure 5. Also, we only considered a change in the C3 component which accounts for $77% of global NPP (section 2.3.3).

R b
$125 Pg C/yr gross terrestrial biosphere-atmosphere return flux. R b 0 $ 31 Pg C/yr ( Figure 5) In Figure 5, we isolate the component of R b that has a residence time in the biosphere that is longer than 1 year, denoted here by mean C3 13 C plant discrimination. We assume in analysis of Figure 5 that the terrestrial sink is C3. Changes in C4 productivity are discussed in section 4.3. _ Á ab 0.01*Á ab = 0.19% anomaly in C3 13 C discrimination. In Figure 5 we explore the implications of a discrimination anomaly that is 1% of the mean C3 value. d ba d ba À Á ab = À7.8%À19.0% d 13 C composition of the terrestrial biosphereatmosphere return flux in Figure 5 (red arrow). " e ao À1.8% mean 13 C discrimination associated with net ocean uptake 33% isotopic disequilibria of the terrestrial biosphereatmosphere flux varying anomaly, Á D ab . In Appendix A, we show that if Á D ab is assumed to equal to 0, then equation (9) simplifies to equation (7). However, when Á D ab is not equal to zero, the product G b Á D ab has the potential to induce an atmospheric forcing. This atmospheric forcing is erased as the isotopic composition of respiration (d ba ) adjusts to (equilibrates with) the new value of discrimination.
[15] With equations (8) and (9), there are now five unknowns, Á D ab , R b , d ba , G b , and N o , and two constraints: the atmospheric rate of change of CO 2 and d 13 CO 2 . Three additional pieces of information are needed to solve this more general system of equations. They are described in the following three subsections.

Anomalies in Plant Discrimination
[16] Variability in Á D ab arises from physiological, ecosystem, and biome-level factors. A model that couples stomatal conductance and photosynthesis may capture some of the response of internal leaf CO 2 concentrations and 13 C discrimination to short-term environmental variability at leaf and canopy levels. When environmental perturbations persist for a period of weeks to seasons, the integrated ecosystem response will also reflect changes in the contribution to the gross flux from species with different intrinsic growth rates, allocation patterns, and water use efficiencies. At the biome level, changes in the disturbance regime can affect the distribution of different plant functional types such as C3 trees and C4 grasses.

Magnitude of the Biosphere-Atmosphere Return Flux
[17] At a given time interval, it may be possible to estimate R b from an ecosystem model that takes into account carbon accumulation (the past history of G b ), allocation, and if available, contemporary information on rates of decomposition, plant respiration, fire, and other loss pathways.

Isotopic Composition of the Biosphere-Atmosphere Return Flux
[18] Accurate estimates of d ba are challenging because they require integration of the past history of G b , Á D ab , " D ab , and d a , and information about biosphere-atmosphere loss pathways. For example, the carbon released by fire has a different isotopic composition than that released from the same ecosystem by decomposition because of differences in the chemical composition and ages of substrates consumed by the two processes [Schuur et al., 2002]. 2.3.3. Vector Diagram of 13 C and CO 2 Budgets for 1990 [19] To illustrate the difference between conventional (equations (6) and (7)) and general (equations (8) and (9)) descriptions of terrestrial biosphere isotopic exchange, we constructed a mean annual 13 C and CO 2 atmospheric budget for a single year (1990) using data from multiple sources [Andres et al., 2000;Francey et al., 1995;Fung et al., 1997;Gruber and Keeling, 2000;Tans et al., 1993] as summarized in Table 1. We present these budgets using the vector diagram approach developed by Enting et al. [1993]. For the conventional budget diagram, terrestrial fluxes were partitioned into net sink and disequilibrium components. For the general budget diagram, we estimated the component of the gross biosphere-atmosphere flux that would retain a discrimination anomaly after 1 year. With both the general and conventional budgets, we considered the effect of a 1% change in global terrestrial C3 discrimination over a year (i.e., a global change from 19.0% to 19.19%).
[20] We assumed that 77% of GPP is C3 [Still et al., 2002], NPP is $50% of GPP [Ryan, 1991;Woodwell and Whittaker, 1968], and that all plant respiration is respired in the same year as initial fixation so these fluxes have no impact on the annual atmospheric 13 C budget. Of the remaining C3 NPP, we assumed that $1/3 was released via herbivory, microbial decay, and fire during the same year as the initial fixation (see paragraph below). The component of C3 NPP that remained in the ecosystem after 1 year was assumed to consist of carbon allocated to wood (with very long residence times), and the fraction of leaves and fine roots not consumed by herbivores, microbial decay, or fires. With these assumptions, and starting with a global GPP of 125 Pg C/yr, we obtained a C3 NPP flux of 32 Pg C/ yr that retained our hypothetical 1% discrimination anomaly after 1 year. We denote this decimated flux as G b 0 .
[21] A precise estimate of G b 0 is difficult because it requires accurate assessments of plant allocation, longevity estimates of fine roots, leaves, and stems, decomposition rates of these plant tissues, and the distribution of plant functional types at regional and global scales. Compared with other biosphere models used in isotope studies, our estimate of G b 0 is conservative (low). For example, with the five-box model of the biosphere constructed by Emanuel et al. [1981] and employed by Francey et al. [1995], approximately 70 Pg C (out of an original GPP flux of 120 Pg C/ yr) remained in the biosphere after 1 year [Thompson and Randerson, 1999]. Similarly, with the CASA model approximately 45 Pg C (out of an original NPP flux of 55 Pg C/yr) remained in the biosphere after 1 year [Thompson and Randerson, 1999].
[22] A significant fraction of G b 0 consists of NPP allocated to wood (that remains intact within plants typically for a period far in excess of 1 year). Using the CASA biosphere model as a tool to help with scaling, we found that between 7.8 Pg C/yr and 15.1 Pg C/yr of global NPP was allocated to wood, with the variation depending largely on whether a fixed biome vegetation map (upper estimate [DeFries and Townsend, 1994]) or a map of fractional woody vegetation cover (lower estimate [DeFries et al., 2000]) was used with the model allocation scheme (G. van der Werf, personal communication, 2002).
[23] The physiological and ecosystem studies that provide the underpinnings for these models also provide direct evidence that our estimate of G b 0 is conservative. Tropical dry forest trees, temperate (cold) deciduous trees, and perennial grasses store starches and sugars in stems and boles at the end of one growing season for the purpose of rapid canopy construction at the onset of the following growing season [Aerts and Chapin, 2000]. Thus, some small component of plant respired C (autotrophic respiration) has a residence time of a year or more. Of the remaining carbon that contributes to a cohort of NPP, rapid loss pathways are likely to be leaf and fine root decomposition, herbivory, and fire. For the case of litter decomposition losses to the atmosphere, the C residence time must also include the lifetime of leaves and fine roots. For tropical and boreal evergreen species, leaf lifespan typically varies between 1 and 2 years [Aerts and Chapin, 2000]. Even for deciduous trees, forbs, and grasses, the mean leaf lifespan is about 0.5 years [Aerts and Chapin, 2000]. Decomposition of leaf litter adds another 0.3 to 0.75 years to the biosphere residence time of carbon allocated to leaves in tropical ecosystems and from 1 to 10 years in temperate and boreal ecosystems [Aerts and Chapin, 2000]. Root longevity and rates of decomposition are more uncertain, and occur on a timescale of weeks to years. A considerable fraction of fine roots appear to live for periods greater than 1 year in temperate forests [Gaudinski et al., 2000]. Finally, fires release no more than 5 Pg C/yr to the atmosphere [Crutzen and Andreae, 1990], and it is likely that a large component of fire emissions have a residence time in the terrestrial biosphere longer than 1 year.

A Simplified Double Deconvolution Inversion With Variable Discrimination
[24] We made two simplifications that allowed us to solve for the 13 C and total CO 2 budgets using equations (8) and Figure 1. We solved a variable discrimination inversion with CO 2 and d 13 C using an iterative approach with equations (8), (9), and (10). The steps in the iteration are outlined with this flowchart. Prior to 1750, the atmosphere was assumed to be in steady state isotopically with the terrestrial biosphere (no disequilibria).
(9). First, if GPP and discriminations do covary at regional and global scales in response to drought stress, then it may be possible to write a linear equation relating discrimination anomalies to GPP anomalies: where S A D is a linear scaling factor that describes how an annual anomaly in GPP translates into an anomaly in discrimination, i.e., a sensitivity factor. S A D will vary over the land surface depending on what environmental variables are regulating productivity and the timescales that are considered (see sections 2.2, 3.2, and 4.2). Where concurrent productivity and isotope data are (or become) available, direct measurement of S A D provides a means for testing the tropical and global scale discrimination-GPP covariance hypothesis presented here.
[25] In our preliminary analysis presented here, we used S A D values of 0.0 and 0.50 in an iterative double deconvolution process described below (Figure 1). A S A D value of 0.0 represents invariant discrimination, i.e., the conventional double deconvolution approach (equations (6) and (7)). A S A D value of 0.5 corresponds to the following: a 1% change in annual GPP causes a 0.5% change in annual (flux-weighted) discrimination. The 0.5 S A D value we adopted is conservatively consistent with the response of C3 tropical forest in SiB2 to the sum of VPD and soil moisture effects (sections 2.2, 3.2, and 4.2) and is broadly consistent with many (but not all) tree ring studies that measured both growth increments and wood isotopic composition (see section 4.2 for a more detailed discussion).
[26] The second simplification is that we used a 1-D pulse response (Green's function) model of the terrestrial biosphere with a monthly time step to estimate R b as a function of the past history of GPP and d ba as a function of the past history of G b , Á D ab , " D ab , and d a . Our pulse model only considered the C3 terrestrial biosphere, and started with C3 GPP of $90 Pg C/yr in 1750. From 1750 to 1981, ice core records of d a  were used as a boundary condition for the pulse model, to properly capture the isotopic disequilibria of the changing atmosphere (or Suess Effect). Similarly, we increased NPP from 1750 to 1982, following atmospheric CO 2 (with a b factor of 0.4 [Wullschleger et al., 1995]), to generate a difference between G b and R b required by the double deconvolution approach at the start of the 1981-1994 period.
[27] With this simple pulse-response approximation, interannual variation in climate did not affect the magnitude of R b , and we assumed that C4 ecosystems did not contribute to interannual variation in gross or net terrestrial exchange. These assumptions allow for a simplistic representation of the terrestrial biosphere, allow us to demonstrate the consequences of C3 discrimination variability for the global carbon budget, but are insufficient for a comprehensive assessment of changes in atmospheric d 13 C and CO 2 (see section 4.4).
[28] For this analysis, we used a smoothed atmospheric record of d 13 C and CO 2 record from Cape Grim archived air from 1982 to 1994 . We solved Figure 2. (a) Annual NPP from the CASA model is greatest in equatorial regions of South America, Africa, and Asia (g C m À2 yr À1 ). This NPP product is derived primarily from satellite-derived products of normalized difference vegetation index (NDVI) and solar radiation (see text for details). (b) Global precipitation (PPT) from the GPCP Version 2 monthly product (mm yr À1 ) has a very similar global pattern to NPP. (c) Pearson correlation coefficient, r, between annual PPT and a normalized Southern Oscillation Index (SOI) (unitless) for the 1980 -2000 period. In positive correlation areas (denoted by red), rainfall decreases during El Niño events. Note the relatively high correlation coefficients in northern South America, Southeast Asia, and western equatorial Africa, where both mean annual NPP and PPT are large. (d) The fractional difference in mean annual PPT between a La Niña period (1999 and 2000) and an El Niño period (1997 and 1998) (unitless). equations (8), (9), and (10) using an iterative process (Figure 1).

Climate Variability and Global Ecosystem Function
[29] Mean annual NPP derived from satellite-derived estimates of canopy light absorption was highest in equatorial regions of South America, Africa, and Asia ( Figure  2a) and had a similar spatial pattern to mean annual PPT (Figure 2b). Interannual PPT was strongly correlated with the Southern Oscillation Index in northern South America, southern Africa, southeast Asia, and eastern and northern Australia (Figure 2c). Weak positive correlations occurred over western equatorial Africa and western boreal North America. Negative correlations occurred over the central part of North America, the southeastern part of South America, and central Asia. Decreases in PPT during the 1997 -1998 El Niño event were substantial over many regions (in excess of 35%) when compared with PPT during the 1999 to 2000 La Niña period (Figure 2d).
[30] Over the entire terrestrial biosphere (equation (1)), PPT changed substantially from year-to-year, with a 13% a) b) c) Figure 3. (a) NPP-weighted precipitation, (b) surface air temperature, and (c) vapor pressure deficit time series were calculated according to equation (1) (solid line with circles). The normalized Southern Oscillation Index (SOI) was highly correlated with rainfall over the terrestrial biosphere (dashed lines with squares) over the 1980 to 2000 period. NPP-weighted precipitation varied by over 13% during this period. The global terrestrial biosphere experienced drought stress during the 1982-1983, 1987, and 1997 -1998 El Niño events, as measured by low rainfall, high air temperatures, and high vapor pressure deficits.
difference between the maximum rate in 1989 and the minimum rate in 1992 (Figure 3a). The relatively strong 1997 -1998 El Niño event caused a global decrease in biosphere-weighted PPT of 5% in 1997 and in 1998, relative to the long-term mean. Biosphere-weighted PPT was also below average during other El Niño events of varying intensity: 1983, 1987, and 1991 -1992. Biosphere-weighted PPT had a significant positive correlation with the Southern Oscillation Index from 1980 to 2000 (r = 0.75, n = 19, P < 0.05) (Figure 3a).
[31] Biosphere-weighted temperatures reached a maximum in 1998 during the 1997 -1998 El Niño event. The lowest biosphere-weighted temperature during the 1980 -2000 period occurred in 1992, probably as a result of stratospheric aerosol loading from the Pinatubo volcanic eruption in June of 1991 (Figure 3b). Biosphere-weighted VPD obtained from NCEP data varied by over 11% from 1980 to 2000. The maximum VPD (0.58 kPa) occurred in 1998, while the minimum VPD (0.51 kPa) occurred in 1993 ( Figure 3c).

GPP and Discrimination Responses to Environment and Physiology
[32] The influences of atmospheric vapor pressure deficit (VPD) and FPAR on D ab and A c are illustrated in Figure 4. A c increased by more than 20 mmol CO 2 m À2 s À1 as a result of increasing FPAR from 0.2 to 0.9 (a change in leaf area index from 0.5 to 5) causing D ab to increase by about 0.7% to 1.1% depending on VPD. Positive feedback between A c , g c , and h s caused C c /C a and D ab to increase with FPAR. A 33% increase in VPD (from 16.4 hPa to 22.9 hPa) reduced D ab by about 0.4% to 0.8% depending on A c .
[33] Table 2 shows the sensitivities (S ) to environmental conditions likely to vary on interannual timescales, specifically, VPD, photosynthetically active radiation (PAR) and soil water availability. Increasing VPD causes stomata to close, decreasing C i /C a and C c /C a . Since both A c and D ab are functions of C c there is a near proportional change in D ab relative to A c as VPD changes. Reducing water stress has a smaller effect because in our model parameterization water stress affects the flux A c , and the conductances, g c and g m linearly which tends to maintain C c constant. The fact that there is moderate sensitivity is in part a result of feedback between A c , g c , h s and C i . That is, water stress reduces A c and g c causing reduced transpiration and increased leaf temperatures. The resulting lower h s reduces g c further relative to A c such that C i and thus C c and D ab decrease. Increasing PAR increases A c and g c but does not influence g m causing C c and D ab to decrease slightly.
[34] Model parameters reflect the characteristics of the vegetation. Specifically, FPAR indicates the amount of green canopy, V max the photosynthetic capacity of the leaves, m bb the reciprocal of the intrinsic water use efficiency of the leaves and n, the mesophyll capacity for CO 2 transport. These characteristics vary among species and for an individual plant depending on environmental conditions. Though in SiB2 m bb is constant for specific vegetation types, previous studies show that plants adapted or exposed to water stress conditions can have higher intrinsic water use efficiencies (lower m bb ) which would have the net effect of increasing D ab relative to GPP [Farquhar et al., 1989;Ehleringer, 1993]. FPAR and V max have only small impacts on S A D and n tends to be fairly conservative among various types of plants [Evans and Loreto, 2000]. The high sensitivity for changes in n does illustrate the potentially strong influence that the parameterization of g m can have on predictions of S A D .
3.3. Implications for the Global Atmospheric 13 C Budget 3.3.1. Vector Diagram of 13 C and CO 2 Budgets for 1990 [35] Different amounts of discrimination against 13 C by C3 photosynthesis ($19%) and by air-sea gas exchange Figure 4. C3 discrimination from the SiB2 model is highly sensitive to VPD. Each model trajectory was created by varying canopy LAI between 0.5 and 5.0. The upper trajectory (solid line) was obtained by using meteorological conditions for a tropical forest site from Sellers et al. [1989]. The lower trajectory (dashed line) was created by increasing mean VPD by 33% relative to the meteorological conditions used in the upper trajectory. VPD is the vapor pressure deficit of the atmosphere, PAR is the incident photosynthetically active radiation, FPAR is the fraction of PAR absorbed by the canopy, V max is the maximum rubisco capacity, and n is the mesophyll conductance scalar.
($2% as a global ocean average) uniquely determine the slope of the net land and ocean carbon sink vectors shown in Figure 5a. A 1% change in global 13 C discrimination (the slope of the dark green arrow) with the conventional double deconvolution representation (Figure 5a) will have a negligible impact on land/ocean sink partitioning: less than 0.02 Pg C/yr shift between land and ocean sinks.
[36] With the more general description of the terrestrial biosphere, a 1% change in annual C3 13 C discrimination (0.19%) has a substantial impact on the land/ocean sink partitioning: a 0.37 Pg C/yr shift between land and ocean sinks (equations (8) and (9)). This increased impact on land/ ocean sink partitioning occurs because the discrimination anomaly changes the slope of the relatively large gross atmosphere-biosphere terrestrial flux vector (G b 0 , note the change in scale between Figures 5a and 5b). We assumed that 2/3 of C3 NPP retained the 1% discrimination anomaly after 1 year (see section 2.3.1).

A Double Deconvolution Inversion With Variable Discrimination
[37] During the 1980s and early 1990s, smoothed trends in atmospheric CO 2 and d 13 C were highly variable ( Figures  6a and 6b). For example, the El Niño event of 1982-1983 was characterized by large growth rates of atmospheric CO 2 and large decreases in d 13 C. In contrast, during the early 1990s, CO 2 growth rates were smaller and d 13 C decreases almost ceased.
[38] In our inversion analysis, an almost neutral terrestrial biosphere in the early 1980s and a large terrestrial carbon sink in the early 1990s provided the best fit with atmospheric observations (Figure 6c). When discrimination was held constant, the terrestrial carbon sink ranged between À0.3 Pg C/yr and À2.7 Pg C/yr over the 13-year period. When discrimination was allowed to covary with GPP, the terrestrial carbon sink ranged between À0.4 Pg C/yr and À2.1 Pg C/yr. The effect of allowing discrimination to covary with GPP in our simplified inversion was to decrease the range of the terrestrial carbon sink by 0.7 Pg C/yr or 28%. The impact on the range of the ocean sink was more modest: a decrease in the range by 0.4 Pg C/yr or 14%. With the constant discrimination inversion, the mean difference between the ocean and land sink in any given year was 1.0 Pg C/yr, whereas with variable discrimination, the mean difference was 0.6 Pg C/yr (this represents over a 40% reduction).
[39] The range of GPP required to sustain a terrestrial carbon sink over this period (to keep photosynthesis outpacing respiration) was substantially less when discrimination was allowed to covary with GPP ( Figure 6d). In the inversion with constant discrimination, C3 discrimination was fixed at 19%. When GPP and discrimination were allowed to covary, discrimination of the C3 terrestrial biosphere ranged between 18.9% and 19.3% (Figure 6e). This range represents less than a 2% change in the discrimination of the C3 terrestrial biosphere.

Climate Variability and Global Ecosystem Function
[40] The PPT, temperature, and VPD data presented in Figures 2 and 3 suggest that the terrestrial biosphere, as a Figure 5. (a) Conventional double deconvolution inversion approach for estimating land and ocean carbon sinks for 1990. The brown vector represents fossil fuel emissions (6 Pg C/yr, À168 Pg C %/yr), the light green vector represents the land isotopic disequilibrium forcing (0 Pg C/ yr, 20 Pg C %/yr), and the light blue represents the ocean disequilibrium forcing (0 Pg C/yr, À55 Pg C %/yr). The black vector represents the observed change in atmospheric composition (3.0 Pg C/yr, À38 Pg C %/yr). Only one combination of land (dark green vector) and ocean (dark blue vector) net carbon sinks can close both the 13 C and total C budgets. (b) The land flux can be alternately represented as the difference between the two one-way gross exchanges of GPP (G 0 b , dark green vector) and the biosphere-atmosphere return flux (R 0 b , red vector). Note the change in axis scale between Figures 5a and 5b. The land isotopic disequilibria forcing term (arising from the Suess effect) is the same in Figures 5a and 5b, but in Figure 5b it is included as a part of R 0 b . The two approaches, Figures 5a and 5b, are equivalent when discrimination is invariant (Appendix A), but rapidly diverge when discrimination anomalies change the slope of G 0 b .
whole, undergoes drought stress during major El Niño events like the one in 1997 -1998. Biosphere-weighted PPT decreased, while biosphere-weighted temperature increased (Figure 3). Anomalously high air temperatures enhanced atmospheric driving gradients for water loss. The primary regions of the terrestrial biosphere that experienced below normal rainfall levels were Central America, northern South America, southern Africa, Southeast Asia, and northern and eastern Australia. These areas account for a substantial fraction of global terrestrial NPP (Figure 2a). Other Figure 6. Driving data and results from a 1-dimensional double deconvolution inversion with variable discrimination. (a) The atmospheric CO 2 content (left axis) and its annual rate of change (right axis) . (b) The atmospheric d 13 C composition (left axis) and its annual rate of change (right axis) . (c) An inversion with constant discrimination (solid lines) and with variable discrimination that varies linearly with global GPP (dotted lines). Ocean fluxes are denoted by blue and terrestrial biosphere fluxes are denoted by green. (d) Changes in GPP required to sustain a terrestrial carbon sink for the constant (solid) and variable discrimination (dotted) inversion cases. Plant respiration and decomposition were calculated using a pulse-response model of the terrestrial biosphere [Thompson and Randerson, 1999]. (e) Discrimination for the constant (solid) and variable (dotted) inversion cases.
analyses of PPT using station data from the last century suggest that the magnitude of PPT decreases in these regions is substantial during El Niño events [Dai et al., 1997;Janicot et al., 1996;Kripalani and Kulkarni, 1997;Ropelewski and Halpert, 1996] and that when PPT anomalies are integrated over all land areas, the signal remains strong [Dai et al., 1998]. Ropelewski and Halpert [1996] found that for large areas of central America, northeastern South America, India, southeastern Asia, and Australia, PPT decreased by about 20% to 40% during El Niño warm events relative to the long-term mean, and by twice that amount relative to the mean of La Niña cold events. Satellite observations during the last two decades confirm that leaf area (or FPAR) decreased during major El Niño events in these regions (and globally) and that terrestrial vegetation growth is linked with interannual patterns of sea surface temperature and PPT [Asner et al., 2000;Kogan, 2000;Los et al., 2001].
[41] Together, this climate and satellite evidence suggests that at the global scale, ENSO-induced change in water availability is the single most important factor that regulates interannual fluxes from the terrestrial biosphere. During the 1982During the -1983During the , 1987During the , and 1997During the -1998 El Niño events, the growth rate of atmospheric CO 2 significantly increased in the mid to latter stages of the event [Conway et al., 2001], although a quantitative partitioning of mechanisms remains elusive [Yang and Wang, 2000]. The terrestrial biosphere component of the anomalies probably involved some combination of decreased GPP [Tian et al., 1998], enhanced respiration [Braswell et al., 1997], and greater fire losses [Nepstad et al., 1999].

Environmental and Physiological Controls on 13 C Discrimination
[42] The SiB2 model simulations illustrate the short timescale physiological responses to drought that determine GPP and D ab in tropical C3 vegetation. Changes in D ab relative to A c (or GPP) are particularly significant in response to VPD, soil water availability, and intrinsic water use efficiency. Large-scale drought conditions brought on by an El Niño would include soil water stress, higher VPD and higher PAR all of which would tend to decrease D ab . In addition, increases in water use efficiency brought on by plant response to drought or by a reduction in herbaceous plant productivity (a shift toward woody species) can lead to lower D ab . At the ecosystem scale, plants with lower C i /C a may do proportionally better than plants with higher C i /C a in response to drought; i.e. drought stress over a period of months to years may amplify flux contributions from plants with higher water use efficiencies and lower discrimination values. The combination of these factors is likely to cause GPP and D ab to covary at interannual timescales.
[43] Interannual changes in precipitation water availability associated with ENSO can impose large changes in stand-level GPP and canopy conductance [Goldstein et al., 2000]. What is more difficult to predict is the response of ecosystem respiration relative to GPP, given that increasing temperatures will increase metabolic activity of plants and microbes [Ryan, 1991;Lloyd and Taylor, 1994], but decreased water availability will decrease growth respiration and the flow of dissolved organic compounds through the soil, a flux that is required for microbial metabolism [Neff and Asner, 2001;Raich and Potter, 1995]. In a temperate deciduous forest, decreases in respiration following a drought outpaced the decline in photosynthesis, leading to net CO 2 uptake by the ecosystem [Goulden et al., 1996]. Although interannual ecosystem flux data from tropical ecosystems is relatively sparse, evidence of deep taproots and access to deep soil water reserves in seasonally dry forests [Nepstad et al., 1994] suggests that in the tropics, ecosystem respiration may also be more sensitive to drought stress as compared with GPP.
[44] While the impacts of El Niño on ecosystem respiration are ambiguous, it is clear that one of the primary impacts of ENSO-induced drought stress in tropical forests is to increase fire frequency and intensity. Global estimates of biomass burning are between 2 and 4.5 Pg C/yr, and largely occur in seasonally dry tropical forests and savannas [Crutzen and Andreae, 1990]. While interannual variability in biomass burning is not well quantified, aircraft trace gas measurements (and in particular CO:CO 2 ratios) provide evidence that anomalously high levels of biomass burning occurred during the 1997 -1998 El Niño [Matsueda and Inoue, 1999]. Other impacts of El Niño in tropical forests include increased tree mortality [Villalba and Veblen, 1998] and forest impoverishment [Nepstad et al., 1999].
[45] Experimental evidence at the ecosystem scale for a covariance between GPP and discrimination in ecosystems that undergo drought stress comes from several sources. Saplings grown under low levels of water availability exhibit decreased instantaneous photosynthesis rates and water use efficiencies, long-term plant growth, and plant discrimination [Ares and Fownes, 1999;Zhang et al., 1997]. Measurement of phloem d 13 C increased by as much 8% from the wet to dry season in Australian eucalyptus forests, whereas irrigation fed stands exhibited almost no seasonal response and had an annual mean value that was depleted by over 3% as compared with the nonirrigated stands [Pate and Arthur, 1998]. Many tree ring studies report a negative correlation between tree ring growth increments (a surrogate for productivity) and wood isotopic composition in regions that undergo interannual variability in drought [Barber et al., 2000;Livingston and Spittlehouse, 1993;McNulty and Swank, 1995;Saurer et al., 1995]. It is important to note, however, not all tree ring analyses show a significant negative relations between growth indices and isotope composition [Flanagan et al., 1997;Stuiver et al., 1984], and in the tropics, where drought-stress is likely to be a major driver of productivity changes, there are few published interannual time series of productivity and plant 13 C discrimination from tree rings or other sources.
[46] Other isotopic work at the ecosystem scale suggests that the effects of water availability on discrimination are substantial and immediate. Over the course of a growing season in temperate forest ecosystems, the isotopic composition of soil respiration closely tracks water availability and can vary by up to 4 or 5% [Bowling et al., 2002;Ekblad and Hogberg, 2001; J. E. Fessenden and J. Ehleringer, Temporal variation in d 13 C of ecosystem respiration in the Pacific Northwest: Links to moisture stress, submitted to Oecologia, 2002]. Bowling et al. [2002] provide evidence that the link between humidity and the isotopic composition of soil respiration occurs rapidly (with a lag only of a few days) implying rapid changes in ecosystem water use efficiency and rapid C turnover within the soil.

Implications for the Global Atmospheric 13 C Budget
[47] If GPP and discrimination covary at regional and global scales, then less interannual variability in net land and ocean carbon fluxes may be needed to explain atmospheric d 13 C anomalies. During the later stages of the 1982 -1983 and 1997-1998 El Niño events, the rate of decline in atmospheric 13 C accelerated [Battle et al., 2000;. Our analysis suggests that part of this acceleration may be simply a result of some combination of vegetation VPD response and increased water use efficiency and thus lower than average discrimination in tropical C3 ecosystems. A negative anomaly in discrimination would provide less resistance to the dilution of atmospheric d 13 C by fossil fuels, independent of any changes in the net terrestrial flux.
[48] A negative anomaly in discrimination, while accelerating the decline in atmospheric d 13 C during the initial El Niño event, may have the opposite effect on the atmosphere in subsequent years because of feedbacks through the terrestrial biosphere. Although these feedbacks may not have been fully captured with our simple pulse-response biosphere model, they would work in the following way. Organic matter that was fixed during the initial drought stress event will be relatively enriched in d 13 C. As this carbon is released back to the atmosphere in subsequent years, it will tend to enrich the atmosphere as compared with respiration of an average isotopic composition.
[49] The initial response of global discrimination to El Niño along with feedbacks in subsequent years, have the combined effect in double deconvolution analyses of requiring less interannual variability in the net land carbon sink to account for atmospheric anomalies d 13 C. Because ocean and land carbon sinks are solved for simultaneously in equations (8), (9), and (10), ocean sink variability also decreases when terrestrial discrimination is allowed to vary with GPP. Thus, the mechanism described here may help to reconcile differences between the double deconvolution approach and other (less variable) methods used to estimate ocean fluxes [Lee et al., 1998].
[50] While changes in atmospheric O 2 /N 2 levels are consistent with large year-to-year changes in the terrestrial biosphere net sink derived from d 13 C during the 1990s [Battle et al., 2000], O 2 /N 2 levels may also respond to the large perturbations in ocean NPP and circulation that occurred during the decade. In the Battle et al. [2000] inversion of carbon sources and sinks, retrieval of land and ocean sinks with O 2 /N 2 assumed that all changes in O 2 / N 2 originated from the terrestrial biosphere and fossil fuels; ocean contributions were assumed to be invariant. Given the large changes in ocean circulation and biogeochemistry that occurred during the 1997/1998 El Niño, it is possible that oceans contributed to some of the atmospheric O 2 /N 2 variability. For example, Behrenfeld et al. [2001] found that ocean NPP increased by 5 Pg C/yr from 1997 to 2000, during an El Niño/La Niña transition, primarily in the eastern tropical Pacific.

Next Steps
[51] We designed our double deconvolution analysis to demonstrate the consequences of time-varying terrestrial C3 discrimination for the global atmospheric carbon budget. It was not meant to be a comprehensive treatment of d 13 C and CO 2 . We ignored several crucial processes. Specifically, shifts in C3 and C4 plant may also induce disequilibria [Ciais et al., 1999;Townsend et al., 2002;Still et al., 2002]. These shifts must be treated in the same way C3 fluxes are considered here, in terms of the isotopic forcing of gross photosynthesis and the gross biosphere-atmosphere return flux (mostly respiration, volatile organic carbon emissions, and fires). Interannual changes of a percent or two in the relative contribution of C3 and C4 GPP would have a comparable effect to the anomalies in C3 discrimination described here, and could either amplify or cancel variability in C3 discrimination caused by drought associated with ENSO. Indeed, many of the large drought stress anomalies associated with ENSO ( Figure 2c) are in areas with predominantly C3 vegetation. A decline in total C3 GPP relative to C4 GPP would have a similar impact on the global atmospheric 13 C budget as the drought in C3 ecosystems described here.
[52] In addition, we neglected variation in ocean disequilibria: by using a single average value from the work of Gruber et al. [1999]. As the ocean to atmosphere 13 C flux depends strongly on the d 13 C of pCO 2 in ocean surface layers [Gruber et al., 1999;Zhang et al., 1995], it will vary with changes in ocean biology [Behrenfeld et al., 2001], sea surface temperature [Reynolds and Smith, 1994], and salinity that accompany El Niño events.

Conclusions
[53] Small changes in terrestrial discrimination have significant impacts on regional and global 13 C budgets. While this complicates the direct use of d 13 C and CO 2 for land ocean partitioning, recent success in characterizing discrimination responses to climate at the ecosystem scale suggest that these responses are predictable, and in the future it will be possible to characterize and isolate their contribution. Three primary points of the paper are summarized below: 1. The terrestrial biosphere, as a whole, undergoes drought stress during major El Niño events.
2. Field observations, tree ring analyses, and ecophysiological models suggest that for many C3 ecosystems, GPP and discrimination simultaneously decrease in response to drought stress.
3. Relatively small changes in global 13 C discrimination (<1%) can have a significant impact on atmospheric 13 C content because they are associated with the one-way GPP flux. Equations describing isotopic exchange with the terrestrial biosphere need to allow for this possibility. If discrimination and the net land sink covary in response to El Niño (consistent with a decrease in C3 GPP), then less interannual variability in land and ocean sinks is required to simultaneously reconcile 13 C and CO 2 atmospheric budgets as compared with past applications of the double deconvolution approach that have assumed a constant discrimination.

Appendix A
[54] Considering only the terrestrial biosphere, 13 C fluxes can be represented in d notation as the sum of the isotopic forcing from the one-way fluxes of GPP and R b .
where D ab can be represented as the sum of an anomaly component,Ḋ ab , and a mean, " D ab : Assuming that G bḊab is equal to zero (that discrimination is invariant), A2 simplifies to: [55] Now, let the net land carbon sink, N b , equal the sum of GPP (G b ) and the biosphere -atmosphere return flux (R b ): Also, if plant discrimination is constant, then the isotopic composition of ecosystem respiration (d ba ) in equation A3 is given by: where d a e À d a is known as the isotopic disequilibrium, and is equal to the difference between the contemporary atmosphere, d a , and a past atmosphere, d a e , that represents the composition of the atmosphere when R b was fixed.
[56] Substituting equations A4 and A5 into A3, we obtain the conventional double deconvolution representation of the terrestrial biosphere: