Climate of the Past Simulation of the last glacial cycle with a coupled climate icesheet model of intermediate complexity

A new version of the Earth system model of intermediate complexity, CLIMBER-2, which includes the three-dimensional polythermal ice-sheet model SICOPOLIS, is used to simulate the last glacial cycle forced by variations of the Earth’s orbital parameters and atmospheric concentration of major greenhouse gases. The climate and icesheet components of the model are coupled bi-directionally through a physically-based surface energy and mass balance interface. The model accounts for the time-dependent effect of aeolian dust on planetary and snow albedo. The model successfully simulates the temporal and spatial dynamics of the major Northern Hemisphere (NH) ice sheets, including rapid glacial inception and strong asymmetry between the ice-sheet growth phase and glacial termination. Spatial extent and elevation of the ice sheets during the last glacial maximum agree reasonably well with palaeoclimate reconstructions. A suite of sensitivity experiments demonstrates that simulated ice-sheet evolution during the last glacial cycle is very sensitive to some parameters of the surface energy and mass-balance interface and dust module. The possibility of a considerable acceleration of the climate ice-sheet model is discussed.


Introduction
Simulation and understanding of glacial cycles, which dominated climate variability over the past several million years, still remain a major scientific challenge.Although a large body of palaeoclimate evidences supports the hypothesis of the Earth's orbital variations as pacemaker for glacial cycles, a number of scientific questions remain unsolved.Among Correspondence to: A. Ganopolski (andrey@pik-potsdam.de) them is the nature of the dominant 100 kyr cyclicity over the past 1 million years, which is not present in the frequency spectrum of "Milankovitch forcing" (summer insolation in the high latitudes of the Northern Hemisphere, NH).Numerous hypotheses addressing the nature of glacial cycles have been proposed in recent decades, but testing these hypotheses with comprehensive Earth system models represents a challenge due to very long orbital time scales.So-called state-ofthe-art comprehensive coupled models of the general circulation of the atmosphere and the ocean (GCMs) are still too computationally expensive to simulate even a single glacial cycle.First simulations of glacial cycles were performed in 80th and 90th with rather simple climate-cryosphere models such as zonally averaged or two-dimensional energy-balance climate models coupled to simplified ice-sheet models (Pollard, 1982;Deblonde et al., 1992;Gallée et al. 1991, Berger et al., 1999).These experiments demonstrated that, when forced by variations of the Earth's orbital parameters, simulated ice sheets experience large variations on all major orbital frequencies (of precessional angle, obliquity and eccentricity) with a clearly asymmetric temporal dynamics consistent with palaeoclimate data.It was also shown by Loutre and Berger (2000) that glacial-interglacial variations in CO 2 concentration alone cannot drive the glacial cycles.Only when these variations are added to the orbital forcing does agreement between simulated and reconstructed glacial cycles improve considerably.This suggests that climate-carbon cycle feedback plays an important role in shaping glacial cycles.More recently, a number of simulations of glacial cycles were performed using three-dimensional ice-sheet models forced by patterns of glacial climate changes (e.g.Zweck and Huybrechts, 2003;Charbit et al. 2007) or individual climate forcing components (Abbe-Ouchi et al., 2007) obtained separately from the GCM experiments and scaled usage of palaeoclimate data or output of the ice-sheet models.Petoukhov et al. (2000).
Abbreviations are explained in the Table 1.
In recent years, a new class of models, the so-called Earth system model of intermediate complexity (EMICs; Claussen et al., 2002) became available for the study of glacial cycles.These models are far less computationally demanding than coupled GCMs but incorporate substantially more physical processes than simple models.In particular, several EMICs simulate the most recent glacial inception (around 120-115 kyr BP) with reasonable agreement with palaeodata, at least, in respect of the total ice volume (Wang and Mysak, 2002;Calov et al. 2005).A stability analysis of the climate-cryosphere system performed by Calov and Ganopolski (2005) has shown that the glacial inception represents a bifurcation transition in the climate-cryosphere system and that the climate system possesses multiple equilibria within the range of Milankovitch forcing, thus supporting the nonlinear paradigm proposed by Paillard (1998Paillard ( , 2001)).Recently, the entire last glacial cycle was simulated with the same climate model used in this study, but coupled to another ice sheet model and using a different coupling technique (Bonelli et al., 2009).
The purpose of our paper is to present the results of our simulations of the last glacial cycle with the CLIMBER-2 model, where the climate and ice-sheet components are coupled bi-directionally using a physically-based surface energy and mass balance interface.This work represents a continuation of our previous work (Calov et al., 2005), but with an improved version of the model.It is important to note that in this work, similar to previous modelling studies, we prescribed time-dependent radiative forcing of the major greenhouse gases, in addition to orbital forcing.Such a modelling experiment represents an important test for the climate-cryosphere model, but cannot be considered as a decisive test for the Milankovitch theory, since, in the real world, changes in the greenhouse gas concentration are not an external forcing, but rather an important internal feedback in the Earth system.

Model description
The model used for this study is the newest version of CLIMBER-2.It is similar to that described by Calov et al. (2005) (hereafter C05).Since the model and its performance has been described in a number of previous publications, here we give only a short, general description and will discuss in more detail only the essential improvements compared to the version used in C05.
The CLIMBER-2 model includes six components of the Earth system: atmosphere, ocean, sea ice, land surface, terrestrial vegetation and ice sheets.The first five components are represented by coarse-resolution modules of intermediate complexity and were described in detail by Petoukhov et al. (2000) and Brovkin et al. (2002).The ice-sheet component is represented by the relatively high-resolution three-dimensional polythermal ice-sheet model SICOPOLIS (Greve, 1997).The ice sheet model is applied only to the Northern Hemisphere (NH).Unlike the majority of previous studies of glacial cycles, which have employed a simple empirical-based parameterisation (the so-called Positive Degree Day method) for the simulation of the mass balance of the ice sheets, in our model the coupling between climate and ice-sheet components is provided via the high-resolution physically-based Surface Energy and Mass-balance Interface (SEMI) described in C05 (see Fig. 1).Using daily fields of temperature, humidity, cloudiness, precipitation, wind, longwave and shortwave radiation fluxes simulated by the coarse resolution climate component, SEMI simulates annual surface mass balance and temperature of the ice sheets on the 1.5 • ×0.75 • grid of the SICOPOLIS model.Compared to the commonly used Positive Degree Day method, the physically based SEMI has the important advantage that it is equally applicable to any region of the Earth and to any climate conditions -irrespective of how different they are from the present one.Moreover, it explicitly accounts for the variations of the shortwave radiation associated with orbital forcing and it directly accounts for the effect of dust deposition on snow albedo.In C05, we have shown that this effect is important for the mass balance of the ice sheets.In the current study, similar to C05, we applied the dust deposition rate based on results of simulations for the present-day and LGM (Last Glacial Maximum) conditions by Mahowald et al. (1999).In addition, we implement a parameterisation for the dust deposition originating from glacial erosion in the model.This is, according to Mahowald et al. (2006a), an important additional source of dust in the vicinity of the ice sheets.We also explicitly included the direct radiative forcing of atmospheric dust in the same way as Schneider et al. (2006).As in C05, the cryosphere component affects climate via changes in the fraction of grid cell covered by ice, which affect surface albedo, elevation and land fraction for each grid cell of the coarse resolution climate component of the CLIMBER-2 model.Note that although the fractions of grid cells covered by the ice sheets change only in the NH; changes in elevation and land fractions caused by sea-level variations were applied to all model grid cells.In the present work, we additionally accounted for the freshwater flux into the ocean originating from the ice sheets with explicit coupling of the ice-sheet module to the ocean module.In the following sections, we will describe the most important model improvements compared to the model version used in C05.

