# Simple Stochastic Modeling of Snowball Probability Throughout Earth History

###### Abstract

Over its multibillion-year history, Earth has exhibited a wide range of
climates. Its history ranges from snowball episodes where the surface was
mostly or entirely covered by ice to periods much warmer than today, where the
cryosphere was virtually absent. Our understanding of greenhouse gas evolution
over this long history, specifically carbon dioxide, is mainly informed by
deterministic models. However, the complexity of the carbon cycle and its
uncertainty over time motivates the study of non-deterministic models, where
key elements of the cycle are represented by inherently stochastic processes.
By doing so, we can learn what models of variability are compatible with
Earth’s climate record instead of how exactly this variability is produced. In
this study, we address why there were snowballs in the Proterozoic, but not
the Phanerozoic by discussing two simple stochastic models of long-term
carbon-cycle variability. The first, which is the most simple and represents
CO2 concentration directly as a stochastic process, is instructive and perhaps
intuitive, but is incompatible with the absence of snowballs in the
Phanerozoic. The second, which separates carbon source from sink and
represents CO2 outgassing as a stochastic process instead of concentration, is
more flexible. When outgassing fluctuates over longer periods, as opposed to
brief and rapid excursions from a mean state, this model is more compatible
with the snowball record, showing only modest increases in the probability of
snowball events over Earth history. The contrast between these models
illustrates what kind of variability may have characterized the long-term
carbon cycle.

Geochemistry, Geophysics, Geosystems

Department of Earth and Planetary Sciences, Harvard University Department of
Earth and Planetary Sciences, Yale University

Mark<EMAIL_ADDRESS>
We explore models of randomness in the long-term climate and carbon cycle

Symmetric noise in the CO2 concentration is instructive, but incompatible with
the snowball record

Modeling noise in the CO2 outgassing rate, instead, is more flexible and can
be qualitatively reconciled with the snowballs

## Plain Language Summary

The Earth’s long-term climate record is complicated and variable. At some
points, our planet has been extremely cold and almost entirely covered in ice
(snowball periods). At others, the climate was much warmer than today and the
surface was ice-free. These changes may have many overlapping causes and are
quite complex. In light of this, we use simple models to understand how
randomness in the carbon cycle affects the long-term climate. In one model,
the atmospheric carbon dioxide concentration is treated as a randomly
fluctuating quantity. We can learn from this model, but it is at odds with the
timing of Earth’s snowball periods. In a second model, the CO2 _outgassing
rate_ is represented by a randomly fluctuating quantity instead of the
concentration itself, and is more consistent with the geological record. These
models offer some qualitative lessons about the kind of randomness that is
compatible with Earth’s snowball episodes.

##  1 Introduction

Earth’s climate has varied dramatically over the past 2.5 billion years,
ranging from “hothouse” climates such as the Eocene (55.8 – 33.9 Ma) [Berner
(1990), Greenwood Wing (1995)] to “Snowball Earth” episodes, when virtually
all of the Earth’s surface was encased in ice [Hoffman Schrag (2002),
Pierrehumbert . (2011)]. Geological evidence points to at least three snowball
episodes, with the first event occurring at the beginning of the
Paleoproterozoic era, at approximately 2.5 Ga [Evans . (1997), Kirschvink .
(2000)]. The next two occurred during the Cryogenian period (720 – 635 Ma),
with inception times about 50 Myr apart [Hoffman . (1998), Rooney . (2015),
Prave . (2016)].

What caused snowball climates at some points, but warm climates at others?
Over geologic timescales, the atmospheric carbon dioxide (CO2 ) concentration,
and therefore global temperature, is set by the balance between sources
(volcanic outgassing) and sinks (silicate weathering) [Urey (1952), Walker .
(1981), Marshall . (1988), Berner Lasaga (1989)]. Silicate weathering and
marine carbonate burial sequester CO2 from the ocean and atmosphere until it
is recycled back by the subduction of oceanic crust [Berner . (1983)]. Indeed,
variations in volcanic outgassing and availability of weatherable materials
have both been suggested to be dominant controls on Earth’s long term climate
evolution [<]e.g.,¿McKenzie-Horton-Loomis-et-
al-2016:continental,macdonald_arc-continent_2019.

Earth’s temperature and CO2 concentration are not well known deep into the
planet’s history. Although the Sun’s luminosity has slowly increased since the
birth of the solar system, there is no evidence for secular warming of the
climate [Feulner (2012)]. To compensate for the increase in absorbed solar
radiation, the mean CO2 concentration has very likely decreased over time
[Halevy Bachan (2017), Krissansen-Totton . (2018)]. This decrease has not been
smooth, but rather accompanied by significant variability, the cause of which
may be multifactorial and is still an active area of research [Franks .
(2014), Montañez . (2016), Lenardic . (2016), Macdonald . (2019), Park .
(2021), Baum . (2022)].

Past studies have generally attempted to understand variability in atmospheric
CO2 concentration in a deterministic fashion, by investigating specific
components of the carbon cycle and their potential effect on climate. For
example, we can attempt to constrain the independent effects of continental
configuration, topography, biochemistry, ocean circulation, and a wide range
of other processes. Here we take a different approach, investigating the
effects of explicitly random or stochastic processes in the long term carbon
cycle. This approach views the climate system as an inherently stochastic,
complex system and attempts to understand what models of randomness agree with
our sparse and incomplete record of Earth’s climate on the longest timescales.
Specifically, we anchor our modeling and discussion to the aforementioned
Snowball Earth episodes, as they represent very important excursions from the
mean state. Our conception of long-term climate variability must be compatible
with the observed timing of snowball events.

That is not to say that stochastic models are absent from the climate
literature. There is a long history of stochastic climate modeling touching
virtually all components of the Earth system [Hasselmann (1976), Imkeller Von
Storch (2001), Imkeller Monahan (2002), Farrell Abbot (2012), Franzke .
(2015)]. Some recent studies have used stochastic methods to model Earth’s
carbon cycle [<]e.g.,¿Zeebe-Westerhold-Littler-et-
al-2017:orbital,wordsworth_how_2021,Arnscheidt-Rothman-2021:asymmetry. In this
paper, we take a similar approach by setting the complexity of the three-
dimensional, chaotic climate system aside to examine the consequences of two
simple stochastic models of the carbon cycle over the past $\sim$2 Gyr. We do
not suggest that our modeling is a precise representation of Earth’s long-term
carbon record, but use it to understand what kinds of variability (what
models) are compatible with the snowball record.

