Quantifying the influence of the terrestrial biosphere on glacial–interglacial climate dynamics

. The terrestrial biosphere is thought to be a key component in the climatic variability seen in the palaeo-record. It has a direct impact on surface temperature through changes in surface albedo and evapotranspiration (so-called biogeophysical effects) and, in addition, has an important indirect effect through changes in vegetation and soil carbon storage (biogeochemical effects) and hence modulates the concentrations of greenhouse gases in the atmosphere. The biogeochemical and biogeophysical effects generally have opposite signs, meaning that the terrestrial biosphere could potentially have played only a very minor role in the dynamics of the glacial–interglacial cycles of the late Quaternary. Here we use a fully coupled dynamic atmosphere– ocean–vegetation general circulation model (GCM) to generate a set of 62 equilibrium simulations spanning the last 120 kyr. The analysis of these simulations elucidates the relative importance of the biogeophysical versus biogeochemical terrestrial biosphere interactions with climate. We ﬁnd that the biogeophysical effects of vegetation account for up to an additional − 0.91 ◦ C global mean cooling, with regional cooling as large as − 5 ◦ C, but with considerable variability across the glacial–interglacial cycle. By comparison, while opposite in sign, our model estimates of the biogeochemical impacts are substantially smaller in magnitude. Ofﬂine simulations show a maximum of + 0.33 ◦ C warming due to an increase of 25 ppm above our (pre-industrial) baseline atmospheric CO 2 mixing ratio. In contrast to shorter (century) timescale projections of future terrestrial biosphere response where direct and indirect responses may at times cancel out, we ﬁnd that the biogeophysical effects consistently and strongly dominate the biogeochemical effect over the inter-glacial cycle. On average across the period, the terrestrial biosphere has a − 0.26 ◦ C effect on temperature, with − 0.58 ◦ C at the Last Glacial Maximum. Depending on assumptions made about the destination of terrestrial carbon under ice sheets and where sea level has changed, the average terrestrial biosphere contribution over the last 120 kyr could be as much as − 50 ◦ C and − 0.83 ◦ C at the Last Glacial Maximum.


