INTRODUCTION

To describe the Earth climate system (ECS) in more detail, ice sheets (IS) should be included in numerical models [1–4]. Taking into account the inertial (and principally nonlinear) components may lead to both the formation of new feedback [5] and the development of multistability in the ECS [1]. The IS regime is a significant integral indicator for a critical level in the climate changes [6]. IS melting during the continued anthropogenic climate changes caused the ocean level to rise: 21 ± 2 mm in 1992‒2020; by 2150, their melting may increase this level additionally by 5 m [7]. Freshening of the World Ocean related to IS melting in the regions where deep-water convection is formed contributes to weakening of an oceanic conveyor with potentially catastrophic consequences for the Earth system [7, 9]. The regional consequences related to the IS formation and melting in the ECS may manifest themselves even tens of thousands of years later. In recent years, the formation of conical depressions (craters) on Yamal and adjacent regions was revealed. In [8], their formation is associated with the decomposition of methane hydrates at a shallow depth, which were formed under high pressure when these regions were covered by an ice sheet for tens of thousands of years.

The realistic reproduction of ISs depends very heavily on the uncertainty of paleodata for climate changes in the Pleistocene. In particular, according to [10], this uncertainty is significant, especially at a regional level compared to strong cooling at a global level during the Last Glacial Maximum (LGM) [10, 11].

The ice sheet models are compared to the direct data about the boundaries of distribution of the main ISs in the LGM and during other time intervals in the Pleistocene, as well as to independent data about the decrease in the ocean level. This comparison makes it possible to formulate the problem about the direct calculation of the influence exerted by the uncertainty in the paleoreconstructions on the IS dynamics followed by the estimation of the upper boundary of the indicated uncertainty.

The purpose of this work is to analyze the results of the ensemble calculations with the IS dynamics for the last glacial cycle and to estimate the influence exerted by uncertainty of paleoreconstructions on the ice sheet dynamics.

THE MODEL OF ICE SHEETS IN USE AND SIMULATIONS

We used the modification of the IceBern2D vertically isothermal ice sheet model described in [12]. The initial variant of the model is characterized by extremely strong growth of ice sheets near the mountain massifs, including at the boundary of the calculated domain. This is caused by the impossible transport of an ice mass across the boundary in the code of the model; therefore, even a small excess of ablation over melting can accumulate a significant ice mass per unit area for tens of thousands of years. In this context, the accumulation of the ice mass was not allowed in Himalaya and in the first two cells at the boundary of the calculated domain (in the direction perpendicular to this boundary). The more intense growth of the ice sheets can be associated with the thin ice approximation used in the model [13] (which is similar to the shallow water approximation in hydrodynamics); consequently, prohibiting accumulation of an ice mass in the indicated regions is an attempt to describe roughly the steep boundaries of the ice sheet and the input of precipitates for the formation of enclosed flow. The equations of the model are solved on the Arakawa C stereographic computational grid with a 40 km horizontal resolution (this resolution is at the upper boundary of the allowable value for the typical IS models [3]) and the time step of one year.

