A model-based constraint on CO2 fertilisation

We derive a constraint on the strength of CO 2 fertilisation of the terrestrial biosphere through a “top-down” approach, calibrating Earth system model parameters constrained by the post-industrial increase of atmospheric CO 2 concentration. We derive a probabilistic prediction for the globally averaged strength of CO 2 fertilisation in nature, for the period 1850 to 2000 AD, implicitly net of other limiting factors such as nutrient availability. The approach yields an estimate that is independent of CO 2 enrichment experiments. To achieve this, an essential requirement was the incorporation of a land use change (LUC) scheme into the GENIE Earth system model. Using output from a 671-member ensemble of transient GENIE simulations, we build an emulator of the change in atmospheric CO 2 concentration change since the preindustrial period. We use this emulator to sample the 28-dimensional input parameter space. A Bayesian calibration of the emulator output suggests that the increase in gross primary productivity (GPP) in response to a doubling of CO2 from preindustrial values is very likely (90 % confidence) to exceed 20 %, with a most likely value of 40– 60 %. It is important to note that we do not represent all of the possible contributing mechanisms to the terrestrial sink. The missing processes are subsumed into our calibration of CO2 fertilisation, which therefore represents the combined effect of CO2 fertilisation and additional missing processes. If the missing processes are a net sink then our estimate represents an upper bound. We derive calibrated estimates of carbon fluxes that are consistent with existing estimates. The present-day land–atmosphere flux (1990–2000) is estimated at −0.7 GTC yr−1 (likely, 66 % confidence, in the range 0.4 to −1.7 GTC yr−1). The present-day ocean–atmosphere flux (1990–2000) is estimated to be −2.3 GTC yr−1 (likely in the range−1.8 to−2.7 GTC yr−1). We estimate cumulative net land emissions over the post-industrial period (land use change emissions net of the CO 2 fertilisation and climate sinks) to be 66 GTC, likely to lie in the range 0 to 128 GTC.


