The 2018–2019 seismo-volcanic crisis east of Mayotte, Comoros islands: seismicity and ground deformation markers of an exceptional submarine eruption


 On 10 May 2018, an unprecedented long and intense seismic crisis started offshore, east of Mayotte, the easternmost of the Comoros volcanic islands. The population felt hundreds of events. Over the course of 1 yr, 32 earthquakes with magnitude greater than 5 occurred, including the largest event ever recorded in the Comoros (Mw = 5.9 on 15 May 2018). Earthquakes are clustered in space and time. Unusual intense long lasting monochromatic very long period events were also registered. From early July 2018, Global Navigation Satellite System (GNSS) stations and Interferometric Synthetic Aperture Radar (InSAR) registered a large drift, testimony of a large offshore deflation. We describe the onset and the evolution of a large magmatic event thanks to the analysis of the seismicity from the initiation of the crisis through its first year, compared to the ground deformation observation (GNSS and InSAR) and modelling. We discriminate and characterize the initial fracturing phase, the phase of magma intrusion and dyke propagation from depth to the subsurface, and the eruptive phase that starts on 3 July 2018, around 50 d after the first seismic events. The eruption is not terminated 2 yr after its initiation, with the persistence of an unusual seismicity, whose pattern has been similar since summer 2018, including episodic very low frequency events presenting a harmonic oscillation with a period of ∼16 s. From July 2018, the whole Mayotte Island drifted eastward and downward at a slightly increasing rate until reaching a peak in late 2018. At the apex, the mean deformation rate was 224 mm yr−1 eastward and 186 mm yr−1 downward. During 2019, the deformation smoothly decreased and in January 2020, it was less than 20 per cent of its peak value. A deflation model of a magma reservoir buried in a homogenous half space fits well the data. The modelled reservoir is located 45 ± 5 km east of Mayotte, at a depth of 28 ± 3 km and the inferred magma extraction at the apex was ∼94 m3 s−1. The introduction of a small secondary source located beneath Mayotte Island at the same depth as the main one improves the fit by 20 per cent. While the rate of the main source drops by a factor of 5 during 2019, the rate of the secondary source remains stable. This might be a clue of the occurrence of relaxation at depth that may continue for some time after the end of the eruption. According to our model, the total volume extracted from the deep reservoir was ∼2.65 km3 in January 2020. This is the largest offshore volcanic event ever quantitatively documented. This seismo-volcanic crisis is consistent with the trans-tensional regime along Comoros archipelago.


I N T RO D U C T I O N
The volcanic Comoros archipelago ( Fig. 1) is composed of four main volcanic islands, from west to east: Grande Comore, Mohéli, Anjouan and Mayotte, and further east, the Geiser and Leven volcanic banks (Daniel et al. 1972). Volcanism spans between the Miocene and Holocene. In Mayotte Island, the earliest phase of magmatic activity is dated 11 Ma (Debeuf 2004;Pelleter et al. 2014), whereas the latest reported volcanic event occurred within a few 1000 yr in Petite Terre, East of Mayotte (Zinke et al. 2003a,b). At the present time, there is a diffuse and moderate seismicity in the archipelago and few historically felt events.
The Mayotte volcano-seismic sequence began on 10 May 2018. Over 1 yr, we recorded 32 earthquakes with magnitudes greater than 5, the largest one reaching M w = 5.9 on 15 May 2018, the largest event ever recorded in Comoros. Such intense and long seismic activity had never been recorded before in the archipelago. More than 1 yr after the onset of the crisis, and after ∼2 months of intense seismic unrest at the beginning of the sequence, the seismicity was persisting with a large number of small earthquakes ) and still ongoing ground deformations are observed by Global Navigation Satellite System (GNSS) and Interferometric Synthetic Aperture Radar (InSAR). There were also frequent episodic monochromatic very long period (VLP) earthquakes (Poli et al. 2019;Satriano et al. 2019;Cesca et al. 2020; in review) typical of active magmatic or hydrothermal processes as described in other cases (e.g. Chouet 2003;McNutt 2005;Shapiro et al. 2017).
In an active volcanic context, magma movement can cause seismic crisis and volcano-tectonic earthquakes migrations can lead to potential eruptions (e.g. Toda et al. 2002;Battaglia et al. 2005;López et al. 2012;Martí et al. 2013;Gudmundsson et al. 2014;Sigmundsson et al. 2015;Ágústsdóttir et al. 2016). There are several past examples of intense seismic activity associated with deformation in volcanic contexts. Such observations can be markers of magma intrusion especially through propagating and opening dykes, as in the Gulf of Aden, between November 2010 and March 2011 (Ahmed et al. 2016), where the spatio-temporal evolution of intense seismicity, with 29 M ≥ 5 earthquakes, was interpreted as the sign of magma ascent, followed by the propagation of dykes. Regarding the intensity of the largest events, in Izu islands (Japan) in 2000 (Toda et al. 2002), five M ≥ 6 events occurred within a few weeks associated with the eruptions of Miyakejima volcano and dyke intrusion.
During seismo-volcanic crises and magma intrusion into dykes processes, swarms of small to moderate earthquakes (M < 6) are commonly observed. As for example in Afar (Ethiopia) where a phenomenon that began in 2005 involved seismicity organized in clusters, including M > 5 earthquakes: between 2005 and 2010, 14 dyke intrusions were identified along a 60-km-long rift segment and earthquakes are organized in clusters, that is distributed at spatially reduced areas close to the dykes (e.g. Ayele et al. 2007;Grandin et al. 2011). Another example illustrating both seismicity and ground deformation patterns as markers of magmatic phenomenon is the rifting event near the Bárdarbunga volcanic system (Iceland) that lasted from 16 August to 6 September 2014. As revealed by seismicity and ground deformation, it involved a 45-km-long segmented dyke intrusion delineated by earthquakes, magma source deflation and slow collapse of a caldera, magnitude M > 5 earthquakes and effusive fissure eruption (Gudmundsson et al. 2014;Sigmundsson et al. 2015;Ágústsdóttir et al. 2016. Finally, the El Hierro eruption (Canary Islands) initiated on 10 October 2011 and had been preceded since July 2011 by intense and fluctuant seismicity (including acceleration and migration, not exceeding M L = 4.3) and by the concomitant evolution of deformation patterns. Pre-eruptive phases were distinguished, especially according to seismicity and displacement monitoring (López et al. 2012). Analysis of earthquake migration and surface deformation have shown that, after the initial fracturing and intrusion episode, a migration phase appeared and extended for ∼20 km, preceding a submarine eruption (Martí et al. 2013).
The 2018-2019 seismo-volcanic activity at Mayotte is unusual due to its duration and the intensity of its seismic sequence and deformations, marking the end of a volcanic quiescence period in the area. The crisis has been monitored since the beginning of the unrest with an evolving seismic network, originally not designed for accurate monitoring volcano seismicity. From 10 May 2018 to 15 May 2019, the seismic monitoring network remained sparse, but allowed for the analysis of the evolution of the seismicity pattern. Permanent GNSS stations and InSAR data registered different displacement patterns that we identified as different eruptive phases.
Here we present an analysis of the first year of the Mayotte 2018-2019 seismo-volcanic episode combining both seismological analysis and deformation observations (GNSS and InSAR) before the significant improvement of the initial seismic monitoring network in May 2019 operated through the installation of ocean bottom seismometers (OBS). We propose a model of scenario that explains, with different phases, the spatio-temporal evolution of the seismicity and the geodetic observations.

G E O L O G I C A L S E T T I N G A N D R E G I O N A L C O N T E X T
The Comoros Islands form an overall east-west (E-W) trending archipelago located north of the Mozambique Channel, between the eastern coast of Mozambique and the northern tip of Madagascar ( Fig. 1). North of the Comoros Archipelago, most authors agree that the crust of the Somali basin is oceanic dating to the Mesozoic time (Ségoufin & Patriat 1981;Rabinowitz et al. 1983;Malod et al. 1991;Sauter et al. 2018) from 152 to 120 Ma (Davis et al. 2016). South to the archipelago, the nature of the Comoros basin is debated between a thinned continental crust (e.g. Flower & Strong 1969;Bassias & Leclaire 1990;Roach et al. 2017) and an oceanic crust setting during the Jurassic Magnetic Quiet Zone (Talwani 1962;Recq 1982;Rabinowitz et al. 1983;Klimke et al. 2016). The present-day morphology of this area arose from the Permian-Triassic Karro northwest-southeast (NW-SE) rifting, which resulted in the separation between Gondwana and Madagascar continental blocks (e.g. Malod et al. 1991;Davis et al. 2016). The Mozambique and Somalia basins opened while Madagascar Island drifted southward along the Davie ridge (Fig. 1). The tectonic regime in the Mozambique Channel is dominated by overall east-northeast by west-southwest (ENE-WSW) extension, which is also identified with the East African Rift System and Madagascar (e.g. Bertil & Regnoult 1998;Piqué 1999;Delvaux & Barth 2010).
The Comoros archipelago began emplacement during the Cenozoic in the northern Mozambique Channel. Michon (2016) proposed that magmatic activity first appeared in Mayotte, ∼20 Ma and then in Anjouan, Mohéli, and Grande Comore ∼10 Ma. Comoros volcanism seems synchronous with volcanism in surrounding areas, that is in Madagascar and in some provinces linked to the East-African Rift System (e.g. Nougier et al. 1986;Debeuf 2004;Michon 2016). According to rock geochronology (e.g. Nougier et al. 1986;  (Blewitt et al. 2018). The values are found in Table S1. A clockwise correction of 2.1 • Myr −1 was applied to cancel the overall rotation around Mayotte. Main features in the area of the East African Rift are plotted (Davie Ridge, Eastern Rift, Main Ethiopian Rift, Western Rift, Malawi Rift and Assoua Fracture Zone). The plate boundaries in white lines are modified from Saria et al. (2014) and Stamps et al. (2018). Elevation grid is from Gebco 2014 (Weatherall et al. 2015). Inset: The Comoros islands surrounded by the Somali and Comoros basins. Bathymetry from Gebco 2014 completed around Mayotte by SHOM (2015) data. J. represents Les Jumelles: bathymetric edifices lying NE from Mayotte (visible on maps from Fig. 4). Debeuf 2004;Pelleter et al. 2014), the youngest volcanic activity is in Grande Comore (0.13 ± 0.02 Myr to present) with the active Karthala volcano, and the oldest was reported in Mayotte (∼11 Myr), with a similar age reported in Anjouan. In Mayotte, these authors identify different eruptive phases separated by quiescent periods. The northeast area of Mayotte is composed by a more recent volcanic complex dated from 500-150 kyr in Grande-Terre and up to 4 kyr in Petite-Terre (Nehlig et al. 2013;Pelleter et al. 2014). According to Zinke et al. 2003a,b), the latest volcanic events documented are some ash deposits in the barrier reef between 4 and 7 kyr. Offshore, on the insular slope of Mayotte, tens of small volcanic seamounts are observed in the north, northwest and northeast parts (Audru et al. 2006) and is mainly distributed off Petite Terre along a WNW-ESE ridge where a new active volcanic structure has recently been documented (Paquet et al. 2019;Feuillet et al. in rev.).
The origin of the Comoros volcanism is not well understood yet. It could be due to the presence of a hot spot, based on the westward migration of volcanism age (e.g. Emerick & Duncan 1982), or to the influence of lithospheric fractures (Nougier et al. 1986). Recent work (Debeuf 2004;Michon 2016) suggests a mixed solution; regional extensive tectonics interacting with asthenospheric processes.
A catalogue of instrumental seismicity of the Comoros was extracted for 1964-2018 from the International Seismological Centre catalogue (ISC 2016). For 1978-1995 was completed using the Malagasy catalogue (Bertil & Regnoult 1998), which lowers the detection threshold, depending on the distance from the Tananarive observatory (from m b ≈ 4.0 eastward from Mayotte to m b ≈ 4.3 around Grande Comore). The location accuracy is poor, up to 30 km for many events in this poorly instrumented region. The Mozambique Channel seismicity is marked by a northsouth (N-S) band of activity (Fig. 2), associated with the Davie ridge ( Fig. 1) showing normal faulting mechanisms (Grimison & Chen 1988). A zone of diffuse seismicity, ∼100 km wide, runs E-W from 48 • E and 42.5 • E along the volcanic line of the Comoros islands (Fig. 2). From 1982 to 2016, eight moderate earthquakes occurred (M w 5.0-5.3), spread throughout this band presenting normal faulting and strike-slip earthquakes, especially east of Mayotte (Fig. 2). In Madagascar, a moderate and distributed seismicity is associated with inherited structures and an overall E-W extension (Bertil & Regnoult 1998;Rindraharisaona et al. 2013).
The overall actual trend of E-W extension over the East African region is illustrated by seismicity pattern (Fig. 2), geodetic observations, and the expression of extensive tectonics from EARS to Madagascar, including Mozambique Channel and Comoros volcanic axis (e.g. Bertil & Regnoult 1998;Delvaux & Barth 2010;Rindraharisaona et al. 2013;Stamps et al. 2018). Deville et al. (2018) reported offshore extensional and trans-tensional fault systems, respectively, around the Davie Ridge and South of it in the Mozambique Channel. They linked these active fault zones associated with volcanic activity with the southernmost part of the eastern branch of the EARS. Along the Comoros axis, the few available focal mechanisms show normal faulting and strike slip, compatible in terms of orientation to a NE-SW tensional axis. From onshore field observations on Mayotte, Anjouan and Mohéli islands ( Fig. 1), Famin et al. (2020) proposed that the Comores archipelago delineates the northern boundary between Somalia and Lwandle plates, undergoing strike-slip stress field associated with a NW-SE shortening (where reverse faulting, strike slip and enéchelon extensional fractures coexist). Whereas, from offshore data, Feuillet et al. (in rev.) suggested that this diffuse E-W striking area accommodates the deformation between offshore East African rifts system and Madagascar grabens, thanks to major extensional features in a trans-tensional area.
The NE-SW extension across the line formed by the Comoros islands is also revealed by GNSS data. Fig. 1 shows the velocities of the GNSS stations in the area (Table S1) Lemoine et al. is bending,both GNSS (Fig. 1) and focal mechanisms (Fig. 2) confirm the presence of a NE-SW tensional axis. The GNSS station distribution is uneven and particularly sparse around Comoros. The relative velocity of SEYG (Seychelles) with respect to the stations PMBA and NMPL located in Mozambique, on the Rovuma Plate, is of 2.6 ± 0.5 mm yr −1 in the azimuth N50 • . As the Comoros axis is associated with volcanism since 20 Ma and seismically active, while there is no identified seismicity within the Somalia basin and between Seychelles and the Comoros, a significant part of this NE-SW extension between Seychelles and Mozambique could well be accommodated across the volcanic line of the Comoros, through oblique extension along the volcanic axis.
This relative motion possibly accommodated in the Comoros would fit with the focal mechanisms reported in Fig. 2 along the Comoros archipelago: normal faulting and strike-slip events, W and SSE of the Comoros that could indicate a trans-tensional stress regime.
Furthermore, VHMR and MTVE GNSS stations (Fig. 1), respectively, at the northern tip of Madagascar and in the Tanzania coastal area, point in opposite directions (relatively to Mayotte station MAYG), which suggest an overall and significant counterclockwise rotation in the area. South of the Davie ridge is dominated by left lateral shear at a rate of 1 mm yr −1 , according to the difference of velocities between PMBA and NMPL on the Rovuma Plate and ABPO, VOIM and ZMBT on the Lwandle Plate ( Fig. 1).
Of note, between the Neogene and Quaternary, tectonic and volcanic activities occurred in several branches in EARS (e.g. Déprez et al. 2013), but also in several part of Madagascar, especially in its northern part (e.g. Rindraharisaona et al. 2013;Michon 2016) and in Comoros (Nougier et al. 1986;Debeuf 2004;Michon 2016). Michon (2016) proposed that the Comoros archipelago could potentially coincide with the eastern border of the EARS. Debeuf (2004) discussed the influence of extension linked to the East African Rift on Comoros volcanism. She mentioned the presence of an ancient regional lithospheric structure: the Assoua fracture zone (east Africa), a Precambrian lineament of lithospheric dimension that, despite its age, could have localized the ascent of asthenospheric plumes in the recent geological past.

