Modelling satellite-derived magma discharge to explain caldera collapse

On 29 August 2014, one of largest effusive eruptions occurring in historic times in Europe began ∼45 km northeast of the center of the Bárðarbunga volcanic system, Iceland (Sigmundsson et al., 2015; Fig. 1A). The Holuhraun eruption followed 15 d of sustained seismicity that accompanied the propagation of a 45-km-long segmented dike, initiated at 10–12 km beneath the caldera of Bárðarbunga volcano (Sigmundsson et al., 2015; Ágústsdóttir et al., 2016; Gudmundsson et al., 2016). The effusive activity persisted for ∼180 d and was accompanied by the slow collapse of the ice-covered summit caldera of Bárðarbunga which began a few days after the beginning of the seismicity (Gudmundsson et al., 2016). The collapse was coeval with a series of magnitude M 5 earthquakes whose distribution correlates with the margins of a subsiding piston ∼7 km in diameter (Riel et al., 2015). Caldera floor subsidence decelerated exponentially with time and, by the end of the eruption on 27 February 2015, had formed a bedrock depression ∼65 ± 3 m deep, ∼1.8 ± 0.2 km³ in volume (Gudmundsson et al., 2016).

Using moderate resolution imaging spectroradiometer (MODIS) thermal data (see the GSA Data Repository¹), we estimated a total volume of erupted lava (∼1.8 ± 0.9 km³) similar to that of the final caldera depression. We also show that lava effusion associated with the distal fissure eruption was characterized by an exponential decay in effusion rates until the end of January 2015 (Fig. 2). Notably, 1 mo before the end of the eruptive event, the trend in effusion rates deviated drastically from a typical exponential decay and declined more rapidly until the complete ceasing of surface activity. This overall pattern very closely reproduces the subsidence of the caldera floor measured by GPS (Gudmundsson et al., 2016) thus suggesting an inelastic link between lateral magma withdrawal and caldera collapse (Fig. 3).

Satellite-based thermal data have been routinely used to estimate lava discharge rates during effusive eruptions (Harris and Baloga, 2009, and references therein). Here we used the radiant flux provided by the Middle Infrared Observation of Volcanic Activity system (MIROVA; Coppola et al., 2016b), which allows quantification of the effusive activity (Fig. 2) by using MODIS infrared data (see the Data Repository).

The highest effusion rate (360 ± 180 m³ s^–1) was recorded on 11 September 2014 (Fig. 2) and likely reflected the “buffered” thermal response related to the initial emplacement of a 17-km-long lava flow that moved toward the northeast (Fig. DR1 in the Data Repository). Successively, two other main flow units (emplaced on the southern side of the above) produced two other thermally buffered pulses, with decreasing peaks, recorded on 1 October (270 ± 135 m³ s^–1) and 26 October (175 ± 87 m³ s^–1). The inferred effusion rate declined gradually throughout December and mid-January when it was 100 ± 50 m³ (Fig. 2). However, in early February 2015, the effusion rate began to decrease more rapidly to 10 ± 5 m³ s^–1 by mid-month, and <0.5 m³ s^–1 on 26 February 2015, 1 d before the eruption was over (Fig. 2). Based on the MODIS data, we calculate a total volume of erupted lava of 1.8 ± 0.9 km³ (Fig. 3A), overlapping with the 1.3–1.7 km³ estimated independently on the basis of topographic surveys (Gudmundsson et al., 2016).

Besides the large quantity of lava erupted, the most significant feature of this eruption was the overall exponential decay of the lava flux (R² = 0.81) recorded during the first 5 mo of activity (blackline in Fig. 2). This decay is also evident in the cumulative volume of lava discharged (R² = 0.99; Fig. 3A) and suggests that the effusion rate resulted from the systematic decrease of the pressure driving the magma flow within the dike.

Although different relationships are obtained for different cross-sections of the pipe (dike) or for turbulent flow (White, 1981), Equations 2 and 3 state that for larger R reached, the smaller the eruption rate: Q(t) = P_f(t)/R.

