An improved land biosphere module for use in reduced complexity Earth System Models with application to the last glacial termination

Interactions between the land biosphere and the atmosphere play an important role for the Earth’s carbon cycle and thus should be considered in studies of global carbon cycling and climate. Simple approaches are a useful first step in this direction but may not be applicable for certain climatic conditions. To improve the ability of the reduced-complexity Danish Center for Earth System Science (DCESS) Earth System Model DCESS to address cold climate conditions, we reformulated the 5 model’s land biosphere module by extending it to include three dynamically varying vegetation zones as well as a permafrost component. The vegetation zones are formulated by emulating the behavior of a complex land biosphere model. We show that with the new module, the size and timing of carbon exchanges between atmosphere and land are represented more realistically in cooling and warming experiments. In particular, we use the new module to address carbon cycling and climate change across the last glacial transition. Within the constraints provided by various proxy data records, we tune the DCESS model to a Last 10 Glacial Maximum state and then conduct transient sensitivity experiments across the transition under the application of explicit transition functions for high latitude ocean exchange, atmospheric dust, and the land ice sheet extent. We compare simulated time evolutions of global mean temperature, pCO2, atmospheric and oceanic carbon isotopes as well as ocean dissolved oxygen concentrations with proxy data records. In this way we estimate the importance of different processes across the transition with emphasis on the role of land biosphere variations. 15

ocean vertical diffusion profile function in the high latitude model ocean that is applied for generating Last Glacial Maximum (LGM) conditions in the DCESS model and the formulation of the resumption of deep ocean mixing.Furthermore, we present some additional information on the generation of model LGM conditions in atmosphere and ocean, a literature review on the Mystery Interval (MI) and an overview for model and proxy data changes during the MI.The Supplement ends with additional information on the various 14 C production rate time series that were applied in the model simulations, an analysis of the DCESS 10 model calculated 14 C production rate time series and the isotope ratio definitions.

S1
In the DCESS model, albedo, α, is taken to be constant and equal to 0.62 for all snow or ice covered areas.For non-snow/ice covered areas.For non-snow/ice covered areas, α is expressed as where Θ is the latitude, b = 0.175 and a = 0.3 for present day conditions.This functional form and these constant values have been based on present day observations (Hartmann, 1994).The albedo of non-snow/ice covered areas should vary with vegetation type since forested areas have lower albedo than non-forested areas (Bonan, 2008).As seen in Fig. 2, as the Earth cools from present day, both forested model areas (EF and TF zones) contract while the non-forested model area (GSD zone) expands slightly, in part in response to dryer conditions (Gerber et al., 2004).This would lead to higher albedo and a positive feedback on the cooling.For completeness in our new treatment of the role of the land biosphere in climate and to capture such albedo variations within the context of our new land biosphere module, we assume that a in Eq.S1 may be related to vegetation type such that where the factor 0.3 is the present day value of a, γ is a multiplier, the value of which is determined by calibration (see below), f rac is the ratio of the area of the GSD zone to the total non-snow/ice covered area (i.e. the sum of the areas of the EF, GSD and TF zones) and f rac 0 is this ratio for present day.Note that f rac(δT glob ) can be taken from Fig. 2 or calculated explicitly using Eqs.S1 and S2 and the snowline/ice sheet dependency on δT glob (see Fig. S4 below).Fig. S1a shows a plot of f rac(δT glob )/f rac 0 .

S2
The vegetation albedo forcing for the LGM (δT glob = −3.5 • C) relative to present day has been determined in more complex models from which we choose the value of −0.7 W/m 2 as being representative (Köhler et al., 2010).Together with Eq.S2 and the model latitudinal distribution of solar forcing, we find that this LGM vegetation albedo forcing anomaly is obtained in our model simulation for a γ value of 0.02, a value we adopt here.Fig. S1b illustrates new albedo distributions with latitude for the specific cases of δT glob = −4, 0 and 4 • C for which a = 0.3027, 0.3 and 0.2976, respectively.

