Our study area in the central and eastern parts of the EVZ is two-thirds covered by the Vatnajökull ice cap (Figures 1A,B), explaining why the knowledge of its geology is much less detailed than that of ice-free regions of volcanic zones. Four volcanic systems have been identified in the EVZ: Grimsvötn-Laki, Bárðarbunga-Veiðivötn, Vonarskard-Hagöngur and Kverkfjöll. (1) The Grimsvötn-Laki volcanic system extends from the Laki volcanic fissure (ice-free domain, south of the Vatnajökull) to the Þórðarhyrna and Grimsvötn central volcanoes (ice-covered). In the two last decades, the Grimsvötn volcano has been the most active volcano in Iceland (Larsen et al., 1998; Larsen, 2002); the most recent eruptions occurred in 1998, 2004, and 2011. (2) Bárðarbunga-Veiðivötn is a large volcanic system roughly 180 km long. It extends from the edge of the Torfajökull central volcano at its southwestern part to the Dyngjufjöll-Ytri hyaloclastite ridge at its northeastern part. It comprises the Bárðarbunga, a large ice-covered central volcano, together with the Hamarinn central volcano that lies 15 km southwest of Bárðarbunga. This volcanic system is characterized by numerous fissure swarms (eruptive and non-eruptive) trending NE-SW, presumably extending through the margins of the glacier (Jóhannesson and Sæmundsson, 1998). The Bárðarbunga central volcano was considerably active during historical times. (3) The Vonarskard-Hagöngur, situated northwest of the Bárðarbunga-Veiðivötn volcanic system, is a volcano-tectonic system in which the main central volcano is Vonarskard-Tungnafellsjökull. (4) The Kverkfjöll volcanic system, around the Kverfjöll central volcano, belongs to the North Volcanic Zone.

Only a few seismic profile images show the crustal structure beneath the Bárðarbunga-Grimsvötn volcanic complex (Darbyshire et al., 1998, 2000; Allen et al., 2002). The upper crust is believed to be thin below this area, being only about 4–5 km thick, while the Moho reaches a depth of roughly 40 km directly above the plume center and thins rapidly away from it (Darbyshire et al., 1998, 2000; Allen et al., 2002). A recent study based on gravity data has shown that the Vonarskard, Bárðarbunga, Kverkfjöll and Grimsvötn volcanoes are associated with major high density anomalies located in the upper crust (Guðmundsson and Hognadóttir, 2007), interpreted as predominantly basic intrusions. Note that, in contrast, no gravity anomaly is observed for the Hamarinn central volcano.

More than 240,000 earthquakes with seismic moment magnitude Mw ranging from −2.4 to 6.4 were recorded locally between 1990 and 2010 by the IMO seismological network; they cluster mainly around the transform zones (Figure 1A). Other large clusters of earthquakes are observed around some of the central volcanoes, particularly in the EVZ, which includes the youngest and most active volcanic systems of Iceland (Figure 1B).

The checkerboard test is a very common technique in geophysics to evaluate the resolution and sensitivity of inversion results to the data. More specifically, we have designed a synthetic stress field as the model m containing a stress tensor at each node ij. We generated synthetic data d' based on the locations and fault planes of the earthquakes in the catalog and computed the new predicted slip vector for each focal mechanism k given by the model m. By inverting the data d', we obtained the model m^ and we inspected the result by comparing the model m with m^. The aim of the following tests is to evaluate the grid size and random choice of the actual fault plane proposed by Hardebeck and Michael (2006). No noise is added to the synthetic data.

