Tidal triggering of earthquake at Axial Volcano

Axial Volcano, which is at the intersection of a mid-ocean ridge with a hotspot (Fig. 1), erupts on a decadal time scale. Each eruption is followed by a caldera collapse accompanied by thrusting on outwardly dipping ring faults, followed by a re-inflation period, at the latter stages of which the ring faults become reactivated in normal faulting11,12,13. The best observations of tidal triggering were for the normal faulting earthquakes in the months prior to the 2015 eruption7.

Fig. 1 Cross-section of seismicity at Axial Volcano. East–west cross-section of seismicity in 3 months preceding the 2015 eruption at Axial Volcano. Red curve is the roof of the axial magma chamber15. The bathymetry is from a compilation of latest cruises, the most recent in 201062. Inset, location map for Axial Volcano Full size image

At Axial Volcano, the ocean tides are very large (3 m) so that ocean loading dominates the solid earth tides and the vertical tidal stress dominates and is in phase with the ocean tides (Supplementary Fig. 1), so we need only consider the vertical component in our analysis. Tension is taken as positive for tidal stresses, so the maximum tidal stress corresponds to the minimum water depth. To avoid ambiguity, in this paper we will refer to high and low tides in the conventional way as high and low water, recalling that low water produces tension and high water produces compression.

Figure 1 shows a cross-section view of the seismicity for the three months prior to the 2015 eruption, which illuminates the ring faults. This data set contains ~60,000 earthquakes with a magnitude of completeness (M c ) = 0.114. Figure 2 shows a histogram of the seismicity plotted as a function of tidal phase, in which 0° is the maximum low tide. The correlation is obvious and requires no statistical treatment. It was first proposed that this was also a case of fault unclamping6,8,10, but when it was established that these earthquakes were dominated by normal faulting12 this viewpoint became untenable. Both the seismicity trends in Fig. 1 and the focal mechanisms12 indicate a mean fault dip of 67°. The reduction of vertical stress brought about by low tide will produce a Coulomb-stress change on such steeply dipping normal faults that inhibits their slip. It is, rather, the high tides that will produce a Coulomb stress on the faults that encourages slip. This seeming paradox is resolved by including the effect of the axial magma chamber on the distribution of stress.

Fig. 2 Seismicity versus tidal phase. Histogram of earthquakes plotted vs. the phase of the vertical component of the tidal stress, in which 0° is the peak stress (tension is positive), which corresponds to the lowest ocean height Full size image

The response of the magma chamber

The red curve in Fig. 1 delineates the roof of the axial magma chamber obtained from seismic imaging15. Inflation of the magma chamber drives the normal faulting on the ring faults. This is demonstrated in Fig. 3, where we show east–west cross-sections of a 3D model containing a magma chamber with dimensions defined by seismic imaging15,16. Figure 3a shows the Coulomb failure stress change, ΔCFS = Δτ + μΔσ, on 67° dipping faults that results from a magma chamber overpressure of 1 MPa (Δτ is the change in shear stress resolved on the fault in the slip direction, Δσ is the change in normal stress on the fault plane, and µ is the friction coefficient). Positive (red) ΔCFS values encourage normal fault slip, negative ones (blue) inhibit it. The primary features in Fig. 3a are the zones of positive ΔCFS that correspond to the seismicity shown in Fig. 1. See Methods for details about the model.

Fig. 3 Coulomb stress on normal faults above the magma chamber. Distribution of Coulomb stress changes on 67° dipping normal faults near the axial magma chamber. Positive values favor fault slip, negative values inhibit it. a For an overpressure of 1 MPa within the magma chamber with friction coefficient of the faults µ = 0.8. b For a decrease in vertical stress equivalent to a reduction in water level of 1 m. In b an effective friction µ’ = 0.4 is used. The bulk modulus of the rock is assumed to be K r = 55 GPa, and the bulk modulus of the magma chamber is assumed to be K m = 1 GPa. The heavy line in Fig. 3b indicates the fault upon which ΔCFS was measured Full size image

