Tidal Influence on Seismic Activity During the 2011–2013 El Hierro Volcanic Unrest

The El Hierro volcanic unrest started in July 2011, with an increase in observed seismicity rates and surface deformation. After the initial onset, hypocenters migrated southward through September 2011, culminating in a submarine eruption beginning on October 10, 2011 and finishing in February 2012. The seismic activity continued, with remarkable periods of unrest through 2012 and 2013. The most significant episodes of seismic activity during this unrest are related to magma migration at depth. In this work, we compute tidal stress for each earthquake, at its hypocenter depth, and assign them a tidal stress phase angle. We have found statistically significant correlations between the occurrence of earthquakes and tidal stress phase angles, corresponding mainly to increasing tidal stress change rates. We found primarily that the magnitude of vertical and E‐W horizontal tidal stress values and their changing rates with time were correlated with earthquake occurrence times. We also found that there is no correlation between tides and seismicity at times with no observed surface displacements, suggesting that tidal modulation might be related to overpressure during migration of magma. Tidal modulation changes with depth and the influence of ocean‐loading tides is stronger than the influence of solid Earth tides. Our results support the hypothesis that tidal stress may modulate the seismicity during volcanic unrest, particularly during shallow depth magma migration.

Plain Language Summary Soli-lunar tides have long been postulated to influence tectonic and volcanic activity. However, the exact mechanisms remain elusive. Here, we study the influence of tides during a period of volcanic unrest that produced a significant number of earthquakes. This large number of earthquakes over a long period of time provided a unique opportunity to study evidence of a robust statistical significant correlation between tides and volcano-tectonic activity. The El Hierro volcanic crisis began in July 2011, with an increase in observed seismicity rates and surface displacements in and around the island, culminating in a submarine eruption that started the 10th of October and finished in February 2012. The seismic activity continued, with remarkable periods of unrest through 2012 and 2013. We have found significant correlations between earthquakes and tides. Tidal correlation was found statistically significant only at periods with observed surface displacements, suggesting that tidal modulation might be related to overpressure during migration of magma. The influence of ocean-loading tides is stronger than solid Earth tides, at least for this island setting. Our results support the hypothesis that tides may modulate the seismicity during volcanic unrest, particularly during magma migration occurring at shallow depths.

MIGUELSANZ ET AL.
In this work, we study the correlation between tidal stress and seismicity rates during the volcanic crisis that occurred in El Hierro Island (27.7° N;18.0° W) between 2011 and 2013. El Hierro is both the youngest and the smallest of the Canary Islands and features the highest concentration of volcanoes in the archipelago . The most recent eruption prior to 2011 was in 1793 (Carracedo et al., 2001). The structure of the island is based on three volcanic rifts in directions NE, NW, and S (Figure 1a), separated by three important gravitational landslides: El Golfo to the north, Las Playas to the east, and El Julan to the south (Masson et al., 2002). Gorbatikov et al. (2013) note the existence of a solidified magma reservoir at the northwest part of the island at a depth between 15 and 25 km. The boundary between the crust and the mantle is set at 12-16 km (Martí et al., 2013;Watts, 1994). An intricate array of magma pockets at depths between 12 and 30 km may represent the main magma storage system at El Hierro (Stroncik et al., 2009).
The volcanic crisis started in July 2011, with an important increase in seismic activity and ground deformation in El Hierro Island (González et al., 2013;Ibáñez et al., 2012;López et al., 2012). In the first stage, most of the hypocenters concentrated under Tanganasoga volcano and El Golfo valley in the center and north of the island (Figure 2a). Starting on mid-September, there was a southwards migration of the hypocenters until a submarine eruption occurred offshore of the southern part of the island on October 10-12, 2011 ( Figure 2b). After the eruption onset, seismic activity was located in the north around the El Golfo area (Figure 2c). Correlation between gravity variations and tidal vertical strain was observed (Sainz-Maza et al., 2014), suggesting potential tidal triggering of earthquakes during the first days of the eruption. This eruption ended on March 5, 2012 (González et al., 2013). However, the seismic activity continued, with periods of high seismicity and deformation through 2012 and 2013  related to magma migration at depth (see Figures 2d and 2e).
Because the rates of diffuse CO 2 emissions may vary before events of volcanic activity (Carapezza et al., 2004), estimating the output of this gas has implications in early warning systems applied to volcanoes. In the case of El Hierro Island, diffuse CO 2 emission have been monitored since 1998 (Melián et al., 2014), and especially during the seismo-volcanic crisis that started in July 2011. Through the analysis of these measurements of CO 2 emission, Melián et al. (2014) found two remarkable periods of efflux increase that may be considered as precursory signals. The first one occurred between 29 September and 12 October, that is, starting almost 2 weeks prior the beginning of the submarine eruption and 10 days before an event of magnitude M = 4.3. The second significant increase in the rate of diffuse CO 2 emission occurred during the period 24 October to 27 November, preceding both the bubbling observed at the sea between 3 and 7 November and the M = 4.6 event recorded on November 11. However, the time series of diffuse CO 2 emission in El Hierro Island does not show any correlation to meteorological parameters such as rain, soil moisture or barometric pressure (see Figure 7 by Melián et al., 2014). MIGUELSANZ  , with the location of the bestfitting spherical point sources: orange, deep crustal source (crust-mantle reservoir, or CMR following the notation by González et al., 2013); and dark red, the shallower crustal reservoir. Seismicity flux (events/km 2 ), which represents the 2-D clustering of background seismicity, is shown the background. Inset shows the vertical cross section b-b' (modified from González et al., 2013).