We gridded the study area with a 2.5 × 2.5 km mesh to obtain on average 30 focal mechanisms per cell (Figure 3). Because earthquakes cluster under each central volcano and are located in a limited range of depth typical for each volcano, we adopted a 2D grid as a first order stress field (see Figures 1B, 3). We also performed experiments using different grid sizes such as 5 × 5 km and 1 × 1 km. The grid size of 2.5 × 2.5 km has the advantage of giving a high number of cells associated with a stress tensor to resolve kilometric stress deviation, while keeping an average of 25–30 focal mechanisms per cell, which is statistically relevant for a representative stress tensor. On the one hand, reducing the grid size significantly diminishes the average number of focal mechanisms per cell. On the other hand, taking a larger grid size increases the number of focal mechanisms per cell and decreases the number of computed stress tensors. In this case, the stress field is resolved but small stress variations would not be revealed.

We produced a simple synthetic model m as the strike slip regime with σ→₁ alternatively heading to N 000° and N 045° for each adjacent node ij. Then, we ran two inversion tests with (1) no damping (e = 0) and (2) the optimal e-value of damping chosen at the lower left corner of the trade-off curve, which yields the best compromise between misfit and model variances (Figure 2A). Stress tensors are calculated if at least 8 focal mechanisms are contained per cell for both tests.

For inversion with no damping, the results show that the stress field m^ issued from the inversion of the synthetic data gave a good reproduction of the synthetic stress field m in terms of stress orientations in all cells (Figure 3). The values of the fourth unknown Φ are on average also well-retrieved from the inversion. Examining the plunge of σ→₁, σ→₂, and σ→₃ on a triangular diagram (Figure 3D) shows that all stress tensors are strike-slip in type with σ→₂ within 20° of the vertical. This first test reveals that a random choice of the nodal planes does not significantly affect the resulting stress field. Thus, the size of the grid is well-adapted to the inversion.

For the next inversion, we used a value for the damping parameter e = 0.6, corresponding to the optimal value given by the trade-off curve depicted in Figure 2A. Figure 3 shows that the stress orientation of the principal axes is well reproduced with a damped stress inversion. Examining the plunge of σ→₁, σ→₂, and σ→₃ on a triangular diagram (Figure 3H) reveals that all stress tensors remain strike-slip in type, although more stress tensors show a σ→₂ away from the vertical than do the ones from the model with no damping. These results also reveal that cells containing a low number of data are first influenced by cells with a higher number of data when the damping parameter increases (Figure 3E).

These two checkerboard tests show that the sensitivity of the damped inversion is dependent on the number of data per cell and thus on the grid size. The synthetic stress field is well reproduced, despite the random choice among the two nodal planes. One should note that the heterogeneity of data within a cell is not evaluated with these tests.

We performed a damped stress inversion on the selected earthquake catalog described in section The Earthquake Catalog Used for Inversion. We used the grid size chosen from the checkerboard tests, as well as the same minimum number of focal mechanisms per cell to compute a stress tensor (Figure 3). Because the solution is not unique for such an inversion method, we computed four damped inversion models as a function of the scalar damping parameter e to evaluate the stability of the stress field in each cell. The e-values were chosen from the optimal value at the lower left corner of the trade-off curve (e = 0.6) compared to an a priori over-damped model (e = 1.4; Figure 2B).

Considering all the models, the effect of damping can clearly be seen in Figure 4, where a smoother model was obtained with high damping. All type of faulting were revealed, but normal and strike-slip faulting dominated (Figures 4, 5). In detail, the ternary diagrams show a continuum regime between strike-slip and normal faulting, revealing the occurrence of stress permutations between σ→₁ and σ→₂. This continuum is in agreement with the high stress shape ratio (e.g., Φ = 1) determined in the models (see Table 1). Nevertheless, stress permutation between σ→₂ and σ→₃ associated with a low stress shape ratio also occurs and this undergoes a slight decrease with the damped model e = 1.4. In a first approximation, the stress trajectory is rather scattered for the model e = 0.6. However, according to the geographical location, some trends remain stable irrespective of the e-value (Figure 4). Thus, the stress field can be divided into one of four bins as a function of the central volcano: Bárðarbunga (B), Grimsvötn (G), Kverkfjöll (K), or Hamarinn (H). For convenience, we calculated the minimum horizontal stress σ_hmin as defined by Lund and Townend (2007) to read the trend of the stress trajectory (Figure 6). (1) The σ_hmin in the H bin trends NW-SE with two main peaks at N135E and N090 for damping of e = 0.6. Upon increasing the damping e-value, the preferred orientation reaches N120 to N130. (2) The G bin shows more scattered σ→_hmin trends. However, the two main peaks are in the range of N060-90 and N000-N010. Damping tends to decrease the height of the first peak and to converge the σ→_hmin orientations toward a more narrow range between N070 and N080. (3) The stress pattern of the B bin is the most constant, showing a main peak in the range of N090 to N100, with few σ→_hmin axes tending perpendicular to this direction. When damping is strengthened, the main peak lies from N070 to N080. (4) The K bin constantly shows the same direction of σ→_hmin, which ranges from N010 to N030.

