ESD Ideas: The Peclet number is a cornerstone of the orbital and millennial Pleistocene variability

We demonstrate here that a single physical phenomenon, specifically, a naturally changing balance between intensities of temperature advection and diffusion in the viscous ice media, may influence the entire spectrum of the Pleistocene variability from orbital to millennial timescales.

1 Introduction

About 1 Myr ago, the dynamics of glacial–interglacial cycles experienced a major change, transitioning from the predominant 40 kyr periodicity to approximately 100 kyr variability (e.g., Ruddiman et al., 1986; Lisiecki and Raymo, 2005). This change in the glaciation rhythmicity is called the Middle Pleistocene transition (MPT). On the other hand, Hodell and Channell (2016) presented the analysis of 3.2 Myr records of stable isotopes and sediment physical properties and observed coherence between millennial and orbital-scale variabilities: after the MPT, millennial events are more frequent and more intense than before. In our previous work (Verbitsky et al., 2019), we have shown that millennial variability may penetrate upscale and influence the slower pace dynamics of glacial–interglacial cycles. Here, we will demonstrate that the same ice sheet non-linearity that is responsible for the penetration of the millennial oscillations into the orbital domain may also contribute – without invoking any additional physics – to the millennial variability and explain its coherence with orbital-scale dynamics. To underline the significance of this proposal, we briefly summarize the current understanding of millennial variability.

Classically, millennial variability during the glacial periods is explained by reference to two families of mechanisms acting in concert. The first one relates to iceberg calving, which occurs when ice sheets have reached the ocean margins and ice shelves developed. Episodes of intense ice rafting generate a halocline and reduce ocean ventilation (the Heinrich events); on the other hand, changes in ocean temperature and sea-level influence the stability of ice shelves and may pre-condition calving events. The second family of mechanisms involves the dynamics of ocean circulation, deep-water mixing, and sea-ice formation at high latitudes. The ocean circulation effectively causes horizontal and vertical transports of buoyancy, which have a role in maintaining the circulation itself. Because of the non-linear character of this feedback, the geographical structure of convection in the North Atlantic realm can shift very rapidly (within years) between different configurations. There is converging evidence that the Dansgaard–Oeschger events identified in Greenland correspond to such circulation changes. Millennial variability associated with all these mechanisms is most likely to occur during cold periods (ice sheets must be large enough) but not too cold (sea ice must not be locked into a deep freeze): a so-called “sweet spot” (Mitsui et al., 2018; Pedro et al., 2018; Li and Born, 2019) during which iceberg calving and ocean-circulation changes can entertain a rich pattern of oscillations (Schulz et al., 2002). The mechanism that we suggest here is not subject to the same sweet spot and can actually explain a pervasive variability during both glacial and interglacial periods.

2 Methods

We consider the non-linear dynamical model of the global climate system presented in Verbitsky et al. (2018). It is derived from the scaled conservation equations of the non-Newtonian ice flow and combined with a linear feedback equation of the climate temperature:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{4}}{\mathrm{5}}}{\mathit{\zeta }}^{-\mathrm{1}}{S}^{\mathrm{3}/\mathrm{4}}(a-\mathit{\epsilon }{F}_{\mathrm{S}}-\mathit{\kappa }\mathit{\omega }-c\mathit{\theta }),\text{(2)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}t}}={\mathit{\zeta }}^{-\mathrm{1}}{S}^{-\mathrm{1}/\mathrm{4}}(a-\mathit{\epsilon }{F}_{\mathrm{S}}-\mathit{\kappa }\mathit{\omega })\mathit{\{}\mathit{\alpha }\mathit{\omega }+\mathit{\beta }[S-S_{\mathrm{0}}]-\mathit{\theta }\mathit{\}},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}\mathit{\omega }}{\mathrm{d}t}}={\mathit{\gamma }}_{\mathrm{1}}-{\mathit{\gamma }}_{\mathrm{2}}[S-S_{\mathrm{0}}]-{\mathit{\gamma }}_{\mathrm{3}}\mathit{\omega }.\end{array}$$

