A range of methods have been used to estimate the volumes of volcanoes and individual eruptive units. For whole edifices these methods include finding the area between a pre-volcanic datum and a digital elevation model of the volcano’s surface (often with a GIS program: e.g., Samaniego et al., 2016), and using a 3-dimensional shape to approximate the geometry of the edifice (e.g., Godoy et al., 2018). For individual eruptive units, volume estimation methods include taking field thicknesses and multiplying by the areal distribution of lava flows that are either observable, or buried and inferred (e.g., Mazama, Cascades: Bacon and Lanphere, 2006; Three Sisters, Cascades; Fiertsein et al., 2011; Calvert et al., 2018; Tongariro, New Zealand; Pure et al., 2020). Methods for individual eruptive units also include subdividing sectors of a volcano using 3-dimensional shape approximations to estimate volume, based on where each major eruptive sequence occurs (e.g., Parinacota, Chile: Hora et al., 2007; Ruapehu, New Zealand; Conway et al., 2016).

To accurately reconstruct the long-term shifts in eruptive rates at a Pleistocene-age volcano, the volumes of lost materials need to be estimated. In undertaking this exercise, the most common approaches are (1) to assume that each eruptive vent had a perfect cone geometry (Thouret et al., 2001) and/or (2) to assume that valleys were previously filled with erupted deposits and were subsequently eroded (Singer et al., 1997). However, sometimes the surviving edifice volume is simply assumed to represent the total amount of erupted material, even where glacial unconformities are reported (e.g., Samaniego et al., 2012; Samaniego et al., 2016), which may be inaccurate at volcanoes that have experienced glacial advances and sector collapses.

Approaches (1) and (2) noted above are justified for vents and flanks that grew in ice-free environments (e.g., Figure 6C). In such cases, the existence of laterally excavated sequences of planar lava flows (e.g., Figure 4A; green units in Figure 6C) and valley-bounding moraines (e.g., see “m” in Figures 4A, C, D, F) may be used to infer distinct periods of edifice growth and erosion. In applying the ice-free growth model, some studies have estimated the amount of denudation by glaciers within valleys to infer eroded volumes (e.g., Mixon et al., 2021). Using this framework, the estimated volume of eroded material at Paniri (Chile) is 4.4 km³, relative to 79.9 km³ of surviving edifice material (Godoy et al., 2018). Elsewhere, estimates of the down-cutting rates from glacial erosion have been estimated by assuming that lava ridges perched above glaciers are remnants surviving after glacial erosion on surrounding pre-existing lava. Using this approach, Hildreth et al. (2003b) inferred down-cutting rates of 0.4–1.75 m/kyr at Mount Baker. However, lava ridges may also be primary (constructional) landforms rather than secondary (erosional) landforms because many glaciated volcanoes demonstrably grew in the presence of ice (e.g., Conway et al., 2016). Furthermore, although valleys at most volcanoes were occupied by glaciers, in many cases it is difficult to differentiate whether such valleys were built around ice or were eroded by ice, or both (e.g., Singer et al., 1997; Singer et al., 2008).

The estimated volumes of Holocene lava flows and consolidated pyroclastic deposits are relatively likely to represent the actual amount of erupted magma because, in most cases, there has been minimal opportunity for erosion to occur via glaciation or through sector collapse events. Evidence for this includes lavas and pyroclastic deposits that have intact carapaces and no glacial striations, and are fully traceable to their source vents, as illustrated by the dark red units in Figure 6C (e.g., the Ngāuruhoe cone of Tongariro, New Zealand: Hobden et al., 1999). In the example shown in Figure 6, the volcano may have had a steady eruptive rate over 200 kyr, but the volume of syn-glacial eruptive products may be under-represented compared with post-glacial products, due to the more complete preservation of the latter.

Tephras and unconsolidated pyroclastic deposits are much less likely to be preserved on a volcano’s edifice over its long-term record than lavas and some pyroclastic deposits. A study of tephra record completeness in ice-free regions suggests that <1% of products from eruptions with a volcanic explosivity index of ≥2 have been preserved in post-LGM stratigraphic records, although the total volume of preserved, erupted material may be >1% (Watt et al., 2013). Unconsolidated pyroclastic deposits that are not capped by lavas are rare in the edifice records of most glaciated volcanoes (e.g., Hackett and Houghton, 1989). Typically, explosive eruption deposits on volcanoes comprise the minority of preserved edifice material (e.g., p. 732 of Hildreth et al., 2003b; p. 364 of Fierstein and Hildreth, 2008; p. 405 of Martínez et al., 2018). Collectively, these observations show that edifice records on volcanoes are biased towards preserving lavas and pyroclastic deposits capped by lavas, relative to tephras and uncapped pyroclastic deposits.