At the beginning of Section 2, we explain the primary physical equations and
assumptions governing both of these models. In Section 2.1, we interrogate the
stochastic model recently published by wordsworth_how_2021 and discuss its
compatibility with Earth’s record of snowball periods. In 2.2 we describe an
alternative model with a deterministic weathering sink but stochastic CO2
outgassing. In Section 3 we explain how we solve the model equations
numerically and produce a large ensemble of these solutions. In Section 4 we
present the results of these simulations, compare with the model of Section
2.1, and discuss compatibility with Earth’s snowball record. Finally, in
Section 5, we offer some concluding remarks and discuss possible avenues of
future research.

##  2 Simple Stochastic Climate Modeling

We use a simple, zero-dimensional model of stochastic, long-term climate
variability representing the long-term evolution Earth’s carbon cycle. The
time-evolution of Earth’s temperature is governed by the balance between
outgoing longwave radiation (OLR) and absorbed solar radiation,

$C\dot{T}=\frac{F(t)}{4}\left[1-\alpha\right]-\textrm{OLR}(T,f)\,,$ (1)

where $C$ is the system’s heat capacity, $T$ is the global mean surface
temperature, $t$ is time, $F$ is the time-dependent incident solar flux,
$\alpha$ is the planetary albedo, OLR is the temperature-dependent and
CO2-dependent outgoing longwave radiation, and $f$ is the atmospheric CO2
concentration in ppm. We set the albedo to a constant 30 %.

In reality, due to the response of clouds and ice to surface temperature,
albedo is a complex function of temperature. However, the goal of this model
is to simulate the likelihood of different temperatures over time, not
represent the full-complexity of the ice-albedo feedback. Instead, we assume a
hypothetical snowball threshold temperature of $T_{\textrm{snow}}=280$ K
[Pierrehumbert . (2011)] where the ice-albedo feedback runs away. The precise
value of this threshold does not affect our conclusions.

The solar flux increases over time according to the standard approximation
[Gough (1981)],

$F(t)=F_{0}\left[1+\frac{2}{5}\left(1-t/t_{0}\right)\right]^{-1}\,,$ (2)

where $F_{0}$ is the modern value of 1366 W/m2 and $t_{0}$ is 4.5 Gyr.

The heat capacity of the ocean-atmosphere system is such that the temperature
approaches equilibrium over $\sim$103 yr, much more rapidly than the 106-109
yr timescales of interest in our simulations. Accordingly, we assume
instantaneous temperature equilibrium and set $\dot{T}$ to zero. We also
linearize OLR around Earth’s preindustrial surface temperature and CO2
concentration of $T_{e}=288$ K and $f_{0}=285$ ppm, respectively.

$\textrm{OLR}=\textrm{OLR}_{0}+a\left(T-T_{e}\right)-b\log\left(f/f_{0}\right)\,,$
(3)

where $a=2$ W/m2/K [Abbot (2016)] and $b=5.35$ W/m2 [Myhre . (1998)]. To
balance Equation 1 at preindustrial conditions and $t=4.5$ Gyr, OLR0 must be
equal to $F_{0}(1-A)/4$. These assumptions, plugged into Equation 1, yield

$0=\frac{F(t)}{4}\left[1-\alpha\right]-\textrm{OLR}_{0}-a\left(T-T_{e}\right)+b\log\left(f/f_{0}\right)$
(4)

which can be solved directly for temperature as a function of time and CO2
concentration,

$T(t,f)=T_{e}+\frac{1}{a}\left[\frac{F(t)}{4}\left(1-\alpha\right)-\textrm{OLR}_{0}+b\log\left(f/f_{0}\right)\right]\,.$
(5)

This is the governing equation of all our subsequent modeling. We note that
the linearization in Equation 3 is a fairly good approximation when $T$ is
within 210–310 K [Koll Cronin (2018)].

###  2.1 Stochastic CO2 Concentration

Equation 5 defines temperature in terms of time $t$ (capturing the brightening
sun) and the carbon dioxide concentration $f$. To introduce randomness into
this system, one approach is to model $f$ as a stochastic process. We analyze
and discuss this approach here, along with its primary assumptions and
implications, before continuing with the development of our new model
equations.

Recently, wordsworth_how_2021 used an Ornstein-Uhlenbeck process [Jacobs
(2010), Dobrow (2016)] to represent randomness in $f$. This approach models
$f$ with Gaussian noise that relaxes to a prescribed value $\chi$. They used a
nondimensional representation, $y=f/f_{0}$, and the model equation is

$dy=\frac{1}{\tau}\left[y-\frac{\chi(t)}{f_{0}}\right]dt+\sigma dB\,,$ (6)

where $\tau$ is the relaxation timescale, $\chi$ is the time-dependent value
of $f$ that would achieve perfect temperature equilibrium (at $T_{e}$),
$\sigma$ scales the noise, and $dB$ is a Wiener process (Gaussian or Brownian
noise). In this case, $\tau$ is analogous to the timescale required for
chemical weathering to restore temperature equilibrium.

In their model, randomness occurs _directly_ in the CO2 concentration and $f$
is the sole prognostic variable. The stochastic differential equation above is
integrated numerically over 108–109 yr and the resulting sequence of discrete
$f$ values determines the simulated temperature history via Equation 5. As
time proceeds, $f$ varies randomly but is restored to equilibrium over a
timescale of $\tau\approx 3$ Ma, which is short compared to the integration
timescale. Figure 1 shows one example simulation using this equation, where
$\sigma$ was chosen for a reasonable chance that the temperature would reach
$T_{\textrm{snow}}$.

Figure 1: A single realization of the stochastic CO2 concentration model
described by Equations 5 and 6, over 250 Myr, with ${\tau=3}$ Myr and
${\sigma=10}$. Panels (a) and (b) show the CO2 concentration over time, with
the equilibrium value $\chi$ drawn in gray. Panel (a) shows the entire time
series and panel (b) shows the same simulation, but zoomed into the $\sim$5
Myr surrounding the temperature minimum. Panels (c) and (d) show temperature
evolution for the same simulation, again zoomed in for panel (d), with
temperatures below 280 K highlighted in blue. Panels (b) and (d), in
particular, highlight the exquisite temperature sensitivity when CO2
concentrations are low. A relatively small change in $f$ near the minimum
between 13 and 11 Ma produces a very dramatic drop in temperature. Note that
sharp negative excursions are a consequence of the logarithmic dependence of
temperature on CO2, and not an ice-albedo feedback.

Relaxation of $f$ toward $\chi$ happens considerably faster than the slow,
insolation driven drift in $\chi$ itself. As a result, over many random
realizations of this model, the CO2 concentration is approximately normally
distributed around $\chi$ at all points in time.

$f\sim\mathcal{N}(\chi(t),\sigma^{2})\,,$ (7)

where $\mathcal{N}$ is the normal distribution with mean of $\chi(t)$ and
variance of $\sigma^{2}$. It is relatively straightforward to examine the
resulting temperature probabilities over time, using Equation 5.

