A rigorous detection of Milankovitch periodicities in volcanic output across the Pleistocene-Holocene ice age has remained elusive. We report on a spectral analysis of a large number of well-preserved ash plume deposits recorded in marine sediments along the Pacific Ring of Fire. Our analysis yields a statistically significant detection of a spectral peak at the obliquity period. We propose that this variability in volcanic activity results from crustal stress changes associated with ice age mass redistribution. In particular, increased volcanism lags behind the highest rate of increasing eustatic sea level (decreasing global ice volume) by 4.0 ± 3.6 k.y. and correlates with numerical predictions of stress changes at volcanically active sites. These results support the presence of a causal link between variations in ice age climate, continental stress field, and volcanism.


Volcanic activity varies over a wide range of temporal scales, from cycles with less than one-year periods in single volcanic systems, to intervals extending out to plate tectonic time scales (Cambray and Cadet, 1994; Kennett et al., 1977; Mason et al., 2004; Paterne et al., 1990). Connections between volcanism, the global carbon cycle and climate appear to be well established at the longest of these time scales (Walker et al., 1981). Moreover, on seasonal-to-decadal time scales, volcanic eruptions are known to influence climate (e.g., Hansen et al., 1992; Robock, 2000), and perhaps vice versa (Novell et al., 2006; Rampino et al., 1979). We focus on the intermediate time scales (103–105 yr) relevant to variations during the ice age. There is evidence that subaerial volcanism increased significantly during the last deglaciation phase of the ice age (Huybers and Langmuir, 2009; Jull and McKenzie, 1996), and a connection between climate and volcanism over broader ice age time scales has been established in regional studies (e.g., Jellinek et al., 2004; McGuire et al., 1997; Novell et al., 2006; Paterne et al., 1990). However, a general link between glacial cycles and global volcanic activity has remained elusive.


We analyze marine records of widely dispersed, and well-preserved, tephra layers associated with the Circum-Pacific volcanic chain referred to as the “Ring of Fire” (ROF). The ROF accounts for about half the global length of active plate subduction. Ash plumes of highly explosive Plinian eruptions travel through the atmosphere and deposit volcanic tephras downwind from the eruption sites; sites included in this study lie 200–500 km from their respective sources in the prevailing stratospheric wind directions (Fig. 1).

Our data set is comprised of two subsets. The first involves records from Central America and consists, in part, of 49 eruptions preserved in tephra layers from 56 cross-correlated gravity cores. These data provide an almost complete record of large explosive eruptions (VEI >5 [Volcanic Explosivity Index]) over the past 490 k.y. in the Central American Volcanic Arc (CAVA; Fig. 1). We augment these data with 42 tephra layers extending over ∼1 m.y. found in Deep Sea Drilling Project (DSDP) and Ocean Drilling Program (ODP) Legs offshore of Central America. The marine tephra records in Central America are dated using estimated sedimentation rates and/or through correlation with radiometrically dated on-land deposits (e.g., Kutterolf et al., 2008; also see the GSA Data Repository1).

The second data subset includes DSDP, ODP, and IODP (Integrated Ocean Drilling Program) drill core records from other sites along the ROF (Fig. 1; Table DR1 in the Data Repository). The records extend back to ca. 1 Ma, but rarely include the shallowest sediments with ages up to ca. 100 ka, and have thus been complemented by data from gravity cores. These tephra layers are also dated using estimated sedimentation rates.

The time series of ash deposits at a specific location is based on the individual drill site with the most continuous and complete tephra record. Preference is given to drill sites on the incoming plate since these are generally less affected by erosion than sites on the continental slope. The record is complemented by core data from adjacent sites if possible. Tephra layers of the same age and depositional depth obtained from sites in close proximity (<50 km) are generally counted as one eruptive event. Tephras of similar age recorded at sites at greater distance from one another are counted as separate events. In total, we have identified 408 tephra layers at sites along the ROF. The estimated ages of the tephra layers have uncertainties ranging to 14% of their age. These uncertainties are accounted for in subsequent analyses. For a detailed discussion and explanation of the dating and associated errors see the Data Repository.