Although explosive eruptions can represent a significant proportion of a given volcano’s eruptive output, the pre-LGM record is usually poor and only preserves the largest-volume events (Kiyosugi et al., 2015). In regions affected by continental glaciation, the ring plain tephra record is often obliterated by glacial erosion (e.g., Hildreth, 1996). Accurate reconstructions of past eruptive rates at glaciated volcanoes may only be possible if ring plain records, which comprise lavas and pyroclastic debris (eroded via fluvial or glacial action), tephras, and sector collapse deposits, can be correlated with specific source volcanoes based on radiometric dating, compositional affinity, drainage patterns and detailed stratigraphic mapping (e.g., Tost and Cronin, 2015). Such exercises are extremely resource- and time-intensive and are also undermined by long run-outs of some sector collapse debris flows that can travel hundreds of kilometres from their source volcano (Vallance and Scott, 1997; Waythomas and Wallace, 2002; Tost and Cronin, 2015).

There are also other biases that appear to affect the preservation of erupted materials in the edifice record. The most notable is when eruptives are emplaced onto ice masses but never make contact with bedrock, which prevents them from being preserved because they are flushed away as debris during deglaciation. This has been inferred from eruptive records (e.g., Jicha and Singer, 2006; Conway et al., 2016; Pure et al., 2020) and observed for modern eruptions (e.g., Edwards et al., 2015; Loewen et al., 2021).

The summary in Table 2 shows that of 33 studies, 20 considered lava-ice interaction and 10 of these studies considered or reported evidence of supra-glacial emplacement of erupted materials. Notably, the only studies that considered supra-glacial emplacement were those that also reported evidence or likelihood of lava-ice interaction (Table 2: Akutan, Katmai, Middle Sister, Baker, Rainier, St. Helens, Ruapehu, and Tongariro volcanoes). No studies of Andean volcanoes from Ecuador, Peru, Argentina or Chile considered supra-glacial emplacement of erupted materials, even though the ages of these volcanoes and extents of ice coverage would have made supra-glacial emplacement of eruptives highly likely (Table 2). Lava flows that were emplaced in ice-free valleys during past interglacial periods may also have low preservation and/or exposure potential, because they can be eroded in subsequent glacial advances and buried by younger till and volcanic deposits (e.g., Conway et al., 2016; Figure 4C).

Post-glacial sector collapses of asymmetric landforms that grew in the presence of ice have been reported widely in both northern and southern hemispheres. Capra (2006) identified over 20 post-LGM sector collapse events from Argentina, Bolivia, Chile, Ecuador, Mexico, New Zealand, and the United States, which is not a comprehensive list (cf. Osceola Mudflow at Rainier, Cascades: Vallance and Scott, 1997; Murimotu Formation at Ruapehu, New Zealand; Palmer and Neall, 1989). Such sector collapse events remove edifice material of a variety of ages depending on which materials were part of the collapsed landform. The existence of hydrothermally altered material in sector collapse deposits, which may form in subglacial environments via interactions with meltwater, show that subglacial materials are commonly removed from the edifice in collapse events (e.g., Vallance and Scott, 1997; Waythomas and Wallace, 2002). Subglacial and syn-glacial eruptives are thus likely to have low preservation potential because of their propensity to be fragmented, altered, buried, collapsed, and eroded. The ability to accurately reconstruct eruptive rates in materials that pre-date sector collapse events may be hampered by mass wasting that removes significant volumes of older material from the edifice, which could therefore lead to underestimates of pre-collapse (i.e., syn-glacial) eruptive rates.

Estimates of erupted volumes are a major source of uncertainty in studies that seek to address the question: does deglaciation cause increased eruptive rates? There is then a marked contrast between studies that interpret variations in time-volume trends as arising from erosion and eruptive quiescence (e.g., Singer et al., 1997; Singer et al., 2008) and other studies that infer a preservation bias whereby erupted materials were flushed off the edifice to the ring plain (e.g., Pure et al., 2020). In many cases, field records provide strong evidence that 1) glacial and fluvial erosion occurs, 2) greater volumes of erupted materials are preserved in periods of reduced ice coverage because of a preservation bias, 3) volcanoes erupt both explosively and effusively, and 4) Pleistocene-age volcanoes contain significant volumes of unobservable buried material. On the contrary, field records struggle to accurately show 1) how much erosion has occurred because assumptions about pre-erosional geometry cannot be verified, 2) what volume of material that was erupted syn-glacially has been flushed to the ring plain because clast provenance can be difficult to establish and continental glaciation can obliterate ring plain records, 3) how much volume has been erupted effusively versus explosively because different lithologies have different likelihoods of being preserved, and 4) the ages and provenance of completely buried edifice materials.