Freshwater flux from the ice sheets
Changes in freshwater flux into the ocean caused by growth and decay of the ice sheets primarily affect climate via changes of the Atlantic meridional overturning circulation (AMOC), which is associated with meridional oceanic heat transport and sea ice cover.In particular, a reduction of the freshwater into the ocean contributes to a strengthening of the AMOC during rapid growth of the ice sheets, while during episodes of massive ice surges into the ocean (Heinrich events) and glacial terminations, an increased freshwater flux into the ocean leads to a weakening, or even a complete cessation of the AMOC, which affects climate in the North Atlantic realm considerably.Surface melting of ice sheets and iceberg calving into the ocean are simulated by the ice-sheet module.Surface ice-sheet melting is now added to the river runoff simulated by the climate component of the model, and calving is treated as surface freshwater flux into the nearest oceanic grid cell.During glacial time, the river routing differed from the present configuration and was evolving in time.Hence, in principle, the river routing scheme should be updated regularly taking into account changes in the icesheet distribution and land elevation.However, for simplicity we use a constant glacial river routine scheme in this study, which due to the coarse CLIMBER-2 spatial resolution differs from the modern one only in few grid cells over the North American and European sectors where modern northward direction of runoff was changed to the eastward and westward directions.In our simulations, we did not found large sensitivity of the simulated glacial cycle to the choice of the routing scheme.The effect of temporal accumulation of meltwater in periglacial lakes has not been taken into consideration and all meltwater from the ice sheets is released into the ocean immediately.Note, that both meltwater and iceberg calving are treated as surface freshwater flux and, hence, the potential effect of penetration of meltwater into the deep ocean with sediment-laden flow (e.g.Roche et al., 2007) is not accounted for.

Radiative forcing of atmospheric dust
Palaeoclimate data indicate that the amount of dust in the atmosphere during glacial times was considerably higher than during interglacials.Modelling studies suggest that enhanced atmospheric dust content during glacial time represented an important additional radiative forcing component (e.g.Claquin et al., 2003).This forcing was extremely inhomogeneous spatially and varied seasonally.Recent estimates with GCMs indicate that the globally averaged radiative forcing of atmospheric dust during the LGM has an order of magnitude of −1 W/m 2 , and according to simulations by Mahowald et al. (2006b) andSchneider et al. (2006), it caused an additional global cooling effect of about 0.5 to 1 • C with a stronger cooling in the NH.Note that this is only radiative (shortwave and longwave) forcing of the dust in the atmosphere.The impact of dust on snow albedo was treated separately (see below).Since the simulations of atmospheric dust radiative forcing are only available so far for the LGM time slice, a parameterisation is needed for the evolution of radiative forcing of the atmospheric dust in the transient experiments.We applied a similar procedure to prescribe the temporal development of radiative dust forcing, as that which was used in C05 for the dust deposition rate.Namely, we assume that the additional (compared to the present one) radiative forcing of dust R d is proportional to the global ice volume as where R d LGM is the radiative forcing of dust at the LGM aggregated over the coarse climate-module grid cells.This relation was used by Schneider et al. (2006).V is the simulated NH ice volume and V LGM =100 msl (metres of sea-level equivalent) is an approximate estimate for the NH ice volume at the LGM.Obviously, this parameterisation is rather crude.First, the assumption that radiative forcing of dust is directly proportional to the global ice volume is somewhat arbitrary and is not held for any location on Earth.Second, the radiative forcing of dust depends not only on the dust loading in www.clim-past.net/6/229/2010/Clim.Past, 6, 229-244, 2010 the atmosphere but also on surface albedo, which changed considerably during the glacial cycle.These problems can only be avoided by the incorporation of the dust cycle directly into the CLIMBER-2 model, work which is now in progress.

