Maximillian VAN WYK de VRIES

Department of Earth Sciences, University of Minnesota, Twin Cities campus, Minneapolis 55455, MN, United States

Abstract Mt Erebus is the most active Antarctic volcano, on the flanks of the world’s largest ice sheet. Despite this, the interactions between its eruptions and the ice cover have not been studied in detail. Focusing on the most recent deglaciation, we build a glacial retreat model and compare this to recent lava geochemistry measurements to investigate the processes involved.This analysis exposes a previously unknown link between Antarctic glaciation and eruptions, of vital importance to the understanding of volcanism in this context. We find that deglaciation led to rapid emptying of the shallow magma plumbing system and a resulting peak in eruption rates synchronous with ice retreat. We also find that the present day lavas do not represent steady state conditions, but originate from a source with up to 30% more partial melting than older ＞4 ka eruptions.This finding that deglaciation affects volcanism both on short and longer timescales may prompt a re-evaluation of eruptions in glaciated and previously glaciated terrains both in Antarctica and beyond.

Keywords glaciovolcanism, Mt Erebus, Antarctica, Earth Systems

1 Introduction

Continental scale geological changes are typically slow processes; mountain belts form and plates collide over tens to hundreds of millions of years. Ice sheets are one important exception to this rule, bodies of ice covering millions of square kilometres can grow and recede over mere 100s to 1000s of years. In glaciated zones this rapid change exhibits a strong control on both the surface and subsurface. Volcanism is sensitive to pressure variations,thus joint glacial and volcanic regions are particularly dynamic (Smellie and Edwards, 2016).

The relationship between glaciation and crustal loading is fairly well understood: when glaciation is initiated the subsurface is temporarily compressed as the ice grows faster than it is isostatically compensated. Conversely when deglaciation occurs the crust and mantle are temporarily decompressed as the crust rebounds and the denser mantle slowly flows back (Smellie and Edwards,2016). One often overlooked consequence of this is that volcanism can be locally inhibited by the growth of ice sheets and magnified as ice retreats. Wherever the mantle is near its melting point, changes in ice coverage affect mantle melt fraction and thus eruption rate (Jull and McKenzie,1996).

Increased eruption rates have been identified in many previously glaciated landscapes (Huybers and Langmuir,2009). For instance, volcanism in the rocky mountains was higher in the early Holocene corresponding to the removal of the Cordilleran and Laurentide ice sheets (Watt et al.,2013). In Iceland, volcanism rates in the thousand years immediately following deglaciation (10–9 ka) were almost 30 times background glacial and present day rates(Schmidt et al., 2013). As well as being apparent in the abundance and thickness of flows erupted at that period,increased melt production in Iceland has been confirmed via numerical modelling (Jull and McKenzie, 1996) and identified in the geochemistry of lavas (Hardarson and Fitton, 1991).

Mt Erebus is the youngest and most active volcano on the heavily glaciated Antarctic plate (LeMasurier and Thomson, 1990). Mt Erebus is a 3792 m high polygenetic stratovolcano, composed of a stratocone built upon a mafic shield base, similar to Mount Etna (Esser and McIntosh,2004). The eruptive centre itself is 1.3 Ma old, and has been erupting at an increased rate for the last 250000 years(Parmelee et al., 2015). To this day Mt Erebus hosts a convecting phonolite lava lake and a number of active ice fumaroles (Kelly et al., 2008; Wardell et al., 2004). Figure 1 shows the location and overall geological context of Mt Erebus- note the abundance of volcanism throughout West Antarctica.

Figure 1 Overall geological, volcanic and glaciological context of Mt Erebus. Note the abundance of volcanism on the Antarctic plate,especially in the deep subglacial basins of West Antarctica. Mt Erebus is located on one branch (Terror Rift) of the West Antarctic Rift System that bisects the continent. Volcano locations from Van Wyk de Vries et al. (2017); local Mt Erebus geology from Rilling et al.(2007).

Volcanism in Antarctica is by no means unusual,particularly in the western regions. A large continental rift, the West Antarctic Rift System, extends from the Ross Sea to the Antarctic Peninsula and is associated with extensive volcanism (Van Wyk de Vries et al., 2017).The region is still strongly data-limited, however around 40 exposed and nearly 100 subglacial volcanoes have been identified (Figure 1). Recent studies suggest that runaway deglaciation of sections of the West Antarctic ice sheet is possible, or even likely (DeConto and Pollard,2016), thus many of these volcanoes may soon experience reductions in pressure similar to around Mt Erebus 10–5 ka (Hall and Denton, 1999). Understanding the effects of deglaciation on Mt Erebus can thus play an important role in forecasting the future of volcanic activity in Antarctica, and may be useful for accurately measuring future ice-melt rates.

