GCM characteristics explain the majority of uncertainty in projected 21 st century terrestrial ecosystem carbon balance

One of the largest sources of uncertainties in modelling of the future global climate is the response of the terrestrial carbon cycle. Studies have shown that it is likely that the extant land sink of carbon will weaken in a warming climate. Should this happen, a larger portion of the annual carbon dioxide emissions will remain in the atmosphere, and further increase global warming, which in turn may further weaken the land sink. We investigate the potential sensitivity of global terrestrial ecosystem carbon balance to differences in future climate simulated by four general circulation models (GCMs) under three different CO 2 concentration scenarios. We find that the response in simulated carbon balance is more influenced by GCMs than CO 2 concentration scenarios. Empirical orthogonal function (EOF) analysis of sea surface temperatures (SSTs) reveals differences among GCMs in simulated SST variability leading to decreased tropical ecosystem productivity in two out of four GCMs. We extract parameters describing GCM characteristics by parameterizing a statistical emulator mimicking the carbon balance response simulated by a full dynamic ecosystem model. By sampling two GCM-specific parameters and global temperatures we create 60 new “artificial” GCMs and investigate the extent to which the GCM characteristics may explain the uncertainty in global carbon balance under future radiative forcing. Differences among GCMs in the representation of SST variability and ENSO and its effect on precipitation and temperature patterns explain the majority of the uncertainty in the future evolution of global terrestrial ecosystem carbon in our analysis. We suggest that the characterisation and evaluation of patterns and trends in simulated SST variability should be a priority for the further d velopment of GCMs, in particular as vegetation dynamics and carbon cycle feedbacks are incorporated.

Abstract. One of the largest sources of uncertainties in modelling of the future global climate is the response of the terrestrial carbon cycle. Studies have shown that it is likely that the extant land sink of carbon will weaken in a warming climate. Should this happen, a larger portion of the annual carbon dioxide emissions will remain in the atmosphere, and further increase global warming, which in turn may further weaken the land sink. We investigate the potential sensitivity of global terrestrial ecosystem carbon balance to differences in future climate simulated by four general circulation models (GCMs) under three different CO 2 concentration scenarios. We find that the response in simulated carbon balance is more influenced by GCMs than CO 2 concentration scenarios. Empirical orthogonal function (EOF) analysis of sea surface temperatures (SSTs) reveals differences among GCMs in simulated SST variability leading to decreased tropical ecosystem productivity in two out of four GCMs. We extract parameters describing GCM characteristics by parameterizing a statistical emulator mimicking the carbon balance response simulated by a full dynamic ecosystem model. By sampling two GCM-specific parameters and global temperatures we create 60 new "artificial" GCMs and investigate the extent to which the GCM characteristics may explain the uncertainty in global carbon balance under future radiative forcing. Differences among GCMs in the representation of SST variability and ENSO and its effect on precipitation and temperature patterns explain the majority of the uncertainty in the future evolution of global terrestrial ecosystem carbon in our analysis. We suggest that the characterisation and evaluation of patterns and trends in simulated SST variabil-ity should be a priority for the further development of GCMs, in particular as vegetation dynamics and carbon cycle feedbacks are incorporated.

Introduction
Discussion about climate change uncertainties has tended to focus on climate sensitivity, i.e. the global mean warming induced by an increased atmospheric CO 2 concentration ([CO 2 ]), which differs among climate models (atmosphereocean general circulation models, AOGCMs, hereafter "GCM") depending on the assumed strength and sign of feedbacks that may dampen or amplify the direct radiative forcing of CO 2 through the greenhouse effect (Knutti and Hegerl, 2008). The terrestrial biosphere and surface layers of the oceans together sequester around 50-60 % of anthropogenic CO 2 emissions. The fate of these sinks in future decades constitutes a feedback that may significantly influence future climate change, but is not accounted for by many current GCMs (Canadell et al., 2007;Denman et al., 2007;Meehl et al., 2007b;Le Quéré et al., 2009).
Studies with dynamic global vegetation models (DGVMs) have shown a considerable spread among models in the future trajectory of terrestrial biosphere carbon balance when forced with identical output fields from climate models (Cramer et al., 2001;Sitch et al., 2008). Studies in which a single DGVM is forced with data from an ensemble of GCMs likewise show a considerable spread in carbon balance, in this case stemming from differences in the climate Published by Copernicus Publications on behalf of the European Geosciences Union.