Introduction
Terrestrial vegetation interacts with the climate in complex ways, both responding to and impacting climate conditions and hence creating an important feedback in the Earth system (e.g. Claussen, 2009;Davies-Barnard et al., 2014b;Harrison and Prentice, 2003;Jahn et al., 2005;Matthews et al., 2003;Pongratz et al., 2010). The influence of the terrestrial biosphere on climate occurs in two distinct ways. Firstly, there are a number of biogeophysical mechanisms such as changes in albedo or evapotranspiration that provide a direct physical influence on surface climate via changes in net solar radiation transfer, infrared loss, roughness length, latent heat loss, and, less directly, via changes in moisture exchange and hence transport. Climate feedbacks driven by these changes in terrestrial vegetation have been hypothesised to be partially responsible for some of the major past climate states (e.g. Bradshaw et al., 2015;Claussen et al., 2006;Crucifix and Loutre, Published by Copernicus Publications on behalf of the European Geosciences Union. T. Davies-Barnard et al.: Terrestrial biosphere carbon emissions 2002;de Noblet et al., 1996;Zhou et al., 2012), with many studies particularly focussing on the biogeophysical effects at the Last Glacial Maximum (LGM) (e.g. Hopcroft and Valdes, 2014;Jahn et al., 2005;Kageyama et al., 2012;O'ishi and Abe-Ouchi, 2013). The second way in which the terrestrial biosphere can influence climate is via variations in the carbon stored in vegetation and soil. This is a crucial component for understanding changes in the carbon cycle through the last glacial-interglacial cycle (Montenegro et al., 2006), and numerous attempts have been made to estimate the total carbon storage using a range of methods, such as inferences from marine and terrestrial carbon isotopes (e.g. Shackleton et al., 1977;Bird et al., 1994), databases of pollen (e.g. Adams and Faure, 1998;Crowley, 1995), and simple and complex modelling (e.g. Prentice et al., 1993;Kaplan et al., 2002;Köhler and Fischer, 2004;Brovkin et al., 2012;O'ishi and Abe-Ouchi, 2013). The resulting range of carbon storage change estimates is from a few hundred to about 1000 PgC (Ciais et al., 2012). One could add to this changes in the weathering of soil minerals and hence CO 2 uptake from the atmosphere, as well as nutrient, particularly phosphate, supply to the ocean and hence changes in the ocean productivity. For simplicity, we will not address these further here (except to include a basic silicate weathering feedback in our model analysis of the impacts of terrestrial carbon storage change).
Simulations of future vegetation changes show that the biogeochemical aspect can globally be around the same magnitude as the biogeophysical effects (e.g. Davies-Barnard et al., 2014b), meaning that there is uncertainty even in the sign of the net feedback with climate change. Both biogeophysical and biogeochemical effects likely also play an important role in past climate change and potentially the same fundamental uncertainty in the sign of the climate feedback might arise. However, model simulations have generally focussed on either the biogeophysical impacts of vegetation changes (e.g. Bradshaw et al., 2015;Claussen et al., 2006;Jahn et al., 2005;O'ishi and Abe-Ouchi, 2013;Shellito and Sloan, 2006) or biogeochemical impacts (e.g. Kaplan et al., 2002;Ciais et al., 2012) and the question of the overall feedback on climate rarely addressed, although Claussen (2009) argues that the net effect at the LGM is dominated by the biogeophysical effects.
One of the few examples where both have been combined and the net effects of vegetation on past climate have been estimated over long time periods is Brovkin et al. (2012). They used an Earth system model of intermediate complexity (EMIC) to suggest that the net effect of vegetation is to decrease global temperatures during the last glacial-interglacial cycle. But the model used is relatively coarse in resolution (10 • in latitude and 51 • in longitude) and reduced in physical process complexity. This may be important because of the local and spatially heterogeneous nature of biogeophysical effects and, depending on the location of the forest, the biogeophysical and biogeochemical effects of forest change can be very different (Bonan, 2008). For this reason, fully coupled general circulation models (GCMs) are commonly used in quantifying future climate changes to vegetated land surface (Brovkin et al., 2013a, b;Davies-Barnard et al., 2015;Davin and de Noblet-Ducoudré, 2010). The importance of considering both biogeophysical and carbon cycle impacts together at a finer scale when assessing the climate impacts of vegetation is illustrated by work quantifying the climate impacts of forest changes. Studies have found that deforestation would cause local high-latitude cooling (Betts, 2000), global warming (Davin and de Noblet-Ducoudré, 2010), or even slight global cooling (Davies-Barnard et al., 2014b, a). These outcomes are not predictable from looking at the biogeophysics or terrestrial biogeochemistry alone at coarse resolutions.
Here we present the first model analysis using a fully coupled dynamic atmosphere-ocean-vegetation GCM over the last 120 kyr that quantifies the net effect of vegetation on climate. (A prior study -Singarayer and Valdes, 2010 -did not have dynamic vegetation and hence could not directly evaluate the biogeophysical effects.) We separate the biogeophysical and biogeochemical effects of vegetation to understand the overall climate effect of vegetation over the last glacial cycle. We show that over the whole period the biogeophysical is the dominant effect, and that the biogeochemical impacts may have a lower possible range than typically estimated. We also highlight how the temporal scale affects the net impact of terrestrial biosphere changes.

Methods
We use the GCM HadCM3 to run a series of simulations with and without dynamic vegetation to provide the biogeophysical changes and the land carbon changes. To look at the climate impact of those vegetation carbon changes, we then use the GCM terrestrial carbon changes as an input to the EMIC cGENIE to calculate the resulting change in atmospheric CO 2 and global temperature.
For future climate changes studies, the response of atmospheric CO 2 concentrations (and hence climate) to changes in terrestrial carbon storage can be calculated using the transient response to cumulative emissions (TRCE) approach (Gillett et al., 2013), which demonstrated proportionality between carbon emissions and temperature rise (Goodwin et al., 2015). We include these estimates for completeness. However, this approach is only valid for relatively rapid changes. On the longer timescales of glacial-interglacial change, we need to take into account the full changes in ocean carbon chemistry and including the interactions of ocean and atmosphere with the solid Earth (e.g. weathering). To do this, we employ the "cGENIE" EMIC.  (Valdes et al., 2017). Though not from the latest generation of climate models, HadCM3 remains an extensively used model for many research applications around the world due to its computational efficiency, which means that long integrations and many ensemble members can be run.
HadCM3 is a three-dimensional, fully coupled, fully dynamic ocean, non-flux-adjusted global climate model (Collins et al., 2001). The atmosphere component, HadAM3, has a Cartesian grid with a horizontal resolution of 2.5 • × 3.75 • , 19 vertical levels, and a time step of 30 min (Pope et al., 2000). The ocean and sea-ice component has the same horizontal resolution as the atmosphere, with 20 vertical ocean levels.
The land surface scheme used for the atmosphere component of HadCM3 is the Met Office Surface Exchange Scheme, MOSES2.1 (Gregory et al., 1994;Cox et al., 1999). MOSES can also use an additional vegetation and terrestrial carbon model, TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics) (Cox, 2001;Cox et al., 1998). TRIFFID predicts the vegetation based on plant functional types using a competitive, hierarchical model. TRIFFID has two modes: equilibrium mode, which quickly brings the vegetation cover into equilibrium by running 50 years of TRIFFID for each 5 years of the climate model run, and dynamic, which runs TRIFFID every 10 days. TRIFFID and MOSES have nine land surface types, five of which are vegetation: broadleaf trees, needleleaf trees, shrubs, C 3 grasses, and C 4 grasses. These are known as plant functional types (PFTs) and have different leaf area index limits and other phenological differences in the model. Soil moisture in the model is represented on four layers of thicknesses (measured from the top) of 0.1, 0.25, 0.65, and 2.5 m (Essery et al., 2001).
The soil carbon is a single pool, increased by litterfall and decreased by respiration (Cox, 2001;Cox et al., 1998). The soil respiration is controlled by moisture and temperature and returns carbon dioxide to the atmosphere unless, as is the case here, the atmospheric carbon dioxide is fixed. The litterfall is an area-weight sum of the litterfall of the five PFTs in each grid cell (Cox, 2001;Cox et al., 1998). There is no permafrost component in the model, and soil in frozen regions is treated the same as in any other.
Assessment of the pre-industrial (PI) vegetation cover of HadCM3 by Valdes et al. (2017) shows good agreement with reconstructions of 1800 vegetation.

