Interdependence of the growth of the Northern Hemisphere ice sheets during the last glaciation: the role of atmospheric circulation

Abstract. The development of large continental-scale ice sheets over Canada and northern Europe during the last glacial cycle likely modified the track of stationary waves and influenced the location of growing ice sheets through changes in accumulation and temperature patterns. Although they are often mentioned in the literature, these feedback mechanisms are poorly constrained and have never been studied throughout an entire glacial–interglacial cycle. Using the climate model of intermediate complexity CLIMBER-2 coupled with the 3-D ice-sheet model GRISLI (GRenoble Ice Shelf and Land Ice model), we investigate the impact of stationary waves on the construction of past Northern Hemisphere ice sheets during the past glaciation. The stationary waves are not explicitly computed in the model but their effect on sea-level pressure is parameterized. We tested different parameterizations to study separately the effect of surface temperature (thermal forcing) and topography (orographic forcing) on sea-level pressure, and therefore on atmospheric circulation and ice-sheet surface mass balance. Our model results suggest that the response of ice sheets to thermal and/or orographic forcings is rather different. At the beginning of the glaciation, the orographic effect favors the growth of the Laurentide ice sheet, whereas Fennoscandia appears rather sensitive to the thermal effect. Using the ablation parameterization as a trigger to artificially modify the size of one ice sheet, the remote influence of one ice sheet on the other is also studied as a function of the stationary wave parameterizations. The sensitivity of remote ice sheets is shown to be highly sensitive to the choice of these parameterizations with a larger response when orographic effect is accounted for. Results presented in this study suggest that the various spatial distributions of ice sheets could be partly explained by the feedback mechanisms occurring between ice sheets and atmospheric circulation.


Introduction
The Quaternary era is characterized by a succession of glacial and interglacial phases.During the last ice age, large land-ice masses covered Canada and northwestern Eurasia (Peltier, 2004;Lambeck et al., 2006;Tarasov, 2012;Clark et al., 1993;Dyke and Prest, 1987;Svendsen, 2004).These large ice sheets represent a crucial element of the climate system (Clark, 1999).Because of their highly reflective surface and their high altitude, they induce zonal anomalies in surface temperature and topography.These anomalies modify large-scale atmospheric circulation by generating zonal asymmetries often referred to as stationary waves (Cook and Held, 1988).
The relations between atmosphere and ice sheets have been previously investigated with several modeling studies.Based on the analysis of simulations carried out with general circulation models (GCM) under Last Glacial Maximum (LGM) conditions, Broccoli and Manabe (1987) found that ice sheets are the main cause of change of stationary waves during ice age climate.In line with these previous findings, Pausata et al. (2011) have shown that the topography of ice sheets is the dominant factor altering the northern extratropical atmospheric large-scale circulation.Moreover, several Published by Copernicus Publications on behalf of the European Geosciences Union.
authors put forward that one of the main effects of an ice sheet on the atmospheric circulation is to change the strength and the position of the subtropical jet, which has also an influence on storm tracks (Kageyama and Valdes, 2000;Hall et al., 1996;Riviẽre et al., 2010;Laîné et al., 2008).This leads to changes in the pattern of precipitation and consequently to changes in the accumulation over ice sheets.
Changes in stationary waves also induce changes in surface temperature.Using a simple ice-sheet model based on an idealized geometry coupled with a stationary-wave model, Roe and Lindzen (2001a, b) highlighted the importance of accounting for the feedbacks between ice sheets and the temperatures induced by stationary waves to properly simulate the evolution of an ice sheet.In the same way, with a threedimensional stationary wave model, Liakka et al. (2011) showed that the southern margin of ice sheets strongly depends on the temperature anomalies due to stationary waves.All these studies illustrate the existence of feedbacks between ice sheets and atmospheric circulation.This suggests that the construction of a given ice sheet (e.g. the North American ice sheet) may influence the growth or the decay of the other one (e.g. the Eurasian ice sheet) through changes in atmospheric circulation that induce modifications in both temperature and precipitation patterns.In turn, these modifications directly influence the surface mass balance of the ice sheets.
However, up to now, no study has been undertaken to understand the role atmospheric circulation may play on icesheet evolution throughout a glacial-interglacial cycle.The aim of this study is twofold.First, we examine how the planetary waves may influence the evolution of ice sheets over the last glacial period.Secondly, we investigate how past Northern Hemisphere ice sheets (i.e.Laurentide ice sheet (LIS) and Fennoscandia ice sheet (FIS)) interact together through induced changes in planetary waves.To achieve this goal, we use the climate model of intermediate complexity CLIMBER 2.4 (Petoukhov et al., 2000) fully coupled with the 3-D thermo-mechanical ice-sheet model GRISLI (GRenoble Ice Shelf and Land Ice model) (Peyaud et al., 2007) to carry out a series of numerical experiments covering the last glacial-interglacial cycle.

