# Seismicity near Mayotte explained by interacting magma bodies: Insights from numerical modeling

Clément de Sagazan, Lise Retailleau, Muriel Gerbault, Aline Peltier, Nathalie Feuillet, Fabrice J. Fontaine, Wayne C. Crawford

► To cite this version:

Clément de Sagazan, Lise Retailleau, Muriel Gerbault, Aline Peltier, Nathalie Feuillet, et al.. Seismicity near Mayotte explained by interacting magma bodies: Insights from numerical modeling. Journal of Volcanology and Geothermal Research, 2024, 446, 10.1016/j.jvolgeores.2023.107985. insu-04462198

HAL Id: insu-04462198

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

Submitted on 22 Nov 2024

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Seismicity near Mayotte explained by interacting magma bodies: Insights from numerical modeling

de Sagazan Clément <fnref n="1" />,*, Retailleau Lise <fnref n="1" />,<fnref n="2" />, Gerbault Muriel <fnref n="3" />, Peltier Aline <fnref n="1" />,<fnref n="2" />, Feuillet Nathalie <fnref n="1" />, Fontaine Fabrice J. <fnref n="1" />,<fnref n="2" />, Crawford Wayne C. <fnref n="1" />

## Abstract :

Mayotte island experienced a large volcanic eruption 50 km offshore in 2018–2021, creating the submarine volcano “Fani Maoré”. The eruption was accompanied by intense seismicity at mantle depths (20–45 km), divided into a “proximal” and a “distal” cluster centered 10 and 30 km east from the island, respectively. Previous studies suggest that two separate magma reservoirs may lie at the top and bottom of the proximal cluster. Here, we assess whether two reservoirs are a mechanically viable explanation for the proximal cluster's truncated conical shape. We developed finite-element models of pressurized magma reservoirs in a 2D axisymmetric domain, modeling the reservoirs as compliant elastic ellipsoids embedded in an elastoplastic host rock. We find that, at these depths, extremely low friction is required to generate failure at realistically low reservoir pressures. This implies in turn that mechanical weakening must occur at these depths. The weakening could be induced by fractures or pore fluid overpressure in the volcanic system. We find that two superimposed reservoirs can generate a plastic domain between them, if they are spatially close enough. Several reservoir geometries (from spherical to sill-like) are plausible. A conical fracture domain is more likely to appear for reservoirs with opposite pressure loads (i.e. one inflating, one deflating). Given the geometrical match with the proximal seismicity cluster at Mayotte, we suggest that the shallower (Moho-depth) reservoir is inflating, creating a potential hazard for Mayotte island.

## Highlights

* Two interacting magma reservoirs can explain the earthquake distribution near Mayotte. A 20-km deep reservoir near the island may have inflated during the recent eruption. Shear patterns between modeled magma reservoirs depend on the applied pressure loads. The occurrence of mantle-deep earthquakes near Mayotte may indicate mantle weakening.

<footnote n="1">¹ Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, 1 rue de Jussieu, Paris 75005, France</footnote>

<footnote n="2">² Observatoire Volcanologique du Piton de la Fournaise, Institut de Physique du Globe de Paris, 14 RN3 km 27, La Plaine des Cafres, 97418 La Réunion, France</footnote>

<footnote n="3">³ GET (UMR 5563, IRD, Université P. Sabatier, CNRS, CNES, OMP), 14 avenue Edouard Belin, Toulouse 31400, France</footnote>

<footnote>* Corresponding author : Clément de Sagazan, email address : desagazan@ipgp.fr</footnote>
**Keywords** : Mayotte, Fani Maoré, Magma reservoir failure, Numerical modeling, Finite elements, Non-associated plasticity
# 1. Introduction

Volcanic eruptions are surface manifestations of the Earth's internal dynamics. Telltale signs associated with volcanic activity include ground deformation (e.g. Xu et al., 2016; Peltier et al., 2022) or earthquakes (e.g. Levy et al., 2018; Agústsdóttir et al., 2019) linked with host rock fracturing by fluid motion or accumulation at depth. Deformation and earthquake observations can therefore shed light on the subsurface structure and dynamics of magmatic systems, and are key to monitoring volcanic activity.

Since May 2018, the French island of Mayotte, in the southwest Indian Ocean, has experienced seismic activity without known historical precedent (Figure 1, Cesca et al., 2020; Lemoine et al., 2020). In addition to the emplacement of the new submarine Fani Maoré volcanic edifice east of the island in 2018 (Feuillet et al., 2021), thousands of earthquakes have been detected (REVOSIMA, 2023). Despite the apparent end of the eruption in early 2021 (Berthod et al., 2022), seismicity persists at a low but steady level (REVOSIMA, 2023). The seismicity, exceptional for its early intensity, its depth and its longevity, reveals the presence of an active magmatic system near Mayotte (Feuillet et al., 2021; Berthod et al., 2021a), calling for a reevaluation of the regional volcano-seismic hazard.

In the context of this recent seismo-volcanic activity, interpretative models of the magma feeding system have been proposed, based on geophysical and geological observations (details in Section 2). These models suggest that a ≈40-km deep magma storage zone was drained eastward by an ascending dyke or conduit (e.g., Feuillet et al., 2021; Lavayssière et al., 2022; Mercury et al., 2023). Nonetheless, this system's detailed structure and dynamics remain poorly understood, and could involve multiple magma reservoirs (Berthod et al., 2021a,b; Lavayssière and Retailleau, 2023) and a mush zone (Mittal et al., 2022).

The seismicity distribution attracts particular attention. Most earthquakes occur at mantle depths (22-45 km), but close to land (5-15 km), beneath a possible former caldera ringed with geologically recent eruptive products (Puzenat et al., 2022; REVOSIMA, 2023). A horseshoe-shaped volcanic structure on its southern rim (the “Fer à Cheval”) is currently experiencing degassing (Feuillet et al., 2021). This highlights the threat of a
connection between deep magmatic structures and outlets near or on the island, especially since degassing also occurs onshore (Liuzzo et al., 2022). Understanding the structures and processes generating this seismicity should help to assess the volcanic hazard affecting nearby populated areas.

In this study, we use elasto-plastic numerical models to investigate the origin of the seismicity near Mayotte. Our goal is to gain first-order mechanical understanding of magmatic structures and their interaction which could explain the observed seismicity distribution.

## 2. Context

### 2.1. Geodynamic setting

The Comoros archipelago is a NW-SE trending alignment of four volcanic islands: Grande Comore, Mohéli, Anjouan and Mayotte, lying north of the Mozambique channel and south of the Somali basin (Figure 1, Michon, 2016). The area features seamounts, volcanic ridges, lava flows and hundreds of volcanic cones (Feuillet et al., 2021; Thinon et al., 2022).

The Comoros mark an immature transtensional plate boundary between the Lwandle and Somali plates (Famin et al., 2020; Feuillet et al., 2021; Thinon et al., 2022). Together with a negative S-wave velocity anomaly in the upper mantle (Mazzullo et al., 2017), this suggests that volcanism in the region may result from an interaction between lithospheric-scale tectonics and a mantle plume.

Volcanism at Mayotte initiated 26-27 Ma (Masquelet et al., 2022). The youngest traces of onshore volcanism in Mayotte are 4 ka old (Zinke et al., 2003). Before 2018, the only historically known active volcano in the Comoros was Karthala on Grande Comore. Hence, the region was considered seismically quiet. Only about 30 earthquakes of $M_w \ge 4$ were recorded in the Comoros between 1976 and 2018 (Famin et al., 2020).
## 2.2. Seismo-volcanic activity since 2018

The Mayotte seismo-volcanic crisis began on May 10, 2018. Over the first month, hundreds of earthquakes were felt, with magnitudes up to 5.9 (Cesca et al., 2020;

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

Figure 1: Geophysical observations of the 2018 eruption near Mayotte. Top: seismicity distribution (Lavayssière et al., 2022) and deformation sources (Peltier et al., 2022). The arrows show the orientation of the cross-sections below. Red symbol: location of Fani Maoré volcano. $M_L$: local magnitude. Contours: 0.2, 0.4, 0.6, 0.8 and 1 isovalues of the probability density function of GNSS-derived deflation sources (Peltier et al., 2022) for 2018 (red), 2019 (green) and 2020 (cyan), with the most opaque areas representing the highest densities. Red frame on top right globe: location of the Comoros archipelago. The square frame beneath the globe details this region, showing 1: Grande Comore. 2: Mohéli. 3: Anjouan. The red box in this plot outlines the bounds of the top left map. Bottom: 2-km-wide cross-sections of the proximal cluster. The horizontal distances are the distances along the cross-section axes. Contours: Vp/Vs ratio from seismic tomography (Foix et al., 2021). Areas below 45 km depth (in gray) are unresolved.
Lemoine et al., 2020; Mercury et al., 2023). In mid-June 2018, the number and magnitude of earthquakes decreased, and VLP (Very Long Period) events began to be observed (Cesca et al., 2020; Laurent, 2023). VLP events, thought to be generated by resonating fluid-filled cavities, are relatively common during volcanic episodes (Chouet, 1996). In July 2018, the island began to subside eastward. This subsidence, reaching in places 10-20 cm, is best modeled by a ≈ 40-km-deep deflating source east of Mayotte (Figure 1, Lemoine et al., 2020; Peltier et al., 2022; REVOSIMA, 2023). The VLP events and the subsidence are thought to mark the beginning of the offshore eruption (Cesca et al., 2020; Lemoine et al., 2020).

The new Fani Maoré volcanic edifice was discovered during the first MAYOBS cruise in May 2019 (Rinnert et al., 2019; Feuillet et al., 2021). The edifice lies 50 km east of Mayotte, at 3500 m depth, (Figure 1). The 820-m high volcano and associated lava flows have a volume of 6.55 km³ (REVOSIMA, 2023), as of 2021, making the 2018-2021 episode one of the largest basaltic eruptions ever observed and documented (Feuillet et al., 2021). Since the first MAYOBS cruise, repeated oceanographic cruises continue to monitor the volcanic edifice, with the last active lava flow detected in January 2021 (Rinnert et al., 2019; Berthod et al., 2022).

The abundant seismicity since 2018 sheds light on the architecture of the Fani Maoré feeding system. The seismicity is divided between two clusters, named “distal” and “proximal”, based on their distances from Mayotte (Saurel et al., 2021). Figure 1 shows the relocated seismicity catalog from 25 February 2019 to 9 May 2020, built by Lavayssière et al. (2022).

Most if not all of the first month's earthquakes occurred in the distal cluster (Mercury et al., 2023), whose base lies at 40 km depth, 30 km east of Mayotte, and which extends upwards and southeast towards the Fani Maoré edifice (Figure 1). Most active in May-June 2018, this cluster is thought to mark host-rock fracturing prior to magma ascension towards the seafloor (Cesca et al., 2020; Mercury et al., 2023).

The proximal cluster initiated in late September 2018 (Mercury et al., 2023). Mostly located between depths of ≈ 22 km and 45 km, it is centered 10 km east of
Mayotte's Petite-Terre island (Figure 1). This cluster is less understood than the distal one. To a first order, it is roughly axisymmetrical. In cross-section, it has an “hourglass” shape (Lavayssière et al., 2022), with a conical bottom part (below ≈ 32 km depth) and an ellipsoidal top part. Most of the proximal earthquakes are volcano-tectonic (VT), but VLP and Long Period (LP) events are also located in this cluster (Laurent, 2023; Retailleau et al., 2022).

## 2.3. Possible interpretations of the proximal earthquake distribution

At first glance, the conical shape of the proximal seismicity cluster evokes outward dipping ring faults such as those observed during caldera collapse events at Axial Seamount on the Juan de Fuca Ridge (Levy et al., 2018) and at Bárðarbunga in Iceland (Agústsdóttir et al., 2019). However, there is no seafloor geomorphological evidence of a caldera collapse associated with the Fani Maoré eruption.

The volcanic contexts mentioned above feature seismicity at much shallower depths than in the context of Mayotte (0-2 km for Axial Seamount, 0-6 km for Bárðarbunga). McTigue (1987) linked, based on analytical work, caldera collapse above an inflating spherical cavity to circumferential stress acting along its walls and the presence of a free surface. More generally, many numerical and analogical studies have modeled the development of ring fracture zones linking shallow magma reservoirs to the surface (e.g., Gudmundsson, 2007; Holohan et al., 2013; Gerbault et al., 2012, 2018; Cabaniss et al., 2020).

Deeper reservoirs have received less attention. Grosfils et al. (2015) show that, when a reservoir's depth is large enough compared to its size (usually depth/size ratios of 2-3), circumferential (i.e. “ring”) roof failure becomes less favorable than vertical dyking from the reservoir's apex or lateral sill emplacement from its sides (depending on reservoir geometry). Interpreting the proximal cluster as a marker of “ring-fault”-associated collapse is thus challenging because (1) the necessary failure pattern is unlikely at such depths and (2) there is little to no seismicity between the seafloor and 20 km depth (Figure 1).