Introduction
Experimental evidence almost without exception shows a stimulation of leaf photosynthesis when plants are exposed to elevated CO 2 (Koerner, 2006). In addition to this direct effect on photosynthesis, the short timescale physiological effect of reduced stomatal opening increases water-use efficiency and additionally increases the efficiency of photosynthesis (Field et al., 1995). The combined effects, which we group under the label of CO 2 fertilisation, act as a negative feedback for anthropogenic carbon emissions. Increased concentrations of atmospheric CO 2 lead to increased photosynthesis and more efficient drawdown, transferring some fraction of these emissions from the atmosphere into terrestrial carbon pools. However, the strength of the fertilisation effect is poorly quantified, especially under natural conditions. Some studies have failed to detect a measurable effect in nature, while others suggest that any effects may be short-term as CO 2 is only one of a number of potentially limiting factors on plant growth (Koerner, 2006). In addition to the implications for future vegetation and crop growth, improved quantification of the globally integrated effects of CO 2 fertilisation in nature is crucial to reduce the uncertainties in carbon-cycle projections from process-based models. In response to SRES A2 forcing, projections of 2100 CO 2 from 11 C 4 MIP models ranged from 740 to 1030 ppm, the largest source of uncertainty coming from the terrestrial response to elevated CO 2 (Friedlingstein et al., 2006). The concept of CO 2 fertilisation has been well demonstrated in many studies under controlled conditions. The most extensive of these experiments applied the FACE (freeair CO 2 enrichment) approach to both coniferous and deciduous trees in four separate sites (Norby et al., 2005). This study measured a 23 % increase of NPP following a doubling from preindustrial CO 2 concentrations. However, there remains doubt as to whether the effect is realised under natural conditions as a result, for instance, of nitrogen limitation (Norby et al., 2010) or temperature limitation in boreal forests . Girardin et al. (2011) failed to measure an effect in boreal forests, setting an upper limit on the CO 2 fertilisation effect at the detection limit, estimated at 14 %. In a model-based study of the period 1973-2004, the inclusion of nitrogen limitation was found to reduce the strength of the global carbon-cycle feedback to just 27 % of that with only carbon cycling (Bonan and Levis, 2010). Quantification is made even more difficult because elevated CO 2 levels are associated with increased temperatures, increasing respiration rates and opposing any change due to fertilisation. Hickler et al. (2008) validated the dynamic vegetation model LPJ- GUESS (Smith et al., 2001) against the FACE experiments of Norby et al. (2005) and applied the model to a global simulation that projected substantial regional differences in the response of forests to a CO 2 doubling, ranging from a 15 % increase in NPP in boreal forests to a 35 % increase in tropical forests. Over a thousand scientific articles have been published on the topic, yet there is no clear consensus (Koerner, 2006). In short, it is far from straightforward to extrapolate empirical evidence from controlled experiments to the global scale.
We here derive, for the first time, a probabilistic calibration of the CO 2 fertilisation effect (implicitly net of nitrogen limitation) through a "top-down" approach. In essence, we force the intermediate complexity Earth system model GE-NIE with historical emissions and land use changes and perform a Bayesian calibration of the modelled CO 2 fertilisation strength, constrained by the post-industrial change in the atmospheric CO 2 reservoir. A great deal of research has focused on the "missing carbon sink" (Broecker et al., 1979), the presumed uptake of carbon by the terrestrial biosphere that is required to reconcile the accumulation of CO 2 in the atmosphere with historical anthropogenic emissions and oceanic uptake. Earth system models implicitly assume that climate effects, reforestation, CO 2 fertilisation and (in some models) nitrogen deposition are the dominant processes controlling the terrestrial sink. However, it is important to note that some of the possible contributors to the terrestrial sink (see Denman et al., 2007) are not modelled, such as woody encroachment (which may itself be CO 2 driven, Bond et al., 2003) or changes in fire management and agricultural practices. In the following analysis, these un-modelled processes are effectively subsumed into the calibration of CO 2 fertilisation. It is also essential to understand that even if the rate of CO 2 fertilisation in the past could be determined precisely, its value in the future will be different as a result of structural errors that cannot be determined on the basis of past observations. It is therefore essential to account for structural error, but its magnitude must be based on expert judgement.
A number of approaches have been applied to constrain the residual terrestrial sink, all of which are associated with significant uncertainty (Denman et al., 2007). The approach developed here is most closely related to the top-down approach of single deconvolution (Siegenthaler and Oeschger, 1987) in that both use an ocean model to constrain carbon uptake and solve for the residual terrestrial source/sink. The fundamental difference here is that we are using the carbon budget constraint to calibrate a terrestrial carbon model, thus constraining the dynamics of the vegetation rather than just the residual flux.
In order to perform this calibration it was first necessary to incorporate a representation of land use change (LUC) into GENIE, described in Sect. 2 and validated in Sect. 4. Emissions from land use change constitute a substantial (and poorly quantified) source of anthropogenic emissions. We perform a 671-member ensemble of simulations with GE-NIE (Sect. 3). We use these simulations to build an emulator of atmospheric CO 2 concentration that we apply to probe the high dimensional parameter space more fully (Sect. 5). We then apply a Bayesian analysis to emulator outputs in order to constrain CO 2 fertilisation, calibrating vegetation parameters against observed CO 2 (Sect. 6.1), and apply the parameter calibrations to the simulated ensemble to generate probability density functions (pdfs) of modelled outputs (Sect. 6.2).
The reduced complexity of GENIE is ideal for performing the large number of simulations required for such an analysis. All vegetation is treated as a single plant functional type that responds to the changing CO 2 concentration according to the same functional dependence. In effect we are constraining a globally averaged rate parameter using global average observations and large ensembles of simulations that are nevertheless constrained to generate reasonable spatial simulations of climate and vegetation states. An equivalent procedure using a model that resolved multiple vegetation classes, while in principle being subject to reduced structural error variance, would require elicitation of combined prior estimates and uncertainty ranges for a large number of input parameters and demand substantial computation resource.
The GENIE calibrations are supplemented with three simulations of LPJmL (Bondeau et al., 2007). The purpose of these simulations is to quantify errors arising from simplifying assumptions made in the GENIE LUC implementation, and to validate the GENIE simulations with results from a complex, bottom-up modelling framework. P. B. Holden et al.: A model-based constraint on CO 2 fertilisation 341 2011) comprises the 3-D frictional geostrophic ocean model GOLDSTEIN (at 36 × 36 × 16 resolution) coupled to a 2-D Energy Moisture Balance Atmosphere and a thermodynamic-dynamic sea-ice model. Vegetation is simulated with ENTS, a dynamic model of terrestrial carbon storage . Ocean chemistry is modelled with BIOGEM  and is coupled to the sediment model SEDGEM at a 36 × 36 resolution . More detail of the specific model configuration can be found in Holden et al. (2012). The principal equations of ENTS are summarised in the Appendix.
The major extension required to ENTS was an implementation of the effects of land use change (LUC), induced by human management of land, on the carbon cycle and climate. We call this revised model ENTSML. For the modelling of LUC, each grid cell is apportioned between natural vegetation and cultivated vegetation. No distinction is made between crop and pasture. In the ENTS formulation , the vegetation state within a grid cell is defined by two state variables, vegetative carbon density C v and soil carbon density C s . In ENTSML no new state variables are added, instead we reinterpret the existing variables to be C vN (the natural, or "potential", vegetation carbon density) and C sT (the total soil carbon density, averaged over both natural and cropped regions). While this allows an adequate representation of deforestation, a weakness of the approach is that by only modelling potential vegetation, we cannot properly capture the timescales of reforestation (discussed in the text that follows Eq. 6). However, the alternative would be a substantially more complicated scheme that separately modelled the vegetation dynamics of natural and cultivated regions. The vegetative carbon density in cultivated regions is defined to be zero, equivalent to the simplifying assumption that cultivated vegetation is instantaneously harvested (crops) or grazed (pastures) and released into the atmosphere.
In ENTSML, as in ENTS, natural vegetative carbon density is given by Appendix A provides the expressions for the rates P (photosynthesis, Eq. A1), R vN (vegetation respiration, Eq. A2) and L (leaf litter, Eq. A3). These are unchanged from the ENTS formulation , but are expressed in terms of the potential vegetation density C vN . The dependence of photosynthesis on CO 2 is given by the Michaelis-Menten relationship where k 19 normalises the fertilisation response to unity at CO 2 = 278 ppm. The parameter k 14 is varied in the ensemble analysis that follows. The parameter k 13 is held constant at 29 ppm throughout. Although the expression for CO 2 fertilisation parameterises the uncertain saturating increase in gross primary productivity (GPP) under elevated CO 2 , the simplified moisture balance formulation does not capture the changes in evapotranspiration, soil moisture and run-off due to the physiological effect (Leipprand and Gerten, 2006;Betts et al., 2007).
In ENTSML, soil carbon density is defined as the average over both natural and cultivated regions. In naturally vegetated regions the rate of change of soil carbon density is given, as in ENTS, by where the temperature-dependent soil respiration rate R s is given in Appendix A (Eq. A4).
Harvesting reduces the input of carbon to the soil in arable regions. Agricultural processes such as tillage, which increases oxidation rates, may further aggravate this loss of soil carbon. Arable soils are estimated to be losing carbon in all European countries at an estimated average rate of 70 gC m −2 yr −1 (Janssens et al., 2005). To capture this effect in ENTSML, the leaf litter input to the soil in cultivated regions is derived from the natural leaf litter rate, reduced by a fraction k C (or the "land management parameter"). Previous modelling studies have applied a similar approach: Olofsson and Hickler (2008) applied a 33 % reduction in the fraction of the litter decomposition flux that is transferred to the soil; Stocker et al. (2011) applied a 43 % reduction. Here, k c is a variable parameter in the ensemble, distributed uniformly between 0 (leaf litter input to the soil is unchanged from the natural rate) and 1 (no leaf litter input to soils in cultivated regions). The rate of change of soil carbon density in cultivated regions is thus given by ENTSML does not distinguish between crops and pasture because at the level of detail of this parameterisation, and with a single plant functional type PFT, the principal difference between crop and pasture is the value of k c . This varies between crop and pasture but equally has substantial spatial variation between different cropland areas (and also pastures). In our parametric approach to the quantification and reduction of uncertainty , the most appropriate way to represent this variation is to consider a range of values for k c . Further discussion, including a description of what uncertainty in k c is likely to represent, can be found in Sect. 4.
The rate of change of soil carbon overall is the weighted average of the value in natural and cultivated regions (Eqs. 3 and 4): where f c is the grid cell fraction under cultivation (a prescribed forcing that is both temporally and spatially variable, see Sect. 3). This simplifies to When cultivated areas are expanded we assume that 50 % of the cleared carbon is released immediately to the atmosphere. The remaining 50 % is added to the soil carbon reservoir where it respires on timescales that are dependent on the local climate. The model of Stocker et al. (2011) releases 25 % directly to the atmosphere, with the remainder split evenly into two pools that decay on timescales of 2 and 20 yr. These authors concluded that the assessment of LUC emissions on decadal or longer timescales was not sensitive to these values. Although ENTS has only one soil carbon pool, we capture the uncertainty in the decay rate of cleared carbon through the ensemble parameter k32  that controls the temperature dependence of respiration. To illustrate the range of variability, at 15 • C, the approximate global average temperature, the ensemble spread in k32 produces respiration decay timescales of 11 to 25 yr.
Reforestation is modelled analogously to deforestation. Vegetative carbon is instantaneously returned to natural values, 50 % of the carbon coming from the soil and 50 % from the atmosphere. Effectively this balances an unphysical instantaneous recovery of vegetation carbon by an unphysical instantaneous loss of soil carbon, the effect of which is to allow a more realistic gradual recovery of the total terrestrial carbon pool, at a timescale dependent on soil respiration. This seems preferable to the instantaneous removal of 100 % natural vegetation carbon from the atmosphere. The approximation arises from the simplifying decision to only model potential vegetative carbon in ENTSML, so that timescales of reforestation from LUC-driven bare soil are not captured. It is worth noting that this is behaviour of ENTSML, not ENTS itself which does describe the evolution of potential vegetation from bare soil.
We model the climatic impacts of LUC through its effect on surface albedo and roughness length. (The soil moisture bucket implementation is unchanged but bucket capacity is affected by LUC due to the change in soil carbon.) We note that although soil carbon is reduced by the direct effect of LUC via the k c land management parameter, which is assumed to be positive, in general the climatic impacts of LUC act to oppose this change (reduced respiration rates due to albedo-driven cooling). This is discussed further in Sect. 4.
In cultivated regions that are not snow covered, albedo is assumed equal to α C , an ensemble parameter that is allowed to vary in the range 0.12 to 0.18 (parameter α c in Table 1). The selection of this range is based on values estimated for various crops (Hansen, 1993). The albedo in natural regions, α N , is calculated as a function of local potential vegetative and soil carbon density, using the standard ENTS expression . The albedo of the grid cell is a weighted average of the two: Snow-covered albedo and local climate are strongly dependent on the nature of the local vegetation (Betts, 2000). Rather than introducing an additional parameter, we capture the increased albedo of snow-covered deforested regions with the ENTS expression applied to the (reduced) average vegetative carbon density across the gridcell C v = (1 − f c ) C vN (see Appendix A, Eq. A5). Roughness length is derived similarly, by applying the average vegetative carbon density to the standard ENTS expression.
In summary, ENTSML is a minimal spatial model of the carbon cycle and surface physics feedbacks to land use change, with an appropriate level of complexity for the single PFT ENTS model on which it is based. The simple parameterisations and single PFT formulation are well suited to investigate the parametric uncertainty associated with global scale quantities. ENTSML is not a robust tool at subcontinental spatial scales. This reflects an appropriate level of complexity for GENIE, especially in the EMBM-atmosphere configuration applied here.