The climate model CLIMBER 2.4
The CLIMBER 2.4 model is a revised version of the CLIMBER 2.3 model extensively described in Petoukhov et al. (2000).It is based on simplified representations of the ocean with three zonally averaged ocean basins for the Atlantic, Indian and Pacific oceans (2.5 • × 20 uneven layers), of the atmosphere (with resolution 10 • in latitude, 51 • in longitude), of the vegetation and of the mutual interactions between these three components.The main advantage of this model is its fast computational time allowing long-term simulations at glacial-interglacial timescale.
The atmospheric module is designed to only resolve largescale processes, the effect of synoptic weather systems on heat and moisture transports being parameterized.The atmospheric variables such as humidity and temperature are calculated in 2-D, and their 3-D profiles are then computed from hypotheses on the vertical structure of the atmosphere.Stationary waves are not explicitly resolved but their effect on sea-level pressure (hereafter SLP) is parameterized.In the standard version of the model, it is described by relationship between long-term large-scale azonal temperature and pressure fields in quasi-stationary planetary-scale waves adapted for the CLIMBER model: the longitudinal shift between the pressure and the temperature is neglected because of the coarse spatial resolution of CLIMBER.The azonal sea-level pressure (p 0 ) is expressed as a function of the azonal component of sea-level temperature T 0 (Petoukhov et al., 2000): with with where p 0 is the sea-level pressure, P 0 = 1012 hPa the mean of sea-level pressure.R = 287.058J kg −1 K −1 is the specific gas constant for dry air, g = 9.81 m s −1 is the gravity constant and H T is the computed tropopause height.Equation (2) shows that the influence of topography on the azonal SLP component is not explicitly accounted for.Here, the impact of orographic changes on stationary waves is studied through a new parameterization (described in Sect.3).The use of parameterizations to account for the effect of stationary waves on sea-level pressure offers the opportunity to test the influence of thermal and orographic effects separately.

The ice-sheet model GRISLI
GRISLI is a three-dimensional thermo-mechanical model which simulates the evolution of ice-sheet geometry (extension and thickness) and the coupled temperature-velocity fields in response to climate forcing.A comprehensive description of the model can be found in Ritz et al. (2001) and Peyaud et al. (2007).Here, we only summarize the main characteristics of this model.The equations are solved on a cartesian grid (40 km × 40 km).Over the grounded part of the ice sheet, the ice flow resulting from internal deformation is governed by the shallow-ice approximation (Morland, 1984;Hutter, 1983).The model also deals with ice flow through ice shelves using the shallow-shelf approximation  (MacAyeal, 1989) and predict the large-scale characteristics of ice streams using criteria based on the effective pressure and hydraulic load.At each time step, the velocity and vertical profiles of temperature in the ice are computed as well as the new geometry of the ice sheet.The temperature field is computed both in the ice and in the bedrock by solving a time-dependent heat equation.The surface mass balance is defined as the sum between accumulation and ablation computed by the positive degree day (PDD) method (e.g.Reeh, 1991;Fausto et al., 2009).This method assumes that melt rates of snow and ice are linearly related to the number of PDD through constant degree-day factors (Braithwaite, 1984;Braithwaite, 1995).The PDD method relies on the assumption that the annual surface melting is proportional to the sum over one year of the excess of temperature above the melting point.The number of positive degree days is calculated from the normal probability distribution around the mean daily temperatures and the variations around the daily mean (Reeh, 1991;Braithwaite, 1984).This approach has been widely used in ice-sheet models and is expressed as

Clim
where T is the near-surface temperature, T a (t) is the mean daily near-surface temperature, and σ the standard deviation around the daily mean.The mean daily temperature T a (t) follows a sine cycle during the year.This formulation allows positive temperatures in a given day even if the mean daily temperature is negative: the higher the σ value, the larger the probability of having positive temperatures is.Therefore, high σ values favor ablation.In the standard PDD formulation (Reeh, 1991), σ is fixed to 5 • C.However, the daily temperature variability is strongly dependent on the altitude.Based on measurements from automatic weather stations in Greenland, Fausto et al. (2009) derived a parameterization of σ expressed as where z s is the altitude of the ice-sheet surface and σ 0 the sea-level value; in this parameterization, σ 0 = 1.574 In the present study we used the same type of relationship between σ and the altitude but with different numerical values.We used different couples of σ 0 and α values (see Sect. 3) in order to modulate the amount of ablation and, thus, to simulate ice sheets of different sizes.

Coupling procedure between CLIMBER-2.4 and GRISLI
The mean annual and summer surface temperatures computed by CLIMBER, as well as snowfall, are used as inputs to GRISLI to compute surface mass balance.To account for the resolution difference between both models, we apply a specific downscaling procedure: for each CLIMBER grid box and each surface type, the temperature is computed on five vertical levels using the free atmospheric lapse rate to account for the dependency with the altitude.This vertical temperature profile is used to compute the vertical humidity and the resulting vertical precipitation profile.For temperature fields, these calculations are performed for each CLIMBER surface type and then averaged over the GRISLI surface types (land ice, ice-free land and ocean).The three climatic fields used as forcing are then tri-linearly interpolated on the ice-sheet model grid.GRISLI returns to CLIMBER new boundary conditions for each CLIMBER grid box: new surface type, altitude of each grid point and fresh water fluxes coming from the potential melting of ice sheets.A more detailed explanation of the coupling procedure is described in Charbit et al. (2013).
3 Experimental set-up

