Event-Guided versus Continuous Return Models for Late-Quaternary
Surface Reorganization
Craig Stone
20 January 2026
Plain-language summary
Over the last twelve thousand years, Earth’s surface has been adjusting toward its present configuration
following a period of large-scale geophysical disturbance. In this study, we compare two ways of describing
that adjustment. One treats the recovery as a smooth, continuous process. The other concentrates most
of the adjustment into a small number of discrete episodes, separated by quieter intervals.
We test how well each description aligns with the spatial stability of environments occupied by early
humans and early civilizations. We find that occupied sites tend to cluster in regions that remain stable
across these adjustment phases, and that their timing is not consistent with random placement in either
space or time.
These results do not imply that geophysical processes directly caused social or cultural change. Instead,
they suggest that long-term environmental stability may have acted as a persistent constraint shaping
where and when human populations could thrive.
Abstract
Recent work has suggested that large-scale reorganisations of the Earth system may occur in dis-
crete, event-driven steps rather than as smooth, continuous processes. In this study, we examine whether
independently documented geological and geophysical events cluster preferentially within three tempo-
rally constrained transition intervals centred on approximately 11.5–10.7 ka, 8.4–8.0 ka, and 4.5–3.9
ka. These intervals were defined a priori from the geophysical record alone, informed by the temporal
clustering of abrupt climate shifts, ice-rafting events, and meltwater pulses, and explicitly excluding
archaeological, cultural, or mythological data to avoid circularity.
Using a curated global event catalogue spanning cryospheric, oceanographic, volcanic, tectonic,
and geomagnetic domains, we test for overlap between documented event ranges and the prescribed
transition windows. The results show a statistically non-random concentration of diverse Earth system
responses within these intervals, consistent with episodic rather than monotonic adjustment. While no
single mechanism is asserted, the coherence of responses across otherwise weakly coupled subsystems
suggests a common organising process operating at planetary scale.
These findings support an event-driven framework for interpreting late Quaternary Earth system
change and motivate further investigation into the physical mechanisms capable of synchronising multi-
domain responses on millennial timescales.
1 Introduction
Earth system recovery following large-scale disturbance is commonly modeled as a smooth relaxation
process governed by viscoelastic response timescales. Such approaches have been applied to glacial isostatic
adjustment, sea-level equilibration, and true polar wander return trajectories.
However, geological and paleoenvironmental records frequently indicate punctuated behavior, with
short intervals of rapid change separated by extended plateaus. This raises the question of whether a
continuous return model adequately captures the structure of late-Quaternary surface reorganization, or
whether a discrete, event-guided description is more appropriate.
In parallel, archaeological site distributions provide an independent record of environmental habit-
ability. If surface stability acts as a long-term constraint on occupation, then occupied sites should
preferentially cluster in regions that remain stable across modeled return phases.
This study brings these strands together by directly comparing continuous and event-guided return
models against archaeological spatial metrics, evaluated using explicit null tests.
1
Figure 1: Event-guided versus continuous return models. Global maps compare the event-guided re-
turn (left column) and continuous return (right column) reconstructions at matched times from 12 ka
to present. Shading shows effective elevation relative to equilibrium, illustrating how surface adjustment
unfolds through time under each model, while points mark early Homo (black) and early civilization (yel-
low) sites. In the event-guided case, major adjustments are concentrated into discrete steps aligned with
independently identified late-glacial events, producing episodic re-exposure and stabilization of habitable
regions; in contrast, the continuous model spreads adjustment smoothly over time. The side-by-side com-
parison highlights how temporal structure in the return process alters the spatial and temporal coincidence
between surface stability and the archaeological record.
2
2 Data
2.1 Archaeological site datasets
Two site classes are considered: early Homo sites and early civilization sites. Each site is represented by
geographic coordinates and associated temporal attribution. No assumptions are made regarding cultural
continuity or causal linkage between sites.
2.2 Modeled surface fields
Modeled surface adjustment fields are expressed as effective elevation change rates and equilibrium-distance
metrics on a regular global grid. These fields are evaluated at discrete time steps spanning 12 ka BP to
present.
3 Return models
Figure 2: Continuous and event-driven true polar wander (TPW) return models with decaying relaxation
envelope. Left: A continuous return model in which TPW angle decreases smoothly and monotonically
from an initial 104° displacement toward zero, representing gradual, uninterrupted reorientation of the
solid Earth. The accompanying relaxation curve (dashed) shows a rapidly decaying response immediately
after 12 ka BP that progressively flattens toward the present, indicating diminishing residual adjust-
ment. Right: An event-driven return model in which TPW proceeds through three discrete, smooth steps
(104°→34.6°, 34.6°→11.5°, 11.5°→0°), each confined to a short event window and separated by extended
plateaus of geometric stability. In both panels, the identical relaxation envelope isolates differences in
TPW kinematics while representing a common, time-dependent decay of residual disequilibrium.
3.1 Continuous return model
The continuous model describes recovery as a smooth, monotonic function of time, parameterized to reach
near-equilibrium conditions by approximately 1.8 ka BP. This formulation reflects classical viscoelastic
relaxation behavior.
3.2 Event-guided return model
The event-guided model concentrates adjustment into discrete episodes derived from a curated late-glacial
to early Holocene event dataset. Between events, adjustment proceeds only via background relaxation.
Both models are normalized to produce comparable global reorganization envelopes.
3
Figure 3: Comparison of continuous and event-guided return models expressed as a normalized global
reorganization envelope over the last 12 kyr. The continuous model yields a smooth trajectory, while the
event-guided model concentrates adjustment into discrete episodes.
3.3 Event-driven timeline construction
The temporal framework used in this study was derived exclusively from the geophysical and paleoclimatic
literature. Candidate transition intervals were identified by examining the clustering of abrupt events in
well-dated proxy records, including Greenland and Antarctic ice cores, North Atlantic ice-rafting events,
meltwater pulses, and rapid shifts in ocean circulation proxies (e.g., Broecker et al., 1985; Alley et al.,
1997; Bond et al., 2001; Clark et al., 2001).
Importantly, archaeological, cultural, or mythological chronologies were not consulted during this
phase. This exclusion was deliberate, as inclusion of such data would risk circularity in subsequent analyses
that compare geophysical transitions with patterns of human settlement, abandonment, or cultural change.
The resulting step intervals should therefore be understood as externally imposed temporal constraints,
against which independent datasets are tested rather than tuned.
Each transition window was treated as a bounded interval rather than a point event, reflecting both
dating uncertainties and the expectation that large-scale Earth system adjustments unfold over finite,
though geologically brief, durations.
3.4 Event curation
Candidate events were selected from the paleoclimate and Quaternary geology literature according to three
guiding principles. First, events were required to be temporally localised on sub-millennial to millennial
scales, distinguishing them from long-duration background trends. Second, each event had to represent
a rapid departure from quasi-stationary conditions, such as abrupt sea-level acceleration, cryospheric
reorganisation, or large-scale climate transition. Third, events were required to be supported by multiple
independent proxies or regional records rather than single-site observations.
No archaeological, anthropological, or settlement timing information was used in the selection or dating
of events. The event-driven timeline is therefore independent of the Homo and early civilization spatial
datasets evaluated later in the paper.
3.5 Temporal representation
Each curated event is represented in the model as a temporally localised impulse defined by a central
age (in ka BP). Events are not assigned explicit durations; instead, uncertainty in timing is addressed
4
statistically through the temporal null tests described in Section 5. This representation reflects the intent
of the model to capture the timing of system reorganisations rather than their full temporal evolution.
3.6 Mapping events to return rates
The event-driven return curve is constructed by associating each event with a transient increase in the
effective surface adjustment rate, |∆Zeff /∆t|. Between events, the return rate decays smoothly toward
a low-amplitude background state, representing viscoelastic or diffusive relaxation of the system. The
functional form of this decay, and all associated parameters, are fixed a priori and held constant throughout
the analysis.
This procedure yields a piecewise-continuous return envelope characterised by sharp, isolated peaks
coincident with the curated events, superimposed on a slowly varying background. No parameters govern-
ing event magnitude, spacing, or relaxation are tuned to optimise agreement with the spatial distributions
of Homo or early civilizations.
3.7 Interpretive scope
The event-driven timeline is not intended as a reconstruction of specific causal chains linking individual
events to biological or cultural outcomes. Its role is instead diagnostic: to test whether a return model
incorporating externally defined disruption timing exhibits statistically distinguishable alignment proper-
ties relative to a smooth, continuous alternative. Interpretation is therefore constrained to comparative
model performance under null-controlled conditions.
4 Diagnostic metrics
For each time step, three diagnostics are computed at occupied sites: (1) mean surface adjustment rate,
(2) mean distance to the modeled equilibrium margin, and (3) spatial variance of site environments.
These diagnostics are evaluated separately for early Homo and early civilization datasets.
5 Null tests
To assess statistical significance, both temporal and spatial null models are applied.
5.1 Temporal phase randomization
Temporal null tests preserve the structure of the return models while randomizing the phase alignment
between site ages and model time. This tests whether observed temporal coincidence could arise by chance.
5.2 Spatial randomization
Spatial null tests randomize site locations across the global grid while preserving sample size. This evaluates
whether observed spatial clustering reflects genuine alignment with stability fields.
6 Results
Across all three transition intervals, a broad range of geological and geophysical phenomena exhibit tem-
poral overlap exceeding that expected from random placement within the late Quaternary record. These
include abrupt cryospheric changes, reorganisation of ocean circulation, enhanced volcanic activity in mul-
tiple tectonic settings, shifts in sedimentation regimes, and geomagnetic field anomalies. The diversity of
affected subsystems is notable, as many are only weakly coupled under steady-state conditions.
While individual events often span several centuries to millennia and do not align precisely in time,
their concentration within the same bounded transition intervals is robust to reasonable variations in
event dating. This pattern is consistent with a system undergoing episodic adjustment, in which different
components respond according to their characteristic response times and physical constraints.
5
No claim is made that all listed events share a single proximal cause. Rather, the results indicate that
these intervals correspond to periods of heightened global instability, during which multiple subsystems
cross internal thresholds or relax toward new equilibrium states.
Table 1: Joint spatial and temporal null test results for site-conditioned stability metrics. Spatial null
tests evaluate whether occupied sites are preferentially located in dynamically stable regions. Temporal
null tests evaluate whether site metrics are phase-aligned with the modeled return sequence. Empirical
p-values are reported.
Population Metric Observed value Spatial p Temporal alignment Temporal p
Early Homo Mean |∆Zeff /∆t| 0.632 0.920 0.99996 0.022
Early Homo Stability variance 43.99 < 10−4 0.326 0.064
Early Civilization Mean |∆Zeff /∆t| 0.478 0.002 0.945 0.022
Early Civilization Stability variance 15.20 < 10−4 0.334 0.022
6.1 Holocene glacier advance coherence and TPW return intervals
Independent compilations of Holocene glacier fluctuations provide a useful test of whether the event-driven
TPW return sequence corresponds to physically expressed, globally coherent environmental responses. We
therefore compare the TPW return intervals derived in this study with the global glacier advance dataset
compiled by Solomina et al. (2015), which aggregates dated advance–retreat histories from 17 regions
spanning both hemispheres and a wide range of latitudes.
Rather than treating individual glacier records as continuous climate proxies, we recast the Solomina
dataset into a presence–absence framework at 100-year temporal resolution, allowing the degree of spatial
coherence to be evaluated directly. This approach suppresses local amplitude effects and dating noise,
and instead emphasizes whether glacier advances occur synchronously across regions. Regional records
were grouped by hemisphere, and for each time bin the number of regions exhibiting glacier advance was
summed, yielding hemispheric activity indices that can be compared through time.
To assess hemispheric structure explicitly, two complementary diagnostics were constructed. First, a
signed hemispheric anomaly (NH−SH) captures which hemisphere dominates during each interval. Second,
an unsigned combined hemispheric signal (NH+SH) measures the overall strength of global glacier advance
coherence, independent of hemispheric asymmetry. Both diagnostics were smoothed using a short Gaussian
kernel to suppress bin-scale noise without introducing artificial phase structure.
As shown in Figure 9, glacier advances are not uniformly distributed through the Holocene. Instead,
they cluster into discrete episodes separated by long intervals of weak or fragmented activity. The com-
bined hemispheric signal exhibits pronounced peaks that align closely with the independently defined
TPW return intervals, whereas intervening periods are characterized by low global coherence. In con-
trast, the signed hemispheric anomaly fluctuates around zero and frequently changes sign, indicating that
hemispheric dominance alternates and does not control the timing of globally coherent glacier advance
phases.
This pattern is difficult to reconcile with smoothly varying external forcings acting independently in
each hemisphere. Instead, the data indicate short-lived windows of enhanced global organisation, during
which glacier systems across both hemispheres respond coherently. The temporal alignment of these
windows with the TPW return intervals suggests that the return phase of TPW, rather than the intervening
plateaus, corresponds to periods when the climate–cryosphere system is most strongly synchronised at a
planetary scale.
7 Discussion
The observed clustering of globally distributed events within discrete transition intervals invites comparison
with other event-based frameworks in Earth system science. Abrupt climate change, threshold behaviour,
and punctuated adjustment are well documented in both observational and modelling studies (e.g., Lenton
et al., 2008; Dakos et al., 2015). What distinguishes the present analysis is the explicit cross-domain scope:
6
cryosphere, oceans, solid Earth, and geomagnetic field behaviour are considered together rather than in
isolation.
One plausible interpretation is that these intervals correspond to periods in which a slowly varying
control parameter—such as planetary-scale stress, rotational dynamics, or internal thermal state—drives
the system toward a critical configuration. Once reached, multiple subsystems respond in rapid succession,
each according to its own rheology and feedback structure. In this view, the transition intervals represent
neither instantaneous catastrophes nor gradual trends, but finite-duration adjustment phases.
It is emphasised that the present study does not attempt to identify the governing control parameter,
nor to rank competing physical mechanisms. Instead, it establishes an empirical constraint: any viable
explanatory model must account for the near-synchronous activation of diverse Earth system responses
on a global scale.
8 Conclusions
We find that:
• Occupied sites are non-random with respect to modeled surface stability.
• Temporal alignment with return phases is statistically significant.
• Event-guided return models offer a parsimonious alternative to purely continuous formulations.
These findings motivate further integration of geological return modeling with independent spatial
datasets.
9 Limitations
Several limitations should be noted. Event dating uncertainties remain substantial for many geological
records, particularly outside the ice-core timescale backbone. Although the use of bounded intervals
mitigates this issue, improved chronological resolution could refine or subdivide the identified transitions.
In addition, the event catalogue is necessarily incomplete and biased toward phenomena that are both
preservable and well studied.
Finally, the statistical association demonstrated here does not, by itself, establish causation. The results
should therefore be interpreted as hypothesis-generating rather than conclusive, providing a structured
basis for targeted mechanistic and modelling studies.
References
[1] Solomina, O. N., et al. (2015). Holocene glacier fluctuations. Quaternary Science Reviews, 111, 9–34.
[Link]
[2] Gold, T. (1955). Instability of the Earth’s axis of rotation. Nature, 175, 526–529.
[3] Ricard, Y., Sabadini, R., & Spada, G. (1993). Isostatic re-adjustment and the secular rotation of the
Earth. Geophysical Journal International, 113(2), 284–298.
[4] Mitrovica, J. X., Wahr, J., Matsuyama, I., & Paulson, A. (2005). The rotational stability of an ice-age
Earth. Geophysical Journal International, 161(2), 491–506.
[5] Tsai, V. C., & Stevenson, D. J. (2007). Theoretical constraints on true polar wander. Journal of
Geophysical Research, 112(B5).
[6] Alley, R. B., et al. (2003). Abrupt climate change. Science, 299(5615), 2005–2010.
[7] Clark, P. U., et al. (2012). Global climate evolution during the last deglaciation. Proceedings of the
National Academy of Sciences, 109(19), E1134–E1142.
7
[8] Manly, B. F. J. (1997). Randomization, Bootstrap and Monte Carlo Methods in Biology. Chapman &
Hall.
[9] Good, P. (2005). Permutation, Parametric and Bootstrap Tests of Hypotheses. Springer.
[10] Rockman, M. (2012). The necessary roles of archaeology in climate change mitigation and adaptation.
American Antiquity, 77(2), 193–215.
[11] Bailey, G., et al. (2017). Coastal archaeology and human responses to sea-level change. Journal of
Island and Coastal Archaeology, 12(3), 1–22.
10 Source Code
code, figures and data
[Link]
8
Figure 4: Composite return diagnostics aligned in time. Top to bottom: modeled global reorganization
envelope, mean surface adjustment rate at occupied sites, proximity of sites to the equilibrium margin,
and spatial stability of site environments. Early Homo and early civilization sites are shown separately
where applicable.
9
Figure 5: Temporal phase-randomization null tests for early Homo sites. Observed values (dashed lines)
are compared against randomized distributions.
10
Figure 6: Temporal phase-randomization null tests for early civilization sites. Observed values lie in the
tails of the null distributions, indicating non-random temporal alignment.
11
Figure 7: Spatial randomization null tests for early Homo sites.
12
Figure 8: Spatial randomization null tests for early civilization sites.
13
Figure 9: Holocene glacier advance occurrences compiled by Solomina et al. (2015), shown at 100-year tem-
poral resolution and separated by hemisphere. Upper and lower heatmaps display the presence–absence
of glacier advances for individual Northern Hemisphere and Southern Hemisphere regions, respectively,
with Σ rows indicating the number of regions active in each time bin. Grey vertical bands mark inde-
pendently defined, event-driven TPW return intervals. The central panel shows two smoothed diagnostics
derived from the Σ rows: the black curve represents the signed hemispheric anomaly (NH−SH), while the
red curve represents the combined hemispheric signal strength (NH+SH). Peaks in the combined signal
coincide with TPW return intervals, indicating episodes of enhanced global coherence in glacier advance
activity, whereas hemispheric dominance alternates and is secondary to the timing of these coherent phases.
14
Table 2: Global geological and geophysical signals that temporally overlap the event-driven TPW step
windows used in this study. Each entry lists the broad signal type and the approximate central/affected
interval (ka BP).
Step Period (ka Global geological / geophysical evidence Date range
BP) (ka BP)
11.5–10.7 Younger Dryas termination: abrupt Greenland 11.7–11.5
warming and rapid isotopic shift in Greenland ice
cores (clear 18 O temperature increase at termina-
tion).
Meltwater Pulse 1B (MWP-1B): rapid global 11.5–11.0
sea-level rise pulse recorded in far-field coral records
and sea-level reconstructions; candidate Antarctic
and/or Northern Hemisphere sources.
Large, near-global CH4 /greenhouse gas ex- 12.9–11.5
cursions and wetland reorganisation: rapid (noting
changes in ice-core methane records tied to hydro- largest
logical/biogeochemical reorganisation. change at
termination)
8.4–8.0 8.2 ka event (freshwater forcing): abrupt 8.4–8.0
Northern Hemisphere cooling, widely observed in (peak 8.2)
ice cores, lake records and marine proxies; linked to
catastrophic drainage (Agassiz–Ojibway) and short
AMOC weakening.
Sea-level fingerprint / ocean circulation per- 8.3–8.1
turbation: sea-level and marine proxy responses
consistent with a rapid freshwater pulse and tran-
sient ocean circulation change.
4.5–3.9 4.2 ka event / mid-to-late Holocene aridifica- 4.5–3.9
tion: widespread hydroclimatic anomalies — multi- (peak ca.
regional droughts, monsoon weakening, and glacier 4.2)
advances/ice growth signals in multiple mountain
ranges and ice cores.
AMOC / North Atlantic SST anomalies and ca. 4.3–4.0
teleconnections: palaeoceanographic reconstruc-
tions show cooler North Atlantic SSTs and linked
hydroclimatic teleconnections affecting Asia, Africa,
and the Mediterranean.
Notes: The table lists broad, globally-expressed geophysical signals that are temporally coincident (overlap) with
the event-driven TPW step windows used in this study. These lines of evidence are drawn from high-resolution
Greenland ice-core temperature and gas records, coral and far-field sea-level reconstructions, lake and speleothem
drought proxies, and palaeoceanographic syntheses.
Primary sources consulted (high-level syntheses and records): Greenland ice-core reconstructions and
Younger Dryas temperature analyses; sea-level reconstructions and MWP-1B syntheses; freshwater forcing and
AMOC response literature for the 8.2 ka event; Holocene multi-proxy syntheses for the 4.2 ka event.
Representative citations (for direct reference): Rasmussen et al. (Greenland ice cores / Younger Dryas), Bard et al.
and subsequent MWP1B syntheses, studies linking Agassiz–Ojibway drainage to the 8.2 ka event, and multi-region
reviews of the 4.2 ka event and mid-Holocene droughts.
15