Dust deposition rate
It has already been shown by Warren and Wiscombe (1980) that even a tiny amount of dust (relative concentrations of about 100 parts per million of weight, ppmw) mixed with snow decreases the albedo of fresh snow by 0.1.For old snow, this reduction can even reach 0.3.In their initial growth phase, the major NH ice sheets occurred in areas with rather small dust deposition rates.But during the latter stages of the glacial cycle, the ice sheets spread into areas where dust deposition rates, especially under glacial conditions, became sufficient to appreciably affect the albedo of snow.As was shown in C05 and confirmed using a more sophisticated model by Krinner et al. (2006), accounting for the effect of dust on snow albedo prevents the growth of the ice sheets in areas with high dust deposition rates, such as Eastern Siberia.
Probably even more importantly, the impact of dust on snow albedo is related to the dust produced due to glacial erosion.According to empirical data (Mahowald et al., 2006a), the dust deposition rate associated with this mechanism was as high as 10 to 100 g/m 2 /yr just south of the North American and northern European ice sheets.If a typical annual precipitation of about 500 mm/yr is assumed in these areas, such dust deposition rates would result in average concentration of dust in snow of 20 to 200 ppmw, assuming that dust is uniformly mixed with snow over the whole year.In reality, most of the dust associated with glacial erosion is deposited during summer; hence, the concentration of dust in the upper snow layer during snowmelt is expected to be much higher than it would be in case of uniform mixing over the year.Moreover, when snow melts, only a fraction of dust is removed by meltwater and therefore the concentration of dust in snow increases with time.Although the accurate modelling of all of these processes is problematic, related uncertainties are not crucial, because a saturation effect occurs for dust concentration in snow of more than 1000 ppmw and the albedo of snow reaches a value comparable to that for dirty ice.
In C05, we implemented a simple technique similar to Eq. ( 1), which produces a spatial distribution of dust deposition rates at any location using results of present-day and LGM simulations from Mahowald et al. (1999), proportionally weighted to the simulated global ice volume.Such an approach, however, is not justified for the glaciogenic dust, because the sources of glaciogenic dust (and therefore its deposition) are strongly related to the extent of the ice sheets and their internal dynamics.Therefore, it is not justified to apply glaciogenic dust simulated for the LGM ice sheets to other periods with very different configurations of the ice sheets.To resolve this problem, we introduced, in our model, a rather simple parameterisation of dust deposition produced from the glaciogenic sources.This parameterisation is based on the assumption that the emission of glaciogenic dust is proportional to the delivery of glacial sediments to the edge of an ice sheet (see Appendix A for details).Most of the glaciogenic dust originates from the southern flanks of the ice sheets and this source is significant only for mature ice sheets, which reach well into areas covered by thick terrestrial sediments.Parameters of the glaciogenic dust module were tuned to reproduce the reconstructed rate of dust deposition during the LGM.Simulated glaciogenic dust was added to the dust deposition rate taken from Mahowald et al. (1999).Note that Mahowald et al. (1999) did not account for glaciogenic dust sources.

Sliding parameterisation
As shown by modelling studies by Calov et al. (2002) and Calov and Ganopolski (2005), fast sliding processes over areas covered by deformable sediments play an important role in ice-sheet dynamics, both on millennial and longer time scales.In the current work, we use three types of surfaces underlying the ice sheets: rocks, deformable terrestrial sediments and marine sediments.The first two types correspond to present-day land and were distinguished by prescribing a minimum sediment layer thickness in the data source by Laske and Masters (1997), while all areas which lay below present-day sea level were considered to be covered by marine sediments.Basal sliding appears only if the base of an ice sheet is at the pressure melting point, in which case, the sediment type of underlying surface becomes important.Two different sliding laws were applied for rock (Calov and Hutter, 1996) and sediment (C05), while the difference between terrestrial and marine sediment was represented by the value of the sliding parameter, which is an order of magnitude higher for marine sediment.

Temperature correction for North America
As shown by Petoukhov et al. (2000), CLIMBER-2 has reasonable skill in simulating spatial and seasonal variability of different climatological fields relevant for ice sheets.In particular, biases in simulated present-day summer temperatures, the key climatological factor determining ablation of ice sheets, are within 1-2 • C and annual precipitation is typically within 30% of observed values.However, with its coarse spatial grid the CLIMBER-2 model cannot resolve some regional features of climate characteristics which are important for the evolution of ice sheets.This is especially true for North America, which is represented in the model by a single grid column.At the same time, due to orographic and geographic peculiarities, there is a strong zonal temperature gradient over North America, which cannot be resolved by the model, but is essential for the development of the Northern American glaciation.Preliminary experiments demonstrated that the lack of this zonal gradient causes a problem for the realistic simulations of the North American ice sheets (see also the discussion in Sect.6).Thereby, we implemented a sub-grid correction for the atmospheric temperature in the North American sector, which has the shape of a dipole and closely resembles the deviations of the observed summer air temperature (corrected for elevation effect) from its zonal mean over North America (Fig. 2).The maximum of this applied temperature correction is 3 • C.Although this temperature structure is more representative for the summer season, we kept the temperature correction constant in time, since it is the summer temperature that is most important for the mass balance of an ice sheet.The temperature correction was added to the surface air temperature interpolated from the coarse climate module grid.This temperature correction directly affects the surface sensible heat flux and, in addition, the downward longwave radiation by assuming a linear relationship between temperature and downward longwave radiation.

