Journal of Geophysical Research: Solid Earth Insights Into the Stress Field Around Bárðarbunga Volcano From the 2014/2015 Holuhraun Rifting Event

The two weeklong rifting event at Bár ð arbunga volcano in 2014 led to the Holuhraun eruption, which produced 1.5 km 3 of lava and was the largest in Iceland in over 200 years. Predicting when and where an intrusion will lead to eruption requires detailed knowledge of the underlying stress ﬁeld. Previous studies have explained the dike propagation path with a model that includes a tectonically induced stress ﬁeld set up by a uniform amount of plate spreading across a straight rift axis. Here we test this hypothesis by modeling the tractions acting on the dike walls, constrained by data from Global Navigation Satellite System and Interferometric Synthetic Aperture Radar. Our results show that the majority of the opening and shearing in the ﬁnal two dike segments is due to stresses built up by plate spreading since the last eruption at Holuhraun, as expected, but that the tectonically induced stress magnitude must be much lower to explain the movement of the dike walls further south. This result implies that most of the tectonically induced stress beneath the ice cap has been released, presumably due to intrusions associated with the Bár ð arbunga volcanic system and the nearby Grímsvötn volcanic system, which have not been detected due to their subglacial nature. Modeling of the 2014 Bár ð arbunga rifting event therefore not only yields insights into the event but also provides a window into undetected volcanic activity in the past.


Introduction
The subglacial Bárðarbunga Volcano in Iceland lies beneath the Vatnajökull ice cap (Figure 1). The associated volcanic system consists of a central volcano, a caldera, and fissure swarms, which extend to the south-southwest and north-northeast for a total length of 170 km. Bárðarbunga has had 23 confirmed eruptions since the settlement of Iceland 1,100 years ago (Thordarson & Larsen, 2007). It also produced the largest Holocene lava flow in the world by both volume and area (25 km 3 and 967 km 2 , respectively), the Þjórsá lava field, which erupted around 8,600 years before present (Hjartarson, 2003).
The 2014 Holuhraun rifting episode began with a seismic swarm on 16 August . Seismic activity and Global Navigation Satellite System (GNSS) observations showed a dike moving initially radially away from the caldera toward the east-southeast and then turning toward the north-northeast ; see seismicity in Figure 1. Dike progress continued for 20 km until 19 August, after which propagation stopped for 80 hr. On 23 August, the dike briefly turned hard left to propagate in a north-northwesterly direction. The final change of direction left the dike to propagate in a north-northeasterly direction, in which it continued until 27 August, with the tip of the dike located approximately 10-km north of the Vatnajökull ice cap. The diking event was accompanied by significant subsidence of the caldera of tens of meters in the first few weeks after the onset of the eruption, showing subsidence rates exceeding 50 cm/day (Gudmundsson et al., 2016;Riel et al., 2015).
The first of two eruptive events commenced on 29 August 2014, a minor event lasting only 4 hr. On 31 August 2014, a second eruption started from the same fissure, which continued for 6 months, until 27 February 2015 (Gudmundsson et al., 2016). The eruptive fissure is located in the older Holuhraun lava field, which is thought to have been emplaced sometime between 1794 and 1864 (Hartley & Thordarson, 2013). From this fissure, a lava field developed that covered 85 km 2 by the end of the 2014/2015 eruption, with a volume of 1.5 km 3 of lava, making it the largest eruption in Iceland in over 200 years (Schmidt et al., 2015).  Sigmundsson et al. (2015). The color of each dot represents the day of the earthquake. and thus the processes behind them. Is it magma pressure that causes the dike to open? Or is the opening the result of external stresses caused by tectonics or topography? Mechanical models, which solve for the stresses on the dike walls, are required to answer these kind of questions. Finite element modeling could be used for this, but these methods require a large volume of the Earth to be modeled, much larger than the area of interest, making it computationally inefficient. Boundary element modeling (Cayol & Cornet, 1997) on the other hand only considers the relevant boundaries, like dike or magma chamber walls, vastly reducing the problem size, resulting in an efficient method that can be used in Markov Chain Monte Carlo inversions (Hooper et al., 2011).
Here we use stress-driven boundary element modeling to constrain GNSS and Interferometric Synthetic Aperture Radar (InSAR) measurements of the dike propagation and the early stages of the eruptive event. We model several time steps, starting the day before the onset of the event on 16 August 2014 until 4 September, when the dike had finished its migration northward and the fissure eruption had been ongoing for several days. Our modeling investigates how the internal magma overpressure and the regional stress field interacted to create the opening and shearing pattern of the dike, shedding light on the dominant process and the stress field at the time of the eruption.