1518
A. Ahlström et al.: GCM based carbon balance uncertainty simulated by different GCMs under the same CO 2 emission scenarios (Berthelot et al., 2005;Schaphoff et al., 2006). The spread among models -either DGVMs or the GCMs used to force them -may often be traced to contrasting simulated dynamics in a few specific regions, such as the Amazon Basin, where vegetation dieback and associated depletion of ecosystem carbon pools (Cox et al., 2000;Malhi et al., 2008) is simulated by some model combinations, but not by others (Berthelot et al., 2005;Schaphoff et al., 2006;Sitch et al., 2008). The impacts may be particularly dramatic when feedbacks of vegetation changes to the atmosphere through the carbon and hydrological cycles are accounted for by coupling the models to each other, simulated climate forcing the evolution of vegetation patterns and carbon pools, and viceversa (Friedlingstein et al., 2006). In short, the available evidence suggests that the uncertainties induced by knowledgeand data gaps on both climate sensitivity (as encapsulated by different GCMs) and carbon cycle response (as encapsulated by different DGVMs) are large and of similar importance in terms of the limitations they impose on our current ability to project changes in the global climate.
While studies have demonstrated and quantified the uncertainty in the evolution of global terrestrial carbon balance over the coming century stemming from differences among GCMs (Berthelot et al., 2005;Schaphoff et al., 2006;Morales et al., 2007), the model characteristics and simulated mechanisms underlying this uncertainty have not been systematically analysed. In this paper, we employ a statistical emulator of a global terrestrial carbon cycle model (DGVM) as a tool to characterize GCM-based uncertainty in projected 21st century ecosystem-atmosphere carbon exchange, identifying major model characteristics underlying such uncertainty, and suggesting priorities for the further development and improvement of the GCMs.

Overview of approach
Here we summarise the methodological steps used in this paper to investigate and quantify the uncertainty induced by GCM characteristics on simulated carbon fluxes. 2. DGVM simulations often result in widely different results in terms of carbon balance depending on the choice of forcing GCM. Based on the initial results from Step 1, we hypothesised that differences in SST variability account in part for the observed discrepancies in simulated carbon fluxes among DGVM simulations.
To evaluate this hypothesis, we applied EOF (empirical orthogonal function) analysis to simulated SSTs from each of the GCMs to separate the long-term trend (mainly warming) seen in all the GCM projections from other aspects of variability (Sect. 2.4). SST variation in space and time emerges as the most important component of this residual variability in all the GCM projections analysed. By correlating the dominating patterns of SST variability, as characterised by the second EOF mode in our analysis to carbon fluxes simulated by LPJ-GUESS, we analysed the impact of differences between the GCMs originating from differences in simulated SSTs and their influence -via circulation patterns -on local land surface climate.
3. Having established GCM-simulated warming and SST variability as important determinants (factors) of change in global terrestrial ecosystem carbon balance, as simulated by LPJ-GUESS, we wished to quantify the relative contribution of, and the residual variation not explained by, each factor. As global simulations with a full DGVM and the subsequent analysis of the considerable output dataset generated are expensive, we parameterised a statistical emulator mimicking the LPJ-GUESS results when forced by global temperature and [CO 2 ] as sole drivers (Sect. 2.5). To account for GCM-dependent differences in carbon balance response, proxies for two GCM-dependent factors were included as parameters in the emulator model. The first of these parameters, α, is a modifier on global gross primary production (GPP) and reflects the character or intensity of SST variability as simulated by different GCMs. The second parameter, γ , is a scalar which translates the anomaly of global mean surface temperature in a GCM climate projection relative to a modern baseline to the corresponding anomaly in global land temperatures.
4. In the last step we employed a factorial simulation approach using the statistical carbon cycle model emulator to attribute and partition variability (uncertainty) in future carbon balance changes to variability in the main global drivers and uncertainty stemming from differences in GCM behaviour, as encapsulated by α and γ (Sect. 2.5). An analysis of variance (ANOVA) was applied to the output of 192 synthetic simulations with the emulator model, spanning the observed space in the drivers (CO 2 concentration and associated temperature change from the original ensemble of GCM projections), α and γ .