Lavayssière et al. (2022) propose an alternative explanation, in which the proximal cluster represents a conical fault zone connecting two superimposed reservoirs. They
suggest that its top outlines an ovoid aseismic zone. Seismic tomography by Foix et al. (2021) also shows two roughly ellipsoidal areas with high Vp/Vs at the cluster's top and bottom boundaries (Figure 1). Further constraints on reservoir depths and dimensions are detailed in Section 3.2.

Few studies have investigated the mechanical relationships between neighboring reservoirs. Gonnermann et al. (2012) explain the apparent coupling between Mauna Loa and Kilauea volcanoes (Hawaii) by modeling pore pressure diffusion within the melt accumulation layer. Pascal et al. (2014) model magma bodies as pressurized cavities within an elastic domain to assess the impact of magma reservoir interactions on surface displacement. Albino and Sigmundsson (2014) use a similar framework but study how a neighboring magma body locally alters the failure threshold on a reservoir's wall. Karaoğlu et al. (2020) present elastic models of the Karliova Triple Junction (eastern Turkey) featuring three neighboring magma bodies, and predict likely areas of failure and magma transfer based on stress distributions within the host rock. Because of the complexity of this area, the authors choose a multi-layered and highly heterogeneous domain with numerous preexisting faults. Their modeled magma bodies are also much closer to one another (≈ 2 km) than the anomalies imaged by Foix et al. (2021) (10-15 km) near Mayotte. Their results thus seem difficult to relate directly to our case study.

Another challenge concerning the proximal cluster is to explain the occurrence of seismic failure at such great depths. As depth increases, so does the confining pressure, which tends to inhibit tensile or shear rock failure (Grosfils, 2007; Gerbault, 2012). The role of internal fluid pressures or weak chemical phases is often invoked to reduce effective yield strengths (e.g. Rubey and King Hubbert, 1959; Brantut et al., 2016).

Our models therefore address two main questions: (1) are two interacting magma reservoirs a mechanically viable explanation for the occurrence and shape of the proximal seismic cluster observed in Mayotte?; and (2) what conditions can lead to reservoir wall failure at these depths?
<figure><img src="image_3.png" /><figcaption>Figure 2: Numerical setup of our models. A: Boundary conditions and parameter definitions (see Table 1 for more details). B: Overview of the used mesh. Minimum element size is 150 meters on the reservoir walls.</figcaption></figure>

# 3. Methods

## 3.1. *Numerical setup*

We perform finite element calculations using the COMSOL Multiphysics® software (v6.1). We use the Structural Mechanics module of COMSOL which solves for the conservation of mass and momentum on an elasto-plastic material; the equations and the method are detailed in Appendix A and Appendix B. We choose a cylindrical coordinate system, because we assume that the structures generating the axisymmetric proximal cluster are also axisymmetric. We define a 200-km-by-200-km 2D axisymmetric domain, with $r = 0$ corresponding to the symmetry axis (Figure 2). The top of the modeled domain is a free surface ($z = 0$, representing the seafloor). Considering the depths of interest, we neglect the presence of seawater, which is discussed in Section 3.3. The bottom boundary is a no-slip boundary and the lateral boundaries have a roller (or “free-slip”) boundary condition. Gravity is applied ($g = 9.81 \text{ m.s}^{-2}$) and the domain is lithostatically loaded.

The domain includes the host rock as well as one or two spheroidal subdomains at the symmetry axis, representing magma reservoirs (Figure 2). Along the reservoir walls, we apply a normal stress $\Delta P$, representing the reservoir's pressure difference with respect to lithostatic pressure. A negative value represents an underpressure (deflation) and a
positive value an overpressure (inflation). Reservoir pressurizations are noted $ΔP_{up}$ and $ΔP_{low}$ for the upper and lower reservoir, respectively.

We build our mesh using triangular elements. The resolution is high near the reservoirs (150 m) and coarsens progressively towards the domain's edges (up to 7500 m). This allows us to capture fine structures while keeping computation times reasonable over a domain that is large enough to avoid border effects.

## 3.2. Geometrical parameters

We define magma reservoir depths and dimensions based on geophysical and geological data acquired near Mayotte. Table 1 provides numerical values.

In addition to the proximal cluster's shape, there are several lines of evidence for magma reservoirs near its top and bottom. LP and VLP earthquakes, generally associated with resonating fluid-filled areas (Chouet, 1996), seem to be located within a VT event gap at the axis of the proximal cluster (Retailleau et al., 2022; Laurent, 2023). Analysis of rock samples from the Fani Maoré and Fer a` Cheval areas suggest that a main storage reservoir lies between 37 and 48 km depths and that several differentiated melt lenses may be stagnating near the Moho, at 15-20 km depths (Berthod et al., 2021a,b). Foix et al. (2021) image two areas with high Vp/Vs (Figure 1) above and below the proximal cluster, which they interpret as melt-rich zones. These anomalies range from 22 to 30 km depth, and from 38 to 50 km depth, with respective volumes of 205 km³ and 234 km³.

To account for these observations, our simulations feature two magma reservoirs. We also provide a numerical benchmark with a single reservoir in Section 4.1. Based on seismic tomography and geochemical data, we choose a central depth of 40 km for the deep reservoir, and a range of depths from 20 to 30 km for the shallower one. When testing spherical reservoirs, we define radii of 3 km and 4 km for the upper and lower reservoir, respectively. This yields initial volumes of 113 km³ and 268 km³, respectively. We examine the effects of reservoir geometries using thickness/width ratios of 1 (spheres), 1/2 (oblate) and 1/4 (sill-like).
<table><caption>Table 1: List of numerical parameters. For C and φ, parametric tests are presented in Figure S1.</caption><thead><tr><th>Parameter description</th><th>Parameter name</th><th>Value(s)</th></tr></thead><tbody><tr><td>Upper reservoir wall pressure</td><td>ΔP<sub>up</sub></td><td>—</td></tr><tr><td>Lower reservoir wall pressure</td><td>ΔP<sub>low</sub></td><td>—</td></tr><tr><td>Pressure ratio</td><td>ΔP<sub>up</sub>/ΔP<sub>low</sub></td><td>-1 to 1 (increment of 0.1)</td></tr><tr><td>Upper reservoir central depth*</td><td>d<sub>up</sub></td><td>20, 22.5, 25, 27.5, 30 km</td></tr><tr><td>Lower reservoir central depth*</td><td>d<sub>low</sub></td><td>40 km</td></tr><tr><td>Upper reservoir width</td><td>w<sub>up</sub></td><td>6 km</td></tr><tr><td>Lower reservoir width</td><td>w<sub>low</sub></td><td>8 km</td></tr><tr><td>Upper reservoir thickness</td><td>h<sub>up</sub></td><td>6, 3, 1.5 km</td></tr><tr><td>Lower reservoir thickness</td><td>h<sub>low</sub></td><td>8, 4, 2 km</td></tr><tr><td>Reservoir separation*</td><td>(d<sub>low</sub> - h<sub>low</sub>/2) - (d<sub>up</sub> + h<sub>up</sub>/2)</td><td>13, 10.5, 8, 5.5, 3 km</td></tr><tr><td>Bedrock Young's modulus</td><td>E<sub>b</sub></td><td>150 GPa</td></tr><tr><td>Reservoir Young's modulus</td><td>E<sub>m</sub></td><td>0.1 GPa</td></tr><tr><td>Density</td><td>ρ</td><td>3300 kg.m<sup>-3</sup></td></tr><tr><td>Poisson's ratio</td><td>ν</td><td>0.25</td></tr><tr><td>Cohesion</td><td>C</td><td>10 MPa</td></tr><tr><td>Friction angle</td><td>φ</td><td>1°</td></tr><tr><td>Dilatancy angle</td><td>ψ</td><td>0</td></tr></tbody></table>

*“Reservoir separation” is the distance from the upper reservoir’s base to the lower reservoir’s apex. The central depths shown here are those tested for the spherical reservoir cases. For non-spherical cases, they were modified so that reservoir separations remained identical.

### 3.3. Mechanical parameters

We define the elastic parameters throughout the modeled domain (Young’s modulus E, Poisson’s ratio ν and density ρ). We approximate the reservoirs as purely elastic (following Gerbault et al., 2012; Got et al., 2013) and the host rock as elasto-plastic with prescribed failure parameters (cohesion C and friction angle φ).
To simplify our analysis, we choose a uniform density $ρ = 3300$ kg.m$^{-3}$, corresponding to mantle peridotites. This “replaces” 17 km of crust ($ρ = 2800$ kg.m$^{-3}$) with mantle, adding 83 MPa to the confining pressure. We do not consider the buoyancy effects of less dense magma reservoirs on the stress field (see Discussion) and neglect the weight of seawater, which subtracts ≈ 35 MPa on top of the domain. These approximations yield a ≈ 50 MPa excess in confining pressure at mantle depths. Compared to the confining pressures of 800 MPa at 25 km depth and 1.3 GPa at 40 km depth, these simplifications yield relative errors of 6% and 3.5%, respectively.

Based on our chosen density and on the velocity model of Lavayssière et al. (2022), we calculate a bedrock Young’s modulus $E_b = 150$ GPa, which is a reasonable value for the upper mantle (e.g., Afonso et al., 2005). Poisson’s ratio $v$ is set to 0.25 for the reservoirs and the host rock. Estimating an average magma reservoir Young’s modulus $E_m$ is difficult due to internal heterogeneities linked with the presence of several physical phases (crystal layering, volatiles) that can evolve transiently before and after a volcanic eruption. Several studies argue for complex poro-visco-elastic behavior (e.g., Mittal and Richards, 2019; Liao et al., 2021; Alshembari et al., 2022, 2023). Here, we neglect these internal variations and choose an homogeneous constant elastic behavior.

We test different values of $E_m$ for a deflating reservoir in a fully elastic model (Figure 3A) and then fix it for subsequent model cases, based on how well the resulting surface displacement fits analytical solutions for pressurized spherical reservoirs (McTigue, 1987). We find that more compliant magma reservoirs produce better fits. Ultimately, we choose $E_m = 0.1$ GPa: smaller values require excessively small numerical time-steps for the elasto-plastic simulations. Alternatively, we could have chosen empty (not meshed) magma reservoirs, similarly to other studies (e.g., Karaoğlu et al., 2020; Zhan and Gregg, 2019); our choice stems from faster numerical convergence when considering a meshed magma reservoir. Gerbault (2012) compared both approaches, and a numerical benchmark is shown in Appendix C.

Elasto-plastic behavior in COMSOL is dealt with in the following manner (see Appendix B). Elastic materials follow Hooke’s law (e.g., Jaeger et al., 2007). Plastic behavior is prescribed with a Drucker-Prager failure criterion, a smoothed approximation
of the Mohr-Coulomb criterion in 3D, and follows a non-associated plastic flow rule
appropriate for rocks and granular materials (dilatancy angle ψ = 0, Vermeer and de Borst,
1984). Plastic behavior is irreversible and non-associated plastic flow rules are difficult to
compute numerically, because they require drastically small time-steps (Gerbault et al.,
1998). Reservoir wall pressures must therefore be imposed with small increments, as a
succession of quasi-static steps (e.g., Gerbault, 2012; Novoa et al., 2022). We choose an
initial pressure step of 1 MPa. When the solver fails to converge, it automatically reduces
the pressure step by down to seven orders of magnitude. Nevertheless, in our simulations,
the solver fails to converge shortly after a specific pressure step, due to the large amount of
deformation (e.g., mesh distorsion) induced by plastic behavior (see Results). Note that our
modeled domain does not present initial strength heterogeneities; the role of preexisting
structures and other factors that would favor rupture at specific locations are discussed later
in Section 5.2.

Ground displacement inversions on Mayotte suggest that a 40-km-deep deflation
source lies near the distal cluster (Figure 1, Lemoine et al., 2020; Feuillet et al., 2021;
Peltier et al., 2022). However, the source's position is poorly constrained due to the non-
optimal GNSS network disposition and the simplicity of the purely elastic rheology.
Multiple sources and deformation occurring at the proximal cluster cannot be excluded.
Magma withdrawn from the area of the proximal cluster could have replenished a drained
reservoir further east (e.g., Mittal et al., 2022). We thus assume that the ≈ 40-km-deep
lower reservoir is deflating, i.e. ΔP<sub>low</sub> < 0. For the upper reservoir, we test positive, null
and negative wall pressures (ΔP<sub>up</sub>). Without constraint on the actual evolution of pressures
within the two hypothetical reservoirs, we assume simultaneous variations of magma
pressure in both reservoirs.

$\Delta P_{up}$ and $\Delta P_{low}$ are imposed progressively such that the $\Delta P_{up}/\Delta P_{low}$ ratio remains
<figure><img src="image_4.png" /><figcaption>Figure 3: Benchmark with a single reservoir. A: vertical displacement generated by a 40-km deep spherical reservoir of radius 4 km, underpressurized by 1.5 GPa (hence deflating by 5 km³). Simulations are performed in an elastic bedrock with E<sub>b</sub> = 150 GPa, and for different values of E<sub>m</sub>. The inset shows misfits with the analytical solution by McTigue (1987), in log scale for both axes. B: Shear stress field obtained with a deflating reservoir in an elastic medium. The color levels show our numerical solution. The contours show the analytical solution by McTigue (1987), in GPa. C: plastic shear strain obtained with a deflating reservoir in an elasto-plastic medium, with C = 10 MPa and φ = 1°. σ<sub>1</sub>, σ<sub>2</sub> and σ<sub>3</sub>: principal stresses (from most to least compressive).</figcaption></figure>