Equations (1)–(3) describe dynamical properties of a large continental (e.g., Laurentide) ice sheet. Here S (m²) is the glaciation area, θ (^∘C) is the basal ice sheet temperature, and ω (^∘C) is the global climate temperature. The profile factor ζ (m^1∕2) is determined by ice viscosity, density, and acceleration of gravity; a (m s⁻¹) is snow precipitation rate; F_S is normalized mid-July insolation at 65^∘ N (Berger and Loutre, 1991); ε (m s⁻¹) is the amplitude of F_S; κ ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$) and c ($\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{}^{\circ }{\mathrm{C}}^{-\mathrm{1}}$) define ice mass balance sensitivity to ω and θ correspondingly; the adimensional coefficient α is the basal temperature sensitivity to ω changes, β (${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$) defines basal temperature sensitivity to the changes in the ice sheet area S, and S₀ (m²) is a reference glaciation area; γ₁ (^∘C s⁻¹), γ₂ (${}^{\circ }\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) and γ₃ (s⁻¹) define the relaxation dynamics of ω.

We now supplement this system with Eqs. (4) and (5) that describe the thermodynamical properties of a smaller-size (e.g., Greenland) ice sheet.

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{G}}}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{4}}{\mathrm{5}}}{\mathit{\zeta }}^{-\mathrm{1}}{S}_{\mathrm{G}}^{\mathrm{3}/\mathrm{4}}\left(a^{\prime }-{\mathit{\epsilon }}^{\prime }{F}_{\mathrm{S}}-{\mathit{\kappa }}^{\prime }\mathit{\omega }-c^{\prime }{\mathit{\theta }}_{\mathrm{G}}\right),{\displaystyle }& {\displaystyle \frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{G}}}{\mathrm{d}t}}={\mathit{\zeta }}^{-\mathrm{1}}{S}_{\mathrm{G}}^{-\mathrm{1}/\mathrm{4}}\left(a^{\prime }-{\mathit{\epsilon }}^{\prime }{F}_{\mathrm{S}}-{\mathit{\kappa }}^{\prime }\mathit{\omega }\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }\times \mathit{\{}{\mathit{\alpha }}^{\prime }\mathit{\omega }+{\mathit{\beta }}^{\prime }(S_{\mathrm{G}}-{S^{\prime }}_{\mathrm{0}})-{\mathit{\theta }}_{\mathrm{G}}\mathit{\}}.\end{array}$$

Equations (4) and (5) are identical to Eqs. (1) and (2). Here S_G (m²) and θ_G (^∘C) are the area and the basal temperature of the Greenland ice sheet; all other parameters have the same dimensions and meaning as the corresponding parameters in Eqs. (1) and (2). Some (but not all) of them may have different numerical values, and for this reason we mark them with an apostrophe. The area of the Greenland ice sheet is limited by the size of the island and its variations have limited impacts on the global climate. Therefore, for reasoning on scaling relationships, we can neglect its contribution to the ω-evolution and leave Eq. (3) unchanged. Hence, the Greenland ice sheet model may be reduced to a system formed by Eqs. (4) and (5), forced by the astronomical forcing and the global climate temperature ω. The dynamical system formed by Eqs. (1)–(5) is represented graphically in Fig. 1a. Without external forcing, the evolution of S_G and θ_G is fully determined by five parameters (namely a^′, ς, $S_{\mathrm{0}}^{\prime }$, β^′, and c^′). Some combinations of these parameters generate a damped-oscillation behavior, and, as we show now, it is reasonably straightforward to estimate the period P of this oscillation. Indeed, if we take parameters a^′, $S_{\mathrm{0}}^{\prime }$, and β^′, as parameters with independent dimensions and consider parameters ς and $S_{\mathrm{0}}^{\prime }$ as constant, then, using the generalized π theorem (Sonin, 2004), the time scale P must obey the following:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }P=\left(a^{\prime }{)}^{-\mathrm{1}}{S^{\prime }}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}\mathrm{\Psi }\right({\mathrm{\Pi }}_{\mathrm{1}}),\end{array}$$

where Ψ is a dimensionless function of ${\mathrm{\Pi }}_{\mathrm{1}}=a^{\prime }/({\mathit{\beta }}^{\prime }{c}^{\prime }{S^{\prime }}_{\mathrm{0}}$). It can be determined experimentally that the period P is smaller when Π₁ is smaller. Following Verbitsky and Chalikov (1986), we now introduce the Peclet number as $Pe=\widehat{a}H/k$, where $\widehat{a}$ is a characteristic mass influx, i.e., accumulation minus ablation, H is ice thickness, and k is temperature diffusivity. The Peclet number measures the balance between temperature advection and diffusion. Since the Greenland ice sheet is relatively thin, for comparable $\widehat{a}$, its Peclet number is smaller than that of a big ice sheet, and therefore the effect of the geothermal heat flux on the basal temperature is, in relative terms, stronger. Verbitsky et al. (2018) showed that the parameter β, which determines the basal temperature response to the changes in the ice sheet size, emerges as a delicate balance between vertical advection, internal friction, and geothermal heat flux, and that it is proportional to $P{e}^{-\mathrm{1}/\mathrm{2}}$. Thus, the relatively small size of the Greenland ice sheet implies a higher β^′ and, according to Eq. (6), its damped oscillations have a higher frequency than, say, the Laurentide ice sheet. In fact, numerical experiments with Eqs. (4)–(5) and with a reasonable set of parameters produce millennial-range oscillations of the Greenland ice sheet.