S3 Transition function for dust concentration
As described in the manuscript, the dustier atmosphere during LGM climate conditions (e.g.Maher et al., 2010) leads to cooling as forced by an additional radiative forcing (A Dust ) of −1 W m 2 (see Mahowald et al., 2006) and increased iron-fertilization, addressed in the model by an increased efficiency factor (f F e−lim ) in the high latitude ocean (Lambert et al., 2013(Lambert et al., , 2015)).For the transition of these parameters between LGM and PI conditions in the transient simulations, we applied functions derived from the correlation between dust flux and global mean atmosphere temperatures in proxy data records.Fig. S2 shows proxy data records from Antarctica of dust flux and temperature (Jouzel et al., 2007;Lambert et al., 2012).As stated in Lambert et al. (2012), the major changes in dust flux occurred before 15 kaBP.An analysis of the correlation of the two quantities yields an exponential dependency of dust with temperature for both high and low latitudes that can be described by where a and b are constants to be determined..The same climate-dust exponential relationship also holds true for a remote site from the central Pacific (Winckler et al., 2008).
Figure S2.Dust fluxes and Antarctic temperature anomaly across the last glacial termination (Jouzel et al., 2007;Lambert et al., 2012).for the additional dust radiative forcing can be established.In Fig. S3, these functions are plotted for the range of T glob considered here.

S4
Fig. S4 shows the ice sheet line (converted as described in the manuscript from the data by Shakun et al., 2012) and the snowline (latitude of T glob = 0 • C) as calculated by the transient model simulation with all transition functions (the All_TF simulation) from Sect.3.2.The thick parts of the lines mark the current minimum of the two functions, which is the value that is applied in the calculations at the actual time step (see Sect. 2.3).This means that in this simulation, the ice sheet line has no influence on the results between 0 and 16 kaBP and the same is true for the snowline between 16 and 25 kaBP.Due to the transition of the minimum between the functions at around 16 kaBP, some results may show a kink there.
The equatorward displacement of the model snowline and the expansion of terrestrial ice sheets lead to three important vegetation-climate feedbacks in this DCESS model version.Firstly, they influence the calculation of the global land albedo, leading to a positive cooling feedback.Secondly, they limit the poleward expansion of the EF vegetation zone.This decreases In the article we present a profile for the high latitude ocean vertical diffusion with depth, which is applied to generate an analogy to isolated Southern Ocean (SO) deep and intermediate waters (see e.g. Watson and Naveira Garabato, 2006).In this section, we present details about this profile and explain how it was established.The mathematical formulations of the profile and of the resumption to full vertical diffusion in the transient simulations are given.

S5.1 Description of the profile
For PI climate conditions, the ocean vertical diffusion parameter is set to D P I = 2.3 • 10 −3 m 2 s −1 evenly throughout the high latitude DCESS model ocean.To generate isolated deep water and thereby LGM climate conditions, we impose the function on this parameter.Here, d denotes the ocean depth (in m) and the factor 0.3 was chosen as a best-guess to set the steepness of the transition. and are used to calculate the maximum and the minimum of the diffusion profiles in relation to the standard value for PI conditions Here, p min represents the relative minimum of the diffusion, which is set to 0.03, and R C determines the change of diffusivity above the transition.However, the minimum of the vertical diffusion by depth is determined by the low latitude ocean value (Shaffer et al., 2008)

S5.2 Resumption of deep ocean mixing
In order to simulate the resumption of deep ocean mixing during the MI, we impose a steady deepening of the transition depth of the imposed vertical diffusion profile.Starting at year 17.5 kaBP we modify p a in Eq.S6 by where n denotes the time step after the start of the resumption.This steadily decreases the transition depth of the vertical diffusion ocean profile over time.Effectively, it leads to the resumption of the high latitude ocean to PI mixing within around 3 ka.In other words, this simulates time-dependent upwelling of the deep ocean waters that had been isolated through the imposed vertical mixing profile.Thus, it approximates the sudden event of deep ocean mixing during the MI, that has been hypothesised in a number of studies (e.g.Burke and Robinson, 2012).