The model was used for simulations of IS evolution for the last 128 kyr. In the simulations, the initial sea level is +7 m relative to the modern level and ISs are not present in the Northern Hemisphere, which corresponds to the state of the Earth system 128 ka (the so-termed the Eemian Interglacial). Similarly to [12], the  climate for this time interval was reconstructed by  using the data of the δ18O content in the ice cores   at the EPICA Dome C Antarctic Station (https://www.ncei.noaa.gov/pub/data/paleo/icecore/ antarctica/epica_domec/edc3deuttemp2007.txt):

$${{Y}_{{{\text{ctrl}}}}}(t)\, = \,{{Y}_{0}}\, + {{\delta }^{{18}}}{\text{O}}(t)({{Y}_{{{\text{LGM}}}}}\,-{{Y}_{0}}){\text{/}}({{\delta }^{{18}}}{{{\text{O}}}_{{{\text{LGM}}}}}\,-{{\delta }^{{18}}}{{{\text{O}}}_{0}}).$$
(1)

Here, Y is the temperature or precipitates, subscript “0” corresponds to the modern regime, subscript “LGM” corresponds to the LGM regime 21 ka. The modern regime was characterized by the long-term mean values based on the ERA5 reanalysis data (https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-pressure-levels-preliminary-back-extension?tab=form) for 1950‒1978. The LGM regime was characterized by the long-term mean values based on the results of the relevant equilibrium calculation with the IPSL-CM4-V1-MR climate model, which was performed under PMIP3 International Project (Paleoclimate Models Intercomparison Project, phase 3) (https://pmip3.lsce.ipsl.fr/). This simulation with the IS model is called the control (CS) and is marked by subscript “ctrl.” We note that, in the general case, the suggested approach may lead to deviations of initial conditions for the model integration from reconstructions for the Mikulino Interglacial.

In addition to control simulation with the model, we performed ensemble simulations with perturbation of temperature and precipitation in accordance with

$$Y{\kern 1pt} '\left( t \right) = {{Y}_{{{\text{ctrl}}}}}\left( t \right) + A\sin [t{\text{/}}(2\pi T) + {{\Phi }_{0}}],$$
(2)

with the period T in the range from 120 to 2400 yr and amplitude A, which varies from 1 to 20% of the maximum difference in the temperature (precipitation) between the modern period and the period of the Last Glacial Maximum. According to (2), the maximum temperature perturbation is about 2°С, which agrees with the estimates of uncertainty of paleoreconstructions for the Last Glacial Maximum [10]. Along with the period and amplitude of perturbations in (2), the initial phase Φ0 varied in the range from 0 to 2π as well.

RESULTS

Reconstruction of the last glacial cycle. Based on the results of our simulations for the modern period in the presence of only the Greenland IS in the Northern Hemisphere, the largest thickness Н of the ice sheet in terms of space is close to 3 km. This agrees with the actual data (up to 3.4 km). In this case, the propagation area of the ice sheets calculated as the total area of the regions with H > 0 equaled 1.85 million km2, which is close to the observed area (1.83 million km2).

Figure 1 presents the changes in the global ocean level in the control simulation with the IS model compared to the paleoreconstruction data [14] for the last glacial cycle. According to Fig. 1, the most noticeable differences in the model simulations and the results of the paleoreconstructions manifest themselves at the beginning of the period under analysis, i.e., the first third of the last glacial cycle, while for the rest of the last glacial cycle, the conformity is much better. This is naturally associated with the influence of initial conditions that do not conform sufficiently well to the conditions of the Eemian Interglacial. A similar effect is recorded in particular for thermophysical processes in the bottom sediments of the Arctic shelf, another inertial component of ECS [15].

Fig. 1.
figure 1

Changes in the global ocean level in the control simulation (black curve) with the model of the ice sheets compared to the paleoreconstruction data [14] (red curve).

Full size image
Fig. 2.
figure 2

Thickness of the ice sheet in the control simulation (a) for the modern conditions and (b) for the conditions of the Last Glacial Maximum 20 ka as well as (c) in the GLAC-1 paleoreconstruction.

Full size image
Fig. 3.
figure 3

Thickness of the ice sheets in the Northern Hemisphere averaged for the time interval of 58‒51 ka (a) in the control simulation and in the simulations at (b) А = 0.2, Т = 2400 yr, Φ0 = π, and (c) А = 0.2, Т = 2400 yr, Φ0 = 0.

Full size image

During the Last Glacial Maximum (~20 ka), according to the model simulations, the World Ocean level decreased by 113 m compared to the modern regime (Fig. 1). This value falls into the range of uncertainty for the paleoestimates of the total contribution of the Laurentide, Scandinavian, and Greenland Ice Sheets, as well as smaller-sized glaciers in the LGM (from ‒95 to ‒113 m) [16]. The contribution of the Antarctic IS and the change in the ocean level due to the increasing density of the sea water is not considered in this work. The calculated values of the heights of ice sheet domes during the Last Glacial Maximum (to 5.3 km for the Laurentide IS and 3.5 km for the Scandinavian IS) agree in general with the PaleoMIST 1.0 data (4.2 and 3 km, respectively; see https://doi.pangaea.de/10.1594/PANGAEA.905800). The calculated height of the dome for Eurasia is also consistent in general with the GLAC-1 model reconstruction (3.1 km; one of the variants of the ICE-6G reconstruction; in detail https://pmip3.lsce.ipsl.fr/), but the values for North America are slightly smaller according to the data of the same reconstruction (3.5 km).

Compared to the modern regime, the model extension in the propagation area of the ice sheets in the LGM (23‒18 ka) is 17.5 million km2, including 13.9 million km2 in the Western Hemisphere (the Laurentide IS and the Greenland IS) and 3.6 million km2 in the Eastern one (the Scandinavian IS). This area is smaller than the area derived according to the area obtained based on the GLAC-1 data (a total of 21.8 million km2, including 16.5 million km2 in the Western Hemisphere and 5.3 million km2 in the Eastern Hemisphere, respectively).

For the interval of 51‒58 ka (the closest to the first third of MIS3 (marine isotope stage 3)), the model heights of the domes on the Laurentide IS, the Greenland IS, and the Scandinavian IS were equal to about 4.5 km, 2.5 km, and 1.5 km, respectively. For the first Laurentide and Greenland ice sheets, this is consistent with the results in [17], while the value obtained in this work for the height of the dome on the Scandinavian ice sheet is much lower than the corresponding value in [17]. We note that, according to the ICESHEET 2.0 data, in the Western Hemisphere, the maximum height of the ice sheet for the time interval of 51‒58 ka is lower than that obtained by the model simulation (3.6 km), while in the Eastern Hemisphere, it is greater (2.6 km). The propagation area of these ice sheets in the Northern Hemisphere was estimated as 16.6 million km2 in the model simulations, including 13.7 million km2 in the Western Hemisphere and 2.9 million km2 in the Eastern Hemisphere.

The influence of climate variations on the ice sheet dynamics in the Northern Hemisphere. The climate variations with amplitude А = 0.1 (which corresponds to uncertainty of the near-surface global temperature of ~1°С) do not significantly change the dynamics of the ice sheets of the Northern Hemisphere in the last glacial cycle. At a larger amplitude А = 0.2 (which corresponds to uncertainty of the near-surface global temperature of ~2°С), the corresponding climate variations lead to quality changes in the ice sheet dynamics and the ocean level in the past 120 kyr (Fig. 4). In particular, at the initial phase Φ0 = 0, during the Last Glacial Maximum, the depth of the ocean level fall changes by 30 m with respect to the period Т, which is approximately by a fourth of the minimum sea level depth in the control simulation. In this respect, the IS propagation area in the Northern Hemisphere during the LGM is not changed significantly compared to the results of the control simulation.

Fig. 4.
figure 4

Change in the ocean level in the numerical calculations with the model of the ice sheets (a) versus the period T (at the initial phase Φ0 = 0) and versus the initial phase Φ0 (at (b) T = 1200 yr and (c) T = 2400 yr) of temperature perturbations with amplitude А = 0.2.

Full size image

The most noticeable changes were recorded for the interval of 51‒58 ka, when the orbit configuration of the Earth was close to the threshold to leave the glaciation regime [2]. For this time interval, the positive temperature variations may lead almost to total melting of the ice sheets in the Northern Hemisphere, except for the Greenland IS (Figs. 3, 4). In the numerical simulation with А = 0.2, Т = 2400 yr, and Φ0 = π with a minimum spread of the ice sheets in the indicated time interval, their area is 14 million km2, including 12 million km2 in the Western Hemisphere and 2 million km2 in the Eastern Hemisphere with a maximum height of the domes equal to 3.2 and 1 km, respectively. In this case, the average thicknesses of the Laurentide and Scandinavian ISs are 0.7 and 0.3 km, respectively.

The changes in the initial phase Φ0 of climate variations may lead to the development of the cold anomaly of the temperature at this time interval followed by the intensification of glaciation in the Northern Hemisphere and the ocean level depth compared to the level reached in the Last Glacial Maximum. For example, in the numerical simulation at А = 0.2, Т = 2400 yr, Φ0 = 0 with the maximum spread of the ice sheets on the time interval indicated, their area is 18.6 million km2, including 15.2 million km2 in the Western Hemisphere and 3.6 million km2 in the Eastern Hemisphere with the maximum height of the domes equal to 5.2 and 3.1 km, respectively. In this case, the average thicknesses of the Laurentide and Scandinavian ice sheets are close to 2 and 1 km, respectively. The thickness and spatial structure of the Greenland IS do not change significantly in these simulations compared to the control simulation.

When the period of climate variations Т = 1200 yr, according to the model results with respect to the initial phase Ф0, either the interglacial is formed or the glaciation intensifies, up to the level comparable to that reached in the LGM. Thus, the model results obtained have low sensitivity to the period of climate variations, while their amplitude and the initial phase play a much bigger role.

The formation of the model interglacial at MIS3 and excessive glaciation on the same time interval are examples of potential anomalous dynamics. This does not manifest itself in the model simulations at the global amplitude of temperature variations (perturbations) ≤1°С but is likely to occur at the relevant amplitude ≥2°С. The latter estimate may serve as the upper estimate for the conditions of manifestation of the significant uncertainty in the regimes of ice sheets with respect to climate variations. In general, this estimate is in agreement with the corresponding estimate [10].

CONCLUSIONS

The ensemble simulations were carried out with the model of thte ice sheet dynamics for the last glacial cycle (128 kyr) in the Northern Hemisphere. According to the results of the numerical calculations, in general, the model realistically reproduces the spatial structure of the ice sheets in the Northern Hemisphere and the heights of their domes, as well as the associated changes in the ocean level. The most significant differences in the model simulations and the results of the paleoreconstructions were recorded at the beginning of the study period, which is naturally related to the influence of the initial conditions. The MIS3 and LGM regimes were simulated realistically.

When possible paleoclimate variations (perturbations) with quite large amplitude are taken into account, the results of simulating the ice sheets in the Northern Hemisphere show noticeable sensitivity. In particular, the depth of the decrease in the ocean level during the LGM changes by approximately a fourth of its value in the control simulation; the Laurentide IS and Scandinavian IS either melted almost completely at the MIS3 or their areas and volumes increased to the values comparable to those reached during the LGM. The analysis of the conditions for implementing such regimes of ice sheets dynamics makes it possible to estimate the possible range of uncertainty of the paleoclimate reconstructions. According to our results, the significant uncertainty about the ice sheet regimes for the last glacial cycle is observed at temperature variations of about 2°С, which is consistent with the available estimates. Except for the cases indicated in this work, the data uncertainty does not lead to qualitative change in the system.

To refine the estimates obtained, the model of ice sheets should further include a module of thermophysical processes within the ice with respect to, among others, the vertical temperature profile inside the ice sheet and its interactive calculation.

A significant limitation of this work is the use of the vertically isothermal ice sheet model. In particular, the consideration of the vertical temperature profile inside the ice sheet (and its interactive calculation in the model) can affect quantification in this work, including for the upper limit of uncertainty of temperature reconstructions in the Pleistocene. Further, it is planned to expand the model by using the block of thermophysical processes within the ice with an option to refine the estimates obtained.