Mt Erebus has been active for much of the last 250 ka,so the underlying mantle necessarily oversteps its solidus to some degree. Figure 2 illustrates how the growth and retreat of ice sheets can affect the mantle geotherm and generate additional melt. Thickening and thinning of ice sheets affects the pressure in the underlying mantle and can raise (deglaciation) or depress (glaciation) the geotherm from its average position. In Figure 2a, a typically non-volcanic region becomes active during deglaciation periods. Note that this does nevertheless require a mantle close to its melting point to be effective;most non-volcanic regions will not melt even during deglaciation. In Figure 2b, a moderately volcanic region usually generates small amounts of melt, and experiences peaks in melt production during deglaciation. This same region may also become volcanically inactive during the glacial advances due to depression of the geotherm. In Figure 2c, a strongly volcanic region naturally generates large amounts of melt. These conditions make it particularly sensitive to changes in glacial loading-large peaks in volcanic activity would be expected during deglaciation. These peaks may however be rapidly overprinted as volcanism rates remain high in steady state conditions and moderate even during glacial advances.

Figure 2 Schematic illustration of how both glaciation and deglaciation can affect the generation of melt in the subsurface. See Huybers and Langmuir (2009) and Jull and McKenzie (1996) for more details.

Another potential source of lava is the emptying of already established shallow magma chambers due to rapid loss of overburden pressure. Recent advances in both the dating of Erebus’s ＜10 ka lava flows (Parmelee et al., 2015)and the timing of glacial retreat in the region (Anderson et al., 2014; Bentley et al., 2014; Ingólfsson, 2004) makesa study of these processes more reliable. In this study we aim to:

(i) Determine whether ice loss during the last interglacial period stimulated higher eruption rates at Mt Erebus.

(ii) If so, constrain what igneous, glacial or tectonic processes generated this increase in activity.

(iii) Establish whether present day Mt Erebus activity is steady state, or whether it is still adapting to past glacial changes.

First we build up a model to evaluate the magnitude of changes around Mt Erebus and what effect this could have on the rate of volcanism. We will then follow that up with an analysis of the geochemistry of Mt Erebus’s eruptive products. Finally, these results will be considered within the context of West Antarctica as a whole to assess their relevance to other glaciovolcanic systems.

2 Coupled deglaciation model

The timing of glacial retreat and maximum ice thickness in the Ross Sea has been the focus of many studies (Anderson et al.,2014; Bentley et al., 2014; Ingólfsson, 2004; Hall and Denton,1999). Bathymetric studies show glacial erosion marks, glacial lineations and moraines up to the edge of the continental shelf,confirming the extent of ice at the glacial maximum (Anderson et al., 2014; Bentley et al., 2014). The last glacial maximum grounding line is in places over 1000 km beyond its current position, which is a measure of the volume of ice lost. The inland reaches of the Antarctic ice sheet are not thought to have changed much in the last 30 ka, and greatest ice loss occurred in the Ross and Weddell seas (Bentley et al., 2014; Ingólfsson,2004; Hall and Denton, 1999).

The Ross Sea’s current bathymetry is made up of deep glacial troughs and wide plateaus due to differential glacial erosion and graben-horst style rifting (Paulsen et al., 2014).There is some lateral variation in depth, but the average is close to 500 m below sea level (Fretwell et al., 2013). Marine based ice sheets are often prone to rapid retreat, especially if strong grounding lines are breached. Recent studies in Greenland have shown that bases of glaciers exposed to warm sea waters can be rapidly eroded at the contact (Rignot et al.,2010). The large ice shelves present in Antarctica’s Ross and Weddell seas are remnants of former ice sheets covering the area that have been undercut by warm sea waters with only floating ice remaining. Ice shelves are composed of floating ice, so do not exert any extra load on the crust.

Sea based ice is thinned mainly by rapid grounding line retreat and basal melting by seawater, whereas land based ice in cold regions thins mostly through excess ice outflow(dynamic thinning). The two components are linked as land based ice retreat rate is related to marine grounding line position. Figure 3 shows the two main hypotheses for ice retreat in the Ross Sea, one with a gradual retreat perpendicular to the Victoria Land coast, and the other with rapid growth of a central embayment (Lee et al., 2017). The discrepancy between the models stems from an absence of data from the central part of the Ross Ice Shelf, however the retreat is much better constrained around Mt Erebus. Note that both models involve the grounding line retreating past Mt Erebus between 8 ka and 7 ka, and involve similar timings of inland ice retreat on Victoria Land.