When the scalar damping parameter e increases, the main stress trajectories are more clearly highlighted. This could suggest that the main tendency observed from the stress trajectories is statistically strongly required by the data, although some local perturbations are maintained even with a high damping as shown in Figure 6. Here, we only consider the trends that remain stable irrespective of the e-values. These correspond to a stress field strongly required by the data.

To evaluate the stability of stress inversion in each cell, we generated a data vector that includes random noise and inverted this to evaluate the propagation of the data errors. This experiment is similar to the “Monte Carlo” tests commonly used in seismic tomography.

To carry out error propagation tests with the damped stress field models determined in section Stress Field Results, we first added random perturbations on the strike, dip and rake of our dataset in the range of the uncertainties of the fault plane solution (±15°, section The Earthquake Catalog Used for Inversion), following a normal distribution (Figure 7D). We then inverted the noisy data and finally evaluated the difference between the noisy models and the optimum damped model (e = 0.6) determined in section Stress Field Results. A convenient way to evaluate error propagation is to compute the scalar product of the principal stresses between a perturbed stress tensor and a non-perturbed one at the ij cell, giving the angular deviations. For convenience, the angular deviation for each of the principal stresses is noted herein as σ→₁ − σ→′₁, σ→₂ − σ→′₂, and σ→₃ − σ→′₃.

The results are presented for each principal stress axis in Figure 7. The histograms in Figures 7E–G reveal the mean angular deviations associated with the mean errors in the range of 10 ± 13°, 14 ± 16°, and 9 ± 12° for σ→₁ − σ→′₁, σ→₂ − σ→′₂, and σ→₃ − σ→′₃ respectively, showing relatively small deviations. Interestingly, σ→₃ is the most stable stress axis while σ→₂ is the least stable one. The number of cells with a deviation of under 10° for σ→₃ − σ→′₃ is 59 (73% of the cells). The surface map deviations depicted in Figures 7A–C show that cells with high deviations (more than 25°) are systematically found on either a and c or b and c maps. For instance, cell 53 for σ→₂ − σ→′₂ is associated with a high angular deviation that can also be found in the same cell for σ→₁ − σ→′₁ but not for σ→₃ − σ→′₃ (Figure 7). This suggests that stress rotation occurs around the σ→₃ stress axis, which remains the constant axis. Nevertheless, cell 16 is found to be deviated for all the principal stress axes (Figure 7). Stress deviations are therefore mainly expressed by rotations along one of the principal stress axes. It is not surprising that σ→₂ is the least stable stress axis since it is the common denominator for stress permutations such as σ→₁/σ→₂ and σ→₂/σ→₃. However, compared to the damped stress field models, we already pointed out the capacity of the stress field to permute mainly between σ→₁ and σ→₂ as shown on the ternary diagrams (Figure 5). So this stress permutation revealed through the stress field models is in good agreement with the results of the error propagation test. Applying tests with different damping parameter values (i.e., e = 0.8, 1, 1.4) yielded lower mean values concerning the deviation for each principal stress axis, keeping σ→₃ as the most stable axis. Adding larger noise, of the order of ±30° following a normal distribution, led to an increase in the mean angular deviation for each stress component of the order of 3°. Although more cells are disturbed, 44 cells for σ→₃ − σ→′₃ have a deviation of under 10°, which represents 58% of all cells. Overall the error propagation tests show that the computed stress fields are stable.

