Modelling satellite-derived magma discharge to explain caldera collapse

Many effusive eruptions are characterized by effusion rates that decay exponentially with time, a trend which is generally ascribed to elastic relaxation of a deep magma chamber. Thermal emissions, detected by satellite during the A.D. 2014–2015 BárðarbungaHoluhraun eruption (Iceland), indicate that the volume of the erupted magma and effusion rates followed an overall exponential trend that fits the observed major subsidence of the Bárðarbunga caldera floor. This trend continued until a critical flow rate was reached. Hence, the subsidence slowed down and the eruption rapidly ceased, reflecting the ultimate closure of the magma path. We present a model of inelastic magma withdrawal that very closely reproduces all the observed phenomena and provides new insights into the caldera collapses and the driving pressure of other effusive eruptions. INTRODUCTION 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 km3 in volume (Gudmundsson et al., 2016). Using moderate resolution imaging spectroradiometer (MODIS) thermal data (see the GSA Data Repository1), we estimated a total volume of erupted lava (~1.8 ± 0.9 km3) 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 1 GSA Data Repository item 2017165, materials and methods and Figures DR1– DR3, is available online at http://www.geosociety.org/datarepository/2017/ or on request from editing@geosociety.org. *E-mail: diego.coppola@unito.it GEOLOGY, June 2017; v. 45; no. 6; p. 523–526 | Data Repository item 2017165 | doi:10.1130/G38866.1 | Published online 16 March 2017 © 2017 The Authors. Gold Open Access: This paper is published under the terms of the CC-BY license. B A C Figure 1. A: Map of Bárðarbunga-Holuhraun volcanic system (Iceland) showing subsiding caldera, lateral dike, and lava flow field associated with the A.D. 2014–2015 eruption. B: Sketch of modeled plumbing system under Bárðarbunga caldera. Cyan, gray, and red colors refer to ice, rock, and magma levels, respectively. C: Subsidence at Bárðarbunga caldera as measured on 5 September 2014 by aerial radar profiling (modified after Sigmundsson et al., 2015). Black dashed circle outlines extension of modeled collapsing piston with radius r = 3500 m.


INTRODUCTION
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 3 in volume (Gudmundsson et al., 2016).
Using moderate resolution imaging spectroradiometer (MODIS) thermal data (see the GSA Data Repository 1 ), we estimated a total volume of erupted lava (~1.8 ± 0.9 km 3 ) 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-DERIVED EFFUSION RATES DURING THE HOLUHRAUN 2014-2015 ERUPTION
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 3 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 3 s -1 ) and 26 October (175 ± 87 m 3 s -1 ). The inferred effusion rate declined gradually throughout December and mid-January when it was 100 ± 50 m 3 (Fig. 2). However, in early February 2015, the effusion rate began to decrease more rapidly to 10 ± 5 m 3 s -1 by mid-month, and <0.5 m 3 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 3 (Fig. 3A), overlapping with the 1.3-1.7 km 3 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 2 = 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 2 = 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.