We mapped the ROF data set into a binary time series in which the occurrence of a tephra layer is given a value of unity. To display low-frequency variability, we applied a stable-phase, running average low-pass filter with a width of 5 k.y. to the time series (red line in Fig. DR3A in the Data Repository). It is instructive to consider the spectral estimates associated with a sequence of progressively older, overlapping 400 k.y. segments of the tephra time series (Fig. DR4), which we obtain using a multitaper spectral analysis (Thomson, 1982; Percival and Walden, 1998). A peak at a frequency of 1 per 41 k.y., the Milankovitch obliquity frequency, is present in the most recent segment, which is expected to have greatest age accuracy. Examination of older segments shows that the amplitude and frequency of the spectral peak nearest to the obliquity frequency varies considerably in other portions of the time series, possibly because of timing errors that lead to distortion of the spectral estimate (Huybers and Wunsch, 2004).

To assess the statistical significance of the peak nearest to obliquity, we use the full 1.2 m.y. interval of tephra observations. Timing errors are allowed for by searching for the largest peak within ±14% of the 1 per 41 k.y. obliquity frequency, which is identified at a frequency of 1 per 44 k.y. Fitting an exponential function to the spectral estimate between 1 per 10 k.y. and 1 per 80 k.y. indicates that the peak at 1 per 44 k.y. rises 6.3× above the background continuum. The significance of this departure relative to the local background continuum is assessed using a Monte Carlo method, wherein synthetic ash layers are realized according to a Poisson distribution so as to have roughly similar numbers of tephra layers as observed. Spectral analysis is performed on each realization of the synthetic ash sequence, and the peak within 14% of the obliquity frequency that rises furthest above the background continuum in a fractional sense is recorded. To build up a null distribution, this procedure is repeated 100,000 times, and indicates that the 6.3 ratio identified for the actual 1 per 44 k.y. peak would be realized less than 1% of the time, according to our null model. Because no other process is expected to concentrate energy in the band of frequencies that we consider, these results constitute a statistically significant detection of variations in tephra related to the quasi-periodic 1 per 41 k.y. changes in the obliquity of Earth’s spin axis.

To further evaluate the spectral peak that we associate with obliquity, we compress all tephra ages by a factor of 10%. Compression has the largest effect on the absolute ages of the oldest samples, for which uncertainties are the largest. The eruption frequency computed on the basis of this compressed time series is given in Figure 2A. The associated power spectrum of this time series is shown in Figure 2B (and by the black line in Fig. DR3D). In addition to the peak within the obliquity band, the spectrum exhibits some additional concentration of energy near the 1 per 100 k.y. and 1 per 23 k.y. frequencies. Although these concentrations contain no significant peaks, the presence of excess energy is in accord with other indicators of late Pleistocene climate.


The foregoing spectral results strongly suggest that ice-age climate variations induce volcanism, and the physical link may be the changes in crustal stress associated with ice-ocean mass redistributions during the glacial cycles (Glazner et al., 1999; Huybers and Langmuir, 2009; Jull and McKenzie, 1996; McGuire et al., 1997; Sigvaldason et al., 1992). Volcanic output would be expected to increase during periods in which the local confining pressure decreases, lowering the integrity of the enclosing rock. In addition to a depressurization-induced increase in magma production (McKenzie, 1984), glacial unloading also facilitates dike formation and propagation from magma chambers to the surface (Jellinek et al., 2004).

To further explore this connection, we consider phase relationships between tephra layers and variations in ice-age climate. We use a time series of benthic δ18O as a proxy for eustatic sea-level variations or, equivalently, variations in ice volume. Figure 2C shows the power spectrum of the time rate of change of a Pleistocene stack of globally distributed benthic δ18O time series values (Lisiecki and Raymo, 2005). The δ18O rate has three distinct peaks centered at periods of 100 k.y., 41 k.y., and 21 k.y., corresponding to the average time scale associated with glacial/interglacial cycles, the period of obliquity variations, and the average periods of climatic precession (Hays et al., 1976; Berger and Loutre, 1992).

Notably, the obliquity bands contain the greatest concentration of spectral power in both the volcanic tephra records and the rate-of-change of δ18O records. Examination of the coherence and phase relationship within the obliquity band shows that volcanism is significantly coherent with, and slightly lags, obliquity by 24 ± 22° or, equivalently, 2.7 ± 2.5 k.y. (see Fig. 3, and the Data Repository). This relationship can be seen, to some extent, through comparing a band-pass filtered version of the tephra time series to the variations in obliquity (Fig. 2A, inset). As follows from this relationship, and the fact that changes in obliquity are roughly anti-phased with the rate of change of δ18O (Roe, 2006), we find that volcanism peaks roughly when the rate of change of δ18O is most negative, with the latter lagging by 169 ± 9° or 19.3 ± 1.1 k.y. Thus, volcanism peaks 35 ± 31° or 4.0 ± 3.6 k.y. after the greatest rate of sea-level rise and ice volume decline (Fig. 3). This is consistent with the regional study of Jellinek et al. (2004), who found that volcanism lagged glacial unloading by 3.2 ± 4.2 k.y. in eastern California (United States).