The top row of Figure 2 shows distributions in the CO2 concentration $f$ at
three different times over the most recent 500 Myr. Earlier in time, the sun
is fainter so a higher $\chi$ value is required for a mean temperature of 288
K. At later times, the brighter sun results in lower $\chi$ values, shifting
the distribution down. The bottom row of this figure shows the same
distributions, but plotted against the corresponding temperatures. Because of
the logarithmic relationship between $T$ and $f$, the probable temperature
space widens as time proceeds and cold temperatures become significantly more
likely than warm ones. This happens because, for a single value of $\sigma$,
the variation in $f$ represents a larger fractional change when the mean is
lower. In this model, snowball temperatures are dramatically more likely at
later times. This phenomena was also recently discussed by Graham-2021:high.

Figure 2: In the top row, normal distributions in the atmospheric CO2
concentration ($f$), for several $\sigma$ values, at three different snapshots
in time. The $\chi$ values are the distribution means required to produce a
mean temperature of ${T_{e}=288}$ K using Equation 5. The bottom row shows the
exact same distributions, but plotted in temperature space instead. These
plots show how, even though the mean temperature does not change (by design),
the space of probable temperatures widens over time and cold temperatures
become strongly favored. The logarithmic relationship between $T$ and $f$
produces this effect. The hypothetical snowball temperature of 280 K becomes
dramatically more likely at later times.

Figure 3 shows exactly how much more likely cold excursions become at later
times in this model: orders of magnitude. The left panel shows the cumulative
probability of the CO2 concentration required to achieve a hypothetical
snowball temperature of 280 K, for several values of $\sigma$, over time. Note
the logarithmic vertical axis. Beyond the most recent periods of time, a
snowball is effectively impossible. The right panel shows the temperature
corresponding to the first percentile of the $f$ distributions over time,
showing a precipitous drop in this temperature as time proceeds.

Figure 3: On the left, the cumulative probability of snowball temperatures for
different values of $\sigma$ over time. Each line represents the probability
that random variation in $f$, as defined by Equation 6, would produce $T\leq
280$ K. As expected from Figure 2, the snowball probability plummets at
earlier times. On the right, the temperatures corresponding to the first
percentile of the $f$ distributions. Although the first percentile is an
arbitrary metric, it illustrates the very strong bias against cold
temperatures at earlier times.

To summarize, this model represents randomness in the climate system by making
the CO2 concentration approximately equivalent to a normally distributed
random variable with a mean value that drifts over time to maintain stable
temperature. Outside the most recent period, stochastically induced cold
temperatures are extremely low probability—effectively prohibited. For higher
values of $\sigma$, as wordsworth_how_2021 showed, the same transition occurs
as time proceeds, but it simply occurs at earlier times.

However, it is worth considering whether randomness in the climate system is
well represented by Gaussian noise in the CO2 concentration directly. First,
the principal results of this model strongly contradict Earth’s geologic
record, which indicates snowball episodes occurred roughly 2.5 Ga and 700 Ma.
To be fair, wordsworth_how_2021 focuses mainly on application to exoplanets,
but in that case assuming Earth-like weathering sinks and Sun-like stellar
properties considerably limits the model’s applicability.

Second, it is questionable to assume that CO2 concentrations are well
represented by symmetric noise. To first order, atmospheric CO2 is controlled
over long periods by volcanic outgassing (source) and chemical weathering
(sink). It is intuitive to assume that volcanic outgassing is noisy, as it is
comprised of eruptions of various magnitudes and timescales. This justifies a
noisy, non-negative, _source_ of carbon dioxide. While the weatherability of
Earth’s crust has probably varied considerably over Earth history [Macdonald .
(2019)], these fluctuations occur on longer timescales and are unlikely to
resemble the character of noisy volcanic outgassing.

This model yields conceptual insights. It demonstrates the basic point that
any change in atmospheric CO2 concentration produces a larger temperature
response when the initial concentration is lower. This response is strongly
skewed toward cold temperatures because of the temperature’s logarithmic
dependence on CO2. Figure 2 demonstrates these effects. However, this model
contradicts Earth’s snowball record and may be of limited relevance to
exoplanets because of its strong assumptions about Earth-like weathering
behavior and boundary conditions. It also makes no explicit distinction
between weathering and outgassing, limiting its physical interpretability.

###  2.2 Stochastic CO2 Outgassing

We develop a model with a more physically consistent representation of long-
term CO2 outgassing and removal. In our model, the temperature is still driven
by CO2 concentration and the secular brightening of the sun defined by
Equation 5. However, in this case the CO2 concentration is not a stochastic
variable itself, but is governed by the balance between deterministic
weathering (sink) and stochastic volcanic outgassing (source).

$\dot{C}(t)=V(t)-W(C)\,,$ (8)

where $C$ is the total ocean atmosphere-ocean reservoir of CO2 , $V$ is the
volcanic CO2 outgassing rate, and $W$ is the rate of CO2 sequestration via
silicate weathering. Weathering is prescribed by the traditional WHAK
formulation [Walker . (1981)], where the weathering rate is exponentially
dependent on temperature changes,

$W(T)=k\exp\left(\frac{T-T_{e}}{T_{s}}\right)\,.$ (9)

Here, $k$ is a calibration constant and $T_{s}$ is the exponential scaling
factor. Note that we exclude direct dependence on carbon dioxide concentration
because, with it, weathering can no longer be calibrated to the same
equilibrium temperature at different solar forcing [Pierrehumbert (2010)].

To evaluate the weathering equation, we convert the total carbon $C$ to
temperature $T$ by first converting to concentration $f$ in ppm. To do this,
we use the simple ocean-atmosphere partitioning model of mills_timing_2011,

$p=0.78\frac{C}{C+h}\,,$ (10)

where $p$ is the CO2 partial pressure and $h$ is a constant 2.33$\times
10^{8}$ Tmole CO2 . From there, we assume a background partial pressure of
approximately 1 bar to compute the CO2 concentration $f$, which is used to
evaluate Equation 5. Finally, the resulting temperature is used to compute the
weathering and sequestration rate.

The outgassing rate $V$ fluctuates stochastically and is defined by an
Ornstein-Uhlenbech process,

$dV=\frac{1}{\tau}\left(\mu-V\right)dt+\sigma dB\,,$ (11)

where $\tau$ is a prescribed relaxation time scale, $\mu$ is a prescribed mean
outgassing rate, $\sigma$ scales the noise term, and $dB$ is a Weiner process
(Gaussian or Brownian noise). Note that this equation and its parameters are
distinct from Equation 6. Here, $\tau$ no longer represents a weathering time
scale because this stochastic process represents outgassing, not CO2
concentration directly.

All together, our model consists of two differential equations. The first is a
straightforward definition of how total ocean-atmosphere carbon changes in
response to the balance of weathering and outgassing. The second defines
stochastic variation in the CO2 outgassing rate.