Prior to the eruption, no particular seismic activity was observed. A significant seismic activity started on the 29th of September 1996 with a main shock of 4.8 moment magnitude (Figures 8B,C). The eruption onset lasted from the 30th of September to the 13th of October and took place on a 6 km long Gjálp fissure (Guðmundsson et al., 2004, Figure 1B), while dyke intrusion is believed to have occurred around the Barðarbunga central volcano (Pagli et al., 2007a). The selected cluster includes 184 earthquakes starting from the main shock on the 29th of September 1996 until the 25th of October 1996 (Figures 8A,C).

Figure 9 shows the general distribution of the strike, dip and rake of the cluster. Two main strike orientations are noted in the range of N095 to N115° and N195 to N210° (Figure 9A). Both nodal planes show dips that can be either relatively gentle (from 45 to 70°) or steep (from 70° to 90°) (Figure 9B). The rake presents two main peaks in the range of −45 to 30° and −170 to 150° (Figure 9C). This cluster denotes a major tendency for strike-slip mechanisms during eruption.

To detail the state of stress, we used another approach to the one that served for the regional stress inversion. A random stress search from the Fsa package (Célérier, 2013, see section Random Stress Search) allowed interactive detection of the mechanically homogenous data subset.

A first stress inversion was applied to the 184 focal mechanisms. It yielded a reduced stress tensor which barely explains all the data. A second stress inversion was run but this time 30% of data were randomly selected and used for inversion to detect different datasets that are mechanically consistent. The resulting solution shows a reduced stress tensor with σ→₃ orientation of N065° associated with a plunge axis of 22° and σ→₁ orientation of N169° associated with a plunge axis of 30° (Φ = 0.8), thus defining a strike-slip regime. This solution explained 78 of the focal mechanisms with an angular α misfit lower than 30°. So this reduced stress tensor was used to reject data that yields an angular α misfit higher than 90°, to obtain a mechanically consistent data subset. Two subsets were thus sorted out, a first one containing 104 focal mechanisms and a second one containing 80 focal mechanisms. From this stage onwards, stress inversions are run on each subset.

Before inverting the data, a P, B, and T-axis analysis for 104 focal mechanisms is performed and presented in Figure 10. Interestingly, the P and T axes have a preferred orientation and the T-axis ranges between N055 and N075° (Figures 10A,B). According to the ternary diagram, most focal solutions are strike-slip in type although reverse and normal ones also occurred. A random stress inversion was applied to this data set. The 10 best reduced stress tensors are retained and these are shown in Figures 10D–H. The orientation of σ→₃ for the best solution is N073° with a plunge of 28° while the orientation of σ→₁ is N333° with a plunge of 17° for a stress shape ratio of 0.72, characterizing a strike-slip regime away from the Andersonian position (Figure 10E). This solution explains 83% of the data (with an angular misfit inferior to 30°, Figure 10G). The average angular misfit of the best solution was 15.9°. As for the regional stress field, an error propagation test was run, showing that the parameters of the reduced stress tensors remain very stable (see section Error Propagation Tests). The 10 best solutions obtained are in agreement with the regional stress field obtained over the 20-year span of this study (Figure 6).

As for the first subset, a P, B, and T-axis analysis was performed with the remaining 80 focal mechanisms and the results are presented in Figure 11. A preferred orientation of the T-axis is found; T-axes are mainly orientated in the range of N150–170° (Figures 11A,B). In terms of mechanisms, strike-slip regimes are dominant, although many focal solutions range between reverse and strike-slip positions in the ternary diagrams (Figure 11C). A first stress inversion was performed, yielding an average angular misfit of 20°. This misfit is acceptable, but it is preferable to refine it and run a second stress inversion, to remove the focal solutions with an angular misfit of higher than 90°. The second stress inversion yields an angular misfit of 14°, removing 5 focal solutions. These solutions are satisfactory and very stable, according to the estimation of the 10 best reduced stress tensors. An error propagation test was performed and this confirmed that the reduced stress tensor parameters are stable. Moreover, the best estimated stress tensor explains more than 90% of the data. This stress tensor solution shows a σ→₃ axis of N322° with a stress shape ratio Φ = 0.6 and is characterized by a strike-slip regime away from the Andersonian positions (Figure 11).