Parameterization of stationary waves in CLIMBER
As mentioned in Sect.2.1, in the standard version of CLIMBER, the azonal component of sea-level pressure only depends on sea-level temperature, so the topography does not exert any direct influence on sea-level pressure.The orographic impact on stationary waves depends on the topography and on the strength of the jet zonal wind (Held et al., 2002;Vallis, 2006;Holton, 1979).We parameterize the impact of the orography on stationary waves following the same basic principle as the one described by Eq. (2).Since the thermal contribution is expressed as a linear relation between the azonal component of sea-level pressure (p th ) and the azonal component of temperature (T 0 ), we expressed the orographic contribution (p oro ) as a function of the deviation of orography from its zonal mean (h 0 ).This accounts for the deviation of the zonal wind by an anomaly of topography.This new parameterization is expressed as follows: where P 0 = 1012 hPa is a constant, h topo is the altitude of the surface and h topo the mean zonal altitude.Since the zonal wind is driven by the equator-to-pole temperature gradient, we added a dependency to this term ( T E/P ).Adding this dependency allows to account for seasonal changes of the influence of the orographic effect on the sea-level pressure.
The larger this gradient is, the stronger the zonal wind and the influence of topography are.On the other hand, if T E/P is too small (i.e.below the T limit ), the zonal wind and its deviation due to orography will be negligible.Therefore, the T limit represents the limit below which the orographic contribution is negligible.
The numerical value of this parameter has been determined by comparing the present-day sea-level pressure obtained with the new parameterization (p 0 = p th + p oro ) with the NCEP reanalysis (Kalnay et al., 1996).The best correlation has been obtained for T limit = 25 • C.
Using this approach, we can study separately orographic and thermal effects, as well as a combination of both.We tested four SLP parameterizations: a. p 0 = p 0 (without any waves), b. p 0 = p 0 + p th (thermal forcing only), c. p 0 = p 0 + p oro (orographic forcing only), d. p 0 = p 0 + p oro + p th (thermal and orographic effect).
Figure 1 shows the azonal sea-level pressure for the northern winter for NCEP reanalysis (Fig. 1a) and for the three parameterizations B, C and D (Fig. 1b-d).The most striking feature is that the patterns of azonal sea-level pressure simulated by CLIMBER are smoother than the NCEP ones due to the spatial resolution.The thermal effect (Fig. 1b) allows to account for the main structures linked to the landsea temperature difference.Low pressure over North Pacific and North Atlantic and high pressure over the continental regions are clearly represented.However, their meridional extent is too large.This leads to a negative SLP anomaly over Greenland and a positive one over the Scandinavian and the Barents-Kara regions, in contradiction with NCEP reanalysis.Moreover, the amplitudes of the anomalies over the Northern Hemisphere are weaker than those of the NCEP database.This implies that the anti-cyclonic structure over the North American continent is almost absent.The amplitudes of sea-level pressure anomalies are smaller in the Southern Hemisphere (Fig. 1a) than in the Northern Hemisphere and occur over a smaller spatial scale.Due to the coarse horizontal resolution of CLIMBER, these structures are poorly resolved whatever the azonal sea-level parameterization is.With the thermal effect (Fig. 1b) this translates into a large anticyclone centered over South Atlantic that expands over the Indian Ocean, western Pacific and the Antarctic ice sheet.As a result, the bipolar structure over Antarctica observed in NCEP is not represented in the simulations carried out with the thermal forcing.Nevertheless, the Southern Hemisphere regions are beyond the regions of interest for the purpose of the present study.
With the orographic parameterization (Fig. 1c), a high pressure appears over Greenland, in agreement with Pausata et al. (2011).Subsequently this leads to a negative anomaly over the other regions of the same latitudinal band, especially over the Barents-Kara seas and over Scandinavia.The spatial structures are in a better agreement with NCEP but the amplitude of the negative anomaly over the Barents-Kara sector is still too weak.A high pressure is also simulated centered over the Caspian sea and over Central and North America.The high pressure over North America is weaker in our simulations than in NCEP.This is due to the Rocky Mountains poorly resolved in CLIMBER because of the zonal structure of the model.Finally, the bipolar structure over Antarctica is represented, although the amplitudes of low and high pressures over western and eastern parts respectively are weaker than the NCEP ones.This is also likely due to the spatial resolution.
Combining both effects (thermal and orographic, Fig. 1d) leads to amplitudes of SLP anomalies in better agreement with NCEP, especially over North Atlantic, North Pacific and Eurasia.Over the Barents-Kara seas, the negative anomaly due to the orographic effect clearly appears, but its amplitude is too small and overtaken by the influence of the thermal effect.As expected the high pressure over North America has a larger extent than that simulated under the thermal forcing alone, but its amplitude remains too weak with respect to the NCEP data set.