In Figure 3, volcanic growth histories are compared against the stacked benthic δ¹⁸O curve from Lisiecki and Raymo (2005), which largely reflects global ice volume, to assess whether eruption rates at arc volcanoes have been influenced by climate change. Within the reviewed papers, the distinctions between reduced and advanced terrestrial ice cover have generally been taken as those that define each Marine Isotope Stage (MIS). These are sometimes integrated with regional paleoclimate information to confirm the timings of glacial advance and retreat, however, few detailed glacier reconstructions exist for arc volcanoes (cf. Eaves et al., 2016; Eaves et al., 2019). When testing whether sea-level changes have affected volcanism (e.g., Wallmann et al., 1988), it can be assumed that rates of lithostatic (un)loading are constant at the scale of an individual volcano or region. On the other hand, ice advance or retreat can occur at different rates within separate valleys of the same volcano due to local morphological and climatological effects (Eaves et al., 2019). The following discussion focusses on the general timing of ice retreat during the Last Termination for the purpose of comparing that timeframe to the precision of eruption age measurements for lava flows, because links between long-term eruptive rates and small-scale glacier mass fluctuations are not well-studied.

The rapid rise in sea-level from 18 to 11.7 ka reflects the wholesale melting of ice sheets and glaciers during the Last Glacial Termination (Denton et al., 2010), which is supported by the stacked record of benthic foraminiferal δ¹⁸O compositions (Lisiecki and Raymo, 2005). For the sea-level record, the average 2σ error on corrected coral ages calculated from published U and Th isotope ratio data over that period is ±0.7 ka (Thompson and Goldstein, 2006). On land, moraine records at or near the volcanoes assessed in this review reflect a period of approximately 10 kyr for the retreat of glaciers from LGM to near-historical extents. Davies et al. (2020) summarized that the Patagonian Ice Sheet in South America started to retreat from greatest extent at ∼28 ka and then stabilised at 21–18 ka. Subsequently, the ice sheet underwent rapid deglaciation and separated into disconnected ice masses by 15 ka, followed by readvances or stabilisations of glaciers at 14–13 ka and 11 ka. These conclusions are supported by ¹⁰Be exposure ages from moraines near Lago Palena on the border of Chile and Argentina (43.9° S; 71.5° W), which have 2σ errors on mean moraine ages of 0.8–1.4 ka (Sotores et al., 2022). This level of precision is similar to that (2σ values of 0.4–1.8 ka) reported for ³He exposure dating of LGM and late-glacial moraines at Tongariro and Ruapehu volcanoes in New Zealand reported by Eaves et al. (2016), Eaves et al. (2019).

The ⁴⁰Ar/³⁹Ar and K-Ar geochronometers have been applied extensively to define the eruptive histories of Quaternary volcanoes worldwide, and they are vital tools for constraining ages for late Pleistocene arc volcanic rocks. One or both of the dating methods have been used for each of the studies that are covered in this review. They are best suited to dating the microcrystalline (non-glassy) groundmass material of lava flows (Gamble et al., 2003; Fleck et al., 2014; Calvert et al., 2018; Sisson et al., 2019). In particular, lavas that were emplaced subaerially next to glaciers provide ideal products for dating because they pond into thick flows with slowly cooled interiors that produce coarsely crystalline groundmass textures (groundmass microlites >20 μm; Conway et al., 2016). In contrast, subglacial volcanism can produce glassy (and hydrothermally altered) volcanic products that are difficult to date using these methods (Flude et al., 2010; Guillou et al., 2010). Eruptive activity during interglacial stages at volcanoes may form rubbly, blocky, or thin lava flows over steep flanks or within ice-free valleys that are also characterized by glassy groundmass textures (Pure et al., 2020). Low radiogenic ⁴⁰Ar yields and high atmospheric ⁴⁰Ar in low-K and/or young (<20 ka) lavas also presents a major challenge for defining the recent eruption histories at arc volcanoes (Jicha, 2009).

The precision of geochronological data in studies assessed in this review is compared with the duration of the Last Termination in Figure 7. Out of 44 ages between 30 and 10 ka, five have 2σ values of 6%–9% relative to the age (i.e., 2σ values of 1.1–2.2 ka). A further ten age measurements have 2σ values of >9–20% relative to the preferred age. Sixteen lavas have 2σ values between 21% and 42% relative to the preferred age, and the remaining thirteen have 2σ values of 52%–93%. An age of 17.8 ± 2.2 ka reported by Conway et al. (2016) indicates, with 95% confidence, that the lava erupted between 20.0 and 15.6 ka. Despite being a relatively precise age, this ± 2σ range spans the period from full glacial conditions to major ice retreat at Ruapehu volcano. The 17.8 ± 2.2 ka age is however distinct from a re-advance at ∼14–11 ka corresponding to the Antarctic Cold Reversal (Eaves et al., 2019). At Calbuco, an age of 11.5 ± 8.0 ka reported by Mixon et al. (2021) spans the period from LGM extent of the Patagonian Ice Sheet at 19.5 ka, through rapid retreat at 17 ka, to re-advances of glaciers at 14–13 ka, 11 ka, and 6–5 ka (Moreno et al., 2015). Combined with field observations, these age constraints are suitable for broadly identifying syn-versus post-glacial eruptive stages, but cannot be directly compared with climatic events on sub-kyr-scales (Figure 7). As such, it is difficult with currently available analytical precision to test the sensitivity of the volcanic response to deglaciation. For example, it would be valuable to know whether eruptions are initiated during large-scale glacier retreat, when a threshold of 70% of ice was lost, or 3 kyr after near-modern glacier extents were reached. Similarly, more precise age data may allow testing of whether minor stabilizations or re-advances of glaciers paused or slowed any potential increases in eruptive rates after the LGM. This has significant implications for understanding whether future glacial retreat under current warming will lead to increased rates of volcanism, and our ability to provide meaningful forecasts for such volcanic activity.