constant throughout the simulation.

GNSS inversion by Peltier et al. (2022) propose 5 km³ of volume loss at depth between 2018 and 2019, while Berthod et al. (2021a) suggest 4.8 km³ based on petrology. Throughout the simulations, we also survey the volume variation of our modeled magma reservoirs (noted ΔV<sub>up</sub> and ΔV<sub>low</sub> for the upper and lower reservoir, respectively) by integrating the normal displacement at the reservoir walls. We use the 5-km³ deflation as an upper bound for the volume loss at the proximal cluster (see Discussion).

We represent our results graphically using the “Von Mises” or “shear” stress (second invariant of the deviatoric stress tensor, here noted τ) and plastic shear strain (here noted ε<sup>pl</sup><sub>II</sub>) distributions (more details in Appendix A.2 and Appendix A.5). Our goal is to produce modeled distributions of τ and/or ε<sup>pl</sup><sub>II</sub> that match the conical shape of the proximal cluster at the considered depths.
# 4. Results

## 4.1. Preliminary tests with a single reservoir

We first perform numerical benchmarks using a single spherical reservoir in an elastic bedrock. This framework, commonly used in volcanology (e.g., McTigue, 1987; Peltier et al., 2022), is a first-order approximation of magma reservoir behavior. Figure 3B shows the result for a 40-km-deep spherical reservoir with a 4-km radius. Segall (2010) gives the equation linking reservoir wall pressure and volume variation for a spherical cavity of radius a:

$$ \Delta V = \frac{\pi a^3}{G} \Delta P, \quad \text{with} \quad G = \frac{E}{2(1 + \nu)} \qquad (1) $$

Assuming $\Delta V = -5$ km³ for Mayotte (Berthod et al., 2021a; Peltier et al., 2022), we obtain $\Delta P = -1.5$ GPa, greater than the 1.3 GPa confining pressure at 40 km depth. The reservoir dimensions are uncertain and, although our chosen radius seems to match the observations, Equation 1 shows that the relationship between $\Delta P$ and $\Delta V$ strongly depends on reservoir size. It also implies that the decrease of Young's modulus at high temperatures (e.g., Cabaniss et al., 2020) yields a smaller $\Delta P$ to match a similar volume loss. Nevertheless, a wall pressure of at least several hundred MPa is still needed to deflate the reservoir by 5 km³. While our modeled distribution of $\tau$ matches the solution by McTigue (1987) (Figure 3B), the elastic framework thus requires an unrealistic pressure in order to reproduce the observed volume variation.

Furthermore, given the reservoir's small radius/depth ratio, the distribution of $\tau$ is concentric. Predicting conical "ring"-like failure would require stress concentrations on the sides of the reservoir (McTigue, 1987; Grosfils et al., 2015), which do not appear here. Purely elastic simulations are thus insufficient to predict both the proximal earthquake cluster's occurrence and its spatial distribution.

While elastic simulations can predict where host rock failure may initiate, they generally fail to compute its evolution (Gerbault et al., 1998). Therefore, we perform a second test using the same reservoir dimensions but in an elasto-plastic bedrock. We consider self-consistent plasticity to assess whether a single reservoir can reproduce propagating fractures or shear zones that match the shape of the proximal cluster.
We first explore different values of the failure parameters (Figure S1). At the considered depths, cohesion (C) has little impact, so we use a standard value of C = 10 MPa. The critical pressure for reservoir wall failure is mainly driven by the friction angle (φ). Figure S1 shows that, unless we set φ below 5° (as opposed to the standard 30° for undamaged rocks), shear failure at 40 km depth requires a ΔP of several hundred MPa. Although we lack constraints on the actual reservoir pressure, such large values are likely unrealistic (see Section 5.2). Therefore, in the rest of the study, we show simulations with φ = 1°, representing an end-member case with extremely low friction; several studies argue for such a low value, referred to as an “effective” friction (e.g., Gerbault et al., 2022, and references therein).

Figure 3C shows the $\epsilon^{pl}_{II}$ field generated by a deflating spherical reservoir with C = 10 MPa and φ = 1°. Volume loss of 5 km³ is achieved for ΔP = -180 MPa, ≈ 8 times less than for an elastic bedrock. Once the failure threshold is reached, the amount of strain generated for a given pressure is larger than for the elastic case, in agreement with, e.g., Gerbault et al. (2012). The modeled plastic domain is concentric around the reservoir wall, consistent with Gerbault et al. (2018) for a spherical reservoir in 3D. At the considered pressure, no conical fracture zone appears.

We also perform tests with an oblate and a sill-like reservoir (Figure S2). The nonsphericity of the magma source generates stress and strain concentrations at the reservoir’s tips, i.e. the areas of greatest curvature, in agreement with the literature (Gudmundsson, 2007; Grosfils et al., 2015; Gerbault et al., 2018; Zhan and Gregg, 2019; Sigmundsson et al., 2020). Nonetheless, we observe that the generated plastic domains remain mostly diffuse, and do not match the shape of the observed proximal cluster near Mayotte. Even in an elasto-plastic framework, a single reservoir thus seems unlikely to generate the observed earthquake distribution.

## 4.2. A reference model case

We now present a reference model with two spherical reservoirs within an elasto-plastic bedrock. In this reference case, only the deepest reservoir is pressurized ($\Delta P_{low} < 0$ and $\Delta P_{up}/\Delta P_{low} = 0$).
<figure><img src="image_5.png" /><figcaption>Figure 4: Evolution of the strain field during an elasto-plastic simulation with two spherical reservoirs. Only the lower reservoir has an applied underpressure. A-C: key snapshots of the simulation. $σ_1$, $σ_2$ and $σ_3$: principal stress orientations (from most to least compressive). D: volume variation for both reservoirs as a function of lower reservoir wall pressure during the simulation. A, B and C correspond to the snapshots on the left. The dashed black line highlights the 5 km³ volume loss. Stage 1-4: see main text.</figcaption></figure>

Figure 4 shows the evolution of $e^{p_{II}}$ and of the reservoirs' change in volume during the simulation. This evolution can be divided into four stages:

* Stage 1: the lower reservoir deflates elastically. $ΔV$ is proportional to $ΔP$, consistent with Equation 1.

* Stage 2: the lower reservoir wall reaches failure at $ΔP_{low} \approx -40$ MPa, in agreement with our analytical predictions (Appendix A.3). A concentric plastic domain forms around it and progressively expands (Figure 4A).

* Stage 3: the plastic shear domain reaches the upper reservoir. A zone of localized plasticity progressively appears in between them (Figure 4B).

* Stage 4: the plastic domain connecting both reservoirs is now a well-defined conical fracture zone that concentrates most of the plastic strain (Figure 4C). Volume variation sharply increases for both reservoirs (Figure 4D).
This reference case shows that adding an area of greater compliance above a pressurized reservoir can enhance stress localization and favor the development of conical shear bands. Based on the orientation of the principal stresses and on the displacement field (Figure S3), the produced plastic domain accommodates inverse displacement, and delimits a conical subsiding area between the reservoirs. In this case, the establishment of a well-defined conical shear zone happens largely after the lower reservoir has deflated by 5 km³ (Figure 4D), which is discussed in Section 5.1.

### 4.3. The influence of differential reservoir pressure on shear zone geometries

We now explore the effects of applying a differential pressure at the walls of both reservoirs. We run numerical simulations varying $\Delta P_{up}/\Delta P_{low}$ from -1.0 to 1.0 (for example, for $\Delta P_{up}/\Delta P_{low}$ = -1.0, when the lower reservoir deflates by 100 MPa, the upper reservoir inflates by 100 MPa). $\Delta P_{low}$ is negative in these simulations. Figure 5 shows the effect of $\Delta P_{up}/\Delta P_{low}$ on the critical $\Delta P_{low}$ needed to reach the different stages described in the reference model case above.

For values of $\Delta P_{up}/\Delta P_{low}$ ranging from -1 to ≈ 0.3, the critical pressure for stage 4 increases with increasing $\Delta P_{up}/\Delta P_{low}$. In other words, rapid volume variation
<figure><img src="image_6.png" /></figure>

Figure 5: Effect of magma reservoir pressurization, for two spherical reservoirs. A: critical pressures for the different stages on a set of runs with varying $\Delta P_{up}/\Delta P_{low}$. Solid blue curve: applied critical magma pressure at which plastic deformation initiates at the lower reservoir wall. Dashed blue curve: same for the upper reservoir. Dashed black curve: critical applied magma pressure for which the lower reservoir has deflated by 5 km³. Red curve: onset of stage 4 (Figure 4). Black arrow: case described in Section 4.2 ($\Delta P_{up}/\Delta P_{low} = 0$); its direction indicates the progression through the run. The limit between “axial” and “conical” fracture zones is based on the observation of the final stress fields shown in Figure S4. B-D: final stages of the end-member simulations (see Figure S4 for intermediate cases). Principal stress orientations are shown similarly to Figure 4.
(accompanying brittle failure at the reservoir walls) begins earlier when the two reservoirs'
differential pressures have opposite signs. This effect reverses for $\Delta P_{up}/\Delta P_{low}$ greater than 0.3
(Fig. 5A, red curve).

The $\Delta P_{up}/\Delta P_{low}$ ratio also influences the geometry of the plastic domain (Figure S4). For negative ratios, a well-localized conical shear zone appears, as seen in the reference case (Figure 4). As this ratio increases, this conical shear zone becomes less and less localized, and strain tends to concentrate around the reservoirs’ walls. For $\Delta P_{up}/\Delta P_{low}$ above 0.5, the shear failure area shifts towards the center of the domain, which we refer to as an “axial fracture zone”.

In summary, applying opposite wall pressures on two reservoirs (1) tends to generate better-defined conical failure zones and (2) requires smaller applied pressure to establish a plastic domain connecting the reservoirs.

### 4.4. The effect of reservoir depths

Next, we investigate the effect of reservoir depths on the geometry of the plastic shear zone. Figure 6 shows simulations with the lower reservoir set at a central depth of 40 km and the upper reservoir set at varying depths (Table 1). We define the reservoirs’ separation height as the distance between the upper reservoir’s base and the lower reservoir’s apex (1). We examine three pressure regimes: $\Delta P_{up}/\Delta P_{low}$ equal to 0.3 (top row), 0 (middle row) and -0.3 (bottom row).

With $\Delta P_{up}/\Delta P_{low} = -0.3$, a conical shear zone appears for all tested upper reservoir depths. When the reservoirs are farther apart, the shear zone appears to be less well defined. The case with a separation of 13 km (Figure 6C) shows poorly localized strain, with $\epsilon^{pl}_{II}$ seemingly more concentrated on the reservoirs’ walls. Figure 6F (separation = 10.5 km) shows a ≈5-km wide network of curved faults. Cases with a smaller separation display more sharply defined conical fracture zones.

For the other two pressure regimes, the plastic domain connecting both reservoirs is hardly visible beyond a certain reservoir separation. For $\Delta P_{up}/\Delta P_{low} = 0$, it is poorly developed at 10.5 km separation (Figure 6E) and disappears at 13 km separation (Figure 6B). For $\Delta P_{up}/\Delta P_{low} = 0.3$, it is poorly defined at 8 km separation (Figure 6G) and disappears beyond 10.5 km separation (Figure 6A, D).
<figure><img src="image_7.png" /></figure>

Figure 6: Effect of the separation distance between reservoirs on the plastic strain field, for spherical reservoirs. For each case, the presented step is the last before solver divergence. $\Delta P_{up}/\Delta P_{low}$ ratios are 0.3 for the top row, 0 for the middle row and -0.3 for the bottom row. Horizontal and vertical axes are the radial distance and depth, in km. Reservoir separations are indicated by the black arrows. Principal stress orientations are shown similarly to Figure 4.

In summary, there seems to be a critical separation height below which a localized plastic connection can develop between the reservoirs. We find that opposite reservoir wall pressure signs raise this critical height.

### 4.5. The effect of reservoir geometries

Figures S5 and 7 show models featuring oblate and sill-like magma reservoirs. Reservoir depths are slightly modified from Table 1 in order to keep the same reservoir separation heights as in Figure 6.
With non-spherical reservoirs, the critical separation below which a well-defined plastic connection forms decreases. For example, for opposite pressure signs and a 13-km separation, the plastic connection is hardly visible with oblate reservoirs (Figure S5C) and it is absent with sill-like reservoirs (Figure 7C), whereas it is clearly visible in the spherical case (Figure 6C).

Grosfils et al. (2015) predicted that, at great depths, inflating oblate and sill-like reservoirs tend to generate laterally propagating fractures. With deflating reservoirs, the most compressive principal stress $σ_1$ shifts by 90° from wall-normal to wall-parallel, as can be seen on our results. But slip line plasticity theory (Gerbault, 2012, and references therein) predicts that shear bands form at an angle $45° ± φ/2$ from $σ_1$. Thus, fault orientation will not change significantly with the wall pressure sign (since we choose $φ = 1°$), which is why we observe laterally propagating faults. With oblate or sill-like reservoirs, these generated faults need to cover a broader curved pattern prior to flipping back and connecting with the upper reservoir, which could explain why it seems more difficult for pressurized sills to generate a connecting shear zone.