Experimental setup
Similar to C05, the model was forced by variations in orbital parameters computed following Berger (1978) and greenhouse gas (GHG) concentrations derived from the Vostok ice core.Unlike C05, we now account for two other important GHGs: CH 4 and N 2 O. Since the radiative scheme of the CLIMBER-2 model does not include CH 4 and N 2 O, the radiative effect of these gases was incorporated via the socalled equivalent CO 2 concentration, which is determined as the CO 2 concentration which has the same radiative forcing as the combined radiative forcing of all major greenhouse gases (see Appendix B for details).All model runs started from the equilibrium state corresponding to the orbital configuration and the concentration of GHGs at 126 kyr BP.The model then was run fully interactively and synchronously for 126 000 years until present-day, thus covering the whole last glacial cycle.
We will start the discussion of modelling results with a so-called Baseline Experiment (BE).This experiment represents a "suboptimal" subjective tuning of the model parameters to achieve the best agreement between modelling results and palaeoclimate data.Obviously, even with a model of intermediate complexity it is not possible to test all possible combinations of important model parameters which can be considered as free (tunable) parameters.In fact, the BE was selected from hundred model simulations of the last glacial cycle with different combinations of key model parameters.Note, that we consider "tunable" parameters only for the ice-sheet model and the SEMI interface, while the utilized climate component of CLIMBER-2 is the same in previous studies, such as those used by C05.In the next section, we will discuss the results of a set of sensitivity experiments, which show that our modelling results are rather sensitive to the choice of the model parameters.Therefore, the simulation of a realistic glacial cycle would represent a too ambitious task if the model was too computationally expensive to perform a large set of experiments.The selection of the best fit to the palaeodata is based on several criteria, among which are an accurate simulation of the total ice-volume variations over the whole glacial cycle, and an agreement between simulated and reconstructed ice sheets during LGM.It is important to note that some systematic biases exist in all model simulations; we were not able to eliminate them with any combination of the model parameters.several palaeoclimate reconstructions.When we compared the model results with sea-level change reconstructions, the problem was that our model does not account for the Southern Hemisphere ice sheets, which, in accordance with different estimates, contributed an additional of ca.20 m to the LGM sea-level drop.This contribution, estimated by Huybrechts (2002), was added to the modelled NH ice volume and plotted in the figure together with an estimate of the range of sea-level variations based on recent coral data compilations (Lambeck et al., 2002;Thompson and Goldstein, 2006; supplementary data therein) and an independent estimate of sea-level change based on δ 18 O c and deep ocean temperature (Waelbroeck et al., 2002).As seen in the figure, the model simulates all major aspect of the ice-volume variations during the whole glacial cycle in reasonable agreement with palaeoclimate estimates.As compared to most previous simulations of the last glacial cycle based on different modelling approaches (e.g.Tarasov and Peltier, 1997;Marshall and Clark, 2002;Zweck and Huybrechts, 2004;Abe-Ouchi, 2007;Charbit et al. 2007;Bonelli et al., 2009), our model simulate much more pronounced (and, in generally, more realistic) variations in the ice volume on the precessional time scale (ca.20 kyr).This is likely to be explained by fact that the previous modelling studies employed the positive degree day approach, which only indirectly accounts for insolation changes.In our model, the mass balance of the ice sheets is computed using a surface energy balance approach, which is more sensitive to variations in insolation.
Similarly to C05, initial buildup of the ice sheets during the last glacial inception starts soon after 120 kyr BP in Scandinavia, Arctic Archipelago and Labrador (Fig. 4a).Already around 112 kyr BP (Fig. 4c), the ice volume reaches 50 meters in sea-level equivalent (msl).During warm stages MIS 5c (Fig. 4d) and 5a most of the European ice sheets disappeared, whilst a relatively large ice sheet survived in the eastern part of North America.The magnitude of ice volume variations during MIS 5 as well as the steep rise of the ice sheets during MIS 4 are correctly simulated by the model.During most of MIS 3 the model simulate a relatively stable ice volume and the rise of the ice volume resume around 30 kyr BP, in agreement with palaeoclimate data.The maximum NH ice volume of ca.90 msl was reached at around 20 kyr BP, i.e. one to several thousand years later than indicated by palaeoclimate reconstructions.At the LGM, the North American ice sheets contribute about 70 m to global sea-level drop, while the northern European ice sheet contributes slightly below 20 m.The glaciation in the eastern Siberia, which developed shortly before the LGM, contributes only several meters.The spatial extent of simulated North American ice sheets (Fig. 5a) is close to the palaeoclimate reconstructions.In particular, east-western gradient in the position of the southern margin of North American ice complex is captured well.Compared to ICE-5G reconstruction (Fig. 5b) by Peltier (2004) the elongated ice dome in the middle of the simulated Laurentide ice sheet is too thin, although our model captures the major topology of this structure.Compared to the palaeoclimate reconstructions, simulated glaciation of Alaska is excessive.The extent and ice thickness of the Fennoscandian ice sheet, the south-western part of the European ice complex, are in good agreement with the ICE-5G reconstruction, while the north-eastern part of this ice complex does not reach far enough to the east, i.e. there is only a limited ice cover over Kara sea in the simulations.
Deglaciation of the NH began soon after 20 kyr BP and accelerates significantly after 16 kyr BP.Between 16 kyr BP and 14 kyr BP the averaged rate of ice-sheet melting reaches 15 msl per thousand years, which is comparable, though still smaller, than observed rate of sea-level rise during the MWP-1A event.During termination, the large portion of the European ice sheets in the response to rising summer insolation and GHG concentration disappear before 13 kyr BP (Fig. 4f) and the most of the remaining European ice sheet melted prior to 10 kyr BP (Fig. 4g), which is ca. 1 kyr too early compared to palaeoclimate reconstructions.The initial response of the Laurentide ice sheet is slower and only around 15 kyr BP does its retreat accelerate significantly.At around 13 ky BP (Fig. 4f), the Laurentide and Cordilleran ice sheets separated and the remains of Laurentide ice sheet on Baffin Island (Fig. 4h) completely disappear soon after 8 kyr BP.
Figure 6 shows the evolution of the volume and area of the North American and European ice sheets during the last glacial cycle.The volume of the European ice sheet experiences stronger variability with the precessional period and it almost completely vanishes during MIS 5c and MIS 5a (Fig. 6a).The North American ice sheets are also affected by precessional variations but it survived periods of high summer insolation and has a general tendency of growth through the whole glacial cycle.The area of the NH ice sheets closely follows NH ice volume except for the periods of rapid glacial advances (Fig. 6b).For example, during glacial inception the ice-sheet area reaches 1/3 of its maximum (LGM) value already at 115 kyr BP, while the total ice volume at that time is only 1/10 of the LGM value.This is explained by an abrupt appearance of large snow fields caused by strong positive snow-albedo feedback, as discussed in C05. Figure 6c   melting point (temperate basal area).The relative fraction of the temperate basal area of the North American ice sheets, in spite of large fluctuations, has a tendency to increase towards the end of the glacial cycle, which is qualitatively similar to results by Marshall and Clark (2002).At the onset of glacial termination, almost one third of the ice-sheet area is temperate at the base; this contributes to rapid deglaciation due to fast sliding of the ice at the southern margins of the ice sheets.Figure 7 shows the temporal evolution of three major components of the mass balance of the NH ice sheets: accumulation, ablation and ice calving.The total accumulation relatively closely follows the ice-sheet area.Although the area of the NH ice sheets varies by one order of magnitude during the glacial cycle and precipitation decreases considerably during glacial maximum, the averaged accumulation rate over the ice sheets remains surprisingly stable, about 1 mm/day.In turn, surface ablation of the ice sheets closely follow Milankovitch forcing, even though variations in GHGs also play an important role in shaping the mass balance of the ice sheets.Interesting enough, ablation actually leads Milankovitch forcing by several thousand years, which is explained by the fact that ablation also depends on the area of the ice sheets, which shrinks with increasing insolation.This fact clearly demonstrates the non-linear nature of the ice sheets response to the orbital forcing and implies potential problems of the dating technique based on orbital tuning since the latter is based on an assumption of linear response.The peak in the ablation rate of about 0.3 Sv during deglaciation at around 15 kyr BP is associated with the rapid resumption of the Atlantic thermohaline circulation (see Fig. 8), which causes pronounced warming over the North Atlantic realm.
Simulated ice sheet calving is rather noisy, which is partly attributed to the numerics of the ice sheet model.For this reason, the total calving rate shown in Fig. 7 is smoothed by a 10 year running-mean window.Still, it reveals strong temporal variability; the most pronounced of which is related to the large-scale instability of the Laurentide ice sheet via a mechanism described by Calov et al. (2002) and which resembles in many respects real Heinrich events.During these simulated ice surges, the calving rate exceeds 0.1 Sv.Although ablation remains the major negative contribution to the mass balance of the ice sheets over the whole glacial cycle, ice calving plays an important role in the millennial scale variability of AMOC.

