# Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation

Simone Cesca, Jean Letort, Hoby N. T. Razafindrakoto, Sebastian Heimann, Eleonora Rivalta, Marius P. Isken, Mehdi Nikkhoo, Luigi Passarelli, Gesa M. Petersen, Fabrice Cotton, et al.

<figure><img src="image_1.png" /></figure>

► To cite this version:

Simone Cesca, Jean Letort, Hoby N. T. Razafindrakoto, Sebastian Heimann, Eleonora Rivalta, et al. Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation. Nature Geoscience, 2020, 13, pp.87-93. 10.1038/s41561-019-0505-5. insu-04836354

HAL Id: insu-04836354

https://insu.hal.science/insu-04836354v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation

Simone Cesca<fnref n="1" />*, Jean Letort², Hoby N. T. Razafindrakoto<fnref n="1" />, Sebastian Heimann<fnref n="1" />, Eleonora Rivalta<fnref n="1" />, Marius P. Isken³, Mehdi Nikkhoo<fnref n="1" />, Luigi Passarelli⁴, Gesa M. Petersen<fnref n="1" />,⁵, Fabrice Cotton<fnref n="1" />,⁵ and Torsten Dahm<fnref n="1" />,⁵

The dynamics of magma deep in the Earth's crust are difficult to capture by geophysical monitoring. Since May 2018, a seismically quiet area offshore of Mayotte in the western Indian Ocean has been affected by complex seismic activity, including long-duration, very-long-period signals detected globally. Global Navigation Satellite System stations on Mayotte have also recorded a large surface deflation offshore. Here we analyse regional and global seismic and deformation data to provide a one-year-long detailed picture of a deep, rare magmatic process. We identify about 7,000 volcano-tectonic earthquakes and 407 very-long-period seismic signals. Early earthquakes migrated upward in response to a magmatic dyke propagating from Moho depth to the surface, whereas later events marked the progressive failure of the roof of a magma reservoir, triggering its resonance. An analysis of the very-long-period seismicity and deformation suggests that at least 1.3 km³ of magma drained from a reservoir of 10 to 15 km diameter at 25 to 35 km depth. We demonstrate that such deep offshore magmatic activity can be captured without any on-site monitoring.

In May–June 2018, global earthquake monitoring agencies detected a series of large earthquakes in an otherwise relatively quiet area ~35 km east of the island of Mayotte. A peak $M_W$ 5.9¹ earthquake occurred on 15 May 2018, the largest ever recorded in the region. In November 2018, long-duration (~20 min) very-long-period seismic signals (VLPs) with periods of ~16 s were discovered in global seismic recordings<fnref n="1" />, triggering the curiosity of the scientific community. The origin of the VLPs was traced back in the vicinity of the swarm activity close to Mayotte. Most often, VLP signals² are recorded in volcanic areas³ suggesting a magmatic origin for the sequence.

Mayotte is one of the four principal volcanic islands in the Comoros archipelago and home to a population of ~256,000 (2017). The island belongs to a 250 km long NW–SE chain of basaltic volcanoes located between Africa and Madagascar (Fig. 1). The region has been affected by multiple tectonic processes⁴,⁵ including an episode of NE–SW trend rifting during the Permo–Triassic, associated with the fragmentation of Gondwana⁶,⁷ and the formation of the Somali and Mozambique oceanic basins, during which Madagascar drifted southwards⁸⁻¹². Proposed sources of volcanism include hotspot¹³, passive magma ascent through lithospheric discontinuities¹⁴ or rifting coeval to rifting in Southern East Africa¹⁵. Volcanism at Mayotte started at about 10–20 Ma (ref. ¹⁵,¹⁶) and subsequently migrated to produce the other islands. Mayotte last erupted 4,000 ± 500 years ago¹⁵. The nature of the crust and its thickness beneath the Comoro Islands is debated¹⁴,¹⁵. Only a few $M > 4$ earthquakes have been recorded in this area, including the 1993 $M_b$ 5.2 event, which caused ~1.7 M Euros of damage in Mayotte, and the 2011 $M_b$ 4.9 (ref. ¹⁶). Focal mechanisms and GNSS data support NE–SW transtension¹⁷,¹⁸.

By modelling seismological data at regional and teleseismic distances and ground displacement recordings at Mayotte (Supplementary Fig. 1) collected between May 2018 and April 2019,

we provide evidence of the drainage of 1.7 ± 0.4 km³ of magma from a ~30 km deep sub-Moho magma reservoir by a dyke that propagated to the seafloor. The dyke propagation caused a swarm of almost 7,000 volcano-tectonic earthquakes (VTs) and the downsag of the host rock overlying the reservoir, which in turn triggered its resonance, emitting 407 long-duration VLPs, and helped sustain a high outflow rate throughout. We argue that reservoir roof failures at calderas, which are well-studied, can be used as a downscaled analogue to evaluate future scenarios.