Because the magma chamber is a soft inclusion, its presence will profoundly affect the stress field in its vicinity resulting from any external load. We simulate the response to tides by calculating the distribution of ΔCFS on 67° dipping faults resulting from a reduction in vertical stress corresponding to a 1 m drop in the ocean tide. This is shown in Fig. 3b. The pattern is very similar to that of Fig. 3a, demonstrating how a low tide can stimulate activity on these faults. This pattern arises because the reduction of vertical stress causes the magma chamber, owing to its higher compressibility, to inflate relative to the surrounding rock, which produces a stress field congruent with that of Fig. 3a. This is superimposed on a uniform ΔCFS from the tidal stress, which is negative in the case of a low tide. Likewise, high tides cause the magma chamber to deflate, which also produces Coulomb stresses opposite in sign to the tidal ones. Which component is larger determines whether normal faulting earthquakes are stimulated by the low tide or the high tide.

The relative expansion of the magma chamber depends inversely with K m /K r , the bulk modulus of the magma chamber relative to that of the surrounding rock, so this is the critical parameter that determines the behavior of the system. In the calculation of Fig. 3a we used µ = 0.8. The lack of a phase shift between the tidal stress and the seismicity (Fig. 2) indicates that tidal loading must be under undrained conditions (see Supplementary Fig. 2 for a check on this assumption). Therefore, for calculations such as shown in Fig. 3b we use an effective friction µ' = (1 − B)µ, where B is Skempton’s coefficient. The proper value to use for B is not well established: plausible values are between 0.5 and 117, so we explored values of µ' from 0.4 to 0. In Fig. 3b we used µ' = 0.4, K m = 1 GPa and K r = 55 GPa.

The ratio K m /K r and µ' control the ΔCFS on overlying normal faults that results from an applied vertical tidal stress. We use a metric, χ, which normalizes the ΔCFS by the vertical tidal stress and which, therefore, is time invariant. We define χ by the ΔCFS on a 67° dipping fault, indicated by the bold line in Fig. 3b, averaged from the corner of the magma chamber to the surface, normalized by the applied vertical tidal stress. This parameter is shown in Fig. 4 as a function of K m for several values of µ’ at a fixed K r = 55 GPa. Positive χ values indicate that normal faulting earthquakes will be favored by low tides, negative values by high tides and vice versa for thrust earthquakes. All conditions in Fig. 4 within the red region therefore favor normal faulting earthquakes on the low tide and inhibit them on the high tide, and within the blue region, vice versa. Normal faulting may be generated by low tides for values of K m < 8 GPa, depending on the value of µ'. The point indicated by the plus in Fig. 4 is the case illustrated in Fig. 3b. The bulk modulus of gas-free mid-ocean ridge (MOR) magma is 12 GPa18, but at the pressure of the magma chamber (~40 MPa) this value can be reduced by one to two orders of magnitude by the presence of volatiles19. Thus, at this pressure, a magma of K m = 1 GPa would contain 2650 ppm CO 2 by weight18. This is greater than the highest values typically seen for CO 2 content of MOR magma20, but this difference could easily be accounted for by the inclusion of exsolved H 2 O. So, the values of K m that we find would promote normal faulting at low tide are realistic. Our illustrative example (Fig. 3b) indicates χ = 0.176, a figure that will enter into the modeling calculations of the triggered seismicity in the next section.

Fig. 4 Properties of the magma chamber deformation system. The vertical axis χ is the average ΔCFS on a 67° dipping normal fault from the tip of the magma chamber to the surface (bold line in Fig. 3b) resulting from an applied vertical tidal stress, normalized by that stress. The red area defines the conditions in which low tides encourage seismicity on normal faults and high tides discourage it, and the blue area vice versa. The plus sign indicates the conditions shown in Fig. 3b Full size image

Comparison of the earthquake triggering with theory

During the three month period prior to the eruption, the seismic moment release on the eastern ring fault indicated a slip magnitude approximately the same as observed geodetically12, indicating that the faults are seismically coupled. In this case, there are two models that relate change in seismicity rate to a rapid change in driving stress. These are based on earthquake nucleation models21, one derived from the rate-state friction law22 and the other from subcritical crack growth due to stress corrosion23,24. The rate-state friction version is