$\displaystyle\dot{C}$ $\displaystyle=V(t)-W(C)$ (12) $\displaystyle\dot{V}$
$\displaystyle=\frac{1}{\tau}\left(\mu-V\right)+\sigma B(t)$ (13)

The difference between this model and the one discussed in Section 2.1 is
that, here, neither the carbon reservoir $C$ nor the CO2 concentration $f$ are
_themselves_ stochastic processes. The carbon sink and source terms are
separate. Weathering is deterministic and depends implicitly on temperature,
as defined in Equation 9. The outgassing _rate_ $V$ is stochastic and modeled
with an Ornstein-Uhlenbeck process. In contrast to the model in Section 2.1,
these equations explicitly represent imbalances in the carbon system that
drive long-term temperature change.

##  3 Simulations

We integrate the system of Equations 12 and 13 using the simple Euler-Maruyama
method over a range of different parameters and analyze the statistics of the
ensemble. We vary $\tau$ and $\sigma$ over discrete values between 105–109 yr
and 10-5–10-2 Tmole/yr, respectively, performing many simulations with each
available combination. We execute simulations in two stages. First, we perform
a relatively small number of simulations over each combination of $\tau$ and
$\sigma$. Second, we identify which combinations yield simulations with
appreciable, but not unrealistically high, temperature variability and perform
a much larger number of simulations using those combinations. This is a simple
strategy to avoid wasting computation on inconsequential $(\tau,\sigma)$
pairs. All told, the ensemble is comprised of 27,551,040 unique simulations.

Each simulation begins at $t=2$ Ga, when the equilibrium CO2 concentration is
slightly above 105 ppm, and ends at 0 Ga (present day). Each simulation is
also taken through a spin-up period of 0.5 Gyr to avoid starting them at the
exact equilibrium temperature, then integrated using 1 million time steps. To
prevent evaluating the logarithm of a non-positive number, $V$ is restricted
to be above zero and $C$ is restricted to be above a very small fraction of
its equilibrium value at modern insolation. These restrictions have virtually
no impact on our conclusions and only serve to prevent computation errors with
parameter combinations that produce unrealistically high variability. For each
simulation, we store the values of $C$, $V$, $T$, and $W$ at 17 equally spaced
points in time. We also record the maximum and minimum value of each quantity
along with the precise time of the temperature extrema. Table 1 consolidates
all of the model parameters along with their values and units.

Table 1: Static model parameters with values and units Parameter | Description | Value/Range | Units
---|---|---|---
$f_{0}$ | preindustrial CO2 concentration | 285 | ppm
$T_{e}$ | preindustrial temperature | 288 | K
$\alpha$ | planetary albedo | 0.3 | dimensionless
OLR0 | equilibrium outgoing longwave radiation | 239.05 | W/m2
$h$ | ocean-atmosphere partitioning constant | 2.33$\times$108 | Tmole
$\mu$ | mean volcanic outgassing rate | 7 | Tmole/yr
$t_{1}$ | simulation start time | 2.5 | Ga
$t_{2}$ | simulation end time | 0 | Ga
$t_{0}$ | temporal scaling of insolation function | 4.5 | Ga
$t_{\textrm{spinup}}$ | model spinup duration | 0.5 | Ga
$T_{s}$ | weathering temperature scaling | 11.1 | K
$T_{\textrm{snow}}$ | snowball threshold temperature | 280 | K
$k$ | weathering calibration constant | 7 | Tmole/yr
$a$ | OLR temperature scaling | 2 | W/m2/K
$b$ | OLR $\log f$ scaling | 5 | W/m2
$F_{0}$ | modern solar insolation | 1366 | W/m2
$\tau$ | outgassing relaxation timescale | 105 \- 109 | yr
$\sigma$ | outgassing deviation/variability | 10-5 \- 10-2 | Tmole/yr

##  4 Results and Discussion

For orientation, Figure 4 shows the results of a single simulation with
$\tau=30\times 10^{6}$ yr and $\sigma=2\times 10^{-4}$ Tmole/yr. The volcanic
outgassing rate, which is a realization of the Ornstein-Uhlenbeck process
defined in Equation 13, is shown on the top in red. The total atmosphere-ocean
carbon stock is shown on a log scale just below in blue. These are the
prognostic variables of the model. The CO2 concentration, temperature, and
weathering rate are each derived from $C$ as described in Section 2.

Figure 4: An example simulation of the system defined by Equations 12 and 13.
In panel (a), the stochastic outgassing rate $V$. In panel (b), the total
atmosphere-ocean CO2 reservoir $C$. These are the prognostic model variables.
Each of the remaining quantities is panels (c-e) is derived from $C$. The gray
line in each panel shows the equilibrium value or the mean value appropriate
for each quantity. For the outgassing and weathering in panels (a) and (e) it
shows the mean value of 7 Tmole/yr. In panels (b) and (c) it shows the value
required for the equilibrium temperature of $T_{0}=288$ K, over time. Finally,
in panel (d), it shows the equilibrium temperature of 288 K.

Figure 5 shows more realizations of the stochastic outgassing rate for
different combinations of $\tau$ and $\sigma$, stacked on top of each other.
When $\tau$ is small, the outgassing rate quickly relaxes to the mean and the
value of $\sigma$ must be relatively high to produce random fluctuations of
appreciable magnitude. In this regime, changes in the outgassing rate are
dominated by short-term, rapid excursions. In an intermediate regime, with
$\tau=30$ Myr, longer-term fluctuations begin to appear because the relaxation
time scale is longer, but rapid excursions are still evident. Finally, with
$\tau=300$ Myr, the outgassing rate is dominated by longer period drift with
occasional and smaller magnitude short-term shocks.

Figure 5: Realizations of the Ornstein-Uhlenbeck process defined by Equation
13 for three exemplary pairs of $\tau$ and $\sigma$. Five independent
realizations are stacked vertically on top of each other.

Figure 6 shows the principal results of our ensemble. It shows the first
temperature percentile, over time, for different values of $\tau$. The first
percentile is a somewhat arbitrary choice. However, the same results are
evident for whatever small percentile is chosen. Smaller values of $\tau$
produce more rapid fluctuations in the outgassing rate and the likelihood of
cold excursions reaching the snowball temperature increases markedly with
time. With higher values of $\tau$, the outgassing history transitions to
long-term variability, flattening the first temperature percentile curves.
When $\tau$ is small, snowballs are considerably more likely at later times.
When $\tau$ is large, the outgassing rate changes more slowly and snowball
likelihood is more uniform over time.