Calibrated ¹⁴C ages for organic material underlying, overlying, or included in layers of volcanic material can provide precise geochronological constraints. Analytical uncertainties of ±100 years (2σ) are regularly reported (e.g., Bromley et al., 2019), which are significantly smaller than uncertainties on typical ⁴⁰Ar/³⁹Ar and K-Ar age determinations, as summarised above. Radiocarbon dating has been successfully used to establish high-resolution eruptive records from post-glacial tephras in studies that investigate potential causality between deglaciation and apparent eruption rates (Watt et al., 2013; Rawson et al., 2016). Radiocarbon dating is not always feasible, however, particularly for eruptives with ages >50 ka and on the upper flanks of volcanoes in middle to high latitudes that have little organic material interbedded or overlain by lavas and pyroclastic units.

Neither eruptive episodes nor deglaciation are instantaneous events, and establishing causality between them has been, and will continue to be, a challenging endeavour. The limits of analytical precision with the K-Ar and ⁴⁰Ar/³⁹Ar geochronometers provide age estimates within 1–2 kyr ranges for medium-to high-K andesite lavas of Late Glacial age, but often 10 kyr or more for low-K lavas. Such uncertainties hinder meaningful assessments of whether variations in eruption rates were caused by deglaciation, and prevent the identification and measurement of time gaps (i.e., “lag times”; Figure 1C) that may occur between deglaciation and episodes of heightened eruptive activity.

Of the 20 volcanoes shown in Figure 3, 12 show post-LGM (∼20 ka) increases in edifice growth rates (Tables 2, 3). However, only 7 of these 12 volcanoes show increased edifice growth rates after the MIS 6 deglaciation (∼130 ka), which was the most substantial deglaciation event in the ∼150 kyr preceding the LGM (Singer et al., 2000; Kaplan et al., 2004). Furthermore, of the 7 volcanoes that show post-glacial increases in edifice growth rates after MIS 6 and the LGM, 3 volcanoes show relative increases in edifice growth rates in the MIS 4–2 period as well, which is when ice coverage was generally increasing into the MIS 2 glaciation. Thus only 4 of 20 volcanoes examined (Mazama, Puyehue-Cordón Caulle, Antuco-Sierra Velluda, and Tongariro) show repeated increases in apparent eruptive rates following the MIS 6 glaciation and the LGM without contradictory increased growth rates in the MIS 4–2 period (Bacon and Lanphere, 2006; Singer et al., 2008; Martínez et al., 2018; Pure et al., 2020).

At Mazama, the increased eruptive rate in the Holocene is largely due to a large caldera-forming eruption that ejected ∼50 km³ of magma (Bacon and Lanphere, 2006; Figure 3). Eruptions of this size from Quaternary arc volcanoes are rare in the Cascade arc, with only three caldera-forming eruptions having occurred since 1.15 Ma (Hildreth, 2007). On correlations alone, it appears unlikely that the timing of the caldera-forming Crater Lake eruption at Mazama was triggered by deglaciation, noting also that the Kulshan caldera-forming eruption at 1.15 Ma was subglacial (Hildreth, 1996). In comparison, other Quaternary arcs that experienced similar amounts of ice-loading during Pleistocene glaciations have hosted more caldera-forming eruptions than the Cascade arc, such as eight or more in central Chile and >20 in the Alaska-Aleutian arc (Miller and Smith, 1987; Newhall and Dzurisin, 1988).

In the study of Puyehue-Cordón Caulle, the authors interpreted post-glacial increases in edifice growth rates to be the consequence of ice unloading of the crust, which promoted greater melt ascent from the deep crust than in pre-glacial and glacial periods (Singer et al., 2008). However, time versus compositional data from erupted lavas do not provide convincing support for this conclusion because neither edifice growth rates nor δ¹⁸O climate proxy data are correlated with changes in the compositions of erupted lavas. The authors also acknowledge the influence of sector collapses and glacial erosion in modifying the amount of preserved edifice material at Puyehue-Cordón Caulle, which questions whether volcanic records are sufficiently representative of Pleistocene and Holocene eruptive rates.