Deformation Observations
The deformation associated with the dike propagation was observed using a network of continuous and campaign GNSS stations , which we supplement with four InSAR interferograms. We used the GNSS data presented in Sigmundsson et al. (2015), which comprises 24-hr solutions of 31 GNSS stations processed using the GAMIT/GLOBK software (Herring et al., 2010) version 10.4.
The InSAR data are composed of two Cosmo-SkyMed interferograms, namely, descending track 2631 (13 August 2014 to 29 August 2014) and ascending track 2631 (30 July 2014 to 1 September 2014) and two TerraSAR-X interferograms, descending track 140 (13 August 2014 to 4 September 2014) and ascending track 147 (26 July 2012 to 4 September 2014). The descending Cosmo-SkyMed track and the ascending TerraSAR-X track were also used to constrain the models presented in Sigmundsson et al. (2015), but the other two tracks were not. All interferograms were coregistered and interfered using the Doris software package (Kampes, 1999). For Cosmo-SkyMed descending track 2631, a large time series of interferograms was available, allowing us to estimate the coherence using the RapidSAR method (Spaans & Hooper, 2016). For the other three interferograms, we used the boxcar method (Just & Bamler, 1994) with a window size of 11 by 11 to estimate the coherence. We estimated and removed the spatially correlated part of the phase from all four interferograms prior to coherence estimation to avoid the high-frequency deformation fringes biasing the estimation (Spaans & Hooper, 2016). The spatially correlated phase estimation was achieved using multilook filtering. The wrapped phase values of points with sufficient coherence were unwrapped using the Snaphu software (Chen & Zebker, 2001). To reduce the amount of data for modeling, we applied adaptive quadtree resampling to downsample the unwrapped interferograms (Decriem et al., 2010;Jónsson et al., 2002), which allows arbitrarily shaped "quads," depending on the availability of points. The left column in Figure 2 shows the unwrapped, downsampled interferograms.
The interferograms each cover slightly different periods. The master date of all interferograms is before the onset of the event on 16 August, and we assume that there is no significant deformation in the area between the master date of each interferogram and the onset of the unrest. The slave dates of the interferograms, however, range from 29 August to 4 September. As significant deformation associated with dike emplacement continued until 4 September , we cannot treat all interferograms as if they cover the same deformation. To address this, we define five overlapping time intervals, where each subsequent time period is between 3 and 5 days longer than the previous one. The time periods all start on 15 August, and they end on 19, 24, and 29 August and 1 and 4 September, respectively. The black arrows in Figure 3 show the available GNSS deformation vectors in the vicinity of the eruptive site for each of the five intervals. A handful  Figure 1 fall outside of the area covered in Figure 3. These stations are, however, included to constrain the modeling described below, and all show no significant movements. The slave dates of the interferograms coincide with three of the time intervals. Thus, we can constrain models of the deformation in the first two intervals by GNSS measurements alone and the last three periods by a combination of GNSS and InSAR measurements.