## Data analysis and modelling

We use seismic data at regional and teleseismic distances (ABKAR seismic array, Kazakhstan) and test different velocity models (Supplementary Fig. 2, Supplementary Tables 1–4) to perform full waveform moment tensor (MT) inversions and depth phase analysis, thereby retrieving focal mechanisms, centroid locations (Supplementary Figs. 3–5) and depths (Supplementary Figs. 6–8) for the most energetic VTs¹⁹. Additionally, we take advantage of seismic station YTMZ—deployed over the full study period on Mayotte—to produce an enhanced catalogue¹⁹ of relative locations of weaker VTs (Supplementary Fig. 9). We detected 6,990 VTs and located 1,904 of them (Fig. 1) by retrieving their origin direction and distance (Supplementary Figs. 10 and 11; see Supplementary Information). The VTs are classified into families of events (colour-coded in Figs. 1 and 2a–c) with similar waveforms and distance to station YTMZ (Supplementary Fig. 12). We design a detection tool to scan the broadband data at four regional seismic stations for monochromatic, low-frequency signals, thereby detecting 407 long-duration VLPs¹⁹ (Supplementary Fig. 13). We also develop an algorithm to invert for the centroid MT of the VLPs and resolve the damping constant and dominant frequency of a damped linear oscillator as the source time function, obtaining high-quality solutions for

<footnote n="1">¹GFZ German Research Centre for Geosciences, Potsdam, Germany. ²Observatoire Midi Pyrénées, IRAP, CNRS UMR 5277, Université Paul Sabatier, Toulouse, France. ³University of Kiel, Kiel, Germany. ⁴King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia. ⁵University of Potsdam, Potsdam, Germany. *e-mail: simone.cesca@gfz-potsdam.de</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** | Map view and cross sections of seismic and deformation sources. **a**, Weak VT locations (points) and strong VT MTs (focal sphere, lower emisphere projections of the double-couple components in map view and backprojection of full MTs from East or South in the cross sections, respectively) plotted and colour-coded according to similarity of waveform and distance to station YTMZ (yellow triangle), overlay of 22 VLPs MTs (black focal spheres, negative CLVD convention used), best-fit point deformation sources (black symbols, see Supplementary Information), vertical displacements (black lines) and time evolution of horizontal displacements at four GNSS stations (black scatter points, 1 July 2018 to 1 April 2019); uncertainties are within 2σ of vertical displacements (green boxes, horizontal uncertainties are negligible with respect to vertical ones), best fitting modelled displacements (red lines and arrows) and the inferred VLP source location (dashed red ellipses). The southeast and upward migration paths of seismicity and magma in Phases I–II (dashed black arrows in the cross-sections) reach the location of the discovered seamount<sup>20</sup> (orange circle and bars, assuming a 5 km diameter<sup>20</sup>). **b**, Map of the Comoro Islands, major regional tectonic structures and basin configuration<sup>5</sup>, showing fracture zone (thin lines) and earthquake locations, and MTs before the sequence (circle and focal spheres, source GlobalCMT, GEOFON and USGS catalogues); shaded-relief topography and bathymetry are from the ETOPO1 Global Relief Model from NOAA's National Centers for Environmental Information. The location of the study region is highlighted by a red star on the Earth globe.</figcaption></figure>

22 VLPs¹⁹ (Fig. 3). Due to the uncertainties on Moho depth and crustal structure, all our results are tested against different velocity models with varying Moho depth (see Supplementary Information).

## Chronology of the crisis

We have identified four phases of the volcano-seismic crisis. During Phase I (10 May to 7 June 2018) the most energetic VT burst of the sequence ($M_W > 5$ events between 15 and 20 May) occurred ~35 km East of Mayotte (Fig. 2a). We find a depth of 22 km for the largest $M_W$ 5.9 earthquake. Between 31 May and 7 June 2018, the seismicity migrated upward, as confirmed by the centroid depth estimations and array beam analysis (Fig. 2c, Supplementary Figs. 7 and 8). During Phase II (7 June to 18 June 2018) epicentres migrated southwards approaching the recently discovered seamount²⁰ (Fig. 1, Supplementary Fig. 5). Left-lateral strike-slip faulting dominates in Phases I–II (Fig. 1). A consistent expansion-related isotropic component, increasing gradually up to 20% of the total moment, is found during the upward migration and in Phase II (Supplementary Figs. 5 and 6). By mid June, a new type of activity emerged in the form of