Figure 3 Depiction of the two main theories for glacial retreat in the Ross Sea, modified from Lee et al. (2017) and Anderson et al. (2014).The differences in timing are relatively minor in the Mt Erebus area where retreat has been well constrained (e.g. Hall and Denton, 1999).

a) Marine deglaciation

Three main factors influence the retreat of marine grounded ice sheets: calving rate, basal melting and surface melting. In this region, surface melting is negligible (cold polar climate)and calving can only occur on the ice sheet boundary exposed to basal melting. By assuming calving ice was exposed to basal melting long enough to be floating (ice shelf conditions),calving is also discounted. This leaves only one term: basal melting (b) due to warm seawater. This melting only affects ice on the marine side of the grounding line, so the rate of grounding line retreat must be defined. A coordinate system is set up with x parallel to the Victoria Land coast, y perpendicular to this coast and z vertical. The ice sheet is considered uniform in the y direction, so this is a relation between x and t: x=θ(t), with θ(t) a function of time measuring the retreat of the grounding line. To measure θ(t), grounding line position data (Lee et al., 2017; Anderson et al., 2014) was used, the 10 ka position is taken as the starting point (x = 0). A polynomial function was fitted to this data (R2= 0.998)giving x=θ(t)=45t2–7.5t which accurately depicts the ice retreat rate until around 4 ka. Basal melting b is thus set as b=0 for x＞θ(t) and b=k with k a constant describing the melt rate at x＜θ(t). Data suggests the initial ice sheet thickness was 650 m, and the sea bed has an average depth of 500 m(Anderson et al., 2014; Fretwell et al., 2013). A value of 0.1 m·a-1was here used for k, this is an intermediate melt rate for ice interacting with sea water (Depoorter et al., 2013;Jenkins and Doake, 1991). This gives a thickness z: z=650 for x＞θ(t), and z=650–b(t–tg) for x＜θ(t), with t time and tgthe date at which the grounding line retreated past this point.

Since water has a higher density than ice only a small portion of the marine ice melted will actually result in decompression before the ice sheet reaches ice shelf status equivalent to open sea. The densities of water and ice are respectively 1000 kg·m-3and 900 kg·m-3. Using this we can calculate that the mass of the 500 m thick water layer Twpresent today (Anderson et al., 2014; Fretwell et al., 2013) is equivalent to a layer of ice with a thickness Tigiven by:

with ρwand ρithe respective densities of water and ice. The ice overload is 650–555=95 m, only 95 m of the total ice column affect pressure in the underlying mantle when removed(ignoring eustatic changes). After this point the ice sheet transitions into an ice shelf and subsequent melt has no effect.Figure 4 shows the evolution of the marine sections from an ice sheet to a retreating ice shelf from 9 ka to present day (the thickness is simplified to a uniform 650 m at 10 ka). The model was run with a temporal resolution of 500 years. Mt Erebus is a distance of 260 km from the 10 ka grounding line, at this distance the ice sheet–shelf transition occurs between 7 ka and 6 ka. This model also reproduces the present day thickness of the Ross Ice Shelf with an error of only 50 m (within ±20%).

The reduction in pressure generated by this glacial unloading can be calculated (assuming an ice density of 900 kg·m-3). Pressure P is given by P=ρgh with g the acceleration due to gravity and h the thickness of the layer.Here all 95 m of relevant marine ice are lost, so P=9.8×900×95=8.4×105Pa=0.84 MPa.

b) Grounded ice loss

Land based (grounded) ice reacts differently to marine ice as the grounding line will not be undercut by warm seawater.As a result basal and surface melting are both negligible terms in the overall mass balance. However, as the marine ice sheet retreats, the land ice is no longer in equilibrium and outflow will dominate inflow. This results in a steady thinning of the ice without significant melt. The rate at which thinning occurs in a given area should increase once ice is no longer grounded. A thinning term τ (Depoorter et al., 2013; Rignot et al., 2010) can be defined as: τ = 0.01 m·a-1for x＞θ(t) and τ = 0.02 m·a-1for x＜θ(t).

As the land and marine based ice sheets are defined in the same coordinate systems, this θ(t) term is the same as that defined in the previous section. At 10 ka the average thickness of ice on the coast inland of Mt Erebus was 550 m(Anderson et al., 2014) and the current ice thickness is an average of 100 m thick (Fretwell et al., 2013), which means 450 m of ice were lost. In contrast to marine ice, all of this ice column will contribute to depressurisation when removed.

