Evidence of earthquake seasonality in the Azores Triple Junction

This work presents evidence of seasonal and inter-annual variations of the earthquake occurrence rate in the Azores Triple Junction, on the Mid-Atlantic Ridge (MAR). Annual cycles in microearthquakes are relatively common in intraplate continental regions affected by large hydrological loads, but this is the first time that earthquake seasonality is recognized near mid-ocean ridges. First, we benchmark the methodology by matching the published results of earthquake seasonality in the intraplate New Madrid Seismic Zone (USA). Next, we analyze the Azores earthquake catalogue, from 2008 to 2018, separately for oceanic and island regions. The results demonstrate that the seasonal modulation of the seismicity rate is only observed in the ocean, especially in the vicinity of the triple junction, with more earthquakes occurring during the summer months from May to August. Monte Carlo simulations show that the probability of observing such seasonality by chance is less than 1% for the magnitude band from 3.3 to 4.5, well above the detection threshold and magnitude of completeness of the seismic catalogue. The methodology includes a Jack-Knife approach, which shows that the oceanic seasonality is not the consequence of abnormal or extreme events. Although we speculate about possible earthquake triggering processes, it remains a challenge to definitely establish the mechanism responsible for the observed earthquake seasonal modulation in the Azores.


Introduction
Several seismic regions exhibit interesting periodicities in their earthquake rates. Observed periods ranging from hours to decades are especially well documented on continents and islands, where land seismic networks have a good detection capability and yield quite complete catalogues. In continental areas, seasonal periodicities have been associated with variations in water loads arising from the atmospheric and terrestrial components of the hydrological cycle, among other external forcings. Some of the best examples of this seasonal modulation show a good correlation between earthquake rate and amount of precipitation (Bollinger et al. 2007), thickness of snow and ice (Heki 2003), and surface water levels in reservoirs, lakes, and watersheds (Craig et al. 2017;Johnson et al. 2017).
In the ocean, where the earthquake detection capability is lower due to the sparsity of ocean seismic networks, periodic modulations of the seismicity are more challenging to observe, especially in locations where earthquake rates are low. Still, such modulations have been established at intermediate and fast spreading centers in the Pacific Ocean, for different phases of volcanic activity, using data acquired by temporary OBS (Ocean Bottom Seismometer) surveys (Wilcock 2001;Tolstoy et al. 2002;Stroup et al. 2007;Tan et al. 2018;Scholz et al. 2019). Less attention has been paid to slow spreading ridges, such as the Mid-Atlantic Ridge (MAR), where most studies (e.g., Grevemeyer et al. 2013;Horning et al. 2018) have not fully investigated triggering mechanisms and the temporal variability of seismicity.
Variations of earthquake rates with terrestrial and oceanic tides, especially semidiurnal, diurnal, and fortnightly cycles, have been well documented in previous works (Tolstoy et al. 2002;Stroup et al. 2007). Most of these studies have evidenced that at mid-oceanic ridges the number of earthquakes tends to increase during or towards low tides, when the ocean load decreases and normal faults are unclamped (Wilcock 2001;Tolstoy et al. 2002;Stroup et al. 2007; Leptokaropoulos et al. 2021). However, the seasonal and longer-term periodic modulation of oceanic microearthquakes has been so far hindered by the lack of continuous long-term high-resolution datasets. To the best of our knowledge, this study unveils for the first time a seasonalscale variability in the seismicity of an oceanic region.