Boundary Element Modeling
We used a boundary element approach to model the InSAR and GNSS observations. Our model is based on the method described in Hooper et al. (2011). Given the stresses acting on the boundaries, the displacements of each boundary patch are inverted for using a least squares approach, where the design matrix is defined  Sigmundsson et al. (2015). The circles indicate the 95% confidence region. Also displayed are the maximum a posteriori probability model predictions for the uniform (red arrows) and the depth difference models (green arrows). The red trace shows the path of the dike used in the modeling. by stress mapping functions. The stress mapping functions are the stresses at every patch in the three principle directions (normal, along strike, and along dip) to unit displacement at each patch in turn, for every principle direction. We calculate the stress mapping functions using rectangular dislocations (Okada, 1992), assuming a flat Earth geometry, a homogeneous elastic half-space, and infinitesimal strain. With the design matrix and the stresses at each patch defined, the displacements of the dike patches are inverted for, and the resulting surface displacements can be calculated.
We assume four dike segments, as defined in the modeling of Sigmundsson et al. (2015). We investigate two types of stresses acting on the dike walls. First, there is the magma overpressure in the dike, which is the difference between the magma pressure and the mean pressure. Second, there is the tectonically induced stress resulting from plate spreading. The tectonically induced stress is dependent on the distance from the central rift axis. The central axis of plate spreading has been found to go through the center of the Askja caldera (Sturkell & Sigmundsson, 2000), just north of the fissure location, and strikes at an angle of approximately 15 ∘ 10.1002/2017JB015274 . This would place the northernmost segments of the dike closest to the rift, increasing the influence of the plate spreading on the stress field. To model plate spreading, we assume a constant displacement rate, with opposite sign on either side of the rift axis, which is applied to the base of an elastic slab, over a zone of finite width. Heimisson et al. (2015) show that the displacements due to this model can be well approximated by an arctangent: where u(d) is the displacement due to plate spreading as a function of distance from the rift, U is the far-field plate separation, and D is a parameter related to the locking depth. The elastic slab model is favored over the buried dislocation model used extensively (e.g., Lafemina et al., 2005) due to its more realistic behavior closer to the rift axis. Heimisson et al. (2015) showed that the arctangent model closely resembles an elastic slab finite element model of the plate spreading, especially in terms of the strain and stress, which are what we use in this manuscript.
We use a value of 6,500 m for D, as was found for the same area by Heimisson et al. (2015). We take the derivative of equation (1) with respect to d to obtain the strain as a function of distance from the rift: (2) We then assume a value of 75 GPa for Young's modulus, as was estimated by Auriac et al. (2014) for this area, to calculate the resulting stress using Hooke's law. The direction of the tectonic stress is assumed to be perpendicular to the rift axis, which strikes at 15 ∘ . The magnitude of the tectonic stress is then dependent on the distance between the dike segment and the rift axis and the far-field plate separation parameter U, which we solve for in our inversion. Once the stress on the dike patch is calculated, we decompose it into a normal component and an along-strike component, and the normal component is added to the magma overpressure.
We account for the caldera subsidence by including a contracting Mogi source (Mogi, 1958) at a fixed location and depth based on that resolved in Sigmundsson et al. (2015) and solving for the volume change, the depth of the top of the dike, and the depth extent of the dike. We are not attempting to determine the stresses acting on the magma chamber with this approach but are simply using the Mogi model to remove the contribution of magma chamber depressurization from the surface deformation. Thus, we do not consider the interaction of stress between the magma chamber and dike, although we recognize that there will be a small influence on dike opening from the magma chamber depressurization that we ignore.
We assumed that the measurement errors for both the Global Positioning System and the InSAR were drawn from a multivariate Gaussian distribution. To ensure a positive semidefinite variance-covariance matrix for our InSAR measurements, we assumed a 1-D exponential covariance function similar to that used in Sigmundsson et al. (2015): with a sill of 15 mm 2 , a nugget of 5 mm 2 , and a range of 20 km. The covariance function for the InSAR observations and the estimated standard deviations for the GNSS displacements (which we assume to be uncorrelated with each other and the InSAR displacements) are used to generate the variance-covariance matrix for our observations, which automatically provide the weights in the inversion. We sampled the a posteriori probability distribution using a Markov Chain Monte Carlo sampling algorithm (Mosegaard & Tarantola, 1995). Briefly, we assume that the joint prior probability distribution for the model parameters is constant and chose initial values for each. We select a trial model by taking a random step for all model parameters. We then compare the likelihood of the trial model to the initial model. The trial model is accepted if it fulfils one of two conditions: (1) the likelihood of the trial model is higher than the current model or (2) the ratio of the trial model likelihood and the current model likelihood is greater than a random number chosen between 0 and 1. If accepted, the trial model becomes the current model, and we select a new trial model by taking a random step from this model. If rejected, we take the random step from the previous model. We continue this until we obtain a representative sampling of the probability distribution. To ensure fast convergence, we perform a sensitivity test every 500 trial models to set the maximum step size for each model parameter. We ensure that all model parameters contribute approximately equally to the change in likelihood and that approximately half of the trial models are accepted (Hooper et al., 2013). For our initial model, we assumed constant magma pressure along the dike, and therefore uniform overpressure in the dike, assuming that the dike stays at a depth of constant pressure. We further assume uniform tectonic stress, the magnitude of which depends on the distance to the rift, as described above. We solved for 13 model parameters: the far-field separation due to plate spreading, five overpressures, five Mogi volume changes (one for each time period), and finally the depth of the top and depth extent of the dike. A schematic representation of the stresses applied on the dike for the different models of this study can be found in supporting information Figure S1 and an overview of the properties of all models in supporting information Figure S2. The red arrows in Figure 3 show the maximum a posteriori probability (MAP) GNSS model predictions for the model with uniform overpressure and tectonic stress. This uniform model provides a reasonably good fit for the GNSS observations. The line-of-sight (LOS) displacement vectors predicted by the MAP model are shown in the second column of Figure 2. Even though the far-field displacements measured by GNSS are fairly well predicted by the model, the InSAR near-field displacements are clearly underpredicted. This underprediction is also present for some of the GNSS stations that are close to the tip of the dike, where predicted vectors point too far northward, suggesting a lack of model opening in the tip of the dike.
The Holuhraun dike is certainly not the first observation of additional opening in a dike tip. Pollard and Mueller (1976) observed teardrop-shaped dikes at the exposed Walsen dike and Theater Canyon sill. Assuming homogeneous rock properties along the dike, additional opening in the tip of the dike must be the result of an increase in tractions on the dike walls. In our case, the most likely candidate for these increased tractions is a gradient in the tectonically induced stress field or from additional magma overpressure inside the dike. Decreased topography along the propagation path can be expected to increase overpressure, either due to decreasing lithostatic pressure, for a dike at a fixed depth with respect to sea level, or increasing magmastatic head for a dike that remains at a depth of constant mean pressure. To test if this could generate the additional opening required in the tip of the dike, we estimate the depth of the final dike segment independently from the depth of the other three. As both eruptive fissures happened above this final dike segment , it is perhaps to be expected that this dike segment could be shallower than the other three. We assume a mean pressure gradient of 25 kPa/m, corresponding to a density difference of approximately 2,500 kg/m 3 . This increases the number of parameters to solve for by one (the depth difference between the final segment and the remainder of the dike), bringing the total to 14 for this depth difference model.
The predicted GNSS deformation vectors for the depth difference model are shown by the green arrows in Figure 3, and the predicted LOS InSAR deformation for the four interferograms is shown in the third column of Figure 2. The fit of the depth difference model predictions for both the GNSS and the InSAR deformation measurements is much improved compared to the uniform model, especially for the InSAR data. The predicted depth of the top of the final dike segment is 206 221 194 m below the surface, and the depth of the other three segments is 620 660 600 m. In these values, and subsequent values like it, the normal scripted number is the MAP value prediction, and the superscript and subscript values represent the 95% probability range. The predicted extent of the dike in depth is 5710 5760 5600 m. This corresponds well with the depth of the seismicity along the dike (Ágústdóttir et al., 2016).  Sigmundsson et al. (2015). The circles indicate the 95% confidence region. Also displayed are the best fit model predictions for the constant, large imposed tectonic stress (red arrows) and variable tectonic stress models (green arrows). The red trace shows the path of the dike used in the modeling.
The time evolution of the magma overpressure is displayed in Figure 4. For the uniform overpressure model (Figure 4a), the overpressure continues to increase over time, albeit at a slower rate, reaching a maximum of just below 10 MPa. For the variable overpressure model (Figure 4b), the magma overpressure seems to flatten off after 24 August at around 8 MPa. This flattening might be expected once the dike stops propagating, and the flow reaches a steady state. In the final dike segment, however, the magma overpressure rises to almost 19 MPa for the depth difference model, due to the final segment being shallower. This value is very high and likely higher than the host rock can sustain without the dike propagating further. Furthermore, the far-field separation due to plate spreading for the depth difference model is only 0.47 0.65 0.26 m, resulting in almost no strike-slip motion on the final, northernmost dike segment. This contradicts observations on the ground (Hjartardóttir et al., 2015), previous modeling results , InSAR offset tracking results (Ruch et al., 2016), and focal mechanisms (Ágústdóttir et al., 2016), all of which indicate significant strike-slip motion on the dike and fault. . The second column shows the large imposed tectonic stress model prediction converted to the radar line of sight, and the third column shows the variable tectonic stress model predictions converted to the radar line of sight. Both models have uniform overpressure. The large imposed tectonic stress model has a 2-m far-field separation imposed on it, while the variable tectonic stress model has different plate spreading applied to the first two and last two dike segments. Positive displacements are displacements toward the satellite. White area in the background is area covered by the ice cap. The red trace shows the path of the dike in the model. The depth difference model shows that the deformation measurements can be fit well by creating additional opening in the tip of the dike. There are two factors contributing to this improved fit. First of all, there is increased opening due to the additional overpressure. A second, smaller effect is that the lower depth makes the opening due to the imposed stress be more focused toward the top of the dike. This increases the deformation close to the rift and alters the shape of the deformation pattern slightly, especially close to the dike. This model thus has more freedom to trade off far-field and near-field deformation patterns, resulting in the best fit to the observations of all models; see Figure S2. However, the very large magma overpressure required to achieve this fit and the lack of strike-slip motion that this model predicts make the model less plausible. As magma overpressure is not able to explain the increased opening at the dike tip, it likely comes predominantly from the tectonically induced stress. We explore this option further by imposing a large far-field separation. We set the far-field separation (U in equation (1) to 2 m and shift the rift axis to line up with the final dike segment to maximize the relative influence of the tectonically induced stress field on this segment. The red arrows in Figure 5 show the model predictions for the GNSS vectors for this model. The effect of the shearing of the dike along most of its length is clearly visible in the vectors closest to the dike. The left-lateral strike-slip motion on the dike rotates the deformation vectors to the west of the dike southward and deformation vectors to the east of the dike northward, leading to an overall poor fit ( Figure S2). Furthermore, the InSAR LOS deformation predictions ( Figure 6, middle column) still underpredict the observed values.
Given that a tectonically induced stress field due to uniform plate spreading, even with its position optimized for the final dike segment, cannot explain the deformations observed, we therefore hypothesize that the tectonically induced stress field must vary along the dike. The rift axis is well constrained by GNSS in the Askja area, just north of the Holuhraun fissure (Sturkell & Sigmundsson, 2000). As mentioned above, ground observations of shearing of the graben and strike-slip focal mechanisms at the tip of the dike suggest that the deviatoric stress field at the tip of the dike largely follows this direction. However, further south, there are fewer GNSS observations. The magnitude and direction of the stress field could be very different here. This hypothesis is supported by earlier modeling of the dike propagation path , where it was found that the propagation path, influenced by both topographically and tectonically induced stress, followed the predicted path best in the last two, northernmost segments of the dike.
We investigate the possibility of a gradient in deviatoric stress field by applying a variable tectonically induced stress field to different parts of the dike. We fix the position of the rift axis to the original location and solve for two different far-field separation parameters of the arctangent model, one for the first two segments of the dike in the south and another for the last two segments in the north. The green arrows in Figure 5 show the model predictions for the GNSS vectors for the variable tectonically stress model, and the right column of Figure 6 shows the InSAR LOS predictions for this model. Even with a simplified representation of what is likely a more complicated tectonic stress field, the variable tectonic stress model is able to fit the measurements well. The fit of this model is slightly worse than that of the depth difference model (see Figure S2), likely due to the fact that the variable tectonic stress model has less freedom in balancing the far-field and near-field deformation patterns compared to the depth difference model. Figure 7 shows the maximum likelihood opening predicted for each patch during the five time periods for the variable tectonic stress model. The maximum opening is just over 4 m, in the final segment of the dike. The maximum strike-slip motion predicted is 0.8 m. The predicted far-field separation parameters are 0.2 1.0 0.0 and 3.2 3.6 3.0 m for the southern two and northern two segments, respectively. The magma overpressure and volume contained in the dike are shown in Figure 8. The time evolution of the Mogi source volume change is given in supporting information Figure S3.

Discussion
Our modeling results show that the opening required to fit the GNSS and InSAR observations cannot come from the pressure differential between magma and host rock alone. As expected, the tectonically induced stress field must make a significant contribution to the tractions on the dike, especially in the final segments. We also show that the tectonically induced stress field must be significantly different between the early dike segments and the final dike segments, to create the additional opening required further north.
To explain this difference in stress field, we propose that part of the stress caused by plate spreading was released by previous intrusions in the area to the south, leading to a gradient in the tectonically induced stress field from south to north. The Bárðarbunga system has had 23 confirmed eruptions in historic times, and nearby Grímsvötn volcano over 70 (Thordarson & Larsen, 2007), making these systems two of the most active volcanoes in the country. While current monitoring equipment is able to detect shallow activity in the region well, even 25 years ago, this would have been a different story. Shallow intrusive activity, or even small eruption that did not break the ice, could well have gone undetected, releasing part of the tectonic stress to the south. Our modeling provides no insight into the absolute stress field; it only shows how much stress was released but not how much is left to be released. Our results do suggest that if a new dike were to be initiated in the near future, it would not necessarily head in the same direction and could be expected to erupt closer to the caldera.
The far-field separations found for the variable tectonic stress model are 0.2 1.0 0.0 and 3.2 3.6 3.0 m for the early and late dike segments, respectively. Assuming 20 mm/year of plate spreading (Árnadottir et al., 2009), the far-field separation for the northern segment accounts for 150-180 years of plate spreading. The previous Holuhraun eruption is thought to have occurred in the 1860s (Hartley & Thordarson, 2013), consistent with our result. The modeled far-field separation for the southern two dike segments corresponds to 0-50 years, which would imply that at least one undetected intrusion occurred beneath the ice cap in the last 50 years, releasing the stress due to plate spreading in the area.
In our modeling, we have assumed that the rift axis continues in a straight line from north to south. However, although the location of the rift axis north of the ice cap is well determined, for the dike segments further south, there is a lot of uncertainty. In this area, the rift axis may have jumped west, or rotated toward the west, to transition from the Northern Volcanic Zone to the Eastern Volcanic Zone. This potential model error is not captured by the uncertainty of our results. However, both a jump west or a rotation of the rift axis to the west would result in increased opening of the southerly dike segments, which does not fit with the observations.
The variable tectonic stress model yields a maximum of just over 4-m opening, accompanied by 0.8 m of maximum strike-slip motion in the final dike segment. This strike-slip motion is consistent with the findings Journal of Geophysical Research: Solid Earth 10.1002/2017JB015274 by Ágústdóttir et al. (2016), who found exclusively left-lateral strike-slip focal mechanisms ahead of the dike tip during its propagation northward in the final dike segment.
The volume contained within the dike on 4 September for the variable tectonic stress model is 0.53 0.54 0.51 km 3 and is very close to the volume found in Sigmundsson et al. (2015). The time evolution of the intruded volume also follows a similar pattern, with the volume increase of the dike slowing after the first two weeks of the eruption.

Conclusions
We have modeled the evolution of the Bárðarbunga dike using the boundary elements method. Although we use only 14 model parameters in this approach, compared to hundreds in the kinematic approach of Sigmundsson et al. (2015), we fit the data almost as well. Our results show that the dike overpressure rose rapidly in the first 5 days and then remained quasi-static for the remainder of the dike propagation and early eruption period. They also show that there must be a difference in the tectonically induced stress field between the north and south segments of the dike. The modeled stress due to plate spreading in the two segments furthest north agrees well with the amount of stress accumulated since the previous eruption at Holuhraun in the 1860s. Despite the dike propagation path itself being mostly consistent with a deviatoric stress field due to uniform plate spreading along a straight rift axis Sigmundsson et al., 2015), our results imply that this was not the case. In fact, looking closely at the dike propagation predictions by Sigmundsson et al. (2015) and Heimisson et al. (2015), we can see that the dike prediction is accurate in the last two dike segments and deviates slightly at the earlier segments. This is consistent with our results, which indicates that a significant part of the tectonically induced stress field has been released beneath the ice cap by one or more undetected intrusions in the last 50 years. The 2014 event at Bárðarbunga shows that stress-constrained modeling of rifting episodes can not only tell us about the present but also shed light on what might have happened in the past.