CC-BY 4.0
1. Introduction
Mayotte Island is one of the four islands of the Comoros volcanic archipelago. It is located in the Indian Ocean in the Mozambique Channel between Madagascar and Africa. Mayotte Island shows marked volcanic geomorphology. Volcanism in Mayotte started about 10 to 15 My ago [Audru et al. 2010]. This behavior continued during the Quaternary, with the last volcanic eruption occurring 7000 years ago [Zinke et al. 2003]. The source of volcanic activity in Mayotte is still debated. Emerick and Duncan [1982] suggests that the origin of the archipelago is a hotspot, while Nougier et al. [1986] postulates that the volcanism corresponds to the reactivation of old and deep lithospheric fractures. Michon [2016] also rejects the idea of a hotspot and proposes that the Comoros archipelago volcanic activity can be explained by lithospheric deformation related to the southern extension of the East-African rift.
In general, the archipelago of Comoros is considered a moderately seismic region. However, since May 10, 2018, unusual intense seismicity has been observed in the east of Mayotte Island. From May 10, 2018 to July 31, 2019 almost 2000 events with local magnitude ML⩾3.5, were recorded [REVOSIMA/IPGP 2019]. The largest earthquake occurred on May 15, 2018, with a magnitude of Mw = 5.9. After July 2018, the number of earthquakes decreased, showing less than a hundred earthquakes with magnitude ML⩾3.5 per month [REVOSIMA/IPGP 2019; Saurel et al. 2022]. Although the most significant earthquakes occurred at the beginning of the crisis, the region remains active in 2021 with 141 MW⩾1 Volcano-Tectonic (VT) earthquakes located in December 2021 [REVOSIMA/IPGP 2021]. Geodetic data recorded in Mayotte show transient displacements of approximately 17–20 cm to the east and subsidence of 8 to 15 cm from the beginning of the crisis until the end of July 2019 [Figures 5 and 6 of the Bulletin number 1; REVOSIMA/IPGP 2019]. Very long-period events (VLP) with a dominant period of ∼15 s have also been reported [Cesca et al. 2020; Lemoine et al. 2020]. An oceanographic campaign in the area led to the discovery of a new submarine volcano in May 2019, located approximately 50 km east of Mayotte, which forms a seamount of approximately 820 m height on the seafloor [Feuillet 2019]. This campaign and later oceanographic campaigns have revealed several lava flows along with acoustic plumes in the water column. The new volcano is located at the tip of a 50-km-long ridge composed of many recent volcanic edifices, interpreted as an extensional structure within an east–west striking zone between East Africa and Madagascar [Feuillet et al. 2021; Famin et al. 2020]. Assuming that the new volcano edifice began to be built in July 2018, the mean lava flow rate is approximately 180 m3⋅s−1 (according to a survey conducted 11 months after onset of eruption [Feuillet et al. 2021]). The extruded volume of lava could be as large as 5 km3 transforming it into the largest active submarine eruption ever documented [Feuillet et al. 2021].
Following the deployment of additional seismological stations both onland and offshore, two separate seismicity clusters are now identified [Saurel et al. 2022]. Most of the large MW⩾5.0 earthquakes that occurred at the beginning of the crisis were located in a distal cluster connecting deep parts of the plumbing system (at depth of 30–40 km) to the new volcanic edifice on the seafloor [Cesca et al. 2020; Feuillet et al. 2021; Saurel et al. 2022]. This initial seismic sequence was most likely associated with the migration of magma toward the surface. Later, the magnitude of earthquakes decreased progressively and the seismic activity migrated toward Mayotte island. This proximal cluster, located 10 km to the east of Mayotte, has a cylindrical shape that likely corresponds to a ring fault system below an ancient caldera structure [Feuillet et al. 2021].
In this study, we analyze the source of the VT earthquakes to constrain the spatio-temporal properties of magma migration at the onset of Mayotte volcano-seismic activity. The unusual increase of earthquake activity along with the birth of the large volcanic edifice offshore of Mayotte clearly indicate the magmatic origin of the seismic sequence. However, the relationship between the observed large (Mw⩾5.0) earthquakes and the migrating magma pressure source is still elusive. In particular, several VT events in Global CMT [GCMT, Dziewonski et al. 1981; Ekström et al. 2012] and Cesca et al. [2020] catalogs present large non-double-couple (non-DC) components. These non-DC events could potentially correspond to volcanic source processes such as dike opening due to magma migration, resonance of fluid-filled cracks or even complex ruptures on ring-faults [Ekström 1994; Chouet and Matoza 2013; Rodríguez-Cardozo et al. 2021]. These non-DC components could also be associated with complex ruptures involving multiple tectonic subfaults or have a spurious origin and be induced by measurement errors (e.g., noise) or modeling uncertainties (e.g., inaccuracy of the Earth model). Here we focus on trying to understand to what extent these earthquakes correspond to opening/closing dykes or whether they correspond to shear faulting caused by stress transfer caused by a migrating magmatic pressure source.
With this purpose, we use long-period waveforms to invert for the source parameters of Mw⩾5 earthquakes that essentially occurred during the first months of the Mayotte seismic crisis (May and June 2018). Unfortunately, few regional stations were available during this period. Hence, the inversion is based on long-period seismological observations, including body and surface waves. To mitigate the impact of lateral heterogeneities that can be particularly large for fundamental mode surface waves, we estimate the Centroid Moment Tensor (CMT) parameters using Green’s functions computed for a global 3D Earth model combining S40RTS and CRUST2.0 [Ritsema et al. 2011; Bassin et al. 2000]. We compare our solutions with other studies and discuss the magma transfer at the beginning of the seismo-volcanic crisis.
2. Data and methods
We select events from the GCMT catalog located in the vicinity of Mayotte Island in 2018. This selection gives a total of 27 events with Mw⩾4.8, occurring mostly at the beginning of the crisis (i.e. in May–June 2018). We use teleseismic and regional waveforms from stations with an epicentral distance from 0° to 90°. We use a combination of 150 channels from GEOFON [IS, GE, TT; GEOFON Data Centre 1993], GEOSCOPE [G; Institut de Physique du Globe de Paris (IPGP) and Ecole et Observatoire des Sciences de la Terre de Strasbourg (EOST) 1982], ASL/USGS [AU, IU, IC, GT; Albuquerque Seismological Laboratory (ASL)/USGS 1988, 1992, 1993], IRIS/IDA [II; Scripps Institution of Oceanography 1986], MedNet [MN; MedNet Project Partner Institutions 1990], RESIF [FR, RD; RESIF 1962, 2018], NARS [NR; Utrecht University (UU Netherlands) 1983], Karthala [KA; Institut de Physique du Globe de Paris (IPGP) 2006] as well as station GULU belonging to a temporary network installed in Uganda to investigate plume-lithosphere interactions [XW; Nyblade 2017]. Waveform data has been validated at long period by comparing with synthetics from large teleseismic earthquakes. We do not use the YTMZ strong motion station from the French RESIF-RAP network, which was the only station on Mayotte Island at the beginning, due to large errors in the internal clock during the studied time-period.
For each earthquake, we manually reject noisy records and data with dubious metadata. After data screening, moment tensor solutions are obtained using 23 channels in average. The time window used for the waveform inversion starts at the P-wave arrival and includes the first group of surface waves (R1 and L1). For specific stations, we manually change the time window duration to exclude clipped signals in the inversion. A causal bandpass Butterworth filter of order 4 is applied to the data. To mitigate long-period seismic noise, the filter corner periods are adapted as a function of the earthquake magnitude and the epicentral distance of the stations (see Tables 1 and 2). Using the resulting long-period waveforms, we then perform point source inversions to invert for CMT parameters (i.e., moment tensor elements, rupture duration, and centroid location in time and space). We use a modified version of the W-phase inversion method by Duputel and Rivera [2019]. This algorithm relies on a grid-search approach to find the point-source mechanism, time and location that minimizes the root mean square (RMS) waveform misfit. In this study, we assume three different levels of generality for the moment tensor: (1) Full moment tensor (FMT; i.e., six independent elements of the symmetric moment tensor), (2) Deviatoric moment tensor (DMT; i.e., assuming no net volume change during the rupture) and (3) Double-Couple (DC; i.e., assuming a pure shear rupture). While FMT and DMT inversions are linear for a given centroid location, DC solutions are estimated by grid-searching for strike, dip and rake angles corresponding to the minimum RMS misfit. For the sake of consistency, FMT and DC inversions are performed using the same dataset, time shift (Ts) and centroid location as those estimated for the DMT solution.
Bandpass filter used for records with epicentral distance larger than 4°
| Magnitude | Bandpass |
|---|---|
| 5.3 ⩽ Mw ⩽ 5.9 | 50–150 s |
| Mw < 5.3 | 50–100 s |
With the exception of earthquakes 2018-06-04 21:20, Mw = 5.1 and 2018-06-01 07:14, Mw = 4.7 for which we use 50–150 s and 30–80 s passbands respectively.
List of events for which an ad-hoc bandpass filter is applied for stations at epicentral distances smaller than 4°
| Event | Magnitude | Bandpass |
|---|---|---|
| 2018-05-15 15:48 | 5.9 | 50–150 s |
| 2018-06-12 17:17 | 5.4 | 50–150 s |
| 2018-06-25 17:41 | 5.3 | 50–150 s |
| 2018-06-27 06:40 | 5.2 | 50–100 s |
| 2018-06-18 13:43 | 5.1 | 50–100 s |
| 2018-06-23 19:45 | 5.0 | 50–100 s |
| 2018-06-05 23:02 | 5.1 | 30–80 s |
| 2018-05-30 05:54 | 5.1 | 30–80 s |
| 2018-05-25 09:32 | 5.0 | 30–80 s |
| 2018-06-07 13:06 | 5.0 | 30–50 s |
| 2018-06-10 13:04 | 4.8 | 30–50 s |
To account for large-scale 3D heterogeneities, we compute Green’s functions using SPECFEM3D_GLOBE [Komatitsch et al. 2015] for a global 3D Earth model composed of S40RTS [Ritsema et al. 2011] and CRUST2.0 [Bassin et al. 2000]. Green’s functions are calculated for various source locations with depth ranging from 5.0 to 35.0 km (each 2.5 km), latitude ranging from 13.03° S to 12.7° S and longitudes ranging from 45.48° E to 45.83° E (each 0.025°). The definition of this grid was guided by the earthquake centroid locations provided by the GCMT catalog and preliminary tests.
Once we obtain the moment tensor parameters we perform a decomposition of the moment into the double couple (DC), isotropic (ISO), and compensated linear vector dipole (CLVD). For this, we follow the definition of Vavryčuk [2015],
| (1) |
| (2) |
| (3) |
| (4) |
3. Earthquake migration uncovering a deep magma transfer in May–June 2018
The spatio-temporal distribution of earthquakes after CMT inversion is shown in Figure 1. Events are roughly aligned along a NW–SE structure with deeper earthquakes at the NW and shallower events at the SE (Figure 1a). Although our centroid locations are globally consistent with other catalogs (see comparison with GCMT and Cesca et al. [2020] in Supplementary Figure 1), there are some differences that result from the use of different Earth models, frequency bands and choice of seismological stations. For example, Cesca et al. [2020] use 1D velocity models and includes the local strong-motion station YTMZ (mitigating clock drifting issues by inverting for time-shifts at each station). In contrast, our study is based on a 3D global velocity model and does not include YTMZ. As shown in Figure 1b, the seismic events migrate upward along this NW–SE structure from a depth of ∼35 km toward the surface between mid-May and early June 2018. This seismicity migration is a robust feature that is clearly visible in multiple catalogs (GCMT, Cesca et al. [2020] and Lemoine et al. [2020]; cf., Supplementary Figure 1). These migrating earthquakes have been interpreted as the markers of a deep magma transfer at the origin of the new volcano evidenced by Feuillet [2019]. Magma likely originated from one or multiple reservoirs or mushes located at depth ranging from 25 and 60 km depending on the estimates [Lemoine et al. 2020; Cesca et al. 2020; Feuillet et al. 2021]. The absolute centroid depth of earthquakes might depends on the assumed local velocity structure, which is poorly known in the vicinity of Mayotte [Cesca et al. 2020; Dofal et al. 2021].
Spatial distribution of moment tensor solutions for the deviatoric inversion. (a) The (lower hemisphere) beachballs are positioned on the map at the estimated centroid locations. Colors correspond to the centroid depth. The red triangle indicates the position of the volcano [REVOSIMA/IPGP 2019]. The two side panels are NS and EW vertical cross sections including beachball visualized from the W and S respectively (back hemisphere beachballs). (b) Vertical A–B cross section. Back hemisphere beachballs (visualized from the SW) are colored according to the event date. Color scale starts on May 14, 2018 (red) and ends on June 27, 2018 (blue). The position of the volcano at the free surface is indicated with a red triangle.
Earthquake migration associated with magma transfer has been observed in many volcanic areas such as the Bároarbunga volcanic system in Iceland [Sigmundsson et al. 2015; Ágústsdóttir et al. 2016], the Kilauea volcanoes in Hawaii [Neal et al. 2019; Lengliné et al. 2021] and the Piton de la Fournaise volcano in La Réunion [Battaglia et al. 2005; Lengliné et al. 2016; Duputel et al. 2019]. However, as noted by Cesca et al. [2020] and Feuillet et al. [2021], deep crustal earthquakes associated with deep magma transfer is less common. Such seismicity at depths larger than 10 km has been previously observed at Lo’ihi and Kilauea volcanoes in Hawaii [Wolfe et al. 2003; Merz et al. 2019]. Deeper seismicity associated with volcano eruption has also been reported in La Palma [Torres-González et al. 2020], with event magnitudes ranging from 0.9 to 2.7 and a depth range of 12 to 33 km, which is comparable to what is observed offshore of Mayotte.
The occurrence of such an intense upward migrating seismicity along with the birth of a new volcano clearly indicate that the earthquakes are caused by a deep magma transfer offshore Mayotte. This sequence is striking for the occurrence of several events of relatively large magnitude (Mw⩾5), which is uncommon in an area formerly considered to be a moderately seismic region. To further analyze the relationship between earthquakes and the migrating magma, our moment tensor estimates are discussed in the next two sections.
4. Moment tensor solutions consistent with regional stresses
Moment tensor solutions together with centroid depth, station coverage, and the percentage of CLVD, ISO and DC components are listed in Figures 2 and 3. For each earthquake, deviatoric moment tensor (DMT), full moment tensor (FMT) and double couple (DC) solutions are depicted in red, gray and blue, respectively. DMT and FMT solutions are very similar for most earthquakes, with only minor differences. Some shallow events are associated with large non-DC components, hence with significant differences between DC and moment tensor solutions (i.e., FMT and DMT). This is evidenced in the Hudson diagram [Hudson et al. 1989] presented in Figure 4(a), showing that several earthquakes have a non-negligible CLVD and ISO components. The origin of such large non-DC components is discussed in Section 5.
Moment tensor solutions for the events analyzed in this study. In red, gray and blue are the beachball obtained from the deviatoric moment tensor (DMT), full moment tensor (FMT), and double couple (DC) inversion, respectively. For each earthquake we present the number of channels, the azimuthal gap of the station coverage, depth and the beachball along with the corresponding CLVD, ISO and DC percentage.
Moment tensor solutions for the events analyzed in this study. In red, gray and blue are the beachball obtained from the deviatoric moment tensor (DMT), full moment tensor (FMT), and double couple (DC) inversion, respectively. For each earthquake we present the number of channels, the azimuthal gap of the station coverage, depth and the beachball along with the corresponding CLVD, ISO and DC percentage.
Hudson diagram and tension and pressure axis for each earthquake. (a) Hudson diagram Hudson et al. [1989]. Blue and red indicate inversion solutions for FMT and DMT, respectively. (b) Pressure axes (P) are represented in blue. Tension (T) axes are represented in red. Circles and square indicate inversion solutions for FMT and DMT, respectively.
The dominant strike-slip regime is obvious: almost all earthquakes present a right-lateral strike-slip focal mechanism with a relatively consistent strike angle. This is also visible in Figure 4(b), with P and T axes showing NE–SW extension and NW–SE compression, which is similar to the GCMT catalog and Cesca et al. [2020]. Even earthquakes with a large non-DC component are consistent with this orientation of the principal axes. This observation is in good agreement with regional GNSS observations in the region of the Comoros archipelago [Stamps et al. 2018; Lemoine et al. 2020]. Using a combined stress inversion of focal mechanisms, deformation structures and intrusions, Famin et al. [2020] found that the Comoros archipelago experiences a dextral shear deformation with maximum compressive horizontal stress in the NW–SE direction and NE–SW extension. More locally, earthquakes offshore eastern Mayotte are located beneath a volcanic ridge that exhibits multiple NE–SW extensional features consistent with the direction of our T-axes [Feuillet et al. 2021].
VT earthquakes are generally considered to be brittle ruptures within the volcanic edifice triggered by stress perturbations induced by magma activity. However, different relationships have been observed between VT earthquakes and magma intrusions. In particular, Roman and Cashman [2006] showed that basaltic volcanoes are often associated with migrating earthquakes ahead of the dike tip with focal mechanism P-axes parallel to the regional maximum compressive stress [i.e., dike propagation model; Ukawa and Tsukahara 1996], while other sequences on strato-volcanoes depict random hypocenter distributions around inflating dikes with focal mechanisms presenting P-axes rotated ∼90° from the maximum compressive regional stress [i.e., dike inflation model; Roman 2005]. These differences are likely related to a number of factors including the existence of pre-existing faults in the vicinity of the dike, the regional stress-field and the properties of the ascending magma.
Observations at the onset of the Mayotte volcano-seismic crisis with migrating seismicity associated with focal mechanisms parallel to regional stresses are more consistent with the dike propagation model proposed by Ukawa and Tsukahara [1996]. The observation of basanitic pillow lava flows following the eruption [reported by Feuillet et al. 2021] thus confirms the idea that VT seismicity at Mayotte reflect stress changes induced by dike propagation. Moreover, given their large magnitudes (Mw ⩾5.0), the earthquakes observed at the onset of the sequence probably do not correspond to ruptures of previously intact rocks and rather occur on pre-existing faults that are already loaded by regional stresses. This agrees with stress calculations of Rubin and Gillard [1998], showing that VT events larger than magnitude 1.0 likely occur on pre-existing structures already close to failure. This is also confirmed by observations on volcanoes indicating that earthquakes triggered by dike migration usually occur on pre-existing fault systems [Gargani et al. 2006; Lengliné et al. 2016; Duputel et al. 2019].
5. Non-DC component of shallow earthquakes
As pointed out in the previous section, several VT earthquakes are associated with a large non-DC component during the Mayotte sequence (see Figures 2–4). As shown in Figure 5(a), most of earthquakes with large non-DC components (>50%) correspond to events shallower than 15 km. We do not observe any correlation of non-DC components with the event magnitudes. Similar non-DC components have also been reported by Cesca et al. [2020]. To assess the impact of such large non-DC component on waveform fits, we evaluate the ratio of RMS misfits computed for the FMT and DC solutions. Results shown in Figure 5(b) indicate a relatively moderate overall reduction of data misfit when inverting for the full moment tensor (i.e., FMT) compared to a DC parametrization. Here, we recall that FMT and DC solutions are inverted independently (i.e., DC parameters are not estimated by decomposing the FMT solution). Unsurprisingly, the reduction of data misfit increases as a function of the non-DC percentage, showing that a DC source cannot fit data as well as a FMT solution for earthquakes with a large non-DC component. This is illustrated in Figure 6 showing the waveform fit of FMT and DC solutions for the earthquake with the largest non-DC component (event ID 2018-06-05 T23:02 in Figure 2). FMT and DMT solutions being very similar for this earthquake, we only show synthetics of the FMT solution. The deterioration of waveform fit for the DC solution is visible at multiple stations (see ABPO, FOMA, LODK and ATD in Figure 6).
Percentage of non-DC components compared to the depth of the events and waveform misfit. (a) Percentage of non-DC as a function of the depth of the events obtained by a full moment tensor inversion. (b) The logarithm of the root means square (RMS) waveform misfit ratio between the solutions obtained with the full moment tensor and the double couple inversion. The colors indicate the days since 2018-05-14.
Waveform fit for event 2018-06-06 09:37. Observations (black) are compared with predictions created using the same centroid location but different moment tensor solutions. Gray and blue lines correspond to synthetics computed using the full moment tensor (FMT) and double couple (DC) solutions represented at the top of the figure, respectively.
Different factors could explain the large non-DC component of earthquakes: one possibility could be crack opening/closure related to volcanic activity (e.g., caldera collapse, dike inflation/deflation, perturbation of magma circulation, etc.). Another possibility could be related to mismodeling of a complex source with a single point source model. Large non-DC components related to dike inflation/deflation have been previously reported in Iceland [Hrubcová et al. 2021] for smaller magnitude VT earthquakes (1 < M⩽4). Rodríguez-Cardozo et al. [2021] reported 230 VT earthquakes from 3.7⩽Mw < 5.5 and depths between 1–8 km, some of them with considerable non-DC components in the collapse of the Bárðarbunga caldera, Iceland. In particular, they observed that the CLVD component increases with the magnitude, relating it to slip on a curved fault and to the caldera collapse.
However, given the moderate magnitude of earthquakes that we investigate, and the long-period waveforms used in the inversion, the effect of source complexity seems quite unlikely to be the origin of misfit. Finally, there is the possibility of a spurious origin (e.g., contamination by ambient noise or inaccuracies in the Earth model).
In the following, we explore the possibility that the non-DC events observed offshore Mayotte could correspond to the combination of strike-slip motion and dike opening/closing associated with magma propagation. From the equations provided by Dufumier and Rivera [1997, appendix A], we write the surface area of an opening crack as:
| (5) |
| (6) |
| (7) |
| (8) |
| (9) |
We then estimate the average magma propagation velocity for each event to evaluate if the non-DC component could be related to fluid transport. We consider dike widths (𝛥u) ranging from 0.5 m to 10 m. The former value corresponds to a width typically observed on basaltic volcanoes, and the latter is an extremely large value to get a lower bound estimate of the velocity [Rubin 1995]. For the FMT solutions obtained in this study, we estimate source duration assuming Tr = 2ts, where ts is the time-shift between the origin time and the centroid time. Centroid time-shift has been proven to provide reliable rupture duration estimates that are less affected by long tails in moment rate functions and by arbitrary modeling choices [Duputel et al. 2013; Meier et al. 2017]. For solutions provided by Cesca et al. [2020], we assume the scaling relation [Duputel et al. 2012] since no centroid times are inverted for this catalog. Figure 7 shows the magma migration velocities estimated using (7) with both the Cesca et al. [2020] catalog and our solutions. Estimated velocities range from 44 m/s to 1565 m/s and from 9.9 m/s to 350 m/s for dike openings of 0.5 m and 10 m, respectively. Even considering dike opening of 10 m, most non-DC events are associated with unreasonably large magma migration speeds greater than 100 m/s. For comparison, a rough estimate of the upward migration speed of earthquakes suggest that the dike propagated at an average velocity of 0.05 m/s between 2018/05/19 and 2018/06/08 (see Supplementary Figure 2). To remove any dependency on the dike width, we also estimate the magma flow rate in Figure 8 using (9). The estimated flow rates are globally much larger than the average flow rate of 180 m3⋅s−1 estimated by Feuillet et al. [2021].
Magma migration velocity of the fluid as a function of the non-DC percentage. Gray and blue circles represent the solutions obtained by Cesca et al. [2020] and the present study, respectively. Green lines correspond a propagation velocities for mantle-derived dikes 0.01 and 10 m/s [Rubin 1995]. We consider dike widths (𝛥u) of 0.5 m (a) and 10 m (b).
Magma flow rate as a function of the non-DC percentage. Gray and blue circles represent the solutions obtained by [Cesca et al. 2020] and the present study, respectively. Green lines correspond to the estimate magma flow rate of Mayotte eruption 180 m3/s [Feuillet et al. 2021].
Shallow moderate magnitude (M⩾5) events with significant non-DC components have been previously observed around active volcanoes. Such earthquakes have been associated with caldera collapses driven by pressure variations in magma reservoirs [e.g., Ekström 1994; Shuler et al. 2013; Duputel et al. 2019], migration of hydrothermal fluids, perturbation of magma and circulation of gases through shallow conduits [e.g., Chouet and Matoza 2013], resonance of the fluids [e.g., Maeda and Kumagai 2017] or complexities in the source [e.g., curved fault slip Rodríguez-Cardozo et al. 2021]. During caldera collapse events, we expect to have many focal mechanisms with vertical T or P axes [Shuler et al. 2013; Duputel and Rivera 2019]. Even if strike-slip earthquakes occasionally occur on ring fault systems [cf., Rodríguez-Cardozo et al. 2021], there is currently no indication of any caldera collapse at the location of the distal seismic swarm. There are indications of a resonating magma body (at the origin of VLP signals) with a structure possibly outlining a caldera, but this structure is located closer to Mayotte, in the vicinity of the so-called proximal swarm [Feuillet et al. 2021; Cesca et al. 2020]. Dike resonance also appears unlikely as these phenomena are more typically associated with very long-period signals with marked peak frequencies, while earthquakes studied here have a frequency content similar to regular VT earthquakes.
As mentioned above, large non-DC components can also be associated with spurious origin such as ambient noise or inaccuracies in the Earth model [Šílený 2009; Kumar et al. 2015]. This seems to be the case in Mayotte as DC solutions fit that data almost as well as FMT solutions (cf., Figures 5 and 6); suggesting that the studied earthquakes correspond to pure-shear ruptures. Moreover, as shown in Supplementary Figure 3, the non-DC part of moment tensors often differs between catalogs. In particular, Cesca et al. [2020] suggests positive isotropic components and negative CLVD for most events while our catalog is dominated by negative isotropic components and positive CLVD. Such variability among catalogs is consistent with bootstrapping results of Cesca et al. [2020], showing that the non-DC components are generally poorly resolved. Uncertainty in the moment tensor solutions increases at shallow depth [Morales-Yáñez et al. 2020; Fukao et al. 2018; Sandanbata et al. 2021]. The enhanced non-DC component of shallow earthquakes could be an artifact due to the mismodeling of the crustal structure in the source region. Depending on the velocity models, the Moho depth in the Mayotte area varies from 5 to 30 km [e.g., Coffin et al. 1986; Laske et al. 2013; Bassin et al. 2000; Pratt et al. 2017; Cesca et al. 2020]. While our 3D Earth model is based on CRUST2.0, which assumes an oceanic crust of thickness 12.5 km, Dofal et al. [2021] proposed recently that the region corresponds to a 17 km-thick continental crust overlying a magmatic underplated layer defining a second interface at a depth of ∼27 km.
6. Conclusion
Characterizing the source of earthquakes off the east coast of Mayotte is essential to understand the link between magma transport and seismic faulting in the region. Consistently with other catalogs, our CMT solutions indicate an upward migration of seismic events at the onset of the Mayotte seismo-volcanic crisis. Although the region remains seismically active as of January 2022, the largest earthquakes of the sequence with magnitude Mw⩾5.0 occurred during this initial phase in May–June 2018. Most earthquakes have a strike-slip mechanism in agreement with regional stresses in the Comoros archipelago. Our catalog along with solutions by Cesca et al. [2020] and GCMT include some shallow events associated with large non-DC components that could have been interpreted as opening or closing of sub-vertical dikes. However, our analysis indicate that these large non-DC components are probably not physical (i.e., not related to dike opening/closing) and are rather artifacts induced by mismodeling due to inaccuracies in the shallow velocity structure. The migration of earthquakes along with the consistency of focal mechanism with regional stresses suggest that seismic events are triggered by stress perturbations induced by the upward propagation of magma. Given their large magnitudes (i.e., up to Mw = 5.9), the triggered seismic events likely occur on pre-existing strike-slip faults already loaded by ambient stresses. This volcano-seismic sequence thus demonstrate how a magma intrusion can induce relatively large earthquakes even in a moderately seismic region.
Conflicts of interest
The authors declare no competing financial interest.
Dedication
The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.
Acknowledgments
We thank two anonymous reviewers and the editors for their useful comments during the revision process. Dr. Lise Retailleau comments also helped us to improve the manuscript. This work was supported by the Initiative d’Excellence (IDEX) funding framework (Université de Strasbourg). This project has also received funding from the European Research Council (ERC, under the European Union’s Horizon 2020 research and innovation program under grant agreement No. 805256) and from Agence Nationale de la Recherche (project ANR-17-ERC3-0010). CMY acknowledges funding from ANID/FONDECYT 3220307 and the “anillo precursor project” ANID PIA Anillo ACT192169. We wish to thank the network and station operators for their commitment to collect high-quality seismic data and the Computational Infrastructure for Geodynamics (http://geodynamics.org) which is funded by the National Science Foundation under awards EAR-0949446 and EAR-1550901.