GENIE ensemble
A 671-member ensemble of transient simulations was performed, varying 28 key model parameters. In a previous calibration exercise , 24 "base" parameters of the unmodified GENIE model without LUC implementation were varied to produce an ensemble of 471 plausible preindustrial "spin-up" simulations . These 24 base parameters were selected for their importance to the carbon cycle, either directly or indirectly through their role on climate and ocean dynamics. Four additional parameters were added for the transient experiment described here (Table 1). Two of these parameters (the land management parameter k c and crop albedo α c ) represent uncertain processes in the new LUC implementation (Sect. 3). The third parameter, K LW1 describes the uncertain outgoing longwave radiation response to changes in radiative forcing (Holden et al., 2010). The range for K LW1 produced a climate sensitivity distribution of 2.4 to 5.1 • C (Holden et al., 2010). Finally, k 14 (Eq. 2) describes the uncertain response of photosynthesis to changing CO 2 concentrations.
Four ensemble experiments were performed, summarised in Table 1. Experiment 1 is a 100-member ensemble, randomly selected from the 471-member spin-up ensemble, so that for each ensemble member, all 24 base parameters were equal to those of the corresponding model spin-up. A maximin latin hypercube design was used to set up the values for the additional four parameters. This ensemble applied a conservatively narrow range for α c (0.14 to 0.16), a range that produced a globally averaged LUC radiative forcing from −0.2 to −0.5 W m −2 in exploratory experiments. Experiment 2 was performed on the remaining 371 spin-ups, allowing a broader range for α c (0.12 to 0.18). Experiments 3 and 4, both using the same subset as Experiment 1, were performed in order to better calibrate the tails of the k 14 distribution, regions that exert the greatest leverage on the subsequent emulations and calibration (Sect. 4). In order to ensure a uniform prior, only Experiments 1 and 2 were applied to the calibrated outputs in Sect. 6.2. All four experiments were used as training data for the emulator (Sect. 5).
Each transient simulation continues from the relevant spinup simulation and runs from 1 AD to 2005 AD. Boundary conditions are held fixed for the first 850 yr to allow the equilibration of climate, vegetation and surface-ocean, correcting for any drift due to minor inconsistencies with the spin-up boundary conditions. The four additional parameters have a minimal effect on the preindustrial spin-up state. Parameters k 14 and K LW1 only affect simulations that are not in the preindustrial state. Crop albedo α c and the land management parameter k c do have a small effect on the spin-ups due to LUC at 850 AD (the spin-ups applied no LUC forcing). However, the ensemble-averaged CO 2 drift of −3 ppm over the 850-yr spin-on (1 to 850 AD), suggests that any residual drift during the 150-yr calibration interval (1850 to 2000 AD) is not significant.
Simulations were forced with CMIP5 fossil fuel CO 2 emissions (including cement and gas flaring) from 1751 to 2005 AD, http://cmip-pcmdi.llnl.gov/cmip5/forcing.html# CO2-emissions. Other forcings (850 to 2005 AD) -land use change, solar variability, orbital configuration, non-CO 2 trace gases, direct aerosol effects and volcanic forcing -are described in Eby et al. (2012). In the GENIE implementation, a globally uniform perturbation to the radiative forcing is added to reflect non-CO 2 trace gases, the direct aerosol effect and volcanic forcing.
The land use change dataset applied was PMIP3 data (800 to 1699AD, Pongratz et al., 2007 linearly blended over the period 1500 to 1699 AD with the CMIP5 historical RCP land use data (1500 to 2005 AD, Hurtt et al., 2011). This combined dataset (Eby et al., 2012) provides fractional coverage of crops and pasture at annual resolution and at 0.5 • degree spatial resolution. The crop and pasture data were summed together and integrated onto the ∼10 • GENIE grid (retaining the annual resolution).

LPJmL simulations
To supplement the GENIE analysis, we performed three simulations with the global crop and vegetation model LPJmL (Bondeau et al., 2007). These historical transient simulations ran from 1900 to 2002 AD. The experiments were performed to isolate the responses of LPJmL to LUC, climate and CO 2 fertilisation. The three experiments are 1. LUC, climate change and CO 2 forcing (LPJmL LCC); 2. LUC and climate change forcing (LPJmL LC); and 3. LUC forcing (LPJmL L).
The historical climate forcing (CRU TS2.1, Mitchell and Jones, 2005) and LUC forcing are described in Fader et al. (2010). The LUC dataset prescribes the historic and present geographical distribution of 12 individual crop functional types (either irrigated or rainfed) as an annual fraction of the grid cell. The crops can coexist in a grid cell with the nine competing natural PFTs. Carbon fluxes, (gross primary production, autotrophic and heterotrophic respiration), carbon pools (represented as leaves, sapwood, heartwood, storage organs, roots for the plant biomass and litter and soil carbon) and water fluxes are modelled, accounting 33 explicitly for the dynamics of natural and agricultural vegetation. Carbon and water fluxes are linked to vegetation patterns and dynamics through the linkage of transpiration and photosynthesis and thus take plant water stress into account. Rising atmospheric CO 2 concentration affects transpiration and biomass production through physiological and structural plant responses (Sitch et al., 2003). Above ground biomass of the cropland area (present-day 1495 Mha) was averaged over the growing season, with sowing dates calculated after Waha et al. (2012). Biomass on the set aside area, represented by grass, grows between the cropping season. The time-weighted mean of crop and setaside biomass gives the average above ground biomass estimation of the cropland area. Managed grassland (presentday 2716 Mha) is grown over the entire year and is harvested whenever the increment of net primary production (NPP) exceeds 300 gC m −2 .

Validation of ENTSML
The LUC implementation is validated through the application of a 20-member subset of Experiment 1, described in Joos et al. (2012). This ensemble was filtered from Experiment 1, applying the constraint of plausible atmospheric CO 2 . The subset was used in the interests of computational efficiency, though we note that this approach does not allow us to consider the full range of modelled uncertainty. For the purposes of this validation, we apply the ensemble to a simulation with only direct LUC forcing. In order to isolate the direct effects of LUC from the indirect effects of LUC emissions on climate and photosynthesis, atmospheric CO 2 is relaxed to 280 ppm. Aspects of this experiment have been previously discussed in the inter-model comparison of Eby et al. (2012). Figure 1a and b illustrate the ensemble-averaged spatial distributions of preindustrial (850 AD) vegetation and soil carbon density. These are provided in part to demonstrate that the ensemble distributions are similar to those obtained in tuned simulations  and in part because errors in carbon released by LUC can be related to errors in the simulation of potential vegetation and soil carbon. Deserts are less distinct than observed because atmospheric moisture transport is too diffusive. Boreal forest is generally too sparse because of insufficient moisture transport into the continental interior, and too northerly in location. Some apparent weaknesses of the vegetation carbon distribution in Fig. 1a (excessive forestation in the tropics and in Europe) are reconciled with modern observations through the effects of LUC (not shown). We note, however, that the excessive 850 AD forestation in Central USA is likely a result of the over-diffusive atmosphere rather than the absence of LUC at this time. Figure 1c illustrates the change in fractional LUC cover from 1900 to 2000 AD. The comparison between Fig. 1a, b and c reveals the regions from where we expect dominant LUC emissions (both high potential carbon and large LUC change). Figure 1d illustrates the temporally averaged land-atmosphere LUC emissions from 1900 to 2000 AD. The largest fluxes are driven by removal of vegetation carbon. However, soil emissions can be significant, especially in the Himalayas. This results from high soil densities at altitude due to lapse rate cooling, and likely too high because soil weathering is not modelled. The regional distributions of the LUC fluxes are generally consistent with Houghton (2008). A notable exception is the large modelled flux from central USA where, as previously noted, potential vegetation density is overstated. We do not provide a comparison with the LUC flux distributions of Raddatz (2010) as the two approaches apply quite different assumptions. Dominant emissions in GENIE arise from regions associated with densely vegetated, deforested regions. Dominant emissions in Raddatz (2010) are associated with high population density.
In summary, although ENTS is based around a single PFT, the climatic dependencies of this PFT are sufficient to provide a reasonable spatial description of vegetative carbon density. As a result, simulated spatial patterns of LUC emissions also vary reasonably depending on whether the local potential vegetation is "forest" (high vegetation carbon density) or "grassland" (low vegetation carbon density). Figure 2a compares the ensemble-averaged temporal history of LUC emissions with estimates of Houghton (2008) and with the LPJmL L simulation. GENIE underestimates the LUC flux, especially in recent decades. We note that this underestimate of LUC emissions is a common feature of EMICs (Eby et al., 2012). It is, however, also worth noting that neither the analyses of Houghton (2008) nor the LPJmL L simulation include feedbacks due to the direct climatic impact of LUC. These climatic feedbacks (which are represented in the coupled EMICS), most notably the uptake of soil carbon driven by LUC albedo cooling, may explain part of the simulated differences. However, structural issues, most notably errors in the distributions of potential vegetation, are likely to be at least as significant.

The land management parameter k c
The ensemble average in Fig. 2a conceals a substantial variability in the magnitude of the simulated LUC flux. The land management parameter k c , the fractional reduction in leaf litter input to soils under LUC, is the dominant driver of LUC emission uncertainty. This is illustrated in Fig. 2b, which plots the time integrated LUC emissions (1850 to 2000) as a function of k c . We note that very high values of k c appear to have been ruled out by the requirement for plausible modern CO 2 that was used to filter this 20-member ensemble; the unfiltered ensemble (Experiment 1) evenly samples k c in the range 0 to 1. We apply the constraint of total LUC emissions (1850 to 2000 AD) to derive a simple calibration for k c . Houghton (2008) emissions are 149 GTC. LPJmL emissions are 181 GTC. The linear fit in Fig. 2b suggests that appropriate values of k c are 0.45 and 0.54, respectively. We apply these values as centres of two alternative priors in the Bayesian calibration of k 14 that follows. The 1-sigma width of the priors, assumed to be Gaussian, is taken as the RMSE of the simulations with respect to the linear fit (0.11). It is useful to discuss what uncertainty in k c is likely to represent. By design, this parameter was intended to represent the effects of land management on soil carbon. Soil carbon is reduced under cropped land because harvesting reduces the carbon returned to the soil and because tillage increases soil oxidation rates. However, by applying a prior to k c that favours reasonable LUC emissions, we subsume structural deficiencies of ENTSML into the parameterisation, such as the excessive potential vegetation in central USA or the excessive potential soil carbon at high altitude. Notably, the ENTS soil model has a single layer that reacts instantaneously throughout its depth to surface temperature changes. Deforestation results in radiative cooling due to increased albedo, although, as has been demonstrated in more complex models, this is countered by decreased evaporative cooling, especially in the tropics (Arora and Montenegro, 2011). It is likely that an excessive decrease in soil respiration (and increase in soil carbon density) would result from the deforestation-driven cooling, i.e. the relatively high values of k c needed to produce reasonable LUC emissions are likely due, at least in part, to this sensitive soil carbon response to LUC (both crop and pasture). It is important to note that our approach to calibrating k c has the effect of negating this LUC-climate-soil carbon sink. The parameter k c is shifted by our choice of prior to values that result in LUC emissions that are consistent with Houghton (2008), thus counteracting direct LUC-climate-carbon feedbacks. Although the presence of the albedo-driven sink cannot be ruled out in nature, we prefer to assume that the sink is negligible given the single soil layer in ENTS (i.e. by applying this prior). Sensitivity to this assumption is tested by also considering a significantly lower k c prior (centred on 0.2) in the analysis of Sect. 6.1. This choice is influenced by previous modelling studies with tuned LUC models (Olofsson and Hickler, 2008;Stocker et al., 2011), which suggest a ∼40 % reduction of leaf-litter input to arable soils, approximately equivalent to a global value of k c ∼ 0.2 when applied to all cultivated regions (ENTSML does not distinguish between crops and pasture).

Methods
The flow chart (Fig. 3) summarises the methodology. The ensemble displayed a wide range of responses, with 671member ensemble-averaged terrestrial carbon loss over the period 1850 to 2000 AD of 128 ± 139 GTC (1σ uncertainties are provided throughout). This compares to BernCC model estimates of ∼100 to 120 GTC from the range of LUC scenarios considered by Stocker et al. (2011), neglecting their illustrative scenario of a linear scaling of fractional LUC from 10 000 BC to the present. The comparison between GENIE uncertainty (∼100 GTC, driven by uncertain model parameter values) and BernCC uncertainty (∼10 GTC, driven by alternative temporal evolutions of LUC) suggests that our neglect of the latter is reasonable for these purposes. We note that 123 of the 671 completed GENIE simulations exhibited an increase in terrestrial carbon storage over this period. Although the direct consequences of LUC can only reduce terrestrial carbon in the model (through reductions in vegetation density and in the leaf litter-rate to soils), the effects of carbon emissions (CO 2 fertilisation and climate) and indirect (climate-driven) LUC feedbacks, acting on both natural and cultivated vegetation, can combine to increase the modelled global terrestrial carbon reservoir.
In order to calibrate the model response, we apply the constraint of CO 2 , the simulated increase in atmospheric CO 2 from 1850 to 2000 AD. This constraint is ideal for vegetation calibration as it is associated with a negligible observational error, but exhibits a wide range of simulated values across the ensemble (117 ± 41 ppm, cf. the observed increase of 84 ppm). It is very important to note that the ensemble variability is dominated by two of the poorly constrained terrestrial vegetation and land-use change parameters: 41 % of the variance in CO 2 is captured by a linear dependence on the CO 2 fertilisation parameter k 14 and 28 % by a linear 35 dependence on the land management parameter k c . The ensemble variability due to all other 26 parameters together generate an uncertainty of ∼16 ppm, compared to k 14 and k c , which together contribute ∼38 ppm. For this reason, we only consider the calibration of these two parameters. To further illustrate, the dominant correlations between parameters and CO 2 are CO 2 fertilisation k 14 : −0.64; crop management parameter k c : 0.53; ocean wind-stress scaling WSF: −0.16; soil respiration temperature dependence SRT: 0.15; and crop albedo α c : −0.15. WSF is the dominant parameter driving uncertainty of ocean uptake . The weak correlation with SRT may be partly explained through the two opposing drivers of changing soil respiration -global CO 2 -driven warming and local LUC-driven cooling. Given the dominant controls exerted by k 14 and k c , the remaining 26 parameters are not well constrained by CO 2 , although we note that, in general, they have already been constrained by the requirement for preindustrial plausibility. It should also be emphasised, however, that these 26 parameters do contribute to significant uncertainty across the ensemble that, together with inherent structural error, limit the degree to which CO 2 fertilisation can be constrained.
We choose to constrain on the basis of CO 2 rather than present-day CO 2 as we are here concerned with the response of the system to changing CO 2 and this approach avoids additional uncertainty due to an imperfect prediction of the preindustrial state. We prefer to avoid the alternative approach of tightly constraining by preindustrial CO 2 as this would contradict our ensemble design philosophy which is to rule out only those simulations which are uncontroversially implausible, necessary  for the Bayesian calibration which follows. We note that ensemble-averaged preindustrial CO 2 (280 ± 15 ppm) is well centred on observations. The variability can be largely attributed to ∼centennial-scale instabilities that are apparent in ∼30 % of the simulations that appear to be related to the presence of alternative stable states, likely driven by ocean convection and/or stratification-dependent mixing. We choose not to filter out these simulations, in part as it was found to be very difficult to distinguish the instability from natural variability with an objective test. Imposing (somewhat subjective) tests to eliminate simulations that display an instability was not found to have a significant effect on the emulation and calibration that follows, but including all simulations provides a weaker but more reliable constraint as it captures the complete range of ensemble variability.
It is not appropriate to derive marginal probability distributions for the two parameters without first considering their joint probability distribution, especially given the strong dependencies of CO 2 on both of these parameters which leads to strong correlations between the parameters in plausible parameter space. The 671 completed simulations are not sufficient to sample the 28-dimensional input space to properly quantify either the expectation or the parametric uncertainty throughout the 2D subspace of CO 2 fertilisation and land management. We note that the parametric uncertainty of CO 2 within the 2-D subspace (i.e. arising from the other 26 parameters) is not constant, ranging from 14 ppm (low k c , high k 14 ) to 24 ppm (high k c , low k 14 ). In order to generate sufficient data to adequately sample the 28-dimensional input space, we first build an emulator of CO 2 , building on Edwards et al. (2011) who used an emulator to identify implausible regions of input parameter space. This emulator is needed to fully and evenly sample the 28-dimensional input space in order to provide well-quantified estimates of both expectation and uncertainty at all points in the k 14 /k c subspace.
We take the 28 parameters as emulator inputs and construct a cubic emulator following a similar procedure to that described in Holden et al. (2010). Our emulator attempts to fit the simulated values of CO 2 to a polynomial function of the input parameters. A quadratic emulator was first built, allowing cross terms between all 28 parameters, to which we allow the addition of cubic terms before applying the Bayes information criterion (BIC) to reduce the model size. The additional cubic terms considered were the four cubic terms that were generated by an investigative emulation that allowed all possible cubic terms but considered only the 16 significant parameters in the quadratic model. This approach was taken to improve the efficiency of the process relative to an approach that considers all possible cubic terms from 28 model parameters. Two cubic terms were retained after application of BIC (k 3 14 and κ T × k 2 14 , where κ T is atmospheric heat diffusivity). The resulting emulator fits the simulated data well, with an R 2 of 94 %. We then designed a 14 100-member parameter set to apply as input to this emulator. The design reproduces 470 of the 471 simulated parameter sets (omitting a single parameter set which did not complete in Experiment 2) thirty times, replacing the four transient parameters with a 4 × 14 100 matrix of randomly generated inputs across the ranges in Table 1. The output of this emulated ensemble is compared with the simulated ensemble in Fig. 4, illustrated by the marginal dependencies on k 14 and k c . The response to changes in these parameters is well captured by the emulator. Although the most extreme values of CO 2 are underestimated by the smooth polynomial functions, the emulated ensemble distribution of 115 ± 39 ppm compares favourably to the simulated distribution (117 ± 41 ppm), suggesting this reduction in variance is unlikely to be significant. It is worth noting that the simulator cannot reproduce observed CO 2 (84 ppm) with low values of k 14 1 . To derive a joint probability distribution we first generate a two-dimensional matrix of data bins. We subdivide the k c /k 14 parameter space into 10 × 10 bins, linearly spaced for k c (i = 1, 10) and quadratically spaced for k 14 (j = 1, 10), across the ranges of their prior distributions (Table 1). Within each of the 10 × 10 bins we calculate the mean µ ij and standard deviation σ ij of simulated CO 2 . Within each bin the variance is dominated by uncertainty due to the remaining 26 ensemble parameters as, to first order, k 14 and k c are fixed. We apply Bayes' theorem to the bin statistics to derive a posterior distribution for each bin. Analogously to Rougier (2007), we apply p θ ij | CO 2 = cϕ CO 2 , µ ij + µ ε , σ 2 ij + σ 2 ε p θ ij , (8) where ϕ describes a normal distribution, evaluated at CO 2 = 84 ppm with mean µ ij + µ ε and standard deviation √ (σ 2 ij + σ 2 ε ), µ ε and σ ε are model structural bias and structural error, respectively, c is a normalising constant, θ ij are the binned values of the k 14 /k c parameter pair, and p(θ ij ) the prior probability we assign to that combination. In this application, the resulting pdf was found to give similar results to a more rigorous calibration procedure that explicitly integrated over the individual outputs of the emulated ensemble (i.e. rather than approximating emulator output with a binned normal distribution).
It is essential to include estimates of structural error and bias in a model calibration to avoid over-constraining the resulting pdf. We derive an estimate for structural error by applying an approach influenced by Murphy et al. (2007) who suggest the inter-model spread of a quantity can be used to provide a measure of structural error, reflecting, at least in part, different structural choices that can be made. Although such an approach is likely to provide an underestimate because different models likely share some sources of structural error, it is conservative in the sense that a component of the inter-model variability will arise from parametric uncertainty (but be attributed to structural error). C 4 MIP simulations forced by SRES A2 exhibit a standard deviation of ∼90 ppm in the simulated increase in CO 2 in 2100 relative to 2000 (Friedlingstein et al., 2006). We note that Friedlingstein et al. (2006) found no systematic differences between the behaviours of OAGCMs and EMICs. We obtain an estimate of structural error by linearly scaling the C 4 MIP intermodel variability by the relative increase in CO 2 , i.e. by CO 2 /480 ppm, where 480 ppm is the C 4 MIP intermodel average increase (2000 to 2100). This gives a CO 2 structural error of ±17 ppm. The ensemble relationship between CO 2 and total emissions T (the sum of time integrated fossil fuel and net terrestrial carbon emissions) is well described (R 2 = 91 %) by This relationship suggests that the assumed structural error in CO 2 is equivalent to a 2-sigma uncertainty in accumulated historical emissions of ∼120 GTC, sufficient to encompass uncertainties in historical fossil fuel emissions (cf. accumulated emissions of ∼280 GTC).
In our LUC implementation, a potentially significant source of structural bias is the assumption that vegetative carbon is negligible in cultivated regions, leading to excess simulated LUC emissions. The LPJmL LCC simulation quantifies this neglect, yielding a global above-ground LUC biomass of 3 GTC, equivalent to an over-estimation of CO 2 by ∼1 ppm which can reasonably be neglected in the calibration.

Calibration of CO 2 fertilisation parameter
We apply Eq. (8) to derive the joint probability distribution for k c (land management) and k 14 (CO 2 fertilisation), illustrated in Fig. 5. In these plots, k 14 is re-expressed as the percentage increase in photosynthesis in response to a doubling of CO 2 from preindustrial levels in order to facilitate comparison with alternative estimates. Figure 5a applies a uniform prior; Fig. 5b constrains the loss in cultivated soil carbon to more reasonable values by applying the k c prior centred on 0.45 (Sect. 4). In both plots, structural error of σ ε = 17 ppm and structural bias µ ε = 0 ppm are applied (Sect. 5). We integrate over the joint distribution in Fig. 5b to derive the marginal probability distribution for k 14 , plotted as the blue curve in Fig. 6a and b. This analysis demonstrates a GPP increase in response to doubled CO 2 that is most likely 40-60 %, and very likely (90 % confidence) to exceed 20 %. Note that high values of k 14 are not well constrained by the pdf. This is in part due to the low sensitivity of the saturating Michaelis-Menten function at high values of k 14 , although the analysis would clearly have benefited from an ensemble that spanned higher values. Notwithstanding this, the pdf is fertilisation parameter k 14 , applying Eq. 8 with a) a uniform prior distribution p(θ ij ) = 1 and b) 2 a prior assumption for k c , p(θ ij ) = φ(k c , 0.45, 0.11). Note that k 14 is re-expressed as the 3 percentage increase in GPP in response to a doubling of CO 2 from preindustrial levels.  1 and (b) a prior assumption for k c , p(θ ij ) = ϕ(k c , 0.45, 0.11). Note that k 14 is re-expressed as the percentage increase in GPP in response to a doubling of CO 2 from preindustrial levels. sufficient to quantify a lower bound for k 14 (90 % confidence) and to identify the approximate range of most likely values. Figure 6a illustrates sensitivity to the assumed prior distribution for land management k c . The base case k c prior (mean 0.45, standard deviation 0.11) constrains CO 2 fertilisation with temporally integrated LUC emissions of 149 GTC (Houghton, 2008) as described in Sect. 4.1. The sensitivity experiment with k c centred on 0.54 was constrained with the higher LUC emissions of the LPJmL simulation (181 GTC), and requires a slightly stronger CO 2 fertilisation effect to balance the carbon budget. i.e. to match observed CO 2 . The sensitivity experiment with k c centred on 0.2 is equivalent to a constraint based on LUC emissions of 55 GTC. This assumption is equivalent to an LUC-climate feedback sink (the effect of LUC-driven albedo changes on soil respiration) of ∼90 GTC. The additional sink implicit in this analysis shifts the most likely CO 2 fertilisation effect (in response to doubled CO 2 ) to ∼30 %. A further sensitivity experiment applies a uniform prior and produces a similar, though somewhat broader, pdf to the base case. This similarity does not suggest that the calibration is insensitive to the prior, but rather it reflects the fact that the base case prior is approximately centred on the ensemble distribution of k c . Figure 6b illustrates sensitivity to the structural error assumption. Two analyses are performed to investigate the structural bias term. The first of these has a k c prior centred on 0.45 and a structural bias of µ ε = 3.2 ppm. This bias is included for the underestimate of LUC emissions (converted to CO 2 with Eq. 9) compared to Houghton (2008) when this prior is applied to the 20-member ensemble (Sect. 4.1). The second sensitivity analysis has a k c prior centred on 0.2 and a structural bias of µ ε = 22.6 ppm. This is an alternative approach to the calibration: we assume that the k c = 0.2 prior is reasonable, but instead correct for the low LUC emissions relative to Houghton (2008) through the structural bias term. These two analyses produce similar pdfs. A further illustrative analysis applies no structural error (σ ε = 0 ppm). This calculation imposes a stronger constraint on k 14 , but one which cannot be justified as it is equivalent to assuming that the only model uncertainty is due to uncertainty in the input parameters.
The base case pdf for k 14 is generally robust with respect to these various sensitivities, with the exception of the choice of k c prior. The dependency on k c is unsurprising given that both parameters are comparably important in determining CO 2 , necessitating the analysis of the joint probability distribution. As previously noted, the calibration does not quantify the upper bound for CO 2 fertilisation. The most likely GPP response to a doubling of CO 2 is 40-60 %, equivalent to a k 14 range of ∼300-600 ppm (and consistent with the simulated ranges most likely to produce plausible CO 2 plotted in Fig. 4a). We note that these values for k 14 are significantly greater than in the standard ENTS value (Williamson et al., 2006, 145 ppm equivalent to a 24 % increase in GPP under doubled CO 2 ).