Figure 5 presents the coastal ice sheet’s thickness against distance inland throughout time. A sharp inflexion in the thickness profiles shows the retreat of the grounding line and its thinning effect on land based ice. The reduction in pressure ΔP is given by ΔP=ρgΔh with Δh the total ice lost between 10 ka and a given time t. However, given that Mt Erebus is 60 km offshore (Fretwell et al., 2013) it will experience only a portion of this pressure drop. The flexural equation (Turcotte and Schubert, 2002) can be used to calculate what percentage of the maximum loading is experienced at Mt Erebus. This equation relates an increase or decrease in pressure W to distance γ from the load:

with P the horizontal pressure acting on the plate (positive if compressional, negative if extensional), g the acceleration due to gravity, ρmthe density of the mantle (3300 kg·m-3), ρwthe density of the load (here the load is ice so 900 kg·m-3), flexural rigidity D is a measure of how well these loads are distributed over distance (Capra and Dietrich, 2008; Kaufmann et al., 2005)and q()γ is the vertical loading causing the flexure.

Figure 5 Modelled retreat of coastal ice on the Victoria coast/Transantarctic Mountains margin. A minimum ice thickness was set at 100m to reproduce present day thicknesses in the region (Fretwell et al., 2013). Note a rapid increase in ice loss after the grounding line retreats south of this point between 8 and 7 ka.

We can assume that the horizontal pressure gradients are negligible (P=0) and that loading q(γ) is 0 when γ is not 0. This will model the flexure beyond the edge of a load(here beyond the edge of the land based ice sheet) in an area in tectonic equilibrium (over the timescales considered).This simplifies the flexural equation to:

Using the boundary conditionsat γ=0 and limγ→∞w = 0 this equation can be solved to give:

with w0the depression or pressure overload at γ=0 and α the flexural parameter (Turcotte and Schubert, 2002). Finally, to simplify the comparison with values from the glacial model,we define a normalised value β(γ) such as:

Figure 6 Decay in compression away from the edge of a central load ‘broken plate model’ given by the isostatic flexure equation(Turcotte and Schubert, 2002) with Mt Erebus labelled.

β(γ) is a dimensionless quantity that describes what percentage of the central loading w0is induced at distance γ from this load. Figure 6 shows a plot of β(γ) against γ,showing the fraction of the land based ice load pressurisation that is transmitted to Mt Erebus (60 km offshore from the coastline). The bulge and decay rate predicted in the model closely fit observations of flexure near oceanic islands and seamounts, and post-glacial isostatic adjustment uplift values from Scandinavia and North America (Watts, 2002).

The coastal pressure reduction values can be multiplied by the flexural parameter for 60 km to give the pressure reduction at Mt Erebus. β(60)=69.3%, so the pressure reduction at Mt Erebus is 69.3% of that at the coast.Mt Erebus is assumed to be in contact with the ocean, so flexure is not considered for the marine based ice retreat.

Using this result, the combined pressure loss due to both the grounded and marine ice can be calculated. A total pressure drop of 3.66 MPa occurred between 9 and 3.5 ka, with peak decompression at 6 ka. Figure 7 shows a plot of this pressure drop through time along with the ages of Mt Erebus surface lavas calculated by Parmelee et al.(2015). The timing of the pressure reduction coincides well with the peak in lava age (also 9–3.5 ka), although eruptions are more discontinuous. There are two main ways in which depressurisation could affect eruption rate:an increase in the degree of partial melting at depth and the rapid emptying of pre-existing shallow magma chambers.

Figure 7 Comparison between the postglacial decompression pulse predicted from this deglaciation model and the age distribution of summit lavas from Parmelee et al. (2015). There is a close association between the two that suggests they may be linked.

c) Shallow magma chamber release

Whilst deglaciation does affect the pressure state in the mantle, the relative change in overburden pressurethe pressure drop divided by the total post-deglaciation pressure) will be low at depths. At shallow levels the crust’s temperature is far too low for any partial melt to occur,howevercan be very high (and can even exceed [1 at depths of less than 150 m for the pressure] drop calculated above). As a result of this, the local environment of any previously stable magma chambers, sills or dykes will drastically change. Similarly to how major volcanic edifice avalanches can induce eruptions (Siebert, 1984), rapid deglaciation can relieve pressure on previously stagnant bodies of magma and permit them to escape to the surface.Deglaciation is thought to involve not only reduced vertical stresses but also lesser horizontal stresses and increased regional strain rates (Watt et al., 2013). This joint action serves to promote opening of cracks, dyke formation and fault motion, all of which can induce eruptions. Figure 8 describes schematically the effect deglaciation can have on shallow magma chambers.

