# Unexpected large eruptions from buoyant magma bodies within viscoelastic crust

Freysteinn Sigmundsson, Virginie Pinel, Ronni Grapenthin, Andrew Hooper, Sæmundur A Halldórsson, Páll Einarsson, Benedikt G Ófeigsson, Elías R Heimisson, Kristín Jónsdóttir, Magnús T Gudmundsson, et al.

► To cite this version:

Freysteinn Sigmundsson, Virginie Pinel, Ronni Grapenthin, Andrew Hooper, Sæmundur A Halldórsson, et al. Unexpected large eruptions from buoyant magma bodies within viscoelastic crust. Nature Communications, 2020, 11 (1), pp.2403. 10.1038/s41467-020-16054-6 ird-04996275

HAL Id: ird-04996275

https://ird.hal.science/ird-04996275v1

Submitted on 18 Mar 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
ARTICLE

https://doi.org/10.1038/s41467-020-16054-6

OPEN

<figure><img src="image_2.png" /></figure>

Check for updates

# Unexpected large eruptions from buoyant magma bodies within viscoelastic crust

Freysteinn Sigmundsson ID <fnref n="1" />✉, Virginie Pinel ID ², Ronni Grapenthin ID ³, Andrew Hooper ID ⁴,
Sæmundur A. Halldórsson<fnref n="1" />, Páll Einarsson ID <fnref n="1" />, Benedikt G. Ófeigsson⁵, Elías R. Heimisson ID ⁶, Kristín Jónsdóttir⁵,
Magnús T. Gudmundsson ID <fnref n="1" />, Kristín Vogfjörd⁵, Michelle Parks⁵, Siqi Li ID <fnref n="1" />, Vincent Drouin ID ⁷,
Halldór Geirsson ID <fnref n="1" />, Stéphanie Dumont ID ⁸, Hildur M. Fridriksdottir⁵, Gunnar B. Gudmundsson⁵,
Tim J. Wright ID ⁴ & Tadashi Yamasaki ID ⁹

Large volume effusive eruptions with relatively minor observed precursory signals are at odds with widely used models to interpret volcano deformation. Here we propose a new modelling framework that resolves this discrepancy by accounting for magma buoyancy, viscoelastic crustal properties, and sustained magma channels. At low magma accumulation rates, the stability of deep magma bodies is governed by the magma-host rock density contrast and the magma body thickness. During eruptions, inelastic processes including magma mush erosion and thermal effects, can form a sustained channel that supports magma flow, driven by the pressure difference between the magma body and surface vents. At failure onset, it may be difficult to forecast the final eruption volume; pressure in a magma body may drop well below the lithostatic load, create under-pressure and initiate a caldera collapse, despite only modest precursors.

<footnote n="1">¹Nordic Volcanological Center, Institute of Earth Sciences, University of Iceland, IS-101 Reykjavik, Iceland. ²Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTERre, 38000 Grenoble, France. ³Geophysical Institute & Dept. of Geosciences, University of Alaska Fairbanks, 2156 Koyukuk Drive, Fairbanks AK-99775, USA. ⁴COMET, School of Earth and Environment, University of Leeds, Leeds, UK. ⁵Icelandic Meteorological Office, Reykjavik, Iceland. ⁶Seismological Laboratory, California Institute of Technology, Pasadena, CA, USA. ⁷Iceland GeoSurvey, Reykjavik, Iceland. ⁸Instituto Dom Luiz - University of Beira Interior, Covilhã, Portugal. ⁹Geological Survey of Japan, AIST, 1-1-1 Higashi, Tsukuba, Ibaraki 305-8567, Japan. ✉email: fs@hi.is</footnote>
Some of the largest effusive eruptions on Earth in the last 10,000 years were basaltic fissure eruptions in Iceland, with volumes ranging from 15–25 km³ (ref. ¹⁻³). More recently, during the 2014–2015 Holuhraun eruption ~2 km³ of magma were removed from a magma body⁴,⁵. Yet, precursory signals prior to magma body failure were weak, which raises questions about how magma bodies supporting such large eruptions are assembled, how eruptions are initiated and how large volumes of magma are extracted without major precursors.

The 1783–1784 Laki and the 2014–2015 Holuhraun fissure eruptions in Iceland respectively produced 15 km³ and 1.4 km³ of basaltic lava, formed by concurrent mixing and crystallisation of diverse primary mantle-derived melts⁶,⁷. Such moderate to large effusive basaltic eruptions and the preceding long-term magma accumulation, modulated by tectonics and mantle plume activity, provide a challenge for commonly used mechanical volcano models. This is clearly illustrated by the 2014–2015 activity in the Bárðarbunga volcanic system, when gradual caldera collapse occurred over a period of 6 months in response to drainage of

magma from below the caldera, initially along a lateral dike⁵,⁸ that fed a major effusive eruption 45 km away at Holuhraun (Supplementary Note 1). Several lines of evidence (ground deformation, petrology and seismic activity)⁴,⁵ suggest that the magma came from a single magma body located at a depth of 10 ± 3 km beneath the surface of the ice-filled caldera, below the brittle-ductile transition⁴,⁷,⁹. Although deformation and seismicity suggest some magma inflow prior to the eruption (Fig. 1, Supplementary Note 1, Supplementary Fig. 1), the amount of deformation suggests the short-term pre-eruptive magma inflow volume was minor compared to the volume drained during the eruption, unless the magma was highly compressible, damping observable inflation¹⁰. This is unlikely given the storage conditions and magma composition. This implies that overpressure (pressure in excess of lithostatic pressure) due to magma inflow at the onset of the 2014–2015 events was small compared to the co-eruptive pressure drop; a feature not in line with commonly used mechanical volcano models.

The host rock surrounding a magma body is typically modelled as elastic, responding instantaneously to a pressure change in a

<figure><img src="image_3.png" /></figure>

**Fig. 1** Seismicity and deformation at Bárðarbunga. **a** Earthquakes versus time (1 September 2013–15 August 2014; before diking and eruption) plotted as impulses scaled with magnitude (right axis). Earthquakes prior to M3.7 event on 16 May 2014 shown in blue and red afterwards. Also shown is cumulative seismic moment (shaded in grey; left axis), and horizontal displacement in direction N279.4°E (yellow dots) at GPS-station VONC from detrended time series (Supplementary Fig. 1). Error bars (1σ) in grey. **b** Inferred location of the earthquakes (Methods) shown in **a**, with earthquakes prior to M3.7 event on 16 May 2014 shown in blue and red afterwards. Small map of Iceland shows the study area outlined with a red box, fissure swarms⁷⁶ with grey shading, and oval outlines showing central volcanoes and calderas⁷⁶. **c** Location of M > 4.6 earthquakes during the caldera collapse (Methods). Note aseismic segments of the caldera. Background map shows shaded surface (grey) and Vatnajökull icecap (white) topography from the ArcticDEM database⁷⁷, and outlines of the Bárðarbunga central volcano (oval shape) and its caldera⁴². Straight lines show segments of the lateral dike that formed and black open circles are ice cauldrons⁵.
source dependent on geometry and elastic stiffness. The density difference between magma and host rock is usually also neglected. In such models, the erupted volume can only be a small fraction of total magma chamber volume¹¹,¹². As a result, the generation of large eruptive volumes would require very large magma reservoirs, high magma compressibility, or high overpressure that may far exceed the effective tensile strength of the Earth. High overpressure is, however, expected to induce large pre-eruptive surface displacements, not observed at, for instance, Bárðarbunga. Furthermore, in elastic models neglecting buoyancy effects, an eruption should stop before any large depressurisation of a magma reservoir can be reached thus preventing the formation of a caldera¹³,¹⁴. By contrast, in models where the Earth rheology is viscoelastic, a pressure change induces an initial elastic response followed by ductile relaxation, with relaxation rates that depend on the viscosity of the medium. While long-lived magma bodies may develop a local thermal aureole, modelled as a pressure source with a surrounding viscoelastic shell¹⁵,¹⁶, magma is often stored in multiple sills on its way towards the surface, with progressive compositional evolution¹⁷,¹⁸. For deep bodies, embedded within sufficiently hot media, we can expect entire crustal sections to behave in a viscoelastic manner¹⁹. Such long-term ductile behaviour of the Earth also gives rise to glacial isostatic adjustment (GIA) (ref. ²⁰).

While viscoelasticity is important to consider over the long timescales of magmatic recharge, over the short timescales associated with eruption and brittle failure, the host rock can be approximated as elastic. Furthermore, buoyancy of magma is important to consider when modelling eruption onset, as it has been suggested to play a key role in initiating super-eruptions²¹. Pressurisation associated with magma injection, on the other hand, is thought to be responsible for relatively small and frequent eruptions²¹.