S6 Additional information on model LGM and transition
The long wave radiative forcing (RF) in the model follows the approximation RF = A − B • T glob by Budyko (1969), where T glob denotes the global mean atmosphere temperature (in One of the features of the DCESS model is that it is an open system model that considers input from weathering and volcanism and outputs from burial.There is no clear consensus with regard to weathering changes since cooler LGM temperatures may favour less weathering while lower sea level and exposed shelves together with the action of the ice sheets may favour more weathering (e.g.Vance et al., 2009).Furthermore, there is evidence that volcanism may have been affected by lowering sea level and ice sheet loading (Kutterolf et al., 2012).Given this state of affairs and for simplicity, we have chosen to maintain weathering and volcanism constant at the values found for the PI calibration of (Shaffer et al., 2008).
We implemented the transition functions for the ocean phosphate and salinity in the following manner: First, we associate a sea level of −130 m below present day values (as during the LGM) to a phosphate concentration scaled by f P = 1.035.The temporal change of the sea level across the last glacial termination then linearly determines the phosphate concentrations so that f P = 1.0 is reached, when the sea level is at present day values (0 m). Analogously, the ocean salinity is prescribed with a value of 35.9 psu for LGM sea level (−130 m) and 34.7 psu for present day (see Adkins et al., 2002).
An overview of all model parameters that were modified additionally to the ocean vertical diffusion is provided in Tab.S1.

S7 Model LGM ocean profiles
The ocean profiles for LGM conditions of various variables for the high and the low-mid latitude sector are presented in Fig. S6.
For reference, also observed present day low-mid latitude ocean profiles (Shaffer et al., 2008) are included in the figures.The present day calibrated DCESS model ocean profiles can be found in Shaffer et al. (2008).
A prominent feature in these vertical ocean profiles is generated through the transition of the vertical diffusivity to very low values below around 1800 m depth.Below the transition, the ocean shows high DIC concentrations and strongly depleted 13 C isotope ratios (see Curry et al., 1998;Mackensen et al., 2001;Hesse et al., 2011).For reference, we provide all LGM oceanic sediment δ 13 C measurements south of 46 latitude ocean below the vertical diffusion transition also shows low oxygen concentrations with minima of 0.05 mol m −3 in the deep ocean.This indicates that there is no widespread ocean anoxia in accordance with proxy data (Jaccard et al., 2014).
Even though the isolated deep ocean waters are strongly 14 C-depleted compared to the upper ocean at LGM model conditions, deep ocean ∆ 14 C values are still around 150‰ higher than at PI conditions.However, since atmospheric ∆ 14 C has decreased by more than 400‰ from LGM to PI conditions, the difference between atmospheric and oceanic ∆ 14 C is much higher during the LGM.max(∆∆ 14 C) (the maximum of ∆ 14 C differences between atmosphere and ocean) is around −480‰ during the LGM and also agrees well with observations (Burke and Robinson, 2012).Tab.S2 provides an overview of the result of key model variables of atmosphere and ocean for PI and LGM conditions as well as the proxy data LGM values for comparison.

S10
S8 Background information on the Mystery Interval One leading hypothesis for explaining the initial step of the last glacial termination ("Mystery Interval" (MI), from 17.5 to 14.5 kaBP; Broecker and Barker, 2007) is enhanced upward exchange of and subsequent outgassing from SO deep water (Francois et al., 1997;Sigman and Boyle, 2000).This water had been isolated from surface layers and thereby had been accumulating carbon from remineralisation of organic matter during the previous glacial time (Broecker and Barker, 2007).
Such enhanced transport could result from a number of processes.These include a strengthening and/or poleward shift of the Southern Hemisphere west wind belt that would enhance deep SO upwelling (Toggweiler and Russel, 2008;Anderson et al., 2009;d'Orgeville et al., 2010) and changes in ocean stratification through brine-induced effects on glacial-interglacial time scales that would enhance deep SO convection and mixing (Bouttes et al., 2010(Bouttes et al., , 2011;;Mariotti et al., 2013Mariotti et al., , 2016)).In addition to such changes, a number of other processes, such as variations in iron fertilisation of the SO, ocean volume changes, as well as carbon storage in the terrestrial biosphere (e.g.Kohfeld and Ridgwell, 2009) and in permafrost (e.g.Zech, 2012;Köhler et al., 2014;Crichton et al., 2016), have to be considered.
Recently, consistent records of carbon isotopes from ice cores have become available for use in constraining processes that control the global carbon cycle.Carbon isotope ratios are sensitive to variations in carbon exchange between the atmosphere, the terrestrial biosphere and the ocean reservoirs, but also to other effects on geological time scales such as changes in oceanic biogeochemical conditions.Therefore, analyses of the variations of isotope ratios can provide deeper insights into the mechanisms that are responsible for the pCO 2 and hence temperature increase.For example Reimer et al. (2013) and Schmitt et al.
(2012) presented time series of reconstructions of the two carbon isotope ratios ∆ 14 C and δ 13 C in the atmosphere during the past 25 ka.Both ∆ 14 C and δ 13 C show a relatively sharp drop during the MI.While ∆ 14 C continues to decrease after that until present day, δ 13 C rises again during the Holocene after first practically stagnating for around 4 ka.Schmitt et al. (2012) suggest that δ 13 C is significantly influenced by ocean temperature changes only during the stagnating phase.The sharp drop can be related to the outgassing of isolated deep waters and the regrowth of the terrestrial biosphere can account for the increase during the Holocene.∆ 14 C is largely affected by its cosmogenic production rate (see e.g.Lal and Peters, 1967), but estimates of 14 C production rates indicate that air-sea carbon exchange has to be considered to explain the ∆ 14 C drop during the MI (Laj et al., 2004;Muscheler et al., 2004Muscheler et al., , 2005;;Hain et al., 2014).Burke and Robinson (2012) used deep sea coral data to reconstruct concentrations of this radioactive tracer (Half-life: T 1/2 ( 14 C) ≈ 5730 a) in different oceanic water masses.
Isotopically strongly depleted water masses can be found in the deep and intermediate SO (Sikes et al., 2000;Skinner et al., 2010;Thornalley et al., 2011).Due to the slow radioactive decay of 14 C, isolation of these water masses from surface waters can account for this isotopic depletion during the last glacial.At the onset of the last glacial termination, ∆ 14 C in deep ocean waters rapidly increased, which indicates a sudden event of deep ocean mixing, transporting isotopically depleted waters from the deep to the upper ocean regions and hence also supporting the hypothesis of enhanced SO vertical transport (Burke and Robinson, 2012).Thereafter, deep ocean ∆ 14 C does not change significantly until the early Holocene.This indicates that the increase of atmospheric CO 2 in the latter part of the deglaciation is not primarily driven from the deep ocean.
Nevertheless, several studies call into question enhanced oceanic vertical transport as the conclusive reason for atmospheric CO 2 and temperature variations during the MI.Due to limitations of the volume of an abundant isolated ocean reservoir, Broecker and Barker (2007) doubt that the magnitude of such an event could account for the observed atmospheric CO 2 changes.Mainly because of the lack of low ∆ 14 C in their measurements De Pol-Holz et al. (2010), Broecker and Clark (2010) and Cléroux et al. (2011) question the existence of such strongly depleted ocean reservoirs.Moreover, modelling work by Hain et al. (2011) report shortcomings of the isolated reservoir hypothesis.Amongst other reasons the study claims that the isolation causes wide-spread anoxia in the deep ocean (which was not reported by proxy data, see e.g., Jaccard et al., 2014) and that the carbon signal would rapidly be dissipated and diluted in the rest of the ocean.Other studies question the SO to be the origin of the event.Okazaki et al. (2010) argue that North Pacific water masses could help to explain the release of old carbon during the last glacial termination and Kwon et al. (2012) point out the possible influence of deep North Atlantic water masses for this process.In an Earth System Model study, Huiskamp and Meissner (2012) find old water masses in different oceanic reservoirs in response to variations of the mechanisms for deep water formation.

S9 Model Mystery Interval changes
To assess the impacts of the individual transition functions on atmospheric T glob , pCO 2 , δ 13 C and ∆ 14 C changes, four transient simulations have been conducted in the following manner: First, only the prescription of the two minor greenhouse gases methane and nitrous oxide, as well as the variations in ice sheet expansion are activated (simulation: Min_RF, Minor radiative forcing).In the second simulation, the ocean vertical mixing is additionally restored as described above (Add_vdiff, additional resumption of ocean vertical diffusion).In the third simulation, variations in atmospheric dust concentrations are simulated on top of the other two effects through their impact on the radiative forcing and the iron limitation factor (Add_dust, additional application of the dust transition functions).In the final simulation with all changes, the ocean phosphate inventory and the ocean salinity transition functions are applied as well (All_TF, use of all transition functions).The four panels in Fig. S7 show the four atmospheric variables and data-based reconstructions for these four simulations between 20 and 10 kaBP, including the MI.
As expected, changes in only methane, nitrous oxide and the ice sheet extent (Min_RF) have almost no influence on atmospheric pCO 2 , and the temperature increase is also modest.The growth of the terrestrial vegetation through the prescribed ice sheet retreat even leads to a decrease in pCO 2 until about 16 kaBP.Up to this time, the latitude of the ice sheets extends farther equatorward than the snowline (see Fig. S4), note that in the Min_RF simulation the two lines cross earlier because of the slower increase of the snowline latitude).δ 13 C and ∆ 14 C also show only minor changes, which are due to temperature and vegetation changes for δ 13 C and mainly production rate-driven for ∆ 14 C.The additional resumption of the ocean high latitude vertical diffusion (Add_vdiff) generates a pCO 2 change of around 18 ppm during the MI.Furthermore, δ 13 C and ∆ 14 C drop within these 3 ka by an additional 0.14‰ and 60‰, respectively, due to the outgassing of the isotopically depleted deep waters.
The dust-dependent iron limitation factor, as well as radiative forcing changes that are switched on in the third simulation (Add_Dust) again enhance the changes of all four quantities during the MI.However, the changes of all four quantities only reach about half of the proxy data levels.The additions of the ocean phosphate and ocean salinity transition functions (All_TF) do not resolve this apparent shortcoming of the model simulation.These modulations only lead to very small changes.
For the quantification of the impact of individual processes on atmospheric changes during the MI in our modelling concept, Tab.S3 gives an overview of the differences of the four atmospheric variables between 17.5 and 14.5 kaBP in the model simulations presented in the manuscript and in proxy data records.Table S3.Differences of atmospheric variables between the beginning (17.5 kaBP) and the end (14.5 kaBP) of the MI in various DCESS simulations and in data based reconstructions.* Since δ 13 C increases again after a local minimum within the MI in the proxy data, we here provide the minimum value, rather than the value from the end of the MI (−0.24‰).

S10 Details on 14 C production rate determination
The 14 C production rate production rate is determined by the strength of the Earth's magnetic field, which shields the atmosphere from high energy cosmic particles, and the solar modulation of the incidence of high energy cosmic particles (Lal and Peters, 1967;Hain et al., 2014).Here, we present the three cosmogenic 14 C production rate time series from the studies by Muscheler et al. (2004); Laj et al. (2004) and Hain et al. (2014) that were applied in our experiments.Based on the GLobal PaleoIntensity Stack (GLOPIS-75), Laj et al. (2004) determined the 14 C production rate as a function of past changes in geomagnetic field intensity using the conversion of Masarik and Beer (1999).Hain et al. (2014) have further processed these data with a Monte-Carlo simulation with 1000 members of randomised magnetic field strength records within the GLOPIS-75 uncertainty envelope and a more recent 14 C production rate model by Kovaltsov et al. (2012) with updated cosmic ray spectra.Muscheler et al. (2004) have used a different approach to infer the 14 C production rate across the last 25 ka based on high resolution 10 Be measurements from GRIP and GISP2 ice cores.Assuming that the atmospheric processes of 10 Be (attachment on aerosols and deposition after a mean residence time of 1−2 years in the atmosphere) are understood, the 14 C production rate can be estimated rather directly, because 10 Be and 14 C are produced by similar processes in the atmosphere (Lal and Peters, 1967;Masarik and Beer, 1999).This reconstruction takes into account all probable causes of production rate variations, such as variations in geomagnetic field, solar activity and/or in the interstellar galactic cosmic rays flux (Muscheler et al., 2004).
In their studies, Muscheler et al. (2004) as well as Laj et al. (2004) only provide normalised (around 1) values for the 14 C production rate, because the determination of the absolute values is very uncertain.Only the more recent study by Hain et al.

S14
(2014) provides these absolute values (in atoms/cm −2 s).For this reason, we scaled the two other data sets by factors that yield 14 C production rates throughout the last 25 kaBP that are close to the data by Hain et al. (2014) and moreover, lead to atmospheric ∆ 14 C results close to the data based reconstructions at the beginning of the MI in the DCESS simulations.
An evaluation of good fits for these scaling factors yielded 1.5 • 10 4 atoms • cm 2 s −1 for the data by Muscheler et al. (2004) and 1.35 • 10 4 atoms • cm 2 s −1 for the data by Laj et al. (2004).For the box diffusion model-calculated 14 C production rate (Muscheler et al., 2005, see Fig. 6), the most likely absolute values are generated by multiplication with 1.8•10 4 atoms • cm 2 s −1 (R. Muscheler, personal communication, 2015).In Fig. S8, the three 14 C production rate time series are shown.Generally, the three individual 14 C production rates show similar patterns across the last 25 kaBP but with a couple of differences.The data by Hain et al. (2014) and by Laj et al. (2004) are derived using the same approach and the general features and local extrema therefore show the same timing, although the time series by Hain et al. (2014) seems somewhat smoother with a smaller general trend.The data by Muscheler et al. (2004) shows far stronger and more high-frequent fluctuations than the other two data sets, the general trend, however, is similar.Only between 20 and 23 kaBP do the time series of the different approaches show an opposing behaviour.While the data set by Muscheler et al. (2004) shows a minimum during that period, the other two time series possess maxima here.Between the beginning and the end of the MI, all data sets are exhibit similar structure.However, over these 3 ka, a reduction of 0.11 atmos/cm 2 s can be seen in the data by Hain et al. (2014), 0.19 atmos/cm 2 s in the data by Laj et al. (2004) and even a slight increase of 0.02 atmos/cm 2 s in the data by Muscheler et al. (2004).
agrees well with the data sets, but at around 15 kaBP, the DCESS time series drops while both data sets show an increase.
Hence, the calculated production rate lies below both uncertainty ranges at that period, where also a large offset between observed and modelled ∆ 14 C can be seen in Fig. S9.The box diffusion model calculated 14 C production rate by Muscheler et al. (2005) also shows a strong decrease during this period where the observational data rises.After this peak, the observations decrease again and agree well with the DCESS calculations thereafter, particularly the time series by Muscheler et al. (2004).
Across the Holocene, the DCESS calculated 14 C production rate is rather at the low end of the uncertainty ranges of the proxy data sets.
The model-calculated 14 C production rate lies within the proxy data-based uncertainty limits during most of the simulation but shows large deviations at the end of the MI.This discrepancy can also be observed in the box diffusion model-calculated time series by Muscheler et al. (2005) (which assumes a constant carbon cycle).Hence, this approach can not answer the question if it is the lack of an important process in the model, which could possibly also account for the remaining increase of atmospheric pCO 2 , or errors in the production rate data sets that lead to these deviations.

Figure
Figure S1.a) Normalised GSD zone area fraction as function of global mean temperature deviation from PI climate conditions.The yellow bar marks the region between LGM and PI.b) Latitude dependency of albedo for three different deviations from the global mean temperature (−4 • C, 0 • C and 4 • C).Note that poleward of the snow line the albedo is 0.62 (albedo of snow/ice covered area).