Seismic and Volcanic Activity Related to El Hierro (2011-2013)
Following the analysis of the seismicity by Ibáñez et al. (2012), we can identify different phases in the seismic activity of El Hierro during the period 2011-2013, which are related to different spatial settings of the earthquakes or different stages in the magma migration processes.
In this work, we designate all seismicity prior to the eruption occurred on October 10, 2011, as Phase 1. This episode starts with an abnormal increase of earthquakes in July 17, 2011 (Figure 2a). During this period, magma is detected as an intrusion in the base of the oceanic crust (around 9.5 km depth), at the south of the main concentrations of earthquakes (González et al., 2013).
Around mid-September 2011, epicenters start migrating to the south of the El Golfo area (Figure 2b). GPS stations detect north-south deformation, suggesting that magma was migrating away from the magma intrusion site (crust-mantle reservoir, CMR) (González et al., 2013). Hypocenters are (on average) deeper than before, and average earthquake magnitudes also are greater (Ibáñez et al., 2012). At the end of September, magma began an upward migration, accumulating at a shallower magma reservoir at around 4.5 km depth. This reservoir likely failed around October 8 with an event of magnitude M = 4.3 that opened a path for magma toward the surface (Figure 1b). Later, between October 10 and 12, the submarine eruption began (González et al., 2013).
On October 15, 2011, new seismic activity was located in the north around the El Golfo area. We will consider this a separate stage (Phase 2; Figure 2c) in the seismic activity, as these earthquakes happen far from the eruption site, in the opposite side of the island. Hypocenters in Phase 2 are deeper than in Phase 1, with a depth range between 15 and 25 km (Ibáñez et al., 2012). Four of the earthquakes that occurred during  eruption was officially decreed in March 5, 2012. See González et al. (2013) for a detailed description on the seismic and volcanic activity during Phases 1 and 2.
After the ending of the submarine eruption, new episodes of seismic activity occurred, the most relevant (in terms of number of earthquakes) being those recorded primarily in June to July 2012 and March to April 2013. For the purposes of this research, we will refer to these two crises as Phase 3 (2012-06-14 to 2012-08-21; Figure 2d) and Phase 4 (2013-03-17 to 2013-04-30; Figure 2e). Some authors think that there is oceanographic and geochemical evidence that the seismic activity in Phase 3 is coincident in time with submarine volcanic activity offshore the western coast of El Hierro (García-Yeguas et al., 2014;, but this volcanic activity must be interpreted in terms of magma injection processes that do not necessarily culminate in an eruptive process (Blanco et al., 2015;Pérez et al., 2015). There is no evidence of new eruptions near El Hierro island after the one that occurred between October, 2011 and March, 2012, although several authors do suggest that the magmatic processes around El Hierro island continued after the conclusion of that eruption . Hypocenters in Phase 3 are located under the west wing of the island, whereas events in Phase 4 are, on average, the most distant from the island. It is during Phase 4 that the strongest seismic activity took place, in addition to the largest magnitude earthquake (M = 4.9) and the greatest number of landslides .

Earthquake Data
Seismicity data related to the El Hierro volcanic crisis were obtained from the Spanish Geographical Survey, IGN [www.ign.es]. Our main data source is a catalog featuring all earthquakes within the period 17/07/2011 to 31/12/2013, and including event number, date, time, latitude, longitude, depth, intensity, magnitude, magnitude type, and location area. Very low magnitudes are excluded because the location of small earthquakes may be affected by significant errors, especially in a small island like El Hierro . Therefore, we used only earthquakes with magnitude M c > 1.5, where M c is the minimum magnitude of completeness based on a frequency-magnitude distribution analysis using the Maximum Curvature (MAXC) technique (Mignan & Woessner, 2012). In addition, in the case of the earthquakes that occurred during Phase 1, we used the catalog of Domínguez Cerdeña et al. (2014), where all events previous to the submarine eruption with magnitude M > 1.5 were relocated using the double-differences algorithm HypoDD (Waldhauser & Ellsworth, 2000), which improved the precision in the location of earthquakes by as much as 10 times. Díaz-Moreno et al. (2015) identifies an average horizontal error of 9 km and an average vertical error of 5 km in the preliminary earthquake locations reported by IGN. According to Domínguez Cerdeña et al. (2014), location errors in the earthquakes from the IGN catalog prior to the submarine eruption have mean values of 4.3 ± 1.8 km (depth) and 4.7 ± 2.1 km (horizontal). After the relocation provided by the algorithm hyppoDD, Domínguez Cerdeña et al. (2014) tested the reliability of the new results according to three possible causes of error (time delay uncertainties, station distribution, and velocity model) and obtained maximum uncertainties of 400 m in the relocated events. This relocated catalog is accessible from the IGN website [www.ign.es]. In the IGN catalogs obtained for the El Hierro series there is only one event provided for the best-fitting double-couple fault plane solution for a set of observed first motion polarities.
Application of the Reasenberg declustering algorithm (Reasenberg, 1985) to the seismic catalog shows that for each of the four phases there exists one big cluster which includes 90%-95% of the events that occurred during that phase. The list of parameters used in the declustering process (van Stiphout et al., 2012) can be found in the supporting information (Table S1). Figure 3 shows the location of the earthquakes in each cluster C1-C4. A visual comparison between Figures 2 and 3 suggests that, in every Phase i (i = 2, 3, 4), most of the events not belonging to the corresponding cluster Ci may be a remnant or continuation of previous phases, implying that they are not necessarily related to the magma movements or physical processes which originated each phase. As a result, only the events in the clusters C1-C4 have been considered here, as we are going to study the possible modulation effects of tides on seismicity in each one of the four phases independently.
Tidal stress and strain were calculated for the locations of all events in each cluster, in order to estimate statistical correlations for tidal modulation of the occurrence of earthquakes in the different phases of the seismic crisis. We want to note that the polarity and magnitude of tidal stressing depends on the orientation of the receiving faults, but the available focal mechanism information is insufficient to be considered in this study.

Calculation of Tides
In order to obtain the tidal stress at the earthquake locations, it is necessary to consider both the solid tide and the ocean tide. In addition, the calculation must include the depth of the events. Strains derived from the solid tides at depths less than 25 km are not significantly different from surface strains but strains due to the ocean-loading component may change significantly between the surface and a depth of 25 km (Thomas et al., 2012).
Our method uses the software package SPOTL (Agnew, 2012) for the calculation of tidal strains. Body tides are calculated on an elastic Earth following Munk and Cartwright (1966). The ocean tidal model we have used is osu.tpxo72atlas (Agnew, 2012), which is a hydrodynamic model assimilating altimetry data from TOPEX/Poseidon and Jason. It combines three basin-wide solutions (Atlantic Ocean 2011-atlas, Pacific Ocean 2011-atlas and Indian Ocean 2011-atlas), each one of which also incorporates a set of high-resolution local models. Tidal waves considered in the study are M m , M f , Q 1 , O 1 , P 1 , K 1 , N 2 , M 2 , S 2 , and K 2 . The original SPOTL codes where modified to use depth-dependent Green functions in the calculations, following the methodology described in Royer et al. (2015) for the modeling of tidal strains at depth. We added the load MIGUELSANZ ET AL.  tides to the body tides in order to give the total tidal strain. Then tidal stresses are derived from tidal strains following Hooke's law for a tridimensional isotropic body (Fischer et al., 2006), with Poisson's ratio ν = 0.25 and rigidity μ = 30 GPa. We consider positive stress as extension.
To study the correlations between tides and earthquakes, we assign a phase angle, ϕ, to each event, which can be evaluated later in statistical terms. For this purpose, a golden section search routine (Press et al., 1992) is used to estimate the extreme tidal values closest to the time of occurrence of the event. Finally, we assign the phase angle ϕ to the earthquake by linear interpolation, assuming the time distance between two tidal peaks defines a complete cycle of 360° (see Figure 4).
As explained in Section 3.1, earthquake location may be affected by both horizontal and depth errors. Table S2 in the supporting information describes the possible effect of these location uncertainties in the results.

Statistical Methods
The Schuster test was used in order to determine whether there is a nonrandom distribution of the calculated tidal phases (Emter, 1997). In this test, each ith earthquake is assigned a unit length vector in the direction of the phase angle ϕ i . The modulus R of the vector sum R over the total number of earthquakes N tot is given by where C is a factor correcting local asymmetries in tides (Wilcock, 2009). It must be noted that no significant tidal asymmetries where found in El Hierro in the time periods considered in this study, as all C values observed for tidal strain or tidal stress are in the range 1 ± 0.0020.
Then, the probability P s that the phase distribution to be nonrandom is estimated by (2) P s represents the significance level with which to reject the null hypothesis that the earthquake phase angles represent a completely random distribution. Small values of P s (<0.05) indicate the existence of a high correlation between tidal stress and the occurrence of earthquakes.
As we are treating highly clustered seismicity, we are not dealing with independent earthquakes. Given the high sensitivity of the Schuster test (Emter, 1977), that means the catalog may be biased in a way that favors the rejection of the null hypothesis. One way to assess the significance of the observed P s values is the use of a Monte Carlo permutation test (Fortin & Jacquez, 2000;Noreen, 1989). For each of the four clusters C1-C4, we order the earthquakes according to their phase angle ϕ, and we divide them into segments (phase angle ranges) of equal length (different segment lengths are tested, see next section). Then, a number NSIM of simulations (e.g., 10,000 simulations) are performed for each cluster with the segments in random order, obtaining a probability value P s * according to the Schuster test for each simulation. Finally, for each cluster we count the number NSIG of random simulations in which the P s * value obtained is lower than the P s value of the original distribution of events. The ratio provides the significance level of our test (Noreen, 1989).
In addition, the American Statistical Association expresses concern regarding the hypothesis testing based on p-values, stating that these practices are prone to generate false positives (Wasserstein & Lazar, 2016). Bayarri et al. (2016) propose a simple alternative to the use of p-values in testing null hypothesis, based on the Bayes factor B, which is defined as  Average likelihood of the observed data under the alternative hypothesis . Likelihood of the observed data under the null hypothesis B (4) Although B may be difficult to compute, it can be approximated by an upper bound B which is a function of the p-values obtained with the usual statistical tests used in testing hypothesis. According to Bayarri et al. (2016), and considering the Schuster test as the appropriate one for our null hypothesis testing: Here, we use the p-values P s calculated by means of the Schuster test to derive this upper bound B of the Bayes factor B. According to Bayarri et al. (2016), the standard significant threshold for rejecting the null hypothesis when using this Bayesian alternative is B > 16.
The use of this Bayesian method is a new application to the assessment of correlations between tides and earthquakes, and it has implications for earlier studies of tidal triggering or modulation of seismicity where the statistical significance was primarily based on p-values.

Results
Results for the application of the Schuster test to the earthquake clusters defined in Section 3.1 are shown in Table 1. The upper bound B of the Bayes factor B is also shown in Table 1 in order to compare results.
According to the Schuster test, very high correlation (P s < 0.05) appears in three of the four clusters (C1, C3, and C4) between tidal confining stress (i.e., the trace of the tidal stress tensor) and the occurrence time of earthquakes. However, when we consider the value of the upper bound of the Bayes factor, only the correlations for clusters C1 and C4 are statistically significant.
10.1029/2020TC006201 7 of 21  (Wilcock, 2009). The terms Ah 1/2 and L are the amplitude and time length of the half cycle in which the event occurs, respectively. x is the time difference between the occurrence time of the earthquake and phase ϕ = 180°, positive when the event occurs after phase ϕ = 180° and negative when it happens prior to phase ϕ = 180°. The figure also shows the graph of tidal stress rate (blue line). Note that the maximum increasing tidal stress rate occurs around tidal stress phase ϕ = 270°, whereas the maximum decreasing tidal stress rate happens around tidal stress phase ϕ = 90°. Table 2 shows the results of a Monte Carlo permutation test (10,000 simulations) applied to the earthquake clusters C1-C4, following the methodology explained in Section 3.3. For each cluster, the corresponding catalog has been divided into 72 segments of 5° width. Different segment sizes may produce different results, because the greater the segment length is, more patterns from the original distributions are repeated in the simulations and the number of possible permutations diminishes (Bhatnagar et al., 2016). As a result, a catalog divided into 144 segments of 2.5° width was also tested (see Table 2). According to this permutation test, the correlation is statistically significant for clusters C1, C3, and C4, although the value of the significance level, MC sl , provided for cluster C3, is clearly higher (implying a weaker statistical significance) than those obtained for clusters C1 and C4. On the other hand, the random shuffling of segments in cluster C2 distribution produces a great number of P s * values lower than the original P s value for both segment lengths considered here, meaning that the correlation is not statistically significant. Figure 5 shows the histograms of the frequency distributions of tidal confining stress phase angles for clusters C1, C2, C3, and C4, respectively. In these histograms, earthquakes are collected in phase angle ranges of 30° width. The sinusoidal curve applied to the frequency distributions is obtained by least squares fitting of the expression where ϕ is the phase angle, P 0 is the mean frequency, P 1 is the amplitude of the curve, and α is the phase of the curve (Tanaka et al., 2002). A triangle marks the peak of the fitted curve associated to the histogram, which occurs on the phase angle maximizing Equation 6. Figure 5a shows that the histogram corresponding to cluster C1 (coincident with Phase 1 in the seismic crisis, i.e., all the seismicity prior to the starting of seismic tremor signal) is clearly biased toward the right side of the phase angle axis, the phase range where tidal stress is increasing (180° < ϕ < 360°). The fitted curve peaks around phase angle ϕ = 251°. As seen in Tables 1 and 2, the correlation between confining tidal stress and the occurrence of earthquakes in cluster C1 is remarkably high for every statistical test. The histogram corresponding to cluster C2 ( Figure 5b) does not show a statistically significant correlation with any of the applied methods, although the most populated phase angle range is 270° < ϕ < 300°. Cluster C3 corresponds to Phase 3, the seismic activity recorded during June to August 2012 under the northwest rift.
MIGUELSANZ ET AL.  Note. Different segment lengths were applied. P s is the probability that the phase angle distribution in the original catalog to be random, according to Schuster's test. For each permutation, a Schuster probability value P* is obtained. NSIG is the number of permutations where P*≤ P s for each cluster. SL means "significance level." The significance level MC sl of the permutation test is obtained according to expression 3 in Section 3.3. Occurrences where P s < 0.05 or MC sl < 0.05 are in bold. According to the Schuster test and the Monte Carlo permutation test, there also is high correlation between tidal stress and the occurrence of earthquakes, although both the P s statistical and the ratio MC sl are not as conclusive as those obtained for cluster C1. Besides, the correlation is not significant if we consider the application of the upper bound B of the Bayes factor B ( Table 1). The fitted curve for cluster C3 distribution has a peak around ϕ = 209° (Figure 5c). Finally, the histogram obtained from cluster C4 data (Phase 4, which includes the seismicity observed during 2013 March to April and expanding offshore eastwards) is similar to the cluster C1 histogram, with a distribution biased toward the phase angle range where confining tidal stress is increasing (180° < ϕ < 360°), and a fitted curve whose peak lays around ϕ = 288° (Figure 5d). All the statistical tests considered for cluster C4 data suggest very high correlation (Tables 1 and 2).
As some of the distributions (namely the histograms corresponding to clusters C1 and C4) show an apparent bias toward phase range 180° < ϕ < 360°, a simple binomial test (see e.g., Wonnacot & Wonnacot, 1977) was performed in order to evaluate if the events are more prone to happen in the increasing-tide half cycle (180° < ϕ ≤ 360°) rather than in the decreasing-tide half cycle (0° < ϕ ≤ 180°). The results of applying that test to the four clusters appear in Table 3 and a poor correlation (not statistically significant) for clusters C2 and C3.
The upper bound B of the Bayes factor B has also been calculated for the probability value P i derived from the binomial test (Table 3).

Variable Influence of Tides on Seismicity From El Hierro Swarms (C1-C4 Phases)
The initial results from the statistical tests (Tables 1-3) and the analysis of the phase angle distribution ( Figure 5) suggest that the earthquakes belonging to clusters C1 and C4 are highly correlated with tidal confining stress, supporting the hypothesis that tidal stress may modulate the seismicity, especially in the case of events related to volcanic areas (Emter, 1997;Rydelek et al., 1988). The absence of correlation in the case of cluster C2 may be due to the fact that most of the seismicity belonging to that period of time (2011-10-17 to 2012-01-01) was characterized by a different mechanism than for C1 and C4. C2 phase seismicity has been proposed to respond to gravitational compaction of the magma plumbing system as magma withdrew (Martí et al., 2013), instead of being caused primarily by magma migration as in Phase 1 (cluster C1), or the posteruption Phase 4 (cluster C4), which was also associated with significant migration of magma González et al., 2013). The lack of correlation between tides and earthquakes during the period of cluster C2 could be due to the different conditions for tidal sensitivity at depth. Earthquake depths during the C2 phase where located deeper than C1 and slightly deeper than the C3 and C4 phases. Alternatively, the establishment of stable magma paths during the eruption (Oct 2011 to Feb 2012) could have altered the stress change conditions responsible for the seismicity modulation, as seen in other volcanoes (McNutt & Beavan, 1984).
Finally, the correlation between earthquakes and tidal confining stress in the events belonging to cluster C3 is less clear. The probability values corresponding to the Schuster test and the Monte Carlo permutations test are statistically significant based on the commonly used 0.05 threshold, even though they are higher (implying a weaker correlation) than the values obtained for clusters C1 or C4. But if we consider the value of the upper bound B of the Bayes factor B (Table 1), where the threshold for statistical significance is B > 16, the result is no longer significant. Nevertheless, the earthquakes in cluster C3 occurred during the posteruption process we are calling Phase 3, where magma migration also occurred , so a weak correlation should not be ruled out.
It should be noted that there are differences in the frequency distributions obtained for the different clusters ( Figure 5). The peak of the sinusoidal curve lies around tidal stress phase angle ϕ = 288° in the histogram related to cluster C4 (Figure 5d). For the cluster C1 histogram it is centered at ϕ = 251°, while the histogram for cluster C3 is around ϕ = 209°, closer to the minimum tidal stress, although still in the phase angle range where tidal stress is increasing (Figure 5c). These differences may be due to the different types of focal mechanisms controlling the seismicity in the clusters and/or also related to the depth of the earthquakes, as controlled by magma location.