Two states of stress are well-identified: (1) the strike-slip regime with the σ→₃ orientation of N073° (called T1 hereafter), and (2) the strike-slip regime with the σ→₃ orientation of N142° (called T2 hereafter). To illustrate more clearly the difference between the T1 and T2 tectonic regimes, we computed the distance between the 10 best reduced stress tensor solutions of T1 and T2. Lisle and Orife (2002) and Orife and Lisle (2003) defined a measure of the distance between two stress tensors called stress difference, expressed by a single scalar D that takes into account not only the stress orientation but also the stress magnitude (the reader is referred to these references for further details). This parameter D varies from 0 to 2, where 0 means that the two stress tensors are identical and 2 means that they are strictly opposite. The results are shown in Figures 12A,B. It is clear that T1 and T2 are opposite stress tensors since the D values are in the range of 1.8–2. The T1 regime represents 56% of the cluster while the T2 regime represents 41%. 3% (including the main shock) remain unexplained by either T1 or T2.

To observe the occurrence of the T1 and T2 tectonic regimes over time, we simply plotted the focal solutions associated with T1 and T2 as a function of time (Figure 13A). Furthermore, to evaluate rotation of focal mechanisms with time, we consider three-dimensional rotations (3-D) of the geometry of focal mechanisms by an approach such as that described by Kagan (2007). This method considers 3-D rotations by which one double-couple focal mechanism can be turned into another arbitrary focal mechanism. This 3-D rotation is expressed by a single misfit angle noted β. This misfit β angle varies from 0° when the pair of focal mechanisms is identical to a maximum angular misfit of 120°. To calculate β, focal mechanisms need to be compared to an arbitrary frame. We introduced here the frame Q = (p→, b→, τ→) where p→, b→, τ→ are unit vectors parallel to the seismological P, B, and T axes respectively (here we used the same frames as Célérier, 2008, 2010). T1 was defined by the frame S_T1 (see section Stress Inversion Principle) and was used to set up an arbitrary frame Q_T1 such that: (6)ST1=QT1 Yielding a misfit β = 0 between these two frames. Next, each focal mechanism was defined by its frame Q_k (where k = 1, 2, 3…N the number of focal solutions). This was associated with the tectonic regime T1 compared to Q_T1 and the 3D rotation angle β calculated. One should note that Q_T1 (i.e., the arbitrary frame) and Q_k have to be compared in a common reference frame. We use a geographical frame G = (g→₁; g→₂; g→₃) where g→₁, g→₂, and g→₃ are unit vectors pointing north, east and down, respectively. The same process was carried out for T2, where Q_T2 is the second arbitrary frame set up. Rather than the absolute amount of rotation, we considered the relative rotation to qualitatively observe geometrical rotations of the focal solutions over time. This method is simple and has the advantage of not involving stress inversion. The results of the focal solutions with time are presented in Figure 13A. Also, the focal solutions and the rotation β misfit are presented on maps as functions of T1 and T2 (Figure 13B).

From Figure 13 it is obvious that both tectonic regimes occurred within the same time periods, sometime even within a couple of minutes of each other. Such flip of focal mechanisms over a short time period was seen in the Upptyppingar swarm in 2007 in Iceland, northeast of the glacier (Jakobsdóttir et al., 2008; White et al., 2011). Substantial changes in rotation of focal solutions were mainly observed immediately prior to the eruption (approximately 15 h before its onset) and during the first days of the eruption from 30th of September to 2nd of October). Afterwards, the rotation remains almost standard, but the amount of rotation is relatively greater than during the first days of the eruption. Also, both tectonic regimes occurred spatially in close proximity to each other, so no particular distribution was observed. This feature was also highlighted in the Upptyppingar swarm (Jakobsdóttir et al., 2008; White et al., 2011). However, then the location of the maximum amount of rotation seemed to be located northern of the Barðarbunga volcano where the main shock occurred, and in an ice-free domain (Figure 13B).