$$\frac{R}{r} = {\mathrm{exp}}\left[ {\frac{{{\mathrm{\Delta }}{\mathrm{CFS}}}}{{A\sigma }}} \right]{,}$$ (1)

and the stress corrosion version is

$$\frac{R}{r} = \left[ {1 + \left( {\frac{{{\mathrm{\Delta }}{\mathrm{CFS}}}}{{\Delta \tau }}} \right)^n} \right]{,}$$ (2)

where R is the instantaneous seismicity rate, r is the background rate, here taken as the rate when the tidal stress is zero, and ΔCFS = χσ v , σ v being the vertical tidal stress. The control parameters for the rate-state friction version are the effective normal stress σ and the viscous friction term A. For the stress corrosion version, they are the stress corrosion index n and the earthquake stress drop Δτ.

The fit of these equations to the data is shown in Fig. 5, where the solid blue curve and the dashed red curves are Eqs. (1) and (2), respectively. These two formulations cannot be distinguished and fit the data equally well. There is no detectable phase shift between the seismicity and the tides (Fig. 2), nor is there any hysteresis observed—data for rising and falling stresses fit the triggering curves equally well (Supplementary Figs. 4 and 5). We, therefore, conclude that poroelastic relaxation is negligible in the response to the semi-diurnal tides.

Fig. 5 Normalized seismicity rate change vs. change in Coulomb stress. Coulomb stress is converted from tidal vertical stress using χ = 0.176 (from the state indicated by the plus sign in Fig. 4). Blue curve is the rate-state friction version and the red curve is the stress corrosion version Full size image

The agreement with the theories is excellent, and extends them to far smaller stresses than previously seen5, even into the negative stress regime. In the case shown in Fig. 5, the value of χ used was 0.176, from the illustrative example. The goodness of fit to the theories does not depend on the value of χ obtained from the deformation model: that merely determines the scale of the stress axis. The various implications of this result will be deferred to the discussion section.

Applications to other areas

Wilcock6 searched for tidal triggering on the mid-ocean ridge systems of the NE Pacific, using mainly land-based networks. He found a 15% increase in seismicity within 15° of the lowest tides. The focal mechanisms of the earthquakes, however, were undetermined. With an OBS deployment on the Endeavour segment of the Juan de Fuca ridge, some 2° NE of Axial Volcano, the correlation of seismicity with low tides became much better defined9. Most of the triggered seismicity there was near the ridge axis, where the focal mechanisms indicate normal faulting25. This situation is therefore quite similar to Axial Volcano and the same effect of the magma chamber seems necessary to explain these observations.

At the hydrothermal field at 9°50’N on the East Pacific rise, an OBS deployment also showed evidence for tidal triggering8. There the ocean tides are much smaller than at Axial Volcano and a significant contribution to tidal stresses is from the solid earth tides. The seismicity maximum correlates with the maximum extensional tidal stress, which can reach 1.3 kPa. The dependence of the seismicity on stress is similar to that observed at Axial Seamount (compare Fig. 3c in ref. 8 to our Fig. 5). Evidence for the mechanism of the earthquakes is equivocal: scant focal mechanism data has indicated strike-slip, normal faulting, and reverse faulting26,27, and others have proposed that the seismicity is due to hydrothermally induced extension cracking28. There is also a variation in the tidal phase angle of earthquakes along the strike of the ridge axis. This indicates the earthquake triggering is also modulated by pore pressure changes brought about by hydrothermal circulation29. With this degree of ambiguity and complexity, we cannot assess how deformation of the magma chamber may be related to the tidal triggering in this location.

The unloading model used here was initially tested at Katla volcano (Iceland), where earthquakes show an annual cycle with the maximum seismicity rate occurring in the late summer30 when the snow cover of the glacier above the volcano is minimum (annual fluctuation—6 m). The model31 showed that this was also the period of maximum Coulomb stresses in the area above the magma chamber. However, in this case the focal mechanisms were not known30 so it was not possible to determine if the system was in the red or blue regions of Fig. 4.