long-duration VLPs. However, early VLPs had already occurred on 30 January and 2 June.

Few $M_W > 4.5$ VTs (Fig. 2a) occurred in Phase III (28 June to 17 September 2018), while VLPs became dominant. Phase IV (17 September 2018 to March 2019) started with an increasing VLP rate, this time accompanied by a new type of VTs (Fig. 2c,d). The dominant period of VLPs increased smoothly from ~15.2 s (June 2018) to a maximum of ~15.6 s (October 2018), before decreasing again to ~15.3 s (February 2019; Fig. 2e, Supplementary Fig. 13).

The VLP MTs are similar (Fig. 1), with a predominance of alternating positive and negative vertical Compensated Linear Vector Dipole (CLVD) (Fig. 3). The vertically axisymmetric MTs suggest a subhorizontal crack, dipping slightly to the west. Centroid locations are the same within errors and depths are at $37 \pm 11$ km. The largest VLP occurred on 11 November 2018, with an estimated surface wave magnitude of $M_s$ 5.1. The damping of the VLP source (quality factor, Q) is in the range $72 \pm 6$, which is higher than in other cases²¹,²² but not unique²³. VLPs onsets are often quasi-simultaneous to one or multiple weak VTs, as observed previously²³. Due to the emergent
<figure><img src="image_4.png" /><figcaption>**Fig. 2** | Timeline of the seismic sequence. **a**, VT magnitudes (colour-coding as in Fig. 1) and seismicity rate (black line). **b**, VT depth based on MT inversion (circles) and array analysis (crosses). **c**, Differential S-P time of VTs at station YTMZ as a proxy of the distance to Mayotte (see Supplementary Information for uncertainties). **d**, VLP magnitudes ($M_s$, purple circles) and rate (black line). **e**, VLP dominant periods (purple circles). **f**, Demeaned and detrended East, North and vertical GNSS displacements at station MAYG. **g**, Normalized cumulative number of dyke-related VTs (the blue line refers to blue, cyan and purple VTs in **a–c**), sagging-related VTs (the red line refers to red and green VTs in **a–c**), VLPs (indigo line) and length of GNSS displacement vectors at station MAYG (green line). Sequence phases are marked in all panels with cyan vertical bars.</figcaption></figure>

nature of VLP signals, it is difficult to judge whether VTs precede or follow VLPs, but some VLPs respond to VTs with abrupt signal amplitude changes, depending on whether the successive VTs are in or out of phase with the resonating VLP source (Fig. 4a). This observation suggests that VTs act as repeated forcing, contributing to the exceptionally long duration of some VLPs. The VLP spectra (Fig. 4b) reveal higher modes including non-integer ratios between them.

All VTs accompanying VLPs and most VTs in Phases III–IV (red and green dots, Figs. 1 and 2a,c) are located closer to Mayotte and have steep NE–SW striking thrust mechanisms (Fig. 1), inconsistent with local transtension. An isolated burst of VTs at the end of

August 2018 included earthquake repeaters and anti-repeaters, with highly correlated and anti-correlated waveforms, respectively (Fig. 4c). VTs and VLPs have persisted throughout Phase IV until the time of writing. The cumulative moment of VTs in Phases III–IV is $M_0 = 5.42 \cdot 10^{16}$ Nm, corresponding to $M_w = 5.1$, far less than for the dyke-related seismicity in Phases I–II ($M_0 = 2.03 \cdot 10^{18}$ Nm, $M_w = 6.2$).

Geodetic data from Mayotte Island reveal a steady, long-lasting subsidence and eastward displacement at four GNSS stations. The length of the displacement vector from July 2018 through April 2019 is 18 cm. Early ground deformation is weak, and grows clearly from July onward. We constrain the location, depth, aspect ratio and
a

<figure><img src="image_5.png" /></figure>

b

<figure><img src="image_6.png" /></figure>

c

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

**Fig. 3** | MT solution for the 31 July 2018 VLP. **a**, The MT plot overlays the focal spheres of the ensemble of bootstrap solutions (see Supplementary Information). The red line denotes the overall best solution. The negative CLVD convention is applied. **b**, The resolved, oscillating source time function, with a dominant period $T_R = 15.4$ s. **c**, Comparison of observed (red lines) and synthetic displacement traces (black lines), as well as their differences (red areas) for a selection of stations and components. Synthetic traces were computed assuming a continental (model P2017, see Supplementary Information) or oceanic crust, based on station locations (station name, spatial component, epicentral distance, azimuth, time window starting time, respect to the origin time, and time window duration are reported for each subplot).