Modification of the ablation in GRISLI
Since the early PDD formulation (Reeh, 1991), a number of experimental campaigns over Greenland (Oerlemans and Vugts, 1993;Ambach, 1988) and other glacier locations have revealed strong spatial dependency of degree day factors (see Hock, 2003 andBraithwaite andZhang, 2000 for a compilation of degree day values).Moreover, owing to the fact that the daily temperature variability is larger in continental climate regions (e.g.Siberia) than in regions with oceanic climate (e.g.northern Europe), σ is unlikely to be spatially constant.In this study, σ is used as a tuning parameter which modulates the amount of ablation and allows to simulate more or less large ice sheets.To achieve that goal, we defined four different areas characterized by different σ values: North America (Alaska excluded) (σ LIS 0 ), Greenland (σ GIS 0 ), Fennoscandia (including British Islands) (σ FIS 0 ) and the rest of the grid (σ R 0 ).Each region is characterized by a specific σ 0 value, but the slope α is constant over the entire model grid.The amplitude of the daily temperature variability directly affects the number of positive degree days and, thus, the amount of ablation.Various values of σ 0 and α lead therefore to different shapes and sizes of simulated ice sheets.In this study, we choose σ 0 and α values in order to obtain a significant ice volume over North American and Eurasian ice sheets.In doing so, this method leads to differences in the size of both ice sheets which are large enough to investigate their mutual influence.

Description of the experiments
To explore the relationship between the Laurentide and the Fennoscandian ice sheets and to investigate how it is modulated by planetary waves, we used different values of σ FIS 0 and σ LIS 0 to obtain larger or smaller Fennoscandian and Laurentide ice sheets, respectively.This allows us to study the impact of the FIS geometry on the LIS and vice versa.We carried out three different experiments for the four parameterizations of sea-level pressure (see Sect. 3.1).These twelve experiments are summarized in Table 1.The initial climatic state is given by a time-slice CLIMBER experiment carried out for 126 ka conditions and the initial ice-sheet topography is set to that of the present-day Greenland.Transient simulations are forced by variations of insolation (Berger, 1978), atmospheric CO 2 concentration (Petit et al., 1999) and sea level (Waelbroeck et al., 2002).
This experimental set-up allows to study the evolution of the coupled Laurentide-Fennoscandia system when the growth of one or the other ice sheet is favored.This evolution is studied in response to each stationary wave parameterization described above.Ice volume changes are analyzed in terms of changes in temperature and accumulation fields.Over the last glacial period, changes in Greenland ice-sheet geometry are not large enough to induce changes in stationary waves, especially at our model resolution.Therefore, the impact of Greenland on the Laurentide and Fennoscandian is likely to be negligible and is not discussed in this paper.Moreover, with the model version used in the present study, we fail simulating the complete deglaciation from the Last Glacial Maximum (LGM) until 0 ka.Actually, previous works suggest that an increase of climatic variability  (Bonelli et al., 2009) or the impact of dust deposition on snow albedo (Ganopolski et al., 2010;Bonelli et al., 2009) may have played a critical role in the deglaciation process.Here, we did not account for the influence of these factors to better constrain the mutual influence of past Northern Hemisphere ice sheets driven only by changes in atmospheric circulation.We therefore kept σ FIS the impact of water-saturated sediments on basal sliding that favor the ice retreat has not been taken into account.These missing mechanisms may explain why the last glacial termination is not simulated satisfactorily in our experiments.Our following analysis therefore focuses on the 126-15 ka period.
Our aim here is not to simulate a perfect picture of the last glacial-interglacial cycle in full agreement with all available paleo-data, but rather to investigate the impact of one or the other stationary wave parameterization on the evolution of the ice sheets and their interdependence.