When a plastic connection does form in between the reservoirs, it seems to develop a better-defined cone, i.e. strain focuses into a thinner shear band, in contrast to a diffuse plastic domain or to strain concentration along the reservoir walls. For example, the case with two deflating reservoirs and a 5.5-km separation, which produced a rather diffuse conical failure zone in between the reservoirs and significant concentric strain along the lower one (Figure 6J), shifts to a well-defined conical fault zone when the reservoirs are oblate or sill-like (Figure S5J and Figure 7J). Stress and strain concentrations in areas of greater curvature (e.g., Gudmundsson, 2007; Gerbault et al., 2018, Figure S2) may explain why non-spherical reservoirs
<figure><img src="image_8.png" /><figcaption>Figure 7: Effect of the distance between reservoirs on the plastic strain field, for sill-like reservoirs ($h/w = 1/4$). Horizontal and vertical axes are the horizontal distance from the domain center and the depth, in km. For every case, the presented step is the last before solver convergence fails. $\Delta P_{up}/\Delta P_{low}$ ratios are 0.3 for the top row, 0 for the middle row and -0.3 for the bottom row. $d_{up}$ and $d_{low}$ differ from Figure 6 to make sure that the tested reservoir separations remain the same (see Table 1 caption), i.e. identical to Figure 6. Principal stress orientations are shown similarly to Figure 4. The insets on the right column are zooms in between 30 and 40 km depths.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

Figure 8: Interpretative sketch of the proximal earthquake cluster near Mayotte. Gray dots: VT earthquakes from Lavayssière et al. (2022). Background color: Vp/Vs ratio from Foix et al. (2021). The tomography data is cut below 45 km depth due to the lack of resolution. Red arrows: magma transfer. Blue arrows: possible fluid transfer (hydrothermal fluids, gases). Black arrows: magma reservoir pressurization. LP event location from Lavayssière and Retailleau (2023).

produce well-defined conical strained areas in our simulations. These well-defined shear strain areas appear for smaller applied pressures than for spherical reservoirs, which agrees with Currenti and Williams (2014) and Gerbault et al. (2018).

In summary, oblate and sill-like reservoirs (1) require shorter reservoir separation to generate a well-defined plastic connection, but (2) seem more prone to generate well-defined conical shear bands, and (3) require smaller applied pressures and volume variation to achieve this plastic connection.
**5. Discussion**

**5.1. Understanding the distribution of the earthquakes in the proximal cluster**

Although ground displacement inversions place the deflation source of the 2018-2021 eruption east of the proximal cluster, (Lemoine et al., 2020; Feuillet et al., 2021; Peltier et al., 2022), the area may still have undergone magma withdrawal and fed the eruption (Lemoine et al., 2020; Mittal et al., 2022). We thus used the estimated 5 km³ volume loss as an upper bound. For the same reason, we favor cases requiring applied pressures that are as moderate as possible.

Our goal was to produce plastic shear strain patterns that would match a given earthquake distribution, hence we assumed that brittle plasticity is a proxy for seismic activity, e.g. like in Gerbault et al. (2022). However, while plastic strain indicates areas undergoing brittle behavior (Vermeer and de Borst, 1984), it does not necessarily reflect a seismic distribution (Dieterich, 1979). Fault-induced stress changes can cause earthquakes slightly off the major faults (King and Cocco, 2001), making the earthquake distribution more diffuse than the shear strain pattern. Interpreting our strain distributions in terms of seismicity thus requires caution.

In our models, we find several cases where the interaction of two nearby reservoirs can produce a conical plastic domain matching the shape of the observed proximal earthquake cluster. In particular, scenarios with reservoir wall pressures of opposite signs are more prone to generate such shapes (Figures 5, S4), and require the lowest applied magma pressures (Figure 5A). Assuming spherical reservoirs, these scenarios are also the only ones where a well-defined fracture zone (stage 4) is established before the lower magma reservoir deflates by 5 km³.

Alternatively, oblate and sill-like reservoirs can generate conical plastic domains for smaller volume variations, even when $\Delta P_{up}/\Delta P_{low} \ge 0$ (Figures S5 and 7). However, these models raise two issues. First, oblate reservoirs or sills poorly match the geometry of the anomalies imaged by Foix et al. (2021) (Figure 1). This issue is minor, since recent works (e.g., Cashman et al., 2017, and references therein) propose that magma “reservoirs”, classically conceived as spheroidal cavities, actually consist of multi-layered sill complexes. Due to its low resolution, the tomography by Foix et al. (2021) cannot
discriminate between such structures and vaguely spherical melt-rich zones. Second, oblate spheroidal reservoirs must be closer to one another in order to establish a mechanical connection, especially for cases $\Delta P_{up}/\Delta P_{low} \ge 0$, and require reservoir height separations of less than 5.5-8 km (Figures S5 and 7). The proximal cluster extends over a depth interval of $\approx 15$ km (Figure 1, Lavayssière et al., 2022), hence reservoirs less than 8 km apart are too close to represent the Fani Maoré plumbing system. Only cases with $\Delta P_{up}/\Delta P_{low} < 0$ produce plastic connections for height separations of more than 10 km (e.g., Figure S5F for oblate reservoirs). Non-spherical reservoir geometries may contribute to generating the observed proximal earthquake distribution, but opposite reservoir pressure signs are our favored explanation.