Tidal Influence on Seismicity During the Beginning of Eruption Period (C1 Subphases)
To analysis the influence of tides around the submarine eruption at the end of Phase 1, we have calculated tidal volume strain in the setting of the shallow reservoir (latitude 27.6653°, longitude -18.0370°, depth 4.5 km; see González et al., 2013) at 2-h intervals during the time period covered by Phase 1 (Figure 6). We have compared it to the number of earthquakes in cluster C1 recorded every 2 h during the same time period. Figure 6 is divided in two subperiods as before or after the beginning of the southward migration of earthquakes on September 15, 2011. During the first period (Figure 6a; 2011-07-15 to 2011-09-14) no clear correlation appears between the biggest earthquake swarms recorded and the highest (or lowest) tidal amplitudes in the fortnightly period. In fact, most of the swarms seem to occur when the upper envelope MIGUELSANZ ET AL. Note. N tot is the number of events of each subset. N i,obs is the number of events observed in the increasing-tide phase angle half cycle (180° < ϕ ≤ 360°). P i is the probability of getting a number of increasingtide phase angle events greater or equal to the observed number N i,obs , assuming the null hypothesis that earthquakes are not influenced by tidal stress. B is derived from P i according to expression 5. Occurrences where P i < 0.05 or B > 16 are in bold.