Figure 1 (a) The dynamical system (Eqs. 1–5) as it has been used for scaling reasoning. The Greenland ice sheet is coupled one-way via the climate temperature ω acting as an external forcing for it. Red circles mark positive feedback loops and green circles mark negative feedback loops; (b) the dynamical system (Eqs. 1–5) as it has been used in the numerical simulation. The Laurentide and Greenland ice sheets are fully coupled. Pink circles mark weaker positive feedback loops. (c) Evolution of wavelet spectra over the past 3 Myr for the Greenland glaciation area S_G (10⁶ km²). The color scale shows the continuous Morlet wavelet amplitude, the thick line indicates the peaks with 95 % confidence, and the shaded area indicates the cone of influence for wavelet transform.

With these considerations in hand, let us foresee the consequences of the MPT. As the Laurentide ice sheet reached gradually increasing footprints on the large American continent, its Peclet number increased, implying that temperature advection became the dominant process in the moving ice media. Consequently, β became smaller. We previously showed (Verbitsky et al., 2018; Verbitsky and Crucifix, 2020) that the dynamics of the system described by Eqs. (1)–(3) are largely determined by a V number measuring the balance between positive and negative model feedbacks:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }V={\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}\left(\mathit{\alpha }+{\displaystyle \frac{\mathit{\kappa }}{c}}\right)\left({\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{2}}}{{\mathit{\gamma }}_{\mathrm{3}}}}-{\displaystyle \frac{{\mathit{\gamma }}_{\mathrm{1}}}{S_{\mathrm{0}}{\mathit{\gamma }}_{\mathrm{3}}}}\right).\end{array}$$

Thus, high values of the parameter β (V∼0.1) imply a weak positive feedback in the system, and, conversely, low values of the parameter β (V∼0.75) imply a strong positive feedback. As has been previously shown (Verbitsky et al., 2018), when the orbital forcing is large enough compared to the average ice accumulation (high enough ε∕a ratio), an increase in the V number generates a period-doubling bifurcation, that is, an escape from the 40 kyr oscillation regime towards longer glacial–interglacial cycles. The increased period of the system response to the astronomical forcing commands the glaciation area S to become larger, meaning that the amplitude of the oscillation is larger after the MPT than before. This amplitude increase may be understood as a consequence of a scale-invariance property (Verbitsky and Crucifix, 2020). In this case, we see that the natural evolution of the Peclet number generates higher-amplitude oscillations of the larger ice sheets (i.e., Laurentide ice sheet), along with longer glacial–interglacial cycles. On the other hand, for the Greenland ice sheet, the Peclet number cannot grow. It remains bounded by a small value which defines millennial-scale oscillations. The large post-MPT oscillations of the Laurentide ice sheet may then excite more vigorously the millennial oscillations of Greenland. The idea of the present contribution is therefore the following: the Peclet number is a key quantity for explaining the joint emergence of large-amplitude longer-period oscillations of the bigger ice sheets, along with the increase in the occurrence of the millennial oscillations of the smaller-size ice sheets.