Figure S3 .
Figure S3.Transition functions for iron-limitation factor f F e−lim (left panel) and for additional radiative forcing ADust due to dust (right panel) variations between LGM and PI.

Figure S4 .
Figure S4.DCESS model-calculated "snowline" (red, from the All_TF simulation in Fig. 5) and prescribed line of NH ice sheet expansion (blue) converted as described in the text from data by Shakun et al. (2012).The current minimum of the two lines (thick) is used for the calculation of the albedo and the vegetation extent.The yellow shading marks the period of the MI.
S9)and imposed by computing D p,tot = max(D p , D p,min ).We used R C = 1.0 to keep the parameter at standard PI conditions in the shallow ocean for simplicity.The function describes a sharp reduction of the high latitude ocean vertical exchange at an ocean depth depending on p a (see below).This approach isolates deep and intermediate waters in the model high latitude ocean from upper waters and thereby from the atmosphere.This leads to carbon enrichment and isotopic depletion of these deep and intermediate waters.As a consequence of this reduction in ocean vertical exchange, there is a 20 ppm drawdown of pCO 2 in the PI atmosphere and a 0.5 • C reduction in global mean temperature.Furthermore, atmospheric δ 13 C increases by around 0.1‰ and ∆ 14 C by around 40‰.For our efforts to generate an LGM climate state, which constrains all the mentioned variables to proxy data records concurrently, we used the parameter p a to vary the transition depth and thereby the oceanic volume of isolated deep water.This yielded p a = 5 for our representation of LGM climate and a vertical exchange transition at around 1800 m depth.Fig.S5shows this vertical exchange profile as a function of ocean depth.