Eruptive records from Antuco-Sierra Velluda have been interpreted to show potential increases in eruptive rates since ∼150 ka, or after the MIS 6 glaciation, with a small decrease in edifice growth rates during 60–17 ka (∼MIS 3–2), which the authors consider may be due to less severe glaciation during MIS 4 than in MIS 6 (Martínez et al., 2018; Figure 3). The growth rates since ∼150 ka, however, are relatively uniform, and the variations in edifice growth rates since ∼150 ka may be due to supra-glacial emplacement of eruptives in MIS 4–2.

At Tongariro, time-volume patterns were interpreted to arise from a preservation bias whereby lavas and pyroclastic deposits were more commonly emplaced onto ice during glacial periods that were subsequently transported to the ring plain as debris upon deglaciation (Pure et al., 2020). Systematic fluctuations in MgO concentrations through time, however, are uncorrelated with climatic cycles. Thus, the correlations between ice coverage and edifice growth rates at Tongariro are unlikely to accurately represent variable eruptive rates that were modulated by ice loading and unloading of the crust. Taken together, the time-volume-composition data from Mazama, Puyehue-Cordón Caulle, Antuco-Sierra Velluda, and Tongariro provide limited support for the idea that deglaciation causes increased eruptive rates at arc volcanoes.

A number of studies have examined potential volcanic responses to deglaciation over multiple glacial-interglacial cycles (e.g., Singer et al., 1997; Edwards et al., 2002; Hildreth et al., 2003b; Lanphere and Sisson, 2003; Bacon and Lanphere, 2006; Jicha and Singer, 2006; Singer et al., 2008; Schmidt and Grunder, 2009; Jicha et al., 2012; Calvert et al., 2014; Conway et al., 2016; Calvert, 2019; Coombs and Jicha, 2020; Pure et al., 2020; Santamaría et al., 2022) whereas other studies only considered or have data for edifice growth rates during and after the last glacial period (e.g., Fiertsein et al., 2011; Calvert et al., 2018). Of the studies that solely consider the MIS 3–1 interval, some show maximum edifice growth rates during the glacial period when ice coverage was greatest (e.g., Middle Sister: Calvert et al., 2018) whereas others show maximum apparent eruptive rates after the LGM (e.g., Calbuco, Chile: Mixon et al., 2021), as in Figure 3. For the former examples with relatively high eruption rates during MIS 3–2, the ages and distributions of ice-affected volcanic products can provide valuable information about the past extents of glaciers on their edifices. Although moraines record outer ice limits during advances, ice-bounded lavas that erupted independently of ice dynamics can provide unique constraints on the presence and thickness of glaciers throughout stages of prior glacial advance and retreat (e.g., Lescinsky and Fink, 2000; Edwards et al., 2002; Conway et al., 2015). Such data from glaciovolcanic deposits can be used to constrain paleoclimate models.

Several studies that examined multiple glacial-interglacial cycles have concluded that the clustering of radiometric ages into interglacial periods does not equate to increases in eruptive rates because time-volume trends do not show post-glacial increases in edifice growth rates (e.g., Hildreth et al., 2003b). Such observations contrast with approaches elsewhere that use the statistical distribution of radiometric ages as a proxy for edifice growth rates in the absence of detailed time-volume data (Calvert et al., 2014; Calvert, 2019). In these cases, although the number of age determinations in interglacial periods is often higher than in glacial periods, as at Mount Rainier, this correlation is not supported by time-volume records from edifice-forming lavas and pyroclastic deposits, because the apparent correlation is an artifact of non-deposition and non-preservation during the periods of abundant ice coverage (Calvert et al., 2014; Sisson and Calvert, 2023).

Preservation and erosion act to reduce the volume of older edifice material, relative to younger edifice material. A potential consequence of this is that reconstructed time-volume data may show ever-increasing edifice growth rates through time. Examples of this pattern are Baker, Adams, and Parinacota (Figure 3), which have upward-facing curves with broadly increasing km³/kyr rates through time (Hildreth and Lanphere, 1994; Hildreth et al., 2003b; Hora et al., 2007). If arguments are made that apparent eruptive rates increase after the LGM because of deglaciation and magma system depressurization, then the same test should be applied to earlier periods of the volcano’s history. For Baker, Adams, and Parinacota, apparent eruptive rates increase during the last glacial period (i.e., ∼MIS 4–2) relative to MIS 5, which contradicts the conclusion that deglaciation is the cause of increased eruptive rates.

Using a hybrid dataset of tephra records from 40 ka to present, Huybers and Langmuir (2009) reported that subaerial volcanism increased 2 to 6 times against background levels after the post-LGM deglaciation. The authors interpreted the feature to result from mantle decompression following deglaciation and then predicted that sea level rise would suppress eruption rates at mid-oceanic ridges. We note, however, that the tephra dataset reported Bryson et al. (2006), which was used for the analysis of Huybers and Langmuir (2009), acknowledges that terrestrial tephra records may be incomplete because of irregular preservation.