Calibrated ensemble output
We apply the parameter calibrations to weigh individual simulated ensemble members in order to derive calibrated model outputs of global GPP (2000 AD), post-industrial terrestrial carbon change (1850 to 2000), land-atmosphere carbon flux (1990-2000 average) and ocean-atmosphere carbon flux (1990-2000 average). Each of these quantities is of importance to the global carbon cycle and each is associated with significant uncertainty. These analyses were performed in part to provide independent estimates of these widely studied metrics, and in part as a validation of the calibration.
We consider the 470 completed simulations that comprise Experiments 1 and 2. We do not include Experiments 3 and 4, as these simulations were performed to better quantify the simulated response at the tails of the k 14 distribution. Their inclusion would compromise the prior distribution of k 14 by attributing too much weight to the tails.
The temporal evolution of atmospheric CO 2 is plotted in Fig. 7. The ensemble-averaged distribution overstates the rise in CO 2 , unsurprising given that the calibration of Sect. 6.1 suggests that the ensemble-averaged value of k 14 (212 ppm, equivalent to a 32 % GPP increase under doubled CO 2 ) is too low to match the observed change in atmospheric CO 2 . The calibrated curve is better matched, although it underestimates atmospheric CO 2 by up to ∼10 ppm between 1850 and 1950 (and overestimates CO 2 ). The sharp drop in CO 2 at 1815 is associated with the Tambora volcanic eruption and reflects the unreasonably fast response time of the single-layer soil pool to surface cooling (see Sect. 4.1). However, the comparably rapid recovery suggests that this sink is unlikely to 39 explain the low CO 2 estimates that persist through to 1950. A more likely explanation is the underestimate of early LUC emissions (Fig. 2a).
The following discussion considers a series of posterior pdfs (Fig. 8). In order to generate these pdfs, we subdivided the output space into data bins and summed the probability weights of all simulations that fall into each bin. The width of the bins was chosen to remove fine structure from the pdf. We assume that this fine structure is unlikely to be physically meaningful. The binning approach is crudely equivalent to ascribing a structural error to each model output, rather than accepting it as a point estimate. Each figure illustrates the calibrated GENIE pdf, together with the 66 % confidence interval, and the point estimates from the LPJmL LCC and LPJmL LC simulations, respectively, including and neglecting the direct effect of CO 2 on vegetation. Figure 8a plots the pdf of global natural GPP (2000 AD). ENTSML does not calculate GPP in cultivated regions but applies the simplifying assumption of instantaneous harvesting and grazing. The GENIE calibration (likely in the range 107 to 152 GTC yr −1 ) is therefore limited to global GPP in regions unaffected by LUC. The LPJmL LCC simulation quantifies this neglect to be ∼7 GTC yr −1 (crops) and ∼29 GTC yr −1 (pasture). The LPJmL LCC estimate of natural GPP of ∼97 GTC yr −1 (1991 to 2000 average) is somewhat lower than the GENIE pdf, although consistent with the 2-sigma uncertainty. A recent observational based estimate of global GPP (Beer et al., 2010) is 123 ± 8 GTC yr −1 , of which 15 GTC yr −1 was attributed to croplands. It is important to note that the GENIE calibration of GPP is not well constrained by CO 2 (which rather constrains the postindustrial change in GPP). Instead, GPP is primarily (but only weakly) constrained by the requirements for plausible vegetation that were applied to the spin-up design. These produced an ensemble-averaged preindustrial global vegetative carbon of 492 ± 94 GTC . The dominant drivers of GENIE uncertainty in global GPP are atmospheric atmosphere flux (1990 to 2000) and d) annually-averaged ocean-to-atmosphere flux (1990 to 7 2000). Vertical dashed lines illustrate the 66% confidence interval. In a-c, the dark green 8 vertical bar is the LPJmL_LCC simulation (including the direct CO 2 effect) and the light 9 green vertical bar is the LPJmL_LC simulation (neglecting the direct CO 2 effect). moisture diffusivity (R 2 = 44 %), soil respiration temperature dependence (16 %), baseline photosynthesis rate (13 %), atmospheric heat diffusivity (10 %), baseline leaf litter rate (10 %) and, finally, the CO 2 fertilisation parameter k 14 (6 %).
Global GPP in LPJmL (natural and cultivated) increases by 9.3 % due to the direct effect of a 22 % increase in CO 2 (i.e. between 1900 and 2000 AD). This compares to the GE-NIE calibration (Sect. 6.1) in response to a doubling of CO 2 of a most likely increase of 40-60 %, very likely to be more than 20 %. Figure 8b plots the posterior distribution of post-industrial terrestrial carbon change . This analysis predicts a 66 % probability that the change in terrestrial biomass lies between a loss of 0 and 128 GTC. The probabilityweighted mean of 66 GTC provides a best estimate of net post-industrial land emissions, consistent with data-driven estimates (for the period 1800 to 1994) of 39 ± 28 GTC (Sabine et al., 2004). We note that the linear fit of Eq. (9) estimates most likely historical net land emissions to be 13 GTC. The higher terrestrial carbon losses in the probabilistic calibration likely reflect the ensemble prior distribution for k 14 , which favours a weaker land sink, being logarithmically spaced from 30 to 700 ppm (cf. the calibrated estimate of ∼300 to 600 ppm). The LPJmL LCC result (including CO 2 fertilisation) lies within the calibrated pdf. The LPJmL LC result (neglecting CO 2 fertilisation) does not. This suggests that plausible post-industrial net land emissions can only be simulated in LPJmL (at default parameters) if the direct effect of CO 2 is included. Figure 8c illustrates calibrated pdf of the land-atmosphere carbon flux over the recent period (1990 to 2000 average). This has a probability-weighted mean of −0.7 GTC yr −1 and a 66 % confidence interval of uncertain sign, ranging from −1.7 to 0.4 GTC yr −1 . LPJmL can only simulate a plausible land-atmosphere flux under the assumption of a direct CO 2 effect. The GENIE estimates compare with IPCC estimates (Denman et al., 2007) of −1.0 ± 0.6 GTC yr −1 . More recently, Le Quéré et al. (2009) estimated LUC emissions (1990 of 1.5 ± 0.7 GTC yr −1 offset by a terrestrial sink (1990 to 2000) of 2.6 ± 0.7 GTC yr −1 , implying a net land-atmosphere flux of −1.1 ± 1.0 GTC yr −1 . Joos et al. (1999) applied the single deconvolution approach, closely related to the approach taken here as both use an ocean model to constrain carbon uptake and solve for the residual land-atmosphere flux, to estimate the 1990 flux to be −0.8 GTC yr −1 . Figure 8d illustrates the calibrated ocean-atmosphere flux (1990 to 2000 average). This has a probability weighted mean of −2.3 GTC yr −1 and a 66 % confidence interval in the range −1.8 to −2.7 GTC yr −1 . These figures compare with 352 P. B. Holden et al.: A model-based constraint on CO 2 fertilisation estimates of −2.2 ± 0.4 GTC yr −1 (Denman et al., 2007, Le Quéré et al., 2009). The single deconvolution analysis of Joos et al. (1999) estimated the 1990 ocean-atmosphere flux to be −2.1 GTC yr −1 .

Summary and conclusions
We have incorporated a relatively simple implementation of LUC into ENTS, the dynamic terrestrial carbon module of GENIE, and applied the model in an attempt to constrain the response of terrestrial vegetation to increased atmospheric CO 2 . The advantage of our procedure is that it provides an estimate of CO 2 fertilisation strength that is independent of plant and ecosystem manipulation experiments. Combining all sources of information using Bayesian methods would lead to a better estimate of the true value.
The resulting probability distribution for the GPP increase in response to doubled CO 2 suggests a significant effect, very likely to exceed 20 % (90 % confidence), with a most likely value of 40-60 %. LPJmL simulations (9.3 % increase in GPP in response to a 22 % increase in atmospheric CO 2 ) are not inconsistent with the GENIE estimate. Although a probabilistic treatment of LPJmL is beyond the scope of this work, the LPJmL simulated land-atmosphere fluxes fall within the GENIE pdfs, an essential requirement to validate the calibration as a correct treatment of structural error should ensure that this is the case. Importantly, it is only the case when LPJmL accounts for the direct effect of CO 2 .
It is important to note GENIE calibration applies to the increase from pre-industrial to present. The projected increase for CO 2 concentrations greater than present-day levels cannot be directly constrained by observations, only by the expectation that the functional form (Eq. 2) remains approximately valid at higher concentrations. These calibrated ranges suggest that the standard parameterisation in GENIE  underestimates the magnitude of the land sink. This behaviour may not be unique to GENIE. An intercomparison of intermediate complexity models (Eby et al., 2012) concluded that all the models are either overestimating the climate-driven terrestrial carbon source or are underestimating the terrestrial fertilisation sink (either from CO 2 or deposition of nitrogen).
Although a negligible fertilisation effect appears difficult to reconcile with this analysis, it is important to reiterate that some of the possible mechanisms that are likely to contribute to the terrestrial sink are not modelled. These processes are therefore subsumed into our calibration of CO 2 fertilisation. This estimate should, however, be useful for other Earth system models as they, in general, will neglect similar sinks. It is also very important to reiterate that the subsumed processes represent a source of structural error and that caution is required when extrapolating these results into the future.
The probabilistic approach we applied (Rougier, 2007) is well suited to the refinement of the calibration through the incorporation of a bias term. Un-modelled terrestrial sinks could be included through this term to improve the calibration. Additional independent constraints, such as the isotopic composition of atmospheric CO 2 , could also be applied.
Although the uncertainty is substantial, the result is useful for two reasons. Firstly, the analysis is by definition a calibration of the fertilisation effect net of other limiting factors, most notably nitrogen limitation (and its alleviation through nitrogen deposition). The calibration simply provides a dynamical description of the globally averaged vegetation response that reproduces the observed historical change in atmospheric CO 2 . Secondly, the approach does not impose a direct prior constraint on the strength of CO 2 fertilisation. As such, it provides a useful independent constraint that can be used in conjunction with existing estimates from CO 2 enrichment experiments.
Although the approach cannot capture the regionally disparate responses that are apparent in more complex models  it nevertheless provides a useful global constraint. This may be especially relevant for application to complex vegetation models that are not coupled to global models of the carbon cycle. The estimate of postindustrial carbon loss, in particular the quantification of uncertainty (likely in the range 0 to 128 GTC) may provide a useful constraint on post-industrial terrestrial emissions in order to ensure that they are consistent with the relatively well-constrained estimates of fossil fuel emissions and ocean uptake.
The autotrophic respiration rate is given by R vN = k 24 k 25 e −k 20/ RT a C vN (A2) where k 24 = 0.225 year −1 is the base rate of autotrophic respiration, k 25 normalising constant (defined so that R vN = k 24 C v when T a = 298.15 K), R is the universal gas constant and k 20 is varied across the ensemble ("VRA" in Holden et al., 2012). The leaf litter rate is given by where the base turnover rate k 26 is varied across the ensemble ("LLR" in Holden et al., 2012). The second expression is designed to represent self-shading (k 16 = 11.5 kgC m −2 ). When the land temperature T l ≥ 273.15 K, the heterotrophic respiration rate is given by where k 29 = 0.165 yr −1 is the base respiration rate, k 30 is a normalising constant, k 31 = 308.56 K, and k 32 is varied across the ensemble ("SRT" in Holden et al., 2012). A separate expression is applied when T l < 273.15 K to prevent unrealistic blow-up of the soil reservoir as T l approaches k 32 (see Williamson et al., 2006). Snow covered albedo is given by where α snow = 0.8 (the albedo of snow-covered bare soil) α snow v = 0.3 (the albedo of snow-covered vegetation) and, in ENTSML, C v = (1 − f c ) C vN . The parameter k 7 = 0.461 m 2 kgC −1 is fixed across the ensemble, with uncertainties in snow-covered albedo captured through the parameters that drive the uncertainty in potential vegetation.
We note that parameter values above for k 24 and k 29 differ from the values applied in Williamson et al. (2006). The revised values are the plausibility-filtered ensemble-average parameter values of Holden et al. (2010), in order to facilitate the centring of the ensemble distribution of global terrestrial carbon pools on plausible ranges.