the Creative Commons Attribution 4.0 License.

Special issue: ESD Ideas

**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

^{1,2}

^{2}

**Mikhail Y. Verbitsky and Michel Crucifix**Mikhail Y. Verbitsky and Michel Crucifix

^{1,2}

^{2}

^{1}Gen5 Group, LLC, Newton, MA, USA^{2}Earth and Life Institute, UCLouvain, Louvain-la-Neuve, Belgium

^{1}Gen5 Group, LLC, Newton, MA, USA^{2}Earth and Life Institute, UCLouvain, Louvain-la-Neuve, Belgium

**Correspondence**: Mikhail Y. Verbitsky (verbitskys@gmail.com)

**Correspondence**: Mikhail Y. Verbitsky (verbitskys@gmail.com)

Received: 08 Jul 2020 – Discussion started: 17 Jul 2020 – Revised: 18 Nov 2020 – Accepted: 27 Nov 2020 – Published: 12 Jan 2021

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.

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.

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* (m^{2}) 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^{−1}) is snow
precipitation rate; *F*_{S} is normalized mid-July insolation at
65^{∘} N (Berger and Loutre, 1991); *ε* (m s^{−1}) is
the amplitude of *F*_{S}; *κ*
($\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{}^{\circ}{\mathrm{C}}^{-\mathrm{1}}$) and *c*
($\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{}^{\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}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}$) defines basal temperature sensitivity to the
changes in the ice sheet area *S*, and *S*_{0} (m^{2}) is a reference
glaciation area; *γ*_{1} (^{∘}C s^{−1}), *γ*_{2}
(${}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{-\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{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 *S*_{G} (m^{2}) 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:

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 Π_{1} 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.

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:

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*_{0} from ${S}_{\mathrm{0}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ km^{2} to ${S}_{\mathrm{0}}=\mathrm{12}\times {\mathrm{10}}^{\mathrm{6}}$ km^{2}, thus invoking a hypothetical trend in processes that
control long-term CO_{2} 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}\phantom{\rule{0.125em}{0ex}}{\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\mathrm{km}}^{-\mathrm{2}}$ to *β*=1.9 ${}^{\circ}\mathrm{C}\phantom{\rule{0.125em}{0ex}}{\mathrm{10}}^{-\mathrm{6}}\phantom{\rule{0.125em}{0ex}}{\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
Π_{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
*S*_{G}. We can see that though the Greenland ice sheet itself has a
limited influence on the global climate temperature (*γ*_{2}*S*_{G}≪*γ*_{2}*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^{2}, corresponding to $\sim \mathrm{0.15}\times {\mathrm{10}}^{\mathrm{6}}$ km^{3} 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.

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*_{0}, *γ*_{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.

The MATLAB R2015b code and data to reproduce the results of the numerical experiment as they are presented in Fig. 1c are available at https://doi.org/10.5281/zenodo.4292357 (Verbitsky, 2020).

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

The authors declare that they have no conflict of interest.

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

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

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, https://doi.org/10.5194/cp-12-1805-2016, 2016.

Li, C. and Born, A.: Coupled atmosphere-ice-ocean dynamics in Dansgaard-Oeschger events, Quaternary Sci. Rev., 203, 1–20, https://doi.org/10.1016/j.quascirev.2018.10.031, 2019.

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57
globally distributed benthic *δ*^{18}O records, Paleoceanography, 20, PA1003, https://doi.org/10.1029/2004PA001071, 2005.

Mitsui, T., Lenoir, G., and Crucifix, M.: Is the glacial climate scale invariant?, Dynamics and Statistics of the Climate System, 3, 1–22, https://doi.org/10.1093/climsys/dzy011, 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, https://doi.org/10.1016/j.quascirev.2018.05.005, 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, https://doi.org/10.1016/0012-821X(86)90024-5, 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, https://doi.org/10.1029/2002GL016144, 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., https://doi.org/10.5281/zenodo.4292357, 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, https://doi.org/10.5194/esd-11-281-2020, 2020.

Verbitsky, M. Y., Crucifix, M., and Volobuev, D. M.: A theory of Pleistocene glacial rhythmicity, Earth Syst. Dynam., 9, 1025–1043, https://doi.org/10.5194/esd-9-1025-2018, 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, https://doi.org/10.5194/esd-10-257-2019, 2019.