Since we assume a deflating lower reservoir, we propose that an inflating reservoir may lie atop the proximal cluster (Figure 8). It would be located beneath a possible former caldera on the seafloor, where geologically recent volcanic features have been identified (Puzenat et al., 2022; REVOSIMA, 2023). There is currently no evidence for lava emission during the 2018-2021 period in this supposed caldera, despite active degassing being recorded in the Fer a` Cheval, on its southern rim, only 10 km east of Mayotte (Feuillet et al., 2021; REVOSIMA, 2023). A recharging reservoir below this location highlights the possibility of a future eruption closer to Mayotte than the Fani Maoré edifice. However, our study cannot predict the likelihood of such a scenario, let alone propose a related timescale.

## 5.2. Pressure conditions necessary for failure at the considered depths

The proximal cluster's depth (25-40 km) is a challenging point. The corresponding confining pressures, of the order of a GPa, make both tensile and shear failure especially hard to attain (Grosfils, 2007; Gerbault, 2012). Here, we produce shear failure by lowering the friction angle $\phi$ to an extremely low value (1°) with respect to usual standards in rock mechanics ($\approx 30°$). The case without friction is an end-member scenario tested in several modeling studies (e.g., Novoa et al., 2022, Gerbault et al., 2022), and mimics the influence of elevated bedrock pore fluid pressures (Rubey and King Hubbert, 1959).

In the case of Mayotte, our quasi-null friction hypothesis indicates that the mantle is likely weakened by preexisting sets of fractures. Lithospheric-scale fault zones likely
affect the mantle, given Mayotte's location on a plate boundary. (Famin et al., 2020; Thinon et al., 2022). Besides, the proximal cluster lies beneath volcanic structures on the seafloor (Feuillet et al., 2021). Whether these structures are directly connected to the proximal cluster is unclear, but their presence suggests that pre-existing magma-induced fractures may have been reactivated by the recent eruption. Emission of magma-originated gases on the seafloor near and on Mayotte (Feuillet et al., 2021; Liuzzo et al., 2022) also hint at a hydrothermal system reaching or approaching the identified magma reservoir depths; thus, serpentinization induced by hydrothermal fluids could contribute to mantle weakening (Brantut et al., 2016).

While the critical pressures for the onset of reservoir failure do not exceed 50 MPa in our simulations, the pressures required to develop a well-defined conical shear zone range from ≈ 125 MPa to ≈ 250 MPa for spherical reservoirs (Figure 5). For sill-like reservoirs, they are lower but still on the order of 100 MPa (Figure S7). Numerical modeling studies of magma reservoirs generally feature applied pressures below a few dozen MPa (e.g., Gudmundsson, 2007; Currenti and Williams, 2014; Browning et al., 2021; Novoa et al., 2022; Karaoğlu et al., 2020). Some studies apply pressures beyond 100 MPa for theoretical purposes, and in models featuring both gravity and standard friction values (Gerbault et al., 2012, 2018). Soltanmohammadi et al. (2021) model magma ascent atop mantle plumes via compaction waves, and compute pressures up to 150 MPa. Kiss et al. (2023) also present models with over/underpressures nearing 400 MPa due to thermal contraction/dilatation. Our simulations thus fall on the high end of what is usually tested in terms of magma pressures.

Here, we model deeper magma reservoirs than the studies cited above, and although we examine essentially friction-less cases, our simulations do include gravity and lithostatic loading. While reservoir pressures beyond 100 MPa are problematic for shallower or gravity-less cases (due to tensile failure occurring prior to shear failure for example), it is hard to assess how (un)realistic such applied pressures are in our framework. We note, however, that pore fluids are likely present near Mayotte, in the form of hydrothermal fluids (Feuillet et al., 2021; Liuzzo et al., 2022) or of a magma mush (Foix et al., 2021; Mittal et al., 2022). Partial reduction of the confining pressure by pore fluids (Grosfils et al., 2015) is likely at play near Mayotte. In these conditions, it seems plausible that an overpressure of ≈ 100 MPa would suffice to trigger tensile failure at 20-40 km depth.
An alternative to high magma pressures triggering bedrock seismicity would be that the fractures or faults generating the proximal swarm predate the 2018-2021 eruption. Slipping along a preexisting fault plane requires less stress than along new faults (Lockner and Beeler, 2002; Collettini et al., 2009). In our modeled framework, reactivation of a preexisting fault zone connecting the reservoirs would only require the pressure difference between stages 3 and 4 (see Section 4.2), which is several tens of MPa less than that required to go from stage 1 to stage 4. We lack the data needed to assess the state of the tectono-magmatic system before May 2018, but Mittal et al. (2022) proposed that the absence of delay between magma emission and seismicity in early 2019 could indicate that the proximal cluster faults were already close to critical state and therefore reacted rapidly to magma withdrawal. Further studies are required to better characterize regional tectonic inheritance.

Other factors than a reservoir over/underpressure may also contribute to host rock failure. In choosing a uniform domain density, we neglected magma buoyancy effects (a more realistic magma density would be ≈ 2500 kg.m⁻³). Sigmundsson et al. (2020) showed that a buoyant magma reservoir exerts a stress field on its roof that can bring it close to failure without extra magma over/underpressure. An 8-km-high magma column like our lower reservoir, with density 2500 kg.m⁻³ embedded in the mantle (Δρ = 800 kg.m⁻³) will exert an extra pressure Δρgh ≈ 63 MPa on the reservoir roof. Whether this extra pressure would cause host rock failure is difficult to say. Gudmundsson (2007) argues that reservoir pressures overcoming the host rock's tensile strength (generally 10 MPa) will trigger tensile failure at any depth, but their approach assumes that the host rock is saturated by pore fluids at lithostatic pressure, allowing them to neglect the confining pressure (Grosfils et al., 2015). Although pore fluids seem present near Mayotte as mentioned above, total cancellation of the lithostatic pressure over the entire 40-km thick column seems to be a strong hypothesis.

### 5.3. *Perspectives for further studies*

Our modeling could benefit from a more complete consideration of physical effects. We treated our problem via simple solid mechanics, albeit with an elaborate rheology to model fault initiation and propagation. First, we do not account for the high temperature at 25-40 km depth, especially near magma reservoirs. This thermal state could trigger
viscous behavior instead of brittle behavior (Currenti and Williams, 2014; Zhan and Gregg, 2019; Kiss et al., 2023). Viscous flow (Head et al., 2022) and thermal expansion (Karaoğlu et al., 2020; Kiss et al., 2023) also largely impact stress magnitude and distribution.

Future models could also include fluid transfer between the reservoirs and the host rock. The location of LP seismicity inside the VT gap of the proximal cluster (Figure 8, Lavayssière and Retailleau, 2023) suggests a fluid pathway between the two reservoirs, which could be modeled as a conduit in the intervening domain (Reverso et al., 2014), as a poro-elasto-plastic domain within which fluids can be pumped in and out transiently (Sáez-Leiva et al., 2023), or as a viscous fluid-filled propagating crack (Furst et al., 2023). These approaches could help reproduce more consistent volume variations in between the reservoirs with lower pressures than those imposed here.

Our choice of axially symmetrical 2D models is based on the first-order shape of the proximal cluster. Actually, the proximal earthquake distribution is not perfectly symmetrical, nor does it seem to evolve symmetrically with time (Lavayssière and Retailleau, 2023). Besides, onshore and bathymetric volcanic and tectonic features clearly follow a NW-SE tensile regional stress pattern, and Mayotte is located on a transtensional oceanic-to-continental plate boundary (Famin et al., 2020; Feuillet et al., 2021; Thinon et al., 2022). A three-dimensional setup then becomes necessary to account for the influence of such an anisotropic regional stress field.

Finally, a complete and quantitative study of magma transfer during the Fani Maoré eruption should also include the distal cluster area. Concurring changes in seismicity in mid-2020 suggest a connection between the proximal and distal seismic clusters (Lavayssière and Retailleau, 2023). It has been convincingly argued that the erupted magma transited through the area of the distal cluster located further east (Cesca et al., 2020; Lemoine et al., 2020; Mercury et al., 2023). This lateral propagation of magma calls for further mechanical explanation.

## 6. Conclusion

Our work shows that bedrock failure in-between two superimposed magma reservoirs (one lying around 40 km depth, the other around 20-25 km depth, as indicated by geophysical and petrological data) is a plausible explanation for the presence of the proximal earthquake
cluster offshore east Mayotte island. The proximal cluster's earthquake distribution seems to favor opposite pressure signs on the reservoir walls (one inflating, one deflating). Assuming, based on data, that the deeper reservoir has undergone magma withdrawal during the recent eruption suggests that the shallower reservoir near Moho depth may have been recharged during the same period. Inflation of this upper magma reservoir could allow upward magma propagation to the surface in the future; further studies are needed to ascertain whether this phenomenon could cause an eruption closer to the island. In the meantime, such a structure should be closely monitored due to the hazard it represents.

## Acknowledgments

Since June 2019, Mayotte seismo-volcanic activity is monitored by the REVOSIMA (Réseau de surveillance volcanologique et sismologique de Mayotte; DOI: 10.18715/MAYOTTE.REVO and funded by the Ministry for Ecological Transition and Territorial Cohesion (MTECT), the Ministry of Higher Education and Research (MESR), the Ministry of the Interior and Overseas (MIOM), with the support of the French Ministry for Armed Forces (MINARM). Marine operations benefit from the support of the French Oceanographic Fleet (FOF). Seismic data shown in this paper have been acquired and processed in the framework of the REVOSIMA. We thank Océane Foix for providing us with tomography data. This work is supported by a PhD funding from Université Paris-Cité. Numerical modeling was carried out with COMSOL Multiphysics® v6.1 (https://uk.comsol.com), through the MPH Python interface (Hennig et al., 2023). Figure 1 was made using Generic Mapping Tools via the PyGMT Python module (Jones et al., 2022). We are thankful to editor Diana Roman for handling the review process, and two anonymous reviewers for their helpful feedback.

## Declaration of competing interests

The authors have no competing interests to declare.

## References

* Afonso, J.C., Ranalli, G., Fernández, M., 2005. Thermal expansivity and elastic properties of the lithospheric mantle: results from mineral physics of composites. Physics of the
Earth and Planetary Interiors 149, 279-306. URL: https://linkinghub.elsevier.com/retrieve/pii/S0031920104003747, doi:10.1016/j.pepi.2004.10.003.

Agútsdóttir, T., Winder, T., Woods, J., White, R., Greenfield, T., Brandsdóttir, B., 2019. Intense Seismicity During the 2014-2015 Bárðarbunga-Holuhraun Rifting Event, Iceland, Reveals the Nature of Dike-Induced Earthquakes and Caldera Collapse Mechanisms. Journal of Geophysical Research: Solid Earth 124, 8331-8357. doi:https://doi.org/10.1029/2018JB016010.

Albino, F., Sigmundsson, F., 2014. Stress transfer between magma bodies: Influence of intrusions prior to 2010 eruptions at Eyjafjallajökull volcano, Iceland. Journal of Geophysical Research: Solid Earth 119, 2964-2975. URL: https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2013JB010510, doi:10.1002/2013JB010510.

Alshembari, R., Hickey, J., Williamson, B.J., Cashman, K., 2022. Poroelastic Mechanical Behavior of Crystal Mush Reservoirs: Insights Into the Spatio-Temporal Evolution of Volcano Surface Deformation. Journal of Geophysical Research: Solid Earth 127, e2022JB024332. URL: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JB024332, doi:10.1029/2022JB024332.

Alshembari, R., Hickey, J., Williamson, B.J., Cashman, K., 2023. Unveiling the Rheological Control of Magmatic Systems on Volcano Deformation: The Interplay of Poroviscoelastic Magma-Mush and Thermo-Viscoelastic Crust. Journal of Geophysical Research: Solid Earth 128, e2023JB026625. URL: https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023JB026625, doi:10.1029/2023JB026625.

Berthod, C., Komorowski, J.C., Gurioli, L., Médard, E., Bachèlery, P., Besson, P., Verdurme, P., Chevrel, O., Di Muro, A., Peltier, A., Devidal, J.L., Nowak, S., Thinon, I., Burckel, P., Hidalgo, S., Deplus, C., Loubrieu, B., Pierre, D., Bermell, S., Pitel-Roudaut, M., Réaud, Y., Fouchard, S., Bickert, M., Le Friant, A., Paquet, F., Feuillet, N., Jorry, S.L., Fouquet, Y., Rinnert, E., Cathalot, C., Lebas, E., 2022. Temporal magmatic evolution of the Fani Maoré submarine eruption 50 km east of Mayotte
revealed by in situ sampling and petrological monitoring. Comptes Rendus. Géoscience
354, 1-29. URL: https://comptes-rendus.academiesciences.fr/geoscience/articles/10.5802/crgeos.155/, doi:10.5802/crgeos.155.

Berthod, C., Médard, E., Bachèlery, P., Gurioli, L., Di Muro, A., Peltier, A., Komorowski, J.C., Benbakkar, M., Devidal, J.L., Langlade, J., Besson, P., Boudon, G., Rose-Koga, E., Deplus, C., Le Friant, A., Bickert, M., Nowak, S., Thinon, I., Burckel, P., Hidalgo, S., Kaliwoda, M., Jorry, S.J., Fouquet, Y., Feuillet, N., 2021a. The 2018-ongoing Mayotte submarine eruption: Magma migration imaged by petrological monitoring. Earth and Planetary Science Letters 571, 117085. URL: https://linkinghub.elsevier.com/retrieve/pii/S0012821X2100340X, doi:10.1016/j.epsl.2021.117085.

Berthod, C., Médard, E., Di Muro, A., Hassen Ali, T., Gurioli, L., Chauvel, C., Komorowski, J.C., Bachèlery, P., Peltier, A., Benbakkar, M., Devidal, J.L., Besson, P., Le Friant, A., Deplus, C., Nowak, S., Thinon, I., Burckel, P., Hidalgo, S., Feuillet, N., Jorry, S., Fouquet, Y., 2021b. Mantle xenolith-bearing phonolites and basanites feed the active volcanic ridge of Mayotte (Comoros archipelago, SW Indian Ocean). Contributions to Mineralogy and Petrology, 24.

Brantut, N., Passelègue, F.X., Deldicque, D., Rouzaud, J.N., Schubnel, A., 2016. Dynamic weakening and amorphization in serpentinite during laboratory earthquakes. Geology 44, 607-610. URL: https://pubs.geoscienceworld.org/geology/article/44/8/607-610/188173, doi:10.1130/G37932.1.

Browning, J., Karaoğlu, Ö., Bayer, Ö., Turgay, M.B., Acocella, V., 2021. Stress fields around magma chambers influenced by elastic thermo-mechanical deformation: implications for forecasting chamber failure. Bulletin of Volcanology 83, 48. URL: https://link.springer.com/10.1007/s00445-021-01471-2, doi:10.1007/s00445-021-01471-2.

Cabaniss, H.E., Gregg, P.M., Nooner, S.L., Chadwick, W.W., 2020. Triggering of eruptions at Axial Seamount, Juan de Fuca Ridge. Scientific Reports 10, 10219. URL: http://www.nature.com/articles/s41598-020-67043-0, doi:10.1038/s41598-020-67043-0.
Cashman, K.V., Sparks, R.S.J., Blundy, J.D., 2017. Vertically extensive and unstable magmatic systems: A unified view of igneous processes. Science 355, eaag3055. URL: https://www.science.org/doi/10.1126/science.aag3055, doi:10.1126/science.aag3055.

Cesca, S., Letort, J., Razafindrakoto, H.N.T., Heimann, S., Rivalta, E., Isken, M.P., Nikkhoo, M., Passarelli, L., Petersen, G.M., Cotton, F., Dahm, T., 2020. Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation. Nature Geoscience 13, 87–93. URL: http://www.nature.com/articles/s41561-019-0505-5, doi:10.1038/s41561-019-0505-5.

Chouet, B.A., 1996. Long-period volcano seismicity: its source and use in eruption forecasting. Nature 380, 309–316. URL: http://www.nature.com/articles/380309a0, doi:10.1038/380309a0.

Collettini, C., Niemeijer, A., Viti, C., Marone, C., 2009. Fault zone fabric and fault weakness. Nature 462, 907–910. URL: https://www.nature.com/articles/nature08585, doi:10.1038/nature08585.

Currenti, G., Williams, C.A., 2014. Numerical modeling of deformation and stress fields around a magma chamber: Constraints on failure conditions and rheology. Physics of the Earth and Planetary Interiors 226, 14–27. URL: https://linkinghub.elsevier.com/retrieve/pii/S0031920113001623, doi:10.1016/j.pepi.2013.11.003.

Dieterich, J.H., 1979. Modeling of rock friction: 1. experimental results and constitutive equations. Journal of Geophysical Research: Solid Earth 84, 2161–2168. URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB084iB05p02161, doi:https://doi.org/10.1029/JB084iB05p02161, arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/10.1029/JB084iB05p02161.

Famin, V., Michon, L., Bourhane, A., 2020. The Comoros archipelago: a right-lateral transform boundary between the Somalia and Lwandle plates. Tectonophysics 789, 228539. URL: https://linkinghub.elsevier.com/retrieve/pii/S0040195120302225, doi:10.1016/j.tecto.2020.228539.
Feuillet, N., Jorry, S., Crawford, W.C., Deplus, C., Thinon, I., Jacques, E., Saurel, J.M., Lemoine, A., Paquet, F., Satriano, C., Aiken, C., Foix, O., Kowalski, P., Laurent, A., Rinnert, E., Cathalot, C., Donval, J.P., Guyader, V., Gaillot, A., Scalabrin, C., Moreira, M., Peltier, A., Beauducel, F., Grandin, R., Ballu, V., Daniel, R., Pelleau, P., Gomez, J., Besançon, S., Geli, L., Bernard, P., Bachelery, P., Fouquet, Y., Bertil, D., Lemarchand, A., Van der Woerd, J., 2021. Birth of a large volcanic edifice offshore Mayotte via lithosphere-scale dyke intrusion. Nature Geoscience 14, 787–795. URL: https://www.nature.com/articles/s41561-021-00809-x, doi:10.1038/s41561-021-00809-x.

Foix, O., Aiken, C., Saurel, J.M., Feuillet, N., Jorry, S., Rinnert, E., Thinon, I., 2021. Offshore Mayotte volcanic plumbing revealed by local passive tomography. Journal of Volcanology and Geothermal Research, 14 URL: https://www.sciencedirect.com/science/article/pii/S0377027321002249, doi:https://doi.org/10.1016/j.jvolgeores.2021.107395.

Furst, S., Maccaferri, F., Pinel, V., 2023. Modeling the shape and velocity of magmatic intrusions, a new numerical approach. Journal of Geophysical Research: Solid Earth 128, e2022JB025697. URL: https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022JB025697, doi:https://doi.org/10.1029/2022JB025697.

Gerbault, M., 2012. Pressure conditions for shear and tensile failure around a circular magma chamber: insight from elasto-plastic modelling. Geological Society, London, Special Publications 367, 111–130. URL: http://sp.lyellcollection.org/lookup/doi/10.1144/SP367.8, doi:10.1144/SP367.8.

Gerbault, M., Cappa, F., Hassani, R., 2012. Elasto-plastic and hydromechanical models of failure around an infinitely long magma chamber: FAILURE AROUND A MAGMA CHAMBER. Geochemistry, Geophysics, Geosystems 13. URL: http://doi.wiley.com/10.1029/2011GC003917, doi:10.1029/2011GC003917.

Gerbault, M., Fontaine, F.J., Peltier, A., Got, J.L., Hassani, R., Ferrazzini, V., Gailler, L., Duputel, Z., 2022. What causes the persistent seismicity below the eastern flank of
Piton de la Fournaise (la Réunion Island)? Elasto-plastic models of magma inflation.
Journal of Volcanology and Geothermal Research 431, 107628. URL:
https://linkinghub.elsevier.com/retrieve/pii/S0377027322001597,
doi:10.1016/j.jvolgeores.2022.107628.

Gerbault, M., Hassani, R., Novoa Lizama, C., Souche, A., 2018. Three-Dimensional Failure Patterns Around an Inflating Magmatic Chamber. Geochemistry, Geophysics, Geosystems 19, 749–771. URL: http://doi.wiley.com/10.1002/2017GC007174, doi:10.1002/2017GC007174.

Gerbault, M., Poliakov, A.N., Daignieres, M., 1998. Prediction of faulting from the theories of elasticity and plasticity: what are the limits? Journal of Structural Geology 20, 301–320. URL: https://linkinghub.elsevier.com/retrieve/pii/S0191814197000898, doi:10.1016/S0191-8141(97)00089-8.

Gonnermann, H.M., Foster, J.H., Poland, M., Wolfe, C.J., Brooks, B.A., Miklius, A., 2012. Coupling at Mauna Loa and Kilauea by stress transfer in an asthenospheric melt layer. Nature Geoscience 5, 826–829. URL: https://www.nature.com/articles/ngeo1612, doi:10.1038/ngeo1612.

Got, J.L., Peltier, A., Staudacher, T., Kowalski, P., Boissier, P., 2013. Edifice strength and magma transfer modulation at Piton de la Fournaise volcano: EDIFICE STRENGTH AND MAGMA TRANSFER. Journal of Geophysical Research: Solid Earth 118, 5040–5057. URL: http://doi.wiley.com/10.1002/jgrb.50350, doi:10.1002/jgrb.50350.

Grosfils, E.B., 2007. Magma reservoir failure on the terrestrial planets: Assessing the importance of gravitational loading in simple elastic models. Journal of Volcanology and Geothermal Research 166, 47–75. URL: https://linkinghub.elsevier.com/retrieve/pii/S0377027307001941, doi:10.1016/j.jvolgeores.2007.06.007.

Grosfils, E.B., McGovern, P.J., Gregg, P.M., Galgana, G.A., Hurwitz, D.M., Long, S.M., Chestler, S.R., 2015. Elastic models of magma reservoir mechanics: a key tool for investigating planetary volcanism. Geological Society, London, Special Publications
401, 239-267. URL: http://sp.lyellcollection.org/lookup/doi/10.1144/SP401.2, doi:10.1144/SP401.2.

Gudmundsson, A., 2007. Conceptual and numerical models of ring-fault formation. Journal of Volcanology and Geothermal Research 164, 142-160. URL: https://linkinghub.elsevier.com/retrieve/pii/S0377027307001333, doi:10.1016/j.jvolgeores.2007.04.018.

Head, M., Hickey, J., Thompson, J., Gottsmann, J., Fournier, N., 2022. Rheological Controls on Magma Reservoir Failure in a Thermo-Viscoelastic Crust. Journal of Geophysical Research: Solid Earth 127. URL: https://onlinelibrary.wiley.com/doi/10.1029/2021JB023439, doi:10.1029/2021JB023439.

Hennig, J., Elfner, M., Maeder, A., Feder, J., 2023. MPh-py/MPh: 1.2.3. URL: https://zenodo.org/record/7749502, doi:10.5281/zenodo.7749502.

Holohan, E.P., Walter, T.R., Schöpfer, M.P., Walsh, J.J., Wyk de Vries, B., Troll, V.R., 2013. Origins of oblique-slip faulting during caldera subsidence. Journal of Geophysical Research: Solid Earth 118, 1778-1794. URL: https://onlinelibrary.wiley.com/doi/10.1002/jgrb.50057, doi:10.1002/jgrb.50057.

Jaeger, J.C., Cook, N.G.W., Zimmerman, R., 2007. Fundamentals of Rock Mechanics. 4 ed., Champan and Hall, London.

Jones, M., Uieda, L., Leong, W.J., 2022. Geospatial analysis and visualization with PyGMT (SciPy 2022). doi:https://doi.org/10.6084/m9.fishshare.20483793.v1.

Karaoğlu, Ö., Bayer, Ö., Turgay, M.B., Browning, J., 2020. Thermomechanical interactions between crustal magma chambers in complex tectonic environments: Insights from Eastern Turkey. Tectonophysics 793, 228607. URL: https://linkinghub.elsevier.com/retrieve/pii/S0040195120302900, doi:10.1016/j.tecto.2020.228607.

King, G., Cocco, M., 2001. Fault interaction by elastic stress changes: New clues from earthquake sequences, in: Advances in Geophysics. Elsevier. volume 44, pp. 1-VIII.
URL: https://linkinghub.elsevier.com/retrieve/pii/S0065268700800060,
doi:10.1016/S0065-2687(00)80006-0.

Kiss, D., Moulas, E., Kaus, B.J.P., Spang, A., 2023. Decompression and Fracturing
Caused by Magmatically Induced Thermal Stresses. Journal of Geophysical Research:
Solid Earth 128, e2022JB025341. URL:
https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2022JB025341,
doi:10.1029/2022JB025341.

Laurent, A., 2023. Détection, localisation et caractérisation des signaux très basses fréquences dans la crise sismique de Mayotte (2018-2022). Ph.D. thesis. Université Paris-Cité, Institut de Physique du Globe de Paris.

Lavayssière, A., Crawford, W.C., Saurel, J.M., Satriano, C., Feuillet, N., Jacques, E., Komorowski, J.C., 2022. A new 1D velocity model and absolute locations image the Mayotte seismo-volcanic region. Journal of Volcanology and Geothermal Research 421, 107440. URL: https://linkinghub.elsevier.com/retrieve/pii/S0377027321002699, doi:10.1016/j.jvolgeores.2021.107440.

Lavayssière, A., Retailleau, L., 2023. Capturing Mayotte's deep magmatic plumbing system and its spatiotemporal evolution with volcano-tectonic seismicity. Volcanica 6, 331-344. URL: https://www.jvolcanica.org/ojs/index.php/volcanica/article/view/219, doi:10.30909/vol.06.02.331344.

Lemoine, A., Briole, P., Bertil, D., Roullé, A., Foumelis, M., Thinon, I., Raucoules, D., de Michele, M., Valty, P., Hoste Colomer, R., 2020. The 2018-2019 seismo-volcanic crisis east of Mayotte, Comoros islands: seismicity and ground deformation markers of an exceptional submarine eruption. Geophysical Journal International 223, 22-44. URL: https://academic.oup.com/gji/article/223/1/22/5850758, doi:10.1093/gji/ggaa273.

Levy, S., Bohnenstiehl, D., Sprinkle, P., Boettcher, M., Wilcock, W., Tolstoy, M., Waldhauser, F., 2018. Mechanics of fault reactivation before, during, and after the 2015 eruption of Axial Seamount. Geology 46, 447-450. URL: https://pubs.geoscienceworld.org/gsa/geology/article/46/5/447/529571/Mechanics-of-fault-reactivation-before-during-and, doi:10.1130/G39978.1.
Liao, Y., Soule, S.A., Jones, M., Le Mével, H., 2021. The Mechanical Response of a Magma Chamber With Poroviscoelastic Crystal Mush. Journal of Geophysical Research: Solid Earth 126. URL: https://onlinelibrary.wiley.com/doi/10.1029/2020JB019395, doi:10.1029/2020JB019395.

Liuzzo, M., Di Muro, A., Rizzo, A.L., Grassa, F., Coltorti, M., Ader, M., 2022. The composition of gas emissions at Petite Terre (Mayotte, Comoros): inference on magmatic fingerprints. Comptes Rendus. Géoscience 354, 1-24. URL: https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.5802/crgeos.148/, doi:10.5802/crgeos.148.

Lockner, D.A., Beeler, N.M., 2002. 32 Rock failure and earthquakes, in: International Geophysics. Elsevier. volume 81, pp. 505-537. URL: https://linkinghub.elsevier.com/retrieve/pii/S0074614202802352, doi:10.1016/S0074-6142(02)80235-2.

Masquelet, C., Leroy, S., Delescluse, M., Chamot-Rooke, N., Thinon, I., Lemoine, A., Franke, D., Watremez, L., Werner, P., Paquet, F., Berthod, C., Cabiativa Pico, V., Sauter, D., 2022. The East-Mayotte new volcano in the Comoros Archipelago: structure and timing of magmatic phases inferred from seismic reflection data. Comptes Rendus. Géoscience 354, 1-15. URL: https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.5802/crgeos.154/, doi:10.5802/crgeos.154.

Mazzullo, A., Stutzmann, E., Montagner, J.P., Kiselev, S., Maurya, S., Barruol, G., Sigloch, K., 2017. Anisotropic Tomography Around La Réunion Island From Rayleigh Waves: TOMOGRAPHY OF WESTERN INDIAN OCEAN. Journal of Geophysical Research: Solid Earth 122, 9132-9148. URL: http://doi.wiley.com/10.1002/2017JB014354, doi:10.1002/2017JB014354.

McTigue, D.F., 1987. Elastic stress and deformation near a finite spherical magma body: Resolution of the point source paradox. Journal of Geophysical Research 92, 12931. URL: http://doi.wiley.com/10.1029/JB092iB12p12931, doi:10.1029/JB092iB12p12931.
Mercury, N., Lemoine, A., Doubre, C., Bertil, D., Van Der Woerd, J., Hoste-Colomer, R., Battaglia, J., 2023. Onset of a submarine eruption east of Mayotte, Comoros archipelago: the first ten months seismicity of the seismo-volcanic sequence (2018–2019). Comptes Rendus. Géoscience 354, 105–136. URL: https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.5802/crgeos.191/, doi:10.5802/crgeos.191.

Michon, L., 2016. The Volcanism of the Comoros Archipelago Integrated at a Regional Scale, in: Bachelery, P., Lenat, J.F., Di Muro, A., Michon, L. (Eds.), Active Volcanoes of the Southwest Indian Ocean. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 333–344. URL: http://link.springer.com/10.1007/978-3-642-31395-0_21, doi:10.1007/978-3-642-31395-0_21. series Title: Active Volcanoes of the World.

Mittal, T., Jordan, J.S., Retailleau, L., Beauducel, F., Peltier, A., 2022. Mayotte 2018 eruption likely sourced from a magmatic mush. Earth and Planetary Science Letters 590, 117566. URL: https://linkinghub.elsevier.com/retrieve/pii/S0012821X22002023, doi:10.1016/j.epsl.2022.117566.

Mittal, T., Richards, M.A., 2019. Volatile Degassing From Magma Chambers as a Control on Volcanic Eruptions. Journal of Geophysical Research: Solid Earth 124, 7869–7901. URL: https://onlinelibrary.wiley.com/doi/10.1029/2018JB016983, doi:10.1029/2018JB016983.

Novoa, C., Gerbault, M., Remy, D., Cembrano, J., Lara, L., Ruz-Ginouves, J., Tassara, A., Baez, J., Hassani, R., Bonvalot, S., Contreras-Arratia, R., 2022. The 2011 Cordón Caulle eruption triggered by slip on the Liquiñe-Ofqui fault system. Earth and Planetary Science Letters 583, 117386. URL: https://linkinghub.elsevier.com/retrieve/pii/S0012821X2200022X, doi:10.1016/j.epsl.2022.117386.

Pascal, K., Neuberg, J., Rivalta, E., 2014. On precisely modelling surface deformation due to interacting magma chambers and dykes. Geophysical Journal International 196, 253–278. URL: https://academic.oup.com/gji/article/196/1/253/585883, doi:10.1093/gji/ggt343.
Peltier, A., Saur, S., Ballu, V., Beauducel, F., Briole, P., Chanard, K., Dausse, D., De Chabalier, J.B., Grandin, R., Rouffiac, P., Tranchant, Y.T., de Berc, M.B., Besançon, S., Boissier, P., Broucke, C., Brunet, C., Canjamalé, K., Carme, E., Catherine, P., Colombain, A., Crawford, W., Daniel, R., Dectot, G., Desfete, N., Doubre, C., Dumouch, T., Griot, C., Grunberg, M., Jund, H., Kowalski, P., Lauret, F., Lebreton, J., Pesqueira, F., Tronel, F., Valty, P., van der Woerd, J., 2022. Ground deformation monitoring of the eruption offshore Mayotte. Comptes Rendus. Géoscience 354, 1-23.
URL: https://comptes-rendus.academiesciences.fr/geoscience/articles/10.5802/crgeos.176/, doi:10.5802/crgeos.176.

Puzenat, V., Feuillet, N., Komorowski, J.C., Escartín, J., Deplus, C., Bachlery, P., Berthod, C., Gurioli, L., Scalabrin, C., Cathalot, C., Rinnert, E., Loubrieu, B., Pierre, D., Pitel-Roudaut, M., Tanguy, N., Fouquet, Y., Jorry, S.J., Lebas, E., Paquet, F., Thinon, I., 2022. Volcano-tectonic structures of Mayotte's upper submarine slope: insights from high-resolution bathymetry and in-situ imagery from a deep-towed camera. Comptes Rendus. Géoscience 354, 1-24. URL: https://comptes-rendus.academiesciences.fr/geoscience/articles/10.5802/crgeos.175/, doi:10.5802/crgeos.175.

Retailleau, L., Saurel, J.M., Laporte, M., Lavayssière, A., Ferrazzini, V., Zhu, W., Beroza, G.C., Satriano, C., Komorowski, J.C., OVPF Team, 2022. Automatic detection for a comprehensive view of Mayotte seismicity. Comptes Rendus. Géoscience 354, 1-18.
URL: https://comptes-rendus.academiesciences.fr/geoscience/articles/10.5802/crgeos.133/, doi:10.5802/crgeos.133.

Reverso, T., Vandemeulebrouck, J., Jouanne, F., Pinel, V., Villemin, T., Sturkell, E., Bascou, P., 2014. A two-magma chamber model as a source of deformation at Grímsvötn volcano, Iceland. Journal of Geophysical Research: Solid Earth 119, 4666-4683.

REVOSIMA, 2023. Bulletin de l'activité sismo-volcanique à Mayotte du 1 au 30 avril 2023.
URL: https://www.ipgp.fr/wp-content/uploads/2023/05/Revosima_bull_20230505.pdf.

Rinnert, E., Feuillet, N., Fouquet, Y., Jorry, S., Thinon, I., Lebas, E., 2019. MAYOBS doi:10.18142/291.
Rubey, W.W., King Hubbert, M., 1959. ROLE OF FLUID PRESSURE IN MECHANICS OF OVERTHRUST FAULTING: II. OVERTHRUST BELT IN GEOSYNCLINAL AREA OF WESTERN WYOMING IN LIGHT OF FLUID-PRESSURE HYPOTHESIS. GSA Bulletin 70, 167-206. URL: https://doi.org/10.1130/0016-7606(1959)70[167:ROFPIM]2.0.CO;2, doi:10.1130/0016-7606(1959)70[167:ROFPIM]2.0.CO;2.

Saurel, J.M., Jacques, E., Aiken, C., Lemoine, A., Retailleau, L., Lavayssière, A., Foix, O., Dofal, A., Laurent, A., Mercury, N., Crawford, W., Lemarchand, A., Daniel, R., Pelleau, P., Bès de Berc, M., Dectot, G., Bertil, D., Roullé, A., Broucke, C., Colombain, A., Jund, H., Besanc, on, S., Guyavarch, P., Kowalski, P., Roudaut, M., Apprioual, R., Battaglia, J., Bodihar, S., Boissier, P., Bouin, M.P., Brunet, C., Canjamale, K., Catherine, P., Desfete, N., Doubre, C., Dretzen, R., Dumouche, T., Fernagu, P., Ferrazzini, V., Fontaine, F.R., Gaillot, A., Géli, L., Griot, C., Grunberg, M., Guzel, E.C., Hoste-Colomer, R., Lambotte, S., Lauret, F., Léger, F., Maros, E., Peltier, A., Vergne, J., Satriano, C., Tronel, F., Van der Woerd, J., Fouquet, Y., Jorry, S.J., Rinnert, E., Thinon, I., Feuillet, N., 2021. Mayotte seismic crisis: building knowledge in near real-time by combining land and ocean-bottom seismometers, first results. Geophysical Journal International 228, 1281-1293. URL: https://academic.oup.com/gji/article/228/2/1281/6374867, doi:10.1093/gji/ggab392.

Segall, P., 2010. Earthquake and Volcano Deformation. Princeton University Press. URL: https://www.degruyter.com/document/doi/10.1515/9781400833856/html?lang=en, doi:10.1515/9781400833856.

Sigmundsson, F., Pinel, V., Grapenthin, R., Hooper, A., Halldórsson, S.A., Einarsson, P., Ofeigsson, B.G., Heimisson, E.R., Jónsdóttir, K., Gudmundsson, M.T., Vogfjörd, K., Parks, M., Li, S., Drouin, V., Geirsson, H., Dumont, S., Fridriksdottir, H.M., Gudmundsson, G.B., Wright, T.J., Yamasaki, T., 2020. Unexpected large eruptions from buoyant magma bodies within viscoelastic crust. Nature Communications 11, 2403. URL: http://www.nature.com/articles/s41467-020-16054-6, doi:10.1038/s41467-020-16054-6.
Soltanmohammadi, A., Grégoire, M., Fontaine, F.J., Bédard, L.P., Blanchard, M., Rabinowicz, M., 2021. Melt percolation, concentration and dyking in the hawaiian mantle plume and overriding lithosphere: links to the evolution of lava composition along the volcanic chain. Journal of Petrology 63, 1–24.

Sáez-Leiva, F., Hurtado, D.E., Gerbault, M., Ruz-Ginouves, J., Iturrieta, P., Cembrano, J., 2023. Fluid flow migration, rock stress and deformation due to a crustal fault slip in a geothermal system: A poro-elasto-plastic perspective. Earth and Planetary Science Letters 604, 117994. URL: https://linkinghub.elsevier.com/retrieve/pii/S0012821X23000079, doi:10.1016/j.epsl.2023.117994.

Thinon, I., Lemoine, A., Leroy, S., Paquet, F., Berthod, C., Zaragosi, S., Famin, V., Feuillet, N., Boymond, P., Masquelet, C., Mercury, N., Rusquet, A., Scalabrin, C., Van der Woerd, J., Bernard, J., Bignon, J., Clouard, V., Doubre, C., Jacques, E., Jorry, S.J., Rolandone, F., Chamot-Rooke, N., Delescluse, M., Franke, D., Watremez, L., Bachèlery, P., Michon, L., Sauter, D., Bujan, S., Canva, A., Dassie, E., Roche, V., Ali, S., Sitti Allaouia, A.H., Deplus, C., Rad, S., Sadeski, L., 2022. Volcanism and tectonics unveiled in the Comoros Archipelago between Africa and Madagascar. Comptes Rendus. Géoscience 354, 1–28. URL: https://comptes-rendus.academie-sciences.fr/geoscience/articles/10.5802/crgeos.159/, doi:10.5802/crgeos.159.

Vermeer, P.A., de Borst, R., 1984. Non-associated plasticity for soils, concrete and rock, Heron ed., volume 29, pp. 1–75.

Xu, W., Jónsson, S., Ruch, J., Aoki, Y., 2016. The 2015 Wolf volcano (Galápagos) eruption studied using Sentinel-1 and ALOS-2 data: The 2015 Wolf Eruption in the Galápagos. Geophysical Research Letters 43, 9573–9580. URL: http://doi.wiley.com/10.1002/2016GL069820, doi:10.1002/2016GL069820.

Zhan, Y., Gregg, P.M., 2019. How Accurately Can We Model Magma Reservoir Failure With Uncertainties in Host Rock Rheology? Journal of Geophysical Research: Solid Earth 124, 8030–8042. URL:
https://onlinelibrary.wiley.com/doi/10.1029/2019JB018178,
doi:10.1029/2019JB018178.

Zinke, J., Reijmer, J., Thomassin, B., 2003. Systems tracts sedimentology in the lagoon
of Mayotte associated with the Holocene transgression. Sedimentary Geology 160, 57–
79. URL: https://linkinghub.elsevier.com/retrieve/pii/S0037073802003366,
doi:10.1016/S0037-0738(02)00336-6.

**Appendix A. Equations solved, flow laws and physical parameters**

*Appendix A.1. Fundamental equations*

The conservation of mass and momentum is written as follows (Jaeger et al., 2007):

$$
\left\{
\begin{array}{l}
\nabla \cdot \boldsymbol{\sigma} + \boldsymbol{F}_{\nu} = 0 \\
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \boldsymbol{v}) = 0
\end{array}
\right.
\qquad (\text{A.1})
$$

where $\sigma$ is the stress tensor, $\nu$ the velocity and $F_\nu$ is the sum of body forces (in the present work, it is equal to the gravitational body force $\rho g$).

Constitutive laws define the rheology of the domain, i.e. the relationship between the stress tensor **σ** and the strain tensor **ε**. For linear elastic materials (here, we assume the magma reservoirs to behave elastically), we assume Hooke's law (Jaeger et al., 2007):

$$
E\epsilon = (1 + \nu)\sigma - \nu\text{tr}(\sigma)\mathbf{I} \quad (A.2)
$$

with $\epsilon$ the strain tensor, **I** the 3 × 3 identity matrix, $\nu$ Poisson's ratio and *E* Young's modulus.
The notation "tr" refers to a tensor's trace, i.e. the sum of its diagonal components. We assume
that the bedrock surrounding the magma reservoirs (e.g. the lithosphere) is elasto-plastic, which
means that Hooke's law prevails (Equation A.2) until a failure threshold is reached.

Appendix A.2. Shear stress and failure criterion

Shear failure (i.e. “mode II” or “mode III”) is assumed to be triggered by a Drucker-
Prager criterion:

$$
F_{DP} = \sqrt{J_2(\boldsymbol{\sigma})} + \alpha I_1(\boldsymbol{\sigma}) - k = 0 \quad (A.3)
$$
where

$$ \alpha = \frac{2\sin(\varphi)}{\sqrt{3}(3-\sin(\varphi))} \quad \text{and} \quad k = \frac{2\sqrt{3}\cos(\varphi)}{(3-\sin(\varphi))} C \qquad (A.4) $$

with C the cohesion and $\varphi$ the internal friction angle. $F_{DP}$ is the Drucker-Prager yield function. $I_1(\sigma) = \text{tr}(\sigma)$ is the first invariant of the stress tensor, equal to three times the mean normal stress (counted positively for tensile stresses and negatively for compressive stresses). $J_2(\sigma)$ is the second invariant of the stress deviator $\mathbf{s}$, such that

$$ \mathbf{s} = \boldsymbol{\sigma} - \frac{I_1(\boldsymbol{\sigma})}{3} \mathbf{I} \qquad (A.5) $$

and

$$ J_2(\boldsymbol{\sigma}) = \frac{1}{2} \boldsymbol{s} : \boldsymbol{s} = \frac{1}{2} s_{ij} s_{ij} = \frac{1}{2} (s_{11}^2 + s_{22}^2 + s_{33}^2) + s_{12}^2 + s_{13}^2 + s_{23}^2 \qquad (A.6) $$

A Drucker-Prager yield envelope corresponds to a regular cone in the three-dimensional stress space, compared to a Mohr-Coulomb yield envelope which was developed for two-dimensional approaches and displays sharp hexagonal corners (hence discontinuous domains in 3D). A comparison of both criteria can be found in Wojciechowski (2018).

The “Von Mises stress”, here referred to as “shear stress” and denoted $\tau$, is defined as

$$ \tau = \sqrt{3}J_2 = \sqrt{\frac{3}{2}\boldsymbol{s} : \boldsymbol{s}} \qquad (A.7) $$

In elastic-only simulations, the distribution of $\tau$ is often used to determine the areas prone to shear failure (e.g. Gudmundsson, 2007; Browning et al., 2021; Karaoğlu et al., 2020).

Tensile failure (i.e. “mode I”) is classically thought to occur when the mean normal stress $I_1(\sigma)$ exceeds the rock’s tensile strength $T$, which generally ranges between 1 and 10 MPa (e.g. Grosfils et al., 2015, and references therein). In our study, we do not account for tensile failure since the confining pressure at the studied depths is so large that stresses around our modeled reservoirs were never at tensile state. Besides, laboratory experiments, analytical studies, numerical modeling and field observations suggest that, for depths greater than ≈ 2 km, shear failure occurs prior to tensile failure (Gerbault, 2012; Grosfils et al., 2015, and references therein). Let us recall that double-couple earthquakes are associated to shear failure, and are often observed during volcanic eruptions (e.g. Levy et al., 2018). Hence, we consider the absence of tensile failure in our models to be a minor drawback.
## Appendix A.3. Expected critical pressure for wall failure

The Drucker-Prager criterion predicts that plasticity occurs when Equation A.3 is satisfied, i.e. when:

$$ \sqrt{J_2(\boldsymbol{\sigma})} = k - \alpha I_1(\boldsymbol{\sigma}) \qquad (A.8) $$

Assuming that the confining pressure is nearly lithostatic yields:

$$ I_1(\boldsymbol{\sigma}) = -3\rho g z \qquad (A.9) $$

with z counted positively downwards. McTigue (1987) and Segall (2010) give the stress field, in spherical coordinates $(r, \theta, \varphi)$, for a spherical reservoir of radius $a$ in an unbounded elastic region:

$$ \boldsymbol{\sigma} = \begin{pmatrix} -\left(\frac{a}{r}\right)^3 & 0 & 0 \\ 0 & \frac{1}{2}\left(\frac{a}{r}\right)^3 & 0 \\ 0 & 0 & \frac{1}{2}\left(\frac{a}{r}\right)^3 \end{pmatrix} \times \Delta P \qquad (A.10) $$

with the convention that tensile stresses are positive. On the reservoir wall ($r = a$), this yields:

$$ \boldsymbol{\sigma} = \begin{pmatrix} -1 & 0 & 0 \\ 0 & \frac{1}{2} & 0 \\ 0 & 0 & \frac{1}{2} \end{pmatrix} \times \Delta P \qquad (A.11) $$

This tensor is purely deviatoric, and its second invariant is:

$$ \sqrt{J_2(\boldsymbol{\sigma})} = \sqrt{\frac{1}{2}\boldsymbol{\sigma} : \boldsymbol{\sigma}} = \frac{\sqrt{3}}{2} \times \Delta P \qquad (A.12) $$

This solution agrees with our numerical tests. Our modeled reservoirs seem far enough from the free surface for Equation A.10 to be valid.

Substituting Equations A.4, A.9 and A.12 into Equation A.8 gives the value of the critical pressure at which failure occurs:

$$ \Delta P_{\text{crit}} = \frac{4}{(3 - \sin(\varphi))} (\cos(\varphi)C + \sin(\varphi)\rho g z) \qquad (A.13) $$

Given the large reservoir radius-to-depth ratio considered, we expect that, for a spherical reservoir, failure initiates at its apex (Grosfils et al., 2015). Considering a center depth $d$ and a radius $r$, we thus calculate $\Delta P_{\text{crit}}$ at depth $d - r$. Figure S1 summarizes our parameter
exploration. In particular, for C = 10 MPa and φ = 1, d = 40 km and r = 4 km, we expect ΔP$_{crit}$ ≈ 40 MPa. For a reservoir with d = 25 km and r = 3 km, assuming the same C and φ, we expect ΔP$_{crit}$ = 30 MPa. Given its shallower depth, it is not surprising that the upper reservoir reaches failure for smaller applied pressures than the lower reservoir. Analytical works generally examine the case of an overpressure, but we find that applying an underpressure does not significantly affect the failure threshold.

## Appendix A.4. Non-associated plastic flow rule

When the failure criterion is met, brittle behavior is simulated by a non-associated plastic flow rule. For this rule, the strain rate tensor ˙ε is split between an elastic and a plastic component (Vermeer and de Borst, 1984):

$$ \dot{\epsilon} = \dot{\epsilon}^{el} + \dot{\epsilon}^{pl} \qquad (A.14) $$

The elastic part ˙ε$^{el}$ follows Hooke's law while the plastic part ˙ε$^{pl}$ is the derivative of the plastic potential Q such that:

$$ \dot{\epsilon}^{pl} = \lambda \frac{\partial Q}{\partial \sigma} \qquad (A.15) $$

where λ is a positive proportionality factor. Details on its derivation can be found in Vermeer and de Borst (1984)] or Gerbault et al. (1998). The plastic potential Q is defined similarly to the stress yield function (Equations A.3 and A.4), except that the internal friction angle φ is replaced by the dilatancy angle ψ:

$$ Q = \sqrt{J_2(\boldsymbol{\sigma})} + \alpha_Q I_1(\boldsymbol{\sigma}) - k_Q \qquad (A.16) $$

with

$$ \alpha_Q = \frac{2\sin(\psi)}{\sqrt{3(3-\sin(\psi))}} \quad \text{and} \quad k_Q = \frac{2\sqrt{3}\cos(\psi)}{(3-\sin(\psi))} C \qquad (A.17) $$

The case ψ = φ is referred to as “associated plasticity”, while the case ψ ≠ φ corresponds to “non-associated plasticity” (Vermeer and de Borst, 1984). This latter framework is more appropriate for soils, rocks and granular materials in general (Vermeer and de Borst, 1984), and is the most complete approach in predicting propagating fault patterns (Gerbault et al., 1998). Choosing ψ = 0 yields a theoretically incompressible flow within the plastic areas,
and can thus produce localized shear bands without any preexisting strength heterogeneities (Gerbault et al., 1998).

## Appendix A.5. The interpretation of plastic shear strain

In this study, we have chosen to evaluate the occurrence of shear failure by examining the distribution of the second invariant of the plastic strain deviator, here noted $\epsilon_{II}^{pl}$:

$$ \epsilon_{II}^{pl} = \sqrt{\frac{1}{2} \mathbf{dev}(\epsilon^{pl}) : \mathbf{dev}(\epsilon^{pl})} \qquad (A.18) $$

where $\mathbf{dev}(\epsilon^{pl})$ is the deviator of the plastic strain tensor $\epsilon^{pl}$:

$$ \mathbf{dev}(\epsilon^{pl}) = \epsilon^{pl} - \mathrm{tr}(\epsilon^{pl})\mathbf{I} \qquad (A.19) $$

We refer to $\epsilon_{II}^{pl}$ as “plastic shear strain”.

Other components of the plastic strain tensor provide additional information. For example, the volumetric part shows areas of compression and dilatation, which can indicate the occurrence of tensile failure. Here, we study plastic shear strain since it better shows the amount of distortion generated by failure. In fact, tests we made show that the volumetric part of the plastic strain tensor generally followed the same distribution as $\epsilon_{II}^{pl}$, only one or two orders of magnitude below. Most of the plastic strain generated by our simulations is therefore plastic shear strain.

## Appendix B. Numerical method implemented by COMSOL

Here, we present a succinct overview of the technical aspects of the solver functionalities we use in COMSOL. Full details can be found in the official COMSOL documentation.

COMSOL Multiphysics® (v6.1) uses finite elements to solve ordinary differential equations over a user-defined mesh. We use the so called Stationary Study, which implements damped Newton iterations. The equations solved can be expressed as:

$$ f(U) = 0 \qquad (B.1) $$
with $U$ the solution vector and $f(U)$ the residual. The solver makes an initial guess $U_0$, and calculates its Jacobian $f'(U)$ such that:

$$f'(U_0)\delta U = -f(U_0) \qquad (B.2)$$

The next iteration $U_1$ is then calculated as follows:

$$U_1 = U_0 + \lambda \delta U \qquad (B.3)$$

with $\lambda \in [0,1]$ the damping factor. A new residual is calculated; if it is lower than that of $U_0$, $U_1$ is accepted and iterations continue. If it is greater, $U_1$ is computed again with a smaller damping factor. Iterations are terminated either when the minimal set value of the residual is reached (i.e. the iterations stop since a satisfying solution has been found) or when the maximum set number of iterations is reached (i.e. the solver has not converged after a given number of iterations, hence it stops).

We impose reservoir wall pressures progressively, as a succession of quasi-static steps. To do so, we select in the Stationary Study, the “auxiliary parametric sweep” option. Each pressure increment is solved with damped Newton iterations as described above. We define a maximum reservoir pressure, select an initial pressure step size, and allow the solver to automatically adjust the step size if necessary. In other words, when convergence cannot be achieved with the initial pressure step size, the solver tries again with a smaller pressure step. In our simulations, a drastic reduction in pressure step size typically characterizes the onset of the large deformation phase (see Results and Figure 4).

The calculations stop either (1) when all pressure increments have been successfully computed (i.e. the prescribed maximum value has been reached), or (2) when no solution with a low enough residual was found even when using the smallest defined pressure step size.

## Appendix C. Comparison of results with meshed and unmeshed magma reservoirs

In this study, we have chosen to model magma reservoirs as more compliant elastic subdomains, following the approach of, e.g., Gerbault et al. (2012) and Got et al. (2013). As described in Section 3.3, we have chosen our mechanical parameters so that the surface displacement in our models satisfyingly matched analytical predictions of calculations based on empty reservoirs (McTigue, 1987). In fact, other numerical modeling studies present models of
empty (i.e. unmeshed) magma reservoirs in an elastic (e.g., Grosfils et al., 2015; Cabaniss et al., 2020) or elastoplastic (Gerbault et al., 2018; Ruz-Ginouves et al., 2021; Novoa et al., 2022) domain.

The supplement of Gerbault (2012) shows a benchmark comparing both approaches. Below, we also propose a comparison between models with meshed and unmeshed reservoirs. We perform our tests with two spherical reservoirs, in the same framework used in Section 4.3 (central depths 25 and 40 km, radii 3 and 4 km). We present the evolution of the vertical surface displacement and reservoir volume variation for $\frac{\Delta P_{up}}{\Delta P_{low}} = -1, -0.3, 0, 0.3$ and $1$ (Figure S8). We also show the evolution of the plastic shear strain for the two extreme cases (Figures S9 and S10).

Figure S8 shows that models with meshed and unmeshed magma reservoirs produce very similar results in the early simulation stages. The difference between them progressively increases, since the unmeshed models seem to reach the rapid reservoir volume variation stage for slightly lower applied pressures. Without a meshed reservoir, the solver is also able to converge until larger values of surface displacement and volume variation are reached.

Figures S9 and S10 also show that, for the displayed cases, the plastic shear strain fields produced are very similar between the two approaches. In particular, the generation of an axial shear zone for two deflating reservoirs and a conical shear zone with reservoirs of opposite signs remains valid.

The conclusions of our study therefore seem uninfluenced by the choice of meshed or unmeshed magma reservoirs. In the latter approach, the solver seems able to converge even when very large amounts of deformation are produced (several tens of meters in ground displacement and tens of km³ in reservoir volume variation). This results in a longer simulation runtime, mostly to compute steps that are likely unrealistic when studying the Mayotte 2018-2021 eruption. For that reason, we have opted for the approach presented in our study.

**Supplementary Reference**

Wojciechowski, M., 2018. A note on the differences between Drucker-Prager and Mohr-Coulomb shear strength criteria. Studia Geotechnica et Mechanica 40, 163–169. URL:
Appendix D. Supplementary figures

<figure><img src="image_10.png" /><figcaption>Figure S1: Parametric exploration of the critical pressure required for reservoir wall failure in an axisymmetric setting. A, B: tests with a spherical reservoir, varying the friction angle $\phi$, corresponding to our modeled dimensions for the upper (A) and lower reservoir (B). C: tests for a spherical reservoir varying the cohesion $C$. D: tests with varying reservoir aspect ratios, from oblate to prolate shapes. Analytical solutions for A-C are calculated following Equation A.13.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

Figure S2: Effect of the reservoir aspect ratio on $\tau$ and $e_{II}^{pl}$ distributions with a single deflating reservoir. A-B: oblate reservoir ($h/w = 1/2$). C-D: sill-like reservoir ($h/w = 1/4$). The host rock is elastic (A, C), or elasto-plastic with $C = 10$ MPa and $\phi = 1^\circ$ (B, D). For elastic simulations, the applied $\Delta P$ is -250 MPa, whereas for elasto-plastic simulations, the last pressure step before solver outbreak is represented (hence the different values of $\Delta P$).

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

Figure S3: Displacement field produced by the reference simulation with two reservoirs (Section 4.2, Figure 4). A: plastic shear strain (Figure 4C). B: total displacement. C: horizontal displacement (positive outwards). D: vertical displacement (positive upwards). The conical plastic strain area delimits an area of subsidence between the reservoirs and accommodates inverse movement.
<figure><img src="image_13.png" /><figcaption>Figure S4: Plastic shear strain generated by different $\Delta P_{up}/\Delta P_{low}$ ratios, for two spherical reservoirs. Each subplot represents the final stage of a run shown as dots in Figure 5A. Subplots A, K and U correspond to Figures 5B, C and D, respectively.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

Figure S5: Effect of the distance between reservoirs on the plastic strain field, for oblate reservoirs ($h/w = 1/2$). For every case, the presented step is the last before the solver fails to converge (hence the run stops). $\Delta P_{up}/\Delta P_{low}$ ratios are 0.3 for the top row, 0 for the middle row and -0.3 for the bottom row. $d_{up}$ and $d_{low}$ differ from those of Figure 6 to ensure that reservoir separations remain the same (see the Table 1 caption). Horizontal and vertical axes display radial distance and depth, in km. Reservoir separations are identical to Figure 6. Principal stress orientations are shown similarly to Figure 4.
<figure><img src="image_15.png" /></figure>

Figure S6: Plastic shear strain generated by different $\Delta P_{up}/\Delta P_{low}$ ratios, for two sill-like reservoirs. Each subplot represents the final stage of a run shown as dots in Figure S7. Subplots A, K and U correspond to Figures S7B, C and D, respectively. Reservoir separations (height between the roof of the lower reservoir and the base of the upper reservoir) are identical to Figure S4.
<figure><img src="image_16.png" /></figure>

Figure S7: Effect of magma reservoir pressurization, for two sill-like reservoirs. A: critical pressures for the different stages on a set of runs with varying $\Delta P_{up}/\Delta P_{low}$. Solid blue curve: applied critical magma pressure at which plastic deformation initiates at the lower reservoir wall. Dashed blue curve: same for the upper reservoir. Red curve: onset of stage 4 (Figure 4). Here, unlike with spherical reservoirs (Figure 5), the 5-km$^3$ reservoir deflation is never reached. B-D: final stages of the end-member simulations (all intermediate cases are presented on Figure S6).
<figure><img src="image_17.png" /></figure>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Figure S8: Vertical surface displacement and reservoir volume variation for simulations with meshed and unmeshed reservoirs. Left column: evolution of the vertical surface displacement. Dark red curves: last pressure increment (hence the different values between meshed and unmeshed reservoirs). Middle column: evolution of the vertical surface displacement at radial distances $r$ 10 km and 20 km (since Mayotte lies between 10 and 20 km from the proximal cluster). The gray patch highlights the 10-20 cm subsidence
observed on Mayotte. The circles show the displacement obtained when $\Delta V_{low} = -5 \text{ km}^3$. Right column: reservoir volume variations during the simulations.

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

Figure S9: Evolution of the plastic shear strain throughout a simulation with $\frac{\Delta P_{up}}{\Delta P_{low}} = 1$, with (A) and without (B) a meshed subdomain in the magma reservoirs. The pressure steps are chosen so that $\Delta V_{low}$ are similar between the top and bottom row (approximately -1 km$^3$, -3 km$^3$, -5 km$^3$)
<figure><img src="image_33.png" /><figcaption>Figure S10: Same as Figure S9, but with $\frac{\Delta P_{up}}{\Delta P_{low}} = 1$.</figcaption></figure>