To test our theoretical perspective, we reproduce the Pleistocene glaciation history, solving the system (Eqs. 1–5) jointly and thus calculating a coevolution of the Laurentide and Greenland ice sheets. For the Laurentide ice sheet, we change linearly the ε∕a ratio over the last 3 Myr from $\mathit{\epsilon }/a=\mathrm{0.3}$ to $\mathit{\epsilon }/a=\mathrm{1.7}$ and change the reference glaciation area S₀ from $S_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ km² to $S_{\mathrm{0}}=\mathrm{12}\times {\mathrm{10}}^{\mathrm{6}}$ km², thus invoking a hypothetical trend in processes that control long-term CO₂ levels. The ε∕a ratio is reduced by changing the a component only, following the assumption that the Pleistocene cooling trend goes along with a decrease in the snow precipitation rate. The parameter β is modified as being proportional to $P{e}^{-\mathrm{1}/\mathrm{2}}\sim H^{-\mathrm{1}/\mathrm{2}}\sim S_{\mathrm{0}}^{-\mathrm{1}/\mathrm{8}}$ (since $H\sim S^{\mathrm{1}/\mathrm{4}}$; Verbitsky et al., 2018), changing from β=2.4 ${}^{\circ }\mathrm{C}{\mathrm{10}}^{-\mathrm{6}}{\mathrm{km}}^{-\mathrm{2}}$ to β=1.9 ${}^{\circ }\mathrm{C}{\mathrm{10}}^{-\mathrm{6}}{\mathrm{km}}^{-\mathrm{2}}$. To account for the relatively small effect of Greenland changes on climate, which we have neglected for our scaling reasoning, we add a term $-{\mathit{\gamma }}_{\mathrm{2}}(S_{\mathrm{G}}-{S^{\prime }}_{\mathrm{0}})$ to the right side of Eq. (3). The dynamical system (Eqs. 1–5) as it has been used in the numerical experiment is shown in Fig. 1b. The Laurentide and Greenland ice sheets are fully coupled here. For the Greenland ice sheet, we assume that the reference glaciation area depends on climate temperature, i.e., $S_{\mathrm{0}}^{\prime }\sim -\mathit{\omega }$, ${\mathit{\beta }}^{\prime }\sim {S^{\prime }}_{\mathrm{0}}^{-\mathrm{1}/\mathrm{8}}$, and the ratio Π₁ is an order of magnitude less than the corresponding ratio of the Laurentide ice sheet. In Fig. 1c we present the results as a wavelet spectrum over the past 3 Myr for the Greenland glaciation area S_G. We can see that though the Greenland ice sheet itself has a limited influence on the global climate temperature (γ₂S_G≪γ₂S), the changes in Greenland's geometry reflect all major events of the Pleistocene history – a transition from the double-precession 40 kyr variability to increased amplitudes of double-obliquity oscillations. The global climate temperature ω, driven by the Laurentide ice sheet (taken here as a proxy for all the large ice sheets of the Northern Hemisphere), shifts the equilibrium state of the Greenland ice sheet, and the latter responds with millennial-period damped oscillatory adjustments. Such shifts are more prominent in the Late Pleistocene, and, therefore, consistently with Hodell and Channell (2016), the millennial events of the Late Pleistocene are more frequent and more intense. The amplitude of S_G millennial variability is $\sim \mathrm{0.1}\times {\mathrm{10}}^{\mathrm{6}}$ km², corresponding to $\sim \mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}$ km³ of ice volume or ∼0.4 m of the sea-level change. Talking about the sea level, we must advise our readers that the coupling of the Greenland and Laurentide ice sheets has been done here via global albedo and global temperature only. In other situations (e.g., consider the Antarctic ice sheet) ice sheet volume changes may cause stronger sea-level variations, which are instantly distributed worldwide and would immediately affect the Laurentide ice sheet.

3 Discussion

Until recently, the timescale separation approach (Saltzman, 1990) has been implicitly or explicitly used to consider separately glacial–interglacial cycles and millennial variability. Accordingly, it is generally accepted that the spectrum of Pleistocene variability in the orbital domain is defined by continental ice sheets, while the millennial part of the spectrum is to be attributed to ocean–atmosphere interactions. We first challenged this approach by demonstrating that ice sheet non-linearity allows millennial variability to propagate upscale and influence ice-age dynamics (Verbitsky et al., 2019). Here, we show that the same non-linear ice-flow dynamics can also affect the millennial part of the spectrum.

From Eq. (7), it appears that different mechanisms would explain an increase in the V number throughout the Pleistocene. Therefore, several scenarios can be invoked as the origin of the MPT. For example, a steady decrease in the background climate temperature (S₀, γ₁), an increase in climate sensitivity to the ice volume (γ₂), a change in the sensitivity of ice mass balance and its temperature to global climate temperature (κ and α), or even a decrease in the intensity of ice sliding (c). At the same time, all these mechanisms produce the changes in the Peclet number which we have described. For this reason, we claim that the Peclet number, changing in concert with the glaciation size, is the key similarity parameter connecting the millennial and orbital Pleistocene variations. We find it remarkable that the same physical phenomenon can explain all major events of the Pleistocene in the orbital domain if the evolution of an ice sheet is not restricted and, in addition, define variability in the millennial domain when the size of the ice sheet is bounded.