GCM simulations and experimental methodology
The simulations used here are revised versions of those described in Singarayer and Valdes (2010), who used HadCM3 version HadCM3-M1, which has an older surface scheme (MOSES1) than the MOSES2.1 used here and no dynamic vegetation. Two sets of 62 simulations were performed, covering the time period 120-0 ka: -The first set of 62 simulations used TRIFFID to predict vegetation changes. Each individual simulation was initialised from the previous MOSES1 simulations (which were run for 600 years) and were then run for a further 300 years with "equilibrium" TRIFFID and a final 300 years with fully dynamic vegetation. This set will be referred to as the Dynamic set.
-A second set of simulations uses static vegetation based on the pre-industrial simulation of the dynamic set (extrapolated to new land areas using a simple nearestneighbour algorithm). They are otherwise identical to the "Dynamic" set (see details below). These will be referred to as the "Static" set.
The differences between Dynamic and Static allow us to evaluate the biogeophysical and biogeochemical responses of terrestrial carbon cycle change.
Both sets of simulations are forced with the same changes in orbit, greenhouse gases (CO 2 , CH 4 , and N 2 O) and ice sheets, as in Singarayer and Valdes (2010), except that we use a revised ice sheet extent and elevation, as discussed in Singarayer et al. (2011).
We have also added a parameterisation of water transport from ocean to ice sheet in order to ensure that ocean salinity is conserved during each simulation. In the normal configuration of HadCM3, salinity is conserved by the numerical scheme but water that accumulates as snow on ice sheets is not interactively considered. A predefined (spatially varying) flux of water is prescribed into the model which minimises the salinity drift for the pre-industrial simulation, but this is not normally changed for other time periods. In our new parameterisation, we continue to add the predefined flux but also add an additional flux which is spatially uniform but temporally variable to ensure that the volume integral ocean salinity is relaxed back to its initial value, with a relaxation timescale of 10 years. This prevents any spurious long-term drifts in ocean salinity. Note that this model does not have a closed carbon cycle. There is no representation of carbon in the ocean and terrestrial carbon changes do not feedback to the atmosphere (since the greenhouse gas forcings are prescribed). However, the carbon that would have returned to the atmosphere can be inferred from the change in the carbon stores in the soil and vegetation, allowing the biogeochemical impact of vegetation to be understood, as well as the biogeophysical.  final climatologies are based on the last 30 years of each simulation.

EMIC description
The cGENIE Earth system model is used to calculate the impacts on atmospheric CO 2 over the glacial cycle and hence make a time-varying estimate of the contribution of biogeochemical changes to glacial-interglacial climate change. The model is based around a fast energy-balance-based atmosphere model coupled to a 3-D ocean circulation component and dynamic-thermodynamic sea ice (Edwards and Marsh, 2005), plus representations of ocean-atmosphere , ocean-sediment , and atmosphere-land (terrestrial weathering) (Colbourn et al., 2013) carbon cycling. As employed here, the non-seasonally forced ocean has eight levels and the configuration and selection of model parameterisations and parameter values is identical to that described in Lord et al. (2016). These choices are made to minimise experiment runtime and provide maximum traceability (to a previously used and indepth analysed configuration), respectively.