volume change by inverting these GNSS data for the time interval 1 July 2018 to 31 March 2019. An isotropic source of pressure is unable to simultaneously fit the ratio of horizontal to vertical displacement and the vectors' orientation¹. We use generalized point-source deformation models²⁴,²⁵, first assuming a vertical ellipsoidal source²⁶ and investigating the effect of layering on the inferred source shape (Supplementary Figs. 14 and 15). The best fit is given by a large negative isotropic plus a small positive vertical CLVD source $12 \pm 4$ km East of Petit Terre at a depth of $32 \pm 3$ km, volume change is here $-1.7 \pm 0.4$ km³. We also test non-axisymmetric generalized point sources in a homogenous medium²⁴, obtaining a laterally contracting vertical dislocation as the best-fit deformation mechanism at a

depth of 23 km and expelled magma volume (potency) of $-1.3$ km³ with 99% confidence bounds of (20, 27) km and ($-2.8, -1.2$) km³, respectively. This is a lower bound of the source potency since all GNSS stations are grouped on one side of the source, so that the network has no sensitivity to NS source contraction. Both analyses indicate that vertical shrinking of the source is negligible, suggesting that the shrinking source does not coincide with the source of the VLPs.

## Interpretation of chain of events

In summary, we have identified two main stages of the unrest (Fig. 5). The first stage (Phases I–II), with the migration of an energetic seismic swarm from deep to shallow depth, is consistent
<figure><img src="image_19.png" /><figcaption>**Fig. 4** | Example of seismic signals and spectra recorded at the YTMZ station. **a**, 11 November 2018 VLP. Waveforms correspond to normalized vertical acceleration filtered in different frequency bands. **b**, Stack of normalized vertical spectra for 51 VLPs (indigo bars represent integer multiples of the dominant frequency (0.065 Hz, or 15.4 s) with a red bar for the 2.5 multiple). **c**, Anti-correlated waveforms for two VTs occurring 1.5 min apart on 20 August 2018 (time windows starting at 21:34:15.5 and 21:35:39.6 respectively). Waveforms of the first event (indigo) are compared to flipped waveforms of the second event (blue).</figcaption></figure>

with the propagation of a magmatic dyke through the whole crust. The second stage (Phases III–IV), with simultaneous VLP and VT activity and large deflation, suggests the evacuation of a sub-Moho, large reservoir and a multifaceted interaction between magma and host rock.

The energetic seismic swarms in Phase I marks, to our knowledge, the first case of vertical propagation of a dyke traced for >25 km from a deep magma reservoir to the surface, although dykes have been observed to propagate horizontally for longer distances²⁷,²⁸. Dyke-induced swarms of this magnitude are rare and attributed to large intruded volumes, which scale with the moment of induced events as a power law²⁹. According to this model, the cumulative moment of Phases I–II, $M_0 = 2 \cdot 10^{18}$ Nm, would correspond to a dyke volume of $5 \cdot 10^{-2}$ km³. Strike-slip focal mechanisms have been observed for other propagating dykes²⁸,³⁰⁻³³. Phase II involves a lateral dyke propagation, with VTs migrating 10–20 km towards South and reaching the seamount location²⁰. The relative VT quiescence

in Phase III is consistent with the establishment of an open pathway to the surface. High-rate crustal deformation sets on a week after the beginning of Phase III, consistent with the time needed for magma to widen its pathway by conduit erosion and establish a high flow rate³⁴. Taken together, these observations suggest that the end of Phase II may mark the onset of the submarine eruption. During Phase IV, VTs, VLPs and ground deformation appear interlinked mechanically. VTs plausibly trigger VLPs and accompany the slow evacuation of the magma reservoir (Fig. 2d,f).

VLPs may have different causes, but ringing events such as those observed at Mayotte have previously been explained by the resonance caused by slow standing waves trapped at the fluid-solid interface of a fluid-filled crack or conduit²,³⁵,³⁶. Here, their period and duration are especially long, which may be explained by the large size of the crack²⁰, the stronger host rock because of the depth and basalt as the fluid filling the crack³. We use analytical formulas³⁷ to constrain crack geometry based on the dominant
<figure><img src="image_20.png" /></figure>

**Fig. 5** | Summary schematic. A thin, subhorizontal reservoir feeds a dyke propagating first upward (Phase I) and later laterally southward (Phase II). The drainage of the reservoir's western edge triggers failure of the overlying rock, whose pressure pulses in turn trigger resonance in the reservoir (Phases III and IV).