A comparable but more cautionary result was reported in a separate assessment of tephra records from the Andean southern volcanic zone, the Cascades, and Kamchatka by Watt et al. (2013). This study noted that inferred eruptive rates from the LGM to the late-Holocene varied by no more than a factor of two, and time-volume variations were not statistically significant. Watt et al. (2013) also reported that the extent of tephra record incompleteness in post-LGM records is substantial, and may be as little as 0.005% (1 in 20,000 events preserved), and that under-sampling of small (<0.1 km³) eruptions prevents meaningful quantification of LGM vs. post-LGM eruption rates. For these reasons and those outlined in the ‘Biased preservation and exposure of eruptive products’ section above, it is unclear whether tephra records in the Holocene are sufficiently complete or free of preservation bias to make interpretations about whether global volcanic activity levels increased after the LGM.

Out of 20 time-volume datasets from arc volcanoes affected by glaciation (Figure 3), 4 show reproducible increases in edifice growth rates after deglaciation events and 12 provide unclear evidence for enhanced post-glacial volcanic activity. Interpretations for why edifice growth rates increase after deglaciation events vary widely and do not always consider eroded or unpreserved materials. For these reasons, periods of heightened erosion or reduced preservation may be misinterpreted as eruptive hiatuses. Some edifice records show progressive increases in growth rates through time (e.g., Baker, Adams, and Parinacota in Figure 3), but it is difficult to determine whether such trends are real or if they are artefacts of erosion, mass-wasting events, and burial; each of these phenomena may disproportionately reduce the exposed volumes of older edifice materials relative to younger materials. Terrestrial tephra records are likely to be less useful than lavas for testing the reproducibility of potential volcanism-deglaciation feedbacks because pre-Holocene tephra records are highly incomplete (Watt et al., 2013). Even for the last glacial cycle alone, it is disputed whether tephra records provide credible evidence for increased post-glacial eruptive rates because volcanic records are incomplete and time-volume variations are statistically insignificant (e.g., Huybers and Langmuir, 2009; Watt et al., 2013; Rawson et al., 2016). The response of volcanoes to deglaciation is not expected to be uniform in time or space, considering the differences in ice volume for each glacial stage and each arc. The response is also likely to differ based on which stage of a volcano’s magma production history is intersected by any given deglaciation event. Therefore, although reproducible correlations between increased ice volumes and reduced eruption rates, across multiple glacial-interglacial cycles, would therefore support a causal relationship, the opposite is untrue. That is, increases in ice volume that are uncorrelated with changes in eruptive rate would not support a causal connection. Furthermore, if increased edifice growth rates in the Holocene are cited as evidence for deglaciation causing increased eruptive rates, then this same standard of evidence should apply to previous glacial-interglacial cycles. These points may be useful for (re)considering whether feedback mechanisms between climate and volcanism often or seldom impacted eruption rates throughout the past and whether many or few volcanoes will be “primed” to respond to future ice retreat.

Several studies have modelled crustal stress responses to deglaciation. Jellinek et al. (2004) estimated that a ∼1–3 MPa decrease in loading pressure during deglaciation of eastern California may have caused stalled magma to ascend, and earlier estimates for silicic magmas were higher (∼10–30 MPa decrease required for magma ascent: Jellinek and DePaolo, 2003). In the southern Chilean Andes at Mocho-Choshuenco, Rawson et al. (2016) estimated that loading pressure on crustal magmas at ∼16 km depth decreased by ∼6 MPa when continental glaciers melted following the LGM. Martínez et al. (2018) also inferred that a stress decrease of ∼3–5 MPa following deglaciation that resulted in a basaltic eruption at 17 ka after the LGM at Antuco-Sierra Velluda. Together, these estimates indicate that pressure changes of about <10 MPa are linked to apparent changes in eruptive rates in arc systems.

Estimates of glaciostatic loads on the magmatic system of a representative, glacierized, arc volcano are given in Table 4, based on the conditions for Tongariro volcano during MIS 4 (Eaves et al., 2016; Cole et al., 2018). For a 70 km³ volcano that is 20% covered by glaciers (200 m thick, in valleys) and 10% covered by a summit ice cap (200 m thick: e.g., Smellie and Skilling, 1994; Magnússon et al., 2012; Cole et al., 2018), the ice forms 7% of the mass load on the pre-volcanic datum, with the remainder coming from the volcano itself. The lithostatic load on a magmatic system, which includes the mass of a 2,000 m-high volcano, increases with its depth in the crust: 151 MPa at 5 km depth, 259 MPa at 10 km depth, and 475 MPa at 20 km depth. When an ice cap and glaciers of uniform thickness cover the volcano, the ice adds a further 1.0–9.8 MPa of pressure onto the magmatic system for ice thicknesses of 100–1,000 m, respectively. For magma systems at 5–20 km depth, this ice contributes 0.2%–5.7% of the total pressure on the magmatic system. The volcanic responses predicted by Jellinek et al. (2004) and Rawson et al. (2016), above, correspond to the pressure decreases that could be caused by removing ∼100–600 m of ice from a glaciated volcano. As a percentage, these changes in glaciostatic loading amount to <3% of the total pressure on a typical continental arc magmatic system (Table 4).