Impact of sea-level pressure on temperature and accumulation patterns
The calculation of sea-level pressure leads to modifications in winds which in turn influence humidity and heat transports, altering cloud formation, precipitation and air temperature.The relation between precipitation, temperature and sea-level pressure is complex due to numerous feedbacks involved in the climate system.Therefore, it is difficult to predict from SLP patterns what the temperature and accumulation patterns will be.However, we compared the CLIMBER response of the sea-level pressure to a change in topography to the response of the IPSL-CM5 model.We used two experiments carried out with both models (with the OTH SLP parameterization in CLIMBER) with exactly the same experimental set-up.The first one was run under preindustrial conditions (PI) and the second one under LGM ones (except that the only ice sheet considered as boundary condition was the Laurentide ice sheet).Both model responses exhibit similar SLP patterns (not shown), despite slight disagreements in some places as for the amplitude of high and low pressures.In the present section, we examine the impact of different SLP parameterizations on snow accumulation and surface air temperature and thus on the icesheet surface mass balance.This analysis is made at 125 ka because at this time period there is no large ice sheet in the Northern Hemisphere except Greenland.Moreover, the 125 ka climate state has likely influenced the early phase of glacial inception and thus the construction and the further development of both Laurentide and Fennoscandia.
Figure 2 displays the simulated summer temperature (Fig. 2a) and accumulation (Fig. 2e) when the parameterization of sea-level pressure is removed (i.e.NONE experiment) as well as the difference of summer temperature and annual accumulation between each SLP parameterization (i.e.TH, ORO and OTH experiments, Figs.2b-d and 2f-h).The azonal component in the SLP field produces a warming over a great part of the GRISLI domain (TH and OTH experiments) and over eastern longitudes and North Atlantic Ocean (ORO experiment), especially over Eurasia.Since the effect of the orography is weak in summer (see Sect. 3.1), its impact on summer temperature is less pronounced compared to the impact of the thermal effect (Fig. 2b and d).Over North America, where the orographic effects are accounted for, temperature patterns are dominated by the thermal effect (Fig. 2d).An excess of accumulation with respect to the NONE experiment is produced over western Eurasia with the thermal forcing whereas accumulation is reduced over Canada, Beringia and eastern Eurasia (Fig. 2f).The orographic parameterization has almost the opposite effect (Fig. 2g).The response of accumulation in OTH experiment is a combination of the two previous responses in which the regions of large accumulation rates are northwestern Canada, Barents and Kara seas and southern Europe.
These accumulation and temperature patterns have a direct impact on the construction of ice sheets during glacial inception.The spatial distributions of the simulated ice sheets are displayed in Fig. 3 for 115 ka.At this period, the glacial inception has already started as suggested by the sea-level reconstructions (Waelbroeck et al., 2002).This figure shows that the ice sheet response strongly depends on the SLP parameterization.
Compared to the NONE experiment, results from the ORO simulation show that the orographic effect favors the growth of the Laurentide ice sheet (Fig. 3c).This is due to the large accumulation rate over Canada (Fig. 2g) combined with a cooling effect (Fig. 2c).The Fennoscandian ice sheet is rather sensitive to the thermal forcing (Fig. 3b) which produces a large accumulation rate over Scandinavia (Fig. 2f).However, the simulated Fennoscandian ice sheet is smaller than that simulated by the NONE experiment because the large accumulation rate is widely counterbalanced by the temperature warming effect (Fig. 2b).Combining both SLP perturbations (thermal and orographic) (Fig. 3d) leads to a much smaller ice volume over Scandinavia (with respect to TH experiment) due to less accumulation, and more ice over the Kara Sea (increased accumulation).Likewise, the Laurentide ice volume simulated on the OTH experiment is in between the TH and the ORO ones.
These first results highlight how the representation of stationary waves in CLIMBER (through SLP parameterizations) may influence accumulation and temperature patterns and therefore the construction of the ice sheets.They show that the simulated ice sheets are heavily different depending on the forcing effect that we take into account (i.e.thermal and/or orographic).In the following sections, we examine the mutual interactions of ice sheets under different SLP parameterizations.

Glacial inception: from 126 to 110 ka
Figure 4 displays the evolution of the simulated LIS and FIS ice volumes from 126 to 110 ka for the twelve experiments described in Sect.3.For a given parameterization, the comparison between solid and dotted lines (Fig. 4a-b) and between solid and dashed lines (Fig. 4c-d) illustrates the difference between the standard simulations (referred to as X-REF in the following) and the simulations where σ FIS 0 (Fig. 4a-b) or σ LIS 0 (Fig. 4c-d) is reduced.