Seismo-tectonic setting
The Azores archipelago consists of 9 islands roughly aligned along the Terceira rift (TR), an incipient spreading center located in the triple junction between the Eurasian, North American and Nubian plates (Fig. 1). The region forms a shallow, tectonically, and volcanically active plateau, shaped by the interaction between a mantle plume and divergent plate boundary forces (Neves et al. 2013). The triple junction itself is at present near 38.5°N in the intersection of the MAR with the Pico-Faial islands alignment (Luis and Neves 2006). The tectonic regime is transtensional, and the seismicity is concentrated in discrete clusters aligned with the islands (Fig. 1). The earthquakes are mainly shallow, with depths less than 20 km (Instituto Português Do Mar E Da Atmosfera, I.P (2006)). Based on regional and global earthquake catalogues, Custódio et al. (2016) analyzed focal mechanisms solutions in the Azores region and found that they can be grouped into 4 regional clusters (Fig. 1, inset). The pressure axes of all clusters are vertical, indicating horizontal extension and normal faulting, although individual earthquakes with strike-slip mechanisms also exist. Cluster A is the closest to the triple junction and the one showing the highest heterogeneity of focal mechanisms within the cluster. This is in agreement with morphotectonic interpretations indicating that in this transtensional domain, rifting processes in pure tension coexist with the development of NNW-SSE left lateral shear faults and WSW-ENE relay zones (Lourenço et al. 1998;Madeira et al. 2015). Notable earthquakes include the January 1, 1980, earthquake, which reached a magnitude Ms 7.0 (Borges et al. 2007) and caused more than seventy deaths in the Terceira, Faial, and Graciosa islands, and the more recent July 9, 1998, Faial earthquake with Ms 6 (Matias et al. 2007). The recognition of a seasonal earthquake modulation in this region demonstrates its sensitivity to external forcings, which may affect magmatic systems near critical state and faults on the verge of rupture, thus improving our knowledge of the physical mechanisms of earthquake generation and contributing to better earthquake hazard assessment.