Here we consider the consequences of buoyant magma accumulating in a viscoelastic medium, followed by magma draining through sustained magma conduits that allow magma to flow, despite development of under-pressure (pressure less than lithostatic) in the magma body. We develop a physical model that explains this paradox, which is not addressed by widely used volcano deformation models²² that only incorporate elastic crustal behaviour without considering density contrasts between magma and crust. But if magma accumulation occurs over a long timescale compared to the viscoelastic relaxation time of the host rock, and the magma is less dense than the host rock, the dominant factor driving magma bodies towards failure is magma buoyancy. Our model can be used to determine the conditions under which caldera collapse may initiate without significant precursors and drive large volumes of magma from depth to the surface. The model can account for the sequence of events that led to the 2014-2015 Bárðarbunga caldera collapse⁴ and its associated unrest and eruption⁵.

## Results

**Modelling framework.** When magma injects into a new or pre-existing magma body embedded in viscoelastic host rock, characterised by a linear viscoelastic rheology (Maxwell material), the stress relaxation timescale ($\tau$) is

$$ \tau = \gamma \frac{\eta}{\mu}. \qquad (1) $$

Here the ratio of host rock viscosity ($\eta$) and shear modulus ($\mu$) is the so-called Maxwell time, and $\gamma$ is a factor dependent on both the geometry of the magma body and the viscoelastic structure¹⁵,¹⁶,¹⁹,²³. In general, the viscosity is highly temperature-dependent and the uppermost few kilometres of the crust behave effectively as an elastic material. Appropriate regional values for

the lower crust in Iceland are $\eta = (2-8) \times 10^{18}$ Pa s, as inferred from GIA in response to present day ice thinning²⁰, and $\mu = 10-30$ GPa (ref. ²⁴), giving a regional Maxwell time on the order of 2–25 years (for $\gamma = 1$). While host rock viscosity at active volcanic plumbing systems is expected to be drastically lower than average regional values, due to higher temperature, $\gamma$ is on the order of 10 or higher for magmatic sills emplaced in a viscoelastic medium under an elastic layer¹⁹. Considering these competing effects, the relaxation time of host rock next to sill-like bodies due to magma inflow may be comparable to, or shorter than, the regional Maxwell time inferred from GIA; short compared to the time of magma accumulation over inter-eruptive periods lasting decades to centuries. At low emplacement rates, it follows that viscoelastic host rock rheology allows large magma volumes to accumulate over long periods without developing large overpressure, as stresses generated around a magma body due to inflow of new magma are relaxed.

In a viscoelastic material the density difference between magma and crustal rocks, $\Delta\rho$, generates a buoyancy force that is directed vertically upward for a magma body as a whole (Figs. 2 and 3), with magnitude given by Archimedes law. At the boundary of the magma body, the buoyancy pressure is locally $\Delta\rho gh$, where $g$ is acceleration of gravity, and $h$ is the height above the base of the magma body (Methods). Above a magma body, the resulting buoyancy-induced stresses are similar to those created by magma inflow (Fig. 4). Pressure at the boundary of a magma body results from the combined effects of buoyancy and pressure changes arising from magma inflow or outflow ($\Delta P_{magma body,flow}$).

A failure criterion for a buoyant magma body of uniform density, impacted by external stress, can be written as:

$$ \Delta\rho gh + \Delta P_{magma body,flow} \geq \sigma_{failure} + \sigma_{external} \qquad (2) $$

where $\sigma_{failure}$ is the failure limit. The external stress ($\sigma_{external}$; negative for tension) is the component of deviatoric stress at the boundary of the magma body due to all external processes acting perpendicular to the plane of failure. This can result from a combination of processes, including tectonic, short- to long-term surface loading such as GIA, and topographic effects. We refer to Eq. (2) as the general failure criterion for a magma body. Failure is approached when either (i) a buoyant magma body increases in thickness, (ii) its density contrast increases through magma evolution, (iii) its pressure increases due to influx of magma, (iv) external stresses bring the deviatoric state of the host rock closer to failure, or (v) a combination of any of these. If the rupture is governed by tensile failure¹¹ and in agreement with magma-filled dikes and hydrofracture mechanical models²⁵, then $\sigma_{failure}$ is independent of depth and equals the effective large-scale tensile strength ($T_{eff}$) multiplied by a geometrical factor $\alpha$ ($\sigma_{failure} = \alpha T_{eff}$). The geometrical factor $\alpha$ depends on the magma body geometry and takes a value close to 1 or less for a horizontally elongated shape¹². The effective large-scale tensile strength ($T_{eff}$) is on the order of few MPa (refs. ²⁶,²⁷), up to an order of magnitude lower than laboratory values for intact rocks²⁸. Alternatively, if failure occurs by shear faulting, as sometimes proposed in numerical gravity-loaded elasto-plastic models, the failure limit is up to an order of magnitude higher (several tens of MPa) (refs. ²⁹,³⁰). Buoyancy effects limit the maximum magma body thickness. For example, in the case of tensile failure, if $\Delta\rho = 270$ kg/m³, $T_{eff} = 2.5$ MPa and in the case $\alpha = 1$ and $\sigma_{external} = 0$, the thickness has to be <1 km (Fig. 2) in order to ensure stability. The tectonic setting is also important; at divergent plate boundaries the stress field supports less internal pressure accumulation before failure of the host rock surrounding the magma body than at convergent margins.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Model and stress field due to buoyancy. **a** Schematic model of a magma body residing in viscoelastic material below an elastic layer. The magma is less dense than the host rock and therefore buoyant. The density difference between magma and host rock is Δρ, and h is the elevation above the bottom of the magma body. **b** Stress due to buoyancy around a horizontal ellipsoid, with a 4.3-km semi-major axis and a 0.4-km semi-minor axis, embedded in an elastic medium with a Poisson's ratio of 0.25. The density contrast is set to 270 kg/m³. Pressure equal to Δρgh, with g being the acceleration of gravity, is applied at its boundary (Fig. 3). The colour scale shows the amplitude of the minimum compressive stress (negative values for tension) and dashes represent the direction of maximum compressive stress (σ₁). The crustal volume immediately above the magma body is under compression, but the edge of the body is affected by a large tensile stress allowing for dike initiation. Above the body, the stress field is similar to the one induced by overpressure resulting from magma inflow (Fig. 4).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Buoyancy modelled as traction on boundary of a magma body.</figcaption></figure>

Overpressure is calculated with reference to the deepest point of the magma body. This pressure is transferred to a surface force at the boundary, in a direction normal to it. The vertically directed upward force on the column shown in red is ΔAΔρg (h₂-h₁). Here ΔA is a horizontal cross-sectional area element, equal to (cosφ)δA, where φ is the local angle that the boundary of the magma body makes with respect to horizontal. The integrated upward force due to all columns in the magma body gives an overall upward force following Archimedes law. For this geometry, the integrated horizontal force on the entire boundary will be zero (cancels out for this geometry because of symmetry).

Following magma body failure, a channel will form enabling magma outflow. Its shape depends on several processes and factors, such as tectonic and topographic stresses acting on the channel³¹, zones of weaknesses such as caldera bounding faults, and the buoyancy pressure³². If the channel has vertical extent from depth D₁ to depth D₂, then the driving pressure (ΔP<sub>driving</sub>) for upward flow in the channel can be written³³:

$$ \Delta P_{\text{driving}} = \Delta \rho g h + \Delta P_{\text{magmabody, flow}} + \int_{D_1}^{D_2} \Delta \rho g dz. \quad (3) $$

The first two terms are internal to the magma body and the same as in Eq. (2): the magma body buoyancy, and the effects of magma flow to and from it. The third term captures effects of magma buoyancy in the channel. If crustal and magma density are constant with depth, then it simplifies to ΔρgH, where H is the height of the channel equal to D₂-D₁. Equation (3) applies to