Climate evolution
As shown with the CLIMBER-2 model in previous studies (Ganopolski, 2003;Schneider et al., 2006), approximately half of the global cooling (compared to present) during the last glacial maximum is explained by the presence of the large ice sheets in the NH and the rest by a lowering of the concentrations of the major GHGs (primarily CO 2 ), an increase in atmospheric dust content and a shrinking of vegetation cover.Since all of these factors are closely related and have a similar temporal dynamics, it is not surprising that the temporal evolution of the globally averaged annual surface air temperature during the glacial cycle follows the same pattern as the ice volume and CO 2 .On regional scale, however, the temperature variations are more complicated and diverse.The northern North Atlantic temperature, which we consider as a proxy for Greenland temperature, apart from strong variations on orbital time scales, experienced numerous abrupt changes with a magnitude of up to 10 • C during a large portion of the glacial cycle (Fig. 8).In many respects, these tooth-shaped fluctuations resemble Dansgaard-Oeschger events recorded in numerous locations over the NH.Such fluctuations are attributed to rapid reorganisations of the Atlantic Meridional Overturning Circulation (AMOC).The largest perturbation of the AMOC, when circulation almost stalled, coincides with the quasi-regular ice surges from the Laurentide ice sheet discussed by Calov et al. (2002).Due to the random nature of millennial-scale variability, one cannot expect one-to-one correspondence between simulated and reconstructed temperatures.For example, the model equivalent of Bølling-Allerød occurs thousand years early than in reality.However, the overall qualitative agreement between model and data is quite instructive.
Antarctic temperature follows the concentration of GHGs more closely, but is also strongly affected by orbital variations (Ganopolski and Roche, 2009).It also exhibits pronounced millennial scale variability, which is the counterpart of that in the Northern Atlantic and which results from the seesaw mechanism (Crowley, 1992).Similar to Ganopolski and Rahmstorf (2001), the strongest warming events in Antarctica are associated with the largest disturbances of the AMOC caused by the model's equivalent of Heinrich events (Fig. 8).As discussed above and illustrated by simulations described below, dust deposition plays an important role in the simulation of the glacial cycle with our model.As shown in Fig. 9a, d, dust deposition over the NH extra-tropics was significantly higher at the LGM as compared to the interglacial state.While outside of the ice sheets, the dust deposition is prescribed by the weighted mean of modern and LGM dust deposition rates simulated in Mahowald (1999) (Fig. 9b), the glaciogenic sources of dust play the dominant role in the vicinity of the ice sheets (Fig. 9c).At the LGM, the simulated values of dust deposition rates just south of both major NH ice complexes exceed 50 g/(m 2 yr), which is compatible with empirical data presented in Mahowald et al. (2006a).As shown in Fig. 10, the dust deposition rate south of the European and North American ice sheets experienced variations with precessional period.However, in Europe the peaks in dust deposition are of similar magnitude.The dust deposition rate south of Northern American ice sheet is much higher during LGM, which is explained by the fact that only prior to glacial termination a large portion of the Laurentide ice sheet spreads into the area covered by terrestrial sediments.Simulated dust deposition rate reduced snow albedo significantly and facilitated the retreat of the ice sheets during deglaciation.

Sensitivity experiments
The ice sheet model and the ice sheet-climate interface contain a number of parameters which are not derived from first principles.They can be considered as "tunable" parameters.As stated above, the BE was subjectively selected from a large suite of experiments as the best fit to empirical data.Below we will discuss results of a number of additional experiments illustrating the sensitivity of simulated glacial cycle to several model parameters.These results show that the model is rather sensitive to a number of poorly constrained parameters and parameterisations, demonstrating the challenges to realistic simulations of glacial cycles with a comprehensive Earth system model.

Sensitivity to the ice-sheet model parameters
The ice-sheet model SICOPOLIS used for our palaeoclimate simulations has two major "free" parameters controlling icesheet dynamics.One is the so-called enhancement factor, which appears as a coefficient in the flow law of ice to account for dust impurities in the ice (Paterson, 1994).A higher enhancement factor increases the deformation velocities of the ice.A constant value of 3 for the enhancement factor is often used for the Greenland ice sheet to broadly capture the contribution of dust load deposited on the ice sheet during glacial time and transported deeper into the ice sheet by ice advection.Considering the dust deposition data by Mahowald et al. (1999Mahowald et al. ( , 2006a)), it becomes obvious that the enhancement factor could be much higher for the NH ice sheets, in particular, during glacial times.The other free parameter determines the bottom sliding over the areas where the base of the ice sheet is temperate.
Experiments where the enhancement factor ranges between 3 and 12 show a rather weak sensitivity of simulated glacial cycles to this parameter (Fig. 11a).This is explained by the relatively minor contribution of the deformation velocities, among other factors controlled by the enhancement factor, to the ice-sheet movement, as compared to basal sliding.Another reason for the low sensitivity of the ice dynamics to the enhancement factor lies in the ice-flow temperature coupling: if modelled ice velocities increase due to a higher enhancement factor, the ice sheet becomes cooler overall, because colder and less deformable ice is advected downward, which, in turn, counteracts the velocity increase due to the enhancement factor.Doubling of the bedrock relaxation time scale (not shown) has a similarly small effect on global ice volume.In contrast, the sensitivity to the parameterisation of sliding processes is rather strong.The experiment where all continental grid points are treated as rocks differs considerably to the experiment where all continental grid points are treated as terrestrial sediments (Fig. 11a).In the first case, even when basal temperature is at the pressure melting point, the bottom sliding is rather small and the ice sheet reached at the LGM a much larger thickness than in the BE.In addition, due to much smaller rate of dust deposition simulated in this experiment, complete deglaciation of the NH is delayed by ca. 5 kyr compare to BE.In the experiment where all continents are covered by terrestrial sediments, the ice sheets are much more mobile and considerably thinner than in the BE.As a result, during periods of high summer insolation, the ice sheets vanished completely and regrow only during periods when precessional variations of summer insolation are reinforced by low obliquity.Hence, the simulated glacial cycle is rather sensitive to the parameterisation of bottom sliding, a process which is not yet properly understood.(a) Sensitivity to ice-sheet-model parameterisations.The black line corresponds to the enhancement-factor value 12 (the value used in BE is 3), the blue line corresponds to the experiment without terrestrial sediments, the red line to the experiment where all land is covered by terrestrial sediments.(b) Sensitivity to parameterisations of the surface energy and mass balance interface.The black line corresponds to the experiment without parameterisation of elevationdesert effect, the blue line corresponds to the experiment without effect of elevation slope on precipitation, and the red line represents the experiment without the North American temperature dipole correction.(c) Sensitivity to the effect of the glaciogenic dust deposition.The black line represents the experiment without the effect of dust deposition on snow albedo, the blue line with a halved and the red line with doubled amounts of the glaciogenic dust deposition.