VLP period, which we assume is the fundamental resonance mode. Fundamental periods of 15.2–15.6 s are possible for a wide range of lengths and thicknesses (Fig. 6). However, cracks with $L<8$ km have too little volume to expel more than 1.5 km³ of magma. A crack with $L>15$ km is unlikely, given the distribution of seismicity (Fig. 1). Smooth dominant period changes from 15.2 to 15.6 s and back may result first from the crack predominantly thinning and next shortening (Fig. 6). For example, path for $L=12$ km carries a volume loss of 1.6 km³. Thus, the change of dominant period may be entirely explained in terms of geometry changes consistent with the observed volume loss; variations in other parameters such as melt compressional wave velocity are also possible.

VLPs in combination with VTs have often been observed before and during caldera collapses¹⁸,³⁸⁻⁴². For example, VLPs occurred before and during the 2000 caldera collapse at Miyakejima, Japan, both with oscillating signals³⁹ and 20–50 s single pulse source time functions³⁸,⁴⁰. With a maximum duration of ~60 s they were considerably shorter than those at Mayotte; their modelled source mechanism also had a different geometry³⁸⁻⁴⁰. They were interpreted as the resonance of an axially symmetric structure in response to a shallower trigger³⁹ or as resulting from the intermittent sinking of a vertical piston into the magmatic chamber, causing its sudden volumetric change³⁸,⁴⁰. At Piton de la Fournaise, La Reunion island, 0.02–0.50 Hz VLPs have been attributed to repeated piston-like collapse⁴². The thrust mechanisms of typical VTs in Phase IV, inconsistent with the NE–SW transtensional background stress regime, require a strong stress perturbation, probably provided by the evacuation of the magma reservoir. Steep outward dipping faults are typically formed in the early stage of reservoir roof failure upon depletion⁴³,⁴⁴. They occur on newly formed, distributed faults (as suggested by their small magnitude and existence of several families), rather than on a ring fault, confirming that the faulting is weakening the overlying material rather than representing the slip of a coherent block. As they are close above the reservoir, they may exert an efficient pressure pulse and trigger waves at the fluid–solid interface traveling towards the opposite edge of the reservoir and back, thereby producing a ringing signal³. The observed anti-correlated VT pairs (Fig. 4c) may be explained by the reverse motion on a fault, similar

<figure><img src="image_21.png" /></figure>

**Fig. 6** | Resonance period of a magma-filled crack as a function of crack length and thickness. Thick black lines are periods of 15.2 and 15.6 s that comprise the dominant periods of VLPs analysed here. Contours of reservoir volumes are shown as thin black lines. Zoomed inset shows an example geometry pathway involving first thinning at constant $L$ and then shortening at constant $d$ for the example length of 12 km. $L$, crack length; $d$, crack thickness.

to observations during collapse events at calderas⁴⁵. The January 2018 VLP, observed before the dyke had formed, suggests that the oscillator is a portion of the deep reservoir rather than the feeder dyke, and that the reservoir then had a similar size and shape.

Our deformation models locate a volume loss of at least 1.3 km³ at 25–35 km depth below the point where the four GNSS displacement vectors converge, which coincides spatially with the downsag seismicity and whose size is consistent with the variations of the VLP dominant period. A simple hypothesis is that the easternmost edge of the reservoir corresponds with the location of the earliest seismicity burst, later seen to migrate upwards. Thus, two questions remain: assuming a 15 km long reservoir, why is volume loss at the western edge and not the centre of the crack, and why is the shrinking horizontal rather than vertical? Reconciling the evidence suggest that a dipping sill-shaped reservoir loosing buoyant magma would shrink at its deeper portion, here to the West, and one-sided drainage coupled mechanically with the reservoir's failing roof may explain the observed horizontal shrinking pattern. Alternatively, a more complex sill shape, with a vertical westernmost portion, could account for the crustal deformation pattern. Superposition of outward dipping thrust faulting are equivalent to a vertical positive CLVD⁴⁶, which may explain the model derived from the geodetic data.

In conclusion, our analysis suggests that a subhorizontal magma reservoir of up to 15 km in length lies between the newly discovered seamount²⁰ and Mayotte. Its exceptional depth of 30 ± 5 km makes it the deepest reservoir whose evacuation has been observed in ground displacement data. Preliminary estimates suggest that more than 3.4 km³ of magma effused at the seafloor²⁰, making it also the largest geophysically monitored submarine eruption to date. The mechanical processes activated offshore Mayotte represent a scaled-up version of a caldera formation process at its nucleating, downsag stage⁴³. Phase IV, currently in progress, involves the growth of deep faults at the western edge of the emptying reservoir. Observations and modelling of reservoir depletion⁴³ suggest
that, should it continue, fault growth will affect a progressively wider area and become shallower. A particular hazard is posed by a scenario where the outward dipping faults reach the ocean floor causing the entire block to collapse abruptly, and new normal faults begin to propagate upward, reaching closer to Mayotte. We estimate the volume evacuation threshold needed to trigger such a collapse, $V_{min}$, based on an equation developed for calderas⁴⁷:

$$V_{min} = \frac{10f \rho g h^2 r^2}{\kappa} \qquad (1)$$

where f is rock friction coefficient, ρ is average medium density, g is gravitational acceleration, h is reservoir depth, κ is magma bulk modulus, r is reservoir radius. Using f=0.6, ρ=2700 Kg m⁻³, h=25 km, r=7.5 km, κ=25 GPa we obtain a threshold volume of 230 km³, which is ~50 times the volume of the seamount on May 2019. This scenario appears remote at this stage, but critical to monitor any migration of seismicity or change of focal mechanisms, as well as better constrain reservoir and faults geometry and crustal properties.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-019-0505-5.

Received: 29 May 2019; Accepted: 14 November 2019; Published online: 6 January 2020

## References

1. Lemoine, A. et al. The volcano-tectonic crisis of 2018 east of Mayotte, Comoros Islands. Preprint at https://doi.org/10.31223/osf.io/d46xj (2019).

2. Ohiminato, T. et al. Waveform inversion of very-long-period impulsive signals associated with magmatic injection beneath Kilauea Volcano, Hawaii. *J. Geophys. Res.* **103**, 839–862 (1998).

3. Kumagai, H. & Chouet, B. A. Acoustic properties of a crack containing magmatic or hydrothermal fluids. *J. Geophys. Res. Solid Earth* **105**, 25493–25512 (2000).

4. Scrutton, R. A. et al. Constraints on the motion of Madagascar with respect to Africa. *Mar. Geol.* **43**, 1–20 (1981).

5. Phethean, J. et al. Madagascar's escape from Africa: a high-resolution plate reconstruction for the western Somali Basin and implications for supercontinent dispersal. *Geochem. Geophys. Geosyst.* **17**, 5036–5055 (2016).

6. Reeves, C. The position of Madagascar within Gondwana and its movements during Gondwana dispersal. *J. Afr. Earth Sci.* **94**, 45–57 (2014).

7. Geiger, M. et al. Reappraisal of the timing of the breakup of Gondwana based on sedimentological and seismic evidence from the Morondava Basin, Madagascar. *J. Afr. Earth Sci.* **38**, 363–381 (2004).

8. Coffin, M. & Rabinowitz, P. D. Reconstruction of Madagascar and Africa: evidence from the Davie Fracture Zone and western Somali Basin. *J. Geophys. Res.* **92**, 9385–9406 (1987).

9. Jokat, W. et al. Timing and geometry of early Gondwana breakup. *J. Geophys. Res.* https://doi.org/10.1029/2002JB001802 (2003).

10. Leinweber, V. T. & Jokat, W. The jurassic history of the Africa—Antarctica corridor—new constraints from magnetic data on the conjugate continental margins. *Tectonophysics* **530**, 87–101 (2012).

11. Mahanjane, E. S. A geotectonic history of the northern Mozambique Basin including the Beira High—a contribution for the understanding of its development. *Mar. Petrol. Geol.* **36**, 1–12 (2012).

12. Mahanjane, E. S. The Davie Fracture Zone and adjacent basins in the offshore Mozambique Margin—a new insights for the hydrocarbon potential. *Mar. Petrol. Geol.* **57**, 561–571 (2014).

13. Emerick, C. & Duncan, R. Age progressive volcanism in the Comores Archipelago, western Indian Ocean and implications for Somali plate tectonics. *Earth Planet. Sci. Lett.* **60**, 415–428 (1982).

14. Nougier, J. et al. The Comores archipelago in the western Indian Ocean: volcanology geochronology and geodynamic setting. *J. Afr. Earth Sci.* **5**, 135–145 (1986).

15. Michon, L. in *Active Volcanoes of the Southwest Indian Ocean* (eds Bachelery, P., Lénat, J.-F., Di Muro, A. & Michon, L.) 233–244 (Springer, 2016).

16. Audru, J. et al. Major natural hazards in a tropical volcanic island: a review for Mayotte Island, Comores archipelago, Indian Ocean. *Eng. Geol.* **114**, 364–381 (2010).

17. Delvaux, D. & Barth, A. African stress pattern from formal inversion of focal mechanism data. *Tectonophysics* **482**, 105–128 (2010).

18. Stamps, D. S. et al. Geodetic strain rate model for the east African rift system. *Sci. Rep.* **8**, 732 (2018).