Figure 6: The value of the first temperature percentile over time in our
ensemble. Moving to the right, each panel shows results for a higher/longer
relaxation timescale $\tau$. Within each panel, results for different values
of $\sigma$ are delineated by color. The hypothetical snowball temperature
${T_{\textrm{snow}}=280}$ K is shown in each panel. When $\tau$ is small, the
outgassing rate is characterized by rapid deviations (see Figure 5) and the
likelihood of cold temperatures increases dramatically as time progresses. As
$\tau$ becomes larger, outgassing rates vary over longer periods and the
likelihood of cold temperatures becomes much more uniform over time. The first
percentile is an arbitrary choice, but other low quantiles demonstrate the
same basic relationships between $\tau$, $\sigma$, and the likelihood of low
temperatures.

Comparing the leftmost panel of Figure 6 to Figure 3 is revealing. Low values
of $\tau$ in Figure 6 recover similar results to wordsworth_how_2021, who
found that the likelihood of snowball events dramatically increases over time
along with the secular brightening of the Sun (as explained in Section 2.1).
For larger values of $\tau$, which permit longer-duration excursions in $V$
(Figure 5), the variation of the 1st percentile temperature over time is much
smaller. This is more consistent with the Earth’s geological record of
snowball episodes, which appear around 2.5 Ga, 720 Ma, and 560 Ma.

The results of this simple model suggest some constraints on the variability
of volcanic outgassing over Earth’s history. For temperatures to generally
remain above $\sim$280 K, yet not to preclude early snowball events, $\sigma$
may have been fairly small with a larger value of $\tau$. Qualitatively
speaking, variability in the carbon cycle and atmospheric CO2 must have
considerable long-term drift, not simply short-term excursions. This is
because a model with only short-term variability exhibits a very dramatic
preference for snowball events later in Earth’s history, contradicting the
geologic record.

Although it is difficult to reconstruct the history of volcanic outgassing,
geological evidence indicates values of $\tau>100$ Ma may be consistent with
the proxy record. McKenzie-Horton-Loomis-et-al-2016:continental used detrital
zircon ages as a proxy for continental arc-volcanism over the past 720 Myr and
found fluctuations in arc activity over the timescale of hundreds of Myr.
Furthermore, studies reconstructing paleogeographic volcanic arc distributions
show that variations in arc length of up to a factor of two over may be
responsible for icehouse–greenhouse transitions over the past billion years
[<]e.g.,¿Lee-Shen-Slotnick-et-al-2013:continental,Mills-Scotese-Walding-et-
al-2017:elevated. These long-term variations in outgassing rate require
$\tau>100$ Ma to be realized in our stochastic model.

##  5 Conclusion

In this paper, we discuss simple models of carbon variability deep into
Earth’s history. One model treats the atmospheric CO2 concentration directly
as a stochastic variable. While instructive, this model is limited in its
flexibility and can only produce results that are not compatible with the
Earth’s history of snowball events. Another model, introduced here, treats
carbon sources and sinks independently, and considers the effect of randomly
varying CO2 outgassing with deterministic weathering. This model is slightly
more complex, allowing for a range of different behavior that can be compared
with the first model and with the snowball record. When outgassing varies over
longer periods, specifically with a relaxation time scale of roughly
$\tau>100$ Myr, snowballs are possible early in Earth history, which is more
compatible with observations.

We note here have been several suggestions regarding the absence of snowball
events during the Phanerozoic, including the evolution of land plants, which
have been suggested to bring down atmospheric CO2 while simultaneously
increasing the strength of the silicate weathering thermostat
[<]e.g.,¿Berner-1994:geocarb, pierrehumbert_principles_2010, Dahl-
Arens-2020:impacts, Another suggestion involves the evolution of fungi and
lichens, which are also expected to increase weathering rates [Evans (2003),
Knauth Kennedy (2009)]. The increasing strength of biological regulatory
feedbacks may be have been responsible for the change in character of Earth’s
glacial cycles (from global-scale events to frequent glaciations that do not
extend to the low latitudes). Yet, at face value, the model of
wordsworth_how_2021 predicts that snowball episodes will become increasingly
likely in Earth’s future with increasing solar luminosity unless stabilizing
feedbacks, whether biological or otherwise, strengthen accordingly.

A recent study by arnscheidt_routes_2020 discussed different mechanisms for
transition into snowball, including fluctuations in radiative forcing and
changes in volcanic outgassing or weathering rate. They distinguish between
two routes to global glaciation: “route A” involving the climate system
entering a limit cycle regime, and “route B” involving a rapid fluctuation in
solar luminosity leading to a solitary “transient” glaciation. The negative
temperature excursions caused by variable volcanic outgassing in our model can
only explain “route A” glaciations, and are more consistent with the idea that
the Sturtian, Marinoan, and Gaskiers glaciations represent three iterations of
a limit cycle [<]as suggested by¿mills_timing_2011.

Variability in the carbon cycle has often been characterized as either
stochastic _or_ deterministic [Westerhold . (2020), Arnscheidt Rothman
(2021)]. The usual conceptual division is between rapid, unpredictable
volcanic forcing (stochastic) and a slow change in the background
“weatherability” of the continental land masses. Although this is sometimes a
useful construct, it is important to remember that both of these processes
can, and probably do, influence the climate simultaneously. For example, a
snowball may be triggered by the combination of slowly increasing
weatherability and a relatively sudden change to volcanic outgassing. Our
model results with long $\tau$ represent a simple demonstration of this
scenario, where outgassing undergoes both long-term drift and occasional
shorter-term shocks.

An important caveat of our new model is that we have assumed stochastically
varying volcanic outgassing to be the only driver of long-term climate
variability. In reality, there are many potential triggers for Snowball events
and long-term climate change that are unrelated to outgassing, including
changing continental configuration and composition [Donnadieu . (2004), Cox .
(2016)], evolution of the biosphere [Hedges (2004), Tziperman . (2011)],
sudden radiative forcing from volcanic aerosols [Macdonald Wordsworth (2017)],
or the collapse of a methane greenhouse [Schrag . (2002), Kopp . (2005)].
Indeed, if the Paleoproterozoic event was a transient glaciation caused by
aerosol forcing, the timing of the Cryogenian Snowballs late in Earth history
would be more consistent with the model of wordsworth_how_2021.

We note that there is a wide range of opportunities to explore randomness in
simple climate models. Our model has only considered one form of stochastic
outgassing with a deterministic weathering mechanism that is dependent
exclusively on temperature. This is, of course, an oversimplification
[Macdonald . (2019), Park . (2021), Baum . (2022)] and future work could
investigate model equations that incorporate variability in other terms in Eq.
5 such as the planetary albedo and weatherability (e.g., to represent the
emplacement of highly weatherable materials and changing continental
configuration). For example, it would be straightforward to include a randomly
varying weatherability factor that is also represented by an Ornstein-
Uhlenbeck process, but this term need not be stochastic at all.