The exact additional volumes of lava that this process would provide are uncertain as they depend on the volume of magma stored in sensitive shallow magma chambers during deglaciation. Present day volcanological studies of Mt Erebus suggest that it is underlain by a fairly extensive magma plumbing system (e.g. Oppenheimer et al. 2011;Esser et al. 2004) that could have been tapped during deglaciation. As the stable fractionation of magma in these shallow chambers was interrupted by deglaciation, magmas erupted as a result would be more primitive (less fractionated) than non-affected lavas. In particular compatible element concentrations (Co, Ni, V, possibly Eu,Ca, Mg, etc; all largely insensitive to differences in degree of partial melting) would be higher in affected magmas as the relevant mineral phases would not have fully crystallised out of the melt. Incompatible elements (and even more so incompatible element ratios) would be much less affected by this shorter fractionation.

d) Mantle decompression melting

Figure 8 Schematic bloc diagram of Mt Erebus and the underlying crust-mantle, including the shallow magma plumbing system. This figure shows how the removal of superficial ice (and associated drop in overburden pressure) can cause rapid drainage of shallow magma chambers and increased eruption rates.

We used the MELTS software to quantify the effect of this decrease in pressure on mantle melting under Mt Erebus(pMELTS version 5.6.1: Ghiorsio and Sack (1995), Asimow and Ghiorsio (1998)). Defining a series of pressure and temperature steps allows MELTS to calculate a change in composition over this range. In this case we do not need to know the absolute melt fraction at a specific depth, only the sensitivity to changes in pressure:, with ΔF the change in melt fraction and ΔP the change in pressure. This ratio describes how much an increase or decrease in pressure will affect the amount of melting. Questions remain about the exact composition and temperature of the lithosphere below Mt Erebus, so appropriate lithospheric averages were used (Xia and Hao, 2013; Green et al., 2010;McKenzie and O’Nions, 1991). Bulk composition and temperature are considered constant at a given depth. For each depth,was calculated for three temperatures(within ±10% of the average geotherm) and averaged.

This function goes to 0 at 35 km depth (1 GPa), meaning no melt produced in the upper lithosphere, in line with observations(Esser et al., 2004). At 150 km depth (4.5 GPa)reaches around 0.012% per MPa, consistent with deep melting at a hotspot/continental rift. Huybers and Langmuir (2009) also suggest a comparablevalue of 0.01% per MPa.

The pressure reductions can be multiplied by the volume of a melting cone beneath Mt Erebus, weighted according to the depth vsrelation, to obtain the total melt production (as in Jull and McKenzie, 1996).Accordingly, the 3.66 MPa pressure drop linked to deglaciation generates a total of 280 km3of magma in the conical melt region. The largest part of this melt occurs at great depths (＞100 km), so only a small fraction may reach the surface. Studies have shown that in Iceland 5%–10% of the total melt may reach the surface, erupting on average 1 ka later (Slater et al., 1998; Jull and McKenzie, 1996). Melt at Mt Erebus is deeper and melt fraction is over an order of magnitude lower than in Iceland, so we could hypothesize that 0.5%–1% of the magma reaches the surface on average 10 ka after deglaciation. Figure 9 shows the processes involved in this deglaciation-partial melt increase cycle.

Overall this model predicts that up to 280 km3of additional melt may form at depth, with 1–3 km3of this being erupted. The exact volume erupted and timing of eruption are poorly constrained, however it would likely reach the surface spread out several thousand years after deglaciation. This process would be expected to result in a distinctive trace element pattern—with the higher than usual melt fraction resulting in different incompatible element ratios. Ratios of highly incompatible to moderately incompatible elements(Ce/Y, Nb/Zr, La/Ho, etc; insensitive to fractional crystallisation) would likely be lower during the eruption of these lavas (Hardarson and Fitton, 1991). Compatible element concentrations would be largely unaffected by these processes.

Figure 9 Schematic bloc diagram of Mt Erebus and the underlying crust-mantle. This figure shows how the removal of superficial ice can generate increased degrees of partial melting at depth and induce higher eruption rates.

3 Geochemical analysis

Basic conceptual modelling suggests that a link between ice coverage and volcanism should exist; postglacial depressurization should indeed induce eruptions. However,while the model is useful for understanding the processes occurring at Mt Erebus, it cannot alone be used as a guide for understanding what is occurring in reality. As the previous section explores, there are two main mechanisms by which glacial retreat could generate eruptions: increased deep partial melting and draining of shallow magma chambers (Jull and McKenzie, 1996). In each case the erupted magma has a very different origin which allows us to make specific predictions about its chemistry. As the geochemistry of both modern lavas (Kelly et al., 2008), lava flows from the late glacial period (8.5–4 ka, Parmelee et al., 2015) and older lavas from the glacial period (60–20 ka; Iverson et al., 2014) has been collected we can test these predictions. In this section this geochemical data is analysed in order to better understand the postglacial evolution of Mt Erebus.