Figure S5 .
Figure S5.Depth profile of the vertical diffusivity in the high latitude model ocean for generating isolated intermediate and deep waters (blue line, bounded by red lines).The left red line denotes the background vertical mixing from the low latitude model ocean, which serves as a lower boundary and the right red line is the standard PI diffusion and describes the upper boundary for the LGM profile.

Figure S7 .
Figure S7.Atmospheric values (a: pCO2, b: global mean temperature, c: δ 13 C and d: ∆ 14 C) for transient DCESS simulations from 20 to 10 kaBP and data-based reconstructions.Black: Data-based reconstructions; pCO2 by Lüthi et al. (2008), temperatures by Shakun et al. (2012), δ 13 C by Schmitt et al. (2012) and ∆ 14 C by Reimer et al. (2013); cyan: DCESS simulation only with prescribed ice sheet extent and minor greenhouse gas forcing (Min_RF); magenta: as cyan, additionally with restored high latitude vertical exchange (Add_vdiff); blue: as magenta, additionally with prescribed temperature-dependent dust concentration (for radiative forcing and iron-limitation) (Add_dust); red:as blue, additionally with sea level-dependent phosphate and salinity variation (All_TF).The yellow shading indicates the MI from 17.5 to 14.5 kaBP.Note that at many stages the blue Add_dust line is covered by the red All_TF line.
Schilt et al. (2010)011, forers adjusted to achieve desired global mean temperature and climate sensitivity.In order to generate a climate sensitivity of 2.5 • C for a doubling of atmospheric CO 2 concentration (as proposed bySchmittner et al., 2011, forLGM simulations) and to achieve a PI global mean temperature of 15 • C, we choice a value for the longwave radiation temperature sensitivity, B, of 2.21 W m −2 K −1 and a value for the mean incoming solar radiation at 0 • C, Ao, to be 207.16Wm−2 .In the above radiative forcing relation, A = A o −A t , whereby A t is the sum of the radiative effects of the deviations of carbon dioxide, methane and nitrous oxide from PI levels.Data reconstructions presented bySchilt et al. (2010)show atmospheric partial pressures of around pCH 4 = 380 ppm and pN 2 O= 200 ppm during the LGM.For the radiation calculations, we adopt these values and prescribe the time series in the transient simulations.

Table S1 .
Parameter values used to generate the DCESS model LGM climate state.
(Shaffer et al., 2008)014, and citations therein)in Fig.S6d.These data show quite good agreement with the model simulations considering the sparse measurements and the coarse model resolution.The high Figure S6.Profiles of the DCESS model ocean at LGM conditions after 80 ka of integration.Blue: high latitude model ocean profiles; red:low-mid latitude ocean profiles; black: observed present day low-mid latitude ocean profiles(Shaffer et al., 2008).

Table S2 .
Atmospheric and oceanic model variables for PI and LGM climate conditions and LGM proxy data for comparison.The range of the proxy data values express their variability between 18 and 25 kaBP.