The different stress fields obtained by the Hardebeck and Michael (2006) approach revealed a dominance of normal and strike-slip faulting associated with high stress shape ratio (e.i., σ₁ ≈ σ₂). In detail, the ternary diagrams (Figure 5) show rather a continuum between the normal and strike-slip faulting. These features suggest the occurrence of stress permutation between σ₁ and σ₂ in agreement with the stress shape ratio characterizing seismotectonic behavior. Stress permutations between σ₂ and σ₃ also occur but they remain a minor feature. It is well-established that stress permutations such as switching between normal to strike-slip faulting occurred widely throughout Iceland for the present-day stress field (Bergerat et al., 1998; Garcia et al., 2002; Bergerat and Plateaux, 2012) and for paleo-stress fields (Bergerat et al., 1990; Villemin et al., 1994; Guðmundsson et al., 1996; Bergerat and Angelier, 1998; Plateaux et al., 2012). Local conditions such as the uniaxial extension (σ₁ = σ₂ > σ₃), fluid transport and dyking can effectively affect the stress field locally. However, with the database and the method used here the question arises as to whether this permutation reflects oblique mechanisms (i.e., oblique slips) or rather two kinds of pure mechanisms. To answer this question, we selected the cell with the greatest amount of data. One hundred and four focal mechanisms are located around Bárðarbunga volcanoes as shown in Figure 4 at node 15–16 (stress tensor solution n° 48 in Table 1) as a representative cell. We plotted the P, B, and T axes of the 104 focal mechanisms in this cell.

This test does not imply any stress inversion and gives a “direct” observation of the focal mechanisms. The results are presented in Figure 14 and reveal a continuum between normal and strike-slip types. This strongly supports the hypothesis that the stress regime in the studied area is a continuum between normal and strike-slip faulting.

In terms of stress orientations, two preferred directions of σ_hmin were determined: (1) a direction in agreement with the spreading direction N°105 E according to DeMets et al. (1990, 1994) and (2) a second direction deviated and perpendicular in some places. The damped regional stress field shows that the deviated stress direction is sustained through time mainly around the Bárðarbunga and Grímsvötn central volcanoes while the second direction prevails around the Hamarinn volcano. Interestingly, a gravimetric study of Guðmundsson and Hognadóttir (2007) highlighted that, in contrast to the large central volcanoes such as Grimsvötn and Bárðarbunga, no gravity anomaly is associated with the central volcano of Hamarinn. This volcano is lacking collapse caldera and shares a fissure swarm with the larger volcano of Bárðarbunga. Guðmundsson and Hognadóttir (2007) thus proposed that one of the two central volcanoes within the Veidivötn-Bárðarbunga volcanic system (Figure 1) becomes the dominant one (here the Bárðarbunga) and serves as the main conduit of magma supply, while the Hamarinn is supplied sporadically through sharing-dyke with the dominant one. This might explain why no deviation is observed in the Hamarinn volcano over the period of the study (20 years).

The same preferred directions of σ_hmin were determined for the present day stress and the state of stress during the 1996 eruption. However, it appears that the deviated stress direction prevails. Deviated states of stress from both the regional stress field and during the volcanic crisis show similar stress orientations and stress shape ratio (i.e., σ₁ ≈ σ₂).