Geochemical (major, minor and trace element concentrations) data for Mt Erebus lavas is available from several studies of the summit plateau lavas and volcanic tephra (Parmelee et al., 2015; Kelly et al., 2008), recent eruptive episodes are well sampled. Older flows are mostly buried under recent flows and ice so are less well studied,however englacial tephras from older eruptions have also been analysed (Iverson et al., 2014; samples from EIT008 and 034). For modern (＜ 4 ka) and late glacial lavas(8.5–4 ka), analyses of over 10 individual samples are available (Parmelee et al., 2015). Modern values come mostly from analyses of glass and crystals from tephra scatter which has been the main mode of eruption over the last 4 ka (Kelly et al., 2008).

The mean values and standard deviations of each element concentration were calculated for each time period: modern (＜4 ka), late glacial (8.5–4 ka) and glacial(60–20 ka). This data is presented in Table 1. Incompatible elements, such as K, Zr, Ce, Nb, Rb and Y (McKenzie and O’Nions, 1991) are very sensitive to changes in degree of melting, but less so to changes in fractional crystallisation(particularly if the ratios of these elements are taken,Hardarson and Fitton, 1991). Compatible elements such as Ni, V (and in most basic magmas Ca, Mg, Eu, etc.) are on the contrary largely insensitive to changes in melt fraction,but strongly affected by differences in fractionation(Ghiorsio and Sack, 1995).

Figures 10 and 11 present two plots of compatible elements versus SiO2. Silica concentration is only lightly affected by changes in fractionation and almost unaffected by differences in melt fraction. Figure 10 shows that concentrations of CaO and Eu, both highly compatible in common basic mineral plagioclase (and to a lesser extent clinopyroxene) were considerably higher in late glacial times than both modern and glacial equivalents. Figure 11 shows a similar trend for MgO (particularly compatible in forsteritic olivine and pyroxenes) and trace element V(small, low charge cation compatible in most mineral phases). Overall these differences reveal that late glacial lavas fractionated less mafic minerals than their modern and glacial counterparts, most likely because they spent less time in the magma chambers. This is consistent with the predictions of the shallow magma chamber release model discussed in the previous section.

Table 1 Geochemical data from three different time periods: glacial (60–20 ka), late glacial (8.5–4 ka) and present day(0 ka). In each case the mean and standard deviation are given for major and minor elements shown, along with a selection of relevant trace elements. Data synthesized from Parmelee et al. (2015), Iverson et al. (2014) and Kelly et al. (2008).

Figure 10 Plot of two compatible elements, Calcium (shown in green) and Europium (shown in blue) against Silica in glacial (60–20 ka),late glacial (8.5–4 ka) and present day lavas from Mt Erebus. Both elements are highly compatible in plagioclase and moderately compatible in clinopyroxenes two common early minerals to crystallise. Data from Parmelee et al. (2015), Iverson et al. (2014) and Kelly et al. (2008).

Figure 11 Plot of two different compatible elements, Magnesium (shown in blue) and Vanadium (shown in red) against Silica in glacial(60–20 ka) late glacial (8.5–4 ka) and present day lavas from Mt Erebus. Mg is highly compatible in pyroxenes and olivine, and Vanadium is compatible in most minerals. Data from Parmelee et al. (2015), Iverson et al. (2014) and Kelly et al. (2008).

Figure 12 shows a plot of two incompatible element ratios: Ce/Y and Zr/Nb. These ratios are unaffected by differences in fractionation and the elements involved are not fluid mobile. Cerium and Nb are both highly incompatible elements, compared to Y and Zr that are only moderately incompatible (McKenzie and O’Nions, 1991).This means that with increasing partial melt fraction Ce/Y would decrease and Zr/Nb would increase along a predictable trend line (depending only on the starting composition of the mantle). A non-modal fractional partial melting curve for a primitive garnet lherzolite mantle (a likely source composition for the mantle beneath Mt Erebus and the closest fit; Hardarson and Fitton, 1991) was plotted alongside the element ratios. Modern lavas form a distinctively separate cluster to late glacial and glacial lavas consistent with 0.2%–0.3% more deep partial melt(Hardarson and Fitton, 1991). This is in line with the model prediction that deglaciation would cause an increase in the degree of partial melting at depth, that would then make its way to the surface over several thousand years. This means that present day lavas formed from a mantle source melting up to 30% more than background rates (Figure 12).