magma flow in a channel only after the establishment of a stable geometry and before its final phase of closure, when flux decreases and the magma cools and solidifies. We refer to this as a sustained channel. Within it, upward magma flow will continue as long as the driving pressure remains positive. At the channel top, the overpressure compared to lithostatic pressure (Fig. 5) equals the magma flow driving pressure, neglecting effects of magma velocity and friction (Methods). Near a deflating magma body, the channel may experience high compressive stress resulting in its closure. Exceptions may include situations influenced by inelastic processes. One such process is addition of magma mush (consisting of both crystals and some interstitial liquid) to the magmatic liquid flowing out of a magma body<sup>7,34</sup>, evidenced in some cases by the presence of mineral assemblages, including macrocrysts, that are not in equilibrium with the host lava (Methods). A continuous entrainment process is inferred if the residence time of macrocrysts in the magma, prior to eruption, is short compared to the duration of eruption, and macrocrysts are relatively uniformly distributed in lava, as found for example, in the 2014–2015 Holuhraun lava<sup>7,35</sup>. Relatively uniform content of macrocrysts, despite variable magma flow rates during an eruption, indicates furthermore that the amount of material added to the lava can scale with the magma flow rate. We refer to this process as “magma mush erosion”, which generates space for a magma channel via a process not included in elastic models, supporting the formation of a sustained channel. Thermal effects may also contribute: solidification of thin parts of a magma-filled crack due to cooling along the host rock interface will help to keep open the wider parts of a dike through bridging effects. Conversely, the wider parts of a dike may grow even wider through thermal erosion. This is the same process by which eruptive fissures that initially produce “curtains of fire” eruptions evolve with time; some parts of fissures close while others may widen and focus activity³⁶. Tectonic effects, e.g., the release of tectonic stress in extensional environments may also contribute to the formation of sustained magma channels.

The deep roots of volcanoes may host a complex system of sills and dikes forming transcrustal magmatic systems with variable amounts of magma mush<sup>18,34</sup>. This can be referred to as magma
<figure><img src="image_6.png" /><figcaption>Fig. 4 Comparison of stress field at magma bodies due to buoyancy and magma flow. a Zoom of Fig. 2b showing stress field set up due to buoyancy of magma body. b Stress field due to 2.5 MPa pressure increase caused by magma inflow, for comparison. The colour scale shows the amplitude of the minimum compressive stress (negative values for tension) and the dashes represent the direction of maximum compressive stress (σ₁).</figcaption></figure>

domain³⁷,³⁸ and during an eruption magma extraction may occur from one or more magma bodies within the magma domain. In the instance of magma draining from a single magma body, its size and shape will provide a strong control on the magma flow. Calderas are often several kilometres or more in diameter, implying similar dimensions of underlying magma bodies. When magma is buoyant, the general failure criterion (Eq. 2) requires these bodies to be thin compared to their lateral extent if tensile failure applies. An oblate horizontal ellipsoid represents an appropriate geometric model. On a short timescale, the relation between the volume of magma outflow, ΔV<sub>magma</sub>, and the pressure change within an oblate ellipsoid, filled with a compressible fluid and embedded in an elastic halfspace with a Poisson's ratio of 0.25 and shear modulus μ is (Methods):

$$ \Delta V_{\text{magma}} = \left(\frac{2}{\mu} + \frac{4\pi b}{3 a} \frac{1}{k}\right) a^3 \Delta P_{\text{magmabody, flow}} \quad (4) $$

where 1/k is the magma compressibility, a is the semi-major axis of the ellipsoid and b is its semi-minor axis. The equation shows that for an oblate ellipsoid with thickness much smaller than lateral extent, the effects of magma compressibility may be much smaller than for a more spherical magma body where the semi-minor axis approaches the length of the semi-major axis, as the compressibility is scaled by b/a. After failure, and with the formation of a sustained magma channel from the source, the pressure reduction can be much larger than the initial over-pressure required for failure (Fig. 5).

Volatile content and the chemical composition of magma strongly modulate rheological properties of magma (including its viscosity and density), and in turn, significantly influence subsurface magma flow and eruption dynamics¹⁰,³⁹. For basaltic magma, the most important volatiles to consider are H₂O, CO₂ and S. For example, estimates of depth-dependent changes in magma density require information on magma composition and dissolved volatile content, in addition to an appropriate set of thermodynamic relationships describing how volatile distribution

between gas and melt changes as pressure varies⁴⁰. Given how magma density, ρ<sub>magma</sub>, changes with pressure, the magma compressibility, 1/k, can be inferred, as it equals (1/ρ<sub>magma</sub>) (∂ρ<sub>magma</sub>/∂p) where p is pressure. Silica content and concentration of volatiles have thus also a large influence on the compressibility of magma¹⁰,³⁹. In order to constrain the density difference between magma and the host rock, and thus infer the magma dynamics, a model of the crustal density profile is required.

**Application to Bárðarbunga volcano.** To resolve the discrepancy between minor precursory deformation and large volume eruption, we test our model against the 2014–2015 activity in the Bárðarbunga volcanic system. During this activity, moderately evolved olivine tholeiite magma is inferred to have drained from a single magma body located somewhere in the range of 7–13 km depth below the surface of the ice cover of Bárðarbunga⁴,⁷,⁹,³⁵. In particular, geobarometry suggests that the most probable pressure at which the magma resided was 230 ± 140 MPa (refs. ⁷,³⁵). Crystal-hosted melt inclusions in the erupted lava⁴¹ indicate that the magma contained about 0.5 wt.% H₂O, 500 ppm CO₂ and 1600 ppm S at 230 MPa prior to the eruption. We used two software programs to infer the magma density (Methods); one of which (the D-Compress software⁴⁰) is used to infer variation in magma density and compressibility versus depth resulting from volatile loss (Fig. 6). A generalised model for the crustal density profile in Iceland⁴² is modified to account for the ~1.8-km elevation of Bárðarbunga above sea level, divided into a 0.6-km-thick ice layer and 1.2-km-thick low-density hyaloclastite layer (Supplementary Text and Supplementary Fig. 2). The model gives crustal density versus depth and allows the determination of crustal pressure conditions (Fig. 6). Over the 150–230-MPa pressure range (~6–9 km depth below surface; ~4–7 km depth below sea level) the magma and crustal densities show minor variations; their difference being in the range of 270–300 kg/m³. For conservative modelling of the buoyancy effect, we initially take the lower limit of this range (Δρ = 270 kg/m³). This density
<figure><img src="image_7.png" /></figure>

**Fig. 5** Excess pressure in a magma plumbing system. Pressure history during inflation and deflation. Combined overpressure from magma flow and buoyancy pressure in a magma body (red, solid line considering a constant magma density and dashed line taking into consideration the magma density changes induced by gas exsolution) and in a location higher within a sustained magma channel (blue, solid line considering a constant magma density and dashed line taking into consideration the magma density changes induced by gas exsolution). The two curves are separated by $\Delta\rho gH$ (which is a constant value when volatiles are not considered). Note that volatiles are not expected to influence the overpressure evolution for a magma body located at 8 km or larger depth. Gas exsolution may cause the two curves to slightly diverge as deflation evolves, as such a process will always have a larger influence at shallower depth. Failure of the magma body occurs at $t = 0$ and caldera collapse (failure of caldera faults) begins in this example at day 8 after the failure of the magma body. For a long time prior to that (years), the pressure conditions are close to the magma body failure limit due to magma buoyancy. Overpressure in the magma body is limited according to the general failure criterion, here considered to correspond to be effective tensile strength of 2.5 MPa (black horizontal dotted line). In this example, slip on caldera ring faults begins at day 8 when critical conditions for caldera failure are reached (green horizontal dotted line). Following that, there is small change in pressure as piston collapse drives out a large volume of magma over months. See text for discussion.

difference is also considered in the crustal layers next to the magma body to constrain the magma channel dynamics.

Considering a short relaxation timescale and the location of the magma body within a viscoelastic domain, stresses due to magma accumulation prior to the events may have been continuously relaxed over the 150 years since the previous major Bárðarbunga fissure eruption¹. The Gjálp eruption in 1996 did, however, significantly modify conditions at Bárðarbunga. It was preceded by a 22-year period of ~15 magnitude 5+ earthquakes, the last one leading to a lateral dike propagating from Bárðarbunga. This 1996 dike intersected another magma body outside of the caldera⁴³, which provided the majority of erupted material. Ring fault activity occurred during this event, suggesting weak caldera faults. We infer that in August 2014 the general failure criterion (Eq. 2) was satisfied again at Bárðarbunga, resulting in buoyant magma draining from the underlying magma body (Fig. 7).

Contraction of a magma body under the caldera has been modelled⁴⁴ as a closing sill with lateral dimensions similar to the caldera. We consider a similar but simplified model geometry. We utilise a reference magma body and Earth model in our calculations, and then assess deviations from these. The reference magma body model is an oblate ellipsoid magma body at 10 km