cGENIE carbon cycle simulations
The evolution of terrestrial carbon storage simulated by HadCM3 from 120 ka to the pre-industrial era was used to derive a forcing for cGENIE. In this, we created a continuous time series of the carbon flux from the terrestrial biosphere by calculating the difference in carbon storage calculated at the end of each HadCM3 time slice and then assuming a linear interpolation between these points. For the "Full" simulations, cGENIE was then run for 120 kyr using this forcing and starting from a fully spun-up state of global carbon cycling including an initial balance between the rate of silicate rock weathering and volcanic CO 2 outgassing (see Lord et al., 2016, for details). For the "Carbonate" simulations, the model was run with just carbonate compensation only, as per . For the "Closed" simulations, there was no weathering or sediment response, and hence is just ocean-atmosphere repartitioning. For the "AirOcean" simulations, the carbon remains in the atmosphere. Both the resulting history of atmospheric CO 2 and annual mean global surface air temperature were extracted and calculated as anomalies relative to the late Holocene (pre-industrial).
Using the "Full" setup, cGENIE simulations were run using four different carbon estimations from the GCM simulations (see Table 1). For "Carbonate", "AirOcean", and "Closed" a simulation was run with the GCI_ELE carbon scenario (see Table 1). Therefore, seven transient cGENIE simulations were run in total.
It should be noted that we do not attempt to change the boundary conditions required by the cGENIE model dynamically through the glacial-interglacial cycles, namely orbital parameters, planetary albedo, and sea level (and ocean salinity). These are instead kept fixed at modern (following Lord et al., 2016). Hence, changes in the sensitivity of atmospheric CO 2 to unit CO 2 input (or removal) will not be accounted for. We expect such an effect to exist due to , for example, the dependence of the Revelle factor (the sensitivity of dissolved CO 2(aq) to changes in total dissolved inorganic carbon Zeebe et al., 1999) on both (ocean surface) temperature and atmospheric pCO 2 , changes in ocean circulation and the efficiency of the biological pump, and changes in the carbonate buffering of ocean chemistry. Some of these factors could in theory be imposed (e.g. changes in ocean surface temperatures), but others would require the glacial-interglacial dynamics in both ocean circulation and marine carbon cycling to be sufficiently accurately represented in the model. The latter is far beyond what the current state of understanding of glacial-interglacial global carbon cycling allows for (Kohfeld and ). Hence, our assumption of fixed late Holocene boundary conditions will impart a small bias in our estimates of the atmospheric CO 2 response, but not one that would affect our overall conclusions.
In addition, in making estimates of the mean global air surface temperature change corresponding to the projected change in atmospheric pCO 2 in cGENIE, it is important to also note that the climate sensitivity is effectively prescribed (Edwards and Marsh, 2005). In the Lord et al. (2016) configured used here, only sea-ice cover, via its associated albedo, can provide feedback on climate. In the absence of a dynamical atmosphere, glacial-interglacial changes in climate sensitivity due to changes in atmospheric circulation and clouds are not possible. Nor do we account for the possible influences of changes in total land surface area (from sea-level change) or vegetation cover and distribution. However, the assumption of an effectively fixed climate sensitivity across the glacial-interglacial cycle is unlikely to impart significant bias or unduly affect our overall conclusions.

Results: vegetation dynamics
The changes in climate over time affects the vegetation cover in the Dynamic simulations (shown in Fig. 1 Our results from TRIFFID are consistent with the relative changes, although, since TRIFFID uses the actual climate from the models, the vegetation can have biases (e.g. Australia has a tendency to be too wet in HadCM3 in the present day and hence the coupled model has too much vegetation in this region). However, during glacial times there is a decrease in biomass, consistent with Hoogakker et al. (2016).
Comparison with the BIOME6000 mega-biome maps for LGM (Pickett et al., 2004;Prentice and Jolly, 2000;Bigelow et al., 2003;Harrison et al., 2001) dataset shows general agreement. The model has considerable expansion of grasses in Eurasia, where BIOME6000 has grassland and dry shrubland. Broadly speaking, in both, North America shows little change from the mid-Holocene to LGM. One key weakness of the model is in western Europe, where BIOME6000 shows grassland and dry shrubland, whereas the model has shrubs and needleleaf trees. South-east Asia shows continued warm-temperate, temperate, and tropical forest where our model simulated broadleaf trees, which encompasses all of these biomes. The BIOME6000 reconstructions show around a dozen tundra points on and near the Bering land bridge, and our model simulates this as C 3 grasses, which is the closest PFT to tundra. Over central Asia our model has extensive areas where the dominant land surface type is bare soil, indicating desert or sparse, dry vegetation. BIOME6000 shows a mixture of desert and dry grass/shrubland, which is generally in keeping with the low-productivity, low-density vegetation indicated in our model simulation.
The forest extent in the tropics at the LGM is similar to PI (see Appendix Fig. A2 for shifts in vegetation at 21 ka). This is supported by pollen and other data Anhuf et al., 2006;Wang et al., 2017) and modelling (Cowling et al., 2001) which find that although there is diminished tropical forest, there is still substantial tree cover at the LGM and little sign of widespread grasslands. Conversely, the BIOME6000 data find that the tropical rainforest area was reduced during the LGM (Pickett et al., 2004;Prentice and Jolly, 2000;Bigelow et al., 2003;Harrison et al., 2001) and grasslands expanded, as do some modelling studies (Martin Calvo and Prentice, 2015;Prentice et al., 2011;Hoogakker et al., 2016). It is interesting to note that in the present-day Amazon, BIOME6000 shows three points of tropical forest, two of savanna, two of warm-temperature forest, two of temperate forest, and three of dry grass/shrubland at the LGM. In our simulations the dominant PFT of the same area is broadleaf trees. For comparison, Prentice et al. (2011), using a different dynamic global vegetation model, have tropical forest over the same domain. Therefore, there is little indication that, where TRIFFID may be inconsistent with BIOME6000, another model is necessarily significantly better.
Because of the PFT (rather than biome) approach of TRIF-FID, and the limited number of PFTs, it is difficult to be sure whether trees in the tropics are a tropical rainforest at the LGM, as there are a number of biomes with significant amounts of trees. Although there is little change in PFT in the tropics at the LGM, on the margins there are reductions in vegetation carbon, suggesting a change in vegetation within the large margins of the PFTs used in this model.