Effect of smaller σ FIS 0
As expected, for all SLP parameterizations, decreasing σ FIS 0 favors the growth of the FIS (Fig. 4a).However, the amplitude of this growth is more or less pronounced depending on the parameterization used.When stationary waves are off (i.e.azonal SLP component is set to zero), the simulated NONE-FIS ice volume (i.e.σ FIS 0 reduced) is around 30 % larger at 110 ka than the ice volume simulated in the baseline experiment.This increase is even more pronounced in TH and OTH experiments (270 % of ice volume difference between .In the ORO experiment, the ice volume is increased by 100% compared to the reference case, but due to the small amount of ice, this increase is not really significant.Figure 4b illustrates how a change in the Fennoscandian ice-sheet geometry influences the construction of the Laurentide ice sheet.Compared to the Laurentide ice volume simulated under orographic forcing (15 × 10 15 m 3 ), the difference in Fennoscandian ice volumes between ORO-FIS and ORO-REF is not large enough (0.5 × 10 15 m 3 at 110 ka) to substantially modify the construction of the Laurentide ice sheet.By contrast, when stationary waves are off or when thermal effect is accounted for, the increase of the Fennoscandian ice volume leads to an increase of the volume of the Laurentide ice sheet.This effect may be directly related to the cold temperature anomaly occurring in response to the development of the Fennoscandian ice sheet.The comparison of summer temperature between NONE-FIS and NONE-REF (  progressively spreads all over the GRISLI domain, favoring thereby the growth of the Laurentide ice sheet.The spreading of the cold temperature anomaly can also be observed in OTH and TH experiments.This suggests that the larger the Fennoscandia ice sheet, the greater the effect on the Laurentide ice sheet.However, the differences of Fennoscandian ice volumes between TH-REF and TH-FIS on one hand and between OTH-REF and OTH-FIS on the other hand are similar. However, the effect on the Laurentide ice sheet is twice as large under thermal forcing alone than the effect produced with the OTH experiment.This means that, in addition to the cooling effect (Fig. 5), another mechanism comes into play.Actually, in OTH experiments, the growth of Laurentide ice sheet leads to a high pressure over Canada (Fig. 6ab).This causes a westward shift of the accumulation area (Fig. 6c-d Laurentide ice sheet.Note that this westward shift is also observed in the ORO experiments (not shown).Under thermal forcing, a slight extent of the high pressure area is also simulated with the growth of the Laurentide ice sheet.However, this extent remains insufficient to act on the displacement of the accumulation area.This makes the TH-FIS Laurentide ice sheet more sensitive to the cooling effect induced by a larger Fennoscandia than the ice sheet simulated in OTH-FIS.This also explains why the differences in simulated LIS ice volumes between TH-REF and TH-FIS are larger than those between OTH-REF and OTH-FIS.

Effect of a smaller σ LIS 0
The reduction of σ LIS 0 has a smaller effect on the LIS than the effect of a lowered σ FIS 0 value on the FIS (Fig. 4c).Although, σ LIS of σ LIS 0 is large enough to produce a larger LIS (Fig. 4c).The ice-sheet growth is smaller in ORO-LIS and OTH-LIS experiments due to the westward shift of the accumulation mentioned above.In all experiments, the differences in the ice-sheet size (with respect to the standard simulations) are too small to significantly influence the construction of FIS (Fig. 4d), at least during glacial inception.

Full glacial state: from 80 to 30 ka
The evolution of the simulated FIS and LIS ice volumes over the entire last glaciation are displayed in Fig. 7. Around 70 kyr BP, there is a minimum of insolation in northern latitudes combined to a decrease of the atmospheric CO 2 concentration which allows a rapid decrease of summer temperature that accelerates the glaciation process.

0
A large difference is observed in Fennoscandian ice volume (Fig. 7a) between ORO-REF (red solid line) and ORO-FIS (dotted-red line) from 75 ka to the end of the simulation.Actually, slight differences in simulated ice volumes between ORO-REF and ORO-FIS appear before 80 ka.They induce small differences in accumulation and temperature patterns as soon as 80 ka (Fig. 8).A small excess of accumulation over the Kara Sea is simulated in ORO-FIS (with respect to ORO-REF) as well as slightly colder temperatures over the Eurasian region.These tiny differences are sufficient to trigger a massive glaciation of the FIS in the ORO-FIS experiment, whereas in ORO-REF the ice volume remains at a low level (< 3.0 × 10 15 m 3 ) until 30 ka.These results show how the decrease of insolation and atmospheric CO 2 concentration may be amplified by a small change in the ice-sheet surface mass balance.The development of both LIS and FIS in the ORO-FIS experiment leads to a cooling all over the GRISLI domain.This shifts southward the snow-rain limit (with respect to ORO-REF), increasing the accumulation rate in the North Atlantic and eastern Canada.This effect combined with the previous cooling effect leads to an acceleration of the LIS growth (Fig. 7b, red solid line with respect to red dashed line).
A similar behavior is observed for the FIS ice volume when both forcings are accounted for (OTH experiment) with an acceleration of the FIS growth in OTH-FIS after 75 ka (Fig. 7a, green lines).The underlying mechanism is similar to the one explained in the ORO case.However, it is more efficient because the differences at 80 ka between the LIS ice volumes simulated in the OTH-REF and OTH-FIS experiments are larger than the differences between ORO-REF and ORO-FIS (Fig. 7b).Similarly to the ORO experiments, the acceleration of the FIS growth implies an acceleration of the (with respect to OTH-REF).
Under the thermal parameterization, there is no difference in the FIS ice volumes between TH-REF and TH-FIS throughout the period spanning from 80 to 30 ka (Fig. 7a, blue lines), whereas the LIS is larger (TH-FIS experiment, Fig. 7b).As explained in Sect.4.2.1, at the beginning of the glaciation, the decrease of σ FIS 0 leads to a larger FIS and in turn to a larger LIS.Actually, the LIS has remained larger in TH-FIS (with respect to TH-REF), even during periods where FIS was entirely melted (e.g. 100 ka).This means that the LIS ice volume differences between TH-FIS and TH-REF after 75 ka does not result from a direct effect of a change in FIS geometry, but rather comes from the fact that the LIS has remained glaciated throughout this period even around 80 ka when the FIS ice volume was very small but not zero.
Finally, we obtain a surprising result when the stationary waves are off: after 75 ka, the NONE-FIS experiment simulates a smaller FIS than the NONE-REF experiment (Fig. 7a, grey and black lines), despite a smaller σ FIS 0 .The comparison of accumulation and temperature patterns between both experiments (not shown) does not explain the decrease of the ice volume after 90 ka in the NONE-FIS experiment (Fig. 7a).Instead, a strong decrease of the simulated ice thickness in the southwestern part of the ice sheet (NONE-FIS) is associated with a decrease of ice flow velocities (not shown).This suggests that the mechanisms responsible for this behavior are closely linked to ice dynamical effects.However, an in-depth analysis of the ice dynamics is beyond the scope of this study.

Effect of a smaller σ LIS 0
At the beginning of the glaciation, the decrease of σ LIS 0 does not have a significant influence on the growth of the ice sheets, at least until 75 ka (Fig. 7c and d).After 75 ka, when the LIS is large (NONE and ORO experiments), the decrease of σ LIS 0 has only a poor impact on the ice volume.With the TH and OTH parameterizations, the geometry of simulated ice sheets in both reference runs and runs carried out with a smaller σ LIS 0 value is similar around 80 ka (TH and OTH experiments).However, after 75 ka, the differences of ice volumes between TH-REF and TH-LIS and between OTH-REF and OTH-LIS become significant, in contrast to results obtained in the ORO and NONE experiments.The impact of a larger LIS on the evolution of FIS is clearly visible in the OTH runs but remains very weak with the TH forcing.To understand the origin of this different behavior, we examined the differences between OTH-LIS and OTH-REF and between TH-LIS and TH-REF in terms of summer temperature and accumulation patterns.To be relevant, this comparison has to be made at a time period when the differences of ice volumes between reference runs and small σ LIS 0 runs (hereafter called ice) are roughly of similar magnitude in both OTH and TH experiments.Similar ice values are obtained at 75 (OTH) and 71 ka (TH).Figure 9 displays the differences in SLP, accumulation and summer temperature patterns at these periods.Since we take into account the topography effect on sea-level pressure, the response of a similar LIS ice on sea-level pressure is stronger in the OTH than in the TH experiment (Fig. 9a and b).Changes in accumulation are larger in the OTH case (Fig. 9e and f).A slight excess of accumulation is observed over the Kara Sea region in OTH-LIS experiment (with respect to OTH-REF), which seems to be sufficient to trigger the growth of FIS.Nevertheless, the impact on temperature in the TH experiment is more important than in the OTH experiment, due to a larger LIS ice (Fig. 9c and d).
These results clearly show that the sensitivity of the atmospheric circulation to the presence of an ice sheet becomes significant when the topography effect is accounted for.Moreover, a change in sea-level pressure due to the growth of an ice sheet seems to mainly alter the accumulation with only a little influence on summer temperature.

30 until 15 ka BP
After 30 ka, a new decrease in insolation concomitant with a decrease in atmospheric CO 2 concentration occurs.The response of the ice volume is more or less pronounced depending on the experiment (Fig. 7).As an example, a pronounced response of both ice sheets (Fig. 7a and b, blue dotted line) is simulated in TH-FIS, unlike NONE-REF, NONE-FIS and NONE-LIS.More generally, it should be noted that the response is absent or weak in simulations where both ice sheets have already reached a significant size at 35 ka (more than 18 × 10 15 m 3 for the FIS and more than 28 × 10 15 m 3 for the LIS).This can be explained by the dry cold air above continental-scale ice sheet.In contrast, when one ice sheet is large but the other is small (e.g.ORO-REF) at 35 ka, colder temperatures due to the decreasing of CO 2 atmospheric concentration and insolation favor the growth of the smaller ice sheet (e.g.FIS start to grow at 30 ka).In turn, this leads to even much colder temperatures and a further growth of the larger ice sheet too (LIS starting to grow at 35 ka).Therefore, the response of one ice sheet to an insolation/CO 2 decrease also depends on the size of both ice sheets.

Conclusions
In the present we investigated the atmospheric-based processes that relate the two main Northern Hemisphere ice sheets during the last glaciation.In particular, through appropriate parameterizations, we examined the effect of topography and surface temperature on the sea-level pressure, and we studied the impact of these effects in the relationship between ice sheets.The parameterization of ablation in the icesheet model is used as a trigger to change the size of one ice sheet and to investigate the mutual influence of the ice sheets.
Our aim was to investigate the effect of a change in one ice-sheet geometry on the other one, depending on the sealevel pressure parameterization.First, at the beginning of the glaciation, we showed that the growth of the Laurentide ice sheet is favored by the topographic effect on sea-level pressure, due to a large accumulation rate in this area, whereas there is more accumulation over Fennoscandia with the thermal parameterization.Secondly, we showed that a larger ice sheet (e.g.FIS) leads locally to a cooling that progressively spreads over the northern latitudes and promotes the growth of the other ice sheet (e.g.LIS).This former mechanism is true whatever the sea-level pressure parameterization is.The second process is the shift of snowfall pattern when an ice sheet grows, altering thereby the accumulation over the other ice sheet.This mechanism is dominant when the orographic effect is on, because it increases the sensitivity of sea-level pressure to the presence of an ice sheet.In that case a change in sea-level pressure mainly affects the accumulation patterns.Owing to the fact that stationary waves are not explicitly computed in CLIMBER and due to the low resolution of the model, the importance of various feedbacks discussed in the present study as well as the amplitude of remote effects of ice sheets may be under-or overestimated.However, this study highlights the key role of topography on sea-level pressure and accumulation, and thus on the evolution of the ice sheets themselves.This means that it is necessary to take into account the mutual influence of past Northern Hemisphere ice sheets to properly understand the mechanisms underlying their own evolution.In the same way, our results suggest that feedbacks between ice sheets and stationary waves could be of key importance to understand the various configurations of ice-sheet shapes during different ice ages of the Quaternary era (e.g.Svendsen, 2004).

Fig. 1 .
Fig. 1.Difference between sea-level pressure and its zonal mean in winter (DJF) for (a) NCEP reanalysis, (b) with only thermal parameterization (TH), (c) with only orographic parameterization (ORO) and (d) for the sum of the two parameterizations (OTH) (hPa).

Table 1 .
Experimental set-up.For each experiment, the 1st part of the name (NONE, TH, ORO, OTH) corresponds to the stationary wave parameterization while the second part corresponds to the ablation formulation: REF indicates the couple σ LIS 0 − σ FIS 0 for the baseline experiment, FIS correspond to a PDD parameterization favoring the growth of Fennoscandian ice sheet, and LIS to a PDD parameterization favoring the growth of Laurentide ice sheet.

Fig. 2 .
Fig. 2. Summer surface temperature at 125 ka in NONE-REF experiment (parameterization of SLP are removed) (a), difference of summer surface air temperature between TH-REF (thermal parameterization) and NONE-REF (b); ORO-REF (orographic parameterization) and NONE-REF (c) and OTH-REF (combination of both thermal and orographic parameterization) and NONE-REF (c) at 125 ka.Same for the snow accumulation (e-h).The green line represents the limit where the ice thickness difference exceeds −500 m.

Fig. 4 .
Fig. 4. The left part of the figure shows the evolution of Fennoscandian ice volume (a) for the reference experiments (solid lines, black when the waves are off, blue for the thermal parameterization, red for the orographic parameterization, and green when the two effects are on) and for the experiments with a smaller σ FIS 0 (plotted line, same color code).(b) shows the Laurentide ice volume for the same experiments.The right part of the figure shows the evolution of the Laurentide ice volume (c) for the reference experiments (solid lines) and for the ones with smaller σ LIS 0 (dotted lines), then the Fennoscandian ice volume (d) for the same experiments.For all panels, the evolution of ice volume is shown for the inception.

Fig. 5 .
Fig. 5. Summer surface temperature differences between NONE-FIS and NONE-REF at 125 ka (a) and 120 ka (b).The yellow line represents the limit where the ice thickness difference exceeds 500 m.

Fig. 6 .
Fig. 6.Azonal sea-level pressure (in hPa) at 120 and 115 ka for OTH (a-b) and TH experiments (e-f), respectively, and accumulation (in meter water equivalent) at 120 and 115 ka for OTH (c-d) and TH experiments (g-h).Grey lines indicate the limits of the ice sheets.

Fig. 7 .
Fig. 7.The left part of the figure shows the evolution of Fennoscandian ice volume (a) for the reference experiments (solid lines, black when the waves are off, blue for the thermal parameterization, red for the orographic parameterization, and green when the two effects are on) and for the experiments with a smaller σ FIS 0 (plotted line, same color code).(b) shows the Laurentide ice volume for the same experiments.The right part of the figure shows the evolution of the Laurentide ice volume (c) for the reference experiments (solid lines) and for the ones with smaller σ LIS 0 (dotted lines), then the Fennoscandian ice volume (d) for the same experiments.For all panels, the evolution of ice volume is shown for the all glaciation.
Fig. 8. (a) Accumulation difference (in meter water equivalent) between ORO-FIS and ORO-REF at 80 ka.The yellow line is the limit where the thickness difference exceeds 500 m, the green one where the difference is under −500 m.(b) Summer temperatures differences between ORO-FIS and ORO-REF at 80 ka (in • C).

Fig. 9 .
Fig. 9. Sea-level pressure differences (in hPa) between OTH-LIS and OTH-REF at 75 ka (a) and between TH-LIS and TH-REF at 71 ka (b).Summer temperature difference between OTH-LIS and OTH-REF at 75 ka (c); and TH-LIS and TH-REF at 71 ka (d); accumulation difference between OTH-LIS and OTH-REF at 75 ka (e); and TH-LIS and TH-REF at 71 ka (f).