depth below the surface with semi-major axis $a = 4.3$ km, so that the horizontal width of the ellipsoid is comparable to the dimensions of the Bárðarbunga caldera. We define the semi-minor axis such that tensile stress at the edge of the magma body does not exceed ~2.5 MPa due to buoyancy. This requires the maximum thickness of the magma body to be ~1 km or less. Here we set the semi-minor axis to be $b = 0.4$ km (Fig. 1). We set the compressibility ($1/k$) of the magma in the body to $0.1 \text{ GPa}^{-1}$, similar to the value derived from D-Compress at the depth of the magma body ($0.09 \text{ GPa}^{-1}$). The effective Young's modulus of the Icelandic crust has been estimated ($40 \pm 15$ GPa based on annual cycles in ground displacements due to snow and ice load changes²⁴. The corresponding shear modulus is $\mu = E/(2 + 2\nu) = (10-22)$ GPa. For the reference model we use the lower limit, $\mu = 10$ GPa, to reflect the expectation that crustal rigidity under volcanoes and in a magma domain may be lower than regional averages. Caldera collapse began on day 8 of the activity, heralded by the onset of M5+ earthquakes at the caldera⁴. According to ground deformation modelling⁵, ~0.3 km³ of magma had already flowed into the dike at that time. For this volume change and the parameter values above, the associated pressure change in the magma body, $\Delta P_{\text{magma body,flow}}$, is -15.8 MPa (Fig. 5).

Next, we consider a model for a magma channel that links the magma body under Bárðarbunga and the lateral dike that formed (Fig. 7). During the dike emplacement and eruption, elastic processes dominated, except for inelastic processes that may have created a sustained channel where Eq. (3) applies. Even if the entrained material from the crystal mush is only 1% of the drained magma volume, as we estimate for the 2014–2015 magma (Methods), then the volume of entrained magma mush would still be ~$20 \times 10^6 \text{ m}^3$, corresponding to, for example, a 5-km high conduit segment (dike) with a $4 \times 10^3 \text{ m}^2$ cross-sectional area (e.g., 4 km long and 1 m wide). We suggest a vertical upflow path along the eastern caldera boundary of Bárðarbunga (Supplementary Text). This initial flow path from the magma body may have been a dike that formed along a part of the caldera boundary where it was weak. At a higher level, within the brittle crust, where the magma is less buoyant and topographic and tectonic stress have major influence, a regional dike formed (Fig. 7). The maximum opening of the regional dike is at ~3 km depth (ref. 5) and we consider this to be the average depth to which magma needs to rise in our reference model. If magma flowed vertically upward from 10 km to 3 km depth below the surface (approximate depth of maximum opening in the dike⁵), the overpressure gain in the vertical channel ($\Delta\rho gH$) is 18.5 MPa for $\Delta\rho = 270 \text{ kg/m}^3$. This shows that the driving pressure can remain positive, allowing for flow despite the large pressure drop in the magma body (Fig. 5). Depressurisation in the model produces strong shear stress at the locations of potential ring faults at depth, which favours the initiation of a piston-type collapse (Fig. 8); large under-pressure allowed for the onset of movements on already weak caldera faults following previous events (Supplementary Text). Once the collapse was initiated, release of gravitational energy, associated with the subsiding piston⁴, drove out a large volume of magma along the previously created magma path. During the piston collapse, exponential decay in magma flow-rate occurred with an estimated pressure change of <2 MPa, comparable to the change in lithostatic pressure as the caldera subsided 66 m (ref. 4).

Although the values of the model parameters presented can explain the behaviour of the volcano in the framework of our model, their exact values are uncertain. If the magma – host rock density difference is in the range of 250–300 kg/m³ and the depth extent of the sustained channel is in the range of 4–10 km (the average rise of magma from the magma body to the effective dike inlet), then the difference in pressure between the magma body and the effective dike inlet is in the range of 12–30 MPa (the
<figure><img src="image_8.png" /><figcaption>**Fig. 6** Volatile content, magma and crustal properties versus pressure/depth. a Dissolved volatile contents (in wt.%) of the 2014–2015 Bárðarbunga magma versus pressure (in MPa, left vertical axis). b Bulk magma (black) and crustal (grey hatched) density (in kg/m³) and c bulk magma compressibility (in Pa⁻¹). Conversion from pressure to depth below sea level is shown at the right vertical axis. See also Supplementary text and Supplementary Fig. 2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 7** Schematic cross-section of a model for the 2014–2015 Bárðarbunga collapse. A sill-like magma body resides within a magma domain under Bárðarbunga³⁸,³⁹. Modest inflow of magma (recharge) occurred prior to August 2014, at least in the last few months prior to activity (see Fig. 1). Failure of the magma body occurred because of the combined effects of buoyancy and pressure due to magma flow. After failure, magma flowed up from an exit at the magma body towards an effective inlet to a regional lateral dike at a higher level in the crust. Macrocysts were entrained into the magma via erosion of a magma mush, contributing to the formation of a sustained magma channel. The magma channel may have been in the form of a dike along the caldera boundary (perpendicular to the plane of the figure). At shallower levels tectonic and topographic stresses dominated and a regional dike formed. Despite development of under-pressure in the magma body, the driving pressure for flow remained positive. The magma body involved was located beneath the long-term brittle-ductile transition located in the 6–8-km depth range below sea level, near the lower edge of the dike formed⁹. The structure of the magma domain is shown in a schematic manner. It may have other isolated magma bodies as shown and volumes of crystal mush; at least one volume of magma mush was involved in creating the erupted magma.</figcaption></figure>

offset of the upper set of curves in Fig. 5 relative to the lower set). For the low end of this range the eruption may stop early because the driving pressure for the flow (Eq. 3) falls quickly to zero. If the amount of combined overpressure due to buoyancy and magma flow at the time of failure was higher than the suggested value of 2.5 MPa, then the pressure curves shown in Fig. 5 would all be shifted to higher levels. In turn, the amount of pressure drop needed to create under-pressure and onset of caldera collapse would be larger. The gradient of the two sets of the pressure curves in the time between the onset of magma outflow and when the caldera collapse begins, depends on the all the parameters in Eq. (4). The larger the dimensions of the magma body, and the larger the compressibility, the lower the gradient of pressure drop versus time. However, for an oblate ellipsoid and the model parameters presented above, the dominant factor is the first term in the bracket in Eq. (4), which depends on the crustal shear modulus. The higher the shear modulus, the larger the pressure drop for a fixed magma outflow volume.

## Discussion

The model presented here explains how a magma body can develop an under-pressure that is an order of magnitude higher than the overpressure prior to initial failure, and yet facilitate substantial magma outflow, as we have demonstrated for the case of the 2014–2015 events at Bárðarbunga. This is different from the behaviour of a fully elastic model that does not consider gravity and where stresses never drop below those prior to the onset of magma accumulation, limiting the pressure drop to be equal to or less than the pre-eruptive pressure increase. In our model, a sustained conduit along part of a caldera boundary allows under-pressure to develop in the underlying magma body and, at the same time, reduces the overall resistance on the caldera boundary where it forms. Both factors facilitate a caldera collapse. The under-pressure that develops depends on the size of the magma body: the smaller the body the larger the under-pressure that develops for a fixed outflow volume. The development of under-pressure hinders upward flow of magma and may eventually stop it, but at the same time it increases the likelihood
<figure><img src="image_10.png" /><figcaption>**Fig. 8** Stress fields related to major drainage of magma. **a** Amplitude of the minimum compressive stress, shown by the colour scale with negative values for tension, induced by combined buoyancy and under-pressure due to magma outflow ($\Delta\rho gh + \Delta P_{\text{magma body,flow}}$) equal to -13.3 MPa, corresponding to the model presented for Bárðarbunga. Direction of maximum compressive stress ($\sigma_1$) is represented by dashes. Most of the medium is put under tension with maximum effect immediately above the magma body. **b** Shear stress (amplitude of $\sigma_1 - \sigma_3$; shown by the colour scale) for the same conditions as in **a**. The panels show that the stress field is favourable for caldera collapse initiation induced by the depressurisation.</figcaption></figure>

of a caldera collapse. Thus, the right conditions for initiation of a caldera collapse may rarely be met, consistent with their infrequent occurrence.