Sensitivity to the parameters of the surface energy and mass balance interface
The surface energy and mass balance interface (SEMI) includes a number of parameterisations, which attempt to account for processes that are not explicitly resolved by the coarse-resolution climate component.One such parameterisation is the so-called "elevation-desert effect" on precipitation (Eq.7 in C05).As shown in Fig. 11b, this parameterisation plays an important role.Indeed, in the experiment where this parameterisation was disabled, i.e. precipitation is simply spatially interpolated from the coarse resolution climate grid, the total ice volume is much higher than in the BE www.clim-past.net/6/229/2010/Clim.Past, 6, 229-244, 2010 (40% higher at the LGM).As a result, the NH ice sheets did not disappear completely at the end of the simulation.To the contrary, switching off the parameterisation of slope effect on precipitation (Eq.6 in C05), which leads to a reduction of the precipitation over the marginal parts of the ice sheets, considerably slows down the ice sheet growth.As a result, NH ice sheets vanished completely between 100 kyr BP and 80 kyr BP and reaches only a half of its LGM volume compared to BE.Another attempt to correct unresolved regional climate pattern is the "American temperature dipole" described above.Although, the magnitude of this temperature correction is only 3 • C, which is comparable to the typical temperature biases of state-of-the-art climate models, disabling the temperature correction over North America has a considerably effect on the simulated glacial cycle (Fig. 11b), comparable with the effect of switching off the slope effect on precipitation discussed above.Nonetheless, one can see that using of the "American temperature dipole" parameterisation is not the precondition for our simulation of pronounced ice volume variations at the orbital time scales.

The role of the atmospheric dust
One of the important novelties of the modelling approach used in this study is the explicit treatment of glaciogenic dust.The proposed parameterisation is rather simple and is only weakly constrained by empirical data.To test how sensitive the simulated glacial cycle is to the amount of dust deposition due to glacial erosion, we changed empirical parameter k g in Eq. (A3) by −100%, −50% and +100%, i.e. from zero to a doubling of the glaciogenic dust deposition rate compared to the Baseline Experiment.As shown in Fig. 11c, the simulated glacial cycle is rather sensitive to the amount of glaciogenic dust deposition, especially during the first part of the cycle and glacial termination.In absence of glaciogenic dust, the simulated ice-sheet volume is considerably larger during the entire glacial cycle.During glacial termination, when the deposition of glaciogenic dust reaches its maximum in the BE, in the absence of glaciogenic dust, almost one third of the LGM ice survives the Holocene.In the case of halving and doubling of the standard value of k g , the difference with the BE is not so pronounced, which is explained primarily by the logarithmic dependence of snow albedo under low concentration of dust and the saturation effect under very high concentrations.Hence, at least in our model, accounting for the additional source of dust related to the glacial erosion is crucial for simulating of a complete termination of the glacial cycle, although the simulated ice volume is not too sensitive to the uncertainties in the parameterisation of this dust source.

Acceleration technique
Simulation of even only one glacial cycle with a state-ofthe-art coupled GCM remains computationally extremely demanding.Thereby the possibility of an acceleration of the model runs would be desirable.Since the time scales of the ice sheets are comparable or even longer than periodicity of the orbital forcing, it is not possible to apply any acceleration technique to the ice-sheet model.However, even a relatively high resolution (ca.50 km), modern, three-dimensional thermomechanical ice-sheet model is computationally inexpensive compared to the climate models (about 1% of total CPU time).At the same time, the typical time scales of the atmosphere-ocean system is much shorter compared to orbital time scales and, therefore, it would be justified to accelerate the climate component by artificially stretching the time scales of the external forcing (orbital and GHGs in our case), as proposed by Lorenz and Lohmann (2004).This technique, when applied to the climate ice-sheet model, implies that the ice sheet is simulated in real time, while for the climate component orbital forcing and GHGs concentrations changes N times faster than in reality, where N is the acceleration factor.As a result, the exchange of information between the climate and the ice sheet model components occurs every year for the climate and each Nth year for the ice sheet component.Calov et al. (2009) showed that the climate component can be considerably accelerated without significant loss of accuracy in simulation of the glacial inception.Here, we extend this analysis to the whole glacial cycle.Figure 12a shows the simulated global ice volume in three runs with acceleration factors of 5, 10 and 20 in comparison with the (fully synchronous) BE.With respect to simulated ice volume, the agreement between accelerated and BE remains reasonably good up to an acceleration factor of 20.The same is true for the globally averaged surface air temperature, although, pronounced millennial scale variability is already considerably suppressed for a two-fold acceleration (not shown).Finally, the simulated deep water temperature is considerably affected by acceleration technique (Fig. 12c).Therefore, our experiments indicate that some aspects of glacial variability on orbital time scales can be successfully reproduced, even when using a large acceleration factor but considerable delay is introduced in the deep ocean evolution.The latter problem can be at least partly mitigated by using an additional acceleration scheme for the deep ocean, as proposed by Liu et al. (2004).