Results: biogeophysical feedbacks
The biogeophysical impacts of vegetation are calculated by subtracting the Dynamic simulations from the corresponding Static simulations. We find that vegetation is acting as a positive feedback to the climate, enhancing the cooling (Fig. 2a). Broadly, the Static and Dynamic simulations both agree with an approximation of global temperature over the whole period (the EPICA dataset halved) (Fig. 2a). The Static set generally does better in 70 to 10 ka, whereas the Dynamic set is closer to the EPICA data in the period 110 to 70 ka. The biogeophysical differences between the Static and Dynamic sets alter global, annual mean surface temperature by as much as −0.91 • C (see Fig. 2b). Regionally this temperature cooling is up to 5 • C (Fig. 3).
The albedo changes are in the same location as the vegetation carbon changes, as well as the main temperature changes (see Figs. 3 and 4). These temperature differences are mainly driven by reductions in tree cover and its replacement with bare soil or grasses, which is a result of the vegetation dynamics in the model (see Figs. 1 and 4). Trees have a lower albedo, and when they are replaced by higher albedo grasses, there is a cooling effect. The change in tree fraction between the Static and Dynamic sets is a good predictor of the temperature changes (r 2 = 0.79 using a linear model of the global temperature and tree anomalies). This is exacerbated by the presence of snow cover as the snow-covered visible and nearinfrared albedo of grasses, shrubs, and bare soil is higher than that of trees (Essery et al., 2001). Therefore, when trees are replaced by grasses where there is snow cover for part of the year, there is a larger change in albedo than where there is no snow cover. Thus, the albedo changes can be seen mainly where a change between trees and grasses occurs in an area with snow cover (see Figs. 1 and 5). The exact contribution of the snow as opposed to the no-snow albedo is difficult to disentangle, but the influence of this effect is well established (Betts, 2000).
The land surface albedo changes caused by the vegetation have an even stronger correlation with these biogeophysical temperature changes (r 2 = 0.86). However, we can see that although the forcing is land-based (the dynamic vegetation), significant changes occur in the ocean (see Figs. 3 and 4) that drive the resulting temperature changes. Ocean-only surface albedo anomaly as a determinant of global temperature anomaly has an r 2 of 0.95 -lower only than the r 2 of the global (land and ocean) surface albedo of 0.96. By comparison, the r 2 of the latent heat anomaly as a predictor of temperature anomaly is lower for land, ocean, and globally than surface albedo (0.70, 0.93, and 0.91, respectively). The other parts of the energy balance, in particular the latent heat, sensible heat, and the net shortwave radiation, do not have such a clear relationship with the temperature change (see Appendix Fig. A4 and compare to Fig. 3).
Although the biogeophysical changes cause cooling, there are some minima of biogeophysical temperature change seen at 30, 56, and 100 ka (Fig. 2, filled symbols). These minima have an oceanic source and are caused by vegetation interacting with thermohaline circulation changes. In our new simulations we account for the net transport of water from ocean to the ice sheets by a parameterisation that instantaneously balances any net accumulation of water on ice. This parameterisation results in fresher ocean conditions during times of precession-driven N. Hemisphere summer insolation highs (less water is being used to build the ice sheets). The instantaneous nature of the parameterisation is physically unrealistic but reductions in accumulation and an increase in ablation during precession highs has been seen in fully coupled climate-ice sheet EMIC simulations (e.g. Ganopolski et al., 2010). During weaker accumulation periods, the parameterisation results in a freshening of ocean surface waters and a reduction in Atlantic meridional overturning circulation (AMOC) strength from 16 to 10-12 Sv. Superimposed upon this general behaviour, the addition of interactive vegetation generally does not change the AMOC strength. However, at times of weak AMOC, small changes in runoff and temperature are sufficient to cause some changes in the response. For instance, in the static vegetation simulations there is a relatively weak AMOC in the simulations for 60, 58, and 56 ka. In the interactive vegetation simulations, the weakened AMOC only occurs at 60 ka. Thus, at 60 ka the changes in climate are fairly typical of preceding times but at 58 and 56 ka there is a substantial difference between the static and dynamic vegetation simulations. The cause for this difference is associated with a combination of reduced runoff into the North Atlantic (principally from changes in land surface in North America) and colder temperatures, both of which act to stabilise the AMOC in all three periods, but it is sufficient to prevent the AMOC weakening in the 58 and 56 ka simulations. This threshold-like behaviour of the AMOC is almost certainly highly modeldependent and hence the result is not robust.
The regional patterns of cooling also temporarily affect the precipitation regime (see Fig. 3). This appears to be related to the AMOC weakening. There are some suggestions of similar relationships between the increases in precipitation and the terrestrial changes to previous studies (Gedney and Valdes, 2000;Singarayer et al., 2009). Similar to the temperature changes, it is unclear how model-specific these changes are.