Results of the Binomial Test for the Confining Tidal Stress in the Seismic Crisis of El Hierro
of the tidal strain amplitudes graph is low (or the lower envelope is high), suggesting that the amplitude of stresses does not control seismicity rate. On the other hand, in the second period (Figure 6b; 2011-09-15 to 2011-10-15) it is obvious that the most significant episodes of seismic activity prior to the submarine eruption of October 10, 2011 began around September 26, 2011, when the tidal volume strain graph is close to its extreme values. This also is the time when magma was postulated to begin an upward migration toward the shallower magma reservoir (González et al., 2013), which suggests that tidal strain can modulate more strongly the seismicity associated with magma migration occurring at shallower depths. Similar figures covering Phase 2, Phase 3, and Phase 4 are provided as supporting information (Figures S1-S3). Figure 7 shows the same plot with tidal strain split into three components: north-south, east-west and vertical. The graph shows that east-west horizontal strain and vertical strain are equally important, whereas north-south strain is less relevant. Similar figures analyzing the three components of tidal strain during Phase 2, Phase 3, and Phase 4 are included in the supporting information ( Figures S4-S6). As seen in those figures, vertical strain is the most important component in the last three phases.
The results of the application of the statistical tests to subclusters C1A and C1B are summarized in Table 4   whereas the histograms of the frequency distributions of confining tidal stress phase angles are shown in Figure 8.
As for the cluster C1, both subclusters C1A and C1B show high correlation between confining tidal stress and the occurrence time of earthquakes (Table 4), but there is a significant difference between both distributions. The peak of the fitted curve lies around tidal stress phase angle ϕ = 279° (which approximately coincides with the timing of the maximum increasing stress change rates) in the histogram related to subcluster C1A (Figure 8a). For the subcluster C1B histogram, it lies around ϕ = 216°, with a large number of events concentrating on the tidal phase range 150° < ϕ < 210°, central to the distribution ( Figure 8b). As we consider tidal stress positive in extension, minimum tidal stress is indeed maximum tidal stress in compression. This suggests that, during the period of upward magma migration toward the shallower magma reservoir, tidal compressive stress varies the magma flow and pressure conditions and influence seismic activity. the earthquakes in subclusters C1A and C1B which supports this interpretation. These results, for El Hierro magmatic system, suggest that there might be different tidal modulations of the seismic events related to magma migration and in particular ascent, depending on the depth.