Two distinct models can explain this exponential decay by assuming a different origin for the excess of pressure, P_f, and by involving distinct capacities, C, in the magma plumbing system. These can be ascribed to elastic (e.g., Wadge, 1981) and inelastic (e.g., Ripepe et al., 2015) discharge dynamics.

The best-fit of the exponential decay (Equation 1) to the effusion rates derived by MODIS during the Holuhraun eruption gives, in its logarithmic form (Fig. 2), an initial flow rate, Q₀, of ∼242 (± 121) m³ s^–1 and allows us to calculate the characteristic relaxation time, τ, of 9.6 × 10⁶ s (∼111 d). Using these values in Equation 1, we modeled the lava discharge volumes, which show a very good fit (R² = 0.99) with the volumes measured by MODIS (Fig. 3A).

According to the gravity-driven inelastic model, we consider a cylindrical magma chamber located below Bárðarbunga with radius r = 3500 m, being consistent with the size of the collapsed caldera (Fig. 1B). By means of Equation 8, for a magma density ρ equal to 2650 kg m^–3, the capacity C of the chamber is 1.5 × 10³ m³ Pa^–1. This constrains the bulk flow resistance of the dike, R, to 6.4 × 10³ Pa•s m^–3 and permits an estimation of the range of the dike cross-section as a function of the magma viscosity. Thus, we assume that the viscosity of the Holuhraun magma may range between 10¹ and 10⁴ Pa•s, which is typical for basalts. Given a dike length L equal to 45 km, the radius a of the flow channel, compatible with the calculated resistance of the flow, is constrained between 4 and 20 m (Equation 3). Therefore, the volume of the feeding flow channel is constrained between 2 × 10⁶ and 56 × 10⁶ m³.

Modeling indicates that only 1.8 (± 0.9) × 10⁹ m³ of the total excess of magma volume in the reservoir [V_e(0) = 2.3 (± 1.1) × 10⁹ m³] were finally erupted (Fig. 3), thus suggesting that ∼0.5 (± 0.2) × 10⁹ m³ remained inside the reservoir. This model diverges from the measured effusive rate only after 27 January 2015, when the rate of magma erupted declined more rapidly until the end of the eruption, occurring on 27 February 2015 (Figs. 2 and 3A). By using a best linear fit analysis (see the Data Repository), we calculate that deviation from the model occurred after ∼151 d of activity, when the effusion rate dropped below the critical value Q_c of ∼50 (± 25) m³ s^–1 (Fig. 2). From Equation 2, we estimate (see the Data Repository) that this critical flow rate is related to an excess pressure, P_f, of ∼0.35 (± 0.17) × 10⁶ Pa, which can be ascribed to the critical minimum pressure (P_c = RQ_c) necessary to enable the lateral transport of magma over a length of 45 km.

The volume of lava erupted after 27 January 2015 (∼37 ± 18 × 10⁶ m³; see the Data Repository) falls into the range of values estimated for the feeding flow channel (2–56 × 10⁶ m³) and seems to suggest a relationship with the gradual drainage of the elastic-walled flow path (e.g., Bokhove et al., 2005). Alternatively, the sharp decrease in the effusion rate could be related to a simple thermal effect (freezing) that caused the gradual solidification and closure at the far end of the flow path once the excess pressure fell below a critical threshold (Fialko and Rubin, 1999).

