Transient simulation of the last glacial inception. Part I: glacial inception as a bifurcation in the climate system

Abstract

We study the mechanisms of glacial inception by using the Earth system model of intermediate complexity, CLIMBER-2, which encompasses dynamic modules of the atmosphere, ocean, biosphere and ice sheets. Ice-sheet dynamics are described by the three-dimensional polythermal ice-sheet model SICOPOLIS. We have performed transient experiments starting at the Eemiam interglacial, at 126 ky BP (126,000 years before present). The model runs for 26 kyr with time-dependent orbital and CO2 forcings. The model simulates a rapid expansion of the area covered by inland ice in the Northern Hemisphere, predominantly over Northern America, starting at about 117 kyr BP. During the next 7 kyr, the ice volume grows gradually in the model at a rate which corresponds to a change in sea level of 10 m per millennium. We have shown that the simulated glacial inception represents a bifurcation transition in the climate system from an interglacial to a glacial state caused by the strong snow-albedo feedback. This transition occurs when summer insolation at high latitudes of the Northern Hemisphere drops below a threshold value, which is only slightly lower than modern summer insolation. By performing long-term equilibrium runs, we find that for the present-day orbital parameters at least two different equilibrium states of the climate system exist—the glacial and the interglacial; however, for the low summer insolation corresponding to 115 kyr BP, we find only one, glacial, equilibrium state, while for the high summer insolation corresponding to 126 kyr BP only an interglacial state exists in the model.

where ci is the effective heat capacity of snow or ice surface, ρ w is the water density. Sabs, R↓, R↑, H and LE are the absorbed short-wave solar radiation, downward long-wave, upward long-wave radiation, sensible and latent heat fluxes, respectively. Equation 3 includes snowfall Ps, ablation A and sublimation E. Ablation is computed from the energy balance Eq. 2, if the surface temperature Ts exceeds the melting point. Unlike Wang and Mysak (2002), our model does not include parameterisations for freezing of rain or refreezing of melt water. The amount of absorbed solar radiation is computed by using two-band representation:

where SI↓ and SV↓ are the downward fluxes of near-infrared and visible radiation, respectively. Correspondingly, α I and α V are the surface albedos in near-infrared and visible bands. Downward short- and long-wave fluxes near the surface are computed by a spatial horizontal (bilineal) and linear vertical interpolation of the corresponding fluxes computed by the radiative scheme of the climate module.

Sensible heat flux and sublimation are computed following the standard bulk-formulas as

where CM is the heat and moisture exchange coefficient, which depends on atmospheric stratification. ρ a and cp are the near-surface air density and specific heat, respectively, us is the surface wind speed, Ts and Ta are the surface and near-surface air temperatures, respectively, Qa the surface air specific humidity, and Qsat the saturated specific humidity at the temperature Ts.

The snowfall Ps is computed by using the total precipitation interpolated from the climate module via

where
\(\tilde P\) is the horizontally interpolated total precipitation from the climate module, k1= 0.5 m−1 s and k2= 0.4 are empirical parameters. The fraction of total precipitation in solid form (i.e. snowfall) p1 is a continuous function of air temperature:

where Hb= 1,500 m, and z′ is an averaged spatially interpolated elevation used in the climate component of the model.

The last term in the right-hand-side of Eq. 6 represents the slope effect (e.g. Liljequist and Cehak 1974). In previous parameterisations of the slope effect, enhanced precipitation is assumed over any slope irrespective of the wind direction. Here, we take into account the difference between up- and down-wind effects by including the term u·∇zi (Eq. 6), which represents the vertical component of wind arising from the topography. u is the averaged wind vector, and ∇zi is the gradient of surface elevation. Because the climate component of CLIMBER-2 calculates only a climatologically averaged wind vector, synoptic variability of the wind vector is taken into account by adding the term k2usn|∇zi| to Eq. 6, where usn is the norm of the synoptic component of the wind, computed in the climate module (see Petoukhov et al. 2000). This term, unlike the previous one, is always positive, and thus always increases precipitation over the slope.

Albedo of snow surfaces

The albedo of snow surfaces is calculated in two steps. Firstly, the albedo of snow for diffusive radiation is computed for two spectral bands—visible and near-infrared, taking into account the effects of snow aging (growth of snow grains) and contamination of snow by dust as

$$
\alpha _{dk} = \alpha _{nk} - f_a \beta _k - \delta \alpha _k ,
$$

(8)

where the index k=1 refers to visible radiation, and k=2 to the near-infrared band, respectively. Further, αnk is the albedo of new snow. The non-dimensional effective age of snow fa, ranging between 0 and 1, is a function of surface temperature and precipitation. βk is the maximum reduction of snow albedo due to aging of snow, and δαk is an additional reduction of snow albedo due to contamination by dust. The effect of dust deposition on snow albedo is assumed to take the form

where δαnk and δαok are the reductions of snow albedo due to dust for new and old snow, respectively. Both factors are functions of dust concentration on the snow surface. These factors are derived from theoretical calculations by Warren and Wiscombe (1980). The concentration of dust in snow nd is assumed to be equal to the dust concentration in precipitation
\(n_{\text{d}} = \left( {{d \mathord{\left/
{\vphantom {d P}} \right.
\kern-\nulldelimiterspace} P}} \right),\) where d is the dust-deposition rate and P is the total precipitation.

Mass balance and interval of interchange

After the ice-sheet module has simulated the orography and the distribution of ice sheets (see Sect. 2.2), one year of integration of SEMI (with a one-day time step) elapses, and SEMI provides the mass balances for the ice-sheet module. If a grid cell is covered by ice, then the change in mass balance, either positive or negative, is applied to the ice-sheet model in the subsequent year. If there is no ice in a grid cell, but SEMI predicts positive mass balance, then this grid cell in the ice-sheet model is treated as ice-covered in the succeeding year. If ice does not exist in the given grid cell and mass balance is negative, then the grid cell is considered to be ice-free in the succeeding year. A special case is a grid cell with negative mass balance which has a neighbouring grid cell covered by ice: if the absolute value of mass balance exceeds the amount of lateral transport of ice to the given grid cell, then ice does not appear, but melting is taken into account. In these grid points (“hidden ablation” area) SEMI computes only a potential mass balance, while the actual mass balance is calculated as the amount of melted ice, which is equal to the lateral ice inflow into this grid cell. This actual mass balance in the “hidden ablation” area is used for calculation of the total mass balance of the ice sheets.

Calov R, Marsiat I (1998) Simulations of the Northern Hemisphere through the last glacial-interglacial cycle with a vertically integrated and a three-dimensional thermomechanical ice sheet model coupled to a climate model. Ann Glaciol 27:169–176Google Scholar

Le Meur E, Huybrechts P (1996) A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice sheet during the last glacial cycle. Ann Glaciol 23:309–317Google Scholar