Depth Dependence for the Tidal Influence on Seismicity
In order to test the sensitivity of the correlation to depth, we have focused on tidal confining stress, and we have considered 15 km depth as a threshold, as the boundary between the crust and the mantle lies at 15-16 km in the vicinity of the Canary Islands (Watts, 1994). In subclusters C1A and C1B, where most of the events belong to the depth range 10-15 km, high correlation is found between confining tidal stress and the occurrence time of earthquakes (see Table 4). In cluster C2, where the vast majority of events are >15 km depth, there is no statistically significant correlation with tides, as seen in Tables 1 and 2. In cluster C3, where, again, the great majority of hypocenters are in depth values > 15 km, the correlation is clearly lower than in subcluster C1A or subcluster C1B (see Tables 1 and 2). Finally, in cluster C4 we have considered two subcluster: C4A (earthquakes with hypocenters in the depth range between 0 and 15 km, 526 events), and C4B (earthquakes with depth > 15 km, 1,556 events). Table 5 shows the results of the application of the statistical tests defined in Section 3.3 to both subclusters C4A and C4B. In all cases, subcluster C4A gives a high correlation with tidal confining stress that is better than the one obtained for the whole cluster C4 (as seen in Tables 1 and 2). Figure 8c shows the histogram of the frequency distribution of tidal confining stress phase angles for subcluster C4A. The most populated tidal phase range is 240° < ϕ ≤ 270°, and the peak of the fitted curve lies around tidal stress phase angle ϕ = 300°, a value quite close to the maximum of the fitted curve corresponding to the entire cluster C4 (ϕ = 288°, as seen in Figure 5d). On the other hand, there is no statistically significant correlation with tidal confining stress for subcluster C4B (Table 5). These results suggest that tidal stress correlation with the origin times of the earthquakes is better for events whose hypocenters are at crustal levels (<15 km).