Materials and methods
Our study uses the Azores seismic catalogue of the Instituto Português do Mar e da Atmosfera (IPMA), the national earthquake monitoring agency, which operates a seismic network in the archipelago (Fig. 1). We analyze the earthquake catalogue between January 1, 2008, and December 31, 2018, a recent period during which the seismic network and catalogue are fairly homogeneous. For this period, the catalogue contains a total of 13.492 events. We divide the earthquake catalogue into insular and oceanic regions since the earthquake triggering mechanisms are expected to differ in these two regions. For instance, Martini et al. (2009) detected a seasonal effect in the seismicity occurrence on the island of São Miguel between 2003 and 2004 that was highly correlated to rainfall episodes. However, rainfall can be ruled out as a triggering process in the ocean. The complete data set (including islands and ocean earthquakes) is mostly dominated by the oceanic domain, which comprises 87% of the events (11.719 earthquakes). In addition to the oceanic and island subsets of the catalogue, we also separate the events into winter (November to February, NDJF) and summer (May to August, MJJA) subsets. This group of months is the one that best matches the seasonal patterns of the 1981-2010 climatic normal (temperature and precipitation) of the Azores central and eastern group of islands (Instituto Português Do Mar E Da Atmosfera, I.P 2022).
In order to assess the variability of the earthquake rate in the Azores, we followed the methodology proposed by Craig et al. (2017). To be sure that the methodology was correctly applied, in addition to the Azores data, we also used the seismic catalogue of the New Madrid Seismic Zone (NMSZ) used by Craig et al. (2017) to reproduce their results. The NMSZ catalogue is available from the Center for Earthquake Research and Information (CERI), University of Memphis (http:// www. memph is. edu/ ceri/ seism ic/ index. php).
As a first step in the methodology, the seismic catalogues are first declustered using the Reasenberg (1985) method, in order to prevent foreshock/aftershock clusters from biasing the temporal distribution of seismicity. This process was done in ZMAP software (Wiemer 2001). The magnitude of completeness (Mc) is computed considering the Gutenberg-Richter magnitude-frequency distribution, for the oceanic and islands catalogue (Supplementary material Fig. 1). In the Azores archipelago, the Mc varies noticeably both in space and in time. When comparing land and ocean, the oceanic catalogue has an overall Mc = 2.0 and the insular catalogue has Mc = 1.4 (computed considering the entire 2008-2018 time span). Lower Mc values in the insular catalogue can be explained by the location of the seismic stations within the islands (Fig. 1). When the seasonal subsets are considered, the oceanic catalogue shows Mc = 1.9 in summer (MJJA average) and Mc = 2.6 in winter (NDJF average). Larger Mc values in winter correspond to reduced earthquake detectability caused by an increase in seismic noise produced by ocean storminess, associated with the passage of the North-Atlantic winter storms over the Azores region. When analyzing the spatial distribution of Mc, we found that maximum values (Mc = 2.6) occur for regional cluster C due to the largest distance to the stations (Mc per sub-region considering the entire time span are displayed in Fig. 4).
The second step is to run a Monte Carlo simulation to generate 10,000 random seismic catalogues with the same number of events and frequency-magnitude distribution as the real catalogue. The random catalogues mimic the observed catalogue in a-and b-value and magnitude of completeness, as well as in the range of longitude, latitude, magnitude, and depth (Bollinger et al. 2007;Craig et al. 2017). Next, each of the 10,000 random catalogues is divided into winter (November to February, NDJF) and summer (May to August, MJJA) months. Then, we calculated the ratio of the number of earthquakes in winter vs summer (ratio NDJF/MJJA), considering only earthquakes above a varying threshold magnitude. From these ratios, we calculated 99% and 95% confidence limits that were used to assess the significance of the observed NDJF/MJJA ratios.
The third step of the methodology involves a Jack-Knife approach to identify whether the observed seasonality could be due to abnormal seismicity in one of the years. We carried out the same analysis described above, but now removing one calendar year at a time from the original seismic catalogue. We obtained as residuals the difference between the observed NDJF/MJJA ratio and the 95% confidence interval. This approach allowed us to identify whether a specific year of observation biased the inferred seasonality (Craig et al. 2017).
In the last step of our methodology, we used singular spectrum analysis (SSA) to identify dominant periodicities in the seismicity rate. We computed a time-series of seismicity rate (number of earthquakes per month) for each one of the regions corresponding to the four clusters in Fig. 1 (cluster A to cluster D). The SSA is a form of frequency analysis applied to decompose the time-series into trend, periodic or quasi-periodic components and noise. This study follows the methodology described by Hanson et al. (2004) to extract the reconstructed components (RCs) which take the form of almost-sinusoidal oscillations, and which contribute the most to the total variance of the time-series. The calculations are carried out by diagonalizing a lagged covariance matrix (Vautard et al. 1992) using the US Geological Survey's Hydrologic and Climate Analysis Toolkit (Dickinson et al. 2014). Figure 2 shows the number of earthquakes in our study region recorded throughout time, binned in 1-month intervals, for both the island and oceanic declustered catalogues. A first visual inspection reveals that there is a cyclic variability in the occurrence rate of ocean earthquakes over the 11 years analyzed, which is not so clearly observed in the islands. The cumulative earthquakes rates increase steadily but show no jumps, indicating that the declustering process was effective in removing sudden bursts of earthquake activity. Figure 3 shows the results of the seasonality analysis for the Azores island and oceanic regions, displayed in the first and second columns, respectively, and the results for the NMSZ catalogue (third column), used to benchmark the methodology and only displayed here to demonstrate that our results perfectly match those published by Craig et al. (2017) (see their Fig. 4d, e). The histograms of the bimonthly number of earthquakes show that while the island catalogue displays more earthquakes in May and June, there is not a clear seasonal pattern in the earthquake frequency (Fig. 3a). The observed peak in May-June is actually due to a seismicity crisis that occurred in 2008 (see Fig. 2a), but if this crisis is removed, we still obtain no seasonality in the islands. In contrast, Fig. 3d shows that there is a clear tendency for more earthquakes in the summer in the oceanic region, peaking in July-August. We consider the results robust mainly for M3.3-M4.5, which are well above the detection limit and which are not affected by the summer/winter oscillations of the magnitude of completeness. For M > 5, the number of earthquakes is not sufficient to provide statistically significant results. Figure 3b and e show the winter/summer ratio (NDJF/ MJJA). The relevant results are those above the magnitude of completeness (gray shading marks the region below Mc). The figures also show the 95% and 99% confidence limits obtained from the Monte Carlo analysis. For the islands (3b), the ratio NDJF/MJJA is higher than one above Mc, meaning that more events happen in the winter. However, this seasonality was not confirmed to be genuine since the events lie within the 95% and 99% confidence limits, which means that a similar ratio may be observed by chance. In contrast, in the ocean, the ratio NDJF/MJJA remains lower than 1 for almost all magnitudes, indicating that more events occur in summer than in winter (Fig. 3e). The difference is significant above the 99% level for low-magnitude earthquakes (M2.0-2.5) but also for earthquakes with magnitudes M3.3-4.5. Only for earthquakes with magnitudes in between these two ranges (M2.5-3.3) the significance of the seasonality is not confirmed since it lies inside the confidence limits. While the seasonality of the earthquake rate in the lower magnitude band (M2.0-2.5) is likely related to seasonal variations in seismic noise and resulting detection capability, the seasonality of earthquakes with M3.3-4.5 cannot be explained by observational limitations because these events have larger magnitudes than the detection thresholds indicated by the Mc value (maximum Mc = 2.6 in winter).

Results
Bursts of earthquake activity may influence the seismic catalogue despite the declustering process, and so, the inferred seasonal trend may be influenced by a small number of anomalous events. The results of the Jack-Knife analysis conducted to exclude this possibility and hence further confirm that the seasonality of the ocean catalogue is genuine are presented in supplementary material (Fig. 2). Thus, the exclusion of each calendar year of data separately from the full data series still shows a statistically significant NDJF/ MJJA seasonal variation. The results of this Jack-Knife analysis are summarized in Fig. 3f, which displays the residual (c, f) Residuals of the Jack-Knife analysis showing the difference between the 95% confidence and the observed NDJF/MJJA ratios. Each black line corresponds to a Jack-knife run, obtained by removing one year of data from the analysis. The red line represents the residual using the whole data set between the calculated NDJF/MJJA ratio (red dots in 3e) and the 95% confidence limit. The Jack-Knife analysis for the islands (Supplementary material Fig. 3) confirms the non-existence of statistically significant seasonality as the residuals are always negative.
In order to characterize the spatial variability of the observed seasonality, we decided to further separate the oceanic catalogue into four subsets, coinciding with the 4 earthquake clusters shown in Fig. 1 (Custó dio et al. 2016). This allows us to investigate whether the detected seasonality in the earthquake frequency is spread out spatially or focused along specific sections of the Azores domain. Figure 4 shows the earthquake histograms and Monte Carlo simulations, computed as before, along with the results of the SSA for the seismicity rate (number of earthquakes per month) above Mc computed for each one of these subset regions (A to D). The SSA plots to the left show the seismicity rate above Mc and the first reconstructed component, which corresponds to an annual periodicity. The ratio NDJF/MJJA (middle column) shows that the seasonality is mainly observed in cluster A. According to the Monte Carlo simulations Fig. 4 Analysis of the seasonality for the four clusters in the study region (A to D). First and second columns are similar to Fig. 3. Third column shows the SSA of the seismicity rate for each one of the clusters (colors corresponding to Fig. 1 inset) with black lines showing the first reconstructed component corresponding to the annual periodicity. The percentage represents the contribution of the annual signal to the variance of the time series the seasonality in this cluster is significant above the 99% confidence level even for M > 3.0 events. Here, the annual component contributes to nearly 40% of the total variance of time series of the seismicity rate. Cluster A is also the one where the total number of earthquakes is the largest (4422 events) and the one closest to the Mid-Atlantic Ridge. The region of cluster A is indeed the most tectonically active and corresponds to the present location of the triple junction between the Eurasian, African, and North American plates, as evidenced by morpho-tectonic analysis, magnetic anomaly studies, and geodynamic models (Lourenço et al. 1998;Luis and Neves 2006;Neves et al. 2013).
The earthquakes in the ocean catalogue can also be classified according to epicentral distance (the distance to the closest station that recorded the earthquake). Figure 5 shows the spatial distribution of the winter (NDJF, on the top) and summer (MJJA, on the bottom) subsets of earthquakes considering magnitudes larger than 3. When analyzing the spatial distribution of the M > 3.0 events, the region of cluster A displays the largest difference between summer and winter (92 events in summer and 48 in winter), confirming the results of Fig. 4. Futhermore, most events in the region of cluster A are recorded at epicentral distances < 120 km, with similar spatial distributions in both summer and winter.

Discussion
The methodology followed in this study was developed for earthquakes recorded on land and its application to an oceanic region involves additional limitations. In the Azores, unlike the NMSZ, the ability to detect and locate earthquakes is conditioned by the seismic network configuration, with stations placed on the islands and aligned with the archipelago. The earthquakes in cluster A, in particular, are recorded at distances of less than approximately 120 km by the closest stations, but the range of epicentral distances listed in the seismic catalogue remains similar in summer and winter (Fig. 5). A more important question is whether the detected seasonality is simply the result of variations in the detection capability of the network as a result of varying seismic noise amplitude, due to ocean storminess. Especially in archipelago settings, such as the Azores, the seismic noise level increases considerably during the winter, producing a seasonality in the seismic catalogues that reflects not a natural process but rather an observational limitation. In fact, there is a significant inverse correlation between the significant wave height in the Azores region, which is directly related to seismic noise, and the seismicity rate (not shown here for the sake of conciseness). However, this issue should not affect earthquakes with magnitudes M3.3-4.5, which are well above the average M2.0 completeness magnitude of the ocean catalogue, as well as above the maximum Mc observed when assessing variations of Mc both in space and in time. The positive residuals that result from the Monte Carlo simulations of 10,000 random catalogues combined with the Jack-Knife analysis show that the chance of observing more events in the summer through random chance is less than 5% for the above forementioned magnitude band (Fig. 3f). In addition, the SSA, which does not assume any a priori knowledge of the seismicity time-series, also confirms the existence of a dominant annual periodicity that peaks during summer, especially in region A.
The analysis per region shows that the seasonality is not uniformly present across the whole domain. Regions B, C, and D have a relatively smaller number of events (2422, 597, and 1779, respectively) than regions A (4422). As the seasonal modulation is a process superimposed on the background seismicity, the modulation is expected to be more evident in regions having higher levels of background seismicity, namely regions A and B. The most seismically active regions are those where the faults are closer to rupture and therefore the most prone to be influenced by small external perturbations. Unfortunately, the information concerning the earthquakes (depth and magnitude) does not allow us to discriminate which kind of faults are prone to seasonality, making it difficult to interpret the results in the light of the local tectonics and driving stresses.
To fully understand the evidence presented in this paper, we need to consider possible drivers of oceanic earthquake seasonality. Seasonal hydrological loads commonly invoked to explain the seasonality of the seismicity in continental regions (e.g., rainfall, snow) are not plausible in the ocean. At mid-ocean ridges, enhanced seismicity has been observed during low tide, at both fast and slow spreading ridges (Tanaka et al. 2002;Tan et al. 2018;Leptokaropoulos et al. 2021). These observations have been explained by a decrease in the water overburden pressure, which causes either fault unclamping (Stroup et al. 2009;Wilcock 2009) or magma chamber inflation and subsequent Coulomb stresses that favor slip on steeply dipping normal faults (Scholz et al. 2019). The similar tectonic framework of mid-ocean ridges and the Azores, especially in the region of cluster A, supports the hypothesis that a similar mechanism comprising ocean water mass changes may operate at an annual scale.
Tides in the Azores have small amplitudes (< 1.5 m) and their annual and semi-annual components contribute to less than 5% of the total variance. Indicators of ocean water mass changes such as sea level anomaly (SLA, provided by CMEMS) and estimates of ocean bottom pressure (OBP) provided by global ocean circulation models (ECCO version 4 release 4, ECCO Consortium et al., 2021) and GRACE (Gravity Recovery and Climate Experiment, Wiese et al. 2018) satellite data show annual and semi-annual variations. Preliminary results show fairly good correlations between the seismicity rate and GRACE or OBP (Suplemmentary material Fig. 4). However, in terms of equivalent water thickness, the annual amplitude of tides (5-10 cm), SLA (5-15 cm), GRACE (5 cm), and OBP (10 cm) in the Azores is very small and unlikely to produce Coulomb stress changes capable of triggering earthquakes. Nevertheless, that possibility requires further study and therefore the possible modulating process in the oceanic environment remains an open research question.

Conclusion
This work is the first to present evidence of seasonal variations of oceanic seismicity rate, and we demonstrated it for a midocean ridge triple junction. In the most active regions of the triple junction more earthquakes with magnitudes M3.3-4.5 occur in the summer than in the winter. A robust statistical method involving Monte Carlo simulations and SSA analysis proves that the seasonality is genuine and cannot be attributed to seasonal variations of the seismic noise or variations of the magnitude of completeness of the seismic catalogue. Whether seasonality is observed in other settings near mid-ocean ridges remains to be investigated. Future studies modeling the solid Earth response to plausible earthquake triggering mechanisms involving ocean water mass changes need to be conducted in order to explain the observed modulation of the seismicity by seasonal oscillatory stresses. These studies will continue to improve our understanding of earthquake dynamics, which in turn will contribute to better monitoring of the seismicity in volcano-tectonic active systems such as the Azores.