The hypothesis in Figure 1 states that deglaciation leads to increased eruption rates and that glaciation supresses eruption rates. However, testing causality between deglaciation and increased eruptive rates is challenging because major sector collapses can occur during or shortly after glacial retreat (Capra, 2006). It is therefore essential to determine the relative influence of glacial retreat versus sector collapse on magma systems when interpreting time-volume data for volcanoes because both phenomena cause magma system depressurization. The opposite situation, whereby rock mass is added to a volcano throughout its lifetime, should also be considered alongside glacial ablation when interpreting changes in apparent eruptive rates because these processes may have cancelling effects on magma system pressurization, especially in cases of rapid edifice growth (e.g., Pinel and Jaupart, 2000).

As noted previously (see ‘Biased preservation and exposure of eruptive products’ section), Capra (2006) surmised that numerous sector collapses at ice-capped volcanoes were induced by ice retreat, due to glacial de-buttressing, load discharge and fluid circulation as ice retreated. There are likely to be many more sector collapses associated with deglaciation than the 23 examples cited in that review, including the early and middle Holocene collapses of Ruapehu volcano (Palmer and Neall, 1989; Donoghue and Neall, 2001), and the Osceola Mudflow collapse event at Mount Rainier (Vallance and Scott, 1997). Such sector collapses removed rock (and remaining ice) masses from the summits and flanks of many volcanoes in catastrophic events. Sector collapse events are therefore relevant to interpreting time-volume-composition trends because collapse events remove in situ portions of the eruptive record (discussed in the “Biased preservation and exposure of eruptive products” section, above), and because collapse events immediately depressurize magma plumbing systems concurrently with gradual depressurization caused by deglaciation.

The response of magmatic systems to sector collapse events has been the topic of numerous studies, and the availability of high-resolution datasets for historical and pre-historical examples has inspired several detailed reviews (e.g., Watt, 2019; Romero et al., 2021). Watt (2019) noted that voluminous effusive eruptions of magma with anomalous compositions have often occurred following large sector collapses at arc volcanoes. Short-lived episodes of more frequent mafic eruptions support the idea that edifice unloading facilitates the ascent of dense, mafic magmas. This presents a challenge for interpreting time-composition-volume trends at deglaciated arc volcanoes: sector collapse events and deglaciation can both occur in the same several-millennia period, but both phenomena lead to magma system depressurization. If changes in time-composition-volume trends are related to magma depressurization, then it may be unclear if these trends were driven by sector collapse (over minutes) or deglaciation (over millennia). Whether magma system pressurization and associated changes in time-composition-volume trends were caused by deglaciation versus sector collapse will require high-resolution geochronological constraints that are supported by field evidence. Ultimately, however, if magma systems respond to sector collapse events that were triggered by ice retreat, then glaciation imparts a fundamental control on the evolution of eruptive activity at such volcanoes. This is pertinent for considering future hazards at glacierized edifices that will experience ice retreat in the future. Although magma systems may not directly respond to ice retreat linked to climate change, sector collapses and associated eruptive activity will pose significant threats in the future.

When interpreting time-volume-composition trends at glaciated arc volcanoes (Figure 3), it is pertinent to note that volcanic regions where glaciation has had no (or minor) impact also display variability in the rates of volcanism and geochemistry of erupted magma on time scales of 1–100 kyr. For example, thin glaciers may have been present only on a few volcanoes with summit elevations above 2,000 m a.s.l. during MIS 2 in Japan (Sawagaki et al., 2004), and thus deglaciation is not considered to have impacted their eruption rates. In a review of the volumetric and geochemical evolution of 29 arc volcanoes in Japan, Yamamoto et al. (2018) showed that renewed input of mafic magma and instances of crustal melting led to increased eruptive rates at several volcanoes. In contrast, volcanoes that displayed minor variability in erupted magma compositions were shown to be generally characterized by decreasing eruption rates.