To explore how changes in eustatic sea level and glacial loading can influence volcanism, we calculate ice-age-induced, near-surface stresses during the last glacial cycle. We use a standard normal mode treatment of glacial isostatic adjustment (Wu and Peltier, 1982) that is valid for spherically symmetric, Maxwell viscoelastic Earth models (see the Data Repository). The simulations adopt the ICE5G model for the ice geometry during the last glacial cycle and the VM2 radial profile of mantle viscosity (Peltier, 2004). Sea-level changes were computed using a gravitationally self-consistent theory that accounts for the migration of shorelines due to local onlap and offlap and changes in grounded and marine-based ice, and the feedback into sea level of perturbations in the Earth’s rotation vector (Kendall et al., 2005).

As an illustration of the numerical predictions, Figure 4 (top frame; Central American volcanic arc [CAVA]) and Figure DR5 (global) show the rate of change of radial stress at 10 ka, i.e., during the final deglaciation phase of the ice age, at a depth of 20 km. This depth is appropriate for primitive magma reservoirs that feed more-evolved magmatic systems at shallower crustal levels. In oceanic regions within the far-field of the late Pleistocene ice sheets, the radial stress increases by 0.1 MPa in 1 k.y., which is consistent with the addition of ∼10 m/k.y. of meltwater during this time period in the ICE5G model. The predictions show more complexity in the near field of the ice sheets, where ice melting, and the associated isostatic adjustment of the solid Earth, including uplift of previously glaciated areas, and subsidence of the peripheral bulges, contribute significantly to the state of stress.

In the far-field of the ice sheets, continental regions show a smaller amplitude change in radial stress at 10 ka than oceanic regions. Continental interiors, in particular, show a drop in radial stress at this time; this stress drop is due to ocean loading, which causes a tilting upward of the lithosphere in a landward direction, and a migration of mantle material from below the oceans to below the continents. The associated gradient in the radial stress rate is evident in the map of Figure 4, where the stress drop increases in a direction away from the oceans, and peaks in the interior of Central American land masses. We conclude that the ice age stress field is characterized by significant geographic complexity, and that changes in the state of stress at oceanic sites may be strongly out of phase with stress changes on continental margins.

Figure 4 (bottom frame) shows the temporal evolution of the rate of change of the radial stress at 20 km depth below the CAVA (Fig. 1). The CAVA site has the densest and most complete record of volcanic eruptions during the last glacial cycle. The timing of these eruptions is superimposed on the plot. It is clear from the figure that these subaerial volcanoes experience the largest change in crustal normal stresses during the deglaciation phase. During this post-LGM (Last Glacial Maximum) period, normal stresses at 20 km depth are primarily falling. Most notably, there is a significant increase in volcanic activity during this period of heightened ice-age-induced changes in the radial stress (Fig. 4).


We have reported on the first Pacific-wide detection of obliquity-forced variability in late Pleistocene records of volcanic eruptions identified using a large number of well-preserved tephra layers in widely distributed marine sedimentary cores. The detection confirms a connection between variations in climate and volcanism during this time period. We suggest that this connection is mediated by changes in surface mass loading and the associated isostatic adjustment of the solid Earth. This loading mechanism is supported by an observed phase relationship, wherein increased volcanism slightly lags behind glacial unloading at the 41 k.y. obliquity band, as well as a model analysis of changes in normal stress. Our numerical predictions of ice-age-induced stress changes suggest that future efforts to correlate such changes with eruption frequency should take into account the geographic variability in the ice-age stress field. Further development of tephra time series, to include greater age control, longer time spans, and more-detailed spatial coverage, would also help to better characterize the nature of the coupling between climate and volcanism.

We gratefully acknowledge comments by Hans Graf, Hans-Ulrich Schmincke, Rick Murray, and Wendy Perez on earlier versions of the manuscript. This publication is contribution no. 138 of the Sonderforschungsbereich 574 “Volatiles and Fluids in Subduction Zones“ at Kiel University. Mitrovica acknowledges support from Harvard University and the Canadian Institute for Advanced Research. Huybers acknowledges support from the Packard Foundation.

1GSA Data Repository item 2013055, data, additional figures explaining the methods, and an extended error analysis of the data set and statistics, is available online at www.geosociety.org/pubs/ft2013.htm, or on request from editing@geosociety.org or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.