The large effusive eruption of Holuhraun was associated with the 65 ± 3-m-deep subsidence of the floor of Bárðarbunga caldera (Gudmundsson et al., 2016). The clear correlation between effusive rates measured by satellite and the ground deformation measured by GPS (Fig. 3B) suggests a close link between these phenomena. The inelastic model well explains the contraction of the magma reservoir in terms of magmastatic load and allows calculation (by Equation 9) of the effect of the effusive trend, Q(t), measured by satellite, on the subsidence, h(t), of the Bárðarbunga caldera (measured on site by GPS and corrected for ice flow; Fig. 3B). The initial effusive rate of 242 (± 121) m³ s^–1 and the relaxation time, τ = RC, indicate (Equation 9) that the effusive process at Holuhraun is consistent with the release of an excess pressure of ∼1.6 (± 0. 8) × 10⁶ Pa (Equation 2) at the beginning of the eruption, being equivalent to a total magmastatic load of ∼64 ± 32 m (Fig. 3B). However, during the 180 d of effusive activity, only ∼49 m of the expected subsidence is reproduced by our modeling, resulting in a lack of ∼15 m of total vertical displacement. By taking into account that the caldera had already subsided by ∼15 m in the early stages of the eruption (Gudmundsson et al., 2016; Fig. 1C), we show that our model matches the vertical displacement of the caldera floor, being consistent with GPS data. Notably, the subsidence at the center of the caldera slowed down and deviated from the modeled exponential trend at the beginning of February 2015 (see the Data Repository), substantially when MODIS-derived effusion rates decelerated (Fig. 3B).

The volume of magma depleted from the reservoir during the Bárðarbunga-Holuhraun eruption was essentially compensated by an irreversible deformation of the volcanic edifice. This feature is not compatible with the elastic recovery of magma chamber walls and suggests an alternative, inelastic process to explain the exponential decay of the effusion rate and the collapse of the Bárðarbunga caldera. The gravity-driven model explains the contraction of the magma chamber in terms of simple magmastatic load changes and provides a new mechanism to explain the origin of excess magma pressure during effusive eruptions. This model was originally developed to explain effusive eruptions at Stromboli (Ripepe et al., 2015) and implies the presence of a magma column just above the eruptive vent. However, our analysis indicates that the gravity-driven dynamic is independent of the magma chamber’s depth because it is exclusively based on the excess pressure changes during the eruption.

According to Gudmundsson (2012), two main types of loading may produce excess pressure in the deep chamber until the point of rupture: (1) influx of magma within the chamber, and (2) external extension (e.g., rifting episodes) that reduces the minimum principal compressive stress (σ₃).

In the case of the Bárðarbunga-Holuhraun eruption, we suggest that the reduction of the minimum principal compressive stress, associated with the initial rifting event before the eruption onset (Gudmundsson et al., 2014; Sigmundsson et al., 2015), induced the initial excess pressure required by the gravity-driven model (P_f ∼1.6 MPa). This pressure is similar to the in situ tensile strength (2–4 MPa) of the Icelandic crust (Gudmundsson, 2012) and likely reflects the compressive stress drop (of σ₃) that accompanied the initial rifting episode.

Previous eruptions of Bárðarbunga could also have been driven by a similar mechanism. Gravity-driven dynamics could link, for example, the formation of the subglacial caldera of Bárðarbunga (∼700 m deep) with the largest (V_e ∼25 km³) fissure eruption ascribed to this volcanic system (i.e., the eruption of the Great Thjórsá Lava on the fissure swarm southwest of Vatnajökull; Larsen and Gudmundsson, 2014). In fact, using the characteristic capacity, C = 1.5 × 10³ m³ Pa^–1, found for Bárðarbunga, we calculate that during this large event, the pressure drop inside the reservoir reached ∼16 MPa (P_f= V_e/C), which corresponds to the lowering of the magma level by ∼650 m, in close agreement with the depth of the subglacial caldera.

Other volcanic eruptions exhibiting a decay in the effusive rates (e.g., Coppola et al., 2016a) associated with caldera subsidence (e.g., Geshi et al., 2002) could be related to the above-described mechanism and may be revisited in the light of gravity-driven processes. The model presented herein can also provide a new strategy to better interpret effusive trends. The deviation from exponential decay has been observed one month before the end of the eruption (Figs. 2 and 3) and could be used, in near real time, to eventually forecast the final stages of an effusive eruption.

We thank editor B. Murphy, L. Wilson and F. Sigmundsson for their valuable suggestions and comments on the manuscript. We acknowledge the NASA LANCE-MODIS data system for providing MODIS Near Real Time products.