None of the 33 volcanoes and volcanic systems covered by our review have exhibited any reproducible correlations between erupted magma composition and ice coverage. In keeping with this, several studies additional to those cited in Table 1 have further explained non-systematic trends for the time-sequenced compositions of erupted material at these volcanoes to be the result of complex open-system processes that operate in crustal magma storage regions (TDungan et al., 2001; Gamble et al., 2003; Mangler et al., 2021). There are very studies that have attributed changes in the composition of post-glacial eruptive products to depressurization via ice retreat. The most commonly cited of these is the case study of Rawson et al. (2016), who identified post-glacial stages of evacuation (eruption of andesite to rhyolite magma), relaxation (basaltic andesite to andesite) and recovery (basaltic andesite to dacite) in the period since 12 ka at Mocho-Choshuenco volcano in Chile. In contrast, Mangler et al. (2022) showed that the sizes, styles, and compositions of eruptions at Popocatépetl since 14 ka can be explained by magma recharge and mixing processes. Conway et al. (2018) showed that a distinct trend towards the eruption of more mafic lavas at Ruapehu commenced at ∼26 ka, prior to deglaciation, and thus was not caused by ice unloading of the magma system. The trend was surmised to reflect changes to the composition of assimilated crustal melts in the plumbing system. The eruption of high-Mg basaltic andesites at ∼45–40 ka at Ruapehu also contradicts the idea of suppressing denser magmas during glacial stages (Conway et al., 2020).

The long-term flux of magma from the mantle to the surface affects the thermal and chemical evolution of the crust that hosts magma systems beneath volcanoes in both non-glaciated and glaciated arcs (Sisson et al., 2013). It is well-documented that this internal evolution (as opposed to the external, surficial evolution of varied ice/edifice load over time) impacts the timescales of magma residence, volumes of magma reservoirs, extents of crystal fractionation, and compositions of crustal assimilants (e.g., Bacon and Lanphere, 2006). It is therefore valid to test alternative hypotheses to the one shown in Figure 1, such as those that explain the evolution of crustal magma systems in non-glaciated regions, when interpreting time-volume and time-composition trends at glaciated volcanoes.

Flare-ups of increased magma production at circum-Pacific arcs have occurred in 5–10 Myr-long pulses since 50 Ma, which has been linked to major tectonic plate reorganizations (Jicha et al., 2009). In addition to variations in the rate, composition, temperature and angle of subducting slabs (e.g., Gill, 1981), mantle processes can also influence the productivity of arc magma systems, even over periods of months to years (e.g., Ruprecht and Plank, 2013). Further examples of mantle and slab control on eruptive activity and magmatic production rates follow.

In young lavas from the Alaska-Aleutian arc, George et al. (2003) reported positive correlations between convergence rate, volcano volume, and ²³⁸U excesses that were interpreted to show that the volume of magma produced depends on the size of the fluid flux supplied to the wedge <10 kyr prior to eruption. Evidence presented by Taniuchi et al. (2020) further supports this conclusion by showing that a ∼7 kyr-period of relatively high eruptive output at Rishiri volcano (northwestern Japan) was driven by the production of hydrous primary magmas via aqueous fluid-fluxed melting of mantle peridotite. Geospeedometry data have also provided evidence for fast migration of mantle-derived magmas to the surface. Ruprecht and Plank (2013) reported Ni diffusion time scales for olivine of weeks to months that indicate magma ascent rates of 55–80 m per day through the ∼35 km-thick crust beneath Irazú volcano in Costa Rica. Fe-Mg interdiffusion time scales for orthopyroxene phenocrysts in high-Mg andesite and dacite lavas from Ruapehu volcano indicate mantle-derived magmas hybridized with mid-crustal felsic magmas less than 10 days before voluminous lava eruptions during MIS 3 (Conway et al., 2020).

The eruption of high-Mg andesite at volcanoes in the northern Cascade arc also indicates slab and mantle processes have a fundamental control on short-term episodes of mafic volcanism with no apparent linkage to deglaciation (Sas et al., 2017). Wall et al. (2018) noted that the waxing and waning eruptive output of andesitic volcanoes in the Cascades over the Pleistocene reflects variability in the thermal energy delivered from the mantle at subduction zones, rather than cyclic surficial loading and unloading imposed by climate variability. Variability in magma compositions and eruptive rates along the Cascades volcanic arc was also shown by Till et al. (2019) to predominantly be related to variations in the flux of basalt into the crust.

These examples from arc volcanoes indicate that slab and mantle processes influence eruption rates and the compositions of erupted materials, and that such processes can cause variability on time scales comparable with global climate cycles. Such alternative phenomena should be considered when evaluating whether deglaciation has caused increased eruption rates or variations in erupted magma compositions, or whether such trends were fundamentally controlled by crustal and/or mantle and slab processes that have not been affected by lithostatic or glaciostatic loading changes.

A number of alternative explanations, beyond the deglaciation hypothesis in Figure 1, should be considered and tested before trends in an arc volcano’s time-volume-composition evolution are interpreted to result from ice loading and unloading. The timescales of magma production and transfer through the crust range from several millennia to several days, which shows that processes other than glaciation and deglaciation can alter eruption rates and magma compositions throughout the life cycles of arc volcanoes.