Results: biogeochemistry
We now calculate the total change in terrestrial carbon stores in the HadCM3 simulations. We consider scenarios of terrestrial carbon change with combinations of including or excluding uncertain aspects of the carbon cycle, specifically depending on the fate of soil carbon under ice and the changes related to the expansion of land. Zeng (2003) suggested that the soil and vegetation carbon formed during the warm last interglacial could simply get covered by ice and is stored there, rather than being released into the rest of the system as is typically assumed in past estimates. Equally, if the carbon were exported from under the ice sheets, it is not certain that it would return to atmosphere, rather than be sequestered in riverine and oceanic systems.
Similarly, the amount of carbon stored on newly emerged land is also uncertain as it depends on both the area of emergent land and the surface properties. It is uncertain whether the vegetation and soil inundated in the deglaciation would have been returned to the atmosphere or remained in situ. Therefore, we calculate the changes in soil and vegetation carbon from these various sources. In Table 1 we focus on the changes between pre-industrial and LGM, which corresponds to the largest overall change through the glacialinterglacial cycle.
In  In reality, glacial systems are known to export carbon in a highly labile form (Lawson et al., 2014); erode soil and bedrock, creating major landscape changes; and release large amounts of methane when they retreat (Wadham et al., 2012). Although the conversion of this terrestrial carbon to atmospheric carbon may be through riverine or oceanic systems, it seems likely the majority would return to the atmosphere within the time periods we consider. We therefore use this largest scenario as a conservative option for our main analysis.
The other major change to soil carbon in the model is newly exposed land, which is revealed when the water in the ice sheets causes lower sea levels (see Figs. 1 and 5). For the new land we use a nearest-neighbour interpolation of basic soil properties (e.g. water-holding capacity) and the model is run for sufficient length of time for the soil and vegetation carbon to reach equilibrium.
This estimate of carbon on expanded lands also has uncertainties. We have reasonable confidence in the sea-level estimates and consequent change in land area, but it is much more uncertain whether carbon could accumulate on that land. Uncertainty in the positions and timings of ice sheets, particularly before the LGM exacerbate this. For instance, in our simulations the East Siberia ice sheet is absent (see Fig. 1), whereas some other ice-sheet reconstructions include it (e.g. Niessen et al., 2013). The area of the ice sheet alone accounts for an average of 56 PgC soil carbon in these simulations when it is absent. However, soil carbon takes a long time to accumulate, especially with low net primary prowww.clim-past.net/13/1381/2017/ Clim. Past, 13, 1381-1401, 2017 ductivity (NPP) and vegetation carbon storage averages just 0.5 PgC over all the expanded lands. If exposed land carbon was included and glacial land soil carbon excluded, the terrestrial carbon is a gain from PI to LGM of +37 PgC (see Table 1). However, as discussed above, we would argue that excluding glacial land soil carbon change is probably unreasonable. Most previous studies have also assumed that all carbon under ice is removed. If we include the loss of carbon, then the range in total amount of terrestrial carbon lost in this model between the PI and the LGM at 21 ka is −440 to −257 PgC.
The change in terrestrial carbon found in our simulations contributes to atmospheric carbon dioxide change. Using the cGENIE model to approximate the carbon uptake by the ocean we therefore calculate the atmospheric carbon dioxide change (see Methods section and Fig. 6).
Selecting the largest change in carbon storage (−440 PgC at the LGM, including glacial soil carbon changes and ex-  Figure 6. Contribution of terrestrial biosphere carbon emissions to atmospheric CO 2 . Run with the cGENIE "Full" configuration and normalised to pre-industrial CO 2 levels. The four scenarios are as detailed in Table 1. cluding expanded lands), we find that the results suggest a peak contribution compared to pre-industrial CO 2 of 25 ppm CO 2 (Fig. 6). In all scenarios except GCE_ELI, the terres-Clim. Past, 13, 1381-1401, 2017 www.clim-past.net/13/1381/2017/ trial carbon contribution to atmospheric CO 2 acts as a negative feedback to the climate, dampening the effect of other climate forcings, including the net contribution of the terrestrial biosphere (Fig. 2b). Within cGENIE, the change in atmospheric CO 2 produces a warming at the LGM of 0.29 • C (equivalent to a climate sensitivity of around 2 W m −2 • C −1 ; see Fig. 2b). This is much smaller than the biogeophysical contribution of −0.84 • C. It is also much less variable. For most of the glacial period, from 100 to 20 ka, the implied biogeochemical warming is around 0.26 • C (Fig. 2b). This results in the dominance of the biogeophysical impacts over biogeochemical feedbacks.
However, with different Earth system processes included, the biogeochemical effects vary substantially (see Fig. 7). In the simulations discussed above, silicate and carbonate weathering are both included and this results in the lowest temperature change from the same carbon emissions. The temperature contribution at the LGM increases relative to the baseline "Full" simulation (0.29 • C) as the silicate weathering is excluded (Carbonate, 0.30 • C); all weathering is excluded (Closed, 0.47 • C); a decadal-to millennial-scale carbon uptake is used (TRCE, 0.86 • C); and if all carbon remains in the atmosphere (AirOcean, 1.92 • C). Note that the TRCE as shown above includes the terrestrial biosphere as a sink, so it will slightly overestimate how much carbon will be removed from the atmosphere when the source is the natural vegetation. Comparing these values to the biogeophysical terrestrial effect in Fig. 2b, we can see that the shorter the timescale, the more likely biogeochemical terrestrial processes will dominate as it weakens over time. On longer timescales the biogeophysics dominate because the scale of the effect does not diminish over time relative to the control.