19. Cesca, S. et al. *Seismic Catalogues of the 2018–2019 Volcano-seismic Crisis Offshore Mayotte* (GFZ Data Services, 2019).

20. Feuillet, N. MAYOBS1 Cruise, RV Marion Dufresne (Institut de Physique du Globe de Paris, 2019); https://doi.org/10.17600/18001217

21. Kumagai, H. et al. Magmatic dike resonances inferred from very-long-period seismic signals. *Science* **299**, 2058–2061 (2003).

22. Aster, R. et al. Very long period oscillations of Mount Erebus volcano. *J. Geophys. Res. Solid Earth* https://doi.org/10.1029/2002JB002101 (2003).

23. Talandier, J. et al. Unusual seismic activity in 2011 and 2013 at the submarine volcano rocard, Society hot spot (French Polynesia). *Geophys. Res. Lett.* **43**, 4247–4254 (2016).

24. Nikkhoo, M. et al. Compound dislocation models (CDMs) for volcano deformation analyses. *Geophys. J. Int.* **208**, 877–894 (2017).

25. Heimann, S. et al. Grond—a probabilistic earthquake source inversion framework. V. 1.0. GFZ Data Services, http://pyrocko.org/grond/docs/current/ (2018).

26. Davis, P. M. Surface deformation due to inflation of an arbitrarily oriented triaxial ellipsoidal cavity in an elastic half-space, with reference to Kilauea Volcano, Hawaii. *J. Geophys. Res. Solid Earth* **91**, 7429–7438 (1986).

27. Uhirai, K. & Toriyama, N. Meeting with the May 29, 2015 eruption of Kuchinoerabujima volcano. *Bull. Volcanol. Soc. Jpn.* **60**, 487–490 (2015).

28. Sigmundsson, F. et al. Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland. *Nature* **517**, 191–195 (2015).

29. White, R. & McCausland, W. Volcano-tectonic earthquakes: a new tool for estimating intrusive volumes and forecasting eruptions. *J. Volcanol. Geotherm. Res.* **309**, 139–155 (2016).

30. Hayashi, Y. & Morita, Y. An image of a magma intrusion process inferred from precise hypocentral migrations of the earthquake swarm east of the Izu peninsula. *Geophys. J. Int* **153**, 159–174 (2003).

31. Passarelli, L. et al. Stress changes, focal mechanisms, and earthquake scaling laws for the 2000 dike at Miyakejima (Japan). *J. Geophys. Res. Solid Earth* **120**, 4130–4145 (2015).

32. Ruch, J. et al. Oblique rift opening revealed by reoccurring magma injection in central Iceland. *Nat. Commun.* **7**, 12352 (2016).

33. Ágútsdóttir, T. et al. Strike-slip faulting during the 2014 Bárðarbunga-Holuhraun dike intrusion, central Iceland. *Geophys. Res. Lett.* **43**, 1495–1503 (2016).

34. Wadge, G. The variation of magma discharge during basaltic eruptions. *J. Volcanol. Geotherm. Res.* **11**, 139–168 (1981).

35. Ferrazzini, V. & Aki, K. Slow waves trapped in a fluid-filled infinite crack: implication for volcanic tremor. *J. Geophys. Res.* **92**, 9215–9223 (1987).

36. Chouet, B. Long-period volcano seismicity: its source and use in eruption forecasting. *Nature* **380**, 6572 (1996).

37. Maeda, Y. & Kumagai, H. A generalized equation for the resonance frequencies of a fluid-filled crack. *Geophys. J. Int.* **209**, 192–201 (2017).

38. Kumagai, H. et al. Very-long period seismic signals and caldera formation at Miyake Island, Japan. *Science* **293**, 687–690 (2001).

39. Kobayashi, T. et al. Very long period seismic signals observed before the caldera formation with the 2000 Miyakejima volcanic activity, Japan. *J. Geophys. Res. Solid Earth* https://doi.org/10.1029/2007JB005557 (2009).

40. Kobayashi, T. et al. Intermittent inflations recorded by broadband seismometers prior to the caldera formation at Miyake-jima volcano in 2000. *Earth Planet. Sci. Lett.* **357–358**, 145–151 (2012).

41. Munekane, H. et al. Mechanisms of step-like tilt changes and very long period seismic signals during the 2000 Miyakejima eruption: insights from kinematic GPS. *J. Geophys. Res. Solid Earth* **121**, 2932–2946 (2016).

42. Fontaine, F. R. et al. Very- and ultra-long-period seismic signals prior to and during caldera formation on La Reunion Island. *Sci. Rep.* **9**, 8068 (2019).