For water-poor basaltic magmas, exsolution of volatile elements from melts stored at mid- to deep-crustal levels (>8 km) has less of an effect on magma compressibility relative to water-rich, silicic magmas stored at shallow levels¹⁰. In addition, the presence of high-density fluid inclusions in the Holuhraun melt suggest that it became CO₂ saturated at mid- to lower-crustal pressures in excess of 400 MPa (ref. ⁴¹). A considerable portion of the primary CO₂ was therefore lost via open-system degassing before the Holuhraun magma accumulated in the inferred magma body at ~230 MPa pressure (ref. ⁴¹). This is consistent with other studies⁴⁵,⁴⁶, indicating CO₂ exsolution from basaltic melt under Iceland and its subsequent loss begins at ~25 km depth. Our depth-dependent volatile exsolution modelling (Fig. 6) further suggests that the loss of CO₂ from the melt was a near-continuous process, while both water (which profoundly affects magma compressibility) and sulfur exsolution only become significant at shallow depths (<500 m). As a result, changes in rheological properties of the Holuhraun magma driven by volatile loss mostly occurred in the topmost three kilometres of the crust, above the level of maximum dike opening. In this depth range, volatile loss has thus a major influence on magma ascent dynamics whereas at deeper levels the effects on dike opening are minor.

Our model is appropriate for deep magma accumulation, with magma flow rates below a threshold value (Methods). For Bárðarbunga, the model provides a direct link between the deep magma body, where accumulation and drainage occurred, and the shallower regional dike. Another example of well-studied deep magma accumulation in extensional tectonics is the Socorro Magma Body in the Rio Grande Rift. It is inferred to be ~150 m thick, residing at ~19 km depth, and has been inflating at low rates for decades⁴⁷⁻⁴⁹. One would expect similar processes to be at work there, even if the lithospheric properties are different from those in Iceland. Individual elements of the modelling framework presented here (viscoelastic host rock behaviour, magma

buoyancy, sustained magma channel, development of under-pressure to initiate caldera collapse) may help to understand large volume eruptions in general. The most recent large volume effusive eruptions on Earth include the rift eruption at Kilauea volcano in 2018 when ~0.8 km³ of magma erupted in relation to summit caldera collapse³⁸,⁵⁰, and the ongoing submarine eruption east of the Mayotte island, Indian Ocean, which began in May 2018 and where more than 5 km³ of magma have erupted⁵¹⁻⁵³. At Kilauea, interpretation of geodetic measurements finds pressure drop prior to onset of caldera collapse to be ~17 MPa (ref. ⁵⁰), similar to what we find for the Bárðarbunga collapse. In the activity east of Mayotte, drainage from a deep magma body (>20 km depth) in the mantle is inferred⁵¹,⁵³ where the host rock is viscoelastic. Studies of the eruptive products show the magma was buoyant⁵². More general, the reawakening of dormant volcanoes by inflow and accumulation of magma to shallow depth may result from failure of deep-seated buoyant magma bodies in viscoelastic host rock, residing in the lower crust or the mantle. If magma buoyancy plays a role in bringing magma bodies towards the general failure criterion, then our model suggests failure of such bodies can occur with little warning, and once initiated, can result in unanticipated large lava production.

## Methods

**Location of earthquakes.** Earthquakes shown in Fig. 1 are from routine locations from the national seismic network of Iceland operated by the Icelandic Meteorological Office, shifted southward⁴, to correct for inadequacy of the velocity model that does not account for lateral variations in velocity structure. The southward shift decreases linearly from ~2.6 km at the northern caldera rim to 0.5 km at the southern rim. This shift is applied both to earthquakes in the pre-eruptive period and the diking/eruptive period shown in Fig. 1b, c, respectively. Formal location errors were significantly improved in second half of August 2014 as six new seismic stations were deployed within 45 km radius of the caldera (seven stations before). After late August 2014 the formal epicentral location errors are typically up to 1 km, whereas before they are up to 5 km.

**Buoyancy pressure.** The effect of magma buoyancy on dynamics of spherical magma bodies has been considered before²¹,⁵⁴⁻⁵⁶, but is here expanded to magma bodies of arbitrary shape. The total upward force acting on a magma body of lower density than the surrounding viscoelastic crust is found, according to Archimedes
law, by considering the integrated body force due to gravity acting on the magma body minus the gravity force of the crustal material that the magma body has replaced. At the magma body boundary this upward directed body force is transferred to a surface force acting on the surrounding crust. Figures 2 and 3 show the case for an ellipsoidal magma body. For the case of uniform density contrast, the total upward force that the magma body exerts on the surrounding viscoelastic crust is $\Delta\rho gV$, where $\Delta\rho$ is the difference in density between magma and crust, g the gravitational acceleration, and $V$ the volume of the magma body. For variable density, $\Delta\rho$ is the difference of the average densities of the magma and the crust. As the magma body is embedded in a material that behaves as a viscous fluid on a long timescale, it will have the tendency to "float" upwards and its shape may be distorted. Stokes flow for bodies emerged in viscous fluids, and its extension to a fluid body emplaced within another fluid body suggest that upward movement velocity⁵⁷,⁵⁸ will be on the order of $(1/3)\beta\Delta\rho g\eta^2/\eta$, where $a$ is the semi-major axis for an ellipsoidal magma body and $\beta$ is a scaling factor that takes into consideration the aspect ratio $b/a$ (b being the semi-minor axis of the magma body). If $a=b$ (spherical body), then $\beta=1$. If $b<a$ (ellipsoid) then $\beta$ is <1. Assuming such a body within the Icelandic lower crust with viscosity of $5 \times 10^{18}$ Pa s, density difference $\Delta\rho$ equal to 270 kg/m³, and $a=4.3$ km, the upward rate of movement can be expected to be <10 cm/yr. This is slow enough that a magma body receiving inflow from depth is likely to transfer magma from deeper levels to shallower levels via feeder dikes, rather than diapiric rise of whole sills through the crust. Using the same values, the growth rate of Rayleigh-Taylor instabilities is expected to remain smaller than a few millimetres per century⁵⁹, confirming that the ascent of magma caused by buoyancy will be effective only when the magma body thickness will be large enough to overcome the surrounding medium resistance to rupture.

**Viscous pressure drop induced by magma flow.** Equation 3 gives the pressure at the top of the vertical magma channel neglecting the viscous pressure drop induced by the magma flow. This pressure drop can be estimated⁶⁰ around $\eta_m v H / w^2$, where $\eta_m$ is the magma viscosity, $v$ the magma velocity and $w$ the channel width. Taking $\eta = 100$ Pa s, $w=1$ m, $v=1$ m/s and $H=7$ km, the viscous pressure drop is ~0.7 MPa or <4% of the buoyancy gain over the channel height (18.5 MPa; see main body of text). The overestimation of the pressure at the horizontal dike inlet resulting from neglecting the magma flow is thus expected to be small.

**Sustained magma channel.** Mantle-derived lavas containing plagioclase and other mineral assemblages, too primitive to be in equilibrium with the host lava, are widely distributed around the globe³,⁶¹. One of the most striking examples of such crystal-melt disequilibria are so called plagioclase-ultra-phyric basalts which are found both within the Neovolcanic zone and throughout the near 16 Ma lava pile of Iceland⁶². These lavas are characterised by an abundance (often >10 vol%) of large crystals (>~0.5 mm; macrocrystals) with high relative proportion of plagioclase to olivine. A widespread view is that these macrocrystals crystallised at middle to deep crustal levels and accumulated in a crystal mush layer from which they were picked up by their host/carrier melt during ascent. Well documented examples of Icelandic Holocene lavas, where prior studies have demonstrated the importance of such a magma mush layer, include the 8600 BP Pjórsárhraun lava³, the 1783 Laki lava⁶,⁶³, and more recently the 2014–15 Holuhraun lava⁷,³⁵. Consideration of elemental diffusion in these crystals and their melt inclusions demonstrates that timescales associated with the crystal entrainment process are likely to be short, or equivalent to timescales of these eruptions, i.e., a few weeks to months³,⁷,⁶⁴. This suggests that the crystals have been brought into the magma only shortly before or during these eruptions. In line with these observations, we suggest that an important process for formation of sustained magma channels is the addition of mush crystals to a magmatic liquid as it leaves a magma body and rises towards the surface. In this case, a near-continuous addition of mush crystals into the host/carrier melt⁶³, in proportion to the rate of magma flow, leads to uniform crystal content throughout the erupted magma, despite variable rate of flow. We suggest that erosion of such a crystal mush layer may ultimately lead to formation of sustained openings, which are required for the extraction of a melt from a magma body. The reduction in volumetric stress which happens above and below a deflating sill (Fig. 8) may indeed make the crystal mush more easily erodible and therefore accommodate this process.

For the 2014–2015 Holuhraun lava, the amount of macrocrystals, defined as minerals with long axes >1 mm in hand specimen, is relatively small compared to some older Bárðarbunga lavas. Halldórrsson et al.³⁵, and later work estimated the range to be from ~1% and up to no more than 5% in terms of volume. Yet, even the lower bound of 1% is sufficient to support the role of magma mush erosion in the formation of a sustained channel in the 2014–2015 events.

