The CarbonTracker Data Assimilation System for CO 2 and δ 13 C (CTDAS-C13 v1.0): retrieving information on land–atmosphere exchange processes

. To improve our understanding of the global carbon balance and its representation in terrestrial biosphere models, we present here a ﬁrst dual-species application of the CarbonTracker Data Assimilation System (CTDAS). The system’s modular design allows for assimilating multiple atmospheric trace gases simultaneously to infer exchange ﬂuxes at the Earth surface. In the prototype discussed here, we interpret signals recorded in observed carbon dioxide (CO 2 ) along with observed ratios of its stable isotopologues 13 CO 2 / 12 CO 2 ( δ 13 C). The latter is in particular a valuable tracer to untangle CO 2 exchange from land and oceans. Potentially, it can also be used as a proxy for continent-wide drought stress in plants, largely because the ratio of 13 CO 2


Introduction
The terrestrial biosphere has absorbed about 25 % of global fossil fuel carbon dioxide (CO 2 ) emissions over the last several decades but the future of this sink is highly uncertain in a warming world (Booth et al., 2012;Rowlands et al., 2012).It depends on the small difference between two large fluxes Published by Copernicus Publications on behalf of the European Geosciences Union.
of the terrestrial carbon cycle: photosynthetic uptake or gross primary production (GPP) and terrestrial ecosystem respiration (TER), and is here referred to as the net ecosystem exchange (NEE = TER − GPP + fire disturbances and land use change and harvesting of crops).All these flux terms respond to changes in local temperature, precipitation, nutrient availability, and other key environmental variables (Friedlingstein et al., 2006).Extreme climate events such as droughts can decrease GPP and increase TER and fire disturbances to a point where regional NEE is turned into a temporary carbon source (Ciais et al., 2005;Gatti et al., 2014;van der Laan-Luijkx et al., 2015).These dynamic responses (and positive feedbacks whereby increased CO 2 may lead to more droughts) are now an integral part of climate models that include fully coupled carbon cycling (Booth et al., 2012;Dai , 2012).Such models give rise to a wide range of climate projections primarily as a result of different simulations of terrestrial carbon exchange (Friedlingstein et al., 2006).It is therefore important to test and improve the representation of the terrestrial biosphere in carbon-climate models.Uncertainty in climate projections can be reduced by evaluating present-day performance of these models to observations (Hoffman et al., 2014).This paper presents a data assimilation system that can be used to evaluate existing terrestrial biosphere models by using an extensive number of atmospheric CO 2 observations in tandem with other trace gases.
Measurements of atmospheric CO 2 have been used to infer carbon fluxes at the Earth's surface using a variety of inversion techniques (e.g., Keeling and Revelle, 1985;Keeling et al., 1989;Tans et al., 1993;Ciais et al., 1995;Rayner et al., 2008;Alden et al., 2010).Unfortunately, a limited number of CO 2 observations, errors in atmospheric transport modeling, and the realism of bottom-up carbon flux estimates are limiting the utility of these techniques.For instance, the representation of subgrid scale vertical motion in (and through the top of) the planetary boundary layer is one of the most uncertain aspects in atmospheric tracer modeling and can hinder the accuracy of CO 2 transport (Kretschmer et al., 2012;Miller et al., 2015).In addition, atmospheric CO 2 as a tracer has its own limitations as it only reflects a small residual of different sources and sinks, such as wild fires, anthropogenic sources, ocean in-and outgassing, and terrestrial GPP and TER.
The CarbonTracker Data Assimilation System (CTDAS) has been developed to estimate global net ocean and terrestrial carbon exchange fluxes, with a focus on North America and Europe (Peters et al., 2005(Peters et al., , 2007(Peters et al., , 2010;;van der Laan-Luijkx et al., 2017).This application uses the ensemble Kalman filter (EnKF) as a Bayesian minimization approach for the estimation of weekly ocean and terrestrial carbon fluxes on a 1 • × 1 • horizontal grid to improve the agreement between modeled and measured atmospheric CO 2 .The versatile object-oriented design of CTDAS allows flexible implementation of different components of the data assimilation system (van der Laan-Luijkx et al., 2017).Such modifications include but are not limited to (1) the configuration of the state vector; (2) the expansion of the monitoring network, such as for the Amazon (van der Laan-Luijkx et al., 2015) and China (Zhang et al., 2014); (3) the use of Lagrangian atmospheric transport (He et al., in preparation); and (4) to monitor other tracer gases like methane (Bruhwiler et al., 2014;Tsuruta et al., 2016).
One aspect that has not yet been explored in CTDAS is the monitoring of multiple trace gases in the atmosphere that are strongly related (i.e., gases with a common chemical or metabolic pathway in the ocean and/or terrestrial biosphere).The main purpose of such an application is to improve the estimation of carbon fluxes and to retrieve new information on the underlying flux exchange processes that would otherwise remain undetected.We are in particular interested in the use of the stable isotope 13 C (in atmospheric CO 2 ) as an additional tracer alongside total CO 2 to estimate carbon sources and sinks and their variability.In earlier studies, 13 C was used to distinguish oceanic from terrestrial carbon exchange, as oceans take up 13 CO 2 more efficiently than land surfaces relative to 12 CO 2 .In so-called double-deconvolution methods this particular trait is used to untangle the global land carbon budget from the ocean carbon budget (Keeling et al., 1989;Tans et al., 1993;Ciais et al., 1995).More recently, the 13 C isotope was used to study the diurnal cycle of GPP and TER (Wehr et al., 2016) and was used as a tracer of water use efficiency to study long-term responses to CO 2 increases in tree rings (van der Sleen et al., 2015), and attempts are underway to do the same based on atmospheric records.On regional scales, variations in the ratio of 13 CO 2 / 12 CO 2 (typically reported as δ 13 C in ‰ relative to the VPDB reference ratio) reflect changes in discrimination processes associated with photosynthetic uptake of carbon by plants (e.g., Farquhar et al., 1989;Fung et al., 1997;Scholze et al., 2003;Rayner et al., 2008).Plants generally take up the heavier 13 CO 2 molecules less efficiently than 12 CO 2 molecules, increasing the 13 CO 2 / 12 CO 2 ratio of CO 2 remaining in the atmosphere.This kind of discrimination against 13 C is much stronger for C 3 plants than for C 4 plants, but also varies as a function of moisture conditions in the canopy air and soil (Farquhar et al., 1980(Farquhar et al., , 1989;;Ekblad and Högberg, 2001;Ometto et al., 2002;Suits et al., 2005).That implies that under the right circumstances, measured atmospheric δ 13 C can be used to recognize land usage, such as C 3 /C 4 photosynthesis, and changes in photosynthetic activity resulting from drought stress (Ballantyne et al., 2010;Raczka et al., 2016).
Such an application could also be beneficial to explore other facets of carbon exchange.Any errors in the fossil fuel emission inventories (although relatively small) are in the current CTDAS releases aliased erroneously on the natural ocean and terrestrial fluxes.Assimilation of the fraction of the radioactive isotope 14 CO 2 in the atmosphere would allow independent verification of the fossil fuel emissions as its old organic carbon is radiocarbon-free (Bozhinova et al., 2014;Basu et al., 2016).Other chemical constituents like carbonyl sulfide (OCS) and solar-induced chlorophyll fluo-rescence (SIF) could also be important additions to CTDAS.Inclusion of these tracers in the assimilation could enhance our understanding of carbon exchange, because variations in photosynthetic carbon uptake are recorded in atmospheric OCS and satellite SIF data (Yang et al., 2014;Commane et al., 2015).
Before we can interpret signals derived from these additional tracers, our aim for this paper is (1) to explain how the first dual-species CTDAS application works, with a specific focus on the use of δ 13 C and CO 2 , henceforth the system named as CTDAS-C13 version 1.0; (2) to demonstrate its accuracy in solving the targeted optimization problem in comparison to observations; (3) to test the sensitivity of the system to the introduced nonlinearity arising from simultaneous optimization of terrestrial total CO 2 and 13 CO 2 fluxes; and (4) to verify our new estimates of carbon and isotope exchange with independent drought index data.

Methodology
We present the atmospheric δ 13 C budget (Sect.2.1) before proceeding to describe the integration of δ 13 C within our new dual-species data assimilation framework CTDAS-C13 (Sect.2.2).We then briefly describe the prior estimates and the observational network used (Sect.2.3).Finally, we give a brief description of the different inversion experiments (Sect.2.4).The methodology presented here is based on Sect.4.2 of the lead author's PhD dissertation (van der Velde, 2015).

Atmospheric δ 13 C budget
The use of δ 13 C observations alongside CO 2 observations constitute a useful change to the traditional CO 2 -only CT-DAS application, as it provide an additional constraint on carbon surface fluxes and isotope exchange processes in plants.The rationale behind this is that the 13 CO 2 and 12 CO 2 contents in the atmosphere are affected through the same CO 2 pathways from land and ocean surfaces.There are, however, specific processes that change the 13 CO 2 exchange fluxes slightly differently from 12 CO 2 fluxes.We can write a global mass balance for atmospheric δ 13 C (δ a ) so that the different isotopic processes are explicitly defined and dependent on total CO 2 fluxes (for the derivation of Eq. 1 see Tans et al., 1993).We can then identify the (1) emission forcing terms, (2) net exchange isotope forcing terms, and (3) grossflux isodisequilibrium forcing terms: [emission forcing terms] where C a is the total carbon content (unit mol or mass) in the atmosphere (in the form of CO 2 ).The subscripts ba and oa denote the direction of the one-way gross fluxes (unit mol or mass per unit time).For example, F ba refers to the respiratory release of CO 2 from terrestrial biosphere to atmosphere.The isotopic ratios of 13 C/ 12 C are expressed as δ xx (‰), where the subscripts refer to the signature in biosphere vegetation and soils (b), in biomass burning flux (fire), or in the fossil fuel emission flux (ff).The signature δ eq a depicts the isotopic ratio of CO 2 that is in equilibrium with the ocean surface and δ eq b depicts the ratio in the terrestrial biosphere that would be in isotopic equilibrium with the current atmosphere, which is more depleted in 13 CO 2 than when the biomass was formed years ago.N b and N o refer to net exchange fluxes (gross release minus gross uptake) of CO 2 , and F ff and F fire are the fossil fuel and biomass burning CO 2 emissions, respectively.
The terrestrial (photosynthetic) isotopic discrimination in Eq. ( 1) is expressed as ph = δ eq b − δ a ≈ − ph (‰), and can be derived from a CO 2 gradient-weighted average of different isotope fractionation effects during the transfer of CO 2 molecules from the canopy air until their reaction with the enzyme Ribulose-1,5-bisphosphate (Rubisco) in the chloroplasts of the plant leaf.There are two main fractioning effects along this pathway: the plant fractionates with s = 4.4 ‰ when CO 2 diffuses from leaf boundary through leaf stomata, and with f = 28 ‰ during carboxylation.Smaller fractionation effects occur during diffusion between canopy air and leaf boundary ( b = 2.9 ‰), and during dissolution of CO 2 in mesophyll water ( diss = 1.1 ‰) and transport to chloroplasts ( aq = 0.7 ‰).The parameterization of ph for C 3 plants has been described by Farquhar et al. (1982) and takes the following form as in Suits et al. (2005) where c a, s, i, c represent CO 2 partial pressures in canopy air space, leaf boundary layer, stomatal cavity, and in the chloroplasts, respectively.The overall discrimination ph value reflects mostly the fractionation step with the highest resistivity (O'Leary, 1988).For example, during a drought when the leaf's stomatal conductance is lowered in an attempt to prevent evaporative water loss, the diffusive s is the most limiting factor resulting in a lower overall ph .The opposite happens under more favorable environmental conditions when stomatal aperture is higher and carboxylation is the limiting factor, resulting in a higher overall ph .
The overall discrimination leaves the atmosphere relatively enriched and plants relatively depleted in 13 C. C 3 I.R. van der Velde et al.: A CO 2 and δ 13 C data assimilation system plants are depleted in 13 C by approximately −20 ‰ relative to the atmosphere and C 4 by approximately −4 ‰ as they can assimilate 13 CO 2 more efficiently with Rubisco.C 4 photosynthesis is essentially a more complex form of carbon fixation than C 3 photosynthesis as it shields Rubisco in the bundle sheath cells from wastefully binding with oxygen rather than carbon dioxide.
In addition to discrimination effects during photosynthetic uptake, we also need to account for isotopic enrichment of the atmosphere through respiratory release of carbon with a heavier isotopic signature after spending from one year to several decades or more in the plant and soil organic matter.This respiratory part will still enrich the atmosphere with 13 CO 2 even if net CO 2 uptake equals zero (Ciais et al., 1995), and we refer to it as the terrestrial isodisequilibrium flux in Eq. ( 1).
Discrimination associated with the dissolution of CO 2 in ocean water (Zhang et al., 1995) is much smaller and spatiotemporally homogeneous ( ao = −2 ‰) than in the terrestrial biosphere.The difference between ocean and land discrimination provide an additional constraint on the net fluxes and has already been demonstrated in previous studies (e.g., Keeling et al., 1989;Tans et al., 1993;Ciais et al., 1995;Fung et al., 1997;Rayner et al., 2008).We also have to account for isotopic disequilibrium that exists between the atmosphere and oceans.This isodisequilibrium flux is associated with the outgassing of CO 2 from the ocean waters, and has globally an enriching tendency on the δ a signatures.
Besides the land and ocean discrimination and disequilibrium forcing terms we have two additional terms in Eq. (1).Firstly, there are CO 2 emissions due to combustion of fossil fuels, which have a distinct isotopic signature depending on the organic fuel type, but globally its signature is approximately δ ff = −30 ‰.Secondly, there are CO 2 emissions due to biomass burning, where δ fire bears the signature of the 13 CO 2 and 12 CO 2 fluxes of F fire , which is typically the signature of burnt leaf foliage, woody tissue, and the aboveground litter (van der Velde et al., 2014).

CTDAS-C13
We followed the method presented by Peters et al. (2005) for designing the joint CO 2 and δ a data assimilation system.The architecture is similar to the CarbonTracker Data Assimilation Shell (CTDAS) v1.0 discussed in detail by van der Laan-Luijkx et al. (2017).Just like the traditional CO 2 -only inversions, we aim to close the CO 2 budget through fluxes from fossil fuel combustion, biomass burning, and net exchange fluxes from the terrestrial biosphere and oceans.In addition, we also intend to simultaneously close the 13 CO 2 ( 13 C a ) budget using the same set of CO 2 fluxes.Isotopic signatures themselves are not conserved quantities, therefore we calculate conserved mole fractions of CO 2 and 13 CO 2 in our transport model, which we can sample at designated locations and times to calculate δ a .The combined set of balance equations (unit mol per unit time) takes the following form: After some manipulation of Eqs. ( 3) and (4), by following Tans et al. (1993), we obtain the following: The 13 C a balance equation is now a close analog of Eq. ( 1), because 13 C a is a function of discrimination, N b and N o , and isodisequilibrium fluxes.The R values depict the isotopic ratio of 13 CO 2 /CO 2 in the atmosphere (R a ), in fossil fuel (R ff ), and biomass burning emissions (R fire ), and their values are approximately 0.011.The isodisequilibrium fluxes from land and ocean surfaces are here simply shown as D b and D o , respectively.The term λ discr ph /1000 + 1 represents the optimized ratio between the isotopic signature in the photosynthetic flux and atmosphere (R ph /R a ), and ranges between 0.980 and 0.996 depending on the prior ph and discrimination scaler λ discr .The term ( ao /1000 + 1) represents the ocean flux ratio and is held constant at 0.998 assuming ao = −2 ‰, and is not optimized.The parameters λ b and λ o represent the linear scaling factors for each week and ecosystem region (ecoregion) to adjust the net carbon exchange over land and ocean surfaces, respectively.For land, the scaling factor is associated with one scalar per ecoregion based on the Olson (1985) land use classification following Peters et al. (2005Peters et al. ( , 2007) ) (Fig. 1).The terrestrial biosphere is further divided into 11 larger geographical areas also known as TransCom regions (Gurney et al., 2002).Like in the early CarbonTracker releases, each of the 11 TransCom land regions contains a maximum of 19 ecoregion types (Fig. 2) and the ocean is divided into 30 large basins encompassing largescale ocean circulation features.This gives a maximum of 239 (= 11 • 19 + 30) different scaling factors each week (Peters et al., 2007).The new parameter is λ discr , which is used to scale a maximum of 209 terrestrial discrimination parameters per week.They are associated with the same 1 • × 1 • ecoregions as the terrestrial fluxes.Note that the maximum number of scalable land parameters is in reality ∼ 130, and not 209, because not each land region contains all 19 ecoregion types.
The terrestrial net exchange term in Eq. (5) (λ b N b λ discr ph /1000 + 1 R a ) includes two multiplicative scaling factors, making the required solution nonlinear.This poses a potential problem where variations in net exchange and discrimination are canceling each other out to such a degree that it leads to low signal-to-noise ratios, especially for discrimination.This is further investigated in Sect.3.2.The fossil fuel combustion, biomass burning, and terrestrial and ocean isodisequilibrium fluxes all remain fixed a priori estimates.We describe in Sect.3.1 the tuning of the latter disequilibrium fluxes to close the long-term mean balance of δ 13 C in our system.
The scaling factors λ b , λ o , and λ discr are the unknowns that are combined in state vector x (with dimension s), for which we will try to find an optimal solution by minimizing a quadratic cost function.In this function, there is a balance between information drawn from the observation vector y (with dimension m) with a covariance R (m × m) and prior knowledge from the state vector x p (s) with a covariance P (s × s): The observation operator H (m) represents the atmospheric transport model that propagates the surface fluxes from Eqs. ( 3) and ( 5) and samples accordingly the mole fractions of CO 2 and 13 CO 2 at the same location and moment as the observations y.
I. R. van der Velde et al.: A CO 2 and δ 13 C data assimilation system The solution for x that minimizes J is (Tarantola, 2005) where K represents the Kalman gain matrix (Peters et al., 2005).Equation ( 7) can be expressed in terms of λ (posterior scaling factor), λ p (prior scaling factor), and separate measurements of CO 2 (c) and δ 13 C (δ) with dimensions (j) and (k), respectively: In state vectors x and x p , the scaling factors for terrestrial discrimination are appended after the flux scaling factors.Similarly in the observation vectors y and H (x p ), the δ 13 C observations are appended after the CO 2 observations.The K matrix determines how much a scaling factor needs to change given a set of CO 2 and δ 13 C measurements.The matrices R and P modulate whether observations or bottom-up estimates are given more weight to the solution.
The P matrix contains 448 × 448 elements in total and is shown in Fig. 3.The first 209 × 209 element block contains the land flux uncertainties per ecoregion and their spatial correlations.The second 30 × 30 element block contains the ocean flux uncertainties per ocean basin.We gave the land scalars and the ocean scalars a maximum uncertainty of 80 and 100 % along the diagonal, respectively, as in earlier CarbonTracker releases.The third 209 × 209 element block contains the terrestrial discrimination scalars with a maximum uncertainty of 20 % along the diagonal with an identical spatial correlation structure as applied to the terrestrial flux uncertainty scalars.This implies that we can scale ph by a factor of 1.0 ± 0.2, and thus for a typical C 3 plant ( ph = −20 ‰) the mean and uncertainty lies around −20 ± 4 ‰.Furthermore, there is covariation between ecoregions of nearby TransCom regions, e.g., between North America boreal and temperate regions and between Europe and Eurasian regions.We did not allow covariances between net exchange and discrimination in order to give the parameters enough freedom in the solution.
The covariance structure of R is similar to CO 2 -only CT-DAS, but is extended with additional uncertainties in δ 13 C observations.These expected uncertainties quantify our ability to simulate observations given the uncertainty in atmospheric transport modeling and measurement errors.Section 2.3.5 gives an overview of the used uncertainties for each observation category.
With this inversion framework in place CTDAS-C13 progresses in a similar manner as the traditional CO 2 -only CT-DAS.For each week, the set of unknowns in the state vector are updated in a cycle that contains two steps.First there is a forecast step, which is driven by our fluxes and current background state vector x p to forecast an ensemble of CO 2 and 13 CO 2 mole fractions 5 weeks ahead in time.This is followed by an analysis step to determine the new state of the system with Eq. ( 8) such that it is consistent with the observations for the current week of the cycle.The analyzed state is propagated to the next cycle using the same model as Peters et al. (2007, Eq. 1 of the Supplement), and with this new state a new cycle begins with another forecast step to forecast a new ensemble of the background state 5 weeks ahead in time, now with an additional set of observations from a new week.The ensemble for each tracer is created from 150 ensemble members to provide a Gaussian probability density function of the state vector.
The simulation of atmospheric transport is provided by the two-way nested global transport model TM5 release 3 (Krol et al., 2005).This application simulates the atmospheric transport of CO 2 and 13 CO 2 at a global 6 • × 4 • resolution, with no nesting.It is driven by 3 hourly meteorological output from ECMWF ERA-Interim reanalysis (Dee et al., 2011).All the CO 2 and 13 CO 2 flux fields provided to the model are in units of mol CO 2 m −2 s −1 and mol 13 CO 2 m −2 s −1 , respectively.Atmospheric concentrations of CO 2 and 13 CO 2 are calculated as mole fractions in mol mol −1 .Signatures of δ 13 C are computed to the relative per mil value using the following conversion formulation in order to facilitate comparison with observations: where R ref is the VPDB reference ratio adopted for 13 CO 2 / ( 12 CO 2 + 13 CO 2 ), which is 0.011112 (Tans et al., 1993).R is the ratio of simulated mole fractions 13 CO 2 /CO 2 .

Terrestrial biosphere fluxes
The terrestrial first-guess net CO 2 exchange (N b ) and fire (F fire ) estimates were calculated in the Simple-Biosphere Carnegie-Ames-Stanford Approach model (SiB-CASA; Schaefer et al., 2008) on a 1 • × 1 • grid on a 10 min time resolution and were further processed into 3 hourly mean fluxes to serve as input for CTDAS-C13.SiBCASA is a biogeochemical model that calculates carbon, isotope, water, and energy exchange fluxes.The model inherited the aerodynamic and surface resistance models from SiB (Sellers et al., 1996) to solve for CO 2 partial pressures in an iterative loop to acquire a balance between net assimilation rate, stomatal conductance, mesophyll conductance, and CO 2 partial pressures.The aerodynamic resistance model describes the turbulent transfer processes using the Monin-Obuhkov similarity theory.The surface and interior resistance models describe the pathway of CO 2 (but also water and heat) through the leaf boundary, the leaf stomata, and ultimately the leaf chloroplasts.The Ball-Berry-Collatz model is used to estimate stomatal conductance (Ball, 1988;Collatz et al., 1991) and is coupled to the Farquhar and Collatz photosynthesis models for C 3 and C 4 vegetation (Farquhar et al., 1980;Collatz et al., 1992).A mesophyll conductance formulation was introduced by Suits et al. (2005) to predict realistic CO 2 partial pressures in the chloroplasts.Mesophyll conductance is suggested to be as important as stomatal conductance in terms of magnitude and variability, and it is shown that δ 13 C correlate more precisely with c c /c a than with c i /c a (Flexas et al., 2008).It is parameterized as a function of the canopy photosynthetic rate and soil water stress factor.SiBCASA is driven by 3 hourly ECMWF ERA-Interim meteorology, designed with a semi-prognostic leaf pool to track seasonal plant phenology, and it uses GFED4 daily burned area disturbances to calculate fire fluxes at a fine temporal resolution (van der Velde et al., 2013Velde et al., , 2014)).The model incorporates 12 different aggregated ecosystems according to Olson (1985)  from the plant and soil is calculated in the CASA part of the model using 13 biogeochemical pools with environmentinfluenced turnover rates (Schaefer et al., 2008).

Ocean fluxes
The ocean first-guess net CO 2 exchange (N o ) estimates derive from ocean inversions from Jacobson et al. (2007).These long-term estimates are combined with the quadratic gastransfer velocity from 3 hourly ECMWF ERA-Interim wind fields (Wanninkhof, 1992) to create fluxes on a 1 • × 1 • grid at a 3 hourly temporal resolution.An additional trend was applied to the fluxes to ensure that increases in anthropogenic uptake are proportional to increases in atmospheric CO 2 levels (see http://www.esrl.noaa.gov/gmd/ccgg/carbontracker).

Isotope and disequilibrium fluxes
To calculate the fluxes of 13 CO 2 from land surfaces we used the photosynthetic discrimination parameterization (Eq.2) for C 3 plants in the SiBCASA model (van der Velde et al., 2014).The weighted leaf-level value for C 3 discrimination is typically 19.0 ‰, and given the more efficient CO 2 bonding with the Rubisco enzyme, C 4 discrimination is 4.4 ‰ (Still et al., 2003;Suits et al., 2005).Given the dominance of C 3 plant growth (70 % of global GPP), the global mean discrimination in SiBCASA has been estimated at ph = 15.2 ‰.SiBCASA's spatial heterogeneity of land discrimination is shown in Fig. 4. It reflects the land use distribution and the environmental forcing.Large discrimination values can be found in the temperate regions, the boreal forests, and in the humid environments such as the tropical rain forests in South America, Africa, and South East Asia.Small discrimination values can be found in the United States corn belt and in the dry climate regions such as the African savannas and Australian grasslands, where there is an abundance of C 4 plant growth.More subtle variations in ph in C 3 dominant regions are driven by differences in environmental conditions (e.g., humidity, groundwater availability, and light intensity).Weekly 1 • × 1 • fields for ph were used to map the regular 3 hourly net CO 2 fluxes to 13 CO 2 fluxes: where ph is derived from SiBCASA's ph output.Their relation is straightforward For the calculation of 13 CO 2 biomass burning flux we assumed R fire to be very close to the signature of newly assimilated photosynthates: 13 F fire = F fire λ discr ph /1000 + 1 R a . (12) The 13 CO 2 fossil fuel emissions are calculated with R ff = 0.0107786, given that the global mean value of δ ff is equal to Note that we did not vary δ ff for different fuel types in this version of CTDAS-C13, but such variability could be included in the future based on the work of Andres et al. (2000).
The ocean discrimination parameter ao is assumed to be constant at −2 ‰ as in many comparable studies (e.g., Tans et al., 1993;Ciais et al., 1995;Alden et al., 2010), and is not optimized.The regular 3 hourly net CO 2 fluxes were mapped to 13 CO 2 fluxes as follows: The isodisequilibrium fluxes (D b and D o , in mol 13 CO 2 m −2 s −1 ) were made available on a monthly 1 • × 1 • resolution.D b is calculated using SiBCASA's gross natural respiratory flux scaled with isotopic disequilibrium of the terrestrial biosphere with the current atmosphere, i.e., F ba δ b − δ eq b .Because fossil fuel emissions add isotopically depleted CO 2 to the atmosphere, the biosphere signature δ b follows with a time lag dependent on the residence time of carbon in the vegetation and soils.That implies δ b is larger than δ eq b , which is the biosphere signature that is in equilibrium with the current atmosphere (Tans et al., 1993).D b has a positive tendency on atmospheric δ 13 C as carbon originating from different SiBCASA pools is older and more enriched in 13 C than the isotopic signature of recently fixed photosynthates.The SiBCASA pool configuration is described in detail in van der Velde et al. (2014).
D o is calculated from the outgassing flux of CO 2 scaled with the isotopic disequilibrium of the ocean surface with the current atmosphere, i.e., F oa δ eq a − δ a .The δ eq a term is determined from a global network of δ 13 C measurements in dissolved inorganic carbon (Gruber et al., 1999).F oa is parameterized as a function of surface ocean partial pressure of CO 2 and wind speed after Takahashi et al. (2009).Wind speed and solubility are assumed to remain constant year-toyear.The disequilibrium fluxes are positive from the equator to approximately 60 • of latitude in both directions and are negative beyond that.

Observations
Observations of CO 2 from a wide range of research laboratories are bundled in Observation Package (ObsPack) version 1.0.products that include the provider's original data and metadata reformatted into the ObsPack framework (Masarie et al., 2014).
From the available CO 2 observations, approximately 24 000 weekly flask measurements were used in the assimilation from a fixed network of 58 surface sites.Another large set of 174 000 measurements came from 23 semi-continuous in situ sites.Most CO 2 measurements are obtained with a nominal precision of ±0.1 ppm.The remainder of sites and measurements (including from aircraft or shipboard) were not used because of double records, and some measurements were kept for independent checks.A small fraction was omitted as our model could not resolve certain locations at a coarse resolution.
For the dual-species inversions we also used 22 000 flask measurements of δ 13 C from 53 different surface sites.A further 5600 measurements from five different sites were obtained using programmable flask packages, which measure δ 13 C at a daily resolution.The isotope ratios are measured by dual inlet mass spectrometry with a precision of ±0.01 ‰.
We determined observation uncertainties (model-data mismatch or MDM) for each of the δ 13 C measurement sites in a heuristic manner based on earlier test inversions.These values are added to the diagonal of R. A too small error would give an unrealistic amount of confidence how well the model is expected to represent the measurement location during sampling but a too large error would give very little confidence to the measurement representation.The δ 13 C measurement sites were divided into different categories each with their own MDM value.As with CO 2 , these categories were land, mixed conditions, marine boundary layer, deep Southern Hemisphere, and a special category for problem sites where forecast performance is poor.For each site, we determined the innovation statistic χ 2 , which is a measure for how apt our applied uncertainty level is given the model-data fit.A χ 2 value of 1.0 indicates that the simulated and expected total uncertainty are equal, lower values indicate overestimation of the uncertainty and higher values underestimation.Table 1 gives a summary of the site categories used, together with the assigned MDM for δ 13 C and the category-average innovation χ 2 determined from an inversion experiment.For the majority of sites the innovation values are between 0.7 and 1.3, i.e., around the ideal value of Table 2. Summary of the four inversion experiments, the observations used, the optimized items (ocean and land fluxes, and ph ), and their linearity.The prefix TRAD− refers to traditional, i.e., experiments that have been performed in the past in any way, shape or form.The prefix NEW− refers to a new type of inversions used in this publication.NEW−CO2C13 used the default CTDAS-C13 model setup as described in the Methodology, while NEW−2STEP solved for ph using only δ 13 C data.

Experiment
Observations Optimization Linear?1.0.For the CO 2 measurement sites we used a similar set of MDM values as in previous CarbonTracker releases.

Experiments
We performed four inversion experiments as summarized in Table 2.The simulation period covered the years 2000 through 2011, but our analyses focused on the period 2001-2011, i.e., we omitted the spin-up year.As a benchmark, we performed a traditional inversion to estimate the net carbon exchange fluxes of the ocean and land using only CO 2 observations, which we call TRAD−CO2.For the second inversion, we added δ 13 C observations alongside CO 2 to constrain only the exchange fluxes; therefore, we call this experiment TRAD−CO2C13.The experiment in which we estimated discrimination and fluxes simultaneously is called NEW−CO2C13.This inversion is nonlinear because the discrimination scaling parameter is in the same multiplication term as the net flux scaling parameter.The fourth experiment was a linear inversion experiment where we estimated only the land discrimination parameter using δ 13 C data.We call this experiment NEW−2STEP because discrimination was solved in a second step after optimization of the net exchange fluxes.This means that ocean and land fluxes were derived from the optimized state vector and its covariance from the TRAD−CO2 inversion.

Comparison to observations of CO 2 and δ 13 C from the global network
We first evaluate the global CO 2 and δ 13 C budgets simulated by our combination of fluxes as described in Sect. 2 to assess where we expect the largest changes in the optimization.As shown in Fig. 5a, the prior net exchange flux estimates and unscaled disequilibrium fluxes were not large enough to close the gap with the observed tracers, CO 2 , and δ 13 C.The sum of the flux arrows overestimated the annual CO 2 growth rate along the x axis and overestimated δ 13 C depletion along the y axis.In a traditional TRAD−CO2 inversion, the estimated ocean and land fluxes closed the CO 2 budget along the x axis.The leverage in the net exchange fluxes was, however, not large enough to close the δ 13 C budget along the y axis as well.In an inversion that includes δ 13 C observations, the gap in δ 13 C would adjust the CO 2 flux magnitudes and ocean/land partitioning to unrealistic magnitudes in an effort to overcome the large offset between the simulated and observed δ 13 C growth rate.Instead we chose to use scaled disequilibrium fluxes in our inversions in order to estimate land and ocean CO 2 flux magnitudes that remain close to the results of other traditional carbon cycle budgeting studies (Alden et al., 2010;van der Velde et al., 2013).
We chose the disequilibrium fluxes to adjust because (1) the exact magnitudes of these terms are still unknown due to uncertainties in the carbon pool turnover, gross carbon fluxes, and isotopic discrimination, and (2) these terms do not affect the CO 2 mass balance.It assured a closed mean δ 13 C budget of our inversions without creating unrealistic carbon sinks over land and oceans (Fig. 5b).Most importantly, closing the climatological (11 year) budget allowed us to focus our study on interannual changes in the net fluxes and photosynthetic discrimination.
We obtained the best fit with δ 13 C data when the land and ocean disequilibrium flux were scaled by a factor of 1.2 without changing either their spatial patterns or time trends.This is consistent with recent double deconvolution studies where the global δ 13 C balance was closed with a factor of 1.3 in land and ocean disequilibrium (Alden et al., 2010).Our value was determined after assessing an ensemble of different sets of scaling numbers (ranging from 1.1 to 1.5) in a forward TM5 simulation, which was driven by the optimized net land and ocean flux estimates from the TRAD−CO2 experiment.This assured a closed multi-year δ 13 C budget together with a closed multi-year CO 2 budget.As selection criteria, we used (1) the 11 year mean root mean square difference (RMSD) of a large selection of δ 13 C sites and (2) the average bias between simulated and observed values.In the non-scaled disequilibrium simulation we obtained a RMSD of 0.165 ‰ and a bias of −0.110 ‰ averaged over all sites.The optimal result was obtained with a scaling factor of 1.2, which reduced the RMSD to 0.079 ‰ and the mean bias to −0.010 ‰.Note that these scaling factors cannot be applied to other inversion studies because the disequilibrium scaling factors are tuned for this particular system and time period.
To demonstrate our procedure in terms of individual data sets, we refer to Fig. 6.After scaling the disequilibrium fluxes and using optimized net carbon exchange from the TRAD−CO2 inversion, time series of δ 13 C at 32 of the 46 Northern Hemisphere sites showed no remaining significant trend (sites where p value > 0.05) in the summer residuals, and the residuals from the trend lines were within or close to the MDM specified for our dual-species inversions.Some of the sites with remaining trends are located at great distances from large continental carbon sources and sinks, and exert little influence on the posterior λ discr parameter (e.g., CHR, GMI).Some of the other sites were assigned a large MDM (e.g., BAL, NWR, TAP) giving them less weight in the estimation of the posterior λ discr parameter.The collection of sites with remaining trends do not seem to have a systematic geographic pattern and are likely reflecting a change in local oceanic or biospheric isotope exchange, such as must be the case for the Bermuda West (BMW, non-significant positive trend) and Bermuda East (BME, significant downward trend) site.
With the long-term trend of δ 13 C appropriately captured, we proceeded to optimize NEE and ph with our new framework (NEW−CO2C13).We show that this inversion further reduced δ 13 C residuals (Fig. 7a), without compromising (nor strongly improving) the fit to CO 2 (Fig. 7b) that we attained from the TRAD−CO2 inversion.In Fig. 7a the ratio of δ 13 C RMSD of NEW−CO2C13 to δ 13 C RMSD of the TRAD−CO2 inversions was at most sites smaller or equal to 0.95 (indicating a significantly higher accuracy of NEW−CO2C13 in form of bias and noise reduction): In Fig. 7b, the ratio of CO 2 RMSD of NEW−CO2C13 to CO 2 RMSD of TRAD−CO2 was at most locations between 0.95 and 1.05.This suggests that the two atmospheric constraints applied are complementary, and there is no indication that the TRAD−CO2 results from CarbonTracker were inconsistent with δ 13 C measurements.This is an important prerequisite for a credible estimate of discrimination in our system.Furthermore, Fig. 7a shows a notable latitudinal divide in the reduction of δ 13 C RMSD, indicating the utility of NEW−CO2C13 in the Northern Hemisphere due to the large availability of measurements and scalable discrimination parameters.
At sites like Alert (Nunavut, Canada) the NEW−CO2C13 inversion provided a better fit to the measured data than the TRAD−CO2 inversion (Fig. 8).The 11 year averaged δ 13 C residuals were close to zero for both inversions, as the disequilibrium flux was tuned specifically to prevent large residuals in a-priori simulated δ 13 C as described in Sect.3.1.The 1σ SD of the δ 13 C residuals at Alert were smaller in the NEW−CO2C13 inversion in comparison to TRAD−CO2 due to the additional optimization of ph alongside net exchange fluxes.The CO 2 residuals for Alert in Fig. 8 were for both inversions almost identical.

Linear and nonlinear estimates of net carbon uptake and land discrimination
Simultaneously optimizing both λ discr and λ bio is inherently nonlinear and thus possibly problematic for our assimilation system; therefore, we tested the validity of our approach in the NEW−CO2C13 inversion.We hypothesized that a region's net carbon uptake and discrimination would change in a similar fashion in the nonlinear inversion, as it would for a linear inversion.The linear inversion experiment consisted of two consecutive steps: (1) the optimization of the net exchange fluxes using only CO 2 observations (TRAD−CO2) followed by (2) the estimation of the land discrimination parameter using only δ 13 C observations (NEW−2STEP).In the nonlinear NEW−CO2C13 inversion the optimization of fluxes and discrimination was done simultaneously.For net carbon uptake by vegetation, we refer to net ecosystem exchange or NEE defined as positive when CO 2 is taken up from the atmosphere.For plant isotope discrimination we refer to ph in per mil, which is defined as positive.
As shown in Fig. 9, the 11 year mean NEE for the 11 land TransCom regions are very similar in the nonlinear NEW−CO2C13 and linear TRAD−CO2 inversions.Deviations are in the order of tens of teragrams, and within 1σ SD of the flux interannual variability (IAV).Figure 9  gregated ph values.In the boreal regions, where there is very little C 4 plant growth, the discrimination is at its maximum (approximately 20 ‰, 5 ‰ above the global average), but in regions where there is C 4 plant growth (e.g., due to agriculture in the United States or savannas in Africa) the mean ph values are lower (approximately 12-15 ‰).These regional patterns of ph imposed by SiBCASA (see also Fig. 4) are maintained by the NEW−CO2C13 inversion framework.Because we aimed to retrieve robust temporal patterns of IAV, the most relevant indicators for the robustness of our nonlinear inversion approach are given by the correlation coefficients (r) between the two types of inversions.We calculated r for NEE and ph between the linear and nonlinear inversions.As the seasonal cycles in uptake and discrimination are largely dictated by the prior estimates, we removed them using a 3 month boxcar mean smooth curve fitting to obtain the anomalies relative to the seasonal trend.
The NEE in NEW−CO2C13 is very similar to the NEE in TRAD−CO2, as indicated by the high r values (> 0.96 for N = 52 • 11 weeks) for all TransCom regions.The r values are lower for ph , but still exceed 0.75 in the Northern Hemi-sphere.The correlation is particularly high over North America boreal, North America temperate, and European regions.Smaller correlations are obtained in tropical and temperate South America and tropical Asia.This is expected, however, as these regions typically suffer from a lack of observational constraints.
The linear NEW−2STEP inversion estimated the same large increase in discrimination IAV as in the nonlinear NEW−CO2C13 inversion for the Northern Hemisphere in comparison to the first-guess estimate of SiBCASA (8-fold increase in SD, see Table 3).In addition, we also found in both inversions a strong positive correlation between ph and NEE on annual time scales (r = 0.79) with a significant slope (p = 0.001, 95 % confidence interval of a two-sided distribution with 9 degrees of freedom).In years when annual mean NEE is low (less carbon uptake), the ph is low too (less discrimination), implying that stomata have partially closed, and vice versa.This correlation did not emerge in the TRAD−CO2 estimate based on atmospheric CO 2 observations alone, and it also did not emerge if δ 13 C observations were additionally used in the TRAD−CO2C13 esti-Figure 6. Summer (JJA) residuals of δ 13 C (‰) in CO 2 for 46 sites (excluding aircraft and ships) situated in the Northern Hemisphere.These sites are ordered based on the their latitudinal location; most northern site is placed at the top left (Alert, Canada) and the site nearest to the equator at the bottom left in the second part of Fig. 6 (Christmas Island, Republic of Kiribati).All residuals (simulated minus observed) are calculated from a traditional TRAD−CO2 inversion with scaled disequilibrium fluxes.Assuming a closed long-term mean budget in δ 13 C, we tested the null hypothesis that the slope of the linear regression line is zero.Sites with a trend where the p value is smaller than the significance level of 5 % are shown in red, whereas the remaining sites without significant trend are shown in green.The sample uncertainty (model-data mismatch) used for the NEW−CO2C13 and NEW−2STEP inversions is displayed by transparent gray areas.Sites marked with ** were not included in the inversions but were used for independent verification.For detailed information of the sites and their locations, we refer to the NOAA website: http://www.esrl.noaa.gov/gmd/ccgg/carbontracker/observations.php.ph and the correlation with NEE were driven by δ 13 C observations, and were not a symptom of the systems inability to separately estimate NEE and ph variations.This suggests that the estimated IAV of ph in the nonlinear inversion is truly a signal retrieved from δ 13 C that would otherwise be aliased erroneously into the carbon fluxes or not retrieved at all.

Independent verification with drought indices
A closer inspection reveals that the reported correlation between the Northern Hemisphere's NEE and ph in Table 3 could indicate a moisture-driven response on the ecosystem level.We identified several moments of severe to extreme drought as characterized by the Standardized Precipitation and Evaporation Index (SPEI; Vicente-Serrano et al., 2010) below −1.0 that covered an extensive area of more than a million square kilometers in the States.These droughts are described in literature as the droughts (or heat waves) of summer 2002 (Seager, 2010;Schwalm et al., 2012) and 2011 (Long et al., 2013).The annual averaged maps of SPEI for 2001-2011 are shown in the top panel of Fig. 10 calculated for the Northern American temperate TransCom domain.Independent of the SPEI drought index, we estimated changes in ph and NEE over the same American domain with the NEW−CO2C13 inversion using atmospheric CO 2 and δ 13 C data (Fig. 10 middle and lower panels).A correlation between ph and SPEI could only be established by applying an area weighting function to the SPEI index to give years that experienced large and severe droughts, the strongest association with reductions in ph .We used the following function for the Weighted Drought Index (WDI): Total-area .
In words, we sum over the product of the SPEI index and the grid cell surface area where SPEI is below −1.0 and subsequently we divide it by the total area of the TransCom domain.Hence, the WDI is an expression of the drought in terms of the surface area that is affected.A larger drought surface area will result in a more negative WDI.Using this function, we see that the lower values for ph correspond strongly with years of low SPEI over large serried areas, indicating a temporal correlation of r = +0.75 between the SPEI variable and ph (see correlation in Fig. 11 with a significant slope: p = 0.008, 95 % confidence interval of a twosided distribution with 9 degrees of freedom).The two largest anomalies (> 1σ of 11 year IAV) in annual mean ph correspond with low SPEI in 2002 and in 2011.A third notable drought as recorded in SPEI happened in 2006; although carbon uptake was reduced, it did not amount to a significant signal in ph .Similar correlations do exist over other parts of the Northern Hemisphere in our inversion solution.For instance, severe droughts in western Europe (2003) and Russia (2010) lowered the discrimination by 1.0 ‰, and exceeded more than 1σ SD of its 11 year IAV (not shown).
In addition, in years when ph is low, the annual mean NEE tends to be low too, possibly as a result of reduced GPP.This implies that leaf stomata have partially closed and are therefore affecting both ph and carbon uptake from photosynthesis.The reduction of the optimized net carbon sink for North America is 100-400 Tg C yr −1 during the drought years of 2002, 2006, and 2011 (in comparison to their surrounding years).
These correlations that are averaged over continent sized areas do, however, breakdown on smaller scales.On regional scales, we observed a partial misallocation of the model adjustments of NEE and ph in comparison to SPEI.This is largely a consequence of our limited capacity to monitor CO 2 and δ 13 C.For example, for North America temperate 2002, where the drought index was negative over the mountain states, the impact on the carbon cycle was strongest over the eastern forests of the United States.In these forest ecosystems, CO 2 exchange is much stronger than over the mountains, and hence their impact on atmospheric δ 13 C as well.
Notice that the prior net carbon sink is underestimated in comparison to the optimization because SiBCASA assumes a near steady state between GPP and TER (Fig. 10).SiBCASA was in fact able to simulate small carbon uptake anomalies during the reported droughts using its own environmental response parameterizations.However, it lacked a substantial amount of interannual variability in NEE and ph and lacked a strong correlation of ph with SPEI (Fig. 11).This suggests a potential absence of an important coupling between the hydrology and carbon discrimination processes in the model.

Discussion and conclusions
We developed a new application of the CarbonTracker data assimilation system that simulates two atmospheric tracers simultaneously: CO 2 and the δ 13 C isotope signature of CO 2 .We used measurements of both tracers to optimize the net ocean and land carbon exchange fluxes and the land discrimination parameter ph .The annual reductions in ph were up 0.75 ‰ and exceeded the 1σ SD of the IAV over 11 years in the North American domain (16.4 ± 0.3 ‰).We interpret these negative anomalies in ph as possible reductions of the intercellular CO 2 levels and relative increases in the intercellular 13 CO 2 / 12 CO 2 ratio, resulting from stomatal closure due to drought stress on the leaf level.This is the most plausible explanation as most other factors that affect ph are (a) included a priori in the SiBCASA biosphere model, such as    Seager, 2010;Schwalm et al., 2012;Long et al., 2013).These droughts correlate with reductions in annual mean ph , and reductions in the estimated carbon sinks as reported in Peters et al. (2007).) show no significant correlation between ph and large-scale droughts, while the simultaneous optimization of carbon sinks and ph with atmospheric CO 2 and δ 13 C observations (red triangles) suggests a highly significant correlation can be derived.slope of the red regression line is 1.61 ‰/WDI (p = 0.008, 95% confidence interval of a two-sided distribution with 9 degrees of freedom).The SiBCASA slope is, however, not significantly different from zero (p 0.05).The integrated ph values are GPP-weighted per grid box as in Fig. 10.WDI is based on the SPEI index but area weighted to give years with large serried areas that experienced severe droughts (with SPEI smaller than −1.2) more leverage.
the effects of IAV on the strength of photosynthesis over C 3 and C 4 vegetation; (b) the variations in mesophyll conductance are not expected to vary much from year-to-year, such as ecosystem composition; or (c) the variations in mesophyll conductance would enhance the intercellular CO 2 levels (and thus ph ) rather than reduce it, such as increased radiance of the leaves under reduced cloud cover.This suggests the possibility that the impact of environmental stress on stomatal conductance and carbon uptake is much larger than currently simulated by the widely used drought parameterizations in terrestrial biosphere models.These parameterizations are often derived from laboratory observations or plot-scale observations that often aggregate poorly over much larger scales.
Our first results suggest that a data assimilation system that uses the global atmospheric δ 13 C record, in concert with the CO 2 record, can offer new insights on large-scale drought dynamics of the coupled vegetation-atmosphere system.It is unlikely that our terrestrial biosphere model will reproduce the new large-scale atmospheric constraints on NEE and ph with a simple adjustment of the currently used drought response parameterizations (such as stomatal conductance and soil water stress inhibition functions).We experimented with a different stomatal conductance model based on vapor pressure deficit (VPD; Leuning, 1995) rather than relative humidity as it was shown to better predict changes in the isotopic composition in tree rings (Ballantyne et al., 2010).This modification, however, did not change the annual covariation between NEE and ph in SiBCASA.In addition, modifications in the soil water stress function of SiBCASA, which impacts ph through mesophyll conductance (Seibt et al., 2008) also had little impact on annual variations in ph .Instead, SiBCASA shows minimal dynamic range in the hydrological drivers of drought stress.This was also concluded using satellite-observed soil moisture in SiB-CASA over boreal EurAsia (van der Molen et al., 2016).That means our model is potentially (a) too homogenous regarding its plant and soil characteristics; (b) suffering from a too simple hydrological formulation for runoff and interception; (c) lacking realism in simulating the latency of ecosystem recovery after a severe drought; (d) misrepresenting effects of root-zone soil-moisture stress, as was also diagnosed for the Amazon by Harper et al. (2010) for the closely related SiB model; or (e) suffering from an even more fundamental problem inside the A-g s model where c i /c a and c c /c a are calculated.In SiBCASA, the soil moisture limitations are applied by first downscaling assimilation rate (A), V max , and mesophyll conductance, after which the balance is calculated between A, stomatal conductance (g s ), and c c /c a .In Egea et al. (2011) this approach was shown to conserve incorrectly intrinsic water-use efficiency (iWUE = A/g s ) and ph during droughts.However, initial tests show that a direct coupling of soil moisture stress to g s would affect SiBCASA's iWUE (and ph ) much stronger and more favorable during droughts (E.van Schaik, personal communication, 2017) There is also evidence that the conventional use of land cover types in biosphere models does not adequately describe the spatial variations in carbon exchange (Bloom et al., 2016).
As with any data assimilation system, the number of available observations largely determines the assimilation system's ability to retrieve meaningful signals.Our current method relies on atmospheric δ 13 C anomalies that affect multiple monitoring sites at the same time due to low signalto-noise ratios at each site, but the network coverage over many parts of the world is still sparse.The increase in the number of measurement sites, the addition of δ 13 C to many existing ones, particularly in sparsely populated areas could benefit CTDAS-C13 greatly.New measurement efforts are currently underway to improve our observational coverage in these sparsely sampled areas.Regular measurements of CO 2 from aircraft vertical profiles have recently commenced at four different sites above the Amazon.These data have provided new insights on the carbon cycle under drought conditions (Gatti et al., 2014).These new measurements were successfully used in an application of CTDAS (van der Laan-Luijkx et al., 2015) and confirmed that the Amazonian CO 2 uptake by vegetation was indeed reduced during the severe 2010 drought.Furthermore, some coauthors are currently involved in a new collaborative effort to provide the first high-precision measurements of δ 13 C and other isotopes in CO 2 from a large number of air samples collected over the Amazon basin.Using an assimilation system similar to that described here, these data would bolster our ability to quantify seasonal to interannual changes in the Amazonian carbon balance and better understand the influence of drought stress on NEE.
The retrieved correlation between NEE and ph in the Northern Hemisphere was derived from atmospheric δ 13 C observations through our new dual-species approach, and thereby provided new insights on the land-atmosphere coupling of water and carbon on continental and hemispheric scales.The unconstrained SiBCASA model does not show a large enough response to drought both in terms of NEE and ph .The correlation between droughts and ph over the North American temperate domain (Fig. 10) can only be demonstrated after optimizing NEE and ph by applying atmospheric δ 13 C and CO 2 constraints together.We emphasize that the reported correlations remain robust and significant even when changing the atmospheric transport characteristics (i.e., convection fields from ECMWF ERA-Interim meteorology vs. default TM5 convection scheme), the optimization method (nonlinear vs. linear 2-step), and when changing the assumed model-data errors of our data assimilation system.
A potential problem with the current framework is that we cannot account for changes in the terrestrial isodisequilibrium flux.In Eq. ( 1), we forced all missing isotopic vari-ability into term N b ph without considering additional variability from the isodisequilibrium term.Photosynthetic discrimination is also responsible for a portion of the variability in the terrestrial isodisequilibrium flux (van der Velde et al., 2013), but the extent is hard to quantify.The δ eq b signature (i.e., the biosphere signature that is in equilibrium with the atmosphere) is a function of the current δ a and ph , two quantities that ultimately exert influence on δ b as the isotopic signal carries through the series of carbon reservoirs (i.e., leaves, stems, roots, and ultimately the soils).The absence of direct adjustments to the disequilibrium flux could mean we aliased erroneously isotopic signals only onto the net flux term of the budget.In light of recent observational evidence, the variability in the disequilibrium term might be of more importance than recently thought.Bowling et al. (2014) showed with δ a measurements that the disequilibrium flux can become negative locally due to humidity-induced changes in ph .We found that these effects on the ph estimate are likely small mainly because gross flux variability is fundamentally limited and dampened by the large reservoir sizes from which it comes.In an experiment with SiBCASA where we allowed extra uncertainty in respiration and ph to drive through the disequilibrium isoflux we indeed found an increase in variability in the disequilibrium term, necessitating 10 % less ph variability to keep a closed δ 13 C budget.This indicates that allowing for errors in the disequilibrium fluxes, the variations in the estimated ph parameter might be slightly smaller or larger than estimated with CTDAS-C13, but nonetheless still twice as large than estimated with SiBCASA.On one hand, using a more simplified but physically consistent set of equations only based on gross fluxes (GPP and TER) to express the rate of change of δ a would eliminate the need for a disequilibrium term.On the other hand, this would complicate the closing of the CO 2 budget as it necessitates a way to effectively separate these two gross fluxes.
It is worth mentioning that the carbon residence time in land ecosystems is highly uncertain, and therefore the gross CO 2 exchange as well.Welp et al. (2011) suggested that the current popular estimate of global GPP of 120 Pg C yr −1 , which is also predicted by SiBCASA, may be a lower limit and could in reality be as large as 175 Pg C yr −1 to reflect faster turnover of carbon in the vegetation and soils.Such uncertainties were also underlined by Carvalhais et al. (2014) who found that higher precipitation rates are associated with faster carbon turnover, but that global modeled turnover is in fact often underestimated.We make a cautious conjecture that if GPP is in fact as large as claimed by Welp et al. (2011), and heterotrophic respiration is large too, it will partly explain the current underestimation in the modeled disequilibrium fluxes, which are a function of TER and ocean CO 2 outgassing.In this study, we closed the gap with a predetermined scaling factor of 1.2 on the disequilibrium fluxes for oceans and land without assuming actual changes in GPP, TER, or ph .We could, therefore, benefit from a more integrated as-Geosci.Model Dev., 11, 283-304, 2018 www.geosci-model-dev.net/11/283/2018/similation system where we are using atmospheric data to simultaneously optimize for terrestrial model parameters that exert an influence on GPP, TER, and carbon turnover.The CTDAS modular design (van der Laan-Luijkx et al., 2017) makes it now more straightforward to develop and implement such additional improvements.
To conclude, this study showed there is significant potential to use atmospheric CO 2 and δ 13 C data as constraints on plant NEE and isotopic discrimination using a dual-species assimilation platform.Signals that would otherwise be lost in a single tracer data assimilation system, such as the possibility of a drought-driven covariation between isotope discrimination and NEE or the separation of GPP from NEE, can potentially be detected in the described dual-species application of CTDAS.Continued and additional measurements of atmospheric δ 13 C and CO 2 , especially in future assimilation systems where biosphere model parameters are directly optimized, should help us better understand the hydrological and biogeochemical interactions between the atmosphere and vegetation.

Figure 3 .
Figure 3.The prior P covariance structure represents squared uncertainty of the dimensionless state vector.The first 209 × 209 element block represents the covariance matrix for land NEE with a maximum diagonal uncertainty of 0.64 (equivalent to 80 %), the second 30 × 30 element block represents the covariance matrix for ocean fluxes with a maximum diagonal uncertainty of 1.0 (equivalent to 100 %), and the third 209 × 209 element block represents the covariance matrix for ph with a maximum diagonal uncertainty of 0.04 (equivalent to 20 %).The matrix is organized according to TransCom ocean basins and land regions, where each land region contains 19 potential ecoregions (see Figs. 1 and 2).

Figure 4 .
Figure 4. Mean (2001-2011) modeled discrimination parameter ph (‰) from SiBCASA.The discrimination is more detailed for ph > 16 ‰ to highlight the more subtle variations in ph in the dominant C 3 regions that experience different environmental forcing.

Figure 5 .
Figure 5. Annual mean carbon (x axis) and δ 13 C (y axis) growth rates for (a) the prior estimates and for (b) the NEW−CO2C13 experiment.Colored arrows represent the different sources and sinks of the carbon cycle.A closed budget for both tracers was accomplished in the NEW−CO2C13 experiment, as indicated by the resultant vector (sum of all colored arrows) returning to the black arrow (observed growth rate in atmosphere).To close the long-term trend we increased the isodisequilibrium fluxes by 20 %.

Figure 7 .
Figure 7.A comparison of the relative performance of inversion techniques for the period 2001 through 2006 based on the ratio of the model-data (a) δ 13 C root mean square difference (RMSD) of NEW−CO2C13 to δ 13 C RMSD of TRAD−CO2, and (b) CO 2 RMSD of NEW−CO2C13 to CO 2 RMSD of the TRAD−CO2 inversion.A ratio lower than 1.0 indicates a higher accuracy of the NEW−CO2C13 inversion technique; green sites indicate a ratio less than or equal to 0.95, red sites indicate a ratio greater than or equal to 1.05, and sites where the difference in respective RMSD is less than 0.05 are given in black.The size of each symbol is a measure of the relative performance of NEW−CO2C13 in comparison to TRAD−CO2.The larger the symbol the more the ratio of RMSD differs from 1.0.

Figure 8 .
Figure 8.Comparison of two different inversion experiments at Alert (ALT, Canada).The top panel displays δ 13 C observations (black circles) together with simulated δ 13 C from NEW−CO2C13 (green circles).The top right panel displays the probability density functions (PDFs) of the residuals between NEW−CO2C13 and observed (green) and between TRAD−CO2 and observed (blue).The lower panel displays independent flask measurements (not used in the assimilation) of CO 2 (black circles) at Alert with simulated CO 2 from NEW−CO2C13 (green circles).Notice the almost identical distribution of the residual PDFs between NEW−CO2C13 and TRAD−CO2 inversion techniques.

Figure 9 .
Figure 9. (a) the 11 year mean land carbon uptake (Pg C yr −1 ) for each TransCom region with estimates from the nonlinear NEW−CO2C13 inversion (red) and estimates from the linear TRAD−CO2 inversion (yellow).Error bars depict 1σ SD of the flux IAV.The 11 year correlation coefficients r between the two inversion methods are given on top of the bars.These correlations are based on the 3 month boxcar mean anomalies after subtracting the seasonal cycle.(b) Comparison of ph (‰) between the NEW−CO2C13 inversion and the linear NEW−2STEP inversion.We again provide IAV error bars and correlation coefficients between inversion methods.(c) The 3 month boxcar mean anomalies in ph for the North America temperate TransCom region to illustrate the high degree of similarity between both inversion methods (r = 0.86).

Figure 10 .
Figure10.Top panels: the annual averaged Standardized Precipitation and Evaporation Index (SPEI) estimated for the North American temperate domain (map inserts).Middle panel: the annual GPP weighted averaged ph (‰) of vegetation against 13 CO 2 from NEW−CO2C13 (red) and SiBCASA (blue) estimated for the same domain.It illustrates the summertime isoforcing of δ 13 C towards the atmosphere (as wintertime ph has no impact on atmospheric δ 13 C).Lower panel: net carbon uptake (TgC yr −1 ) from NEW−CO2C13 (red) and SiBCASA (blue) estimated for the same domain.The yellow shaded years(2002 and 2011)  indicate significant drought conditions as recorded in SPEI and other independent reports (e.g.,Seager, 2010;Schwalm et al., 2012;Long et al., 2013).These droughts correlate with reductions in annual mean ph , and reductions in the estimated carbon sinks as reported inPeters et al. (2007).

Table 1 .
Summary of assigned δ 13 C model-data mismatch (MDM), the category-averaged and 1σ SD of the innovation χ 2 , and number of sites per category.

Table 3 .
Northern Hemisphere land net carbon uptake (NEE, (Pg C yr −1 )) and land discrimination ( ph , (‰)) 11 year mean estimates, and IAV (±1σ SD) from SiBCASA (prior) and the four inversion experiments.The last row gives the correlation coefficient r between 11 annual mean NEE and ph values.
mate, to estimate NEE but not ph .The SiBCASA terrestrial biosphere model that provides the first-guess NEE and ph of our data assimilation framework based on commonly used drought response parameterizations, simulated neither the large IAV in NEE and ph nor their strong correlation.It is evident from the NEW−2STEP inversion that changes in Figure 11.Weighted SPEI drought index (WDI) versus annual mean isotopic discrimination ph integrated over the North American temperate domain.Results from the SiBCASA biosphere model (blue circles . The lack of www.geosci-model-dev.net/11/283/2018/Geosci.Model Dev., 11, 283-304, 2018 variability in simulated atmospheric δ 13 C found in van der Velde et al. (2013) could well be (partially) ascribed to the lack of sensitivity towards soil moisture stress in SiBCASA.