MODELS FOR EXPONENTIAL MAGMA DISCHARGE RATE
The exponential decay of effusion rates has been observed in many eruptions, and it is generally ascribed to the release, during an eruption, of elastic energy stored in the magma chamber (Machado, 1974;Wadge, 1981). This trend is generally described by: where Q(t) is the effusion rate, Q 0 is the initial flow rate, t is time, and τ is the decay time constant, which depends (τ = RC) on the "hydrodynamic resistance", R, of the magma flow path and on the "capacity", C, of the magma reservoir (Aki and Ferrazzini, 2001). In a pipe (dike or sill) of length L with equivalent radius a and for a given excess pressure P f (t) in the magma chamber, the effusion rate is controlled by the resistance R of the Poiseuille laminar flow with magma viscosity η (Fig. 1B): . ( 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: Effusive eruptions can be thus seen as the "strategy" of the volcano to reduce the excess pressure P f within the reservoir by draining an excess volume, V e , of magma. In turn, the excess pressure in the reservoir is  considered to be dependent on the capacity, C, of the volcanic reservoir to induce pressure on the fluid, P f (t) = V e (t)/C, for a given volume V e (t) > 0. The effusion rate can be thus written (Equations 2 and 3) as function of the excess volume of magma discharged by the reservoir: . (4) If the effusive rate Q(t) = -dV e /dt is the outflow discharge rate, Equation 4 can be solved in terms of the excess volume (e.g., Wadge, 1981): , which shows how the excess of volume is drained out of the magmatic system following an exponential decay. 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.

Elastic Relaxation Model
Storage of magma within the crust involves the transfer of elastic strain energy to reservoir rocks and to the magma itself, most of which is released again as magma leaves the system (Wadge, 1981). The exponential decrease of the effusion rate thus results from the contraction of an elastic magma chamber with capacity: , being controlled by the total volume of the magma chamber (V m ) and by the bulk modulus K m of magma (e.g., Wadge, 1981). The outflow of magma is then recovered by progressive decompression of the magmatic system (i.e., Hreinsdóttir et al., 2014), with changes in the reservoir volume much lower than the volume of the erupted lava (Johnson et al., 2000). This process (combined inflation and deflation) does not involve a permanent deformation of the volcano edifice (Parfitt and Wilson, 2008) and seems incompatible with the piston-like collapse of a caldera inside the magma chamber. Therefore, the effusive trend observed during the Holuhraun eruption and the associated collapse of the Bárðarbunga caldera become rather difficult to explain in terms of progressive elastic contraction of the deep magma chamber.

Inelastic Gravity Model
The inelastic gravity-driven model, recently proposed by Ripepe et al. (2015), explains the exponential decay of effusion rates as driven by the decrease of magmastatic pressure (without the elastic recovery of the magma chamber). The model assumes that the pressure, P f , controlling the effusion rate (Equation 2) is equivalent to the magmastatic excess pressure: , induced by an excess of magma volume, V e (t) = πr 2 h(t), stored in a cylindrical-like reservoir with height (h) and equivalent radius (r), located above the outlet to the channel feeding the eruptive vent (Fig. 1b). In this model the capacity of the reservoir is represented by a pure static coefficient: , which depends on gravity (g), magma density (ρ) and on the geometry of the magma reservoir. The gravity-driven model explains at once the exponential decay of effusion rate and the lowering of the magma level inside a cylindrical reservoir: , which makes this model suitable to explain permanent deformations associated with effusive dynamics. This model is somewhat similar to the one that has been recently proposed to explain the subsidence of the Bárðarbunga caldera (Gudmundsson et al., 2016) but assumes that the flow path is fully horizontal. In this case the pressure needed for vertical flow over the eruption site is provided by the weight of the collapsing piston so that the excess pressure, P f , is the unique driving force governing the gravity-driven dynamic.

MODELING THE EFFUSIVE TREND
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 0 , of ~242 (± 121) m 3 s -1 and allows us to calculate the characteristic relaxation time, τ, of 9.6 × 10 6 s (~111 d).
Using these values in Equation 1, we modeled the lava discharge volumes, which show a very good fit (R 2 = 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 3 m 3 Pa -1 . This constrains the bulk flow resistance of the dike, R, to 6.4 × 10 3 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 1 and 10 4 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 6 and 56 × 10 6 m 3 .
Modeling indicates that only 1.8 (± 0.9) × 10 9 m 3 of the total excess of magma volume in the reservoir [V e (0) = 2.3 (± 1.1) × 10 9 m 3 ] were finally erupted (Fig. 3), thus suggesting that ~0.5 (± 0.2) × 10 9 m 3 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 3 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 6 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 6 m 3 ; see the Data Repository) falls into the range of values estimated for the feeding flow channel (2-56 × 10 6 m 3 ) 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).

DISCHARGE MODELS AND CALDERA COLLAPSE
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 3 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 6 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).

DISCUSSION AND CONCLUSIONS
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 (σ 3 ).
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 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 σ 3 ) 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 3 ) 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 3 m 3 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.