White et al. (2011) have demonstrated that the flip of normal to reverse focal mechanisms may occur within a minute and in the same location. The two opposite states of stress characterized by the strike-slip regime that were highlighted during the volcanic crisis in 1996 support the same observation of two opposite fracturing mechanisms. Because the two T1 and T2 states of stress show stress regimes distant from the Andersonian positions, P, B, and T axes should be interpreted in terms of fault and slip geometry rather than in terms of stress (Célérier, 2010). White et al. (2011) proposed that the simultaneous normal and reverse sequence mechanisms may be associated with the progressive melt intrusion of a dyke. We propose that during the 1996 eruption, magma intrusion may have caused such sequence as well.

The stress directions described here have been observed elsewhere in Iceland. For instance, a large swarm of earthquakes occurred in 2007 north of the Vatnajökull in the Upptyppingar region at a depth of 14–17 km and it is believed that this was due to dyke intrusion (Jakobsdóttir et al., 2008). Analyses of P, B, and T-axes from fault plane solutions achieved by the Icelandic national survey showed two preferred directions of the T-axis: one is in agreement with the direction of regional spreading and another one lies perpendicular to this (Jakobsdóttir et al., 2008; Plateaux et al., 2012). Interestingly, such an observation has been described for other volcano-tectonic regions such as Mount Etna (Cocina et al., 1996; Patané and Privitera, 2001), Guagua Pichincha in Ecuador (Legrand et al., 2002), the Soufrière Hills volcano (Roman et al., 2008) or Mount St. Helens in Washington (Letho et al., 2010). These studies propose that intruding magma must be the cause of these major changes in the direction of the pressure or tension axes. Thus, we argue that the occurrence of the stress deviations that we observe must be the consequence of volcanic processes prevailing locally on the tectonic ones, not only during a volcanic crisis but also over several years.

Iceland is currently undergoing glacial isostatic adjustment due to melting of the glaciers (Árnadóttir et al., 2009). Moreover, according to these authors our study area is the fastest rebounding area in Iceland. Therefore, we can expect the regional stress field to be affected by elastic rebound due to glacier thinning. On the base of numerical models of magma transport and isostasic rebound, Hooper et al. (2011), Jull and McKenzie (1996), Maclennan (2002), Sigmundsson et al. (2010a,b), and Pagli and Sigmundsson (2008) showed that decompression due to the recent retreat of the glacier will significantly increase the rate of melt generation in the mantle. However, Sigmundsson et al. (2010a,b) argued that the mode and time scale of transport of the additional magma toward the surface may take longer than decades or centuries. So, we might expect no direct influence of deglaciation on the regional stress field during the short time period of our seismotectonic study (20 years). According to Sigmundsson et al. (2010a,b), the elastic rebound due to the effect of present day ice thinning on the failure of shallow magma chambers should remain small, compared to their natural activity. Therefore, we can discard a notable effect on the expected state of stress (in terms of stress orientation and magnitude) due to the present-day retreat of the glacier, which can be compared to our results.

A number of studies have attempted to model the stress field in terms of stress orientation and magnitude for volcano-tectonic interactions (Guðmundsson et al., 2005, Guðmundsson and Andrew, 2007). Numerical modeling of Andrew and Guðmundsson (2008) considered for their calculations the magma chambers as elastic inclusions. With these simplifications, they highlighted that a contrast of elastic properties between central volcanoes and their host matrix (e.g., lava pile) may significantly change the stress field. More specifically, the maximum principal stress (σ₁) can change in both magnitude and orientation in the host rocks immediately surrounding the elastic inclusion. These numerical models predict local stress deviation from the spreading direction around volcanoes. Pressurization and depressurization of the magma chamber must also change the stress in direction and in magnitude. Vargas-Bracamontes and Neuberg (2012) propose an interaction between the regional and the magma-induced local stress field, based on the Coulomb failure criterion, to explain the 90° shift of pressure axes (P-axes) referred therein as deviated stress direction from the regional stress field.. According to these authors, when magma pressure starts to increase, the slip (in either a regional or rotated sense of slip relative to the regional maximum compression) also increases. A further increase would tend to promote faulting in the opposite sense from faulting induced by regional stress fields.