Trends and drivers of regional sources and sinks of carbon dioxide over the past two decades

Introduction Conclusions References

the global and regional climate and atmospheric CO 2 -driven trends in land and oceanic CO 2 exchanges with the atmosphere over the period 1990-2009, attribute these trends to underlying processes, and quantify the uncertainty and level of model agreement. The models were forced with reconstructed climate fields and observed global atmospheric CO 2 ; Land Use and Land Cover Changes are not included for the DGVMs. Over the period 1990-2009, the DGVMs simulate a mean global land carbon sink of −2.4 ± 0.7 Pg C yr −1 with a small significant trend of −0.06 ± 0.03 Pg C yr −2 (increasing sink). Over the more limited period 1990-2004, the ocean models simulate a mean ocean sink of −2.2 ± 0.2 Pg C yr −1 with a trend in the net C uptake that is indistinguishable from zero (−0.01 ± 0.02 Pg C yr −2 ). The two ocean models that ex- 15 tended the simulations until 2009 suggest a slightly stronger, but still small trend of −0.02 ± 0.01 Pg C yr −2 . Trends from land and ocean models compare favourably to the land greenness trends from remote sensing, atmospheric inversion results, and the residual land sink required to close the global carbon budget. Trends in the land sink are driven by increasing net primary production (NPP) whose statistically significant