43. Acocella, V. Understanding caldera structure and development: an overview of analogue models compared to natural calderas. *Earth-Sci. Rev.* **85**, 125–160 (2007).

44. Levy, S. et al. Mechanics of fault reactivation before, during, and after the 2015 eruption of axial seamount. *Geology* **46**, 447–450 (2018).

45. Jónsdóttir, K. et al. Evaluating changes of the Bárðarbunga caldera using repeating earthquakes. *AGU Fall Meeting Abstracts*, 2017, S22C-07 (2017).

46. Cesca, S. & Heimann, S. in *Moment Tensor Solutions: A Useful Tool for Seismotectonics* (ed. D'Amico, S.) 163–181, (Springer, 2018).

47. Geshi, N. et al. Evaluating volumes for magma chambers and magma withdrawn for caldera collapse. *Earth Planet. Sci. Lett.* **396**, 107–115 (2014).

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© The Author(s), under exclusive licence to Springer Nature Limited 2020
Methods

Details about the methods used in this paper are given in the Supplementary information.

Data availability

Seismic data used in this study pertain to networks II⁴⁸, IU⁴⁹, GE⁵⁰, G⁵¹, PF⁵² and RA⁵³ and are available at IRIS (Incorporated Research Institutions for Seismology), GEOFON (GEO-FOrschungsNetz), ORFEUS EIDA (Observatories and Research Facilities for European Seismology—European Integrated Data Archive) and/or the Réseau Sismologique et Géodésique Français (French seismological and geodetic network, RESIF⁵³) web services. Geodetic data are available at the web facilities of the Nevada Geodetic Laboratory, at the University of Nevada, Reno⁵⁴.

Code availability

All codes used in this work are open source. The codes used to generate individual results are available through the contact information from the original publications. Requests for further materials should be directed to S.C.

References

48. *Global Seismograph Network - IRIS/IDA (GSN)* (Scripps Institution of Oceanography, International Federation of Digital Seismograph Networks, 1986); https://doi.org/10.7914/SN/I

49. *Global Seismograph Network—IRIS/USGS (GSN)* (International Federation of Digital Seismograph Networks, 1988); https://doi.org/10.7914/SN/IU

50. *The GE Seismic Network* (GEOFON Data Centre, 1993); https://doi.org/10.14470/tr560404

51. *GEOSCOPE* (Institut de Physique du Globe de Paris and Ecole et Observatoire des Sciences de la Terre de Strasbourg, 1982); https://doi.org/10.18715/geoscope.g

52. *Piton de la Fournaise Volcano Observatory Network (Reunion Island) (OVPF)* (FDSN, accessed 01 Aug 2019); https://www.fdsn.org/networks/detail/PF/

53. *RESIF-RAP Accelerometric Permanent Network* (RESIF, RESIF-RAP French Accelerometric Network, 1995); https://doi.org/10.15778/RESIF.RA

54. Blewitt, G., Hammond, W. C. & Kreemer, C. Harnessing the GPS data explosion for interdisciplinary science. *Eos* https://doi.org/10.1029/2018EO104623 (2018).

Acknowledgements

The discovery of the underwater volcano is a result of the MAYOBS 1 campaign that took place from 2 to 19 May, 2019 aboard the Marion Dufresne oceanographic vessel. This campaign was conducted by several French research institutions and laboratories (IPGP/CNRS/BRGM/IFREMER/IPGS) as part of a CNRS-INSU programme²⁰. We thank CNRS-INSU for making public the location and size of the discovered volcano. M.I. thanks H. Sudhaus for her valuable contribution and acknowledges funding by the German Research Foundation DFG through an Emmy Noether Young Researcher Grant (no. 276464525). G.P. is funded by the German Research Foundation DFG project (no. 362440331), a subproject of "SPP 2017: Mountain Building Processes in 4D" (project number no. 313806092). We thank T. James and T. Davis for revising the English text.

Author contributions

S.C. coordinated this project, conceived the manuscript and figures, and analysed, modelled and interpreted local accelerometric data. J.L. performed the VLP and array analysis. S.H. and H.R. performed the MT inversion for VT and VLP events. M.I., M.N. and L.P. analysed and modelled deformation data. G.P. assessed the seismic data quality. E.R., F.C. and T.D. contributed to the interpretation of results and discussion section. S.C., E.R., F.C. and T.D. drafted the manuscript. All authors reviewed the manuscript.

Competing interests

The authors declare no competing interests.

Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41561-019-0505-5.

**Correspondence and requests for materials** should be addressed to S.C.

**Peer review information** Primary Handling Editors: Stefan Lachowycz; Melissa Plail.

**Reprints and permissions information** is available at www.nature.com/reprints.