Articles | Volume 12, issue 1
Special issue:
ESD Ideas
12 Jan 2021
ESD Ideas |  | 12 Jan 2021

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

Mikhail Y. Verbitsky and Michel Crucifix

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:


Equations (1)–(3) describe dynamical properties of a large continental (e.g., Laurentide) ice sheet. Here S (m2) is the glaciation area, θ (C) is the basal ice sheet temperature, and ω (C) is the global climate temperature. The profile factor ζ (m1∕2) is determined by ice viscosity, density, and acceleration of gravity; a (m s−1) is snow precipitation rate; FS is normalized mid-July insolation at 65 N (Berger and Loutre, 1991); ε (m s−1) is the amplitude of FS; κ (ms-1C-1) and c (ms-1C-1) define ice mass balance sensitivity to ω and θ correspondingly; the adimensional coefficient α is the basal temperature sensitivity to ω changes, β (Cm-2) defines basal temperature sensitivity to the changes in the ice sheet area S, and S0 (m2) is a reference glaciation area; γ1 (C s−1), γ2 (Cm-2s-1) and γ3 (s−1) 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.


Equations (4) and (5) are identical to Eqs. (1) and (2). Here SG (m2) 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 SG and θG is fully determined by five parameters (namely a, ς, S0, β, 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, S0, and β, as parameters with independent dimensions and consider parameters ς and S0 as constant, then, using the generalized π theorem (Sonin, 2004), the time scale P must obey the following:

(6) P = ( a ) - 1 S 0 1 / 2 Ψ ( Π 1 ) ,

where Ψ is a dimensionless function of Π1=a/(βcS0). It can be determined experimentally that the period P is smaller when Π1 is smaller. Following Verbitsky and Chalikov (1986), we now introduce the Peclet number as Pe=a^H/k, where 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 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 Pe-1/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 SG (106km2). 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:

(7) V = 1 β α + κ c γ 2 γ 3 - γ 1 S 0 γ 3 .

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 ε/a=0.3 to ε/a=1.7 and change the reference glaciation area S0 from S0=2×106km2 to S0=12×106km2, thus invoking a hypothetical trend in processes that control long-term CO2 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 Pe-1/2H-1/2S0-1/8 (since HS1/4; Verbitsky et al., 2018), changing from β=2.4C10-6km-2 to β=1.9C10-6km-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 -γ2(SG-S0) 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., S0-ω, βS0-1/8, and the ratio Π1 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 SG. We can see that though the Greenland ice sheet itself has a limited influence on the global climate temperature (γ2SGγ2S), 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 SG millennial variability is 0.1×106km2, corresponding to 0.15×106km3 of ice volume or ∼0.4m 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 (S0, γ1), an increase in climate sensitivity to the ice volume (γ2), 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.

Code and data availability

The MATLAB R2015b code and data to reproduce the results of the numerical experiment as they are presented in Fig. 1c are available at (Verbitsky, 2020).

Author contributions

MYV conceived the research and developed the formalism. MYV and MC contributed equally to writing the paper.

Competing interests

The authors declare that they have no conflict of interest.


We are grateful to our two anonymous reviewers for their helpful comments.

Financial support

Michel Crucifix is funded as Research Director with the Belgian National Fund of Scientific Research.

Review statement

This paper was edited by Anders Levermann and reviewed by two anonymous referees.


Berger, A. and Loutre, M. F.: Insolation values for the climate of the last 10 million years, Quaternary Sci. Rev., 10, 297–317, 1991. 

Hodell, D. A. and Channell, J. E. T.: Mode transitions in Northern Hemisphere glaciation: co-evolution of millennial and orbital variability in Quaternary climate, Clim. Past, 12, 1805–1828,, 2016. 

Li, C. and Born, A.: Coupled atmosphere-ice-ocean dynamics in Dansgaard-Oeschger events, Quaternary Sci. Rev., 203, 1–20,, 2019. 

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records, Paleoceanography, 20, PA1003,, 2005.  

Mitsui, T., Lenoir, G., and Crucifix, M.: Is the glacial climate scale invariant?, Dynamics and Statistics of the Climate System, 3, 1–22,, 2018. 

Pedro J. B., Jochum, M., Buizert, C., He, F., Barker, S., and Rasmussen, S. O.: Beyond the bipolar seesaw: Toward a process understanding of interhemispheric coupling, Quaternary Sci. Rev., 192, 27–46,, 2018. 

Ruddiman, W. F., Raymo, M., and McIntyre, A.: Matuyama 41,000-year cycles: North Atlantic Ocean and northern hemisphere ice sheets, Earth Planet, Sc, Lett,, 80, 117–129,, 1986. 

Saltzman, B.: Three basic problems of paleoclimatic modeling: A personal perspective and review, Clim. Dynam., 5, 67–78, 1990. 

Schulz M., Paul, A., and Timmermann, A.: Relaxation oscillators in concert: a framework for climate change at millennial timescales during the late Pleistocene, Geophys. Res. Lett., 29, 2193,, 2002. 

Sonin, A. A.: A generalization of the Π-theorem and dimensional analysis, P. Natl. Acad. Sci. USA, 101, 8525–8526, 2004. 

Verbitsky, M. Y.: Supplementary code and data to ESD paper “ESD Ideas: The Peclet number is a cornerstone of the orbital and millennial Pleistocene variability”, by: Verbitsky, M. Y. and Crucifix, M.,, 2020. 

Verbitsky, M. Y. and Chalikov, D. V.: Modelling of the Glaciers-Ocean-Atmosphere System, Gidrometeoizdat, Leningrad, 1986. 

Verbitsky, M. Y. and Crucifix, M.: π-theorem generalization of the ice-age theory, Earth Syst. Dynam., 11, 281–289,, 2020. 

Verbitsky, M. Y., Crucifix, M., and Volobuev, D. M.: A theory of Pleistocene glacial rhythmicity, Earth Syst. Dynam., 9, 1025–1043,, 2018. 

Verbitsky, M. Y., Crucifix, M., and Volobuev, D. M.: ESD Ideas: Propagation of high-frequency forcing to ice age dynamics, Earth Syst. Dynam., 10, 257–260,, 2019. 

Short summary
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.
Special issue
Final-revised paper