Overall the geochemical analyses allow us to distinguish two distinct interactions that have occurred: a release of magma from shallow magma chambers shortly after deglaciation and an increase in the degree of deep partial melting that only makes it to the surface after a lag of several thousand years. These results allow us to build up an integrated model for the evolution of Mt Erebus since the last glacial maximum 20 ka. This model is presented in the following section, and evaluated in the context of other relevant studies.

4 Discussions

Figure 12 Plot of two incompatible element ratios: Ce/Y and Zr/Nb. These element ratios are well suited as they are not prone to fluid remobilisation and are unaffected by changes in degree of fractional crystallisation. A non-modal fractional melting curve for the most likely mantle composition (primitive garnet lherzolite) is also plotted on this diagram revealing that modern lavas derive from a mantle melting 0.2%–0.3% more than late glacial and glacial lavas(red symbols show average for each age bracket). Glacial lavas exhibit higher spread, likely due to the higher age variation and changes in the magmatic system. Data from Parmelee et al. (2015);Kelly et al. (2008) and Hardarson and Fitton (1991).

In Iceland, the best studied glaciovolcanic province, both deglaciation and the volcanic response to it occurred very rapidly (from fully glaciated to present conditions in as little as 1000 years). However Iceland’s unusual climatic and geodynamic environment make it particularly sensitive and quick to respond to changes. Comparisons between volcanoes in different environments are not always possible- while both Iceland and Mt Erebus experienced a postglacial burst in volcanism, this was related to different processes. In Iceland this peak (up to 30 times the background rate; Slater et al., 1998; Jull and McKenzie,1996) was mostly due to rapid flushing of large volumes of decompression melt, whereas at Mt Erebus it is instead due to evacuation of preexisting shallow magma chambers.

Figure 13 Generalised geological evolution sections for the Mt Erebus region from the last glacial maximum (roughly 20 ka) through to a theoretical future stable state. One take home point is that deglaciation affected Mt Erebus’ eruption rate through two entirely distinct processes:one with a rapid response and one with a much slower response still ongoing to the present day. See Figure 14 for further discussion.

Figure 14 Schematic interpretation of the volcanic response to deglaciation over time in three different environments. In high melt regions (e.g. Iceland) the increase in mantle melting dominates and the two processes are quasi indistinguishable (Jull and McKenzie, 1996).In lower melt regions (e.g. Mt Erebus or elsewhere in the Antarctic rift volcanics) the mantle melting response is much slower and drawn out, and in extreme cases may be have a negligible effect on eruption rates. This effect is ongoing at Mt Erebus, with erupted magmas originating from up to 30% more mantle melting than background rates.

Figure 13 shows a generalised history of Mt Erebus from the last glacial maximum (20 ka) to 10 ka in the future.The eruptive mode (and lava chemistry) would in theory be similar whenever the system is in steady state: at the last glacial maximum (stable glacial) and 10 ka in the future(stable non-glacial). The 10 ka future step of course assumes the surface conditions do not further change during this time (no further ice melt or glacial readvances to disrupt the system). The late glacial period (8.5–4 ka) is dominated by eruptions sourced from shallow magma chambers rendered unstable via the loss of glacial overburden pressure. Importantly we note that the present day conditions are not steady state, Mt Erebus is still erupting magmas derived from unusually high melt fractions due to deep mantle decompression.

Deglaciation has previously been hypothesized to have a longer timescale effect on volcanoes than that observed in Iceland (e.g. Huybers and Langmuir, 2009;Jull and McKenzie, 1996), however traces of the last deglaciation have not previously been identified in present day lavas. There are two main reasons why Mt Erebus is an ideal location to observe this: peak deglaciation occurred later than elsewhere (8.5–4 ka,Anderson et al. (2014); as opposed to 15–10 ka in Iceland, Jull and McKenzie (1996); and North America,Booth et al. (2003)) and deep, low degree melting means that magmas take a longer time to reach the surface. The lack of weathering and vegetation coverage also makes analysis of past lava flows easier than in lower latitude regions (Esser and McIntosh, 2004). Figure 14 shows the changes in eruption rate linked to release of shallow magma chambers and changes in melt fraction at depth.In Iceland the two peaks are high amplitude, short period and occur contemporaneously. At Mt Erebus (and likely in other low melt regions, for instance throughout the West Antarctic volcanic province) a first peak related to shallow magma release occurs shortly after deglaciation.The second peak, from an increase in melt fraction at depth is however drawn out and forms a low amplitude-long period peak streching several thousand years following deglaciation. This longer period is due to a slower mantle response, but also linked to longer magma source to surface timescales.