**Horizontal oblate ellipsoid filled with compressible fluid.** The volume change of a penny-shaped crack of radius $a$, embedded in an elastic halfspace with shear modulus $\mu$ and Poisson's ratio 0.25, in response to a pressure change $\Delta P$ is⁶⁵

$$ \Delta V_{\text{penny}} = \frac{2a^3}{\mu} \Delta P. \qquad (5) $$

The volume change of a fluid body in response to pressure change due to its compressibility is $V\Delta P/k$, where $1/k$ is the fluid compressibility. If the fluid is in the

form of an ellipsoid with semi-major axis $a$ and semi-minor axis $b$, then its volume is $(4/3)\pi a^2 b$. In the case of magma outflow from an ellipsoidal body that takes the form of a penny-shaped crack ($b \ll a$), the volume of magma outflow equals the combined volume change of a penny-shaped crack and the volume change due to compressibility. Adding these two terms together gives Eq. (4).

**Numerical modelling.** The conceptual model presented here relies on the ductile behaviour of the lower crust surrounding a magma body. Nevertheless, we use numerical models considering elastic rheology to support it. This approach is justified since the buoyancy force in a viscoelastic medium, contrary to the overpressure induced by magma inflow, cannot be relaxed through time by a change in volume or shape. Its effect can be quantified using an elastic medium and taking the reference state of stress as lithostatic (fully relaxed). The stress field resulting from a pressurised buoyant body is thus determined from the equations for linear elasticity which we solved using the Finite Element Method in axisymmetrical geometry with the COMSOL software. We use a mesh of about 120,000 triangular units that is refined around the magma body (maximum element size of 20 m next to the magma body). No displacement is allowed at the lateral and bottom boundaries, which are set at a distance equal to 20 times the depth of the magma body whereas the upper boundary is considered a free surface. Stress due to buoyancy is similar to the one induced by an overpressurized magma body (Fig. 4). When the combined stress exceeds a given threshold, it favours tensile rupture of the medium consistent with magma propagation through the medium. Within a ductile medium, tensile stress is expected to favour rupture and magma propagation through sheet intrusions the same way as in an elastic medium except for a difference of the tip shape, the tip being less sharp in a ductile medium⁶⁶.

**Density of magma.** The density of the 2014–15 Holuhraun magma was calculated using both the D-Compress as well as the Petrolog3 software⁶⁷,⁶⁸. In the latter case, we assumed that whole-rock composition of the lava flow represents the melt compositions. From the average of 62 whole-rock analyses³⁵, a density value of 2722 ± 5 kg/m³ was obtained for the 2014–15 Holuhraun melt at 230 MPa (corresponding to 2.3 kbar). Pressure dependency of these calculations is negligible and within the error introduced by considering standard deviation of the mean value of 62 whole-rock analyses adopted for the calculations. A fixed $Fe^{3+}/Fe^{total}$ ratio of 0.15 was assumed, which is considered a representative value for evolved rift-associated basalts in Iceland. Water content of 0.5 wt% was selected as it best characterises prior-degassing H₂O levels of the Holuhraun magma at ~230 MPa. The inferred density value is in good agreement with values obtained for aphyric lavas of similar compositions and volume from North Iceland⁴⁵.

**Density of the crust at Bárðarbunga.** The Icelandic crust is 3–4 times thicker than normal oceanic crust⁶⁹. However, the upper crust in Iceland is similar in structure to that of the normal oceanic crust, with high gradients in both seismic P-wave velocities and density in the upmost 5 km, while the gradients are much lower at greater depth⁷⁰. Two density models for basaltic crust, based on empirical data, are those of Carlson and Herrig⁷¹ and Christensen and Wilkens⁷². The latter model is based on velocity-density systematics from a 1.9-km deep drillhole in east Iceland. These models have been used to constrain the sources of gravity anomalies in the upper crust in the Bárðarbunga region⁴². We use these models together with ice cover information to obtain an estimate of the density-lithostatic pressure systematics for Bárðarbunga. The Carlson and Herrig model provides a minimum density estimate and the Christensen and Wilkens model a maximum. We use the mean of these two models to determine the most plausible pressure-density function (Supplementary Fig. 2). The ice layer is assumed to be 600-m thick and to have a density of 900 kg m⁻³. The density in the uppermost kilometre of the crust has considerable uncertainty. We use a mean value of 2300 kg m⁻³ for the rocks at the ice-bedrock interface and a strong gradient down to 4 km depth below sea level (density 2950 kg m⁻³, pressure of ~130 MPa, see Supplementary Fig. 2). The lower panel in Supplementary Fig. 2 shows density as a function of depth. Depth values refer to Bárðarbunga, where zero is sea level.

**Volatiles and magma properties versus depth.** Following Kilbride et al.¹⁰, we adopted the D-Compress thermodynamic model of volatile saturation and partitioning⁴⁰ to predict the composition of a gas phase in equilibrium with the 2014–15 Holuhraun melt as a function of volatile content, melt composition, temperature, pressure and oxidation state as the magma rose from the magma body at ~7.5 km depth b.s.l. (230 MPa) to the surface (based on the crustal density model; Supplementary Fig. 2). On the basis of these output parameters, we calculated magma compressibility considering the exsolved gas phase following $(1/\rho_{\text{magma}})(\partial\rho_{\text{magma}}/\partial p)$ where $\rho$ is density and $p$ is pressure. Assuming an oxygen fugacity buffer of NNO+0.5, a temperature of 1170 °C and pressure of 230 MPa (ref. ³⁵), our modelling demonstrates that CO₂ will continue to be effectively lost from melt (Fig. 6a), reaching a near-zero (~0.001 wt.%) value in erupted and fully degassed melt at 0.1 MPa (1 bar), consistent with observations from the Holuhraun glass⁴¹. However, once vapour saturation is reached, sulfur (we only plot melt-dissolved H₂S in Fig. 6a) exsolves continuously, with a major loss restricted to pressures below 10 MPa. Moreover, water remains largely dissolved in the melt until pressures of about 5 MPa (<200 m). Rapid and efficient exsolution and resulting bubble
growth is therefore expected to occur once the melt has reached such low pressures. To quantitatively estimate the effects on magma buoyancy at mid-crustal levels, it is important to consider the resulting change in magma density. At 230 MPa, D-Compress calculates a bulk density of the magma of 2694 kg/m³, almost identical to the value calculated for melt (2698 kg/m³) at the same depth. These values are somewhat lower than the density inferred from the whole-rock composition of the lava flow described above (2722 kg/m³). D-Compress is used to calculate the density (Fig. 6b) and compressibility (Fig. 6c) changes of the bulk magma as a function of depth.

**Threshold for basal inflow rate.** Magma may accumulate in a deep magma body without inducing its rupture if the mean magma basal inflow rate, $\phi$, is small enough such that the volume of magma entering the body over a time period equal to the relaxation time ($\Delta V = \phi\tau$) does not induce pressure increase beyond the rock tensile strength. Using Eq. (4), we can thus derive an upper threshold value for the mean basal inflow:

$$ \phi < \left( 2 + \frac{4\pi b \mu}{3 a k} \right) \frac{a^3 T}{\tau \mu}. \quad (6) $$

Taking $a = 4.3$ km, $b = 0.4$ km, $T = 2.5$ MPa, $\mu = k = 10$ GPa, and the relaxation time as 10 years, the threshold value for the mean basal flow would be about 5 million m³/yr. This value is smaller than the inferred rate of basal inflow in the nearby Grímsvötn volcano<sup>73,74</sup> in the period from 2004 to 2012, consistent with the higher frequency of eruption observed at Grímsvötn<sup>75</sup>. Considering this threshold value for the mean basal flow, ~6000 years would be required to accumulate the 30 km³ of magma in an ellipsoid of dimensions suggested in this example, which is an approximate estimate of the maximum size of the magma body that fed the 2014–2015 activity under Bárðarbunga.

## Data availability

The source data underlying Figs. 1, 2b, 4, 6, and 8 and Supplementary Figs. 1 and 2 are provided in a Source Data file.

Received: 15 October 2019; Accepted: 8 April 2020;

## References

1. Thordarson, T. & Larsen, G. Volcanism in Iceland in historical time: volcano types, eruption styles and eruptive history. *J. Geodyn.* **43**, 118–152 (2007).

2. Pyle, D. M. *The Encyclopedia of Volcanoes*. p. 263–269 (Elsevier: MA, USA, 2000).

3. Halldorsson, S. A. et al. Isotopic-heterogeneity of the Thjorsa lava - Implications for mantle sources and crustal processes within the Eastern Rift Zone, Iceland. *Chem. Geol.* **255**, 305–316 (2008).