Discussion and conclusions
We have presented simulations of the last glacial cycle using the Earth system model of intermediate complexity CLIMBER-2, which incorporates the three-dimensional polythermal ice-sheet model SICOPOLIS bi-directionally coupled to the climate component via a physically-based interface.Compared to our previous work (C05), we additionally introduced radiative forcing of atmospheric dust and the effect of deposition of dust produced by glacial erosion.
Our experiments demonstrate that the CLIMBER-2 model with an appropriate choice of model parameters simulates the major aspects of the last glacial cycle under orbital and greenhouse gases forcing rather realistically.In the simulations, the glacial cycle begins with a relatively abrupt lateral expansion of the North American ice sheets and parallel growth of the smaller northern European ice sheets.During the initial phase of the glacial cycle (MIS 5), the ice sheets experience large variations on precessional time scales.Later on, due to a decrease in the magnitude of the precessional cycle and a stabilising effect of low CO 2 concentration, the ice sheets remain large and grow consistently before reaching their maximum at around 20 kyr BP.The spatial extent of the simulated ice sheets at the LGM agrees reasonably well with palaeoclimate reconstructions.Only the thickness of the central part of Laurentide ice sheet is underestimated and Alaska appear to be glaciated too strongly.Simulated components of the mass balance are strongly affected by the ice volume and orbital forcing.Simulated total ablation is closely related to Milankovitch forcing, but leads it by several thousand years, while the total accumulation closely follows the area covered by the ice sheets.
From about 19 kyr BP, the ice sheets start to retreat with a maximum rate of sea level rise reaching some 15 m per 1000 years around 15 kyr BP.The northern European ice sheets disappeared first, and the North American ice sheets completely disappeared at around 7 kyr BP.Fast sliding processes and the reduction of surface albedo due to deposition of dust play an important role in rapid deglaciation of the NH.Thus our results strongly support the idea about important role of aeolian dust in the termination of glacial cycles proposed earlier by Peltier and Marshall (1995).
During the second part of the glacial cycle, the Laurentide ice sheet experienced large scale internal oscillations with a typical periodicity of about 6000 years, in many respects resembling observed Heinrich events.During the same period of time, glacial AMOC becomes unstable and experiences numerous transitions between different modes of operation, resulting in abrupt climate shifts over the North Atlantic realm resembling observed Dansgaard-Oeschger events.It is important to stress that the model simulates millennial scale climate variability without any explicit external forcing with the same periodicity.Millennial scale climate variability in our model is the results of internal instability of two components of the climate system: ice sheets and AMOC under glacial climate conditions.
Results from a set of sensitivity experiments demonstrate high sensitivity of simulated glacial cycle to the choice of some modelling parameters, and thus indicate the challenge to perform realistic simulations of glacial cycles with the computationally expensive models.At the same time, our results indicate that the computational cost of the simulation of the glacial cycle can be considerably reduced by using of the appropriate acceleration technique.At least as far as the simulated global ice volume is concerned, acceleration with a factor of up to 20 does not significantly disturb the simulation of the glacial cycle.
where r = V /V LGM , V is the NH ice volume (in msl), V LGM =100 msl is the NH ice volume at LGM, D LGM and D MOD are dust deposition rates computed in Mahowald et al. (1999) for LGM and modern conditions respectively.
The deposition of the glaciogenic dust is computed using a rather crude parameterisation based on the assumption that the glaciogenic source of dust is determined by the sediment transport through the margin of ice sheet.Namely, we assume that a certain fraction of the sediments accumulated in the ice-free grid cells neighbouring to the ice covered grid cells becomes airborne and its transport in the atmosphere can be described as a macroturbulent diffusion.It is naturally to assume that the production of glaciogenic dust is proportional to the sediment flux through the margin of the ice sheets.The later is related to the sliding velocity of the ice sheet in the proximity of the ice margin and the advance/retreat velocity of the ice sheet margin.To overcome numerical problems related to fluctuations in simulated values of both velocities, we introduced a new characteristics, the thickness of the "sediment excess" H e , which represents time-integrated convergence of the sediment transport by the ice sheets.The corresponding equation for H e is written as with the additional constraint H e ≥ 0. Here, (u,v) is the sliding velocity of the ice sheets, which is also assumed to be the velocity of the upper layer of sediments, H a is the thickness of the active (water saturated) sediment layer, i.e. the unconsolidated layer of sediments which is moved by the ice sheet.The value of H a is set to 0.1 m for rocks and 1 m for continental area covered by terrestrial sediments.By keeping H a independent of time, we assume that changes in the area covered by terrestrial sediments can be neglected over one glacial cycle.This would not be an appropriate assumption for the simulation of the whole Quaternary.
The term E in Eq. (A2) represents the decay of accumulated "sediment excess" due to water and wind erosion.This term is only applied outside of the ice sheets and is expressed in the form where the characteristic time scale of erosion τ is set to 1000 years.Assuming a fixed fraction of wind erosion in the total erosion of the "sediment excess", the emission of glaciogenic dust per unit of area is where the values of the empirical parameter k g =10 5 g/m 3 was chosen to match reconstructed glaciogenic dust deposition rate given in Mahowald at al. (2006a).Parameter k g has the units of density and its value implies that less than 10% of the sediments transported by the ice sheet become airborne.
Equation (A2) is formally solved over the whole continental area, but positive "sediment excess" is formed only in the grid cells along the margin of the ice sheets and these grid cells become the source of aeolian dust only if they are ice free.
Assuming an universal vertical profile of dust in the atmosphere, isotropic macroturbulent mixing of the aeolian dust with the constant diffusion coefficient K and the constant residence time of the dust in the atmosphere τ a , the balance equation for the total dust load in the atmospheric column Q (in g/m 2 ) can be described by the equation where is the Laplace operator and Q τ = D g is the dust deposition rate.Although Eq. ( A4) is simple, solving it on the ice sheet model grid is computationally expensive.Therefore, we use instead an equivalent but less expensive approach by computing the deposition rate of glaciogenic dust at any location (x o , y o ) by integrating the emission of dust due to glacial erosion S g over the surrounding area D g (x 0 ,y 0 ) = S g (x,y)f (l) dxdy, where l = (x − x 0 ) 2 + (y − y 0 ) 2 and, due to conservation of mass, the universal function f (l) must satisfy the condition f (l) dxdy = 1.
The function f was obtained by solving Eq. (A4) numerically for a singular dust source and using K=2•10 6 m 2 /s and τ a =5 days.The combination of these two parameters gives the length scale L = √ Kτ =10 6 m, which represent the radius of influence of a singular dust source.Since f (l) decrease rapidly with the distance from the source, the size of the integration domain can be chosen with sufficient accuracy on the order of the magnitude of L.
To take into account the effect of elevated ice sheets on the dust transport towards the ice sheet interior, we imposed an additional reduction on the dust deposition over the ice sheet as a function of the distance from the ice sheet margin.The length scale for the dust deposition decay is set to 500 km.