Carbon cycle model
LPJ-GUESS is an individual-based, globally-applicable model of vegetation dynamics and biogeochemistry. It combines features of the Lund-Potsdam-Jena (LPJ) DGVM (Sitch et al., 2003;Gerten et al., 2004) with a more detailed treatment of plant resource competition and demography based on simulated interactions of plant individuals co-occurring in patches . Vegetation is represented by a mixture of plant functional types (PFTs, 11 in this study), differentiated by bioclimatic limits, growth form, phenology, photosynthetic pathway (C 3 or C 4 ) and life history strategy. LPJ-GUESS has been evaluated extensively and exhibits comparable skill to other approaches and models in reproducing observed temporal and spatial variation in large-scale vegetation patterns and ecosystematmosphere carbon exchange (Piao et al., 2013). Simulated carbon and evaporative fluxes have been compared to ecosystem flux measurements (Morales et al., 2005;, forest inventory data and site-based measurements of NPP, leaf area index and biomass, spanning many of the world's biomes Zaehle et al., 2006;Smith et al., 2008Smith et al., , 2011Tang et al., 2010Tang et al., , 2012. See Smith et al. (2001) for a detailed description of the model. The version used in this study includes the updates detailed in Hickler et al. (2012), and the same PFT set and configuration as described in Ahlström et al. (2012a). Our analysis was based on the simulated ecosystem carbon fluxes output by LPJ-GUESS, specifically gross primary production (GPP, i.e. the annual sum of canopy-level net photosynthesis for the vegetation in a grid cell), net primary production (NPP, i.e. GPP minus autotrophic [plant] respiration), ecosystem respiration (Er, i.e. the sum of, autotrophic and heterotrophic respiration), and biomass burning through wildfires (Fire). Net Biome Production (NBP) is the small difference between GPP and the release fluxes, Er and Fire. The total terrestrial carbon pool (Cpool) is the sum of the carbon residing in standing biomass and the organic layers of the soil. The cumulative sum of NBP over time is analogous to the change in Cpool.

Climate data
The LPJ-GUESS simulations were initialised with data from the CRU TS3.0 (Climatic Research Unit time series) observed climate database (Mitchell and Jones, 2005) covering the period 1901-2000. From 2001, we forced LPJ-GUESS with simulations by four different GCMs, each with three atmospheric [CO 2 ] pathways based on the SRES emission scenarios (A2, A1B, B1) (Nakicenovic et al., 2000), giving in total twelve simulations. The data were obtained from the World Climate Research Programme's (WCRP) Coupled Model Intercomparison Project phase 3 (CMIP3) multimodel dataset (Meehl et al., 2007a). Four representative GCMs were chosen: CCSM3 (Community Climate System Model version 3) (Collins et al., 2006), UKMO-HadCM3 (United Kingdom Meteorological Office Hadley Center Coupled Model version 3) (Gordon et al., 2000), IPSL-CM4 (Institut Pierre Simon Laplace Coupled Model version 4) (Marti et al., 2005) and ECHAM5/MPI-OM (European Centre-Hamburg model version 5/Max Planck Institute Ocean Model) (Roeckner et al., 2003). Further details of the simulation set up and forcing are provided in the Supplement, Text S1.

EOF analysis
We focused on sea surface temperature (SSTs) as an overall indicator of those aspects of GCM-simulated climate of importance in terms of impacts on global ecosystem carbon balance. SST records for recent decades clearly reflect the general trend in global temperatures (Dai, 2013). Additionally, variability in SSTs has been shown to have a large influence on regional precipitation trends (Hoerling et al., 2009) and low latitude precipitation in general (Dai, 2006). ENSO is the dominant process of SST variability. ENSO has been shown to be the strongest controlling factor for global precipitation variability (Dai et al., 1997).
We characterized the relationship between SST patterns from a given GCM and NPP simulated by LPJ-GUESS when forced by that GCM by performing a singular value decomposition (SVD) analysis of the global SST field. SVD is the generalization of the diagonalisation process that EOF analysis does for square matrixes, but since the SST matrix is not square we used SVD to calculate the EOFs (Uvo and Berndtsson, 1996). EOF analysis, sometimes also referred to as principal component analysis (PCA; although definitions vary), can be used to capture and illustrate important spatial and temporal patterns of geophysical fields (e.g. Wallace et al., 1993;Uvo et al., 1998;Quadrelli and Wallace, 2004). Central to EOF analysis methods is the concept of transforming the original data via composite expressions (the EOFs or "modes") ordered by the proportion of the total variability in the original dataset they explain. The information encapsulated in a particular mode may be visualised as a spatial pattern (often a map) indicating the spatial loading of that mode, e.g. where in space the variability is centered, a time series representing variability or a trend in the mode, or a scalar value (eigenvalue, singular value) representing the amount of variability explained by that mode. Hereinafter the term "mode" in this paper denotes the sense outlined above.
The time series of the different modes, which are vectors the same length (here 80 yr) as the original data, are orthogonal (uncorrelated, independent), and often but not always reflect the influence of different physical processes on spatial and temporal variability in the original data.
We derived EOF modes from the standardised annual SST field (transformed to a mean of 0 and standard deviation of 1) of the SRES A2 simulations for each of the four GCMs considered in our study (Sect. 2.3). Due to corrections (see Text S1 in Supplement) of the variability and trend applied to the climate data between 2001 and 2020, we confined the analysis to data from the period 2020-2099. We chose to present the results of the EOF analysis as correlation maps, both for SSTs and for NPP. The correlation maps share features with the spatial loading patterns for the SST modes. Because they are unitless they may be directly compared with the corresponding NPP maps.
For the three first modes, we calculated the correlation between the resultant time series and both the original SSTs and the simulated NPP in every gridcell (only first and second modes are presented in this paper). The resulting spatial pattern of correlations between modes and original GCMsimulated SSTs shows how much of the variability of the GCM-simulated SSTs is related to the interannual variability described by each mode. The pattern of correlations between modes and NPP may be taken to reflect the impact of that aspect of SST variability captured by each mode on the NPP of different regions. The statistical significance of the calculated correlations was estimated using a bootstrap technique with 500 repetitions (Wilks, 2006).
The fraction of the total SST variation explained by a given mode was quantified as the squared covariance fraction (SCF): where l i is the singular value of the i-th mode. The analysis was repeated for NBP, and for precipitation and temperature to better understand the underlying driver of the resulting NPP and NBP patterns.

Partitioning uncertainty in future carbon balance
We wished to partition uncertainty in terrestrial ecosystem carbon balance change due to temperature, [CO 2 ], and GCM characteristics. A model statistically emulating LPJ-GUESS was defined as a tool to enable uncertainty propagating from the four GCMs investigated to be extrapolated to a wider "space" of GCM behaviour by deterministic parameter resampling. To characterise the influence of GCM-dependent SST variability, we defined a parameter, α. To relate α to the GCMs SST variability, we needed a scenario-independent measure of SST variability. The global mean warming trend is well separated from other variability -not directly due to increased radiative forcing -in the EOF analysis of the SSTs of the SRES A2 simulations. The inverse of the SCF of the first mode could therefore be used as a proxy for the relative importance of variability not directly related to climate change. The global warming is, by contrast, not always well separated in the SSTs of the lower emission scenarios, A1B and B1, and the SCF was not found to be a useful proxy of SST variability for these scenarios. Therefore, to be able to compare the α values from the statistical carbon cycle emulator with a scenario-independent measure of SST variability we calculated the global average standard deviation of all the detrended SST time series. For each grid cell we removed the 2015-2085 trend (based on a 30-yr moving average of the global SST time series) using linear regression. We then computed the standard deviation (SD) over the period 2015 through 2085 in the resultant trend-free time series, averaging the SDs among all grid cells to produce a global average SST variability for each GCM simulation used except for CCSM3-B1 (SSTs for the latter were not available in the CMIP3 database).
Simple polynomials were fitted to global sums of GPP, Er, Fire and Cpool simulated by LPJ-GUESS under GCM forcing, with global average land temperature and atmospheric [CO 2 ] as the predictor variables. For GPP (in Pg C), the resultant function has the form (2) +β 6 Cpool t + β 7 Cpool 2 t + α, where t indicates the time step in years, β i are coefficients fitted to the data (see below), T ( • C) is the global land temperature and Cpool (Pg C) is the total terrestrial carbon pool from Eq. (4). The last term in the equation, α, is a simulationspecific parameter (Pg C) that is allowed to vary between the 12 simulations used for the calibration. As outlined above, it represents factors influencing GPP not captured by land temperature and [CO 2 ], e.g. regional patterns of temperature, variability and variables not included in the forcing data for the statistical emulator, such as precipitation and shortwave radiation. α was set to zero during the historical period , it was then linearly interpolated until 2020, similar to the forcing climate data. After 2020 it was held constant until the end of each simulation.
Equation 3 describes the sum of the global carbon fluxes from ecosystem respiration and wild fires (Pg C) as a function of global land temperature, GPP from Eq. (2) and Cpool from Eq. (4).
Er + Fire t = β 8 + β 9 T t + β 10 T 2 t + β 11 GPP t + β 12 Cpool t . (3) In Eq. (4) the total carbon pool of the next time step (year) is calculated as the present year's carbon pool plus the current year's net carbon exchange, NBP.
The coefficients β 1−12 (Eqs. 2 and 3) were estimated simultaneously by minimizing the sum of the square of the errors (SSE) between the resulting GPP, Er + Fire, Cpool and the corresponding LPJ-GUESS simulation results; see Table S1 for parameter values in the Supplement. Cpool was initialised with the total carbon pool of year 1901, the first year after the initial ("spin-up") phase of the relevant LPJ-GUESS simulation (see Supplement, Text S1). We performed a separate validation study by excluding one scenario from each GCM in the calibration of the statistical emulator. This showed (Supplement, Text S2 and Figs. S1 and S2) that the emulator is robust even after exclusion of about one third of the calibration data, and that there are no changes that would affect the conclusions drawn in this paper if only a subset of the data were used in the calibration. Because the estimated parameters (α) were crucial for the analysis in this paper, however, we use all available data in our calibration.
To better separate the CO 2 signal from the temperature signal we used the full length of the climate data series, including the "stabilization" period after 2100 when [CO 2 ] were held constant, but temperatures continued to evolve. The length of the stabilization period varied between the datasets of the A1B scenario, resulting in a slightly uneven weighting of the four GCMs.
The emulator uses the area-weighted average of land temperatures, T , as a predictor variable. However, land temperatures may evolve at a different rate than the global mean surface temperature -which also accounts for temperatures over the oceans, inland water bodies and ice sheets -in GCM simulations. We calculated a factor, γ , that expresses the degree of amplification in the temperature trend over land relative to the global mean surface temperatures, with 1961-1990 mean temperatures as a baseline: where LT t is the mean land temperature change by year t compared with the 1961-1990 average. GT t is the global mean surface temperature change compared with the 1961-1990 average. Now, land temperature T in the replacement model (Eqs. 2 and 3) can be calculated from global mean surface temperature following where T 61−90 represents the average 1961-1990 land temperature in the CRU dataset.
To partition the variation in future terrestrial carbon balance suggested by our simulations to underlying factors in the GCMs and [CO 2 ] pathways providing the forcing, we adopted a permutation approach, using the statistical emulator (Eqs. 2-4) to generate 192 synthetic time series, corresponding to all possible combinations of the main independent forcing factors encompassed by our analysis, namely the three [CO 2 ] pathways, the four corresponding global temperature change time series ( GT) (one simulated by each GCM), four GCM average values of the simulation-specific parameter α in Eq. (2), and four GCM average values of the land-to-global warming ratio, γ , from Eq. (5).
For each [CO 2 ] pathway (A2, A1B, B1) we thus created 64 (4 3 ) combinations of α, γ and GT. Four of these correspond to the original models, while the remaining 60 may be thought of as "artificial" models with different, plausible GCM characteristics. Applying the ensemble of the 60 "artificial" GCMs plus the original four in combination with the three [CO 2 ] pathways produced a suite of 192 plausible carbon balance trajectories.  To attribute variation in the resultant carbon balance time series to underlying forcing factors, we performed analysis of variance (ANOVA; Draper and Smith 1998) with α, γ , and [CO 2 ] pathway (A2, A1B, B1) as independent factors. The ANOVA gives information on which factor induces the largest spread in simulated total carbon pool by the statistical emulator. It separates the total sum of squares in the simulated change in total carbon stock from 2000-2099 into sum of squares due to each of the three factors above and a residual term that contains effects due to interactions between factors and an unexplained component, such as temperature. Since the temperature depends on the specific CO 2 scenario, it is non-trivial to quantify the effect of different temperatures, forcing us to limit the ANOVA to the three factors α, γ , and CO 2 .

Carbon cycle simulations
The LPJ-GUESS simulations forced by different GCMs and [CO 2 ] pathways result in a wide range of trajectories for Cpool (Fig. 1). Simulations forced by climate data from the same GCM under different CO 2 emission scenarios tend to cluster, indicating that differences between GCMs exert a stronger control on Cpool than [CO 2 ] pathways. By 2099, the largest difference between the simulations, 247 Pg C, is found between the simulation forced by CCSM3 under the A1B and HadCM3 under the B1 emissions scenario. Regardless of emissions scenario, the largest contrast is found between the simulations forced by CCSM3 and HadCM3.

GCM climate-carbon cycle relationships
The EOF time series of the first two modes, the SST spatial correlation patterns and the NPP correlation patterns for the time period 2020-2099 are illustrated in Fig. 2. The first mode is closely related to the global temperature trend and variability as indicated by the overall strong correlation over all the world's oceans (Fig. 2a). The temporal variability and trend of the first mode show similarities with the standardized global land temperatures (Fig. 2e). The SCF differs between the GCMs. A high SCF in the first mode implies a relatively uniform pattern of SST variability across the globe, as seen for CCSM3, while a lower SCF in the first mode, as seen for HadCM3 and ECHAM5, implies that differences in SST variability in different areas of the world's oceans more strongly contribute to the total variability globally.
The second mode of variability (Fig. 2b) explains around an order of magnitude less of the total variance of the simulated SSTs for all four GCMs, and is generally characterized by El Niño-Southern Oscillation (ENSO) like patterns, domi-nated by variability in the central tropical Pacific. The ENSO patterns are more prominent in the SSTs of ECHAM5 and HadCM3 compared with the other GCMs, while for CM4 they are more pronounced in the third mode (not shown here) than in the second mode, indicating that other variation, not primarily associated with ENSO, is influencing the second mode of variability for this GCM. Time series for the first two modes (Fig. 2e) show, as expected, a uniform up-going trend in the first mode consistent with global warming response to rising greenhouse gas concentrations, while the second mode exhibits continuous fluctuations with no obvious trends for any of the models.
The third and fourth columns (Fig. 2c-d) show the correlation between the modes and simulated NPP. Global temperature increase, together with the covarying [CO 2 ] (closely related to the first mode of variability of the simulated SSTs) are associated with increasing NPP in most areas, except in some tropical regions, where a negative correlation is seen. America and western North Africa see a decreased NPP resulting from the negative impacts of decreased precipitation and/or increased evapotranspiration under higher temperatures on plant water relations (Supplement Fig. S3). The correlation patterns between the second mode and NPP (Fig. 2d) illustrate the impact of ENSO-related regional climate patterns on land ecosystems, particularly water relations. For HadCM3, ENSO leads to pronounced decrease of NPP in the tropics and Australia and increased NPP in western North America and the Middle East. In ECHAM5 the negative impact is more pronounced in Australia but less pronounced in the tropics. The regions that show increased NPP are similar for both GCMs. The direct cause of the NPP declines in HadCM3 and ECHAM5 is decreased precipitation and increased temperature (Supplement Fig. S3).
CCSM3 shows a similar but considerably weaker correlation pattern compared to ECHAM5, while CM4 shows no strong patterns, probably because ENSO dynamics for that GCM are expressed mainly in the third mode. Correlation patterns between SST modes and NBP, shown in Fig. S4 SST Standard Deviation °C α (Pg C/y) 1 Figure 4. Relation between the simulation-specific parameter α and glo 2 driving role of vegetation productivity in the carbon cycling of ecosystems as a whole.

GCM-dependent uncertainty in carbon balance
GCM-simulated global mean land temperature and the [CO 2 ] pathway of the underlying emission scenario explain part, but not all, of the variation between simulations in the evolution of the terrestrial carbon pool over the 21st century, as simulated by LPJ-GUESS. The remaining, unexplained, source of variation, signified by α in Eq. (2), encapsulates the regional patterns of carbon balance and their evolution over time, specific to each GCM × emissions scenario combination (Fig. 3). The displacement between the black (simulation-specific fit for α) and grey (α set to the average of the twelve simulation-specific values) curve in each frame of Fig. 3 represents the effect of this unexplained variation in each individual simulation.
The EOF analysis reveals differences in the impact on simulated NPP and NBP, induced by differences among GCMs related to SST variability and associated weather patterns (see Fig. 2). We found a strong relationship between α and global SST variability (Fig. 4). The α values cluster according to GCMs and not according to CO 2 emission scenarios, signifying a strong GCM dependence for this parameter. The CCSM3 simulations show the largest uptake of carbon (Fig. 1). The corresponding α values are also the least negative (Fig. 4). The CM4 simulations group in the middle, and show a decrease of about 7 Pg C in GPP per year compared to the historical simulation. ECHAM5 and HadCM3 show larger SST variability compared to CM4 and CCSM3, and the simulations forced with these models also have the most negative α values of all the simulations. All simulations have an α value lower than the historical simulation (for which α = 0).

Global versus land temperature
The ratio of land to global warming γ varied from 1.26 to 1.46 between the 12 datasets for 2001-2009 (Fig. 5). CCSM3 shows the lowest γ (least amplification of warming over land) among all scenarios while HadCM3 shows the largest γ among all scenarios. A single GCM-specific factor seems to be enough to explain most variation in γ as long as the climate forcing is increasing. In the longest dataset used here, the CM4 A1B, top, middle column in Fig. 5, where [CO 2 ] has been held constant since 2100, land temperatures have started to diverge from the constant warming amplification seen during the 21st century, approaching the global temperatures. The range in A2 global temperature (2095-2099) is 0.27 • C, while the range in land temperatures is more than double this at 0.57 • C. The global average temperature change, GT, (2095GT, ( -2099 in the HadCM3 model's A2 scenario is the lowest among the four GCMs (3.9 • C), while the land temperature, LT, for the same period, GCM and scenario is the highest among the GCMs (5.7 • C).

Partitioning uncertainty in carbon balance
Using the statistical carbon balance emulator to sample across 192 combinations of GCM and scenario characteristics (see Methods), we came up with the distribution of carbon balance trajectories shown in Fig. 6. For this ensemble of [CO 2 ] pathways and GCMs, and given the carbon cycle representation in LPJ-GUESS, total terrestrial carbon stocks may increase or decrease by the end of the 21st century, but are more likely to increase. Note that the percentiles and the medians shown in Fig. 6 are partly a result of the sampling. It is apparent that much of the uncertainty represented by the spread among trajectories propagates from different characteristics and behaviour of GCMs. According to the ANOVA analysis, 91 % of the variability among trajectories could be explained by the two GCM-specific factors α (87 %) and γ (5 %). The choice of CO 2 emission scenario explains just 2 % of the variation, while 7 % remains unexplained (Fig. 6).

Discussion
The LPJ-GUESS simulations forced by the ensemble of GCMs and CO 2 emission scenarios chosen for this study result in a considerable spread in the future evolution of carbon balance. Our results demonstrate that the majority of this spread can be traced to a GCM-specific parameter (α) closely related to global SST variability.
As seen in Fig. 4, all the GCMs show a negative α value (α was set to 0 for CRU). One interpretation could be that all GCMs show increasing climate variability or other cli- mate characteristics that negatively influence GPP. However, the offset of the α values in comparison to CRU can also be explained by the offset between the GCMs and CRU as a result of different trends since the 1961-1990 climatology. Additionally, as a result of station data limitations, around 13 % of the CRU gridcells used in this study show at least one consecutive 10-yr period with no interannual variability (see Supplement Fig. S5). Although the importance of these episodes in the CRU data for the simulated carbon balance is unknown, it is difficult to draw conclusions about the α values.
Our results are based on four GCMs. While these were chosen to provide a representative range of climate responses to future emissions, we do not know whether the results or the overall conclusions would have been different, had we chosen a larger sample of CMIP3 GCMs. If the relationships found apply to other GCMs (or later versions of the four we have explicitly considered) remains to be tested. The carbon balance spread shown is this paper is likely to increase if uncertainties surrounding land use and land use change were Fig. 6. Results of the 192 replacement model simulations using "artificial GCM" input. The total spread of the A2 scenario simulations in total terrestrial carbon pool (Cpool) is illustrated in red. Blue shows the A1B scenario and green the B1 scenario. The boxplots show the median, 25th and 75th percentile, and maximum and minimum values of the 2099 total terrestrial carbon pool (Pg C), for each CO 2 scenario and all simulations. The pie chart shows the proportion of 2099 total carbon pool variability explained by, α, γ , CO 2 scenarios and residual unexplained variation (res.). also considered. Different land use projections, as well as different treatment of land use and land use change in DGVMs and carbon cycle models, constitute a large source of uncertainty as to future carbon balance (e.g. Zaehle et al., 2007). Still, our results provide some quantification of GCM discrepancies and their impact on carbon balance as simulated by our model.
By applying different GCMs and scenarios to force a single ecosystem model, LPJ-GUESS, we have focused on the uncertainties arising from differences in the climate as projected by the GCMs. Previous studies show that different DGVMs can project substantially different carbon balance in response to the same forcing (Cramer et al., 2001;Sitch et al., 2008;Piao et al., 2013). This carbon cycle model-dependent aspect of the overall uncertainty has not been addressed in this study. However, LPJ-GUESS shows comparable skill and behaviour to other DGVMs and may be regarded as representative of DGVMs as a group in terms global carbon balance response to climate and [CO 2 ] forcing. For example, a recent study comparing ten DGVMs with each other and with independent datasets (Piao et al., 2013) shows that LPJ-GUESS predicts present day global GPP in the middle of the range of other models, in agreement with an observationbased estimate, and exhibits comparable sensitivity to precipitation as suggested by upscaled ecosystem flux measurements.
When applying the statistical emulator with combinations of the GCM-specific parameters and variables we assume that they are independent. When averaging γ and α across GCMs (n = 4), they show no significant correlation (r = −0.93, n.s.) but when we use all simulations (n = 12), they show a significant correlation (r = −0.92, P < 0.00005). However, this high correlation could be a coincidence result-ing from the specific selection of GCMs. Amongst many factors and processes it is plausible that GCM differences, for example in the partitioning of latent and sensible heat fluxes between the surface and atmosphere can influence both α and γ . Other possible differences between GCMs, such as differences in ocean overturning or the melting of sea ice affecting ocean warming, potentially influencing γ are likely to have smaller effect on α. ENSO is the dominant determinant of global precipitation variability (Dai et al., 1997). A prominent effect of warm ENSO events is negative precipitation anomalies in large parts of the tropics, but also significant precipitation anomalies in other regions can be traced back to ENSO (Dai and Wigley, 2000). Patterns of correlation between pacific SST variability and precipitation reminiscent of ENSO have been found to emerge in GCM simulations. The EOF analysis of sea surface temperatures presented here reveals carbon cycle patterns similar to those found by Dai (2006) when analysing rainfall patterns and tropical SST variability in HadCM3, ECHAM5 and CCSM3, amongst other GCMs. Although the "strength" of representation of ENSO differs between GCMs, a strong dependency of low latitude precipitation and tropical SST variability is typically simulated on interannual and longer timescales (Dai, 2006).
The effects of droughts on carbon fluxes in the Amazon Basin have been debated. Satellite-based studies have reported a "green up" in the Amazon rainforest during dry conditions as a result of increased incoming solar radiation when cloud cover decreases (Huete et al., 2006;Saleska et al., 2007). However, field studies have reported increased tree mortality and carbon loss during severe droughts (Nepstad et al., 2007;Brando et al., 2008;Phillips et al., 2009), resulting from decreased plant available water and/or heat stress (Toomey et al., 2011). Recent studies employing a satellitebased model and a DGVM have suggested a strong NPP dependency on plant available water and susceptibility to droughts in tropical regions over the last decade (Zhao and Running, 2010;Ahlström et al., 2012a). Impacts of droughts on tropical vegetation have been demonstrated to be a major uncertainty in carbon balance in modelling studies applying several GCMs to force ecosystem models (Berthelot et al., 2005;Schaphoff et al., 2006;Ahlström et al., 2012b). The negative NPP and NBP anomalies found have either been explained by decreased precipitation (Schaphoff et al., 2006) or increased temperatures (Berthelot et al., 2005), leading to increased plant respiration and/or exerting an indirect negative effect on plant available water by increasing the atmospheric demand for water vapour. The results presented here concur with previous studies in attributing the majority of the uncertainties in future global terrestrial carbon cycle to decreased plant available water mainly in the tropics, a result of reoccurring droughts, induced by SST variations accompanied by more static regional differences in future climatology. Our results point to a marked dependency of the future development of global terrestrial ecosystem carbon balance on climate characteristics, particularly SST variability and its impact on weather patterns, simulated differently among different GCMs. A further GCM characteristic, the degree of enhancement of warming over land relative to global warming generally, accounts for the second largest proportion of uncertainty in carbon balance response. Uncertainty stemming from the choice of CO 2 emission scenario is much less marked. Finally, our results suggest that improved representations of ENSO dynamics and low-latitude precipitation patterns are important to narrow the uncertainties in future climate change projections using ESMs or uncoupled DGVMs forced by GCM-simulated climate projections.