##  6 Open Research

All of the code developed for this project was written in the Julia language
[Bezanson . (2012), Bezanson . (2017)] and is freely available on GitHub at
github.com/markmbaum/random-volcanic-climate and publicly archived on Zenodo
[Baum (2022)]. Additional code for the WHAK weathering equation can be found
in the GEOCLIM.jl module [Baum Fu (2022)]. Plots were created using the open
source Matplotlib library [Hunter (2007)].

###### Acknowledgements.

The computations in this paper were run on the FASRC Cannon cluster supported
by the FAS Division of Science Research Computing Group at Harvard University.
We thank Kaitlyn Loftus, Camille Hankel, Ned Kleiner, Stephen Bourguet, and
Katherine Keller for helpful and collaborative discussions. We thank Dorian
Abbot and an anonymous reviewer for helpful feedback that improved the
manuscript.

## References

  * Abbot (2016) abbot2016analyticalAbbot, DS.  2016\.  Analytical investigation of the decrease in the size of the habitable zone due to a limited CO2 outgassing rate Analytical investigation of the decrease in the size of the habitable zone due to a limited co2 outgassing rate.  The Astrophysical Journal8272117.
  * Arnscheidt Rothman (2020) arnscheidt_routes_2020Arnscheidt, CW. Rothman, DH.  2020\.  Routes to global glaciation Routes to global glaciation.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences4762239.  10.1098/rspa.2020.0303
  * Arnscheidt Rothman (2021) Arnscheidt-Rothman-2021:asymmetryArnscheidt, CW. Rothman, DH.  2021\.  Asymmetry of extreme Cenozoic climate–carbon cycle events Asymmetry of extreme cenozoic climate–carbon cycle events.  Science Advances733eabg6864.
  * Baum (2022) mark_baum_2022_6799309Baum, M.  202207\.  markmbaum/random-volcanic-climate: release-2. markmbaum/random-volcanic-climate: release-2.  Zenodo.  https://doi.org/10.5281/zenodo.6799309 10.5281/zenodo.6799309
  * Baum Fu (2022) mark_baum_2022_6533814Baum, M. Fu, M.  202205\.  markmbaum/GEOCLIM.jl: v0.1.12. markmbaum/geoclim.jl: v0.1.12.  Zenodo.  https://doi.org/10.5281/zenodo.6533814 10.5281/zenodo.6533814
  * Baum . (2022) baum2022sensitiveBaum, M., Fu, M. Bourguet, S.  2022\.  Sensitive Dependence of Global Climate to Continental Geometry Sensitive dependence of global climate to continental geometry.  Geophysical Research Letterse2022GL098843.
  * Berner (1990) Berner-1990:atmosphericBerner, RA.  1990\.  Atmospheric carbon dioxide levels over Phanerozoic time Atmospheric carbon dioxide levels over Phanerozoic time.  Science24949751382–1386.
  * Berner (1994) Berner-1994:geocarbBerner, RA.  1994\.  GEOCARB II: A revised model of atmospheric CO [sub 2] over phanerozoic time Geocarb ii: A revised model of atmospheric co [sub 2] over phanerozoic time.  American Journal of Science;(United States)2941.
  * Berner Lasaga (1989) Berner-Lasaga-1989:modelingBerner, RA. Lasaga, AC.  198903\.  Modeling the Geochemical Carbon Cycle Modeling the geochemical carbon cycle.  Scientific American22274-82.  10.1038/scientificamerican0389-74
  * Berner . (1983) Berner-Lasaga-Garrels-1983:carbonateBerner, RA., Lasaga, AC. Garrels, RM.  1983\.  The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years.  American Journal of Science2837641–683.  10.2475/ajs.283.7.641
  * Bezanson . (2017) bezanson2017juliaBezanson, J., Edelman, A., Karpinski, S. Shah, VB.  2017\.  Julia: A fresh approach to numerical computing Julia: A fresh approach to numerical computing.  SIAM review59165–98.  https://doi.org/10.1137/141000671
  * Bezanson . (2012) Bezanson2012Bezanson, J., Karpinski, S., Shah, VB. Edelman, A.  2012\.  Julia: A Fast Dynamic Language for Technical Computing Julia: A fast dynamic language for technical computing.  CoRRabs/1209.5145.  http://arxiv.org/abs/1209.5145
  * Cox . (2016) cox_continental_2016Cox, GM., Halverson, GP., Stevenson, RK., Vokaty, M., Poirier, A., Kunzmann, M.Macdonald, FA.  2016\.  Continental flood basalt weathering as a trigger for Neoproterozoic Snowball Earth Continental flood basalt weathering as a trigger for Neoproterozoic Snowball Earth.  Earth and Planetary Science Letters446.  10.1016/j.epsl.2016.04.016
  * Dahl Arens (2020) Dahl-Arens-2020:impactsDahl, TW. Arens, SK.  2020\.  The impacts of land plant evolution on Earth’s climate and oxygenation state–an interdisciplinary review The impacts of land plant evolution on earth’s climate and oxygenation state–an interdisciplinary review.  Chemical Geology547119665.
  * Dobrow (2016) dobrow_introduction_2016Dobrow, RP.  2016\.  Introduction to stochastic processes with R Introduction to stochastic processes with R.  Wiley.
  * Donnadieu . (2004) donnadieu_snowball_2004Donnadieu, Y., Goddéris, Y., Ramstein, G., Nédélec, A. Meert, J.  2004\.  A ‘snowball Earth’ climate triggered by continental break-up through changes in runoff A ‘snowball Earth’ climate triggered by continental break-up through changes in runoff.  Nature4286980.  10.1038/nature02408
  * Evans (2003) Evans-2003:fundamentalEvans, DA.  2003\.  A fundamental Precambrian–Phanerozoic shift in earth’s glacial style? A fundamental precambrian–phanerozoic shift in earth’s glacial style?  Tectonophysics3751-4353–385.
  * Evans . (1997) evans_low-latitude_1997Evans, DA., Beukes, NJ. Kirschvink, JL.  1997\.  Low-latitude glaciation in the Palaeoproterozoic era Low-latitude glaciation in the Palaeoproterozoic era.  Nature3866622.  10.1038/386262a0
  * Farrell Abbot (2012) Farrell-Abbot-2012:mechanismFarrell, B. Abbot, D.  2012\.  A mechanism for dust-induced destabilization of glacial climates A mechanism for dust-induced destabilization of glacial climates.  Climate of the Past862061–2067.
  * Feulner (2012) feulner2012faintFeulner, G.  2012\.  The faint young Sun problem The faint young sun problem.  Reviews of Geophysics502.
  * Franks . (2014) franks2014newFranks, PJ., Royer, DL., Beerling, DJ., Van de Water, PK., Cantrill, DJ., Barbour, MM. Berry, JA.  2014\.  New constraints on atmospheric CO2 concentration for the Phanerozoic New constraints on atmospheric co2 concentration for the phanerozoic.  Geophysical Research Letters41134685–4694.
  * Franzke . (2015) franzke2015stochasticFranzke, CL., O’Kane, TJ., Berner, J., Williams, PD. Lucarini, V.  2015\.  Stochastic climate theory and modeling Stochastic climate theory and modeling.  Wiley Interdisciplinary Reviews: Climate Change6163–78.
  * Gough (1981) gough_solar_1981Gough, DO.  1981\.  Solar interior structure and luminosity variations Solar interior structure and luminosity variations.  Solar Physics74.
  * Graham (2021) Graham-2021:highGraham, RJ.  2021\.  High pCO2 reduces sensitivity to CO2 perturbations on temperate, Earth-like planets throughout most of habitable zone High pco2 reduces sensitivity to co2 perturbations on temperate, earth-like planets throughout most of habitable zone.  Astrobiology21111406–1420.
  * Greenwood Wing (1995) Greenwood-Wing-1995:eoceneGreenwood, DR. Wing, SL.  1995\.  Eocene continental climates and latitudinal temperature gradients Eocene continental climates and latitudinal temperature gradients.  Geology23111044–1048.
  * Halevy Bachan (2017) Halevy-Bachan-2017:geologicHalevy, I. Bachan, A.  2017\.  The geologic history of seawater pH The geologic history of seawater ph.  Science35563291069–1071.
  * Hasselmann (1976) hasselmann_stochastic_1976Hasselmann, K.  1976\.  Stochastic climate models Part I. Theory Stochastic climate models Part I. Theory.  Tellus286.  10.1111/j.2153-3490.1976.tb00696.x
  * Hedges (2004) Hedges-2004:molecularHedges, SB.  2004\.  Molecular clocks and a biological trigger for Neoproterozoic Snowball Earth events and the Cambrian explosion Molecular clocks and a biological trigger for neoproterozoic snowball earth events and the cambrian explosion.  SYSTEMATICS ASSOCIATION SPECIAL VOLUME6627–40.
  * Hoffman . (1998) Hoffman-Kaufman-Halverson-et-al-1998:neoproterozoicHoffman, PF., Kaufman, AJ., Halverson, GP. Schrag, DP.  1998aug 28.  A Neoproterozoic snowball Earth A Neoproterozoic snowball Earth.  Science28153811342-1346.
  * Hoffman Schrag (2002) hoffman_snowball_2002Hoffman, PF. Schrag, DP.  2002\.  The snowball Earth hypothesis: testing the limits of global change The snowball Earth hypothesis: testing the limits of global change.  Terra Nova143.  10.1046/j.1365-3121.2002.00408.x
  * Hunter (2007) Hunter:2007Hunter, JD.  2007\.  Matplotlib: A 2D graphics environment Matplotlib: A 2d graphics environment.  Computing in Science & Engineering9390–95.  10.1109/MCSE.2007.55
  * Imkeller Monahan (2002) imkeller2002conceptualImkeller, P. Monahan, AH.  2002\.  Conceptual stochastic climate models Conceptual stochastic climate models.  Stochastics and Dynamics203311–326.
  * Imkeller Von Storch (2001) imkeller2001stochasticImkeller, P. Von Storch, JS.  2001\.  Stochastic climate models Stochastic climate models ( 49).  Springer Science & Business Media.
  * Jacobs (2010) jacobs_stochastic_2010Jacobs, K.  2010\.  Stochastic Processes for Physicists: Understanding Noisy Systems Stochastic Processes for Physicists: Understanding Noisy Systems.  Cambridge University Press.  10.1017/CBO9780511815980
  * Kirschvink . (2000) kirschvink_paleoproterozoic_2000Kirschvink, JL., Gaidos, EJ., Bertani, LE., Beukes, NJ., Gutzmer, J., Maepa, LN. Steinberger, RE.  2000\.  Paleoproterozoic snowball Earth: Extreme climatic and geochemical global change and its biological consequences Paleoproterozoic snowball Earth: Extreme climatic and geochemical global change and its biological consequences.  Proceedings of the National Academy of Sciences974.  10.1073/pnas.97.4.1400
  * Knauth Kennedy (2009) Knauth-Kennedy-2009:lateKnauth, LP. Kennedy, MJ.  2009\.  The late Precambrian greening of the Earth The late precambrian greening of the earth.  Nature4607256728–732.
  * Koll Cronin (2018) Koll-Cronin-2018:earthKoll, DD. Cronin, TW.  2018\.  Earth’s outgoing longwave radiation linear due to H2O greenhouse effect Earth’s outgoing longwave radiation linear due to h2o greenhouse effect.  Proceedings of the National Academy of Sciences1154110293–10298.
  * Kopp . (2005) Kopp-Kirschvink-Hilburn-et-al-2005:paleoproterozoicKopp, RE., Kirschvink, JL., Hilburn, IA. Nash, CZ.  2005\.  The Paleoproterozoic snowball Earth: a climate disaster triggered by the evolution of oxygenic photosynthesis The paleoproterozoic snowball earth: a climate disaster triggered by the evolution of oxygenic photosynthesis.  Proceedings of the National Academy of Sciences1023211131–11136.
  * Krissansen-Totton . (2018) Krissansen-Totton-Arney-Catling-2018:constrainingKrissansen-Totton, J., Arney, GN. Catling, DC.  2018\.  Constraining the climate and ocean pH of the early Earth with a geological carbon cycle model Constraining the climate and ocean ph of the early earth with a geological carbon cycle model.  Proceedings of the National Academy of Sciences115164105–4110.
  * Lee . (2013) Lee-Shen-Slotnick-et-al-2013:continentalLee, CTA., Shen, B., Slotnick, BS., Liao, K., Dickens, GR., Yokoyama, Y.others  2013\.  Continental arc–island arc fluctuations, growth of crustal carbonates, and long-term climate change Continental arc–island arc fluctuations, growth of crustal carbonates, and long-term climate change.  Geosphere9121–36.
  * Lenardic . (2016) lenardic2016climateLenardic, A., Jellinek, A., Foley, B., O’Neill, C. Moore, W.  2016\.  Climate-tectonic coupling: Variations in the mean, variations about the mean, and variations in mode Climate-tectonic coupling: Variations in the mean, variations about the mean, and variations in mode.  Journal of Geophysical Research: Planets121101831–1864.
  * Macdonald . (2019) macdonald_arc-continent_2019Macdonald, FA., Swanson-Hysell, NL., Park, Y., Lisiecki, L. Jagoutz, O.  2019\.  Arc-continent collisions in the tropics set Earth’s climate state Arc-continent collisions in the tropics set Earth’s climate state.  Science.  10.1126/science.aav5300
  * Macdonald Wordsworth (2017) Macdonald-Wordsworth-2017:initiationMacdonald, FA. Wordsworth, R.  2017\.  Initiation of Snowball Earth with volcanic sulfur aerosol emissions Initiation of snowball earth with volcanic sulfur aerosol emissions.  Geophysical Research Letters4441938–1946.
  * Marshall . (1988) Marshall-Walker-Kuhn-1988:longMarshall, HG., Walker, JCG. Kuhn, WR.  1988\.  Long-term climate change and the geochemical cycle of carbon Long-term climate change and the geochemical cycle of carbon.  Journal of Geophysical Research: Atmospheres93D1791–801.
  * McKenzie . (2016) McKenzie-Horton-Loomis-et-al-2016:continentalMcKenzie, NR., Horton, BK., Loomis, SE., Stockli, DF., Planavsky, NJ. Lee, CTA.  2016\.  Continental arc volcanism as the principal driver of icehouse-greenhouse variability Continental arc volcanism as the principal driver of icehouse-greenhouse variability.  Science3526284444–447.
  * Mills . (2017) Mills-Scotese-Walding-et-al-2017:elevatedMills, BJW., Scotese, CR., Walding, NG., Shields, GA. Lenton, TM.  2017\.  Elevated CO2 degassing rates prevented the return of Snowball Earth during the Phanerozoic Elevated co2 degassing rates prevented the return of snowball earth during the phanerozoic.  Nature communications811–7.
  * Mills . (2011) mills_timing_2011Mills, BJW., Watson, AJ., Goldblatt, C., Boyle, R. Lenton, TM.  2011\.  Timing of Neoproterozoic glaciations linked to transport-limited global weathering Timing of Neoproterozoic glaciations linked to transport-limited global weathering.  Nature Geoscience412.  10.1038/ngeo1305
  * Montañez . (2016) montanez2016climateMontañez, IP., McElwain, JC., Poulsen, CJ., White, JD., DiMichele, WA., Wilson, JP.Hren, MT.  2016\.  Climate, pCO2 and terrestrial carbon cycle linkages during late Palaeozoic glacial–interglacial cycles Climate, pco2 and terrestrial carbon cycle linkages during late palaeozoic glacial–interglacial cycles.  Nature Geoscience911824–828.
  * Myhre . (1998) myhre1998newMyhre, G., Highwood, EJ., Shine, KP. Stordal, F.  1998\.  New estimates of radiative forcing due to well mixed greenhouse gases New estimates of radiative forcing due to well mixed greenhouse gases.  Geophysical research letters25142715–2718.
  * Park . (2021) park2021evaluatingPark, Y., Swanson-Hysell, NL., Lisiecki, LE. Macdonald, FA.  2021\.  Evaluating the relationship between the area and latitude of large igneous provinces and Earth’s long-term climate state Evaluating the relationship between the area and latitude of large igneous provinces and earth’s long-term climate state.  Large igneous provinces: A driver of global environmental and biotic changes153–168.
  * Pierrehumbert (2010) pierrehumbert_principles_2010Pierrehumbert, RT.  2010\.  Principles of Planetary Climate Principles of Planetary Climate (1 ).  Cambridge University Press.  10.1017/CBO9780511780783
  * Pierrehumbert . (2011) pierrehumbert_climate_2011Pierrehumbert, RT., Abbot, DS., Voigt, A. Koll, D.  2011\.  Climate of the Neoproterozoic Climate of the Neoproterozoic.  Annual Review of Earth and Planetary Sciences391.  10.1146/annurev-earth-040809-152447
  * Prave . (2016) Prave-Condon-Hoffmann-et-al-2016:durationPrave, AR., Condon, DJ., Hoffmann, KH., Tapster, S. Fallick, AE.  2016\.  Duration and nature of the end-Cryogenian (Marinoan) glaciation Duration and nature of the end-Cryogenian (Marinoan) glaciation.  Geology448631–634.
  * Rooney . (2015) Rooney-Strauss-Brandon-et-al-2015:cryogenianRooney, AD., Strauss, JV., Brandon, AD. Macdonald, FA.  2015\.  A Cryogenian chronology: Two long-lasting synchronous Neoproterozoic glaciations A Cryogenian chronology: Two long-lasting synchronous Neoproterozoic glaciations.  Geology435459–462.
  * Schrag . (2002) schrag_initiation_2002Schrag, DP., Berner, RA., Hoffman, PF. Halverson, GP.  2002\.  On the initiation of a snowball Earth: Initiation of a Snowball Earth On the initiation of a snowball Earth: Initiation of a Snowball Earth.  Geochemistry, Geophysics, Geosystems36.  10.1029/2001GC000219
  * Tziperman . (2011) tziperman_biologically_2011Tziperman, E., Halevy, I., Johnston, DT., Knoll, AH. Schrag, DP.  2011\.  Biologically induced initiation of Neoproterozoic snowball-Earth events Biologically induced initiation of Neoproterozoic snowball-Earth events.  Proceedings of the National Academy of Sciences10837.  10.1073/pnas.1016361108
  * Urey (1952) Urey-1952:earlyUrey, HC.  1952\.  On the early chemical history of the earth and the origin of life On the early chemical history of the earth and the origin of life.  Proceedings of the National Academy of Sciences384351.
  * Walker . (1981) Walker-Hays-Kasting-1981:negativeWalker, JCG., Hays, PB. Kasting, JF.  1981\.  A negative feedback mechanism for the long-term stabilization of Earth’s surface temperature A negative feedback mechanism for the long-term stabilization of Earth’s surface temperature.  Journal of Geophysical Research: Oceans86C109776–9782.
  * Westerhold . (2020) Westerhold-Marwan-Drury-et-al-2020:astronomicallyWesterhold, T., Marwan, N., Drury, AJ., Liebrand, D., Agnini, C., Anagnostou, E.others  2020\.  An astronomically dated record of Earth’s climate and its predictability over the last 66 million years An astronomically dated record of earth’s climate and its predictability over the last 66 million years.  Science36965091383–1387.
  * Wordsworth (2021) wordsworth_how_2021Wordsworth, R.  2021\.  How likely are Snowball episodes near the inner edge of the habitable zone? How likely are Snowball episodes near the inner edge of the habitable zone?  The Astrophysical Journal Letters9121.  10.3847/2041-8213/abf7c7
  * Zeebe . (2017) Zeebe-Westerhold-Littler-et-al-2017:orbitalZeebe, RE., Westerhold, T., Littler, K. Zachos, JC.  2017\.  Orbital forcing of the Paleocene and Eocene carbon cycle Orbital forcing of the paleocene and eocene carbon cycle.  Paleoceanography325440–465.