Pressure influences both magma production and the failure of magma chambers. Changes in pressure interact with the local tectonic settings and can affect magmatic activity. Present-day reduction in ice load on subglacial volcanoes due to global warming is modifying pressure conditions in magmatic systems. The large pulse in volcanic production at the end of the last glaciation in Iceland suggests a link between unloading and volcanism, and models of that process can help to evaluate future scenarios. A viscoelastic model of glacio-isostatic adjustment that considers melt generation demonstrates how surface unloading may lead to a pulse in magmatic activity. Iceland’s ice caps have been thinning since 1890 and glacial rebound at rates exceeding 20 mm yr−1 is ongoing. Modelling predicts a significant amount of ‘additional’ magma generation under Iceland due to ice retreat. The unloading also influences stress conditions in shallow magma chambers, modifying their failure conditions in a manner that depends critically on ice retreat, the shape and depth of magma chambers as well as the compressibility of the magma. An annual cycle of land elevation in Iceland, due to seasonal variation of ice mass, indicates an annual modulation of failure conditions in subglacial magma chambers.
Volcanic systems are recharged by magma flow from depth, which causes a build-up of pressure until failure conditions are reached, and an intrusion or an eruption occurs. Events that cause stress changes, such as earthquakes and eruptions at other nearby volcanoes, can thus influence the complex inner working of a particular volcano. Another source of pressure changes is variable load on the surface of the Earth. Retreating ice caps due to the warmer climate influence underlying volcanoes. One metre reduction of ice, with a density of 910 kg m−3, corresponds to a normal stress change on the surface of the Earth of (910 kg m−3)×(9.8 m s−2)×(1 m)=9 kPa, an order of magnitude smaller than the stress variation of 1 bar (100 kPa) often referred to in earthquake triggering studies. Here, we model stress induced by variation in surface loads and evaluate how the resulting pressure conditions can modulate magmatic activity. We demonstrate the effects of unloading on volcanic systems, with examples from Iceland (figure 1), which is located at the divergent plate boundary between the North American and Eurasian plates and above an inferred mantle plume (e.g. Sigmundsson 2006). Reduction of pressure due to tectonic extension enables mantle rocks to melt, creating a zone of magma upwelling underneath Iceland. A striking case of a major influence of ice retreat on volcanic activity is the large pulse of volcanic activity in Iceland at the end of the last glaciation. Volcanic activity appears to have been more than 10 times more frequent during the early postglacial history of Iceland than presently (e.g. Sinton et al. 2005), associated with the disappearance of the more than 300 km wide ice sheet that covered Iceland. Pressure decrease at depth leading to mantle melting or flexure of the crust creating ‘easier’ pathways for magma to reach the surface have been suggested as the cause. Modelling of the pressure release effect shows how a large volume of magma was generated (Jull & McKenzie 1996).
The volcanic zone in Iceland experiences various load changes at present. Large ice caps presently cover parts of the plate boundary in Iceland and these have been retreating since 1890 due to the warmer climate. The thinning is responsible for presently ongoing glacio-isostatic adjustment. Widespread uplift occurs, with maximum rates exceeding 20 mm yr−1 (figure 1). The glacio-isostatic adjustment has been measured and modelled, for example, by Pagli et al. (2007) and Árnadóttir et al. (2009). GPS geodetic measurements show furthermore an annual cycle in land elevation that can be attributed to seasonal variation of ice mass (Grapenthin et al. 2006), indicating an annual modulation of failure conditions in magma chambers of subglacial volcanoes. More sudden local changes also occur. The Grímsvötn subglacial volcano at Vatnajökull has a caldera-filled subglacial lake that drains in jökulhlaups, sudden outburst floods. Some of the eruptions of this volcano have taken place after the onset of a jökulhlaup, such as an eruption in 2004 (Sigmundsson & Gudmundsson 2004). Links between the associated load reduction and eruption triggering were initially suggested by Thorarinsson (1953).
Redistribution of surface loads on volcanoes, such as partial destruction of edifices or flank destabilization events, have been linked through models to the onset of eruptions and changes in eruption rate (e.g. Pinel & Jaupart 2005; Manconi et al. 2009). Ice load variations are expected to have a similar effect, with influence depending critically on the size of the load. A retreating ice cap with a radius of only a few kilometres will influence only the shallowest parts of a magmatic system, including a shallow magma chamber. A retreating ice cap with a radius of tens of kilometres or more may, on the other hand, influence conditions in the deepest part of magmatic systems, the melt generation zone within the mantle (figure 2a). Here we consider both of these processes.
2. Pressure variation under a disc load
Magma can be generated at depth in the mantle by decompressional melting when material moves upwards towards lower pressure. Reduction in surface load can similarly reduce pressure at depth in the mantle, causing melting. In addition to details of the surface load, the induced pressure effect at depth in the Earth will depend on rheology. The initial response of the Earth to load variation is elastic, followed by viscoelastic adjustment of the layers below the uppermost brittle elastic crust. The initial elastic response is here compared with that of a viscoelastic Earth model, for loads with axisymmetric geometry. We consider the influence of load variations on variation in stress from a reference stress state existing prior to load variations. Pressure, P, is a parameter of fluids, but within an elastic material it can be defined as the average of three orthogonal normal stress components. In cylindrical geometry it is given as 2.1 where σzz is the vertical normal stress, σrr the radial normal stress and σθθ the tangential normal stress. Ice load variations induce a change in the normal stress at the surface of the Earth equal to ρgh, where ρ is the load density, h the thickness change of the load and g the gravitational acceleration.
(a) Elastic half-space
When modelled as surface load variations on an elastic half-space, the resulting displacement and stress changes at depth in the Earth can be inferred from the solution to a surface point load, known as a solution to the Boussinesq problem. Integration of the effect of a point load over a circular area gives the stress variation under a disc load. At the axis, the vertical normal stress as a function of depth below the surface, z, equals (e.g. Davis & Selvadurai 2001) 2.2 where p0=ρgh is the surface load, and R is the radius of the load. The other two horizontal stress components are equal (e.g. Pinel 2000, Pinel & Jaupart 2004): 2.3 where ν is Poisson’s ratio. From equations (2.1) to (2.3) it is found that the pressure under the centre of a disc load overlying an elastic half-space is given by 2.4 Stress and pressure versus depth are shown in figure 2b.
(b) Elastic layer overlying a viscoelastic half-space
When considering long-term pressure variations, the viscoelastic properties of the mantle and lower crust have to be taken into account. Pagli et al. (2007) and Pagli & Sigmundsson (2008) used an axisymmetric finite element model to study stress and pressure changes below Vatnajökull. The Earth was modelled as an incompressible elastic plate overlying an isotropic, incompressible Maxwell viscoelastic half-space. Vatnajökull was modelled as a circular ice cap. The ice retreat history assumes isostatic equilibrium in 1890 and gradual thinning of the ice cap since 1890. Pagli et al. (2007) showed that a viscosity of 8×1018 Pa s, assuming an elastic plate thickness of 10 km, gives a best fit to their GPS-derived displacement field at the southeast edge of Vatnajökull. Using that model, Pagli & Sigmundsson (2008) estimated the largest rate of pressure decrease of about 1.7 kPa yr−1 in the mantle beneath the ice cap. Árnadóttir et al. (2009) modelled the vertical velocities from a different nationwide GPS data set with respect to glacial isostatic adjustment in a three-dimensional finite element model (after Lund & Näslund 2009), incorporating the four largest ice caps in Iceland. They found that their data is best fitted by an Earth model with a 40 km thick compressible elastic layer overlying a Maxwell viscoelastic half-space of 1×1019 Pa s, or equivalently well by a three-layer model with a 10 km elastic layer and a 30 km, 1×1020 Pa s, viscoelastic layer overlying the same half-space.
Here, we illustrate the pressure decrease below a melting ice cap using an axisymmetric finite element model (figure 3) with a 30 km compressible elastic plate (density 2800 kg m−3, Young’s modulus 40 GPa) overlying a 1×1019 Pa s compressible Maxwell viscoelastic half-space (density 3200 kg m−3, Young’s modulus 130 GPa). The effect of viscosity is shown by considering both the first 10 years (figure 3a,c) and the last 10 years (figure 3b,d) of the simulation. We see that the ice model has a profound influence on the results and that ice caps with a higher deglaciation rate at the edge (the normal situation) produce smaller depressurization magnitudes than would a spatially uniform melting rate. Also, the pressure changes are larger beneath the edge of the ice cap than at the centre in the non-uniform model. However, at depth below the ice cap the pressure changes are similar for the two ice models and more dependent on the size of the ice cap and the total volume loss. The early rate of pressure change is different from the late rate, in both magnitude and spatial distribution. The viscoelastic half-space relaxes the stresses in the lower part of the elastic layer.
(c) Comparison of models
The estimates of pressure induced under a retreating ice cap, inferred from different models, are shown in figure 4. Below 30 km depth and throughout the mantle, the estimates of the pressure changes are similar, and all decay with increasing depth. Even the elastic model predicts values that are comparable with values predicted from viscoelastic models, all mostly depending on the total load. The lowest estimate of the pressure changes is from the model by Pagli & Sigmundsson (2008), which provides therefore a minimum estimate for the pressure release effect. The variation between models depends on the different rheological parameters and assumptions used, as well as the finite element model set up. Within the elastic layer, large variation between models is observed. The pressure within the elastic layer depends strongly on the elastic parameters and the viscosity below it. As the loading or unloading phase of glacio-isostatic adjustment continues, the stress distribution in the elastic layer resembles more and more that of a flexing beam.
3. Influence on deep magma generation
The melt generated by glacial unloading can be obtained from the pressure changes, provided that a relation between pressure change, P, and melt fraction by weight, X, is known. Following Jull & McKenzie (1996), we assume isentropic decompression melting and that the substantive derivative of melt fraction and pressure are related in the following way: 3.1 where is the velocity of the solid. The partial derivative of melt fraction with respect to pressure for constant entropy, (∂X/∂P)S, can be integrated numerically if the relation between X, P and temperature, T, is known, e.g. from laboratory measurements (McKenzie 1984). Pagli & Sigmundsson (2008) integrated eqns (D7) and (D8) of McKenzie (1984), using a fourth-order Runge–Kutta scheme, starting on the solidus (X=0) at a temperature of 1500°C. The partial derivative of melt fraction with respect to pressure is about constant, taking the value of 8.7×10−11 (Pa−1). Multiplication of the rate of pressure change by this factor gives the melt fraction produced per year (figure 5). This is true only where a pre-existing magma upwelling region already exists—like in Iceland or other spreading centres. To infer the total integrated melt volume created by ice decompression, the melting region needs to be specified. We follow Jull & McKenzie (1996) and assume that melting is confined to a triangular melting region centred at the ridge axis. We also assume that melting occurs exclusively in the mantle, between a depth of 25 and 112 km, and consider only the area underneath the Vatnajökull ice cap, where ice thinning and upwelling directly combine. With this model, we find that since 1890 the average melt production rate due to glacio-isostatic decompression is about 0.014 km3 yr−1. Considering that other models predict somewhat larger pressure changes, we take this as a minimum estimate for the melt generated under the Vatnajökull ice cap due to its present-day thinning.
The melt estimate can also be inferred from a three-dimensional model. Lund et al. (2009) used the three-dimensional model of Árnadóttir et al. (2009) to calculate the pressure decrease below Iceland. The maximum pressure change inferred is somewhat larger than the estimate by Pagli & Sigmundsson (2008). As above, they integrated the estimated melting rates, but now in a triangular region below the entire length of the rift axis in Iceland. The resulting melt production rate depends on how the volume of the melting region is estimated and plausible models yield rates of 0.04–0.06 km3 yr−1. This includes not only the Vatnajökull region but also volcanic centres north and south along the rift.
4. Influence on shallow magma chambers
(a) Pressure changes within and around magma chambers
In previous sections, pressure changes induced by surface load variations were estimated considering a laterally homogeneous crust. However, surface load variations above volcanoes usually act on a heterogeneous crust containing localized storage zones. We present axisymmetric numerical models in order to quantify the effect of surface load variations on a shallow magma chamber, considered as a cavity of an idealized shape (sphere or ellipsoid) embedded within an elastic homogeneous crust, filled with an inviscid fluid. The location of a rupture at a magma chamber boundary is sensitive to its shape. Largest hoop stresses are usually generated where the radius of curvature of the boundary is smallest, except for external stress field effects. Here, for modelling purposes, we only considered idealized shapes (spheres or ellipsoids). The presence of protuberances on magma chamber boundaries would modify our results somewhat, with rupture location controlled by the protuberance rather than by the overall reservoir shape. Protuberances on magma bodies can, however, be expected to be a result of previous dyke or sill injections from it. Their distribution, reflecting the state of stress, influences the general shape of a magma chamber. An ellipsoidal-shaped magma chamber model may thus approximate to some extent the influence of repeated sill or dyke intrusions on the shape of magma chambers.
Our model considers only elastic rheology and is most appropriate for short loading/unloading times. The magma has the same density as the surrounding crust and the reference state is lithostatic. Details of the numerical method are given by Albino et al. (in press). We estimate the magma pressure change (ΔPc) occurring within the magma chamber. Magma pressure variation, shown on the grey area in figure 6a, is the result of both the deforming crust and the magma equation of state (Pinel & Jaupart 2005; Albino et al. in press). For an incompressible magma, the chamber volume remains constant and the pressure applied at the surface is compensated by a magma pressure increase. For compressible magma, part of the stress increase is compensated by a decrease of the chamber volume (figure 6b). Magma pressure change is therefore highest for incompressible magma, but decreases towards zero as magma compressibility increases (magma bulk modulus K approaches zero). Stress variations at magma chambers are due to both magma pressure changes as well as cavity deformation.
(b) Dyke initiation at the chamber wall and eruption likelihood
Surface load variation induces both a magma pressure change, ΔPc, and a modification of the threshold magma pressure required for dyke initiation, ΔPr. Following Pinel & Jaupart (2005) we consider that failure in tension of the chamber wall occurs when the minimum compressive deviatoric stress compensates the tensile strength (Ts) of the host rocks. Applying this rupture criterion in three dimensions, one can define a threshold magma pressure Pr required for dyke initiation (Albino et al. in press). An unloading event always induces a magma pressure decrease but can either increase or reduce the threshold pressure for dyke initiation. The difference in the two terms, ΔPr−ΔPc, provides the relative evolution of a volcanic system between the initial and the final state and characterizes the effect of the surface perturbation on the ability of the system to erupt (figure 7). A positive value signifies that the magma chamber state moves away from rupture conditions and an eruption is inhibited. Conversely, if the sign is negative, the magma chamber approaches rupture conditions and an eruption is favoured. Figure 8 shows the consequences of removing a large-scale disc load (50 km radius) above a 10 km3 magma chamber located at 3 km depth. The chamber has an ellipsoidal shape characterized by its ellipticity, e, the ratio between half-width, a, and half-height, b. Pressure variations are normalized by the amplitude of the load removed. The induced magma pressure reduction increases with decreasing magma compressibility. This effect is of a similar magnitude to the load, and larger for laterally elongated chambers. For shallow reservoirs the effect of the threshold for dyke initiation depends on the amplitude of the load and strongly on the reservoir shape. The threshold increases for prolate reservoirs, decreases for spherical and oblate reservoirs, the effect being maximum (two times the load removed) for spherical ones. The pressure difference (ΔPr−ΔPc) shows that an unloading event inhibits rupture initiation for a prolate chamber (positive values) and promotes rupture for a spherical and an oblate chamber (negative values). The enhancement of eruption likelihood is maximum for a spherical chamber. Increasing compressibility of magma inside the chamber increases the ‘enhancement effect’ for spherical and oblate chambers, and reduces the ‘inhibition effect’ for prolate chambers.
Climate warming results in focused peripheral retreat of ice caps rather than uniform thinning. As shown by Albino et al. (in press), for a peripheral unloading, both pressure changes are smaller in amplitude than in the case of a central unloading event, with values always smaller than the surface load removed. The maximum magma pressure change, ΔPc, reaches only 40 per cent of the load removed, and the amplitude of the difference, ΔPr−ΔPc, is always less than half of the load removed. Magma pressure changes strongly depend on the chamber shape as well as on its depth. As expected, for deep chambers the results are similar to the case of a central unloading event. Eruption is then favoured for spherical and oblate chambers, whereas it is inhibited for prolate chambers. For shallow chambers, the magma pressure change ΔPc is highest for prolate ellipsoids, but tends to zero for oblate ellipsoids. Threshold pressure changes ΔPr are negative for prolate shape and positive for oblate shape. In this case, an unloading event prevents rupture initiation, except for a prolate ellipsoidal chamber filled with compressible magma.
(c) Application to Grímsvötn and Katla volcano
The 2004 eruption of Grímsvötn, a subglacial volcano located under the centre of Vatnajökull ice cap, was preceded by a jökulhlaup, a glacial outburst flood of 0.5 km3. Thirty per cent of the total discharge occurred before the eruption, corresponding to a volume of 0.15 km3 or 15 m drop in lake level (H. Björnsson & F. Pálsson 2008, personal communication). Modelling indicates that the sudden removal of a 1.8 km radius, 15 m water load, above a 10 km3 magma chamber at 2.5 km depth, might trigger dyke initiation by reducing the threshold by 10–150 kPa. The effect depends strongly on the chamber shape and magma compressibility, being largest for a spherical chamber filled with compressible magma. We, thus, confirm that small and sudden surface load changes can trigger eruptions as observed at Grímsvötn in 2004 but only if the magma chamber is close to failure conditions.
We also applied our model to Katla volcano, Iceland, which is covered by the Mýrdalsjökull ice cap. Ice load variations include long-term thinning, as well as an annual load cycle, with up to 6 m change in snow thickness from winter to summer (Gudmundsson et al. 2007). Our model predicts that, in the case of a spherical or horizontally elongated magma chamber, eruptions are more probable when the seasonal snow cover is smallest (Pinel et al. 2009; Albino et al. in press). This triggering effect is small, around a few kilopascals, but appears consistent with the fact that all the last nine major historical eruptions at Katla occurred during the summer period (e.g. Eliasson et al. 2006). A long-term ice thinning due to global warming is also occurring, mainly at the periphery of the Mýrdalsjökull ice cap. Our elastic model predicts that the short-term effect of such a process is to reduce the likelihood of eruptions at Katla. However, a model including a viscoelastic upper mantle and lower crust is needed to fully evaluate the influence of long-term ice thinning on magmatic systems.
5. Discussion and conclusions
Pressure drives the thermodynamics of magma generation in the mantle and the mechanics of magmatic ‘hydraulic’ fracturing in the crust. Ice retreat induces pressure release in the mantle capable of magma generation, depending on the size of ice caps. For estimated thinning rates of the Vatnajökull ice cap of approximately 50 cm yr−1, the pressure release is of the order of 0.5–1.5 kPa yr−1 throughout a large part of the melting regime. The velocity of upwelling mantle under Iceland due to tectonic and mantle plume processes is uncertain but has been estimated to be of the order of 1–10 cm yr−1 (Ito 2001), corresponding to a pressure release of up to 3.2 kPa yr−1 for a density of 3300 kg m−3. The pressure release caused by ice retreat is therefore highly significant and new magma is generated due to the warmer climate. We estimate that the induced minimum rate of melt generation under the Vatnajökull ice cap is 0.014 km−3 yr−1. The effect is an order of magnitude smaller than the inferred pressure decrease of up to 19 kPa yr−1 during the last deglaciation of Iceland (Jull & McKenzie 1996). The mode and time scale of transport of the additional magma from approximately 50–100 km depth in the melting regime towards the surface is highly uncertain, and may take longer than decades or centuries. An average melt extraction velocity of more than 50 m yr−1 (Maclennan et al. 2002) would suggest a transport time of less than 1000 years from a depth of 50 km. Although some of the magma may end up as intrusions in the crust, more frequent or voluminous eruptive activity is expected. Such an anomaly in eruptive activity may be difficult to detect because of the widely episodic nature of volcanism.
At shallow depth ice retreat influences both the magma pressure and the threshold pressure required for failure of magma chambers. If the values for the two pressures approach each other, the likelihood of dyke injection or an eruption is enhanced. As the vertical and horizontal stresses induced by ice retreat are different, the shape of the magma chamber plays a critical role. In the case of spherical and oblate-shaped chambers (sill-like) with highly compressible magma, under a retreating surface load that is wide compared with the depth to an underlying magma chamber, the magmatic systems will be driven towards failure at a rate comparable to the pressure release due to surface unloading. This may correspond to a few kilopascals per year, about three orders of magnitude less than the tensile stress of rocks. In general, the effect of present-day ice thinning on the failure of shallow magma chambers should therefore be expected as a small modulation on their natural activity. The model presented here for the influence on shallow magma chambers considers only elastic rheology, which is a limitation for sustained ice retreat. Future work requires consideration of viscous relaxation, as already done when quantifying pressure variation effects in melting areas. Moreover, we only considered pressure variations within magma chambers resulting from surface load changes, but additional melt available in melting zones may induce pressure increase in a magma chamber. Magma routes may be affected. Under long-term surface unloading the behaviour of the elastic layer will increase the pressure near the bottom of the plate, closing conduits, while decreasing them further up, facilitating opening of magma conduits towards the surface.
Support from the University of Iceland Research Fund and the Icelandic Research Fund is acknowledged, as well as funds for EU project VOLUME (Contract 18471). We are very grateful for helpful reviews from Lionel Wilson and Bill McGuire.
One contribution of 15 to a Theme Issue ‘Climate forcing of geological and geomorphological hazards’.
- © 2010 The Royal Society