Historical seismicity
There are only a few testimonies of historical seismicity in the Comoros before the 20th century. Hachim (2004) identified damaging earthquakes in Mayotte in 1606, 1679 and 1788 from oral transmission without being able to constrain testimonies from bibliographical data. Gevrey (1870) and Vienne (1900) quote felt earthquakes in the Comoros Archipelago in 1808, 1829 and 1865, without any detail. Moreover, Lambert (1997) associated the 1829 event not with an earthquake but with a cyclone. From 1910 to 1960, the Malagasy academy published a catalogue of locally felt earthquakes made by the Tananarive observatory. The 16 January 1936 earthquake, west of Mayotte, was felt across the whole Comoros, causing moderate damage in several municipalities (La Dépêche de Madagascar, 19/02/1936, BNF). More recently, regional felt earthquakes have been reported in Sisfrance Océan Indien (2010). No historical chronicle documents any seismic cluster in Comoros.

Instrumental seismicity
Regarding the instrumental period, the catalogue of seismicity we built for the Comoros (see preceding chapter) contains some events within 150 km near the Mayotte edifice. Some events occurred south of Mayotte, including a m b = 5.0 event, on 23 April 1993, located 40 km south of the island. In 2007, at ∼150 km east of Mayotte, a patch of seismicity seems to be located around Zélée bank, including a magnitude 5.1 event on 23 June 2007 showing a strike-slip fault mechanism processed by the Global Centroid Moment Tensor (G-CMT, Dziewonski et al. 1981;Ekström et al. 2012). The most significant event close to Mayotte before 2018 occurred on 1 December 1993, located 40 km west of the island with m b = 5.2 and associated with moderate damage (Lambert 1997). Close to this event, a normal faulting earthquake occurred on 9 September 2011 (M w = 5.0). No significant event was reported near the 2018-2019 seismic sequence during the instrumental period, in a poorly instrumented region associated with high detection magnitude threshold (∼4.5). Moreover, no long lasting seismic sequence is reported in the instrumental catalogue of the Comoros.