4. Gudmundsson, M. T. et al. Gradual caldera collapse at Bárðarbunga volcano, Iceland, regulated by lateral magma outflow. *Science* **353**, 6296 (2016).

5. Sigmundsson, F. et al. Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland. *Nature* **517**, 191–195 (2015).

6. Neave, D. A., Passmore, E., Maclennan, J., Fitton, G. & Thordarson, T. Crystal-melt relationships and the record of deep mixing and crystallization in the AD 1783 Laki eruption, Iceland. *J. Petrol.* **54**, 1661–1690 (2013).

7. Hartley, M. E., Bali, E., Maclennan, J., Neave, D. A. & Halldórsson, S. A. Melt inclusion constraints on petrogenesis of the 2014–2015 Holuhraun eruption, Iceland. *Contrib. Mineral. Petrol.* **173**, 10 (2018).

8. Heimisson, E. R., Segall, P. Physically consistent modeling of dike induced deformation and seismicity: application to the 2014 Bárðarbunga dike, Iceland, *J. Geophys. Res.* https://doi.org/10.1029/2019JB018141 (2020).

9. Hudson, T. S. et al. Deep crustal melt plumbing of Bárðarbunga volcano, Iceland. *Geophys. Res. Lett.* **44**, 8785–8794 (2017).

10. Kilbride, B. M., Edmonds, M. & Biggs, J. Observing eruptions of gas-rich compressible magmas from space. *Nat. Commun.* **7**, 1–8 (2016).

11. Pinel, V. & Jaupart, C. Magma chamber behavior beneath a volcanic edifice. *J. Geophys. Res.* **108**, 2072 (2003).

12. Pinel, V., Jaupart, C. & Albino, F. On the relationship between cycles of eruptive activity and growth of a volcanic edifice. *J. Volcanol. Geotherm. Res.* **194**, 150–164 (2010).

13. McLeod, P. The role of magma buoyancy in caldera-forming eruptions. *Geophys. Res. Lett.* **26**, 2299–2302 (1999).

14. Pinel, V. & Jaupart, C. Caldera formation by magma withdrawal from a reservoir beneath a volcanic edifice. *Earth Planet. Sci. Lett.* **230**, 273–287 (2005).

15. Newman, A. V., Dixon, T. H. & Gourmelen, N. A four-dimensional viscoelastic deformation model for Long Valley Caldera, California, between 1995 and 2000. *J. Volcanol. Geotherm. Res.* **150**, 244–269 (2006).

16. Segall, P. Repressurization following eruption from a magma chamber with a viscoelastic aureole. *J. Geophys. Res. Solid Earth* **121**, 8501–8522 (2016).

17. Kelemen, P. B., Koga, K. & Shimizu, N. Geochemistry of gabbro sills in the crust/mantle transition zone of the Oman ophiolite: implications for the origin of the oceanic lower crust. *Earth Planet. Sci. Lett.* **146**, 475–488 (1997).

18. Maclennan, J. Concurrent mixing and cooling of melts under Iceland. *J. Petrol.* **49**, 1931–1953 (2008).

19. Yamasaki, T., Kobayashi, T., Wright, T. J. & Fukahata, Y. Viscoelastic crustal deformation by magmatic intrusion: a case study in the Kutcharo caldera, eastern Hokkaido, Japan. *J. Volcanol. Geotherm. Res.* **349**, 128–145 (2018).

20. Auriac, A. et al. Iceland rising: Solid Earth response to ice retreat inferred from satellite radar interferometry and viscoelastic modeling. *J. Geophys. Res. Solid Earth* **118**, 1331–1344 (2013).

21. Caricchi, L., Annen, C., Blundy, J., Simpson, G. & Pinel, V. Frequency and magnitude of volcanic eruptions controlled by magma injection and buoyancy. *Nat. Geosci.* https://doi.org/10.1038/NGE02041 (2014).

22. Dzurisin, D. Volcano Deformation, New Geodetic Monitoring Techniques. (Springer-Praxis Books in Geophysical Sciences, London, UK, 2007).

23. Bonafede, M. & Ferrari, C. Analytical models of deformation and residual gravity changes due to a Mogi source in a viscoelastic medium. *Tectonophysics* **471**, 4–13 (2009).

24. Grapenthin, R., Sigmundsson, F., Geirsson, H., Árnadóttir, T. & Pinel, V. Icelandic rhythmics: annual modulation of land elevation and plate spreading by snow load. *Geophys. Res. Lett.* **33**, L24305 (2006).

25. Rivalta, E., Taisne, B., Bunger, A. P. & Katz, R. F. A review of mechanical models of dike propagation: Schools of thought, results and future directions. *Tectonophysics* **638**, 1–42 (2015).

26. Jónsson, S. Tensile rock mass strength estimated using InSAR. *Geophys. Res. Lett.* **39**, L21305 (2012).

27. Haimson, B. C. & Rummel, F. Hydrofracturing stress measurements in the Iceland Research Drilling Project drill hole at Reydarfjordur, Iceland. *J. Geophys. Res. Solid Earth* **87**, 6631–6649 (1982).

28. Benson, P. M. et al. Laboratory simulations of tensile fracture development in a volcanic conduit via cyclic magma pressurisation. *Earth Planet. Sci. Lett.* **349–350**, 231–239 (2012).

29. Gerbault, M. in *Faulting, Fracturing and Igneous Intrusion in the Earth’s Crust*, (ed. Healy, D.) Vol. 367, p. 111–130 (Geological Society Special Publications, 2012).

30. Currenti, G. & Williams, C. A. Numerical modeling of deformation and stress fields around a magma chamber: constraints on failure conditions and rheology. *Phys. Earth Planet. Inter.* **226**, 14–27 (2014).

31. Heimisson, E. R., Hooper, A. & Sigmundsson, F. Forecasting the path of a laterally propagating dike. *J. Geophys. Res. Solid Earth* **120**, 8774–8792 (2015).

32. Roper, S. M. & Lister, J. R. Buoyancy-driven crack propagation from an over-pressured source. *J. Fluid Mech.* **536**, 79–98 (2005).

33. Jaupart, C. *The Encyclopedia of Volcanoes*. p. 237–245 (Elsevier, MA, USA, 2000).

34. Cashman, K., Sparks, S. J. & Blundy, J. D. Vertically extensive and unstable magmatic systems: a unified view of igneous processes. *Science* **355**, eaag3055 (2017).

35. Halldórsson, S. A. et al. Petrology and geochemistry of the 2014–2015 Holuhraun eruption, central Iceland: compositional and mineralogical characteristics, temporal variability and magma storage. *Contrib. Mineral. Petrol.* **173**, 64 (2018).

36. Bruce, P. M. & Huppert, H. E. Thermal control of basaltic fissure eruptions. *Nature* **342**, 665–667 (1989).

37. Sigmundsson, F. New insights into magma plumbing along rift systems from detailed observations of eruptive behavior at Axial volcano. *Geophys. Res. Lett.* **43**, 12423–12427 (2016).

38. Sigmundsson, F. Calderas collapse as magma flows into rifts. *Science* **366**, 1200–1201 (2019).

39. Edmonds, M. & Woods, A. W. Exsolved volatiles in magma reservoirs. *J. Volcanol. Geotherm. Res.* **368**, 13–30 (2018).

40. Burgisser, A., Alletti, M. & Scaillet, B. Simulating the behavior of volatiles belonging to the C–O–H–S system in silicate melts under magmatic conditions with the software D-Compress. *Comput. Geosci.* **79**, 1–14 (2015).

41. Bali, E., Hartley, M. E., Halldórsson, S. A., Gudfinnsson, G. H. & Jakobsson, S. Melt inclusion constraints on volatile systematics and degassing history of the 2014–2015 Holuhraun eruption, Iceland. *Contrib. Mineral. Petrol.* **173**, 9 (2018).

42. Gudmundsson, M. T. & Högnadóttir, T. Volcanic systems and calderas in the Vatnajökull region, central Iceland: constraints on crustal structure from gravity data. *J. Geodyn.* **43**, 153–169 (2007).

43. Einarsson, P. et al. Center of the Iceland hotspot experiences volcanic unrest. *Eos* **78**, 374–375 (1997).
44. Parks, M. M. et al. Evolution of deformation and stress changes during the caldera collapse and dyking at Bårdarbunga, 2014–2015: implication for triggering of seismicity at nearby Tungnafellsjökull volcano. *Earth Planet. Sci. Lett.* **462**, 212–223 (2017).