Tidal Stress Versus Tidal Stressing Rate Dependence on Seismicity at El Hierro
Previous research indicated that the strongest gradients of tidal variations might be more important for the modulation of the seismicity than the absolute value of the tidal stress (McNutt & Beavan, 1981, 1984. Therefore, stressing-rates could be a major controlling factor, rather than the absolute maximum/minimum stress change values, to explain the correlation between tides and the occurrence of earthquakes at El Hierro. With that in mind, we calculate tidal confining stress rates corresponding to all earthquake clusters or subclusters C1A, C1B, C2, C3, C4A, and C4B, and attempt to determine if the events occur more frequently when tidal stress rate is higher. As in the previous case, we use the total tides (load tides plus body tides) in the calculation of the tidal phase angle ϕ. The corresponding histograms appear in Figure 9, and the results for the statistical tests are shown in Table 6. There are three clusters with statistically significant values for all the tests defined in Section 3.3, namely C1A (Figure 9a), C1B (Figure 9b), and C4A (Figure 9e). There are no significant correlations between earthquakes and tidal stress rates in the events belonging to clusters C2 (Figure 9c) and C4B (Figure 9f). Finally, the correlation for cluster C3 appears to be statistically MIGUELSANZ ET AL. Note. Ps is the probability that the phase angle distribution to be random, according to Schuster's test. B is derived from P s according to expression 5. For each permutation, a Schuster probability value P* is obtained. NSIG is the number of permutations where P*≤ P s for each cluster. SL means "significance level." The significance level MC sl of the permutation test is obtained according to expression 3 in Section 3.3. Occurrences where P s < 0.05, B > 16 or MC sl < 0.05 are in bold. The type of tidal stress considered is confining. significant when using the Schuster test or the Monte Carlo permutations test, but it is not significant when considering the upper bound B of the Bayes factor B (Table 6), so this correlation must be considered with caution. In addition, in most of these clusters (Figures 9a, 9c, 9d-9f), the most frequently occurring phase  Note. P s is the probability that the phase angle distribution to be random, according to Schuster's test. B is derived from P s according to expression 5. For each permutation, a Schuster probability value P* is obtained. NSIG is the number of permutations where P*≤ P s for each cluster. SL means "significance level." The significance level MC sl of the permutation test is obtained according to expression 3 in Section 3.3. Occurrences where P s < 0.05, B > 16 or MC sl < 0.05 are in bold. The type of tidal stress considered is confining. angle range is the high-tide half cycle (0° < ϕ ≤ 90° or 270° < ϕ ≤ 360°), where tidal confining stress rates are higher. The exception here is subcluster C1B (Figure 9b), where the highest frequency corresponds to phase angle range 240° < ϕ < 270°. It should be noted that, in the time period covered by subcluster C1B, our previous analysis of tidal strain and stress values suggests that tidal compression plays a key role in seismicity modulation, rather than tidal extension, and this could explain the difference in behavior of subcluster C1B.
Our analysis of tidal modulation of earthquakes in El Hierro based on tidal stress rate agrees with the results from our previous analysis based on tidal stress values (Tables 1, 2, 4, and 5), in the sense that the subsets of events where the statistical significance is poorer are clusters C2, C3, and subcluster C4B, which are mostly composed of earthquakes > 15 km depth, even though in the case of cluster C3 the results are still statistically significant for both Schuster and Monte Carlo tests. If we discard clusters C2 and subcluster C4B (as they obtain the worst results with all the statistical tests considered), there is an interesting similarity between the frequency distributions of clusters C1A and C4A, and between cluster C3 and subcluster C1B. Whereas the maximum of the fitted curve in the distributions of subclusters C1A (Figure 9a) and C4A (Figure 9e) lies around extreme values of the tidal stress rate cycle (ϕ = 356° and ϕ = 11°, respectively), it is centered close to ϕ = 290° in both frequency distributions for subcluster C1B (Figure 9b) and cluster C3 (Figure 9d). Although the correlation found for cluster C3 is not as consistent as the one obtained for subcluster C1B (it is not statistically significant when we apply the Bayesian alternative to p-values) these results may reflect that, in the case of cluster C3 and subcluster C1B, tidal compression plays an influence in magma flow and pressure during migration toward shallower depths. Finally, we have seen that both seismic events belonging to subcluster C1B and cluster C3 can be related to volcanic activity due to magma injection processes, even though these processes only ended with an eruption in the case of the cluster C1B activity.

Solid-Earth and Ocean-Loading Tidal Stress Influence on Seismicity at El Hierro
As suggested elsewhere (Tanaka, 2002;Xu et al., 2011), the statistical significance of the correlation between tidal stress and the occurrence time of earthquakes may vary according to the type of fault. Due to the lack of information about focal mechanisms in the available earthquake catalog, an analysis of tidal stress modulation according to faulting type in the El Hierro volcanic crisis is beyond the scope of this work, but should be considered if additional information becomes available. Nevertheless, an analysis of the correlations based on the three components of tidal stress (north-south, east-west, and vertical) can be found in the supporting information. Histograms of the three components for Clusters C1, C3, C4, C1A, C1B, C4A, and C4B (Figures S10-S16) are almost identical to those obtained for tidal confining stress, also shown for comparison purposes. There are no histograms for cluster C2 because there is no significant correlation for any of the components. Figures S17-S20, supporting information, show which tidal stress components are dominant in each of the Phases 1-4. As in the case of tidal strain, east-west horizontal tidal stress and vertical tidal stress are equally important in Phase 1, whereas vertical tidal stress is the most important during Phase 2, Phase 3, and Phase 4.
In all the time periods considered in this study, ocean-loading tides are stronger than solid Earth tides. According to Cochran et al. (2004), tidal stresses due to body tides can reach values up to 5 × 10 3 Pa, which is 1 order of magnitude less than the values that can be reached by tidal stresses induced by ocean-loading tides (5 × 10 4 Pa) in ocean basins. For all the events in the different clusters in the study, we have calculated half cycle amplitudes (see definition in Figure 4) induced by ocean tidal loading confining stress and half cycle amplitudes due to Earth tidal confining stress. The former is, on average, between 5 and 6 times greater than the latter, depending on the cluster (see data in the repositories zenodo.org and digital.csic.es). Therefore, the ocean tidal signal appears to be more important than the body tidal signal in the modulation of the seismic activity during the El Hierro volcanic crisis. This result could be of general application in many of the volcanic events in oceanic islands and should be confirmed with additional, complementary studies.
If we consider only ocean tidal stress in our study of correlations, the results are almost equal to those obtained using the total tidal stress (see Figures S21 and S22 in the supporting information). Figures S23-S26 in the supporting information show the horizontal component due to the Earth tidal stress only, which is not in-phase with the ocean tidal stress.
In order to show if the magnitude cutoff chosen for the catalogs is relevant in the results obtained, we have performed some of the statistical tests over the catalogs removing all earthquakes with M < 2. As seen in Table S3 in the supporting information, the results are quite similar to the ones obtained in the original analysis, where M = 1.5 was the lower magnitude cutoff. We have not performed the statistical tests over declustered catalogs because, as we have told in Section 3, this seismicity is highly clustered in such a way that replacing each cluster for the biggest or the first event in the cluster series results in catalogs with very few elements. An example is shown in Table S4 in the supporting information, where a declustering process using a cutoff of M = 2.5 is described. Table S5 in the supporting information provides a more detailed description of the declustered catalogs D1-D4 so produced.
We have tested also possible correlations between earthquakes and strong rainfall on the El Hierro island during this period. Meteorological data can be accessed via an open access network from the Spanish Meteorological Agency (AEMET) (http://www.aemet.es/es/datos_abiertos). Figure S8 in the supporting information shows the comparative between earthquakes and rainfall gauges during the period covered by cluster C1. Rainfall can infiltrate into the crust and increase pore-pressure, or cause crustal deformation by changing water storage and hydrological loading (see e.g., González et al., 2012;Johnson et al., 2017). For this analysis, only onshore events have been considered. Although the beginning of all the earthquake sequence appears to occur immediately after the rains occurred on July 20, 2011, there seems to be no other correlation along the series. It should be noted that there is only rainfall data at Hierro Airport Station located at the northeast of the island. There are no figures for clusters C2, C3, and C4 because the events belonging to those clusters are mostly offshore.
Correlation with north-south and east-west tidal tilt also was studied. As seen in Figures S27 and S28, supporting information, we see a repeat of the previous tidal stress analysis, in that the subsets of events where the statistical significance is poorer are clusters C2, C3, and subcluster C4B, even though in the case of cluster C3 the results are still statistically significant for both north-south and east-west tidal tilt components.

A Conceptual Model for the Influence of Tides on the Observed Seismicity at El Hierro
Our interpretation of the volcano dynamics that generated the spatiotemporal variability on tidal correlations is detailed in Figure 10. resulting in a relative decrease in compression and an increase in extension. In particular, correlation is highest at the timing with high values of extensional tidal confining stress rates. This may change magma flow and pressure conditions inside the magma-filled fractures and/or conduits associated to the magma migration from mantle source regions into the storage conditions at the CMR. On the other hand, during the period defined by subcluster C1B (2011-09-26 to 2011-10-10), magma ascended toward a shallower and smaller reservoir than the CMR, and briefly accumulated there until the eruption started on October 10, 2011 ( Figure 10b). As the lateral growth of this reservoir, which affects overpressure, is limited, conditions at the shallow reservoir may be more sensitive to smaller changes in the absolute tidal stress compression. At high tide, sea level rises and the compression exerted on the seafloor by the body of water increases, promoting magma ascent and the build-up of pressure within the magma reservoir. Tidal confining stress compression, together with tectonic and volcanic stresses, would encourage magma to migrate into zones of weakness to ascend, favoring seismic activity in the process. Note that Figure 8b suggests that it is tidal compressive stress which promotes the seismic activity related to this period.
There appears to be no correlation between tides and seismicity during the co-eruption phase, to which the earthquakes in cluster C2 belong (2011-10-17 to 2012-01-01). We attribute this to a change in the triggering mechanism of the observed seismicity, seismicity might not be due to migration of overpressurized magma (González et al., 2013), but due to other processes, such as the postulated gravitational compaction of sections of the magma plumbing system as magma withdrew (Martí et al., 2013). Another explanation could be that the rapid rise and higher flow of magma during the eruption could have generated stress changes stronger than the tidal stress changes, overwhelming that effect.
For the correlation between tidal confining stress and earthquakes in cluster C3, a mechanism based on tidal compression similar to that described in Figure 10a is a possible explanation, but at, probably, lower magma overpressure and involving smaller volume of magma. This magma injection and migration episode did not culminate in an eruption (Blanco et al., 2015). Nevertheless, more research is needed in order to clarify the exact physical mechanisms controlling the correlation because the reported hypocenters in cluster C3 are at deeper depths than those ones of subclusters C1A or C1B. In addition, the correlation found in cluster C3 is poorer than those obtained for subclusters C1A or C1B, and it is not statistically significant when we estimate the upper bound of the Bayes factor (Tables 1 and 6).
Finally, the hypocenter patterns of the events in cluster C4 (see Figure 3d) suggest important magma lateral migration. In this context, we again invoke the mechanism described for Figure 10a, based on extensional tidal confining stress rates controlling changes in magma flow and pressure conditions during magma migration to explain the observed correlation with seismic activity (high for cluster C4 and the subcluster C4A). The lack of correlation in the case of subcluster C4B (even though its hypocenter distribution is similar to the original cluster C4, as seen in Figure S9 in the supporting information) may be due to the greater depth of the events when compared to those in subcluster C4A, or perhaps due to reduced overpressures and flow rates in the conduits at the end of the magma injection and seismic swarm episode.

Conclusions
We have studied the correlation between seismicity and tidal stress during different stages of the 2011-2013 El Hierro volcanic unrest, assigning a tidal stress phase angle to each located earthquake. Our work shows that ocean and body tides modulated the seismic activity during specific phases of the 2011-2013 El Hierro unrest. We found that mainly vertical and E-W horizontal tidal stress values and their rates were correlated with earthquake occurrence times. We determined that the influence of ocean-loading tides is stronger than the influence of solid Earth tides and that tidal modulation seems to be dependent on depth. In addition, results obtained prior to the submarine eruption suggest that tidal compressive stress might modulate the magma migration characteristics, and hence magma static elastic stress changes causing the observed earthquakes. We conclude that overpressure and flow rates in the reservoirs and conduits were modified by the changes in tidal stressing rate, and absolute tidal stresses at shallower depth. There is no correlation between tides and seismicity at times with no observed surface displacements, which strongly suggest that tidal modulation might be related to overpressure during migration of magma.
Our analysis supports the hypothesis that tides may modulate earthquake activity in volcanic areas, particularly at shallow depths, which generally agree with previous studies (e.g., Berrino & Corrado, 1991;Bhatnagar et al., 2016;Dzurisin, 1980;Kasahara et al., 2001;McNutt & Beavan, 1981, 1984Rydelek et al., 1988;Scholz et al., 2019;Sottili et al., 2007). Our results should be considered in the study of future volcanic reactivation at similar ocean island volcanoes, where ocean-loading might be the dominate component of the seismicity modulation. Considering that the seismic monitoring of active volcanoes is one of the primary methods for characterizing activity levels and assessing associated hazard probabilities, insights into the impact of tidal forces on the seismic activity at active volcanoes has the potential to improve our understanding of volcano dynamics and better quantify eruption hazard.

Data Availability Statement
Used seismic catalog can be downloaded from IGN website, www.ign.es. The remainder of the data are available in the repositories zenodo.org (https://zenodo.org/record/3714198#.XnFmI6hKg2w) and digital. csic.es.