The 2018-2019 seismic sequence
At the beginning of the sequence, the closest broad-band stations from international seismic networks operating in this region were more than 650 km away (Fig. 2). In Mayotte, the accelerometric station, YTMZ from the French RESIF-RAP network (Réseau Accélérométrique Permanent, RESIF 1995) was active in real time since 2016. To monitor the crisis, a virtual network was set up, composed of YTMZ and regional broad-band stations. This network evolved in time with the addition of local and regional stations until mid-July 2018 (Supporting Information and Table S2).
Here we present an analysis of the swarm of volcano-tectonic earthquakes from the beginning of the crisis until 15 May 2019, when new onshore and offshore stations allow a good azimuthal coverage of the seismic sequence and a better analysis of seismicity distribution (Saurel et al. 2019;Feuillet et al. in rev.).
The YTMZ station allows finer monitoring of the crisis than global networks. Thanks to it, from 10 May 2018 to 15 May 2019, 1872 events with a local magnitude M L (Supporting Information) greater than or equal to 3.5 could be detected (161 and 32 of them corresponding, respectively, to magnitude greater than or equal to 4.5 and 5, Fig. 3). Only 1163 of these 1872 events could be located with crustal phase data for distance less than 1400 km (Figs 4 and 5).
The location configuration (described in Supporting Information) have been built and tested on earthquakes that occurred after the improvement of the seismic network, that is after the period analysed in this paper. It makes it possible to limit epicentre shifts with the seismicity pattern observed beyond May 2019 thanks to the seismic network improvement Saurel et al. 2019).
Accuracy of the location depends directly on the network configuration that evolved with time and on the magnitude (on azimuthal coverage of the network and on number of manually picked seismic phases available for location). Epicentral location errors range from 10 km in the most unfavourable configuration to 1 km in the most favourable. Depth could not be retrieved from hypo71 (Lee & Valdes 1985) for some events associated to a reduced number of picked phases, especially until latest network improvement (i.e. 15 July 2018). In such cases, location is made fixing the depth at 30 km (depth average in the best network configuration is around 32 km). 30 km is also in the range of depth values from G-CMT locations (Fig. 6) and permit to retrieve the seismicity pattern during the whole crisis (Saurel et al. 2019). In other cases, average vertical location error is estimated to be equivalent to horizontal error (i.e. few km). However, depth distribution depends on the velocity model which is far from being constrained around Mayotte and surrounding areas (e.g. Jacques et al. 2019). Finally, we could state that the locations including more than 10 manually picked seismic phases, are stable, especially horizontally whereas depth distribution should be considered with caution. The crisis started on the morning of 10 May 2018 with small events detected by the YTMZ accelerometer only. In the evening, a M L = 4.3 earthquake, not detected by the global networks, was the first event felt in Mayotte. From May 13 to 15, the activity increased gradually with four events with M L between 4.5 and 5 (Fig. 3). The main event, and largest ever recorded in the Comoros, with magnitude M L = 5.8 and M w = 5.9 (G-CMT), occurred on 15 May The first cluster (light to dark blue) is active mainly from 10 May 2018 to early July 2018, the second (yellow to light orange) from 26 June 2018, and the third (red to pink) from mid July 2018. Within each cluster, colour of epicentres evolves with time (days from the beginning of the crisis, that is 10 May 2018). Green points represent earthquakes registered during the observation period but out of the clustering seismicity (one M L ≥ 3.5 event is reported westward from Mayotte). The bathymetry is from the Homonim project (SHOM 2015) and Gebco 2014 (Weatherall et al. 2015). Bottom panel: seismicity plotted on a vertical cross section between points A and B, only resolved depths are plotted (otherwise, the depth is fixed at 30 km). Yellow triangles indicate the GNSS stations. Epicentres marked by a light gray circle and transparency were localized with less than 10 phases and considered as less relevant (map and cross section).
2018. Then the seismicity was relatively steady for the following ∼15 d with 10 to 30 M L > 3.5 events per day. Two strong events occurred on May 20 and 21 (M w = 5.5). Between June 1 and 7, the seismicity culminated with 52 to 87 M L ≥ 3.5 events daily, with the highest activity on 1 June 2018. After 7 June 2018, the seismicity decreased. It increased again between 19 and 27, June 2018 with several M L ≥ 5.0 events with a new peak of activity on 23 June 2018 (50 M L ≥ 3.5 events). The period between July 10 and August 26 was quiescent with a sharp drop in the number and magnitude of events (Fig. 3). On August 26, eight events with magnitude between M L = 4 and 4.5 occurred (27 M L ≥ 3.5 events), initiating a renewal of activity with two M L = 4.8 events occurring on October 16 and November 7, marking the last day associated with such a high level of activity (afterward, the maximal daily number of M L ≥ 3.5 events were 14). However, their contribution to the total moment release remains low compared to the first weeks of activity (Fig. 3). Since October 2018, a relatively steady level of active seismicity is observed that began to lower at the end of April 2019 (Fig. 3).
The cumulated seismic moment for the whole period corresponds to a M w = 6.5 event, with most of the seismic moment released before early July 2018. The most active period, from 1 to 7 June 2018, corresponds to a M w = 6.2 event, more than the M w = 5.9 of the largest earthquake of May 15, whereas the largest magnitude registered in this period is M L = 5.3. In comparison, the cumulative seismic moment from 1977 to 2017 for the entire Comoros (between 42.5 • E and 48 • E) represents less than M w = 6.0. Between 1 August 2018 and 15 May 2019, the cumulative seismic moment released corresponds to a M w = 6.1 event. After a quiescent phase from 13 July to 26 August 2018, there was a resurgence of seismic activity at the end of August with lower intensity and low contribution to the moment release (Fig. 3).
In order to follow the temporal evolution of the seismicity, we choose to keep S-P time at YTMZ as a reference, as it is the only parameter common to the entire seismic sequence that is independent to location uncertainties or network evolution.
The seismicity is characterized by a long lasting sequence of events localized in a reduced area. Even if the largest earthquake to have been recorded in the Comoros belongs to this seismic sequence, the seismicity cannot be described from a main shock followed by aftershocks. S-P time and magnitude distributions ( . Vertical lines mark principal steps in the seismic network improvements. Events marked by a light gray circle and transparency were localized with less than 10 phases and considered as less relevant. Before mid-July 2018, depth is particularly badly constrained (consequently fixed at 30 km when not retrievable): hypocentral evolution must be considered accordingly. the presence of different clusters. The first seismic cluster began on 10 May 2018 and is associated with an increase of S-P time until 1 July 2018. Two other clusters then appeared associated to distinct values of S-P time. These three main clusters are coherent in time ( Fig. 3) and space, as they are restricted to limited area (Fig. 4).
In the first one, active mainly from May 10 to beginning of July 2018, the events (blue dots) occurred in a ∼20 km diameter area, corresponding to S-P time at YTMZ from 5.2 to 8.1 seconds (i.e. east of YTMZ, ∼30 to 50 km considering hypocentral distance). This first cluster gathers all the largest events of the crisis. During May 2018, most events were concentrated in an around 12 km diameter area located 30 km eastward from the island after the very first earthquakes occurred closest to the island (Figs 4 and 5). The whole period shows a gradual increase of S-P time at YTMZ, especially from the beginning of June 2018 (Fig. 3). From 1 June to 1 July 2018, S-P time at YTMZ increased from 5.8 to 8.1 s, associated with a migration SE, with a total of 20 km shift for the period (Fig. 5). Depths could be determined only for few events during June 2018 (Fig. 5), but they are compatible with an upward migration of the seismicity. SE migration is confirmed by G-CMT locations for the largest events ( Fig. 6) associated in parallel with a depth change from 35 km to 12 km from early June 2018 (depth from G-CMT catalogue is decreasing from June 4, Fig. 6), which is beyond the error bars and confirms the hypothesis of upward migration. Depth migration is also confirmed by Cesca et al. (2020) from moment tensor inversions and analysis of teleseismic depth phases. We distinguished two phases in this first cluster: cluster 1a from May 10 to the beginning of June 2018, and cluster 1b covering early June to beginning of July 2018, that started at the beginning of the seismically very intense first week of June. Cluster 1b is associated with the increase and spreading of S-P time reported at YTMZ and preceded its lowering at initial values at the end of June 2018 (Fig. 3).
From the end of June 2018 (around June 26) to July 2018, S-P time retrieves lower values equivalent to the beginning of the first cluster (between 5.2 and 6.2 s); the second cluster (yellow dots) starts westward from the latest events of the first cluster ( Fig. 4), that is closer to YTMZ (hypocentral distance to YTMZ: ∼40-49 km). This second cluster is close to the initial position of the first cluster, before its migration (i.e. corresponding to cluster 1a). Clusters 1 and 2 overlap in time until the former almost vanishes around 8 July 2018. In July 2018, a third cluster starts (red dots), with S-P time at YTMZ from 3.3 to 5.1 s, overlapping in time with cluster 2, but more energetic (Fig. 3). Clusters 2 and 3 are spatially separated by an aseismic gap around longitude 45.49 • (Figs 4 and 5), and distinction between these two families of events is sharper considering this border than values of S-P time (Figs 3 and 5). After mid July 2018, there is no evidence of overall earthquake migration within clusters 30 A. Lemoine et al.  Table S3. Upward migration is observed from early June 2018 (4 June 2018 is marked by the upper arrow).
2 and 3 (Figs 4 and 5). However, 26 August 2018 marked the end of the period of seismicity quiescence and Fig. 5 illustrates the associated aligned hypocentres distribution (likely along NW-SE axis). This particular pattern marked the beginning of an intense seismicity within the cluster 3. Inside clusters 2 and 3, depth ranged mainly between 25 and 40 km with a sparse seismicity shallower, except on 26 August 2018, when the aligned hypocentres seem to be associated to depth ranging between 5 and 25 km.
The earthquakes of the first cluster seem to be shallower than those of clusters 2 and 3 (Figs 4 and 5); however, the depth distribution of earthquakes must be interpreted with caution, due to the lack of an ad hoc velocity model and to the unfavourable geometry of the network, especially between May and June 2018. As a whole, the three clusters depict an activity migrating westward (Fig. 4) during the onset of the crisis until July 2018 (whereas inside the first cluster, a SE migration is observed, confirmed by G-CMT locations for the largest events), after what a steady seismicity pattern is observed.
Westward from the Mayotte edifice, one M L > 3.5 earthquakes is reported in Fig. 4 (green circle): on 6 April 2019 (M L = 3.7, with 12 phases). It occurred close to the epicentres of both the 1 December 1993 damaging earthquake (Lambert 1997, m b 5.2) and 9 September 2011 event (M w 5, normal faulting, G-CMT).
The 32 largest events (M L ≥ 5.0) are listed in Table S3 with their G-CMT moment when it exists. All together they correspond to 45 per cent of the total moment release from 10 May 2018 to 15 May 2019. All occurred during the first part of the seismic sequence, between 14 May and 27 June 2018, except three events (28 March 2019, M L = 5.0; 5 April 2019, M L = 5.2; and 14 May 2019, M L = 5.1). Except one reverse faulting mechanisms for 14 May 2019 event belonging to cluster 3, most of them have strike slip focal mechanisms and belong to cluster 1. The G-CMT nodal planes are very similar for these events with averages and mean scatters of 352 ± 9 • , 63 ± 9 • and 5 ± 13 • for plane 1 and 260 ± 12 • , 85 ± 11 • and 152 ± 9 • for plane 2 (Fig. 6). The corresponding tension axis is oriented N34 ± 15 • , which corresponds to a most favoured faulting azimuth N304 • . Both nodal planes indicate strike slip on a steep fault, indicating, during May and June 2018, either right lateral shear along ∼E-W structures or left lateral shear along ∼N-S structures (Fig. 6). The former is more consistent with the relative movements of the plates from either side of the archipelago, that is counterclockwise rotation of the Lwandle Plate with respect to the Somalia Plate (Saria et al. 2014), but potential pre-existing active structures that could be responsible for those mechanisms are still not identified (if we consider an influence of the regional tectonics on the Downloaded from https://academic.oup.com/gji/article-abstract/223/1/22/5850758 by Ifremer, Bibliothèque La Pérouse user on 30 July 2020 The Mayotte 2018-2019 seismo-volcanic crisis 31 seismo-volcanic crisis). N-S structures would have the same orientation as ancient transform faults linked to ancient spreading in the Somali basin along the N-S direction (e.g. Rabinowitz et al. 1983). In the surrounding areas, there are other strike-slip focal mechanisms reported in the G-CMT catalogue since 1976 (Fig. 2), making the whole consistent with a present-day trans-tensional regime in the area and a ∼NE-SW regional extension. Moreover, strike-slip focal mechanisms are commonly observed during dyke intrusions (e.g. Toda et al. 2002). 14 May 2019 occurred within cluster 3, at 34 km depth, its reverse focal mechanism is incoherent with regional trans-tensional stress regime and cannot be explained by dyke intrusion. As suggested by Cesca et al. (2020) who reported more events of this type located in our Cluster 3, such thrust mechanisms should be linked to local processes, probably of magmatic origin.

V E RY L O W F R E Q U E N C Y M O N O C H RO M AT I C S E I S M I C E V E N T S
During the crisis, long lasting monochromatic very long period (VLP) events were recorded by both local and international networks (Table S5 describes few of them). After the initial emergent onset of the signal, they show decaying harmonic oscillations comparable to Long-Period (LP) waveform patterns observed, for example, at Kusatsu-Shirane or Galeras volcanoes (Kumagai & Chouet 1999). However, they presented much larger period of oscillations (∼16 s) and long lasting signals. Systematic detection and analysis of such events are reported by Poli et al. (2019), Satriano et al. (2019) or Cesca et al. (2020) who reported more than 400 of such VLP events and located them in the vicinity of the cluster 3. The largest one, on 11 November 2018, well recorded worldwide, awoke the interest of the seismological community and the media (e.g. Wei-Haas 2018; Lacassin et al. 2020). This long lasting monochromatic VLP event was not preceded by a strong earthquake, but there are earthquakes embedded in the signal, during the occurrence of the VLP event, especially at the beginning of the signal (Figs 7 and S2). As observed in Polynesia in 2011 and 2013 by Talandier et al. (2016), such VLP wave train are associated with elliptical particle motion in the vertical plane, in the direction of propagation, between the source and the receiver. It can be attributed to long monochromatic very low frequency Rayleigh waves. As shown in Fig. S3: elliptical particle motion is observed in the vertical plane and can provide information in the back-azimuth, direction of propagation of local station should converge in the source area.
For the largest very low frequency event of 11 November 2018, a first signal appeared at 9:27:27 UTC, followed by a moderate earthquake (M L 3.1 at 9:28:00), and then the VLP signal that lasts more than twenty minutes emerged (Fig. 7), as shown by the E-W raw record of the local MCHI seismometer from G30 (Fig. S2). Others embedded earthquakes appeared in the VLP signal, two of them could be localized within cluster 3, the one that marked the beginning of the very low frequency signal (M L = 3.1 at 9:28:00, 45.3990 • E, 12.7317 • S, depth 21 km, 12 phases) and the smaller one that followed it (M L = 3 at 9:29:33, 45.3683 • E, 12.9208 • S, depth 12 km, 10 phases). It should be mentioned that depth for such events should be considered with caution. A third, smaller earthquake occurs at 9:30:50, followed by two very small ones that we could not localize, but that are well visible in the raw trace shown in Fig. S2(a). The spectrogram allows us to separate the long VLP signal from standard volcano-tectonic earthquakes. The spectrum of the signal (Fig. S2b) shows a peak at 0.062 Hz, thus a period of ∼16 s corresponds to the VLP event integrated in the signal (green dot). Fig. 7 shows the three displacement components for the 11 November 2018 event at seven local and regional stations (Fig. 2), including the YTMZ accelerometer and velocimeters corrected by instrumental responses. After integration, the signal was filtered with a band of 0.05 to 1 Hz. Three distinct phases can be distinguished in the records. At the event onset, the first low frequency phase starts almost at the same time as the small high frequency signal that initiates the sequence (black dots, especially that of MCHI, on Fig. 7). For the local stations, this initial phase lasts around one minute. The second phase is the one with the oscillations and highest amplitude monochromatic waves (green dots). In the decay part of the harmonic oscillations, there are two embedded decay envelopes starting at the green and purple dots and lasting 20-30 min. The duration of the monochromatic very low frequency wave train seems not to depend on distance from the source: on the eastern and vertical component of stations on Fig. 7 very low frequency oscillations from local stations (YTMZ and MCHI) has the same duration than the ones from regional stations (CAB in Grande Comore, or SBV northern Madagascar). This observation suggests that the long duration very low frequency signal is linked to a source process rather than to a propagation effect. 11 November 2018 signal could be related to an oscillating source, in agreement with the model of resonance of a deep magma reservoir proposed by Cesca et al. (2020) Long period seismicity is considered to be related to magmatic and hydrothermal processes in a volcanic context. Their signal can be linked to oscillations within an excited fluid filled oscillator (e.g. Chouet 2003). Physical properties of the resonator can be deduced from the signal characteristics. We modelled the decay envelopes of the very low frequency event of 11 November 2018 (Fig. S4) with a model of a damped oscillator following the equation exp (-t/τ ) with τ = 360 s. From the physics perspective, this is compatible with the excitation followed by the free damping of a viscous fluid oscillating within a reservoir.
Figs 7 and S3 show that, especially at local stations YTMZ and MCHI, the VLP tremor affects mainly the vertical and E-W components, and much less the N-S component, as indicated by particle motion plotted from the decaying part of the signal in order to avoid noise induced by embedded events. The particle motion calculated for YTMZ and MCHI shows how the oscillation is polarized during the monochromatic decay phase. Horizontal polarizations of the closest stations are used in order to estimate the position of the source: in Fig. S3, the green lines from YTMZ and MCHI on the map are converging towards the aseismic gap between clusters 2 and 3, with the azimuth ∼N85 • from MCHI and ∼N100 • from YTMZ. An origin of this long lasting very low frequency oscillation in the vicinity of cluster 3 would be consistent with the location of the small located earthquakes present at the onset of the signal. Those earthquakes could have excited a fluid filled resonator, yet their locations need to be considered with caution.

GNSS data
Six GNSS stations (Table S6)   and data available from the surrounding islands of Anjouan, Mohéli, Grande Comore, Glorieuse. Fig. 8 shows the time-series at the station MAYG. Only this station has long enough time-series to estimate and remove a seasonal term (GNSS data processing is described in the Supporting Information).
The time-series (Fig. S5) show a major component of common mode. Four phases are visible, we call them: A, B, C and D. Phase A spans from the origin of the GNSS observations to the main earthquake of May 15 it shows the background deformation. No precursory deformation is observed at any of the stations. Then, from 15 May 2018, three distinct deformation anomaly phases are observed link to the seismo-volcanic crisis. During phase B, from 15 May to 30 May 2018, there is a small yet well visible deformation on the east component (Fig. S5). During phase C, from May 30 to July 3 there is little horizontal deformation at the stations and a slight global subsidence. During phase D, starting on July 3, all stations show a large and steady drift mostly east and downwards. Table S7 contains the velocity anomalies during phase D, obtained by linear fit. The time-series of the N-S component of the two southern stations, PORO and BDRL, show significant northwards velocity during phase D unlike the other stations. As the E-W velocities at all stations are very similar, we averaged them and, at the first order, we consider that the average vector represents the motion of the entire island of Mayotte at the barycentre of the six points located at 45.164 • E and 12.813 • S. Fig. 9 shows how the pattern of the cumulated seismic moment correlates with the average eastern velocity of Mayotte. Fig. S6 shows the evolution of the east versus north component, and east versus vertical. This figure is dominated by the large drift that occurs during phase D and shows clearly that this drift is linear at the first order during this first investigated period (early July-mid November 2018). Vertical versus east component plots for phase B, C and D are aligned, showing that the subsidence is not restricted to phase D, but occurs also, with lower intensity, during phases B and C. For north versus east, the situation is different and phase B, C and D occupy spaces not aligned together (Fig. S6).
Using the time-series of Fig. S5, we estimate the best-fitting velocity anomalies, thus removing the long term International Terrestrial Reference Frame (ITRF14) velocity and the seasonal component of MAYG-both assumed to be valid at all stations. The velocities anomalies are summarized in Table S7.
We finally analysed the data until mid-January 2020, and we split phase D in six periods of three months (D1 to D6) in order to assess the spatio-temporal evolution of the deformation during the crisis. Fig. S7 shows the time-series at four of the GNSS stations for the periods D1-D4 (July 2018-July 2019) and Fig. 8 MAYG for the six periods. The deformation evolves with time and at a lowering rate during phases D3, D4, D5 and D6 (January 2019-January 2020). Table S8 reports this observation quantitatively.

Interferometric analysis
In order to determine a dense pattern of deformation for the entire island of Mayotte we performed interferometric data analysis. This analysis involved the processing of Copernicus Sentinel-1 mission Synthetic Aperture Radar (SAR) data from both ascending and descending acquisition geometries for the periods of December 2017 to April 2019 (29 scenes) and June 2017 to April 2019 (56 scenes), respectively (the processing is described in the Supporting Information).
InSAR measurements well-outlined the ground displacement pattern of Mayotte Island, showing a ramp subsiding east (Figs 10 and S8), as a part of the overall concentric ground deformation common to volcano-tectonic unrests. The observed maximum rate change from one side of the island to the other (approximately along an E-W direction) reached −115 and −63 mm yr −1 for the ascending and descending tracks, respectively, which are perfectly consistent with the relative motion measured by the GNSS network. Differences of line of sight (LoS) measurements from opposite observation points indicate the presence of significant horizontal motion, as depicted by GNSS vectors. In that case, subsidence combined with motion away from the satellite for the ascending orbit results in higher displacement rates  compared to the descending one. The compatibility between InSAR and GNSS and agreement to the proposed phases of the crisis can also be seen by examining the InSAR LoS displacement time-series (Fig. 10).
It should be noted that as the entire island is deforming, In-SAR can provide only information regarding the relative motions within Mayotte, while the knowledge of the absolute motions is coming from the GNSS. Using the LoS changes estimated from 34 A. Lemoine et al. the GNSS time-series at nearby pixels allow us to overcome this limitation and thus derive absolute LoS time-series at all InSAR pixels.
Of high interest, and not achievable with GNSS, is the observation made by InSAR of the homogeneity of the displacement field within the whole island. This uniformity is an expected consequence of the distance of the source that we calculate in the next section. There is no evidence of landslide or any other land motion that could have been triggered by the repeated moderate seismic events during the crisis and could constitute a potential threat for the inhabitants.

Evidence of three phases of deformation
Based on the characteristics of the seismicity and deformation, we split the crisis in three unrest phases (B, C, D) preceded by the quiescent pre-crisis phase A. The three phases, issued from deformation anomalous phases, are interpreted as follows: phase B, from 15 May to 30 May 2018: seismic crisis and fracturing between a deep reservoir and the surface; phase C, from 30 May to 3 July 2018: magma ascent from the reservoir to the surface and SE migration; phase D, starting on 3 July 2018 and still ongoing in early 2020: eruption.
In the following, we first model phase D and then move backward in time to C and B. Indeed we need first to decipher the dominant phase D accurately as we will use its implications for the correct interpretation of the initial phases C and B.

Model of the phase D (eruption)
For modelling the ground deformation observed during phase D, we use the classical Mogi (1958) model that assumes a deflating point source buried in an elastic half space. Our data well constrain the source azimuth, pitch, and deflation, but not its distance to the stations. This is because all GNSS stations are in the same azimuth. Consequently the GNSS displacements are fit almost equally well with reservoirs in the range of longitude 45.55 • E-45.75 • E. As the ratio vertical/horizontal displacement indicates exactly the direction of the source centre, the modelled depth varies from 25 to 40 km when browsing the range of possible longitudes. On the other hand, if we do not use the vertical information but only the horizontal Downloaded from https://academic.oup.com/gji/article-abstract/223/1/22/5850758 by Ifremer, Bibliothèque La Pérouse user on 30 July 2020 The Mayotte 2018-2019 seismo-volcanic crisis 35 one, the preferred location of the source is at the average crossing place of the horizontal vectors at a longitude at around 45.50 • E thus slightly outside of the range permitted by the 3-D modelling. If now we blend in our inversion 50 per cent of cost of fit of the Mogi model with 50 per cent of cost based on the crossing place of the horizontal vectors (Fig. 11), using as metrics the minimization of the cross product of the displacement vectors and the vectors source-station, we find a best-fitting source at 45.590 • E, 12.777 • S and 28 km depth (depth of the deflation source is discussed in the Supporting Information). The blending of the two criteria in the inversion is interesting because the unconstrained Mogi modelling tends to pull the model further east (to give a better fit at one single observation, the E-W motion of MAYG, which means that these data have to be highly trusted), while the crossing of the horizontal vectors pulls the solution the other way. The combination of both, with similar weight acts as a joint inversion and gives a stable solution with mean residuals of 25, 5, and 20 mm yr −1 in the east, north, and vertical components, respectively (Table S7). With that geometry we infer an almost steady emission rate of 82 m 3 s −1 during the first months of the phase D.

Model of the phase C (magma ascent)
During the phase C, we interpret the deformation (presented in Table  S9, left) as due to the ascent of magma in a conduit connecting a deep reservoir to the subsurface, combined with the deflation of this reservoir of a similar volume. This scenario is supported by the geodetic observations as, during that phase, the E-W motion is nearly zero while there is a slight vertical motion (Fig. S5) of 7 ± 4 mm. Such ground motion can be explained with a magma body migrating upward without emission at the surface, thus without a global volume change. As the above-mentioned subsidence rate is 55 per cent of that of phase D and mostly attributable to the deflation of the reservoir (because a shallow source cannot produce vertical changes at GNSS stations at that distance), we assume that the horizontal rate associated with the deflation is also 55 per cent of that of phase D. Consequently, this allows us to evaluate what we call the C reservoir contribution in Table S9. The remaining part is what we call C conduit. We interpret it as the elastic response measured at the GNSS stations to the opening of the feeding magma conduit (Fig. 11) from the top of the chamber to the subsurface and the further opening of shallow dykes and sills. As the emission rate is 82 m 3 s −1 in phase D, rescaled by 55 per cent, it is 45 m 3 s −1 during phase C, which means that the volume injected in the channel, shallow dykes, and sills is 0.16 10 9 m 3 during phases B and C.
We model the C conduit deformation by using the Okada (1992) formalism and the inversion code developed by Briole et al. (1986) and Briole (2017). As the GNSS stations are far from the source, we cannot resolve the shape of the conduit. We assume that it is an elongated dyke buried at a depth of 22 km, at the top of a deflating chamber that we assume is 12 km in diameter (consistent with the extend of cluster 1a and with dimension of numerous well-known calderas, such as the ones reported by Geyer & Martí (2008), e.g. Campi Flegrei in Italy) with its tip 2 km beneath the seafloor. As a priori azimuth for the conduit, we use N304 • (Section 3.3), which is perpendicular to the tension axis. The best-fitting azimuth found by our inversion is N318 • , though loosely constrained because of the distance of the GNSS stations. The best-fitting solution for the width of the conduit is 4 ± 2 km, and its opening is 1.62 ± 0.7 m. Thus, the best-fitting volume of the modelled conduit is 0.13 × 10 9 m 3 , 20 per cent less than the a priori value. The root mean square (rms) scatter of the residuals is 3 mm for the horizontal components and 6 mm for the vertical. In the uppermost part of the crust, the conduit could extend laterally in a shallow dyke ∼20 km long, along the same axis ∼N318 • . This dyke could accommodate the missing 0.03 × 10 9 m 3 . The less active seismic cluster 2 might correspond to earthquakes induced by redistribution of the stress in the rock surrounding the magma reservoir deflation, as it is located close to its centre (e.g. Gargani et al. 2006).
There is no means with the geodetic data to resolve precisely the characteristics of the shallow channels (dyke and or sill), because their volume is small for the distance and they do not affect the vertical component. From the vertical GNSS data, we can only assess that there cannot be a large horizontal injection at depth (e.g. 10 km long at 10 km depth) because this would give a geodetic signature that is not observed.

Model of the phase B (fracturing of the crust)
For phase B, during which a small yet clear deformation is observed from 15 to 30 May 2018, we used an Okada (1992) model with the following parameters: same location and geometry as phase C, a rectangular fault elongated vertically of 20 km length and 4 km width (Fig. 11, Table S10). The best-fitting slip is 0.85 m of pure strike-slip. Assuming a rigidity of 3.3 10 10 Pa, the corresponding moment is M = 224 10 16 Nm, a value twice as large as the total seismic moment cumulated on May 30, but compatible with the cumulated moment a few days later on June 2. When relaxing the azimuth from its initial value of N318 • (inherited from phase C modelling), the best fit is obtained with N341 • , which corresponds to an azimuth with the tension axis oriented at N26 • , slightly smaller than those found before, but consistent with the G-CMT focal mechanisms. With both geodesy and seismology, we have no means to discriminate between the conjugate fault planes, and the upward fracturing during phase B might well involve ruptures that blend the two families of planes.

Temporal evolution of the deformation during the eruption (phase D)
In the section 6.2, we model the ground deformation by assuming, during phase D, a single source of deformation and assuming a constant rate of deformation during period D from 3 July until 15 November 2018. Below, we use the available data until early January 2020 (Fig. 8 for MAYG), and we split in six periods of three months (D1-D6), the evolution of the deformation during the eighteen months of the deflation signal. The modelling, made by assuming the same location (horizontal and depth) of the deflating reservoir, shows a clear evolution of the effusive rate during the eruption. The rate increases during the first phase of the eruption to reach a maximum around the end of 2018. This is the period during which the largest very low frequency tremor were observed on 11 November 2018, and during which clusters 2 and 3 reach steady activity (Figs 3 and 4). Then, the deformation gradually decreases during the first semester of 2019 (blue line on Fig. 12). Table S11 gives the effusive rate values for the six periods. According to this model, the cumulated volume extracted from the deep reservoir has reached ∼2.65 km 3 in January 2020, 1 yr and a half after the start of the eruption. This is less than the value reported from the estimated volume of the new submarine volcanic cone discovered at the sea floor (∼5 km 3 , REVOSIMA 2019; Feuillet et al. in rev.)

36
A. Lemoine et al. but consistent if one considers a 40 per cent increase of volume between magma at depth and the products at the surface.

Fracturing of the crust and magma injection
The intense seismic sequence begins offshore Mayotte with a first cluster, active from 10 May 2018 to early July 2018. During May 2018 (cluster 1a), intense seismicity (11 M L ≥ 5.0 events) is concentrated in a limited area with a diameter around 12 km. The first week of June represents the most intense week for the seismic sequence (equivalent to M w = 6.2), and the beginning of seismicity migration towards the surface and SE until the end of June 2018 (forming cluster 1b).
The GNSS data provide information on the timing of the creation of the feeding channel that opened the way for the magma to migrate from the reservoir to the surface.
This occurs during phases of the deformation corresponding to fracturing linked to strike-slip faulting, and then to magma ascent, between 15 May and 3 July 2018. Phase B, associated with moderate deformation from 15 May to 30 May 2018, is interpreted as the fracturing phase, and linked to the intense cluster 1a. Moreover, Downloaded from https://academic.oup.com/gji/article-abstract/223/1/22/5850758 by Ifremer, Bibliothèque La Pérouse user on 30 July 2020 the change in the deformation signature suggests that the filling of the conduit started at the end of phase B, during the peak seismicity (and rock fracturing) that occurs between 30 May and 7 June, 2018. During the 2014 Bárðarbunga dyke intrusion, most of the released seismic moment occurred during the rapid propagation phases (Ágústsdóttir et al. 2016). At the same time, cluster 1b shows intense seismicity migrating upward and SE (Figs 5 and 6) that attests to upward fracturing, while the conduit started being filled by magma until early July 2018. The fracturing of phase B can be explained geodetically without the need of filling the ruptured faults, that is it can be modelled thanks to coseismic displacements associated with strike-slip faulting.
We have no strong constraint on the geometry of the feeding conduit. The best-fitting model is an elongated vertical conduit starting arbitrarily at the depth of 22 km on top of the reservoir and ending at a depth of 2 km, yet the latter value is loosely estimated by the inversion. What is well constrained by the inversion, and almost independent on the shape of the source, is the volume of the conduit, 0.13 ± 0.01 10 9 m 3 , and therefore we can infer its section, 6500 m 2 . Assuming a pipe shape, its diameter would be 90 m. It is difficult to figure out how such a cylinder could appear in a few days within 10-20 km of oceanic crust not having been perforated by magma in recent geological time. Therefore, it is much more likely that the feeding conduit is a dyke. The best fit for this dyke is a horizontal width of 4 km, an opening of 1.62 m, and an azimuth of N318 • . The azimuth is loosely constrained because the conduit is narrow, but consistent with the overall SW-NE extension inferred from the focal mechanisms and the orientation of the extension from GNSS regional observations (Fig. 1).
The steadiness of the deformation in June 2018 implies a steady rate (within our observational uncertainties) of the amount of magma released from the reservoir during that month. This suggests that the cracks moved towards the surface and started being filled according to a relatively smooth process.
The G-CMT focal mechanisms are associated with the largest events of May and June 2018. Those events belonging to cluster 1, during which there is an apparent ∼20 km SE epicentral migration (Fig. 5). This migration was confirmed by location of the largest events by G-CMT (Fig. 6). However, there is evidence of upward migration: the event depths determined by G-CMT show a clear upward migration from early to late June 2018 (Fig. 6) that is coherent with depth evolution with time from our location, even if during this period, seismic network sparsity do not permit accurate depth determination for the whole earthquakes. During the 2011-2012 El Hierro unrest, a migration extending over 20 km and ending in an offshore eruption was precisely located (López et al. 2012;Martí et al. 2013). The Bárðarbunga-Holuhraun dyke intrusion (Iceland) in 2014, is associated to magma propagation along 48 km during 2 weeks (Ágústsdóttir et al. 2016, 2019). White and McCausland (2016) compiled data from seismic swarms linked to eruptions and intrusions. They noted that, generally, seismicity began at a distance that can reach tens of km from the eruptive location.
Focal mechanisms (from G-CMT) show the dominance of strike slip events during the fracturing period. In Fig. 6, the strike-slip focal mechanisms represent earthquakes belonging to the period from 14 May to 27 June 2018, that is corresponding to deformation phases B (fracturing) and C (magma ascent). Dyke related seismicity also sometimes produces the so called 'dogbone seismicity distribution', where normal faulting events around the dyke (normal faults parallel to dyke) coexist with strike slip events around it (e.g. Toda et al. 2002). In the case of the Mayotte seismo-volcanic phenomenon, only the largest events are associated with a focal mechanism, leading to a lack of accuracy in the description of the seismicity pattern.
We interpret cluster 1 as corresponding to fracturing of the crust followed by triggering by magma ascent. The seismicity for the first period of cluster 1 (cluster 1a during fracturing phase B) is roughly above the location of the reservoir (Fig. 11), the intensity of the seismicity during this period should be linked to the resistance of the crust probably not associated to a significant fracture density (e.g. Rivalta & Dahm 2004). Then, the upward and SE migration is observed from early June 2018 (cluster 1b, earthquakes triggered by magma ascent during phase C). Magma intrusions phenomena are associated with seismic swarms, as observed for example in 2000 in Izu islands (Japan) were an energetic seismic swarm was accompanied by dyke intrusion and eruptions of Myakejima (Toda et al. 2002;Rivalta & Dahm 2004). As reported during the Bárðarbunga-Holuhraun (

Eruption onset and the intense ground deformation
The seismic and geodetic data suggest that the proper eruption starts associated with a second seismic cluster, much less active than the first one, which appears at the end of June 2018. This cluster is centred NW in comparison with preceding events. It appears during a short intense period of seismicity (end of June 2018), marking the last stage of migration. Then, in July and August 2018, there is a period of relative seismic quiescence. Some days after the beginning of this period, a third seismic cluster appears, located to the west of the second one. Simultaneously, early July 2018 is characterized by the onset of an intense deformation recorded by GNSS, with large eastward and downward drifts of the entire island of Mayotte.
Those two clusters were still active in early 2020, with the seismicity of the third one remaining more intense. This third cluster is closest to Mayotte. Clusters 2 and 3 present b-values significantly higher than the one of cluster 1 (Fig. S1), which is consistent with values observed at some volcanoes (e.g. Wyss et al. 1997).
Because of the poor azimuth coverage of the GNSS stations (Fig. 11), the longitude of the deflation source is not easy to assess. However, as our two criteria (best-fitting Mogi source and horizontal crossing of the vectors) are leading in opposite directions, we believe that our localization of the deflating source is resolved within a few kilometres in a robust way. Subsequently, its depth (28 km) is constrained accurately by the pitch of the deformation and the trade-off between distance and depth in the model (Fig.  S6). Its azimuth is well constrained by the azimuths of the GNSS vectors. Moreover, this reservoir is localized in an area were level of seismicity is very low whereas it is surrounded by earthquakes (Figs 11 and 13). A similarly deep magmatic reservoir can be found in El Hierro (∼25 km, Martí et al. 2013;Klügel et al. 2015) or in the Klyuchevskoy volcano group in Kamchatka (z ∼ 30 km, Levin et al. 2014;Shapiro et al. 2017), linked to shallower sources.

38
A. Lemoine et al. Figure 13. Schematic summary. Yellow full filled circle represents the main deep reservoir inferred by the GNSS data. Best located M L ≥ 3.5 events of clusters 1, 2 and 3 are plotted in blue, yellow and red colour scales. The eruption takes place on the yellow triangle, as observed during marine surveys (Deplus et al. 2019;Feuillet et al. 2019, in review). It is connected to the vertical conduit by a ∼20-km-long dyke. The model includes a main deep reservoir and a secondary source beneath the south end of Mayotte island. A moderate flux of material (∼3 m 3 s −1 ) is inferred from that source to the main reservoir. Large dashed line: proposition for Moho position (e.g. Jacques et al. 2019from Dofal et al. 2018and Coffin et al. 1986).
Our modelling assumes that the source is a point, which is of course not the case, but we cannot resolve the geometry of the source with GNSS network as, at such distance, a source with a radius of several kilometres is indistinguishable from a point source. We assumed that the source diameter is 12 km, because this is a typical value used in literature for deep reservoirs, based on the size of many calderas that can be seen on Earth [e.g. the Campi Flegrei caldera in Italy as described in the database from Geyer & Martí (2008)]. It also ranges with the dimension of the area where earthquakes of cluster 1a were concentrated (interpreted as the fracturing phase preceding any magma ascent, Fig. 11).

Speed of magma ascent during the eruption
With 6500 m 2 for the best-fitting section of the feeding channel and 82 m 3 s −1 for the effusive rate, if we assume a constant emission rate between 3 July and 15 November 2018, the speed of magma ascent is 0.013 m s −1 . The speed is half this rate between 30 May and 30 June 2018, when migration of seismicity revealed the magma migration (during June 2018). Such a low speed of magma ascent has been observed elsewhere and is discussed by Gonnermann & Manga (2013) through the analysis of the crystals in the erupted magmas. This evaluation of the speed could be further compared and tested with the one derived from the analysis of the crystal growth in the fresh rocks sampled on the new volcano (e.g. Bachèlery et al. 2019).

Rates and emitted volume
The rate of magma released is less in June, approximately half, than the one estimated after 3 July 2018. In other words, the flux in the feeding channel doubled in July once the path to the surface opened, and the eruption started. The timing of the start of the eruption, estimated on 3 July 2018, comes from the beginning The Mayotte 2018-2019 seismo-volcanic crisis 39 of the large ground deformations, thus, when there is no longer a balance between the volumes collected at large depth in the reservoir and those stored at shallow depth. Early July corresponds also to the beginning of a seismicity quiescence period, and the onset of a new pattern of seismicity: only few small shallow events around the dyke and the onset of a new and steady seismicity pattern within clusters 2 and 3. In the same way, during the Bárðarbunga-Holuhraun rifting event (Iceland) in 2014, a rapid decrease of the seismicity followed the beginning of the main eruption and a lack of shallow seismicity is observed (Ágústsdóttir et al. 2016, 2019), like in other contexts (e.g. during the 1983 Kilauea dyke intrusion, Rubin et al. 1998). Persistent seismicity is not associated to the propagation of a dyke (as during phase of magma ascent in Mayotte), but to processes of continual fracturing that maintains magma flow (e.g. Ágústsdóttir et al. 2016). Moreover, Ruch et al. (2016) highlighted importance of pre-existing fractures in rifting events as they can control the magma propagation and new distribution of stresses.
With an effusion rate of 45 m 3 s −1 over 42 d (magma ascent phase C) and then, after the estimated beginning of the eruption early July 2018, six phases over 3 months with rates 72.9, 92.4, 77.2, 40.7, 27.0 and 12.2 m 3 s −1 , the total volume extracted from the chamber is 2.65 km 3 , according to our modelling. White & McCausland (2016) proposed an empirical tool for estimating the intruded magma volume from the analysis of a large range of seismicity patterns that preceded past eruptions or intrusions. Considering the cumulative seismic moment between 10 May 2018 and 15 May 2019, the estimated magma volume should be of the order of 1.55 km 3 , lower than what we predict at this date from GNSS data (∼2.21 km 3 ). Moreover, recent marine surveys (R/V Marion Dufresne) performed bathymetric and water column observations that confirm recent submarine volcanic activity, 50 km eastward from Mayotte and differential bathymetry allowed to estimate an important volume emitted at the seafloor of ∼5 km 3 (REVOSIMA 2019; Feuillet et al. in rev.). The volume of released magma from the reservoir that we estimated without a priori assumption is lower than the volume deduced from direct seafloor observations of the eruption. This discrepancy can have different origins, including potential volume change of material between its extraction at 28 km depth and the seafloor (Bachèlery et al. 2019).
Such volume ranks the Mayotte eruption among the major events of the last 1000 yr. It can be compared only to other very large eruptions such as the Lanzarote one between 1730 and 1736, that emitted 3-5 km 3 , associated with an effusion rate five times lower than the one estimated for Mayotte with comparable volumes (Carracedo et al. 1992). The Laki eruption spanned from June 1783 to February 1784 and emitted 14.7 km 3 of magma, with an effusive rate nearly four times higher than the one from Mayotte (Thordarson & Self 1993).

11 November 2018 very low frequency tremor
Only two events reported in 2011 and 2013, close to Rocard submarine volcano, Society hot-spot, Polynesia (Talandier et al. 2016) could be compared to Mayotte unusual 11 November 2018 signal: both 2011 and 2013 intense Rayleigh wave trains were monochromatic (the period observed was 17.0 s) and long-lasting (50-60 min in 2011 and 30-40 min in 2013). Talandier et al. (2016) interpreted these events as due to the strong hydrostatic pressure undergone by the volcanic edifice at 3200-4000 m depth below the sea surface. Under high hydrostatic pressure, magma movement can generate the resonance of a shallow opened conduit considered as a fluid filled reservoir after punctual excitation.
If those unusual events belongs to VLP seismicity (regarding frequency range), their waveforms look like of some LP events associated with the decaying oscillation in a resonator. Among them, some earthquakes are described by decaying harmonic oscillations after a short onset in the period range 0.2-2 s, characterized by one or more frequency peaks related to a source effect (e.g. Chouet 1996;Chouet & Matoza 2013). This pattern is interpreted as a timelocalized pressure excitation mechanism followed by resonance of a fluid filled cavity in response to the excitation. Damped oscillations can be described by its dominant frequency and quality factor that depend on the source properties: fluid types and resonator geometry (e.g. Kumagai & Chouet 2000). Characteristic frequencies seem to be different from one volcano to another. Signatures of signals recorded at Kusatsu-Shirane or Galeras volcanoes present a typical signature of long harmonic codas following a brief onset (Kumagai and Chouet 1999) that exhibit some similarities with the Mayotte 11 November 2018 event, but with a lower range for the characteristic period and duration. Among the family of VLP events (presenting longer periods), signals observed at Hachijo volcano in 2002 present a decaying harmonic oscillation lasting 300 s with period ∼10 s, a pattern that is similar to features of Long Period (LP) events with shorter period ∼1 s. These VLP signals at Hachijo volcano were interpreted as the resonance of a basaltic magma in a dyke (Kumagai et al. 2003). Kumagai (2006) proposed that the variation of frequency and quality factor of the harmonic oscillation are due to the gradual expansion of a crack containing basalt mixed with gas. Kawakatsu et al. (2000) reported that Long Period Tremors (LPT) with a dominant period around 15 s were observed at Aso volcano (Japan), with a short duration (less than one minute), both during periods of volcanic quiescence and unrest. These LPT are interpreted as a response of resonating fluid filled cracks to pressurization within the hydrothermal system at Aso volcano (e.g. Hendriyana & Tsuji 2019).
The 11 November 2018 signal emitted offshore Mayotte (and similar reported VLP events) is associated with a frequency peak included in the range of VLP seismology, but its waveforms are typical of some long period events attributed to a decaying oscillation in a resonator. This event can be considered as an end-member of sources linked to resonator oscillation after excitation, similar to the two events reported in 2011 and 2013 in Polynesia (Talandier et al. 2016).
During the Mayotte seismo-volcanic crisis, the source process of long-lasting VLP events is non-destructive as several hundreds of them were reported since the beginning of the unrest and even few months before (Poli et al. 2019;Satriano et al. 2019;Cesca et al. 2020). The characteristic period of such events increased until late 2018 before decreasing (e.g. Satriano et al. 2019;Cesca et al. 2020). This period of highest period of resonance corresponds to the highest modelled effusion rate (Fig. 12). In parallel, volcano-tectonic seismicity pattern follows a steady pattern for several months from this period: activity is mainly concentrated in cluster 3 and cluster 2 is still active (Fig. 3). Geometry of the volume of the resonator and characteristics of fluid filling should influence VLP signals. The constraints on the geometry of the deep feeding channel, and a possible shallow dyke, can be used as a priori geometric model for the interpretation and modelling of the reported very low frequency tremors. Those could be generated by an oscillation, potentially in a long horizontal shallow dyke or in the magma reservoir, but not in the vertical conduit as the particle motion does not support the hypothesis of a vertical vibrating source. Cesca et al. (2020)

40
A. Lemoine et al. linked the long-lasting VLP events to a 15 km long subhorizontal crack forming a deep reservoir. However, constraining the source of 11 November 2018 event does not belong to the objectives of this work, geodetic network do not allow us to constrain the geometry and the complexity of the plumbing system: future analysis of such unusual events and onshore and offshore acquisitions could help us understand it.
Other tremors are recorded before and after this major tremor, weaker but with very similar characteristics (Poli et al. 2019;Satriano et al. 2019;Cesca et al. 2020), suggesting a common source and a non-destructive process.

Spatial distribution of the source of deformation
The GNSS time-series (Fig. 8) show a change of deformation rate and an increase of the pitch (ratio N-S/E-W velocities) of the deformation with time. To take into account this observation, we refine our model with a dual source. The small secondary source can be regarded as a perturbation of the main one. The localization and evolution of this secondary source is expected to provide clues on the dynamic of the process at depth.
For each of the subphases D1-D6 we now invert for the main source plus the small secondary source. Whatever are the initial conditions set for the secondary source, the inversions bring this source towards a position which is located beneath the south of the Mayotte island around coordinates 45.18 • E and 12.94 • S (Fig.  13). Only the depth of the source cannot be constrained by our model, so we decided to fix it at the same depth as the main source. Table S12 gives the values estimated for this dual source model. With respect to the single source model, the overall volume of magma released is unchanged with most of the magma (100-94 per cent) still extracted from the main reservoir and the secondary reservoir providing a small amount of the total (0.1 per cent in phase D1, 4.1 per cent in phase D2, 4.5 per cent in phase D3, 5.8 per cent in phase D4, 6.6 per cent in phase D5 and 14.5 per cent in phase D6). At the same time, this small perturbation of the model improves the fit to the data by 23, 9, 23, 22 and 14 per cent, respectively, during phases D2-D6, which is remarkable. In terms of emitted volumes, the relative importance of the secondary source is increasing with time.
The modelled secondary source might correspond, at the first order, to a flux of material extracted from beneath Mayotte where the main asthenospheric plume that has built the island should be located. This flux might contribute dynamically to maintain the isostatic equilibrium of the whole system and in that case its time span will be controlled by the average viscosity of the deep rocks involved in the process. This relaxation might continue well after the end of the eruption that could occur by the mid of 2020 based on the trend shown by Fig. 12.
Other authors have studied the connexion between reservoirs at depth. For instance, Bato et al. (2018) proposed deep connections between the two volcanic systems of Grímsvötn and Bárðarbunga (Iceland), ∼25 km away explained by lateral flow hypothesis at depth or shared magma reservoir at depth ∼30 km. From the 2011 to 2012 submarine El Hierro eruption analysis, Klügel et al. (2015) suggested that lateral magma movements at depth are typical of mature oceanic intraplate volcanoes. They proposed that between 2011 and 2014, pre-, syn-and post-eruptive seismic swarms testify to the presence of magma paths in the uppermost mantle and lower crust.

Interpretation of the seismic clusters
Phase D of intense deflation, which has run for more than 1 yr at the time of this paper, is the phase with the lowest seismic rate. This indicates that the fracturing from the reservoir depth to the subsurface is terminated. The second cluster, which spans from the end of June 2018 until now, is located in the vicinity of the deflation centre, likely bellow it, as far as depth distribution is reliable (Fig. 11). During the El Hierro eruption (2011-2012), Martí et al. (2013) interpreted the syn-eruptive seismicity as a response to readjustments of the plumbing system following magma withdrawal. As the plumbing system in the case of Mayotte is not described, it is difficult to propose a reliable explanation for the sources of clusters 2 and 3, both of which are syn-eruptive. A new swarm or geometry in the seismicity pattern can mark changes in magma circulation (e.g. segmented dyke growth in 2014 in the Bárðarbunga volcanic system; Sigmundsson et al. 2015). In the present case, it seems that there is a seismicity gap between clusters 2 and 3, except maybe at larger depths (Figs 4, 5 and 11), however, they are both syn-eruptive and almost steady during months. Cluster 2, due to its vicinity with the centre of deflation, should be linked to the redistribution of stress around the reservoir due to magma withdrawal. Post-eruptive clusters could also reflect regions of magma circulation. Moreover, the third cluster could be linked to fluid movement maybe between two reservoirs (Fig. 13). Cesca et al. (2020) interpreted earthquakes observed offshore Mayotte the followed the onset of the eruption as due to the progressive failure of the roof of a magma reservoir (then inducing its resonance and VLP events).
The three seismic clusters show an overall migration of activity SE and then NW and W, indicating a coaxial complex and extended activity ranging up to ∼40 km depth. All focal mechanisms preceding the eruption (May-June 2018) made by G-CMT indicate an extension axis oriented ∼N50 • . The perpendicular azimuth N310 • is therefore the most favourable axis for a feeding dyke. The centre of deflation, is ∼20 km from the new volcanic edifice (Deplus et al. 2019;Feuillet et al. in rev.) in the azimuth N316 • (Fig. 13), and thus, close to the abovementioned N310 • and to the azimuth of the feeding conduit found by inversion of GNSS data which is N318 • . This azimuth is also consistent with some topographic features observed in the bathymetry east of Mayotte, for example the Jumelles seamounts (NE from Mayotte, Fig. 1), and with that of dykes exposed onshore north of Mayotte (Nehlig et al. 2013). It is also consistent with the overall axis of the Comoros volcanic Archipelago and is roughly perpendicular to the relative motion measured by GNSS (Fig. 1) in azimuth N50 • across the Comoros. So, the episode of magma ascent and dyke intrusion is consistent with regional kinematics and tectonics. Feuillet et al. (in rev.) suggested that the new volcanic edifice belongs to a volcanic ridge, extending between Mayotte up to the new volcanic edifice, which is part of a regional enéchelon system where dextral movement between EARS and Madagascar is accommodated. Famin et al. (2020) described Comores archipelago as the northern boundary of Lwandle Plate were volcanic cones are associated to Riedel shear and enéchelon extension fractures.
The 2018 Mayotte crisis shows some similarities with past telluric crises involving dyke intrusions. For example, Ahmed et al. (2016) suggested a dyking event related to intense seismicity involving 29 M > 5.0 earthquakes in the Gulf of Aden over a few months, where seismic migration is interpreted as magma ascent before dyke propagation. As in the case of Mayotte, cluster 1b (June 2018) may correspond to a phase of shallow injection (upward and lateral The Mayotte 2018-2019 seismo-volcanic crisis 41 migration). Regarding the reservoir depth (28 km) or distance of migration, the Mayotte seismo-volcanic crisis shares some similarities with that of 2011-2012 at El Hierro, Canary Island (deep reservoir depth: 20-25 km), where the activity was characterized by several phases, including phases of migration and submarine eruption (e.g. Martí et al. 2013).

C O N C L U S I O N
We proposed a scenario for the onset and the first 18 months of the volcanic eruption that started off-shore Mayotte in early July 2018 and was eventually confirmed by direct observation 10 months later (Feuillet et al. in rev.). Our scenario is mostly based on the available seismological, GNSS and InSAR data acquired at Mayotte and surroundings during and before the crisis.
The first emergence of the seismic crisis was on 10 May 2018 and the main shock on May 15.
Two phases preceded the proper eruption that initiated around July 3 according to the GNSS data.
First, from May 10 to May 31, the intense seismic activity, with most of the largest events of the crisis (29 among the 32 M L ≥ 5.0 events reported), is interpreted as a phase of fracturing of the crust but without the intrusion of magma in the fractures. This is strongly supported by the pattern of the deformations during those 3 weeks.
Then, during the entire month of June 2018, a second and distinct phase occurred, with the injection of the magma in the fissures. This injection phase started with a week of intense seismicity (equivalent to M w = 6.2) and the migration of the seismicity towards southeast and the surface. The GNSS data are also recording the magma ascent during this period at an average speed that we can estimate around 0.013 m s −1 . This speed might be confirmed by the analysis of the crystals present in the fresh basalts collected at the sea bottom on the slopes of the new volcano, or could be used to constrain the analysis of those crystals.
Based on the clear and abrupt change in the GNSS time-series, we guess that the proper start of the eruption is on 3 July 2010 ± 2 d, with a flux of magma that doubles within the feeding cracks with respect to June, and the beginning of building of the new volcano at the sea floor.
During the 2 months of eruption that follow, in July and August, the seismicity is low, while the slope of the GNSS changes is totally stable, thus the flux of magma perfectly steady ±5 per cent.
The modelling of the ground deformation is robust but with a 20 km range of possible models along the longitude axis, due to the peculiar geometry of the GNSS stations. Our preferred model is with a deflating reservoir located ∼30 km east of Mayotte at 28 km depth.
At the end of June, approximately at the same time the proper eruption starts, a new cluster of seismicity appears, mostly NW from latest events of the first cluster and close to its initial events. Its activity is moderate in comparison with the first cluster, but it remains active during the entire crisis.
A third seismic cluster starts also around the beginning of the deflation, west from the second one. This third cluster is more active than the second one and was still active at the end of the observation period. Those two clusters might correspond to the propagation of part of the magma at depth in dykes or sill, but it might also correspond to cracks triggered by the on-going crisis without the need of invoking the injection of magma in those cracks. There is no signature of those clusters in the GNSS data that are totally dominated by the intense deflation.
An intense long-lasting monochromatic very long period (∼16 sec.) seismic event occurred on 11 November 2018, lasting more than 20 min, corresponding presumably to the oscillation of fluid in the system, either magma in the chamber or in a long dyke. This tremor was roughly located between clusters 2 and 3.
The peak of the effusion rate, 93.8 m 3 s −1 is dated around the end of 2018. Then during the entire year 2019 the effusion rate decreased smoothly. Between October 2018 and April 2019, a sustainable increased seismicity appears, essentially regarding cluster 3 but also cluster 2. Then, at the end of April 2019, seismicity enters a new slowdown phase. The total magma released from depth, according to our model, and for the period July 2018-January 2020, is estimated between 2.5 and 3 km 3 . Whereas the plumbing system and local structure are poorly known, adding a smaller secondary source in our model of deflation improves very significantly the fit of the GNSS data. This second source that we localize, quite robustly, beneath the south of the Mayotte main island might represent the area where material is predominantly collected to partly refill the major depleted area in a process that might be a viscous relaxation flow.
The seismo-volcanic crisis of 2018-2019 offshore Mayotte has shown activity of a WNW-ESE aligned system. It marks the end of a period of quiescence several 1000 yr of the Mayotte island edifice and occurred in a poorly known and poorly instrumented area. However, it confirms the alignment of active seismicity and volcanism around Comoros archipelago, between the Somalian Plate and Lwandle block in an overall extensive regional context, consistent with trans-tensional regime in the area. The ongoing crisis in Mayotte might be important for teaching us more about the dynamics between the Somali Plate and the Lwandle block and the related volcanism. Improving this knowledge and improving the knowledge of the distribution, alignment, and ages of the offshore volcanic features, especially around the main islands, may lead to a better understanding of the behaviour, evolution, and related hazards of this peculiar area.

S U P P O RT I N G I N F O R M AT I O N S
Supplementary data are available at GJ I online. Figure S1. b-Values for the three clusters. Figure S2. The 11 November 2018 very low frequency tremor and embedded events. Figure S3. The 11 November 2018 very low frequency event: particle motions. Figure S4. Modelling of the tremor decay. Figure S5. Time-series for the coordinates of the six GNSS stations during phases A, B, C and D. Figure S6. Ratio between horizontal and vertical motion during phase D. Figure S7. Time-series since January 2018 at stations BDRL, GAMO, KAWE and MAYG. Figure S8. InSAR (Sentinel-1) vertical and E-W components rate maps. Figure S9. Parametric analysis of the depth and source type. Table S1. ITRF2014 horizontal velocity of GNSS stations around Mayotte. Table S2. The stations of the Mayotte monitoring seismic network. Table S3. The 32 events with M L ≥ 5.0. Table S4. Regional velocity model Table S5. Three intense Very Long Period tremors. Table S6. The GNSS stations in Mayotte. Table S7. Velocity anomalies at the Mayotte GNSS stations during deformation phase D. Table S8. Temporal evolution of the velocity anomalies during phase D. Table S9. Deformation during phase C. Table S10. Deformation during phase B. Table S11. Evolution of the effusion rate. Table S12. Evolution of the effusion rate estimated over periods of 3 months each by assuming a dual source. Table S13. Analysis of the fits and effusive volumes as a function of source depth and model Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.