Calculation of the equivalent CO 2 concentration
During glacial cycles variations of CO 2 produce the largest GHG forcing but the role of CH 4 and N 2 O are also not negligible (at LGM their combined radiative forcing is about 0.7 W/m 2 ).The radiative codes of the CLIMBER-2 model does not account for the latter two gases and the simplest way to take them into account is to recalculate CO 2 concentration (the so-called "equivalent CO 2 concentration") which has the radiative forcing equal to the sum of radiative forcings of all three GHGs.However, the CLIMBER-2 model originally was tuned for CO 2 concentration of 280 ppm and introducing of additional radiative forcings of CH 4 and N 2 O will lead to equivalent CO 2 concentration higher than 280 ppm for preindustrial conditions.To avoid necessity of retuning of the model, we calculated equivalent CO 2 concentration in such a way that it is 280 ppm for the preindustrial conditions and gives the right changes in the total radiative forcing for the concentrations of CO 2 , CH 4 and N 2 O different from preindustrial values: where R ECO2 = 5.35ln C e C o is the radiative forcing of the equivalent CO 2 concentration C e , R CO2 = 5.35ln C C o is the radiative forcing of the "real" CO 2 concentration C, and C o =280 ppm is preindustrial CO 2 concentration.Radiative forcings of CH 4 , R CH4 , and N 2 O, R N2O , were computed using corresponding formulas from the IPCC report (Houghton et al., 1996).All radiative forcings are relative to the preindustrial concentrations of corresponding GHGs.By substituting values of individual radiative forcings into Eq.(B1) we calculate the value of equivalent CO 2 concentration, which then was used in model simulations.Note that by definition, C e = C o for preindustrial climate conditions.The concentrations of individual greenhouse gases were averaged over thousand years (in case the resolution of individual records was better than one thousand years).When computing equivalent CO 2 concentration, following Hansen et al. (2005), we took into account, the fact that the resulting radiative forcing of CH 4 is 40% higher than its pure radiative effect due to methane decomposition in the stratosphere and the production of additional water vapour.

Fig. 2 .
Fig. 2. (a) Deviation of observed summer surface air temperature from the zonal average over the CLIMBER-2 sector corresponding to North America.(b) Temperature correction added to the interpolated CLIMBER-2 temperature in the SEMI module.
Figure 3 shows the NH ice-volume variations during the last glacial cycle expressed in terms of global sea-level change relative to present-day values for the BE, in comparison with

Fig. 3 .
Fig. 3. Simulated and reconstructed sea-level evolution during the last glacial cycle, relative to present-day values.The red line shows the sum of the simulated Northern Hemisphere (NH) ice volume and the contribution of the Antarctic Ice Sheet (following Huy-brechts, 2002), expressed in units of sea level change.The grey area represents an interpretation of recent coral data compilations(Lambeck et al., 2002;Thompson and Goldstein, 2006; supplementary data therein).The width of the interval corresponds to two standard deviations.A very small number of well-dated corals for MIS 5d preclude a reliable estimate of the sea level from corals for this period.The solid black line shows the reconstruction byWaelbroeck et al. (2002) based on deep ocean δ 18 O c and temperature records.

Fig. 6 .
Fig. 6.Simulated NH time series.(a) Ice volume, (b) ice area and (c) fraction of the ice sheets temperate basal area.The black line shows the time series of the entire NH, the blue line of the North American and the red line of the European ice sheets.

Fig. 7 .
Fig. 7. Simulated components of the mass balance of the NH ice sheets.(a) Total accumulation (solid blue) and ice-sheet area (dashed black).(b) Total surface melting (blue) and maximum summer insolation at 65 • N (red dashed line).(c) Total ice calving.Units are Sverdrups of water (1 Sv=10 6 m 3 /s).

Fig. 8 .
Fig. 8. Climate forcing and climate evolution.(a) Prescribed radiative forcing of GHGs (green) and orbital forcing (red) represented by changes in maximum summer insolation at 65 • N. (b) Simulated annual anomalies of East Antarctic temperature (red) and anomaly of the deuterium concentration from the EPICA ice core (grey).(c) Simulated annual anomalies of the North Atlantic (60 • -70 • N) temperature (red) and measured anomaly of the δ 18 O concentration from the NGRIP ice core (grey).(d) Simulated rate of Atlantic outflow at 30• S (blue) and total ice calving (green).Note the reverse direction of y-axis for the ice calving.The anomalies in (b) and (c) are relative to present-day.Although the empirical data are arbitrarily scaled to achieve the best fit with the modelled data, the relationships between temperature variations and isotope variations are close to those used for palaeoclimate reconstructions.

Fig. 9 .
Fig. 9. Annual background dust deposition rates for (a) modern conditions and (b) LGM conditions based on Mahowald et al. (1999).(c) Simulated glaciogenic dust deposition rate at the LGM, and (d) total dust deposition rate at the LGM equal to sum of the fields shown in (b) and (c).

Fig. 10 .
Fig. 10.Simulated temporal evolution of annual dust deposition rates in the grid cell south of the Laurentide Ice sheet (blue) and south of Fennoscandian ice sheet (red).The locations of the grid cells are shown in Fig. 9d by rectangles of corresponding colours.

Fig. 11 .
Fig. 11.Simulated NH ice volume in a suite of sensitivity experiments.Grey lines in all figures show the Baseline Experiment.(a)Sensitivity to ice-sheet-model parameterisations.The black line corresponds to the enhancement-factor value 12 (the value used in BE is 3), the blue line corresponds to the experiment without terrestrial sediments, the red line to the experiment where all land is covered by terrestrial sediments.(b) Sensitivity to parameterisations of the surface energy and mass balance interface.The black line corresponds to the experiment without parameterisation of elevationdesert effect, the blue line corresponds to the experiment without effect of elevation slope on precipitation, and the red line represents the experiment without the North American temperature dipole correction.(c) Sensitivity to the effect of the glaciogenic dust deposition.The black line represents the experiment without the effect of dust deposition on snow albedo, the blue line with a halved and the red line with doubled amounts of the glaciogenic dust deposition.

Fig. 12 .
Fig. 12.Effect of acceleration of the climate component on simulated (a) NH ice volume, (b) global annual mean surface air temperature and (c) globally averaged deep ocean (4 km) temperature.The grey lines in all figures show the Baseline Experiment.The black, blue and red lines indicate experiments with an acceleration by a factor of 5, 10 and 20, respectively.
Flow diagram of the model version used in this study.Thick arrows represent new flows of information compared to the "standard" version of CLIMBER-2 described in

Table 1 .
Abbreviations of information fluxes shown in Fig. 1.