45. Hartley, M. & Maclennan, J. Magmatic densities control erupted volumes in Icelandic volcanic systems. *Front. Earth Sci.* **6**, 29 (2018).

46. Hauri, E. H. et al. CO₂ content beneath northern Iceland and the variability of mantle carbon. *Geology* **46**, 55–58 (2018).

47. Fialko, Y. & Simons, M. Evidence for ongoing inflation of the Socorro magma body, New Mexico, from Interferometric synthetic aperture radar imaging. *Geophys. Res. Lett.* **28**, 3549–3552 (2001).

48. Pearse, J. & Fialko, Y. Mechanics of active magmatic intraplate in the Rio Grande Rift near Socorro, New Mexico. *J. Geophys. Res.* **115**, B07413 (2010).

49. Ake, J. P. & Sanford, A. R. New evidence for the existence and internal structure of a thin layer of magma at mid-crustal depths near Socorro, New Mexico. *Bull. Seismol. Soc. Am.* **78**, 1335–1359 (1988).

50. Anderson, K. R. et al. Magma reservoir failure and the onset of caldera collapse at Kilauea Volcano in 2018. *Science* **366**, eaaz1822 (2019).

51. Feuillet, N. et al. Birth of a large volcano offshore Mayotte through lithosphere-scale rifting. Fall Meeting 2019, Abstract #V52D-01 (American Geophysical Union, 2019).

52. Bachelery, P. et al. Petrological and geochemical characterization of the lava from the 2018–2019 Mayotte eruption: first results. Fall Meeting 2020, Abstract V52D-06 (American Geophysical Union, 2019).

53. Cesa, S. et al. Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation. *Nat. Geosci.* **13**, 87–93 (2020).

54. Karlstrom, L., Dufek, J. & Manga, M. Organization of volcanic plumbing through magmatic lensing by magma chambers and volcanic loads. *J. Geophys. Res.* **114**, B10204 (2009).

55. Karlstrom, L., Dufek, J. & Manga, M. Magma chamber stability in arc and continental crust. *J. Volcanol. Geotherm. Res.* **190**, 249–270 (2010).

56. Degruyter, W. & Huber, C. A model for eruption frequency of upper crustal silicic magma chambers. *Earth Planet. Sci. Lett.* **403**, 117–130 (2014).

57. Hadamard, J. Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide viseueux. C. R. Acad. Sci. **152**, 1735–1738 (1911).

58. Weinberg, R. F. & Podladchikov, Y. Diapiric ascent of magma through power law crust and mantle. *J. Geophys. Res.* **99**, 9543–9559 (1994).

59. Seropian, G., Rust, A. C. & Sparks, R. S. J. The gravitational stability of lenses in magma mushes: Confined Rayleigh-Taylor instabilities. *J. Geophys. Res. Solid Earth* **123**, 3593–3607 (2018).

60. Lister, J. R. Buoyancy-driven fluid fracture: similarity solutions for the horizontal and vertical propagation of fluid-filled cracks. *J. Fluid Mech.* **217**, 213–239 (1990).

61. Lange, A. E., Nielsen, R. L., Tepley, F. J. III & Kent, A. J. R. The petrogenesis of plagioclase-phyric basalts at mid-ocean ridges. *Geochem. Geophys. Geosyst.* **14**, 3282–3296 (2013).

62. Oskarsson, B. V., Andersen, C. B., Riishuus, M. S., Sorensen, E. V. & Tegner, C. The mode of emplacement of neogene flood basalts in eastern Iceland: the plagioclase ultraphyric basalts in the Grænavatn group. *J. Volcanol. Geotherm. Res.* **322**, 26–50 (2017).

63. Neave, D. A., Buisman, I. & Maclennan, J. Continuous mush disaggregation during the long-lasting Laki fissure erupton, Iceland. *Am. Mineral.* **102**, 2007–2021 (2017).

64. Hartley, M. E., Morgan, D. J., Maclennan, J., Edmonds, M. & Thordarson, T. Tracking timescales of short-term precursors to large basaltic fissure eruptions through Fe–Mg diffusion in olivine. *Earth Planet. Sci. Lett.* **439**, 58–70 (2016).

65. Rivalta, E. & Segall, S. Magma compressibility and the missing source for some dike intrusions. *Geophys. Res. Lett.* **35**, L04306 (2008).

66. Pollard, D. D. Derivation and evaluation of a mechanical model for sheet intrusions. *Tectonophysics* **19**, 233–269 (1973).

67. Danyushevsky, L. V. & Plechov, P. Petrolog3: integrated software for modeling crystallization processes. *Geochem. Geophys. Geosyst.* **12**, 1–32 (2011).

68. Lange, R. A. & Carmichael, I. S. E. Densities of Na₂O-K₂O-CaO-MgO-FeO-Fe₂O₃-Al₂O₃- TiO₂-SiO₂ liquids: new measurements and derived partial molar properties. *Geochim. Cosmochim. Acta* **51**, 2931–2946 (1987).

69. Brandsdóttir, B. & Menke, W. H. The seismic structure of Iceland. *Jökull* **58**, 17–34 (2008).

70. Flóvenz, Ó. G. & Gunnarsson, K. Seismic crustal structure in Iceland and surrounding area. *Tectonophysics* **189**, 1–17 (1991).

71. Carlson, R. L. & Herrick, C. N. Densities and porosities in the oceanic crust and their variations with depth and age. *J. Geophys. Res.* **95**, 9153–9170 (1990).

72. Christensen, N. I. & Wilkens, R. H. Seismic properties, density, and composition of the Icelandic crust near Reyðarfjörður. *J. Geophys. Res.* **87**, 6389–6395 (1982).

73. Reverso, T. et al. A two-magma chamber model as a source of deformation at Grimsvötn Volcano, Iceland. *J. Geophys. Res. Solid Earth* **119**, 4666–4683 (2014).

74. Bato, M. G., Pinel, V., Yan, Y., Jouanne, F. & Vandemeulebrouck, J. Possible deep connection between volcanic systems evidenced by sequential assimilation of geodetic data. *Sci. Rep.* **8**, 11702 (2018).

75. Larsen, G., Gudmundsson, M. T. & Bjørnsson, H. Eight centuries of periodic volcanism at the center of the Iceland hotspot revealed by glacier tephrostratigraphy. *Geology* **26**, 943–946 (1998).

76. Einarsson, P. Earthquakes and present-day tectonism in Iceland. *Tectonophysics* **189**, 261–279 (1991).

77. Porter, C. et al. ArcticDEM, Harvard Dataverse, V1, https://doi.org/10.7910/DVN/OHHUKH (2018).

## Acknowledgements

The research presented here has benefitted from extended visits of FS during a sabbatical term to, and discussion with scientists at, the University of Leeds, ISTERre University of Savoie Mont-Blanc, USGS Cascades Volcano Observatory, and Geological Survey of Japan. We acknowledge reviews by Philip Benson and Luca Caricchi that helped to significantly improve the paper, as well as reviews of an early version of the paper by two anonymous reviewers. Financial support from the H2020 project EUROVOLC funded by the European Commission is acknowledged (grant number 731070). F.S. acknowledges support from the University of Iceland Research Fund, and R.G. acknowledges partial support through NSF grant EAR-1464546. Fissure swarms, central volcanoes and caldera outlines shown in Fig. 1 are reproduced from publications referred to (refs. 42,76) with permissions from Elsevier, and we acknowledge the use of ArticDEM (ref. 77) to plot surface and ice topography shown in Fig. 1. COMET is the NERC Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics, a partnership between UK Universities and the British Geological Survey.

## Author contributions

F.S. led the development of the ideas and modelling framework presented in this paper, with the participation of V.P., R.G., A.H., S.A.H., P.E., E.R.H., M.T.G., T.W. and T.Y.. Numerical modelling with the COMSOL software was carried out by V.P., and S.A.H. did the D-Compress modelling. The crustal density model was made by M.T.G.; Collection and analyses of seismic and geodetic data was carried out by B.G.Ó., K.J., K.V., M.P., S.D. H.M.F., G.B.G., P.E., E.R.H., H.G., F.S., A.H., S.L. and V.D. All the authors contributed to evaluation of the modelling, discussion of the results and the writing of the paper.

## Competing interests

The authors declare no competing interests.

## Additional information

Supplementary information is available for this paper at https://doi.org/10.1038/s41467-020-16054-6.

Correspondence and requests for materials should be addressed to F.S.

Peer review information *Nature Communications* thanks Philip Benson and Luca Caricchi for their contribution to the peer review of this work. Peer reviewer reports are available.

Reprints and permission information is available at http://www.nature.com/reprints

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_11.png" /></figure>

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020