Discussion
Considering the biogeophysical and biogeochemical effects of the terrestrial biosphere together, we find that the biogeophysics dominate the signal. Over the entire 120 kyr period considered, the average temperature contribution by the terrestrial biosphere is −0.26 • C. At the LGM, the contribution of the terrestrial biosphere to temperature is −0.58 • C. However, if we use the GCE_ELI scenario, where the carbon under the ice sheets and on land exposed by sea-level changes is assumed not to return to the atmosphere, the equivalent values are −0.50 and −0.83 • C. For both the conservative and more uncertain assumptions, it is clear that terrestrial carbon is a much smaller contributor to temperature than biogeophysics over this period.
The biogeophysical results found here broadly concur with comparable model studies of past vegetation biogeophysics. Claussen et al. (2006) found the biogeophysical contribution of vegetation to LGM cooling of around 1 • C in the Northern Hemisphere, whereas Jahn et al. (2005) found around  (2007); "Closed", with no weathering or sediment response, and hence just ocean-atmosphere repartitioning; and "AirOcean", where the carbon remains in the atmosphere. "TRCE" is the simple calculation of the TRCE of HadCM3 (taken from Matthews et al., 2009) for the same carbon inputs into the atmosphere as used for the cGENIE simulations. Note that we include the TRCE for completeness, but it is not a cGENIE simulation (see Methods section).
−0.6 • C, and up to 2 • C locally. Our result of −0.84 • C is in the middle of the other LGM studies. The dominance of the biogeophysical effects found here is contrary to the results found for short timescale problems, which find that biogeochemistry tends to be comparable in magnitude to biogeophysical effects (e.g. Davies-Barnard et al., 2014b;Pongratz et al., 2010). This is because the centennial simulations have a stronger biogeochemical effect since the transient response to cumulative emissions is stronger than the equilibrium response. In climate simulations up to around a century long, more carbon tends to remain in the atmosphere. This creates a strong warming effect that is approximately linearly related to the amount of greenhouse gas emissions (Matthews et al., 2009;Gillett et al., 2013). The transient response to cumulative emissions (TRCE) accounts for the uptake of atmospheric carbon by the ocean and terrestrial biosphere, but only on short timescales. The uptake of atmospheric carbon by the ocean requires hundreds or thousands of years, and it is slower when the increase in carbon into the system is small and staggered (Lord et al., 2016). However, the simulations we use are on a millennial timescale, allowing much of the carbon to be taken up by the ocean (Lord et al., 2016). From a climate sensitivity point of view, this means that, on shorter timescales, the effects of dynamic vegetation can cancel each other out. This provides some rationale for the fact that dynamic vegetation has been generally not included in the majority of state-of-the-art Earth system models used in CMIP5, as it does not significantly affect the climate sensitivity. At longer timescales, it is more important to include dynamic vegetawww.clim-past.net/13/1381/2017/ Clim. Past, 13, 1381-1401, 2017 tion, as without the positive feedback of the biogeophysical effects, the climate sensitivity would be underestimated. For the biogeochemical effects of the terrestrial biosphere, previous estimates of carbon stocks on exposed continental shelves based on models are between 112 and 323 PgC at the LGM (Montenegro et al., 2006). The comparable number in this simulation is 183 PgC, which is on the lower end of the wide range of other models. However, it has good agreement with the vegetation reconstruction (not model) values by Montenegro et al. (2006) of 182 to 220 PgC.
For present day, Hugelius et al. (2014) shows around 75-100 kg C m 2 of soil carbon in far north Siberia and 20-40 further south. Far northern Canada is much more heterogeneous, with values from 20 to 150 kg C m 2 . The modelled PI values of northern high-latitude soil carbon are on the low side, and much more homogeneous, around 15-20 kg C m 2 , but similar to Hugelius et al. (2014), who show far north America to be less consistent, with some higher areas of 35-40 kg C m 2 in the far north. (See Appendix Fig. A4 of the loss of soil carbon under ice sheets at the LGM.) What this suggests is that, while of the correct order of magnitude, the model has a very modest amount of soil carbon that could be considered permafrost. Therefore, we think it is reasonable to include this low estimate of soil carbon in the uncertainties.
The soil carbon change under ice sheets between PI and LGM is modelled as ∼220 PgC. Extrapolating from a comparison with Hugelius et al. (2014), this might be a third too little. If the true value were ∼330 PgC, this would make the total C change PI to LGM 550 PgC, more in line with some previous model estimates. It would affect the global mean annual biogeochemical contribution by ∼ 0.1 K. This would mean the net effect of vegetation was closer to zero, but the biogeophysical effect would still dominate.
However, the exact size of the terrestrial carbon emissions is uncertain. Other carbon stores not accounted for here are potentially important, for example methane during sea-level rises or changes to the wetlands in the tropics. Modelling studies that look at wetlands at the LGM suggest that although the wetland area is larger, the methane emissions are lower compared to modern day (Kaplan, 2002). However, palaeohydrological data indicate a drying in the African tropics (Gasse, 2000). Our model does not have a processbased permafrost or wetlands component, and therefore the changes in methane are not accounted for. This is a particular limitation when considering the carbon stored in deep permafrost soils in northern peatlands. Saito et al. (2013) show that, based on the temperature changes, there is a substantial expansion of permafrost area during glacial times but cannot estimate any changes in carbon storage. Zimov et al. (2006Zimov et al. ( , 2009 have argued that permafrost storage could be a major source of carbon through the deglaciation, and Ciais et al. (2012) argue that there was a large extra pool of inert carbon at the LGM. Similarly, Köhler et al. (2014) have argued that large amounts of carbon were locked into permafrost which were then released rapidly at the Bølling-Allerød.
Research has also suggested that waterlogging and flooding as sea-level rises during the Holocene could cause rapid anaerobic decomposition of vegetation, causing methane emissions . This could account for emissions of as much as 25 PgC for 10 m sea-level rise (ibid). Since our simulations do not account for methane or this effect of inundation, it is likely it there is a slight underestimation of equivalent CO 2 effect of the carbon emissions (as methane is a stronger greenhouse gas than carbon dioxide).
The impacts are mainly determined by the vegetation shifts the dynamic global vegetation model simulates. Each grid box has the potential for five PFTs, but generally the Lotka-Volterra equations used in TRIFFID mean that the grid box is dominated by one PFT. The small number of PFTs means that the range within each PFT is relatively large. Therefore the model probably underestimates the effects of small perturbations in climate, as the large definition of the PFTs allows the PFT to remain the same. Conversely, it makes an abrupt change more likely as the climate tips a grid box from being predominantly one PFT to being predominantly another. Overall, the model could be slightly underestimating the amount of change in vegetation. However, because of the ratio of the biogeophysical to biogeochemical changes, if the vegetation change is underestimated, the sign of the net effect of the terrestrial biosphere is unlikely to change. Similarly, because over the long time periods involved much of the released carbon is taken up by the ocean, the changes in carbon densities of the vegetation would need to be wrong by a lot to change the overall signal.
Our approach here assumes that there is no non-linear interaction between the biogeochemical and biogeophysical effects. Since the biogeochemistry acts as a negative feedback and reduces over time, and the biogeophysics acts as a positive feedback and stays the same over time, there is no strong reason to believe that, at equilibrium, there would be any significant synergy. However, on shorter timescales and on a regional rather than global scale, it is quite possible that there could be some synergies.

Conclusions
Using a fully coupled atmosphere-ocean-vegetation model with static and dynamic vegetation, we find that over the Clim. Past, 13, 1381-1401, 2017 www.clim-past.net/13/1381/2017/ last 120 kyr the net effect of vegetation feedbacks on global, annual mean 1.5 m air temperature is a cooling, which can be as much as −0.66 • C (Fig. 2b). For the vast majority of the last glacial-interglacial cycle, cooling associated with biogeophysical feedbacks dominates over the biogeochemical warming associated with reduced terrestrial carbon storage. The biogeophysical cooling effect is mainly due to the role that vegetation plays in changing surface albedo and particularly related to snow cover and the taiga-tundra transition (Gallimore and Kutzbach, 1996;de Noblet et al., 1996) and is a well-understood climate feedback. The biogeochemical contribution to atmospheric carbon dioxide is small (∼ 20 ppmv) and hence the temperature contribution is small (on average 0.26 • C with a maximum of 0.33 • C).
There are significant uncertainties in this calculation which would further diminish the net temperature impact of the terrestrial biosphere by cancelling out the biogeophysical impact. In this analysis, the only time periods where the effects are comparable are at times when additional mechanisms operate, such as changes in ocean circulation, but these mechanisms may be model-specific.
The key uncertainties in this study originate in the biogeochemistry, especially the soil carbon build-up in newly exposed land, the fate of soil carbon in glacial systems, and the amount of carbon in permafrost (not calculated in this study). Further research is needed to fully understand the functioning of these systems and how they can be best incorporated into climate models. In addition, the technique we use for inferring the biogeochemical effects of terrestrial carbon changes has limitations and is potentially model-dependent. However, the smaller estimate of terrestrial carbon emissions may make the low LGM atmospheric carbon dioxide somewhat easier to reconcile (Montenegro et al., 2006).
Our work confirms previous results using EMICs (Brovkin et al., 2012) that found the net terrestrial biosphere effect to be primarily biogeophysical and that the terrestrial carbon contribution to atmospheric carbon is comparatively small. Our findings also represent a clear illustration of the net climatic effect of vegetation being highly dependent on the timescale, with the biogeophysical response dominating in the longer term in contrast to century-scale future changes.
Code and data availability. The model code is currently available to view at http://cms.ncas.ac.uk/code_browsers/UM4.  LGM and PI in the Dynamic simulations, for the five PFTs: pft1 is broadleaf trees; pft2 is needleleaf trees; pft3 is C 3 grasses; pft4 is C 4 grasses; pft5 is shrubs.