Introduction
Soon after the first high-precision measurements of atmospheric CO 2 started in the late 5 1950s, it became clear that the global-mean CO 2 growth rate is substantially lower than expected if all anthropogenic CO 2 emissions remained in the atmosphere (e.g., Keeling et al., 1976). The search for this "missing" carbon and the identification of the processes driving carbon sinks has been one of the dominating questions for carbon cycle research in the past decades (e.g. Tans et al., 1990;Sarmiento and Gruber, 2002 and 10 others). While much progress has been achieved (e.g., Prentice et al., 2001;Sabine et al., 2004;Denman et al., 2007;Le Quéré et al., 2009), and estimates have converged considerably Khatiwala et al., 2012;Wanninkhof et al., 2013), the recent sink rates for the ocean and land, and particularly their changes through time, remain uncertain. A particularly striking observation is that the combined sinks 15 in the ocean and land have increased approximately in line with the increase in fossil fuel emissions in recent decades (Keeling et al., 1995;Canadell et al., 2007;Raupach et al., 2008;Sarmiento et al., 2010;Gloor et al., 2010;Ballantyne et al., 2012). Thus, for the past two decades alone when fossil fuel emissions increased from 6.2 Pg C yr studies using Dynamic Global Vegetation Models (DGVMs) or Ocean Biogeochemical General Circulation Models (OBGCMs) mechanistically represent many of the key land (Prentice et al., 2007) and ocean processes (Le Quéré et al., 2005), and offer the opportunity to investigate how changes in the structure and functioning of land ecosystems and the ocean in response to changing environmental conditions affect 15 biogeochemical cycles. Therefore DGVMs and OBGCMs allow in principal a more comprehensive analysis of surface carbon trends and provide insight into possible mechanisms behind regional trends in the carbon cycle.
There is a growing literature on regional carbon budgets for different parts of the world (Ciais et al., 1995;Phillips et al., 1998;Fan et al., 1998;Pacala et al., 2001; age, and climate effects on productivity, respiration and climate-caused natural disturbances (see Table S1 for details represented in individual models). A particular focus is on the impacts of climate variation and change on land ecosystems at the regional scale as extreme climate events have occurred over the 1990-2009 period across many regions of the World, e.g. including North America (South West USA, 15 2000-2002), Europe (2003), Amazonia (2005), and eastern Australia (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), raising considerable attention in the ecological community on the consequences of recent climate variability on ecosystem structure and function (Allen et al., 2010) and the carbon cycle (Ciais et al., 2005;Van der Molen et al., 2011;Reichstein et al., 2013). Consideration of Land Use Change (LUC) on regional trends is beyond the scope of 20 the present study. There are large uncertainties in the global LULCC flux (Le Quéré et al., 2009), and in comparison the net land use flux represents only approximately 2-3 % of the large land-atmosphere fluxes associated with net primary productivity and heterotrophic respiration. The response of these large fluxes to climate variability and CO 2 are the focus of this study. In addition the net LU flux for the period 1990- 25 2009 will be influenced by earlier LULCC (i.e. legacy fluxes), confounding the analysis. Other companion papers investigate ecosystem response to interannual and seasonal timescales (Piao et al., 2012), and the carbon balance for individual land and ocean air-sea CO 2 flux arising from the increase in atmospheric CO 2 is often referred to as the flux of anthropogenic CO 2 , while the remainder, induced by changes in the natural cycling of carbon in the ocean-atmosphere system is called the "natural" CO 2 component (e.g., Gruber et al., 2009). Although this conceptual separation has its limits (McNeill and Matear, 2013), it provides for a powerful way to understand how different 10 forcings affect the net ocean sink.
DGVM results are compared with estimates of the Residual Land Sink and with remote sensing products indicating trends of greening and browning in the northern region. Regional sources and sink trends are attributed to processes based on factorial simulations.

Dynamic Global Vegetation Models
Following the studies of Le Quéré et al. (2009) and Sitch et al. (2008), a consortium of Dynamic Global Vegetation Model (DGVM) groups set up a project to investigate further the spatial trends in land-atmosphere flux and agreed to perform a factorial set of 20 DGVM simulations over the historical period, 1901-2009. These simulations have contributed to the RECCAP activity (Canadell et al., , 2013. There are now a variety of DGVMs with origins in different research communities that typically contain alternative parameterisations and a diverse inclusion of processes (Prentice et al., 2007;Piao et al., 2013). DGVMs have emerged from the Land Surface Modelling (LSM), for-A total of four different groups have conducted the factorial simulations over the analysis period with three-dimensional OBGCMs and submitted their results to the RECCAP archive. These are MICOM-HAMOCCv1 (BER) (Assmann et al., 2010), CCSM-WHOI using CCSM3.1 (BEC) (Doney et al., 2009a, b), CCSM-ETH using CCSM3.0 (ETH) (Graven et al., 2012), and NEMO-PlankTOM5 (UEA) (Buitenhuis et al., 2010). Details 15 of the models are given in the Appendix of Wanninkhof et al. (2012). Not all model simulations are independent of each other, as several of them share components. BEC and ETH employ the same BOGCM, but differ in their spinup and surface forcing. The employed models have relatively similar horizontal resolution of the order of 1 • to 3 • in longitude and latitude, i.e., none of them is eddy-permitting or eddy-resolving. The four 20 ecosystem/biogeochemical models are also of comparable complexity, i.e., including explicit descriptions for at least one phytoplankton and zooplankton group, with some models considering up to three explicitly modeled groups for phytoplankton and two for zooplankton. All models use the same gas exchange parameterization of Wanninkhof (1992), although with different parameters. In particular, the ETH model used a lower 25 coefficient than originally proposed, yielding a global mean gas transfer velocity that is more than 25 % lower than those of the other models (Graven et al., 2012). This reduction reflects the mounting evidence based on radiocarbon analyses that the orig-20120

Ocean
Unlike how the land models simulations were set up, no common climatic forcing dataset was used for the ocean model simulations. In fact, some models even provided results with different climatic forcings. Models were forced by the NCEP climatic data (Kalnay et al., 1996)

Atmospheric inversion
Simulated trends are compared with those from version 11.2 of the CO 2 inversion product from the Monitoring Atmospheric Composition and Climate -Interim Implementa-5 tion (MACC-II) service (http://www.gmes-atmosphere.eu/). It covers years 1979-2011 and a previous release has been documented by Chevallier et al. (2010). It uses a climatological prior without inter-annual variability, except for fossil fuel.

10
The land models were forced over the 1901-2009 period with changing CO 2 , climate and fixed present-day land use according to the following simulations: S_L1: changing CO 2 only (time-invariant present-day land use mask, fixed preindustrial climate). S_L2: changing CO 2 and climate (time-invariant present-day land use mask). 15 In both cases, DGVMs spin-up were initially performed, with pre-industrial CO 2 , and climate. For DGVMs including the N cycle, N deposition was a time-variant forcing in both simulations, such that the difference between S_L2 and S_L1 includes the synergistic effects of N deposition on CO 2 fertilisation (Zaehle et al., 2010). Figure 2 shows the historical changes in climate, atmospheric CO 2 concentration, 20 nitrogen deposition over the period 1990-2009 used to force the DGVMs. A summary of DGVM characteristics is given in Table 1. A more detailed description of DGVM process representations is given in Table A1.

Ocean
The ocean models have employed two different strategies for creating the initial conditions for the experiments. The first strategy, followed by CCSM-ETH, CCSM-WHOI and BER, involved first a several century long spinup with climatological forcing and with atmospheric CO 2 held constant at its pre-industrial value, bringing these models 5 very close to a climatological steady-state for preindustrial conditions (in some models ∼ 1750; in others ∼ 1850). In the second step, the models are then integrated forward in time through the historical period until 1948, with atmospheric CO 2 prescribed to follow the observed trend and a climatological forcing. The length of the spinup varies between a few hundred years to several thousand years, resulting in differing global integrated drift fluxes, although their magnitudes are substantially smaller than 0.05 Pg C yr −1 with essentially no rate of change. The second strategy, followed by the University of East Anglia, was to initialize the model with reconstructed initial conditions in 1920, and then also run it forward in time until 1948 with prescribed atmospheric CO 2 , repeating the daily forcing conditions of a single year (1980). The modelled ex-15 port production was tuned to obtain an ocean CO 2 sink of 2.2 Pg C yr −1 in the 1990s.
This second method offers the advantage that the model's carbon fields remain closer to the observations compared to the long spinup approach, but it comes at the cost of generating a drift that affects the mean conditions and to a lesser extent the trend. Tests with the model runs of Le Quéré et al. (2010) suggests the drift in the mean CO 2 S_O2: CO 2 and climate, i.e., as S_O1, but models are forced with "realistic" year-toyear variability in atmospheric boundary conditions (ANTH).
The CCSM-based models performed an additional experiment to better separate between the fluxes of natural and anthropogenic CO 2 : S_O3: pre-industrial and climate, i.e., atmospheric CO 2 is fixed at its pre-industrial 5 level, but atmospheric boundary conditions vary as in S_O2 (PIND). From these simulations, only the results from 1990 through 2009 were analysed. Only the UEA and CCSM-WHOI models made available results for the S_O1 and S_O2 simulations for the entire analysis time. The results for the BER model for 2009 are incomplete and the CCSM-ETH simulations extend only to 2007. In order to maintain a sufficiently large set of models, we decided to focus our analysis primarily on the 1990 through 2004 period, but include occasionally also the results through 2009, with the important caveat that the latter are based only on two models.

15
Output variables include those associated with the carbon, hydrological and energy balance over the period 1901-2009. Variables were considered either "Level 1", and essential, or "Level 2", desirable, for additional analysis/studies. Output variables were typically at either the monthly or annual resolution, and at the spatial resolution of the individual model application. In this study we analyse a sub-set of the DGVM output, 20 and primarily focus on the carbon cycle variables, Net Primary Productivity (NPP), Heterotrophic Respiration (RH), and the Net Biome Productivity (NBP), and Leaf Area Index (LAI), a measure of vegetation greenness. The land to atmosphere net CO 2 flux is equal in magnitude to the NBP but has the opposite sign, i.e. we adopt the sign convention where a negative value for the net CO 2 flux represents a carbon sink. Here Introduction , calculated as the annual anthropogenic CO 2 emissions (fossil fuel, cement manufacture and land use C flux) minus the annual CO 2 growth rate and model mean ocean C sink as given by Friedlingstein et al. (2010). The ocean uptake is from the same OGGCMs as the ones used here, and the land use C flux is 10 based on a book-keeping approach from Houghton (2010).
The regional analysis will focus on 3 large regions: northern, tropical and southern land regions (Fig. 1). Within these regions, trends at a finer spatial resolution, from multi-grid cell to the sub-region are analysed. The Northern land region is divided by continent into sub-regions, North America, Europe and North Asia, and into  The comparison of DGVM simulated trends in the northern growing season against satellite-derived NDVI observations was based on 8 models (JULES, LPJ, LPJ-GUESS, NCAR-CLM4, ORCHIDEE, OCN, SDGVM, VEGAS), which provided Leaf Area Index outputs (LAI). The means and trends in the onset, offset and length of growing season were computed. Growing season variables were calculated using the 25 methodology of Murray-Tortarolo et al. (2013). Leaf onset is defined as the day (t) when LAI begins to increase above a critical threshold (CT), defined as: where LAI min and LAI max represent the minimum and maximum LAI over the annual cycle. Similarly, leaf senescence, or offset, is defined as the day when LAI decreases below the CT. The length of the growing season in days is calculated as the offset minus the onset. This calculation was made for each gridcell above 30 • N (i.e. northern extra-tropics) from the models and the satellite data. In addition, any gridcell where LAI 5 varied by less than 0.5 over the annual cycle from the satellite data was considered to be predominantly evergreen (e.g. Boreal Forest), and not considered in the analysis. We also masked out regions where LAI decreases in the summer (drought deciduous vegetation). In addition, when the growing season spans over the end of year (e.g. Mediterranean and some pixels particularly on the southern margin of the domain), we 10 include the first 3 months of the second year in our analysis. Means and trends were calculated using a linear model over the period 1990-2009.

Ocean
The modelling groups provided output on a monthly basis for the years 1990 through 2004 and 2009, respectively, at two levels of priority. Tier one data included the sur- 15 face ocean fields of the air-sea CO 2 flux, oceanic pCO 2 , dissolved inorganic carbon (DIC), alkalinity (Alk), temperature (T ), salinity (S), and mixed layer depth. The second tier data included the biological export at 100 m, the vertically integrated net primary production, and the surface chlorophyll a concentration. Some models supplied also three-dimensional climatological fields of DIC, Alk, temperature, and salinity. 20 To determine the different factors contributing to the modelled trends and variations, we undertook two (linear) separations: The contribution of climate variability on the ocean carbon cycle: X _var = X (S_O2) −X (S_O1), where X is any state variable or flux, where the expression in parentheses represents the results of the corresponding simulation, and where X_var 25 represents the impact of climate change and variability on the ocean carbon cycle.
The contribution of anthropogenic CO 2 : X _ant = X (S_O2) −X (S_O3). 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. For each of the integrations, but particularly for the changing CO 2 and climate simulation S_O2, we analysed the factors contributing to the temporal change in the air-sea CO 2 flux F by a linear Taylor expansion (see e.g., Lovenduski et al. (2007) and Doney et al., 2009a):

BGD
where ws is the wind speed, ice is the sea-ice fraction, sDIC and sAlk are the salinity normalized DIC and Alk concentrations, and where ∂F/∂FS is the change in the air-sea CO 2 flux in response to freshwater fluxes. This latter term includes not only the sensitivity of oceanic pCO 2 to changes in salinity, but also the dilution effects of freshwater on DIC and Alk (see Doney et al., 2009a for details). The partial derivatives were computed directly from the model equations for the mean conditions in each region. The changes in the driving components were derived from the trend computed via a linear regression of the model results and then multiplied by the length of the 15 timeseries.

Land
The ensemble mean global land to atmosphere net carbon dioxide flux from S_L2 is  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a net flux trend over this period ranging between −0.02 and −0.11 Pg C yr −2 , corresponding to DGVMs, OCN and Hyland, respectively (Table 2). DGVMs simulate an increase in the land C sink with an ensemble mean trend of −0.06 ± 0.03 Pg C yr −2 (P < 0.05) over the period 1990 to 2009 (Table 2), in response to changes in climate and atmospheric CO 2 content. The two DGVMs with a fully coupled carbon and nitro-5 gen cycle (CN), OCN and CLM4CN, still simulate an increase in the land sink, but with a trend below the ensemble mean, at −0.02 (P = 0.6) and −0.05 Pg C yr −2 (P = 0.06), respectively. DGVMs suggest the increase in global land sink between 1990 and 2009 is driven by increases in simulated global NPP (Fig. 2). DGVMs simulate an ensemble mean global NPP of 62.9±8.73 Pg C yr −1 over the pe- 10 riod 1990 and 2009 (Table 2). All DGVMs simulate an increase in NPP over this period, with an ensemble mean DGVM trend in NPP of 0.22 ± 0.08 Pg C yr −2 (P = 0.00) (Table 2). Also models with a higher NPP trend also produce a higher NBP trend (Fig. A2). VEGAS, and the two CN models, CLM4CN and OCN, simulate the smallest positive trends in global NPP of 0.11, 0.15 and 0.16 Pg C yr −2 , respectively ( Table 2). The 15 ensemble mean NPP trend of 0.22 ± 0.08 Pg C yr −2 (P = 0.00) from simulation, S_L2 (CO 2 and climate forcing) contrasts with an ensemble trend of 0.19 ± 0.08 Pg C yr −2 (P = 0.00), and 0.03±0.05 Pg C yr −2 (P = 0.24) over the same period for the S_L1 (CO 2 only) and S_L2-S_L1 (the climate effect), respectively (Tables A2 and A3). These results suggest that the simulated increase in global NPP is mainly in response to in-20 creasing atmospheric CO 2 (direct CO 2 fertilization of photosynthesis, in addition to the indirect benefits from an improved water balance in water-limited ecosystems due to the physiological effects of CO 2 on water use efficiency). Again, VEGAS, CLM4CN and OCN, simulate the smallest positive trends in NPP among the DGVMs in response to elevated CO 2 forcing (Table A2). This suggests that the potential CO 2 fertilization effect 25 may be already strongly limited by present day nitrogen availability in some ecosystems (Vitousek and Howarth, 1991). There is more uncertainty among models on the impact of climate changes on global NPP, with only 2 models simulating a significant positive trend (Table A3). DGVMs simulate an ensemble mean global RH of 57.5 ± 9.8 Pg C yr −1 over the period 1990 and 2009 (Table 2). All DGVMs simulate an increase in RH for S2 (CO 2 and climate), with an ensemble mean trend of 0.16 ± 0.05 Pg C yr −2 (P = 0.00) over the period 1990 to 2009 (Table 2). This is lower than the trend in global NPP, thus explaining the trend towards increasing net land carbon uptake. This is unsurpris-5 ing as there is a lagged response in increases in RH relative to NPP, reflecting the turnover time of the newly incorporated plant material. The ensemble mean trend in RH is 0.12 ± 0.06 Pg C yr −2 (P = 0.00) and 0.04 ± 0.02 Pg C yr −2 (P = 0.09) over the same period for the S1 (CO 2 only) and S_L2-S_L1 (the climate effect), respectively (Tables  A2, A3), implying the dominant effect on RH is increased substrate for microbial respiration, with the additional litter input into soils, as a consequence of enhanced NPP, rather than enhanced rates of microbial decomposition with rising temperatures. Nevertheless, the simulated mean residence time (MRT = Soil Carbon/RH) of soil organic matter decreases, in response to warming, which is especially pronounced in high latitude regions (Fig. A3). The differences in land-atmosphere flux trends for the CN 15 models, OCN (−0.02 Pg C yr −2 ) and CLM4CN (−0.05 Pg C yr −2 ) are largely due their differences in RH trends at 0.14 and 0.11 Pg C yr −2 , respectively, rather than differential responses of simulated NPP to elevated CO 2 (Table 2). Only 4 DGVMs provided wildfire fluxes (CLM4CN, LPJ, LPJ-GUESS, SDGVM). No significant trends in the global wildfire flux were reported by any of the DGVMs. Only 20 two DGVMs (CLM4CN and LPJ) simulated small significant trends of opposite sign in the wildfire flux to the 90 % confidence interval for S_L1 (CO 2 only), possibly reflecting the counteracting effects of CO 2 on plant water use, leading to more moist land surface (reduced flammability), and elevated productivity (increased fuel loads).

25
The global ocean is simulated to have acted as a very substantial sink for atmospheric CO 2 , but one that has increased only slightly over the last two decades (see also BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.  Wanninkhof et al., 2013). The mean ocean sink in the 4 models (CCSM-ETH, CCSM-WHOI, UEA and BER), increased from ∼ −2.0 Pg C yr −1 in the early 1990s to ∼ −2.1 Pg C yr −1 during the first five years of the 21st century (Fig. 4). The overall sink is largely a consequence of the increase in atmospheric CO 2 (i.e., it corresponds mostly to the uptake flux of anthropogenic CO 2 ), but it includes a substantial perturbation flux 5 stemming from the impact of climate variability and change on the ocean carbon cycle. This is confirmed if we separate the mean and variable components by using our factorial experiments, i.e., by using S_O1 results to identify the ocean uptake in the absence of climate variability and change, and the difference between S_O2 and S_O1 as measure of the impact of climate change. This separation reveals that in the ab- is entirely driven by the increase of atmospheric CO 2 and is -integrated globally -numerically equivalent to the ocean uptake flux of anthropogenic CO 2 . Climate variability and change modified these fluxes, and particularly the trend in these models. The four models suggest an enhancement of the uptake in the early 1990s (1990)(1991)(1992)(1993)(1994) of about 20 −0.2 Pg C yr −1 , turning into a reduction of the uptake in the subsequent period (1995)(1996)(1997)(1998)(1999), followed by a further reduction in the 2000-2004 period of ∼ +0.1 Pg C yr −1 . This trend toward reduced uptake in response to climate variability and change of +0.03 Pg C yr −2 nearly completely compensates for the anthropogenic CO 2 driven increase in uptake, causing the overall uptake of CO 2 to have a nearly flat trend over the 25 1990 through 2004 period of < 0.01 Pg C yr −2 The same tendencies are found for the two models that extend over the entire 1990 through 2009 period: climate change and variability reduces in these models the CO 2 driven trend of −0.04 Pg C yr −2 by more than +0.02 Pg C yr −2 , to around −0.02 Pg C yr −2 . 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. Consideration of the different factors affecting the ocean carbon sink following our Taylor expansion, we find increasing sea-surface temperature to be globally one of the most important drivers for the positive trends (reduced sinks) induced by climate change and variability. Over the 1990 through 2004 period, the surface ocean warmed, on average, by 0.004 • C yr −1 (0.005 • C yr −1 from 1990 through 2009). Isochemically, this leads to an increase of the oceanic pCO 2 of ∼ 0.06 µatm yr −1 , which appears small. However, it needs to be compared with the trend in the global mean air-sea pCO 2 difference of about ∼ 0.1 µatm yr −1 that is required in order to generate a trend in the ocean uptake of −0.03 Pg C yr −2 (see e.g., Matsumoto and Gruber, 2005;Sarmiento and Gruber, 2006). 3.2 Regional trends

Land
Spatial correlations between trends in drivers and carbon fluxes (NBP, NPP, RH) were significant for precipitation (0.36, 0.5, 0.48 respectively; P < 0.05), but not for temperature (see Table A11), underlining the sensitivity of carbon cycle models to precipitation 15 changes at the multi-year to decadal timescale.

Northern land
All DGVMs agree on a land C sink over the northern land region, with a mean landatmosphere flux of −1.03±0.30 Pg C yr −1 , over the period 1990-2009 (Fig. A4, Table 2).
The ensemble mean land-atmosphere flux trend is near zero for this region between 20 1990 and 2009 (Fig. A5). Although 4 of 9 DGVMs simulate reductions in the land sink over the northern land region, none of the trends from the 9 DGVMs are significant at the 95 % confidence level (Table 2). Of particular interest are sub-regions with a simulated positive land-atmosphere flux trend (Fig. 5), implying a diminishing sink of atmospheric CO 2 , or an increasing source of CO 2 to the atmosphere. At least 6 mod-BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | els out of 9 agree on a decreasing regional land sink across some areas in Temperate North America, Eastern Europe, North-East China and Mongolia (Fig. 5). These net flux trends largely correspond to regions with negative trends in precipitation (Fig. 6). Over the northern region, DGVMs simulate an ensemble mean NPP of 24.1 ± 4.48 Pg C yr −1 , which represents almost 40 % of the global total (Table 2). All DGVMs 5 simulate an increase in northern NPP over this period, with a trend in NPP of 0.06 ± 0.02 Pg C yr −2 (P = 0.00) ( Table 2). However, enhanced productivity in the northern land region accounts for only around 29 % of the simulated global trend in NPP. The ensemble mean NPP trend of 0.06 ± 0.02 Pg C yr −2 (P = 0.00) from simulation S2 (CO 2 and climate forcing) compares to a trend of 0.07 ± 0.03 Pg C yr −2 (P = 0.00), and 10 −0.00 ± 0.04 Pg C yr −2 (P = 0.85) for the S1 (CO 2 only) and S2-S1 (the climate effect), respectively (Tables A2, A3 The ensemble mean NPP for North America, Europe and North Asia are 7.78 ± 1.38 Pg C yr −1 , 5.08±1.40 Pg C yr −1 , and 11.20±1.99 Pg C yr −1 , respectively (Table A5). Despite the two-fold difference in mean NPP, NPP trends for North America, Europe and North Asia are similar, at 0.021±0.008 Pg C yr −2 (P = 0.02), 0.018±0.006 Pg C yr −2 20 (P = 0.00) and 0.024±0.015 Pg C yr −2 (P = 0.01), respectively (Table A5). DGVMs simulate larger mean NPP in Temperate compared to Boreal regions with mean NPP in Temperate North America and Asia of 4.22 ± 0.83 Pg C yr −1 and 7.09 ± 1.03 Pg C yr −1 , respectively, compared to 3.57±1.21 Pg C yr −1 and 4.12±1.46 Pg C yr −1 in Boreal North America and Asia, respectively (Table A5). DGVMs simulate a significant positive trend 25 in Boreal North America and Boreal Asia of 0.014 ± 0.007 Pg C yr −2 (P = 0.02) and 0.018 ± 0.006 Pg C yr −2 (P = 0.01), respectively. DGVM NPP trends in both Temperate North America and Asia are smaller than those for Boreal regions but are not significant at the 95 % confidence level (Table A5). 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. In response to warming, models simulate an earlier onset (ensemble mean model trend = −0.078 ± 0.131 day yr −1 ) and delayed termination of the growing season (0.217 ± 0.097 day yr −1 ) based on LAI, and thus a trend towards a longer growing season in the northern extra-tropics (0.295±0.228 day yr −1 ) (Fig. 7). This is in broad agreement with observed greening trends: onset = −0.11 day yr −1 ; offset = 0.252 day yr −1 5 and growing season length = 0.361 day yr −1 ). There is less agreement among models on reproducing the observed browning trends in some regions of the boreal forest. DGVMs simulate an ensemble mean RH of 21.8 ± 4.6 Pg C yr −1 across the northern land region ( America and Asia, respectively (Table A6). DGVMs simulate a significant positive trend in Boreal North America and Boreal Asia of 0.012 ± 0.003 Pg C yr −2 (P = 0.02) and 15 0.015 ± 0.004 Pg C yr −2 (P = 0.01), respectively (Table A6). DGVM RH trends in both Temperate North America and Asia are smaller at 0.008 ± 0.005 Pg C yr −2 (P = 0.01) and 0.010 ± 0.009 Pg C yr −2 (P = 0.08), respectively than those for boreal regions and both significant at the 90 % confidence level. This is because of relatively smaller increases in substrate (i.e. NPP), in temperate regions, and greater warming in boreal 20 regions stimulating microbial decomposition, reducing mean residence time of carbon in soils (MRT = Soil Carbon/RH; see Fig. A3). No significant trends in the wildfire flux were reported by any of the DGVMs for the northern land region. However DGVMs agree on simulating a small negative trend in wildfire flux across Boreal North America and Tundra, with ensemble mean trends of 25 −0.002 ± 0.003 Pg C yr −2 (P = 0.01) and −0.001 ± 0.002 Pg C yr −2 (P = 0.08), respectively at the 90 % confidence level. 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.

Tropical Land
All DGVMs simulate a land C sink over the last two decades, in response to climate variability and changes in atmospheric CO 2 concentration, with an ensemble mean land-atmosphere flux of −0.96 ± 0.43 Pg C yr −1 over the Tropical land region (Table 2, Fig. A4). All DGVMs simulate an increasing land sink, with an ensemble mean trend of 5 −0.04 ± 0.01 Pg C yr −2 (P = 0.04) for this region (Table 2, Fig. A5). The ensemble mean NBP trend over the tropical land region is significant (P < 0.05) and represents 65 % of the increase in land sink over the last two decades. DGVMs simulate a significant land-atmosphere flux trend (i.e. increasing sink) of −0.016±0.007 Pg C yr −2 (P = 0.00) and −0.008 ± 0.006 Pg C yr −2 (P = 0.05) across Tropical Asia and Equatorial Africa, 10 respectively (Table A4). Land-atmosphere flux trends for the Tropical South American forest and North Africa Savanna regions are not significant (Table A4). DGVMs simulate an ensemble mean NPP of 27.0 ± 4.29 Pg C yr −1 averaged over the tropical region representing 42 % of the global total (Table 2). All DGVMs simulate a significant increase in tropical NPP over this period, with an ensemble mean 15 trend in NPP of 0.10 ± 0.03 Pg C yr −2 (P = 0.00) for S2 (Table 2). This compares to a trend of 0.09 ± 0.03 Pg C yr −2 (P = 0.00), and 0.02 ± 0.02 Pg C yr −2 (P = 0.33) over the same period for the S1 (CO 2 only) and S2-S1 (the climate effect), respectively (Tables A2 and A3). Again the simulated trend in NPP is dominated by the simulated response of ecosystems to elevated atmospheric CO 2 content. DGVMs simu-20 late positive NPP trends of 0.038 ± 0.015 Pg C yr −2 (P = 0.00), 0.031 ± 0.011 Pg C yr −2 (P = 0.00), 0.024 ± 0.01 Pg C yr −2 (P = 0.00), 0.01 ± 0.008 Pg C yr −2 (P = 0.07), across Tropical South America forests, Tropical Asia, Equatorial Africa, and North Africa savanna, respectively (Table A5). Nevertheless there are some areas in Tropical South America and in the North Africa Savanna regions with negative trends in NPP (Fig. 6). 25 DGVMs simulate an ensemble mean RH of 24.49 ± 4.75 Pg C yr −1 over the Tropical land region (Table 2). All DGVMs simulate an increase in RH over the period 1990-2009, with a significant trend in the ensemble mean RH of 0.065 ± 0.025 Pg C yr −2 (P = 20134 BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 0.00). This can be attributed largely to the response of ecosystems to elevated CO 2 ; the ensemble mean trend in RH from S1 is 0.059 ± 0.022 Pg C yr −2 (P = 0.00) (Table  A2). No significant trends in the wildfire flux were reported by any of the DGVMs for the tropical land region. However DGVMs agree on simulating a negative trend in wildfire 5 flux across Equatorial Africa and Tropical Asia, with an ensemble means trends of −0.004 ± 0.002 Pg C yr −2 (P = 0.00) and −0.003 ± 0.003 Pg C yr −2 (P = 0.03).

Southern land
All DGVMs agree on a net land sink over the southern land during the last two decades, with an ensemble mean land-atmosphere flux of −0.38 ± 0.29 Pg C yr −1 ( Table 2, 10 Fig. A4). The ensemble mean land-atmosphere flux trend is −0.02 ± 0.02 Pg C yr −2 (P = 0.20) for this region (Fig. A5). Although all DGVMs simulate an increase in the land sink over the southern extra-tropics, only trends for HYL and ORC are significant at the 95 % confidence level (Table 2). Ensemble mean land-atmosphere flux trends are significant for Temperate South America, Australia and New Zealand, and Southern Africa 15 regions at 0.005 ± 0.005 Pg C yr −2 (P = 0.05), −0.003 ± 0.004 Pg C yr −2 (P = 0.68), and −0.022 ± 0.011 Pg C yr −2 (P = 0.01), respectively (Table A4). For Southern Africa, all DGVMs simulate an increase in the land sink in response to climate changes (S_L2-S_L1) over this period (5 out of 9 are significant to the 90 % confidence level), with an ensemble mean land-atmosphere flux of −0.018 ± 0.01 Pg C yr −2 (P = 0.02) (Ta-20 ble A7). This compares with a land-atmosphere flux trend of −0.004 ± 0.002 Pg C yr −2 (P = 0.08) in response to CO 2 forcing only (Table A7). This is likely related to an increase in precipitation over this region (Fig. 6).
In contrast, the simulated decrease in land sink for Temperate South America is associated with a decrease in precipitation over 1990-2009. All DGVMs, except HYL, 25 simulate a decrease in land uptake in this region for the S_L2-S_L1 simulation (however only 4 out of 9 models are significant at the 90 % confidence level), and an ensemble mean land-atmosphere flux trend is 0.005 ± 0.005 Pg C yr −2 (P = 0.06) (Table A8). 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.  (Table 2). All DGVMs simulate an increase in NPP over this period, with a significant ensemble mean trend 5 of 0.05 ± 0.03 Pg C yr −2 (P = 0.01), i.e. the southern land region accounts for around 25 % of the simulated global trend in NPP. This compares to ensemble mean trends of 0.03 ± 0.02 Pg C yr −2 (P = 0.00), and 0.02 ± 0.01 Pg C yr −2 (P = 0.34) over the same period for the S1 (CO 2 only) and S2-S1 (the climate effect), respectively (Tables  A2, A3). Southern Africa is the only southern sub-region with a significant trend of 10 0.041 ± 0.018 Pg C yr −2 (P = 0.00) (Table A5). This large positive trend in NPP is due to a positive response of plant production to both CO 2 (NPP S1 = 0.014±0.007 Pg C yr −2 , P = 0.00), and climate (NPP S2-S1 = 0.027±0.012 Pg C yr −2 , P = 0.01) (Table A7). An increase in NPP over Southern Africa is likely in response to increases in precipitation over the last two decades (Fig. 5), whereas a decrease in NPP over Temperate South

15
America is likely due to a decrease in precipitation. DGVMs simulate an ensemble mean RH of 11.20 ± 3.60 Pg C yr −1 over the Southern land region (Table 2). All DGVMs simulate an increase in RH over the period 1990-2009, with a significant trend in the ensemble mean RH of 0.03 ± 0.02 Pg C yr −2 (P = 0.00). This is only partly explained by the response of ecosystems to elevated CO 2 ; 20 over Southern Africa the ensemble mean trend in RH from S1 is 0.02 ± 0.01 Pg C yr −2 (P = 0.00), and a climate induced positive trend in RH, 0.01 ± 0.00 Pg C yr −1 (P = 0.00) (Tables A2 and A7). No significant trends in the wildfire flux were reported by any of the DGVMs for the southern land region. However DGVMs agree on simulating a negative trend in wildfire 25 flux across Southern Africa, with an ensemble mean trend of −0.006 ± 0.002 Pg C yr −2 (P = 0.03). 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.

Qualitative change in processes
A qualitative assessment of the differential responses of the underlying land processes to changes in environmental conditions, and their contribution to the sink-source landatmosphere flux trends are shown in Fig. 8. Many regions are simulated to have a negative land-atmosphere flux trend, with increases in NPP leading increases in RH. How-5 ever the locations with positive trends over the period 1990 to 2009, i.e. the red colours in Fig. 8, generally fall into two categories: First, regions where models simulate a positive trend in NPP, but an even larger positive trend in RH (eastern Europe, South East USA, Amazonia, South China, North America Tundra). Warming is likely to enhance both NPP and RH in high latitude ecosystems, but primarily RH in low latitudes. Re-10 duced precipitation may partially or fully offset the benefits of elevated atmospheric CO 2 abundance on NPP, and the response of RH to changes in precipitation is not obvious, as this depends on the initial soil moisture status. This is because microbial activity increases with increasing soil moisture at low moisture levels, before reaching a maximum activity, and then begins to decline as water fills the soil pore spaces and 15 oxygen becomes more limiting to respiration. Locations in eastern USA, southern Asia, northern boreal China, south-eastern South America, western and southern Australia are simulated to have negative NPP trends over the last two decades, as a result of reduced rainfall, and there is a less negative trend in RH, possibly due to a reduction in microbial respiration rates with increased soil dryness. The warming and drying in 20 Inner Asia (North East China and Mongolia) and southern Australia is simulated to reduce the rate of microbial decomposition in these regions (Fig. A3), which partly opposes the NPP-driven lagged decrease of RH. The source trend in eastern Europe is simulated as a combination of a negative trend in NPP, as a result of a combination of elevated temperatures and reduced precipitation (i.e. soil drying), and a positive trend 25 in RH, despite reduced plant litter input.

20137
BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.

Regional fluxes
The large-scale distribution of the modeled mean surface fluxes consists of strong outgassing in the tropical regions, especially in the Pacific, and broad regions of uptake in the mid latitudes, with a few regions in the high latitudes of particularly high uptake, 5 such as the North Atlantic (Fig. 9). This pattern is largely the result of the exchange flux of natural CO 2 that balances globally to a near zero flux, but exhibits regionally strong variations (Gruber et al., 2009). Superimposed on this natural CO 2 flux pattern is the uptake of anthropogenic CO 2 , which leads to uptake everywhere, but with substantial regional differences. Large anthropogenic CO 2 uptake fluxes occur in the regions of 10 surface ocean divergence, such as the equatorial Pacific and particularly the Southern Ocean (Sarmiento et al., 1992;Mikaloff Fletcher et al., 2006). This is a result of the divergence causing waters to come from below to the surface which have not been exposed to the atmosphere for while, thereby permitting them to take up a substantial amount of anthropogenic CO 2 . This reduces the outgassing that typically characterizes 15 these regions as a result of these upwelling waters bringing with them also high carbon loads from the remineralization of organic matter. Over the analysis period, the air-sea CO 2 fluxes exhibit only a remarkably small trend in most places with some regions increasing in uptake, while others show a positive flux anomaly, i.e., lesser uptake. Thus the small global trend in ocean uptake over 20 the 1990 through 2004 analysis period is a result of also the individual regions having relatively modest trends.

Process-decomposition
The regional flux trends are, however, much smaller than expected from an ocean with constant circulation that is only responding to increasing atmospheric CO 2 and hence 25 would tend to increase its uptake of anthropogenic CO 2 through time (Fig. 10) absence of climate variability and change, all regions would have flux density trends of more than −0.05 g C m −2 yr −2 , with some regions, such as the Southern Ocean exceeding −0.15 g C m −2 yr −2 . But climate variability and change compensate these negative trends in every single region by increasing them by +0.04 g C m −2 yr −2 or more (with the exception of the South Pacific), such that the overall trends fluctuate from region to 5 region around zero (Fig. 10). The largest reductions in trends are simulated to occur in the North and Equatorial Pacific and in the North Atlantic, where they even cause a change in the sign of the overall trend. A similar, although slightly more moderate pattern is seen if the analysis is undertaken for the entire 1990 through 2009 period on the basis of two models only. The most important difference is found in the North 10 Atlantic, where the climate variability impact is substantially smaller, and not offsetting the anthropogenic CO 2 trend when analyzed for 1990-2009. The mechanisms driving the oceanic flux trends differ between the analyzed regions. In some regions, surface ocean warming dominates and hence reduces or even cancels increasing ocean anthropogenic CO 2 uptake, as is the case globally (Roy et al.,  1960-1988and 1989). The ensemble mean land-atmosphere flux increased by −1.11 Pg C yr −1 for the same period, compared to the estimated RLS increase of −0.88 Pg C yr −1 from BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.  [1990][1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005][2006][2007][2008][2009]. Although encouraging, these results should be interpreted with caution because the inversion accounts for any trend in the land use change flux over this period, 5 whereas DGVMs had fixed land use. For this reason we do not compare results at scales finer than the zonal one. There is empirical evidence for a large increase in biomass in intact forest in the tropical South America and Africa (Pan et al., 2011;Baker et al., 2004;Lewis et al., 2009a, b), which is consistent with the DGVM projections presented here. Subsequently, Lewis et al. (2009b) found broad agreement between biomass trends from observations and from a suite of carbon cycle models applied with 20th century forcing of climate and atmospheric CO 2 content, using a similar protocol to the current analysis. This suggests a large possible role of CO 2 fertilization in stimulating plant productivity in tropical ecosystems. DGVMs suggest a large component of the uptake trend is associated 15 with a positive NPP response to elevated CO 2 , which is broadly consistent with the enhancement of forest production due to CO 2 observed in Free Air Carbon Enrichment (FACE) experiments (Norby et al., 2005), albeit they are largely located in temperate forest ecosystems. However, recent studies have highlighted the role of nitrogen in limiting the long term CO 2 response (Canadell et al., 2007;Norby et al., 2010) in these 20 ecosystems. The long-term plant response to elevated CO 2 is likely affected by nutrients and its impact on plant C allocation (Zaehle et al., 2014), however only two out of the nine models used here (CLM4CN and OCN) include interactive nutrient cycling (see DGVM characteristics, Table A1).
In contrast to the large trend in net C uptake across the tropics, DGVMs simulate no 25 statistically significant trend over the northern land. In particular, trends in NPP over Temperate regions are smaller than those in Boreal regions, and are also not significant. Many temperate areas have experienced a decrease in rainfall between 1990 and 2009, and suffered periods of prolonged and severe drought. Examples include BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.  (McDowell et al., 2008Anderreg et al., 2012), and the 2003 summer heatwave in Europe (Ciais et al., 2005). These are likely to have had a detrimental effect on the land C uptake (Ciais et al., 2005). This is particularly relevant as climate models agree on a future warming and reduced summer precipitation over continental mid-latitudes (IPCC, 2007). 5 Modelling ecosystem structure and function in water stressed environments remains a challenge for DGVMs (Morales et al., 2005;Keenan et al., 2009). In general, there is a need for a greater understanding of the mechanisms behind drought-induced plant mortality (Allen et al., 2010;McDowell et al., 2011), and changes in plant water use (De Kauwe et al., 2013). Incorporating hydraulic failure under drought stress improved the ability of the LPJ-GUESS model to simulate the distribution and productivity of xeric and arid vegetation types on a global scale (Hickler et al., 2006); however, the latter is not included in most models. In general, an improvement in the description of vegetation is needed with the adoption of more plant trait information now available (Kattge et al., 2011). 15 Satellite observations suggest a general greening trend in high latitudes with an earlier onset and longer growing season in high latitude ecosystems, which is reproduced by the DGVMs. Observations suggest a greening Tundra and a slower greening and possible browning in some regions of the boreal forest (Tucker et al., 2001;Bhatt et al., 2010), especially in North America (Beck and Goetz, 2011). In tundra ecosystems, an 20 earlier onset is attributed to warming and earlier snowmelt. In these ecosystems radiation melts snow in the early spring, and the start of growing season corresponds to near peak in radiation. Thus any temperature induced early snowmelt (McDonald et al., 2004;Sitch et al., 2007) is likely to enhance plant production. Warming may not have such a great effect on the offset of the growing season in Arctic tundra ecosystems as 25 this may be driven primarily by radiation. DGVMs simulate a significant positive trend in NPP Boreal North America and Boreal Asia and Tundra. Nitrogen limitation is also likely to constrain the carbon cycle at high latitudes. Only 2 out of 9 DGVMs studied BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. here included a fully interactive carbon and nitrogen cycle, and it was not possible to quantify N-limitations effects on regional trends in this study. DGVMs simulate decreasing NPP across North East China and Mongolia, contributing to the overall decreasing land uptake trend, in response to recent climate. In a regional study, Poulter et al. (2013) investigated the differential response of cool semi-5 arid ecosystems to recent warming and drying trends across Mongolia and Northern China, using multiple sources of evidence, including, LPJ DGVM, FPAR remotely sensed data (derived from GIMMS NDVI3g) and tree-ring widths. They found coherent pattern of high-precipitation sensitivity across data sources, which showed some areas with warming-induced springtime greening and drought-induced summertime browning, and limitations to NPP explained mainly by soil moisture.
Browning is a consequence of regional drought, wildfire and insect outbreak, and their interaction, especially in North America (Beck and Goetz, 2011). Disturbance plays a key role in the ecology of many global ecosystems. For example, wildfire plays a dominant role in the carbon balance of boreal forest in central Canada and other re-  (Kurz et al., 2008). In general disturbance and forest management are inadequately represented by the current generation of DGVMs, even though several models include simple prognostic wildfire schemes (Table A1), 20 while some are starting to include other disturbance types such as insect attacks (Jönsson et al., 2012) and windthrow (Lagergren et al., 2012). The improvement of DGVMs to include representations of globally and regionally important disturbance types and their response to changing environmental conditions is a priority.

25
The investigated BOGCMs consistently simulate an ocean characterized by a substantial uptake of CO 2 from the atmosphere, but with a global integrated trend in the two recent decades that is substantially smaller than that expected based on the in-20142 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | crease in atmospheric CO 2 . Results based on the predictions from ocean inversion and ocean Green function methods (Mikaloff Fletcher et al., 2006;Gruber et al., 2009;Khatiwala et al., 2010) suggest an increase in ocean uptake with a trend of the order of −0.04 Pg C yr −2 over the analysis period (see also Wanninkhof et al., 2013). These latter methods assume constant circulation, while our simulations here include the impact 5 of climate variability and change.
Our analyses reveal indeed that recent climate variability and change has caused the ocean carbon cycle to take up less CO 2 from the atmosphere than expected on the basis of the increase in atmospheric CO 2 , i.e., it reduces the efficiency of the ocean carbon sink. Globally, we demonstrated that this efficiency reduction is primarily a re-10 sult of ocean warming, while regionally, many more processes (e.g., wind changes, alkalinity/DIC concentrations changes) are at play.
Is this reduction in uptake efficiency over the analysis period the first sign of a positive feedback between global warming and the ocean carbon cycle -or alternatively, could it be due to natural decadal-scale variability in air-sea CO 2 fluxes? Without a formal 15 attribution study, it is not possible to provide a firm answer. We suspect that the majority of the trend in the efficiency is due to "natural" decadal-scale variability, however, largely based on the results of McKinley et al. (2011) and Fay and McKinley (2013) who showed that whereas trends in oceanic pCO 2 (and air-sea CO 2 fluxes) are variable on a decadal time-scale, they converge towards atmospheric pCO 2 trends when 20 analyzed over a longer 30 yr period for most global regions. Nevertheless, they also show that in the permanently-stratified subtropical gyre of the North Atlantic, warming (partly driven by anthropogenic climate change) has started to reduce ocean uptake in recent years. In the Southern Ocean, where Le Quéré et al. (2007) and Lovenduski et al. (2008) used models to suggest a reduction in ocean carbon uptake efficiency 25 over the past 25 yr in response to increasing Southern Ocean winds, Fay and McKinley (2013) concluded that the data are insufficient to draw any conclusions.
We should note that the associated uncertainties remain large. Of particular concern is the moderate success of the models to simulate the time-mean ocean sinks and BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | their long-term seasonal cycle (e.g. McKinley et al., 2006). Furthermore, some of the models underestimate the oceanic uptake of transient tracers such as anthropogenic radiocarbon (see e.g., Graven et al., 2012). Such a reduction in the oceanic uptake efficiency is also not suggested by independent measures of oceanic CO 2 uptake, such as the atmospheric O 2 /N 2 method (Manning and Keeling, 2006;Ishidoya et al., 2012), 5 albeit the large uncertainties in these estimates make the determination of trends in uptake highly uncertain.

Reducing uncertainty in regional sinks
In order to better quantify the regional carbon cycle and its trends, DGVM and ocean carbon cycle models need to improve both process representations and model evalua-10 tion and benchmarking (Luo et al., 2012). There is a need for up to date global climate and land use and cover change datasets to force the DGVMs, as well as a deeper investigation of the quality and differences between the different reanalysis products used to force ocean carbon cycle models. Also techniques such as detection and attribution can be applied to elucidate trends in the regional carbon cycle and their drivers.

Model benchmarking
There is a critical need for comprehensive model benchmarking, as a first step to attempt to reduce model uncertainty. Several prototype carbon cycle benchmarking schemes have been developed (Randerson et al., 2009;Cadule et al., 2010). A more in depth evaluation and community benchmarking set needs to be agreed and im-20 plemented, which also evaluates models for their implicit land response timescales (especially relevant in the discussion on future tipping elements and non-linear future responses) and for the simulated carbon, water and nutrient cycles. New emerging frameworks now exist (Blyth et al., 2011;Abramowitz, 2012;Luo et al., 2012;Dalmonech and Zaehle, 2013). 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al.

Model resolution
Simulated ocean carbon dynamics may be sensitive to horizontal resolution, particularly as model resolution improves sufficiently to adequately capture mesoscale eddies. Mesoscale turbulence influences the ocean carbon cycle in a variety of ways, and the present eddy parameterizations may not adequately capture the full range of effects 5 and the responses to climate variability and change. For example, mesoscale processes are thought to modulate biological productivity by altering the supply of limiting nutrients (Falkowski et al., 1991;McGillicuddy et al., 1998;Gruber et al., 2011). A particularly crucial issue involves the wind-driven overturning circulation in the Southern Ocean, where non-eddy resolving models indicate a strong sensitivity of the overturn- 10 ing circulation and ocean carbon uptake to surface wind stress (Le Quere et al., 2007;Lovenduski et al., 2008). Some eddy-resolving models in contrast suggest that enhanced wind stress is dissipated by increased eddy activity, leading to only a small increase in overturning (Boning et al., 2008), though more recent results indicate a larger response (Gent and Danabasoglu, 2011).

Model structure
There is a need for improved representation of ecological processes in land and ocean models, e.g. nutrient cycling (N, P), disturbance (wildfire, wind-throw, insects), land use and land cover change in land models and better representation of the key functional diversity in ocean biogeochemical models. DGVMs need to represent land use and 20 land cover changes, forest management and forest age, to improve estimates of the regional and global land carbon budget. There are recent developments to include nutrient dynamics, mostly nitrogen, into global land biosphere models (as reviewed by Zaehle and Dalmonech, 2011). Too few model simulations are available to date to allow for an ensemble model trend assessment. However, a few general trends appear 25 robust: As evident from Table 2, C-N models generally show less of a response to increasing atmospheric CO 2 due to nitrogen limitation of plant production. N dynam-BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. ics further alter the climate-carbon relationship, which tend to reduce the C loss from temperate and boreal terrestrial ecosystems due to warming; however, with a considerable degree of uncertainty (Thornton et al., 2009;Sokolov et al., 2008;Zaehle et al., 2010). Changes in the nitrogen cycle due to anthropogenic reactive nitrogen additions (both fertiliser to croplands and N deposition on forests and natural grasslands) fur-5 ther modify the terrestrial net C balance and contribute with −0.2 to −0.5 Pg C yr −1 to current land sink, with a geographic distribution closely related to that of anthropogenic reactive nitrogen deposition (Zaehle and Dalmonech, 2011). Zaehle et al. (2011), using the OCN model, have estimated the 1995-2005 trend in land uptake due to N deposition as −1.1 ± 1.7 Tg C yr −2 , with strong regional differences, depending on the regional trends in air pollution and reactive N loading of the atmosphere and the nitrogen status of the ecosystems, which are generally lower in less responsive in ecosystems close to nitrogen saturation highly polluted regions.
There are several additional land processes that have not been considered in this current multi-model analysis. These include the effects of aerosols and tropospheric 15 ozone on the carbon cycle. Unlike a global forcing agent such as CO 2 , the effects of air pollutants (aerosols, NO x , and O 3 ) with their shorter atmospheric lifetimes are at the regional scale. Aerosol induced changes in radiation quantity and quality (i.e. the ratio of diffuse to direct) affect plant productivity and the land sink (Mercado et al., 2009). From around 1960 onwards until the 1980s radiation levels declined across 20 industrialized regions, a phenomenon called "global dimming", followed by a recent brightening in Europe and North America with the adoption of air pollution legislation. Reductions in acid rain have been found to greatly influence trends in riverine DOC, vegetation health, and rates of soil organic matter decomposition. Tropospheric ozone is known to be toxic to plants and lead to reductions in plant productivity, and reduce 25 the efficiency of the land carbon sink Anav et al., 2011). Drivers of the land carbon sink related to air pollution, e.g. N deposition, acid precipitation, diffuse and direct radiation, and surface O 3 have varied markedly in space and time over recent decades. Similar gaps need to be addressed in ocean biogeochemical models. The ecosystem modules in the current generation of BOGCMs lack the ability to assess many of the suggested mechanisms by which climate and ocean acidification could alter marine biogeochemistry and ocean carbon storage. Proposed biological processes that could influence ocean carbon uptake and release involve, for example, decoupling of carbon 5 and macro-nutrient cycling, changes in micro-nutrient limitation, variations in elemental stoichiometry in organic matter, and changes in the vertical depth-scale for the respiration of sinking organic carbon particles (e.g., Boyd and Doney, 2003;Sarmiento and Gruber, 2006). Some advances have been made with the incorporation of dynamic iron cycling and iron limitation, multiple plankton groups, calcification, and nitrogen fixation 10 (Le Quéré et al., 2005). However, the evaluation of these aspects of the models is currently hindered by both data and process-level information limitations.

Climate and land use and cover datasets
In addition to model structure, the choice of climate forcing and model initial conditions can also contribute to differences between the simulated terrestrial carbon sink. 15 At regional scales, differences in land cover can introduce ∼ 10 % uncertainty in simulated regional-scale GPP (Jung et al., 2007;Quaife et al., 2008) and about a 3.5 % uncertainty for global NPP. Climate forcing uncertainty tends to have larger effects on carbon flux uncertainty than land cover (Hicke, 2005;Poulter et al., 2011) with up to 25 % differences in GPP reported over Europe (Jung et al., 2007 and a 10 % difference 20 for global NPP (Poulter et al., 2011). Climate forcing uncertainty and land cover (i.e., PFT distributions) can alter long-term trends in NBP and inter-annual variability of carbon fluxes to climate (Poulter et al., 2011). The DGVMs applied here did not consider LULCC. This is an active area of research; models need a consistent implementation of LULCC. Uncertainties in the simulated net land use flux are associated with as- 25 sumptions on the implementation of LULCC gridded maps (e.g. whether conversion to cropland in a grid-cell is taken preferentially from grassland, forest, or from both), simulated biomass estimates, and subsequent decomposition rates. However DGVMs 20147 BGD 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. offer the exciting prospect to disentangle the component fluxes associated with land use (e.g. direct emissions, and legacy fluxes), and to separate the environmental and direct human impacts on the net LU flux.

Conclusions
Land models suggest an increase in the global land net C uptake over the period 5 1990-2009, mainly driven by trends in NPP, in response to changes in climate and atmospheric CO 2 concentration. Over the same period, ocean models suggest a negligible increase in net ocean C uptake; a result of ocean warming counteracting the expected increase in ocean uptake driven by the increase in atmospheric CO 2 . At the sub-regional level, trends vary both in sign and magnitude, particularly over land. Areas 10 in Temperate North America, eastern Europe, North-East China, show a decreasing regional land sink trend, due to regional drying, suggesting a possibility for a transition to a net carbon source in the future if drying continues or drought become more severe and/or frequent. In the ocean, the trends tend to be more homogeneous, but the underlying dynamics differ greatly, ranging from ocean warming, to winds, and to changes in 15 circulation/mixing and ocean productivity, making simple extrapolations into the future difficult.
Our conclusions need to be viewed with several important caveats: Few land models include a prognostic representation of wildfire and no land model represents other disturbances and their interactions, indicating the model response to warming and drought 20 may be conservative in some regions. In addition, only a few models include a fully coupled carbon-nitrogen cycle. Ocean models tend to be too coarse in resolution to properly represent important scales of motions and mixing, such as eddies and other mesoscale processes, and coastal boundary processes. Furthermore, their representation of ocean ecosystem processes and their sensitivity to climate change and other 25 stressors (e.g. ocean acidification, deoxygenation, etc.; Gruber, 2011; Boyd, 2011) is rather simplistic. 10,2013 Trends and drivers of regional sources and sinks of CO 2 over the past two decades S. Sitch et al. There is a need for detailed model evaluation and benchmarking, to reduce the uncertainty in the sinks in the land and ocean and particularly in how these sinks have changed in the past and how they may change in the future. For land ecosystems, a concerted effort is needed in the DGVM community to incorporate nutrient cycling, and land use and land cover change. For the oceans, models need to improve their 5 representation of unresolved physical transport and mixing processes, and ecosystem models need to evolve to better characterize their response to global change.