While the volcanology of Mt Erebus has been studied in quite some detail (Parmlee et al., 2015;Oppenheimer et al., 2011; Kelly et al., 2008), little attention has been paid to the effects of glaciation.Indeed prior studies have suggested that the chemistry of Mt Erebus’s magmatic system has been unchanged over the last 40ka (e.g. Iverson et al., 2014; Kelly et al., 2008).However these conclusions are based almost entirely on40Ar/39Ar dating of flows and tephra with very high uncertainties, for instance Kelly et al. (2008) dated Mt Erebus’s surface flows as between 17±8 ka and 0±4 ka.These same surface flows were constrained by Parmlee et al. (2016) to between 8.50±0.19 ka and 4.52±0.08 ka using exposure age dating. Thus low resolution40Ar/39Ar dates not only failed to recognize an eruption peak during deglaciation, but also smeared both glacial and postglacial chemical analysis data over the last 20 ka creating an artificial equilibrium. These dates are useful for understanding the longer timescale volcanic evolution,however cannot reliably be used to interpret the ＜10 ka evolution of Mt Erebus due to uncertainties inherent in the40Ar/39Ar dating technique.

Given that the region has been glaciated intermittently since as long as 34 Ma (Convey et al., 2018), long before Mt Erebus began erupting, ice played an important role in the formation of this large stratovolcano. For instance if ice re-advances before all of the deglaciation related partial melt pulse is erupted, the increased overburden pressure would lock the magma in subsurface bodies. While this magma would not be apparent on the surface, it would still contribute to the overall structure’s size. Indeed, Mt Erebus’s total volume cannot be explained by its current eruption rate as inferred from exposed lava flows (Esser and McIntosh, 2004). Large subsurface magma emplacement is a way to reconcile the structure’s volume with its eruption rate.

Importantly, the link between deglaciation and volcanism at Mt Erebus may improve our understanding of how the rest of the volcanic province, spanning the length of the West Antarctic Rift System (Van Wyk de Vries et al., 2017), would react to future changes. Rapid grounding line retreat and ice thinning is increasingly occuring, particularly in the Amundsen Sea sector(DeConto and Pollard, 2016). If the volcanoes in that area react similarly to Mt Erebus, a peak in eruption rate would be expected and some currently inactive volcanoes may be reawoken by the pressure drop.

It is difficult to determine whether the increased eruption rate from Mt Erebus had any effect on the rate of deglaciation (Watt et al., 2013). Lava flows only locally melt ice caps and it is unlikely that they had much effect on regional ice coverage beyond removing Mt Erebus’s summit glaciers. In the larger picture however,erupting several cubic kilometres of lava would also emit considerable tephra and greenhouse gasses. Greenhouse gasses warm the climate on a local and global scale by trapping outgoing radiation, while tephra can dirty the ice and decrease its albedo, also warming the climate.This warming can in turn accelerate deglaciation (Watt et al., 2013; Huybers and Langmuir, 2009). This feedback cycle may play a role in promoting runaway deglaciation,and is of particular relevance in today’s rapidly changing world. These results suggest the earth system is deeply coupled and interactions between volcanoes, glaciation and the earth’s interior may be greater than previously expected.

5 Conclusion

In this study the links between deglaciation and volcanism were investigated on Antarctica’s most active volcano, Mt Erebus. We started by building up a deglaciation model to evaluate the effect that deglaciation may have had on the underlying crust and mantle. This revealed that 3.66 MPa of depressurisation occurred between 9 ka and 3.5 ka,coinciding with a burst in magmatism at Mt Erebus. Next,we analysed the geochemistry of a suite of glacial to present day lavas to evaluate model predictions. This led to three main results:

(1) Deglaciation has had a large impact on Mt Erebus’activity over the last 10000 years, and possibly also its long term evolution.

(2) Mt Erebus’ peak in eruption rate between 8.5 and 4 ka is linked to rapid emptying of shallow magma chambers,destabilised by the loss of glacial overburden pressure.

(3) Present day conditions are not steady state, but originate from up to 30% more mantle partial melting than background rates, also linked to the post-glacial pressure drop. Thus Mt Erebus’s eruptive system has not been constant over the last 10 ka, and continues to change today.

With Antarctica currently experiencing the largest and most rapid changes in ice volume since the last deglaciation,other volcanoes throughout the West Antarctic Rift System may also experience similar peaks in activity. While not expected to be a major driver of ice retreat, this question certainly warrants further research in this highly sensitive region.

AcknowledgmentsEarly versions of this manuscript greatly benefitted from discussions with a number of colleagues from the University of Edinburgh and beyond. Valuable and constructive suggestions by anonymous reviewers helped improve and refine the final version of this paper, and provide independent verification of the analytical techniques used.