# Volcano-tectonic crisis of Mayotte (2018-2022): a deformation model for geodetic applications

Raphael Grandin, Xavier Collilieux, Pasquier Inès, Oliver Jamet

► To cite this version:

Raphael Grandin, Xavier Collilieux, Pasquier Inès, Oliver Jamet. Volcano-tectonic crisis of Mayotte (2018-2022): a deformation model for geodetic applications. Journal of Geophysical Research: Solid Earth, In press, 10.1029/2025JB031473. [https://doi.org/10.5295/873]

HAL Id: hal-05295873

https://hal.science/hal-05295873v1

Submitted on 3 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Volcano-tectonic crisis of Mayotte (2018–2022): a deformation model for geodetic applications

Raphaël Grandin<fnref n="1" />, Xavier Collilieux<fnref n="1" />,<fnref n="2" />, Inès Pasquier<fnref n="2" />, Olivier Jamet<fnref n="3" />

*Accepted for publication*
*on 2 October 2025*

<footnote n="1">¹Université de Paris Cité, Institut de physique du globe de Paris, CNRS, IGN, F-75005 Paris, France</footnote>

<footnote n="2">²Univ Gustave Eiffel, ENSG, IGN, F-77454 Marne-la-Vallée, France</footnote>

<footnote n="3">³Institut national de l'information géographique et forestière, IGN, F-94160 Saint-Mandé, France.</footnote>
**Keywords (up to 6)**

Mayotte / volcano deformation / GNSS / InSAR / levelling / terrestrial reference frame

**Index terms (1 to 5)**

• 1207 Transient deformation (6924, 7230, 7240)

• 1229 Reference systems

• 1928 GIS science

• 8155 Plate motions: general (3040)

• 8439 Physics and chemistry of magma bodies

**Key Points:**

• A deformation model is derived for Mayotte from inversion of permanent GNSS and InSAR data, validated against levelling and GNSS campaigns

• Model accuracy is better than 2 cm, allowing for recovering tectonic motion of the island and detecting local deformation at Hamouro

• The model is translated into a semi-kinematic frame for coordinate epoch trans-formation of geodetic data spanning the crisis
**Abstract**

Volcano-tectonic unrest during the 2018–2022 crisis of Mayotte (Comoros archipelago) generated an island-wide surface displacement reaching ~ 20 cm, with a relative (inter- nal) deformation of the island of ~ 10 cm. Correcting for these displacements, which ex- ceed the required accuracy of geodetic network reference coordinates, is necessary for Ter- restrial Reference Frame applications, both for global (e.g. ITRF2020) or local reference frame determinations.

Here, we develop a deformation model defined by inversion of permanent GNSS and InSAR measurements, treating spatial and temporal components of deformation sepa- rately. The spatial pattern of deformation is reproduced by a near-vertical deflating cigar- shaped pressure source situated ~ 20 km to the east of the island and ~ 24 km depth. The temporal evolution of displacement is modeled by a double-exponential function, which likely reflects the progressive drainage of the reservoir feeding the eruption site on the seafloor.

The model effectively removes non-linear variations in permanent GNSS coordi- nates time-series, with residuals on daily positions everywhere smaller than 2 cm. Model- corrected campaign GNSS, levelling and InSAR data acquired independently allow for detecting a previously unidentified local deformation signal, unrelated to the 2018–2022 crisis, in the area of Hamouro.

The deformation model is incorporated in the *PROJ* library, using the formalism introduced by Crook (2020), being the first volcanic deformation model of this kind. We recommend the use of this model to supplement the new local terrestrial reference frame for Mayotte (RGM23) especially for analyzing geodetic data acquired between mid-2018 and 2023, thus realizing a de facto semi-kinematic frame for Mayotte.

**Plain Language Summary**

Starting from May 2018, a major submarine volcanic eruption took place offshore of the island of Mayotte (Comoros archipelago), emitting over 6 cubic kilometers of lava on the seafloor. The vertically-extensive magma reservoir feeding the eruption was par- tially drained, generating a deformation of the crust over distances of tens of kilometers. Over the course of the eruption (2018–2022), the island was progressively displaced down- wards and horizontally by 20 cm. This deformation induced an island-wide change in
the geographic coordinates of all points situated on the island. This motion needs to be
compensated in order to recover the “true” coordinates of geodetic benchmarks, private
properties or public infrastructures. Here, we describe the magma reservoir using a me-
chanical model, with internal parameters constrained using spatially-scattered GNSS data,
which enables us to reconstruct the full deformation field that affected the island. This
model is validated against ground-truth data collected by a geodetic campaign in 2023.
Stability of sites surveyed in 2023 is retrospectively assessed by comparing model pre-
dictions and new measurements, complemented by information provided by satellite radar
interferometry. Model accuracy is estimated to approximately 2 cm at all times, which
is adequate to compensate the motion caused by the volcanic crisis.

Introduction

A geodetic datum, or more specifically a Terrestrial Reference Frame (TRF), must
be adopted to unambiguously and precisely define the coordinates of any object at the
Earth's surface. There is currently a wide set of available TRFs with spatial extents cov-
ering a country, a continent or the entire Earth (e.g. the ITRF2020, Altamimi, Rebis-
chung, et al., 2023). TRFs are generally defined by the coordinates of a network of geode-
tic markers, permanent GNSS antennas or more occasionally through a coordinate trans-
formation. They are most commonly issued by an authoritative agency (typically, a na-
tional land survey or an international working group), with their associated metadata
made available in public databases such as EPSG (2024).

As the Earth is deformable, coordinates of both reference and target markers must
be allowed to change with time (McCarthy et al., 2004). In the stable interior of tectonic
plates, points show negligible relative horizontal coordinate variations with time. Thus,
it is possible to define a “regional” TRF with respect to which coordinate variations are
negligible over a decade (e.g. the ETRF2020, the latest realization of the European Ter-
restrial Reference System 1989 ETRS89, Altamimi & Collilieux, 2024). In such a “plate-
fixed” frame, surveyors can exchange coordinates without necessarily adding a time tag,
hereafter named epoch.

At global scale, the primary contribution to coordinate changes is caused by the
steady motion of tectonic plates. Relative plate motion can be described by angular ve-
locities around a set of Euler poles, usually determined geodetically (D. F. Argus et al.,
2010; Altamimi, Métivier, et al., 2023) or formerly geologically (DeMets et al., 2010). Knowledge of the appropriate Euler rotation parameters suffices to calculate a velocity vector at any point or, equivalently, a displacement between two epochs, the latter being simply proportional to the time difference between measurements. Such a model can be easily incorporated into a 14-parameters transformation by modifying the rotation rate parameters in such a way that a plate-fixed frame can be defined rigorously from a global no-net-rotation TRF (e.g. Altamimi & Collilieux, 2024).

On the other hand, in regions subject to local or regional deformation processes, standard coordinate transformations cannot provide an operational solution. Some deformation processes can still be modeled by a temporally constant velocity field, such as Glacial Isostatic Adjustments (GIA) (D. Argus & Drummond, 2015), or interseismic deformation near plate boundaries (Kreemer et al., 2014). Other processes generating transient deformation, such as seismic cycle processes, volcanic activity, accelerated ice melting or continental hydrology, can cause non-linear coordinate time-variations. Thus, the coordinates measured in the field by surveyors may have a complex dependence on time. Since in most civil engineering, cadastral or cartographic applications, coordinates have to be expressed at the same epoch, coordinate epoch transformations are mandatory to achieve precise geolocation in deformation zones. Performing coordinate epoch transformations requires a deformation model which covers the relevant area and period of interest. Such models may be occasionally available from case studies published in the scientific literature, but their dissemination into the public domain for operational applications is not yet common practice.

Achieving this objective requires the definition and implementation of a standard exchange format. New Zealand was the first country to adopt a deformation model in the definition of its national TRF. Geodesists have established a so-called “semi-kinematic frame” where users are requested to publish all coordinates at the same reference epoch (Donnelly et al., 2017; Stanaway, 2020). Because the New Zealand territory straddles a major plate boundary, all geodetic markers are potentially affected by seismic-cycle processes generating non-linear displacements, often decomposed into static co-seismic displacement and transient post-seismic deformation (see e.g. the case of the 2016 Kaikōura earthquake described by Hamling, 2020). Therefore, in New Zealand, a deformation model that takes into account the contribution of recent earthquake activity needs to be applied by the coordinate producers in order to propagate coordinates at the reference epoch
(Donnelly et al., 2017; Stanaway, 2020). Thanks to this correction, users and software
packages such as Geographic Information Systems (GISs) do not have to rigorously man-
age the coordinate epoch.

The proposed deformation model for New Zealand, NZGD2000, consists of two com-
ponents: (i) a steady velocity field that represents the non-uniform spatial distribution
of interseismic motion, using a grid, and (ii) occasional updates that describe the local
deformation field caused by large tectonic events, commonly referred to as “local patches”
(Blick et al., 2005). These patches are regularly updated, published and quality-checked
(Denys & Pearson, 2022). Similar initiatives for Iceland and Fenno-Scandinavia include,
e.g., those reported in Kierulf et al. (2019) or Häkli et al. (2023, 2016). Motivated by
the growing precision and coverage of modern geodetic measurements, the establishment
of a methodology for defining a deformation model at the national or continental scale
for TRF applications has been an active research topic for two decades (Stanaway et al.,
2024).

In this paper, we use a geodetic dataset acquired over the island of Mayotte (Comoros archipelago) to develop a case-study illustrating a methodology for setting up and qualifying a deformation model for operational geodesy. Geodetic data was acquired before, during and after a transient volcano-tectonic deformation episode that unrolled between 2018 and 2022, which we aim at characterizing and modeling in space and time. The deformation model aims to facilitate the accurate computation of coordinate changes for operational coordinate transformations, and establish the foundations for a semi-kinematic frame applicable to the island. Correction of the deformation transient also allows us to investigate the existence of local deformation unrelated to the crisis, and improve the understanding of long-term motion of the island. The deformation model also gives us the opportunity to study the MAYG GNSS coordinate time series as its enhanced modelling would improve ITRF2020 results (Altamimi, Rebischung, et al., 2023). Finally, even if the procedure involved in modelling and correcting for deformation will inevitably be case-specific, this study aims at paving the way towards a general method that could be applied more systematically to other geological transients, past or future.

The paper is organized as follows. Section 1 provides a background on the 2018–
2022 volcano-tectonic crisis in Mayotte. In Section 2, we present the dataset used in this
study. In order to constrain the deformation field associated with the 2018–2022 volcano-
tectonic episode in Mayotte, we rely on daily positions from nine permanent GNSS sta-
tions installed over the island, complemented by two InSAR time-series that captured
the inception of the deformation transient. This primary geodetic dataset (hereafter re-
ferred to as the “reference dataset”) is assimilated in Section 3 for setting up a defor-
mation model of the whole island. In Sections 4 to 6, we exploit an independent dataset
(hereafter referred to as the “validation dataset”) to estimate the model quality and an-
alyze any residual deformation. In Section 4, uncertainty associated with the deforma-
tion model is estimated by assessing posterior residuals of permanent GNSS position and
baseline time-series. In Section 5, using novel data acquired during a GNSS and level-
ling campaign conducted in 2023, we apply the deformation model to correct for the cu-
mulative displacement induced by the volcano-tectonic crisis, and resolve local deforma-
tion features unrelated to the crisis. In Section 6, residual velocities from model-corrected
permanent GNSS and tide gauge data are used for providing an estimation of the long-
term tectonic motion of the island and the local rate of sea level rise. Finally, in Section 7,
we discuss the limitations of the model, draw perspectives for its improvement and pro-
mote its use in geodesy.

# 1 Geological and geodetic context of the Mayotte crisis

The island of Mayotte (Comoros archipelago) is a 20 km × 40 km tropical island that belongs to an area of diffuse tectonic deformation and active volcanism, at the boundary between the Somalia and Lwandle tectonic plates (Michon, 2016; Famin et al., 2020; Stamps et al., 2021). Prior to 2018, the previous episode of volcanism has been dated to 4–6 kyr BP (Zinke et al., 2003).

From 2018 to 2022, an episode of intense volcano-tectonic activity took place off-
shore of Mayotte, generating a transient deformation pattern that affected the whole is-
land (Figure 1a). Deformation was caused by deflation of a crustal-scale magma stor-
age system, situated ~ 20 km ESE of Mayotte, extending vertically in the 20–50 km depth
range (Lemoine et al., 2020; Cesca et al., 2020), corresponding to a volume loss of sev-
eral cubic kilometers (Peltier et al., 2022). Ascent of magma and subsequent drainage
of the reservoir were likely primed by a rifting event that led to the formation of the ~ 6.5 km³
Fani Maoré submarine volcano at 3000 m bsl., ~ 50 km ESE of Mayotte (Feuillet et al.,
2021), although the precise architecture of the plumbing system remains poorly constrained.
The episode started with an earthquake swarm in May 2018, strongly felt by the population (Lemoine et al., 2020). From July 2018, an accelerating deformation anomaly was reported, characterized by eastward and downward movement of the island. The velocity reached a maximum of ~ 2 cm/month around January 2019, and progressively decayed in the following years, practically falling to an unmeasurable level in 2022. Over 3 years, the episode generated a cumulative absolute displacement of ~ 20 cm, with a tilt and internal deformation of the island exceeding 10 cm (Figure 1c).

The 2018–2022 transient deformation episode in Mayotte represents a challenge for the definition of a local reference frame, since the island can no longer be considered as behaving rigidly. Prior to 2018, the island was equipped with a single permanent GNSS station (MAYG), which only allowed for defining a single velocity vector for the whole

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Cumulative horizontal and vertical GNSS displacement induced by the 2018–2022 volcano-tectonic crisis. Inset shows a simplified representation of magma path. (b) Geodetic dataset used in this study. (c) GNSS displacements 2014–2024.</figcaption></figure>
island. In addition to MAYG, the deformation transient was also detected by a set of
five other GNSS stations installed, by chance, just a few months before the unexpected
inception of the volcano-tectonic episode (Figure 1a). However, due to lack of pre-crisis
geodetic data, the “normal” motion of the island is poorly defined, thus, in turn, mak-
ing it hard to determine the actual date of termination of the deformation. In other words,
it remains unclear when the island can be considered to have regained its stable state.
This information is however mandatory for a correct interpretation of geodetic data ac-
quired since 2022, especially to better constrain tectonic plate models in the East African
region.

Defining an updated local reference frame also has practical implications beyond
the geophysical or geodetic communities. In Mayotte, resolving anomalous ground dis-
placement at < 10 cm accuracy horizontally is essential for surveyors to be able to de-
termine surveyed coordinates at the centimeter-level. Furthermore, a high precision on
the vertical component is crucial for anticipating the impacts of coastal inundation on
infrastructures at high tides, especially in the context of a global sea level rise (Jeanson
et al., 2021; Charroux et al., 2024). Finally, constraining vertical motion of the fringing
reef and lagoon is important for monitoring the evolution of sensitive wildlife, such as
mangroves (Aubry & Tempere, 2022) or breeding sites for sea turtles (Morisseau, 2024).

## 2 Geodetic Data

### 2.1 Reference dataset

#### 2.1.1 *Permanent GNSS daily positions*

We first analyze data from nine permanent GNSS stations installed on the island
of Mayotte, that were operational during the 2018–2022 episode (Figures 1): (i) only one
permanent GNSS station was operating well before the crisis (MAYG, installed in 2014);
(ii) five additional stations were installed just before the onset of the deformation tran-
sient in May 2018 (BDRL, PORO, KAWE, MTSA, GAMO, installed at the beggining
of 2018); (iii) finally, three remaining stations were installed in response to the crisis (PMZI,
KNKL and MTSB, installed in 2019). GNSS position time-series with respect to IGb14
reference frame are computed daily by OVPF with Precise Point Positioning using Gip-
syX software (see Peltier et al., 2022, for details). Stations located outside of Mayotte
(the closest being GLOR, ~ 200 km to the NE of the Fani Maoré volcano) are too far
to detect the deformation transient, and are not used in this study. The time-series are
only considered until March 2022, when the deformation model was developed. Subse-
quent observations are not integrated into the deformation model, which has remained
unchanged since mid-2022.

Daily GNSS position time-series are represented in Figures 1c and S1. From 2014
to 2018, tectonic motion is only clearly visible at station MAYG, which provides the longest
pre-crisis record. After detrending the time-series, the volcano-tectonic transient is de-
tected by all GNSS stations of the network, producing a nearly-uniform motion toward
the east (~22 cm between 2018 and 2022), combined with a North-South contraction
of the network (~12 cm relative northward displacement between KAWE and BDRL).
Subsidence is greater in the eastern part of Mayotte (~19 cm at MAYG) than in the
western part (~11 cm at GAMO). Cumulative displacement vectors clearly show a cen-
tripedal pattern, pointing toward a location situated ~25 km to the east of Petite Terre,
with a pitch angle of ~45° (see inset of Figure 1a). This pattern is consistent with a
source of deflation situated offshore the island. This deflation center is however clearly
distinct from the Fani Maoré volcano, which is located further to the south-east, ~50 km
ESE of Petite Terre.

2.1.2 X- and C-band InSAR time-series

Island-wide deformation is also analyzed by Interferometric Synthetic Aperture Radar (InSAR) using Sentinel-1 data acquired in descending track (incidence angle: 39°, 42 Sentinel-1 images between February 2018 and May 2019). No data was acquired in ascending track by Sentinel-1 over Mayotte prior to December 2018, with exception of an isolated acquisition on 11 January 2018, which makes it unsuitable for mapping ground displacement in ascending pass during the first months of activity. Instead, the Sentinel-1 descending-track time-series is complemented by an ascending-track COSMO-SkyMed dataset (incidence angle: 39°, 50 COSMO-SkyMed images between February 2018 et August 2019). In order to overcome the poor interferometric coherence in X- and C-band, these datasets were processed with a Persistent Scatterer approach using the SqueeSAR method (Ferretti et al., 2011). Low-pass temporal filtering was applied to mitigate atmospheric noise. To reduce costs, and because deformation significantly slowed down after 2020, the processing performed in 2020 was not updated. Hence, this dataset stops in 2020. Over the processed area of 374 km², 18,719 Persistent Scatterers (PS) are detected for Sentinel-1, and
244,430 PSs for COSMO-SkyMed. In addition, 11,382 Distributed Scatterers (DS) are identified for Sentinel-1, and 43,844 DSs for COSMO-SkyMed. An area surrounding GNSS station GAMO is chosen for referencing ascending and descending time-series.

The time-series produced by the SqueeSAR processing are represented in Figures S2 and S3. An East-West displacement gradient is observed across the island, reaching a cumulative ~ 100 mm in the ascending geometry and ~ 50 mm in the descending geometry until mid-2019. This difference is explained by the geometry of acquisition, with ascending data (COSMO-SkyMed) being more sensitive to strain affecting the island.

## 2.2 Validation dataset

### 2.2.1 Permanent GNSS

Reprocessed Precise Point Positioning (PPP) position time series of the 9 permanent stations computed with GipsyX following IGS repro3 standards (Rebischung et al., 2024) are used as a first validation dataset. The time-series cover the period from January 2014 to June 2024, extending the reference dataset by about 2 years. They are also expressed with respect to the IGb14 reference frame.

To assess the dependency of the results to a particular processing strategy, an operational solution computed by the French mapping agency IGN using Bernese 5.2 software and double-difference processing is also considered, but only for assessing baseline variations since May 2018. Indeed, this solution only includes regional sites and is not homogeneously reprocessed (PCO/PCV file modified on 3 May 2018) and makes use of less advanced troposphere modeling (GPT/GMF) than the PPP solution. It has been aligned to IGS14 reference frame.

### 2.2.2 GNSS and levelling campaigns

The latest GNSS surveys prior to the 2018 crisis were conducted in the 2000s. A total of 31 campaign sites covering Mayotte (Petite and Grande Terre) were observed. Two simultaneous GPS sessions of at least 2 hours have been carried out for each site. Collected data were processed by Bernese software 4.2 in double-difference mode, except for a few short baselines computed with the commercial software Sky Pro. These specific baselines were combined with Bernese results using least-squares adjustment. They
were aligned to ITRF2000 at epoch 2004.0 using 9 regional permanent GPS stations out-
side Mayotte. The quality of the georeferencing has been evaluated at the level of 9.0 mm
for each horizontal component and 4.9 mm on the vertical component. These coordinates
define the terrestrial reference frame named RGM04 (IGN, 2004), maintained by IGN.

The vertical reference system in use at the time of writing of this manuscript is named MAYO53 (IGN and Shom, 1953). The height above mean sea level, hereafter referred to as “MSL height”, is provided for 434 benchmarks surveyed by spirit levelling in 2006. However, the quality of the gravity measurements acquired at that time was insufficient to rigorously compute orthometric heights. Thus, gravity was ignored in the processing of 2006 MSL heights.

In order to update the terrestrial and vertical reference from Mayotte after the vol-
canic crisis, a new survey was conducted by IGN from September to December 2023. In
addition to new sites, 16 RGM04 geodetic markers were surveyed by means of 2 GNSS
occupations of at least 2 hours. The GNSS coordinates were processed with Bernese ver-
sion 5.2 and expressed in IGS20 (consistent with ITRF2020) at epoch 2023.75 using 8 re-
gional permanent stations outside Mayotte. The difference between campaign GNSS-derived
coordinates between the 2004 and 2023 surveys, after correction of plate motion (see sec-
tion 5.1), is presented in Figure S4, confirming the centripetal pattern of horizontal vec-
tors, eastward tilt and subsidence of the island.

During the same 2023 campaign, about 260 km of levelling measurements were also
acquired, including loops to check closure. A total of 293 benchmarks in common with
the 2006 survey have been measured. A group of 26 benchmarks were discarded due to
suspicion of local motion following on-site visual inspection by IGN surveyors, leading
to the selection of a subset of 267 common benchmarks. While new absolute and rela-
tive gravity measurements have been collected, they have not been used to process MSL
height in this study for a self-consistent comparison with 2006 MSL heights. The raw
difference of 2006 and 2023 MSL heights is presented in Figure S5a. Heights are com-
puted with respect to the common benchmark N-402-BIS located on the West of the Pe-
tite Terre island before computing the differences: this reference benchmark consequently
displays no height change in Figure S5a.
### 2.2.3 L-band InSAR time-series

A dataset consisting of 56 images acquired by the ALOS-2 satellite is processed for the period between July 2015 and June 2024 (Figure S6). Data is acquired in ascending track, in Stripmap (SM3) mode. Thanks to the good interferometric coherence in L-band, this dataset was analyzed using a standard Small Baseline Subset approach. Processing was carried out after the completion of the model inversion in 2022, hence is not used for constraining its parameters. Interferograms were computed with the ISCE software (Rosen et al., 2018), unwrapping was performed using SNAPHU (Chen & Zebker, 2002), and a time-series of ground displacement was calculated using the NSBAS software (Doin et al., 2011) after referencing all interferograms with respect to an arbitrary reference point situated to the west of Mayotte, and coinciding with GAMO GNSS station. The time-series of interferograms is displayed in Figures S7 and S8.

### 2.2.4 Ocean-bottom pressure gauges

Vertical displacement time-series derived from five ocean-bottom pressure gauge campaigns are used for assessing the reliability of predictions in the offshore part of the domain of interest. Details of the dataset are available from Peltier et al. (2022). Due to logistical complexity of deploying these experimental instruments on the seafloor, no measurements are available before March 2019. Furthermore, since measurements were not telemetered, it was necessary to recover the instruments to retrieve the data, which required two cruises per campaign. Because the pattern and amplitude of offshore ground deformation was largely unknown at the time of deployment, the sensors were moved from campaign to campaign, in order to maximize the chances of detecting some deformation. Hence, these measurements may only provide temporally-discontinuous and spatially-scattered measurements of ground deformation. Due to the short duration of the campaigns (no more than 4 months), a single mean velocity is estimated per campaign and per station. Pressure measurements are detided and corrected from sea surface height variations using the oceanic model ECCO2. Velocity for each station is estimated using a least-squares regression, with uncertainty calculated by taking into account a temporal correlation of 10 days within the time series. A common mode, estimated from principal component analysis, is removed to visually improve the signal-to-noise, with no incidence on estimated rates and uncertainties. The time-series are represented in Figure S9. The estimated vertical velocities are reported in Table S6.
# 3 Deformation model

## 3.1 Mathematical framework of semi-kinematic frame

We adopt the deformation model standard proposed by Crook (2020) for performing coordinate epoch transformations. This standard has been implemented in the coordinate transformation open source software “PROJ” (Evenden, 2024).

In the Crook (2020) formulation, the displacement **d** at the point of geographical coordinates **x** (latitude and longitude, or easting and northing), at any epoch *t* (expressed in a reference time scale or relative to a reference epoch), and in any direction ( conventionally, ENU for East, North, Up), has to take the following form:

$$ \mathbf{d}(t, \mathbf{x}) = \sum_{i=1}^{p} f_i(t) \cdot \mathbf{g}_i(\mathbf{x}) \quad \text{or, equivalently,} \quad \begin{bmatrix} d^E(t, \mathbf{x}) \\ d^N(t, \mathbf{x}) \\ d^U(t, \mathbf{x}) \end{bmatrix} = \sum_{i=1}^{p} f_i(t) \cdot \begin{bmatrix} g_i^E(\mathbf{x}) \\ g_i^N(\mathbf{x}) \\ g_i^U(\mathbf{x}) \end{bmatrix} \quad (1) $$

where $f_i(t)$ are “time functions” (with *i* the source index), $g_i^j(\mathbf{x})$ are the three “spatial grids” describing each component (with $j \in \{E, N, U\}$ the component index) of the displacement field for each source. Boldface notation denotes a vector.

Equation 1 makes the implicit assumption of a separation between temporal and spatial variables. The “strength” of a given source *i* is described by the temporal function $f_i(t)$, which may correspond to the time evolution of slip for a tectonic fault, or pressure change for a magma reservoir. On the other hand, the “spatial signature” of the source is captured by the three spatial grids $g_i^E(\mathbf{x})$, $g_i^N(\mathbf{x})$ and $g_i^U(\mathbf{x})$ (one per coordinate). These spatial patterns, which may be envisioned as Green’s functions, describe the proportionality factor for converting unit slip (or unit pressure change) at the source into a displacement at point coordinate **x**. The grids $g_i^j(\mathbf{x})$ may be constructed by inverting for the geometrical parameters of the source of deformation using e.g. well-known analytic solutions for the elastic response of the Earth for a unit strength of the source, such as Okada (1985) for dislocations or Mogi (1958) for pressure sources. Summation in Equation 1 simply reflects the principle of superposition, where *p* is the number of sources. This mathematical framework is therefore consistent with the notion of linear elasticity of crustal deformation.

Based on prior understanding of the physical mechanisms at play in Mayotte during the 2018–2022 volcano-tectonic crisis, we consider that deformation can be explained
by one or several pressure sources at depth (i.e. $p \ge 1$). In the following, we take advantage of the separation of temporal and spatial variables in Equation 1 to implement a two-step procedure for the determination of the “best” deformation model. The method is schematized in Figure 2a. The first step (left column) consists in the determination of the temporal component, whereas the spatial component is recovered in a second step (right column). A detailed presentation of sub-steps of the method is provided in Supplement Text S1 (temporal component) and Supplement Text S2 (spatial component). Only the key aspects are summarized in the following sections, for brevity.

## 3.2 Principal component analysis

Following an approach similar to Kositsky and Avouac (2010), we first focus on determination of the optimal number of sources necessary to explain the variance of the dataset, as a compromise between misfit to the data and model complexity. To do so, we conduct a principal component analysis (PCA) on a subset of GNSS displacement time-series with near-complete temporal coverage (Supplement Text S1, Section S1.1). PCA relies on a truncated eigendecomposition to capture one or several common modes affecting a set of synchronous time-series (Tiampo et al., 2004). Under the assumption of a low correlation between temporal components of deformation, PCA can be applied for source separation or denoising.

PCA is applied to detrended time-series from the six GNSS stations that operated continuously over the period 2018-02-09 to 2021-05-01 (MTSA, GAMO, KAWE, MAYG, PORO, BDRL), combining ENU components into an 18-row data matrix. Due to the loss of station PORO after May 2021, we cannot extend the temporal window beyond that epoch without degrading the spatial resolution of the analysis. On the other hand, installation of MTSB, PMZI and KNKL between March and May 2019, i.e. well after the onset of the episode, implies that their respective time-series are incomplete, and cannot be included in the PCA.

We find that a single component suffices to decrease the data variance by more than 99.5% (Figure S10). This dominant component captures the progressive acceleration of deformation between May 2018 and December 2018, followed by a longer-lasting deceleration visible up until 2021 at least. Additional components do not lead to a further significative variance reduction which indicates that, within the precision and spatial cov-
erage of the GNSS data, deformation can be explained by a single pressure source that remained in the same place, and kept the same geometry throughout the 2018–2022 period. Hence, we may consider that displacement time-series at the surface are everywhere proportional to a single pressure time-history (i.e. *p* = 1).

<figure><img src="image_3.png" /><figcaption>Figure 2: (a) Procedure for determination of a deformation model of the 2018–2022 volcano-tectonic crisis in Mayotte. (b) Temporal component $f_{\theta}(t)$ sampled on 24-h time intervals. (c) Spatial component $g_{\lambda}(x)$ sampled on a regular latitude-longitude grid. Horizontal components are represented as a vector field and vertical component is represented as a colored and contoured background. Best-fitting source location is indicated by the blue circle.</figcaption></figure>
### 3.3 Trajectory model for the first eigenvector (temporal component)

Following the PCA analysis, we seek to find an analytic function for fitting the eigenvector of the first component of the PCA (hereafter noted “PCA1”). This step is essential to define a smooth trajectory model for all time-series that would eliminate (or at least mitigate) the residual noise in PCA1. We find that PCA1 can be closely matched by the so-called “double-exponential” time function. This time function has been proposed for interpreting geodetically-observed deformation over time spans of several years in other volcanic systems, such as Grímsvötn (Reverso et al., 2014), Laguna del Maule (Le Mével et al., 2016), Lazufre (Henderson et al., 2017), or Yellowstone (Delgado & Grandin, 2021) (see also Walwer et al., 2019).

The “double-exponential” temporal function describes the evolution of pressure within a deflating magma reservoir, with magma flowing viscously through a conduit into an external underpressured source behaving as a sink for the reservoir-conduit system. The main reservoir is initially at magmatic equilibrium before $t = t_0$, after which the external reservoir imposes a linear pressure drop at the end of the conduit between $t = t_0$ and $t = t_1$, followed by a constant underpressure for $t > t_1$. The temporal evolution of underpressure in the main reservoir is expressed $\delta P(t) = P \times f_\theta(t)$, where P is the final pressure difference of the reservoir (relative to its initial pressure) for $t \to +\infty$, and $f_\theta(t)$ is a normalized piecewise function defined as:

$$f_{\theta}(t) = \begin{cases} 0 & \text{for } t < t_0 \\ \frac{t-t_0}{\Delta t} - \frac{\tau_p}{\Delta t} \times \left(1 - \exp\left[-\frac{t-t_0}{\tau_p}\right]\right) & \text{for } t_0 \le t \le t_1 \\ 1 - \frac{\tau_p}{\Delta t} \times \left(\exp\left[-\frac{t-t_1}{\tau_p}\right] - \exp\left[-\frac{t-t_0}{\tau_p}\right]\right) & \text{for } t > t_1 \end{cases} \quad (2)$$

where $\Delta t = t_1 - t_0$ and $\tau_p$ is the characteristic time of the reservoir-conduit system. The function $f_\theta(t)$ is represented in Figure 2b and Figure S12 (see also Supplement Text S1, Section S1.2 for details). The non-linear parameters of Equation 2, gathered in the set $\theta = \{t_0, \Delta t, \tau_p\}$, are fitted on PCA1 by exploration of the parameter space and least-squares minimization of residuals (see Supplement Text S1, Section S1.3 for details).

The inversion yields: $\hat{t}_0 = 25$ May 2018, $\hat{\Delta}_t = 204$ days (or, equivalently, $\hat{t}_1 = 15$ December 2018) and $\hat{\tau}_p = 201$ days, where the “hat” symbol denotes optimal value of a parameter. Uncertainties on all three parameters is of the order of 10 days, but some trade-offs between parameters as visible in Figure S11.
3.4 Source model for PCA1 (spatial component)

Under the hypothesis that only a single source contributes to deformation at May-
otte (Section 3.2), the principle of linear elasticity implies that the predicted displace-
ment is everywhere proportional to δP (t), hence to f̂θ(t).

Using the trajectory model provided by the temporal function $f_{\hat{\theta}}(t)$ (Equation 2), the spatial pattern of displacements $g_{\lambda}(\mathbf{x})$ is recovered by inverting for parameters of a physical source model. The quantity $g_{\lambda}(\mathbf{x})$ describes the displacement for a unit change of the source underpressure. The non-linear parameters of the model (gathered in a set noted $\lambda$) include (i) the geometry of the source (shape, size) and (ii) its location (geographic coordinates, depth). Determination of the “best” source parameters, noted $\hat{\lambda}$, hence $g_{\hat{\lambda}}(\mathbf{x})$, is done in three substeps, as described in Supporting Text S2.

First, in order to complement the subset of six GNSS stations used in the PCA anal-
ysis (for which eigenvalues are already known), we project, for each InSAR pixel and each
of the nine available GNSS station, the corresponding displacement time-series onto the
temporal component $f_{\hat{\theta}}(t)$ from Equation 2. This is achieved by a weighted linear regres-
sion, as described in Supporting Text S2, Section S2.1. This strategy produces a coher-
ent displacement field, which can be interpreted as the best-fitting cumulative displace-
ment over the whole duration of the volcano-tectonic transient, regardless of the respec-
tive temporal coverage of input data. This step is essential to enable incorporation of
GNSS stations installed after the beginning of the crisis. It also provides an appropri-
ate scaling factor for the LOS displacement field derived from the Sentinel-1 and COSMO-
SkyMed InSAR analysis, which otherwise would cover only a fraction of the total du-
ration of the crisis. The resulting uncertainty on cumulative displacement, both for GNSS
and InSAR, is of the order of 10 mm. The projected displacement field is represented
in Figure S13 for the GNSS dataset, and Figures S14 and S15 for the InSAR datasets.

Then, using this internally-consistent InSAR and GNSS dataset, we aim at pre-
dicting the displacement pattern everywhere on the island. To do so, we estimate pa-
rameters of a Compound Dislocation Model (CDM) that best fits the displacement field
sampled at available observation sites (Nikkhoo et al., 2016) (Supporting Text S2, Sec-
tion S2.2). The CDM model is chosen for its ability to emulate the moment tensor of
a triaxial ellipsoidal pressure source of arbitrary shape, with improved computational
efficiency compared to previous numerical methods. A Bayesian inversion is computed
to determine the optimal geometry and location of the CDM that best explains the dis-
tribution of cumulative ground displacements (Supporting Text S2, Section S2.3). Af-
ter testing several source geometries (including a simple Mogi source), we find that the
best-fitting CDM source consists in a degenerate prolate spheroid (cigar-shaped) with
near-vertical orientation, steeply dipping toward the ESE (Figure 3a). Its centroid is lo-
cated at ~ 23.9 km depth, at geographic coordinates (−12.8094°N; 45.5148°E) (see Fig-
ure S16). The actual size of the source is unconstrained, as a source with a finite size
yields the same fit to the data as a point source (Supporting Text S2, Section S2.4). The
CDM source strength represents the cumulative underpressure over the whole crisis (i.e.
up to ~ late 2022), and can be converted into a cumulative volume change of -4.8 km³.
This value is broadly consistent with the bulk volume of the Fani Maoré volcano (–6.5 km³),
noting that the latter estimation is not corrected for vesicularity (Berthod et al., 2022).

Finally, the optimal parameters of the CDM are used in a forward computation of
the predicted ENU displacement over the island and offshore, sampled on a regular grid,
which materializes g<sub>λ</sub>(x). This spatial component, representing the cumulative 2018-
2022 displacement, is represented in Figure 2c.

As a side note, the inversion step arguably adds an extra complexity to the pro-
cedure. A simpler approach, such as polynomial interpolation/extrapolation, may also
perform reasonably well. However, in the specific case of Mayotte, where optimal pre-
cision is sought, the limited spatial sampling of available geodetic data does not allow
to capture spatial gradients of displacement at all scales. Applying an empirical fit on
available data poses the risk of introducing uncontrolled and systematic errors away from
the data points. Instead, a physical model (such as the pCDM used here) has the ca-
pability of predicting realistic gradients, including in areas where no data is available,
i.e. beyond the convex hull of the GNSS network, and in the interior of the island (Fig-
ure S17).

**3.5 PROJ implementation**

The deformation model described above is converted to a format ready for inges-
tion by users through PROJ software command line utilities. Here, we describe the for-
mat conversion procedure, whereas details of the implementation are provided in Sec-
tion Text S4.
In order to comply with the PROJ format, the temporal function of our model, $f_{\hat{\theta}}(t)$, has to be decomposed as a sum of a set of elementary functions (noted $f_i(t)$ in Equation 1). They must be selected from a dictionary that includes step, linear, piece-wise linear or exponential functions. The definition of $f_{\hat{\theta}}(t)$ (Equation 2) makes it straight-

<figure><img src="image_4.png" /><figcaption>Figure 3: (a) Observed and modeled cumulative displacement associated with the 2018–2022 volcano-tectonic transient for permanent GNSS stations. Inset shows a 3D representation of the best-fitting Compound Dislocation Model source. (b) Observed and modeled cumulative line-of-sight displacement from InSAR data for Sentinel-1 descending data. (c) Same as (b) for COSMO-SkyMed ascending InSAR data.</figcaption></figure>
forward to describe the full temporal evolution of deformation in Mayotte with only two piecewise exponential functions and one piecewise linear function. All three functions share the same spatial function $g_i(\mathbf{x})$ provided in a single file (geotiff format).

In our final implementation, the deformation $d(t, \mathbf{x})$ predicted by the model was conventionally modified to predict zero deformation at $t \to +\infty$ so that users working with newly acquired geodetic data do not have to apply the deformation model. The deformation predicted by the obtained model at epoch 2023.75 is lower than 0.03 mm on the three coordinate components for inland sites and can be considered negligible. Conversely, the deformation computed at the start of the crisis in May 2018 by means of the final model will be equal to the opposite of the cumulative deformation.

Note that the deformation model standardization has been recently extended to other grid format and time functions (Crook, 2024).

# 4 Assessment of posterior model uncertainty using a validation dataset

## 4.1 Absolute displacement uncertainty over 2018–2024

To quantify posterior uncertainty of the model, we compare the predicted displacement time-series against those deduced from a reanalysis of permanent GNSS data.

We first form model-corrected coordinate time series by removing the displacements predicted by the model to the reprocessed PPP position time series of the validation dataset. They are transformed from IGS14b to ITRF2020 beforehand. A station-specific parametric model that includes a linear trend, and periodic annual and semi-annual signals, is then fitted to each series. In addition, a time-correlated noise model is estimated to process realistic associated uncertainties. The parameters of a power-law plus variable white noise model are fitted using restricted maximum likelihood (Gobron et al., 2022). Estimated velocities and their associated 95% confidence ellipses are reported in Table S1.

The mean ITRF2020 model-corrected velocity computed over all corrected series is −12.98 ± 0.43, 17.88 ± 0.43 and 11.49 ± 0.27 mm/yr in X, Y and Z component respectively, which corresponds to 21.8 ± 0.3, 12.0 ± 0.2, 0.9 ± 0.5 mm/yr in the East, North and Up components. The maximum differences to that mean are 1.8±0.4, 1.7±0.2 and 3.1±0.6 mm/yr in East, North and Up for PMZI, MAYG and KNKL respectively. Figure 4a shows the residual time series of three permanent stations after remov-
ing residual trend and seasonal variations (all series are displayed in Figure S18). While residual series on the Up component are rather stable, MAYG residuals along North as well as KAWE, GAMO and MTSA in a smaller extent, show a transient anomaly during the first year of the crisis and an apparent velocity change, still subsisting in 2024 (this residual anomaly will be specifically analyzed in Section 4.3). Short-duration anomalies are also visible in the East components for a few stations as BDRL, KAWE, MAYG and PORO. Nevertheless, all time-series remain restricted within a band of ± 1 cm about the mean at any time for North (respectively, ± 2 cm for East and ± 3 cm for Up). These values can be interpreted as a maximum bound for absolute displacement error by the model.

## 4.2 Relative displacement uncertainty over 2018–2024

In order to evaluate the performance of the model to correct for internal deformation of the island induced by the 2018–2022 volcano-tectonic crisis, we further compute all possible baselines between the 9 permanent stations using the GNSS solutions of the validation dataset until June 2024 (Section 2.2.1). Figure 4b shows an example of baseline time-series and the associated deformation predicted by the model. Another example is provided in Figure S19. The mean baseline computed over the last common 300 daily values has been removed so that baseline changes are zero in the recent months. Due to our model conventions (Section 3.5), the model and observed baseline changes are also identical over this period.

Figures 4b and S19 indicate that the model captures well the observed deformation pattern and its magnitude. However, the East component of the KAWE-PORO baseline does not show a perfect fit during the first year of the deformation (as already noted in Section 4.1). To quantify this misfit, we evaluate the mean difference between the modeled and observed baselines at the first 50 days of observations. The maximum differences for the reprocessed GipsyX PPP time series are respectively 2.1, 1.1 and 1.7 cm for East, North and Up components at the baselines KAWE-PORO, MAYG-PORO, KAWE-MTSB. An alternative, independent regional GNSS solution (processed by IGN, not shown here) provides slightly larger maximum values back to May 2018 for the oldest records (2.5, 1.1 and 2.0 cm respectively) due to noisier data at some stations. This independent analysis confirms that the misfits cannot be attributed to GNSS processing errors introduced by the GipsyX PPP processing. This analysis indicates that relative horizontal
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Examples of GNSS coordinate time series (East, North, Up) computed by OVPF corrected for the deformation model and for station-specific linear trends and periodic annual and semi-annual signals. Coordinates have been shifted by 2 cm in horizontal and 4.5 cm in vertical for clarity. (b) Uncorrected KAWE - MTSA baseline changes along East, North and Up components for PPP (GipsyX) and regional double difference (Bernese) solutions. Changes predicted from the deformation model are represented as dashed lines.</figcaption></figure>
displacement between points on the island (i.e. island-wide strain) are corrected with an accuracy of ~1 cm (1-σ).

## 4.3 Residual post-2022 deformation

As observed in Figure 4a, the North component of MAYG model-corrected GNSS time-series shows an apparent velocity decrease, from 14.4±0.2 mm/yr before the crisis, down to 12.2 ± 0.5 mm/yr after May 2018. Given the absence of precursory signs prior to May 2018, we regard the existence of a silent geological process affecting the pre-2018 MAYG velocity as possible, yet entirely speculative, thus do not discuss this hypothesis any further. On the other hand, the velocity computed at MAYG after the crisis does not differ significantly from the mean velocity of the network (Table S1a). To assess this apparent velocity change, we fit station velocities only with data after June 2022. We exclude PORO and PMZI stations from this analysis since they cover less than 2 years of observations. The obtained post-2022 mean velocity for the remaining 7 stations is 21.5 ± 0.3, 12.5 ± 0.3 and -1.9 ± 1.0 mm/yr for East, North and Up components respectively (Table S1b). These values confirm the mean model-corrected 2014–2024 velocity estimated in Section 4.1 on the horizontal components, the vertical component being more uncertain.

Our deformation model predicts a velocity smaller than 0.1 mm/yr for the post-2022 period. Thus, this apparent post-2022 island-wide velocity change likely indicates that a residual deformation process, not taken into account in the model, still contributes to the observed displacement time-series in the 2022–2024 time interval (excluding a pre-crisis anomaly). Since the residual ~2 mm/yr velocity is mainly oriented toward the south, while our model predicts essentially no deformation on this component at MAYG, a modification of the temporal component of the model would not suffice to absorb this residual, and the spatial component of the model would also need to be updated. In other words, an additional source is needed to explain this discrepancy. This task is beyond the scope of this study. In the meantime, we recommend to use the post-2018 model-corrected velocity as representative of motion of the island over the coming years.
4.4 Reliability of offshore extrapolation of modeled deformation

Extending the analysis to offshore sites, estimated sea-bottom vertical velocities deduced from five pressure gauge campaigns are systematically compared to predicted velocities (Figure S20). We recall that these experimental ocean-bottom measurements only provide information on the vertical component of displacement, and are not designed to constrain horizontal movement. Since displacement rates are variable in time, a mean velocity must be extracted from the model over the time span of each campaign. For the first campaign (16 March 2019 – 5 May 2019), at site MOCE, which is situated almost exactly above the inferred source of deformation, the observed subsidence rate (–471 mm/yr) is consistent, at first order, with the model prediction (–627 mm/yr). On the other hand, at site MOSE, located near the submarine volcano Fani Maoré, observed vertical veloc- ity is estimated to –582 mm/yr, whereas the predicted subsidence rate is only –238 mm/yr. A similar underestimation of the subsidence rate by the model is noticed at sites MOAF (third campaign) and MOAG (fourth campaign), both situated close to Fani Maoré. This discrepancy could be explained by the presence of a local source of subsidence near the outlet of the magmatic system, possibly at shallow depth in the crust, which would not be resolved by geodetic observations collected onshore. If these specific three observa- tions points (MOSE, MOAF and MOAG) are excluded, the mean of the residual sub- sidence rates over the remaining 20 sites is reduced to –6 mm/yr, against a mean sub- sidence rate prior to model correction of –147 mmm/yr. On the other hand, the stan- dard deviation of residuals is reduced to 99 mm/yr, compared to a standard deviation of estimated subsidence rates of 148 mm/yr prior to correction. In other words, half of the variance between observations is explained by the model, and the residuals are es- sentially unbiased. We however recall that the deformation model in the offshore part of the domain is an extrapolation of data acquired onland. Nevertheless, this analysis may be conservatively interpreted as an indirect validation of experimental ocean-bottom pressure gauges, which do seem to capture realistic subsidence rates.

5 Quantification of local (intra-island) deformation processes

5.1 Residual horizontal displacement from GNSS campaigns (2002–2023)

The cumulative 2018–2022 displacement field predicted from the deformation model
can be used for bringing into a common reference frame campaign geodetic data acquired
before and after the crisis. Taking into account the absolute and relative uncertainty of
the model (see Section 4), residual displacements spanning the crisis may be interpreted
to identify processes of local deformation. This approach can also be used to retrospec-
tively assess the stability of sites equipped with permanent GNSS stations.

We first compute displacement vectors deduced from two GNSS campaigns in 2003
and 2023, respectively, covering 16 common sites (Section 2.2.2), in a plate-fixed refer-
ence frame, namely ITRF2020 at epoch 2023.75. Data from the 2003 campaigns were
processed in ITRF2000 at epoch 2004.0, while those of the 2023 campaign were expressed
in IGS20 (compatible with ITRF2020) at epoch 2023.75. In order to isolate the defor-
mation related to the volcanic crisis, 2003 campaign coordinates were first transformed
to ITRF2020 at epoch 2004.0 using the official transformation parameters available from
the ITRF web site (IGN, 2024). However, the mean motion of the island as observed in
ITRF2020 has to be corrected for by propagating the coordinates to the epoch 2023.75,
before computing the displacements in a plate-fixed frame. We choose not to use any Eu-
ler pole model due to the uncertain tectonic context of Mayotte. As a consequence, we
first use the model-corrected mean velocity provided in Section 4.1 for this purpose. While
it relies on the model we wish to evaluate, it is the only available estimate based on the
full records of 9 permanent stations. When using it to propagate 2004.0 coordinates to
2023.75, the observed cumulative displacement at the 16 stations differs from our mod-
eled displacement by about 5.0 ± 0.8 cm (1-σ) in the North component. This corresponds
to a velocity difference of about 2.5 mm/yr over 20 years. This discrepancy is coherent
with the result obtained in Section 4.3, where MAYG model-corrected series show an ap-
parent velocity change of 2.2±0.5 mm/yr after the crisis. Alternatively, using pre-2018
velocity of MAYG only to propagate 2003 campaign coordinates improves the fit between
campaign data and the model, with the mean North difference reducing to 0.4 ± 0.8 cm.
However, we previously noted the apparent velocity change at MAYG in GNSS model-
corrected GNSS series. Instead, we use pre-crisis velocity of MAYG and post-crisis model-
corrected mean velocity of all stations to propagate the 2004.0 GPS campaign coordi-
nates to epoch 2023.75 and compute the final displacement field from campaign data.
The mean difference between the modeled cumulated displacement and GPS campaign
displacement increases to 1.6±0.8 cm but this approach is more reliable since it relies
on all available GNSS permanent station data in Mayotte.
The formal errors of the 2004–2023 displacement are derived based on estimated standard deviation as output of the GNSS data adjustment and variance propagation. However, the uncertainty of the 2003 campaign coordinates was likely too optimistic. We therefore rescaled formal errors to reach the same median level of the recent campaign. As the alignment of 2004.0 coordinates to ITRF2000 showed RMS residuals of 9 mm in East and North component and 5 mm in Up component at reference sites, their uncertainties were modified accordingly. The formal errors of the ITRF2000 to ITRF2020 transformation and the coordinate time-propagation were also accounted for, based on their respective standard deviation. The final median standard deviation of the displacements are 1.23, 1.08 and 1.26 cm in East, North and Up respectively.

Figure 5b shows the difference between observed and modeled cumulative displacements at the 16 sites for both horizontal and vertical components. 95% confidence intervals are superimposed. For most of the sites, the differences are not significant, which validates our model implementation. Furthermore, insignificant residual displacement is found at the two permanent GNSS sites colocated with campaign GNSS (PMZI to the east of Mayotte, and GAMO to the west). This confirms that these sites can be considered as stable, at least at local scale. Large horizontal and vertical displacements are however found at sites 12CB and 03CB (respectively to the north-west and the south-east of Mayotte), suggesting the existence of local deformation processes, in addition to the volcanic transient. This issue will be further discussed in the following sections using levelling and InSAR data.

## 5.2 Evidence for local vertical deformation from levelling campaigns (2006–2023)

Data from two levelling campaigns provides an estimate of the cumulative height change over the period 2006-2023. We compare MSL heights from levelling (no gravity correction applied) with modeled height changes for the period 2018–2023. The raw comparison of MSL height differences from levelling and modeled height changes shows an overall good fit and similar patterns (Figure S5a and S5b), although a substantial residual tilt of 1.30±0.10 and 1.54±0.06 mm/km remains along East and North directions (Figure S5c). This tilt translates into a North-South difference of about 55 mm over a distance of 40 km, which is well above the uncertainty of the cumulative vertical displacement from the model (see Section 4.1).
<figure><img src="image_6.png" /><figcaption>Figure 5: (a) Estimated 2018–2024 velocities in ITRF2020 after correction by the deformation model, with 95% confidence ellipses (grey, see inset for full-length vectors, and Table S1). Same velocities, additionally corrected for: Somalia plate motion in ITRF2020 (orange), Lwandle plate motion in ITRF2020 (light blue), weighted averaged mean residual velocity (green). (b) Model-corrected displacements from the two GNSS campaigns of 2003 and 2023 in a “plate-fixed” frame (see text). Horizontal differences are shown in blue and vertical differences in red. (c) Height difference from the two levelling campaigns of 2006 and 2023 (colored disks) corrected from the model and an empirical residual North-South tilt. Residual vertical differences at GNSS campaign sites are shown as red vectors, as in panel (b). (d) Residual line-of-sight velocity from model-corrected 2015–2024 ALOS-2 InSAR time-series, with velocity calculated from averaging in a spatial neighborhood of 10×10 pixels (or 300×300 meters). Note the double-sided colorbar, where LOS velocity is converted to pseudo-vertical displacement cumulated over 17.5 years (time-span between levelling campaigns).</figcaption></figure>
The reason for this tilt is still not fully understood. At first, we evaluated the precision of the levelling height difference at the level of 12 mm for 40 km-distant sites. Indeed, an a posteriori evaluation of the precision of the 2023 campaign from the least squares adjustment of data gives a precision of 0.96 mm/√km. We applied the same uncertainty to the first campaign results since it was surveyed similarly. While the evaluation of the levelling error is probably optimistic, the observed tilt is nearly one order of magnitude larger than the estimated error. A geoid tilt might explain this finding, however, the lack of reliable gravity observations in the 2006 campaign prevents evaluating it with in-situ gravity data. A simple dipole model of the mass transfer caused by the creation of the Fani Maoré Volcano (6-7 km³) can only explain millimeter-level geoid changes for inland sites. Finally, the tilt might be related to the lack of gravity correction. The 2023 gravity measurements at the 9 sites show a variability of about 30 mGal. This would cause the geometric heights to differ from gravity-compensated height by about 3 mm per 100 m difference. Yet, the maximum MSL height of levelling benchmarks is 163 m (except for a few inland sites in the South of Grande-Terre that reach 350 m), which remains insufficient to explain the magnitude of discrepancies. It is difficult to take the analysis any further without a complete model of gravity variations at network points. As a consequence, we decide to remove the tilt from our levelling height change observations to concentrate on small-scale features.

Tilt-corrected height change differences between levelling observations and modelled heights are shown in Figure 5c. The RMS of the difference is 9.7 mm. Levelling observations confirm the large subsidence at the two sites 12CB and 03CA/CB already measured by campaign GNSS (Section 5.1). On one hand, 12CB marker is located on a low wall on a cliff along the coast and is the only benchmark in the area affected by significant motion. We interpreted this anomaly as reflecting instability of the benchmark site at a scale of a few tens of meters. On the other hand, 03CA and 03CB are markers located on a pavement and a rock near the coast, and are separated by a distance of 330 m. They show a consistent subsidence of ~ 45 mm with 4 surrounding markers highlighting a local motion of the area of Hamouro city in the south-east of Grande-Terre. The barely significant uplift in GNSS data at 10AA in the North of Grande-Terre is also confirmed while 04BB subsidence coming from GNSS data is visible with a much smaller magnitude in levelling data as well as the uplift at 09BA in the South of Grande-Terre.
As a conclusion, except for the subsidence in the immediate surrounding area of Hamouro city, no other medium-scale pattern of deformation can be identified.

## 5.3 Confirmation of local deformation in the Hamouro area from ALOS-2 InSAR (2015–2024)

To further investigate the existence of small- or medium-scale deformation patterns, we map island-wide deformation using ALOS-2 SAR interferometry spanning the interval 2015–2024 (Figures S7 and S8). By tracking pixel motion (with respect to an arbitrary reference pixel situated near GAMO), the signature of the volcano-tectonic 2018–2022 transient is clearly detected in the time-series (see Figure S21a). Subtraction of the deformation predicted by the model allows for extracting a residual LOS velocity field over the whole island, after fitting a linear trajectory on residual time-series for every pixel in the map (Figure S21b).

In Figure 5d, the resulting ALOS-2 velocity map is compared to the 2006–2023 levelling difference, converting LOS velocity to pseudo-vertical cumulative displacement by (i) dividing by the cosine of the incidence angle and (ii) multiplying by 17.5 years (2006–2023 time difference). This approach implicitly ignores horizontal motion and assumes that any displacement unrelated to the crisis proceeds steadily between 2015 and 2024. To mitigate potential bias due to vegetation or soil moisture (De Zan et al., 2015), we concentrate our comparison on pixels classified as “Artificialized spaces” in the land cover map of Mayotte (Dupuy et al., 2024, see Figure S22 and Section Text S3 for details).

The site-to-site comparison between InSAR and levelling data confirms the presence of an area of anomalous subsidence at Hamouro, to the south-east of the island. Apparent subsidence occurs at levelling sites NN-308a, NN-308-I and NN-309-I (colocated with GNSS sites 03CA and 03CB as noted in Sections 5.1 and 5.2), coinciding with InSAR pixels affected by motion away from the satellite (Figure 5d). There, the subsidence rate estimated by InSAR, with respect to surrounding areas, reaches ~ 2.5 mm/yr (just above the InSAR noise level), translating into a vertical displacement of ~ 50 mm over 17.5 years. The subsidence rate at Hamouro is consistent between InSAR and levelling, after taking into account the different time spans of these two independent datasets, indicating a steady process unrelated with the 2018–2022 crisis. Attempts to identify other
potentially unstable sites throughout the island were unfruitful, due to the poor signal-to-noise of InSAR data.

# 6 Insights into the long-term motion of the island in the ITRF

The deformation model introduced in this work allows for removing most of the non-linear deformation affecting permanent GNSS time series and can be used to estimate model-corrected velocities from daily positions. The correction is indispensable because the 2018–2022 transient overshadows the slow translation of the island due to plate tectonics. As Mayotte is one of the few continuous GNSS sites available in this region, it is crucial for geodynamic studies to determine this long-term velocity.

Figure 5a shows the estimated ITRF2020 velocities relying on all available GNSS data at Mayotte. These velocities are compared with predictions coming from Euler pole models by Altamimi, Métivier, et al. (2023) for Nubia and Somalia plates. Euler pole models in ITRF2020 are also derived for Lwandle and Rovuma plates using relative Euler pole models from Stamps et al. (2021), relying on the previous Nubia or Somalia absolute poles. Differences between our velocities and the Somalia and Lwandle plate motion model predictions (Somalia absolute + Lwandle relative poles) are plotted on Figure 5a. The Euler pole model that best fits the model-corrected velocities is the one of Lwandle plate, with residual mean velocities of −0.33 and −2.26 mm/yr in East and North component respectively compared to −0.37 and −2.75 mm/yr for Somalia. The residual North velocity decreases by −0.83 mm/yr when not correcting for the ITRF2020 Origin Rate Bias (ORB) which remains insufficient to explain the discrepancy. This bias cannot be explained by an alignment error to ITRF2020. Indeed, our MAYG velocity fitted before May 2018 agrees with the ITRF2020 published velocity for that station at a level better than 0.2mm/yr for each component. We note that MAYG motion shows a good individual fit with the Lwandle plate model (Figure 5a). The MAYG velocity estimated with the records prior to May 2018, namely 21.7 ± 0.2 and 14.4 ± 0.2 mm/yr along East and North directions respectively, even fits better the Lwandle pole model prediction, which is 22.1 and 14.2 mm/yr respectively (22.2 and 14.7 mm/yr for Somalia). However, the fit is not as good for other sites, and yields a systematic difference on the North component. Relying on the deformation model and the full set of permanent stations, we therefore recommend to take our mean velocity as the best estimate of May-
otte displacement for the period 2018–2024 (Table S1), instead of the single-station MAYG rate.

Finally, for the sake of completeness, we here report the model-corrected relative sea-level rise velocity from tide gauge DZAOU (colocated with MAYG GNSS station), which we estimate to +3.4±0.1 mm/yr for the period 1985–2024 (Figure S24). Incorporating the +1.3±0.3 mm/yr ITRF-corrected uplift rate at MAYG, the absolute sea-level rise at DZAOU increases to +4.7 ± 0.3 mm/yr.

# 7 Discussion

## 7.1 Limitations of the deformation model and possible improvements

The seismic swarm that shook the island of Mayotte in May 2018 marked the onset of an unprecedented episode of volcanic and tectonic unrest. A number of studies have already attempted to elucidate the processes at play during this episode using geodetic data, with variable degrees of sophistication. Among the published models, some only considered the spatial component of deformation, individual GNSS time-series were analyzed by fitting a piecewise-linear temporal function in separate time intervals (Lemoine et al., 2020; Cesca et al., 2020; Peltier et al., 2022). On the other hand, Mittal et al. (2022) proposed a porous flow model capable of reproducing the observed time history of displacement, but did not attempt to reproduce the spatial distribution of deformation. Finally, several conceptual models were recently proposed to interpret the distribution of seismic events, but did not incorporate geodetic constraints in the analysis (Lavayssière et al., 2022; Jacques et al., 2024; de Sagazan et al., 2024).

In contrast, our work aimed at bringing all geodetic data into a common reference by describing deformation in space and time using a physical model, while keeping this model as simple as possible. We showed that the deformation pattern can be explained, at first order, by a pressure source draining according to a “double-exponential” time function. This temporal component has been interpreted as reflecting the deflation of a “visible” reservoir, draining into a second “stealth” reservoir, which imposes a piecewise linear underpressure forcing in time and is possibly connected to the surface. On the other hand, the spatial component has been reproduced by a vertically-elongated ellipsoidal cavity, dipping steeply to the ESE, at ~24 km depth and ~20 km off the east of Mayotte. We found that it was possible to analyze the spatial and temporal compo-
nents of deformation separately because they appear to be decoupled. However, in or-
der to shape a coherent picture of the architecture and internal plumbing of the magma
system offshore Mayotte, these temporal and spatial characteristics need to be consid-
ered from a common perspective.

First, regarding the temporal component, we note that the “double-exponential pat-
tern” (Le Mével et al., 2016) is ubiquitous in nature, and can be produced by a broad
variety of processes. Mathematically, our “double-exponential” function represents the
response of a first-order linear system to an input consisting of a finite-duration ramp.
However, it is possible to produce a nearly identical trajectory using a second-order sys-
tem consisting of two reservoirs, combined with a step-like forcing term instead of a ramp
(as shown in Figure S25). Higher-order systems contain an even greater number of in-
ternal degrees of freedom capable of fitting the data equally well. In fact, Mittal et al.
(2022) explain the temporal pattern observed at Mayotte by introducing a draining mush
model, which is essentially equivalent to an infinite-order system (neglecting conduit clo-
sure by freezing and other non-linear processes). Thus, the “double-exponential” tem-
poral evolution of deformation may be compatible with the drainage of a mush.

On the other hand, the spatial component of deformation is well fitted by a nearly-
vertical prolate spheroid, which, at first sight, is harder to interpret. Furthermore, the
source inferred from this analysis is not localized below Fani Maoré, but is instead sit-
uated closer to the island, within a few kilometers of the upper part of a deep (25–45 km)
steeply-eastward-dipping seismically-active structure imaged by Lavayssière et al. (2022)
(Figure 3a). This seismicity has been interpreted as reflecting the rupture of the cone-
shaped roof of a transcrustal reservoir, in response to drainage of a mushy reservoir sit-
uated at even greater depth (40–50 km) (Lavayssière et al., 2022; de Sagazan et al., 2024;
Jacques et al., 2024). Jacques et al. (2024) proposed that magma is transferred from the
lower to the upper reservoir, flowing along the edges of the conical roof. If magma does
move through such a complex network of interconnected pockets, our simple elastic model
is unlikely to be a faithful representation of the actual architecture and internal plumb-
ing of the magma system. In addition, given the inferred depth of the source (> 20 km)
and the time-scale of pressure changes affecting it (several years), host-rock surround-
ing the reservoir certainly does not behaves elastically (e.g. Gregg et al., 2018; Segall, 2019).
Notwithstanding, two independent studies may give a hint for the possible consistency of our results with the “transcrustal mush hypothesis”. First, using discrete element modeling, Holohan et al. (2017) showed that the presence of host-rock fracturing above a deep deflating reservoir may explain (i) the shallow inward-dipping seismicity, (ii) the bias of the centroid of deformation toward the surface and (iii) apparent distortions of the inferred pressure source, if incorrectly modeled as an ellipsoidal cavity embedded in an elastic medium. Second, using finite element modeling, Mullet and Segall (2022) attempted to predict the deformation field induced by a vertically-extended transcrustal mush, capped by a melt lens subject to a rapid pressure drop. The authors showed that “for moderately high levels of mush melt fraction the transcrustal mush model generates surface displacements that have (...) an effective geometry of a vertically elongate prolate spheroid” (Mullet & Segall, 2022).

More work is needed to reconcile available geophysical observations into a physically-consistent model, ideally also incorporating petrological constraints (Verdurme et al., 2024). We however posit that, once combined with our own traditional (elastic) analysis of deformation, the studies of Mullet and Segall (2022) and Holohan et al. (2017) for the spatial component of deformation, Mittal et al. (2022) for the temporal component of deformation, and Lavayssière et al. (2022) for seismicity, are strongly suggestive of the drainage of a transcrustal magma storage system offshore Mayotte starting from 2018.

## 7.2 Relevance of the semi-kinematic frame in the context of a transient deformation event

When a whole territory is deforming, such as in response to post-seismic relaxation or volcanic unrest, reference geodetic benchmarks may be affected by the deformation, and can no longer be considered as stable. Failing to recognize and compensate for this perturbation results in an incorrect determination of coordinates derived from survey networks used in most land planning projects. The temporal and spatial extent of these errors, as well as their magnitude, scale with the size of the causative geologic process. For example, post-seismic relaxation following large (M>8.5) earthquakes can generate > 1 mm/yr velocity anomalies at distances exceeding 1000 km and over durations of nearly 10 years (e.g. Liu et al., 2023; Klein et al., 2016, for 2010 Maule and 2011 Tohoku earthquakes, respectively). Correction of these long-lived transients is often dealt with by implementing *ad hoc* strategies (e.g. Tarrío et al., 2024; Ansari & Park, 2022; Gómez et al., 2016). Volcanic unrest
on the other hand may also induce long-lasting deformation transients, but with a comparatively smaller spatial footprint (e.g. Le Mével et al., 2021; Astort et al., 2024, at Laguna del Maule and Campesino). Thus, it is not commonly considered as an issue for national- or regional-level geodetic network stability.

The case of Mayotte, although specific, illustrates the necessity of establishing a formal procedure, by means of definition of a semi-kinematic frame, to correct for transient deformation induced by volcanic or tectonic activity. In Mayotte, coordinates of points deduced from surveys must be provided, by law, in a certified local reference frame (République Française, 2019). Prior to 2018, the local frame (RGM04) was static, which means that the published coordinates of the reference network were constant and that it was not required from surveyors to associate an epoch to their published coordinates. In practice, GNSS surveys in Mayotte are mostly carried out using standard surveying or differential GNSS (static survey, RTK and N-RTK), which provides coordinates directly in RGM04, assuming that the coordinates of one or several public GNSS permanent stations are known. Conversion between RGM04 and ITRF coordinates is then achieved, if necessary, by means of a 14-parameter transformation (IGN, 2011). Thus, prior to 2018, knowledge of the relative position with respect to MAYG, or any other RGM04 reference point, was enough to obtain an absolute position. Conversely, GNSS surveys using PPP, which have the advantage of directly providing coordinates in the ITRF, could be transformed to RGM04 using the inverse transformation.

Onset of the Mayotte crisis disrupted this traditional approach. Accurate compensation of non-stationary displacements now requires specifying the epoch of measurements. Internal deformation of the island has made all coordinates determined using differential positioning during the crisis biased in RGM04, since each baseline had its own time variation. Furthermore, as RGM04 static coordinates are increasingly biased due to island-wide motion, the 14-parameters transformation has become outdated. Thus, producing a deformation model for the 2018–2022 time interval became essential to georeference surveys in a local reference frame consistent with the shape of the island at the epoch of the measurement, at the required level of accuracy.

Following the 2023 GNSS campaign, IGN has derived an updated “plate fixed” frame for Mayotte, named RGM23 (IGN, 2023). Coordinates (referring here to ITRF2020) are published at a fixed epoch, here 2023.75 expressed in decimal year. The deformation
model presented in this paper has been designed to be used with RGM23 as it predicts no deformation at this reference epoch (Section 3.5). Thus, the model is not required for dealing with any measurement collected after 2023.75, but can be used to predict any RGM23 coordinates back to the past, reconstructing the shape of the island at any required epoch. Such coordinate time propagation is relevant to compute the coordinates of GNSS base stations at the epoch of GNSS surveys. Thanks to these updated base coordinates, collected GNSS data can be rigorously recomputed using differential GNSS at this epoch. To comply with the semi-kinematic frame approach, obtained coordinates have then to be transformed back to the reference epoch 2023.75 using the deformation model.

We also propose a transformation from ITRF2020 (at any epoch) to RGM23 at epoch 2023.75, by combination of a time-dependent translation rate and our deformation model. The translation rate is composed of a different vector before and after the crisis. We recall here the translation rate to be used with respect to our transformation convention (from ITRF2020 to RGM23): $12.98 \pm 0.43, -17.88 \pm 0.43$ and $-11.49 \pm 0.27$ mm/yr after May 2018; $13.5 \pm 0.4, -17.2 \pm 0.5, -14.1 \pm 0.2$ mm/yr before May 2018 for X, Y and Z components respectively. The maximum 3D residual velocity among all available GNSS permanent stations w.r.t. the post-2018 mean velocity is 3.2 mm/yr at PMZI, which represents a 3.2 cm error over 10 years. By computing the position time series in RGM23 for all permanent stations and computing the deviation to the mean coordinate of the last available year, we evaluated the transformation uncertainty to be at the level of 2 cm (95% confidence) for each individual component X, Y, Z. The reason for this time-varying translation is the presence of a residual post-2022 deformation anomaly that our deformation model does not attempt to take into account (Section 4.3). In the future, after data has been collected over a longer observation window, further work will be needed to assess the extent of this residual deformation, which will possibly require an update of the deformation model.

As a final comment, our thorough study of MAYG position time series leads us to propose an improvement of the kinematic model of the MAYG station in ITRF2020. Indeed, in the latest ITRF2020 update (namely, ITRF2020-u2023), no less than 6 position and velocity discontinuities have been introduced to model its motion by means of piece-wise linear functions (Altamimi, Rebischung, et al., 2023; IGN, 2024). To avoid such an over-segmentation of the MAYG coordinate time series and to simplify its kinematic model,
our deformation model could be used for correcting/complementing ITRF coordinates
for MAYG in addition to introducing a velocity discontinuity in May 2018. However, in
the current ITRF2020 kinematic model, no individual forward model is currently used,
although logarithmic and exponential functions have been widely adopted to model post-
seismic motions since ITRF2014. Our alternative recommendation would be to imple-
ment the mathematical functions introduced in the temporal component of our defor-
mation model, namely piece-wise linear model and two exponential functions, to fit a station-
specific kinematic model for MAYG, as currently done for sites affected by post-seismic
deformation.

Conclusion

We have developed a deformation model for the 2018–2022 Mayotte crisis, decom-
posed into a temporal component (modeled as a “double exponential” time function) and
a spatial component (modeled as surface deformation from a buried point compound dis-
location model). The model relies on a number of simplifying assumptions, especially
elasticity, and therefore does not reflect a faithful representation of the actual architec-
ture and internal plumbing of the magmatic system offshore Mayotte, which likely con-
sists in a transcrustal magma storage system. This deformation model can however be
used for correcting geodetic data from the contribution of the transient volcano-tectonic
deformation. After correction by the deformation model, independent pre- and post-crisis
campaign geodetic data, as well as InSAR data, indicate the existence of at least one ad-
ditional source of local subsidence in the Hamouro area, unrelated to the volcano-tectonic
crisis, reaching ~ 2.5 mm/yr over a region a few kilometers across. As a consequence,
our deformation model does not meet accuracy requirements for geodetic applications
in this specific location. Moreover, as our analysis of experimental pressure gauge data
shows that only half of their variance is explained by our model, the model is only rec-
ommended for onshore applications.

By construction, the deformation model does not include a long-term trend. As a
consequence, it can be used to propagate coordinates from an epoch to another in a “plate
fixed” frame at Mayotte. We provide a transformation between ITRF2020 and the new
frame RGM23 for Mayotte, valid at any epoch with a precision of 2 cm (2-σ) for each
component X, Y, Z. This transformation is useful for GNSS PPP users who process co-
ordinates in ITRF. Furthermore, the model makes it possible for local surveyors in May-
otte who collected GNSS data from May 2018 to 2022 to reanalyze their data and ex-
press them in RGM23.

As our model predicts negligible deformation from 2023.75 (RGM23 epoch, which
corresponds approximately to October 1, 2023), it can be ignored for any newly acquired
data. However, analysis of post-2022 permanent GNSS time-series indicates an appar-
ent velocity change w.r.t. single-station (MAYG) pre-2018 velocity, suggesting that some
residual deformation may still subsist at the time of writing. This is in keeping with the
fact that, at the time of writing, detected seismicity remains above the pre-crisis base-
line (REVOSIMA, 2024a). As more data is collected, and especially in the event of an
upturn in activity, the deformation model could be updated.
# Open Research Section

## Data

GNSS data used on this paper have been acquired and processed in the framework of the REVOSIMA (REVOSIMA, 2024b), and are available upon approval by the REVOSIMA Scientific Committee. Facilities of the European Space Agency (ESA) were used for access to Sentinel-1 products (ESA Copernicus, 2024b). ALOS-2 imagery was acquired and made available by JAXA via a cooperation agreement with ESA. COSMO-SkyMed and Sentinel-1 data were acquired and processed by the TRE-ALTAMIRA company, under contract with CNES. Raw tide gauge data is made available by Shom (2024), and is post-processed and validated by SONEL (2024). Ground displacement data deduced from differencing between GNSS and levelling campaigns is from IGN-SGM, and is available from Grandin and Collilieux (2025) at https://doi.org/10.57932/688c6c06-d553-4000-b80c-f95eb92a3b4b.

## Softwares

GNSS time-series were computed using Gipsy-X (Bertiger et al., 2020) on the We-bObs instance of the OVPF (Beauducel et al., 2020; Beauducel, 2024). ALOS-2 data was processed using the codes developed by Curen Liang and Eric Fielding (Liang et al., 2018), and implemented in the ISCE software (Rosen et al., 2012; ESA Copernicus, 2024a). We used the PCA package by Marianne Métois (Métois, 2022), based on the EMPCA by Bailey (2012a, 2016), for PCA analysis of GNSS time-series. We used the CSI package for In-SAR data subsampling and covariance estimation (Jolivet et al., 2012, 2014). The CDM code (Nikkhoo, 2016; Nikkhoo et al., 2016), vectorized by François Beauducel and translated in Python by Tara Shreve, was used for forward calculation of ground displacement. The Python library pyTRF developed by Paul Rebischung was used to estimate the GNSS permanent station velocities using a noise model fitted by restricted maximum likelihood.

## PROJ implementation

The deformation model files to be used with the PROJ library are available from Collilieux (2025) at https://doi.org/10.57932/41cf41f5-5a0c-47d0-b88d-a774e06e9012.
Acknowledgments

Since June 2019, Mayotte seismo-volcanic activity is monitored by REVOSIMA (Réseau de surveillance volcanologique et sismologique de Mayotte) (REVOSIMA, 2024b) 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).

Financial support from CNES to cover the costs of the SqueeSAR processing by the TRE-Altamira company is acknowledged. We thank Catherine Proy (CNES) for facilitating access to COSMO-SkyMed data, and Fifamè Koudogbo and Anne Urdiroz (TRE-Altamira) for managing the SqueeSAR processing. We also thank Pierre Potin (ESA) for approving the request to conduct additional Sentinel-1 acquisitions after the onset of the crisis in 2018. Finally, we thank Julia Kubanek (ESA), Pierric Ferrier (CNES), Osamu Ochiai and Shin-ichi Sobue (JAXA) for devoting time and resources that allowed for densifying acquisitions by ALOS-2 since 2018 over Mayotte.

We thank Valérie Ballu for sharing ocean-bottom pressure measurements and Pierre Sakic for GNSS data processing at IPGP.

We thank the teams of the Service de Géodésie et de Métrologie (SGM) of IGN for sharing data quality information following the 2004 and 2023 field campaign as well as Paul Rebischung for guiding us through the use of the pyTRF library.

Francisco Delgado, Manon Dalaison and Valérie Cayol are thanked for pointing to pertinent bibliographic references. We thank Federico Galetto and an anonymous reviewer for their constructive suggestions.

Digital bathymetry from Shom (2016) was used as background map for some figures.
References

Altamimi, Z., & Collilieux, X. (2024). *EUREF technical note 1: relationship and transformation between the international and the European terrestrial reference systems*. (Release: March 4, 2024)

Altamimi, Z., Métivier, L., Rebischung, P., Collilieux, X., Chanard, K., & Barnéoud, J. (2023). ITRF2020 plate motion model. *Geophysical Research Letters*, 50(24), e2023GL106373.

Altamimi, Z., Rebischung, P., Collilieux, X., Métivier, L., & Chanard, K. (2023). ITRF2020: an augmented reference frame refining the modeling of nonlinear station motions. *Journal of Geodesy*, 97(5), 47.

Amoruso, A., & Crescentini, L. (2009). Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data. *Journal of Geophysical Research: Solid Earth*, 114 (B2).

Ansari, K., & Park, K.-D. (2022). Geodetic analysis inside the south korean peninsula and impact of the 2011 tohoku-oki (to) earthquake. *Acta Geologica Sinica-English Edition*, 96(2), 631–647.

Argus, D., & Drummond, R. (2015). Space geodesy constrains ice age terminal deglaciation: The global ICE-6G_C (VM5a) model. *Journal of Geophysical Research: Solid Earth*, 120(1), 450–487.

Argus, D. F., Gordon, R. G., Heflin, M. B., Ma, C., Eanes, R. J., Willis, P., ...
Owen, S. E. (2010). The angular velocities of the plates and the velocity of
earth's centre from space geodesy. *Geophysical Journal International*, 180(3),
913–960.

Astort, A., Trasatti, E., Caricchi, L., Polcari, M., De Martino, P., Acocella, V., & Di Vito, M. A. (2024). Tracking the 2007–2023 magma-driven unrest at campi flegrei caldera (italy). *Communications Earth & Environment*, 5(1), 506.

Aubry, A., & Tempere, A. (2022). Apport de l'imagerie satellitaire pour l'identification et la cartographie des habitats littoraux de mayotte soumis à une subsidence rapide. *Revue Paralia*.

Bailey, S. (2012). *EMPCA: Weighted Expectation Maximization Principal Component Analysis [Software]*. https://github.com/sbailey/empca. (Accessed: 2024-04-14)

Bailey, S. (2016). Weighted EMPCA: Weighted Expectation Maximization Principal
Component Analysis. *Astrophysics Source Code Library*, ascl-1609.

Beauducel, F. (2024). *WebObs: an integrated web-based system for data monitoring and networks management [Software]*. https://ipgp.github.io/webobs/. (Accessed: 2024-04-19)

Beauducel, F., Lafon, D., Béguin, X., Saurel, J.-M., Bosson, A., Mallarino, D., ... others (2020). *Webobs: The volcano observatories missing link between research and real-time monitoring*. *Frontiers in Earth Science*, 8, 48.

Berthod, C., Komorowski, J.-C., Gurioli, L., Médard, E., Bachlery, P., Besson, P., ... others (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 (S2), 1-29.

Bertiger, W., Bar-Sever, Y., Dorsey, A., Haines, B., Harvey, N., Hemberger, D., ... others (2020). Gipsyx/rtgx, a new tool set for space geodetic operations and research. *Advances in space research*, 66(3), 469-489.

Blick, G., Crook, C., Grant, D., & Beavan, J. (2005). Implementation of a semi-
dynamic datum for new zealand.
In *A window on the future of geodesy:
Proceedings of the international association of geodesy iag general assembly
sapporo, japan june 30-july 11, 2003* (pp. 38-43).

Cesca, S., Letort, J., Razafindrakoto, H. N., Heimann, S., Rivalta, E., Isken, M. P.,
... others (2020). Drainage of a deep magma reservoir near Mayotte inferred
from seismicity and deformation. *Nature Geoscience*, 13(1), 87–93.

Charroux, S., Jeanson, M., Morisseau, S., & Pennober, G. (2024). Seasonal mor-
phodynamics of fringing reef pocket beaches and responses to rapid vertical
tectonic movements. *Journal of Coastal Research*, 40(1), 31–50.

Chen, C. W., & Zebker, H. A. (2002). Phase unwrapping for large sar interferograms: Statistical segmentation and generalized network models. *IEEE Transactions on Geoscience and Remote Sensing*, 40(8), 1709–1719.

Collilieux, X. (2025). *PROJ implementation of the Mayotte deformation model* (*2018-2022*). Retrieved from https://doi.org/10.57932/41cf41f5-5a0c-47d0-b88d-a774e06e9012 doi: 10.57932/41cf41f5-5a0c-47d0-b88d-a774e06e9012

Crook, C. (2020). Proposal for encoding of a deformation model (Tech. Rep.). Land Information New Zealand. Retrieved from https://
github.com/linz/deformation-model-format/blob/master/doc/
JsonGeotiffDeformationModelFormat_20200501.pdf

Crook, C. (2024). *Topic 24 - functional model for crustal deformation* (Tech. Rep.). OGC DOCUMENT: 22-010r4, version 1.0.4. Retrieved from https://docs.ogc.org/as/22-010r4/22-010r4.html

Delgado, F., & Grandin, R. (2021). Dynamics of episodic magma injection and migration at yellowstone caldera: Revisiting the 2004–2009 episode of caldera uplift with insar and gps data. *Journal of Geophysical Research: Solid Earth, 126*(8), e2021JB022341.

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions. *Geophysical journal international*, *181*(1), 1–80.

Denys, P., & Pearson, C. (2022). Post-earthquake accuracy and precision assessment of nzgd2000. *New Zealand Surveyor*, 306.

de Sagazan, C., Retailleau, L., Gerbault, M., Peltier, A., Feuillet, N., Fontaine, F. J., & Crawford, W. C. (2024). Seismicity near mayotte explained by interacting magma bodies: Insights from numerical modeling. *Journal of Volcanology and Geothermal Research*, 446, 107985.

De Zan, F., Zonno, M., & Lopez-Dekker, P. (2015). Phase inconsistencies and multiple scattering in sar interferometry. *IEEE Transactions on Geoscience and Remote Sensing*, 53(12), 6608–6616.

Doin, M.-P., Lodge, F., Guillaso, S., Jolivet, R., Lasserre, C., Ducret, G., ... Pinel, V. (2011). Presentation of the Small Baselin NSBAS Processing Chain on a Case Example: The Etan Deformation Monitoring from 2003 to 2010 Using Envisat Data. In *Fringe symposium*.

Donnelly, N., Crook, C., Stanaway, R., Roberts, C., Rizos, C., & Haasdyk, J. (2017). A two-frame national geospatial reference system accounting for geodynamics. In *REFAG 2014: Proceedings of the IAG Commission 1 Symposium Kirchberg, Luxembourg, 13–17 October, 2014* (pp. 235–242).

Dupuy, S., Lelong, C., Gaetano, R., & Villers, A. (2024). *Occupation du sol - Mayotte - Carte de 2023*. CIRAD Dataverse. Retrieved from https://doi.org/10.18167/DVN1/CTNRHW doi: 10.18167/DVN1/CTNRHW

EPSG. (2024). *Geodetic Parameter Dataset* (Database). https://epsg.org. (Accessed: 2024-04-14)
ESA Copernicus. (2024a). *ISCE: Interferometric synthetic aperture radar Scientific Computing Environment [Software]*. Retrieved from https://github.com/isce-framework/isce2 (Accessed: 2024-04-19)

ESA Copernicus. (2024b). *Sentinel-1 Level-1 SLC Products [Dataset]*. Retrieved from https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-1-sar/products-algorithms/level-1-algorithms/single-look-complex

Evenden, G. I. (2024). *PROJ - PROJ 9.4.0 documentation* (Software). https://proj.org. (Accessed: 2024-04-14)

Famin, V., Michon, L., & Bourhane, A. (2020). The comoros archipelago: a right-lateral transform boundary between the somalia and lwandle plates. *Tectono-physics*, 789, 228539.

Ferretti, A., Fumagalli, A., Novali, F., Prati, C., Rocca, F., & Rucci, A. (2011). A new algorithm for processing interferometric data-stacks: Squeesar. *IEEE transactions on geoscience and remote sensing*, 49(9), 3460-3470.

Feuillet, N., Jorry, S., Crawford, W. C., Deplus, C., Thinon, I., Jacques, E., ...
others (2021). Birth of a large volcanic edifice offshore mayotte via lithosphere-scale dyke intrusion. *Nature Geoscience*, 14 (10), 787-795.

Gobron, K., Rebischung, P., de Viron, O., Demoulin, A., & Van Camp, M. (2022). Impact of offsets on assessing the low-frequency stochastic properties of geodetic time series. *Journal of Geodesy*, 96(7), 46.

Gómez, D., Piñón, D. A., Smalley, R., Bevis, M., Cimbaro, S., Lenzano, L., & Barón, J. (2016). Reference frame access under the effects of great earthquakes: a least squares collocation approach for non-secular post-seismic evolution. *Journal of Geodesy*, 90, 263-273.

Grandin, R., & Collilieux, X. (2025). *Ground displacement at Mayotte from GNSS campaigns (2004-2023) and levelling campaigns (2006-2023)*. Retrieved from https://doi.org/10.57932/688c6c06-d553-4000-b80c-f95eb92a3b4b doi: 10.57932/688c6c06-d553-4000-b80c-f95eb92a3b4b

Gregg, P., Le Mével, H., Zhan, Y., Dufek, J., Geist, D., & Chadwick Jr, W. (2018). Stress triggering of the 2005 eruption of sierra negra volcano, galápagos. *Geo-physical Research Letters*, 45(24), 13-288.

Gualandi, A., Serpelloni, E., & Belardinelli, M. E. (2016). Blind source separation
problem in gps time series. *Journal of Geodesy*, 90(4), 323-341.

Häkli, P., Evers, K., Jivall, L., Nilsson, T., Himle, S., Kollo, K., ... Lidberg, M. (2023). NKG2020 transformation: An updated transformation between dynamic and static reference frames in the Nordic and Baltic countries. *Journal of Geodetic Science*, 13(1), 20220155.

Häkli, P., Lidberg, M., Jivall, L., Nørbech, T., Tangen, O., Weber, M., ...
Paršeliunas, E. (2016). The NKG2008 GPS campaign–final transformation
results and a new common Nordic reference frame. *Journal of Geodetic Sci-
ence*, 6(1).

Hamling, I. J. (2020). A review of the 2016 kaikōura earthquake: insights from the first 3 years. *Journal of the Royal Society of New Zealand, 50*(2), 226-244.

Henderson, S. T., Delgado, F., Elliott, J., Pritchard, M. E., & Lundgren, P. R. (2017). Decelerating uplift at lazufre volcanic center, central andes, from ad 2010 to 2016, and implications for geodetic models. *Geosphere*, *13*(5), 1489-1505.

Holohan, E. P., Sudhaus, H., Walter, T. R., Schöpfer, M. P., & Walsh, J. J. (2017). Effects of host-rock fracturing on elastic-deformation source models of volcano deflation. *Scientific reports*, 7(1), 10970.

IGN. (2004). RGM04 Coordinate Reference System for Mayotte, EPSG code: 4468 (CRS). https://epsg.org/crs_4468/RGM04.html.

IGN. (2011). *Paramètres de transformation ITRF2014 vers les références nationales*. https://geodesie.ign.fr/files/geodesie/2024-11/ITRF2014_vers_references_nationales(1).pdf. (Accessed: 2025-07-11)

IGN. (2023). *RGM23 Coordinate Reference System for Mayotte*, EPSG code: 10671 (CRS). https://epsg.org/crs_10671/RGM23.html.

IGN. (2024). *MAYG station coordinate time-series from ITRF web site*. https://itrf.ign.fr/plots/ITRF2020-U2023/gnss/plot/TMAYG_90101M001.png.
(Accessed: 2024-07-30)

IGN and Shom. (1953). *MAYO53 Vertical Datum for Mayotte* (Vertical CRS). https://cnig.gouv.fr/IMG/documents_wordpress/2015/03/IGN_MAYOTTE.pdf.

Jacques, E., Hoste-Colomer, R., Feuillet, N., Lemoine, A., van der Woerd, J., Crawford, W., ... Bachlery, P. (2024). Ring faulting and piston collapse in the
mantle sustained the largest submarine eruption ever documented. *Earth and Planetary Science Letters*, **647**, 119026.

Jeanson, M., Anthony, E. J., Charroux, S., Aubry, A., & Dolique, F. (2021). Detecting the effects of rapid tectonically induced subsidence on mayotte island since 2018 on beach and reef morphology, and implications for coastal vulnerability to marine flooding. *Geo-Marine Letters*, **41**(4), 51.

Jolivet, R., Duputel, Z., Riel, B., Simons, M., Rivera, L., Minson, S., ... others (2014). The 2013 m w 7.7 balochistan earthquake: seismic potential of an accretionary wedge. *Bulletin of the Seismological Society of America*, **104**(2), 1020–1030.

Jolivet, R., Duputel, Z., & Simons, M. (2012). *CSI: Classic Slip Inversion [Software]*. https://github.com/jolivetr/csi. (Accessed: 2024-04-14)

Kierulf, H. P., Valsson, G., Evers, K., Lidberg, M., Hakli, P., Prizginiene, D., ... others (2019). Towards a dynamic reference frame in iceland. *Geophysica*, **54**(1), 3–17.

Klein, E., Fleitout, L., Vigny, C., & Garaud, J.-D. (2016). Afterslip and viscoelastic relaxation model inferred from the large-scale post-seismic deformation following the 2010 m w 8.8 maule earthquake (chile). *Geophysical Journal International*, **205**(3), 1455–1472.

Kositsky, A., & Avouac, J.-P. (2010). Inverting geodetic time series with a principal component analysis-based inversion method. *Journal of Geophysical Research: Solid Earth*, **115** (B3).

Kreemer, C., Blewitt, G., & Klein, E. C. (2014). A geodetic plate motion and global strain rate model. *Geochemistry, Geophysics, Geosystems*, **15**(10), 3849–3889.

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.

Le Mével, H., Córdova, L., Cardona, C., & Feigl, K. (2021). Unrest at the laguna del maule volcanic field 2005–2020: renewed acceleration of deformation. *Bulletin of Volcanology*, **83**(6), 39.

Le Mével, H., Gregg, P. M., & Feigl, K. L. (2016). Magma injection into a long-lived reservoir to explain geodetically measured uplift: Application to the 2007–2014
unrest episode at laguna del maule volcanic field, chile. *Journal of Geophysical Research: Solid Earth*, 121(8), 6092–6108.

Lemoine, A., Briole, P., Bertil, D., Roullé, A., Foumelis, M., Thinon, I., ...
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(1),
22–44.

Liang, C., Liu, Z., Fielding, E. J., & Bürgmann, R. (2018). Insar time series analysis of l-band wide-swath sar data acquired by alos-2. *IEEE Transactions on Geoscience and Remote Sensing*, *56*(8), 4492–4506.

Liu, T., Fu, G., She, Y., Meng, G., Zou, Z., Wu, W., ... Pupatenko, V. V. (2023). Post-seismic deformation following the 2011 m w9.0 tohoku-oki earthquake and its impact on northeast asia. *Geophysical Journal International*, 235(2), 1479–1492.

McCarthy, D. D., Petit, G., et al. (2004). *IERS conventions (2003)* (Vol. 32). Bundesamt für Kartographie und Geodäsie Frankfurt am Main, Germany.

Métois, M. (2022). *GPSPCA: Weighted Expectation Maximization Principal Component Analysis for GNSS time-series [Software]*. https://bitbucket.org/mmetois/geodepy. (Accessed: 2022-01-01)

Michon, L. (2016). The volcanism of the comoros archipelago integrated at a regional scale. In *Active volcanoes of the southwest indian ocean* (pp. 333–344). Springer.

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.

Mogi, K. (1958). Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. *Earthq Res Inst*, **36**, 99–134.

Morisseau, S. (2024). *Approche biogéomorphologique des plages de nidification pour la conservation des tortues vertes (Chelonia mydas) de Mayotte* (Thebes, Université de Montpellier). Retrieved from https://hal.science/tel-04776384

Mullet, B., & Segall, P. (2022). The surface deformation signature of a transcrustal, crystal mush-dominant magma system. *Journal of Geophysical Research: Solid
Earth, 127(5), e2022JB024178.

Nikkhoo, M. (2016). *CDM: Compound Dislocation Models [Software]*. https://volcanodeformation.com/onewebmedia/CDM_ECM RD.zip. (Accessed: 2022-10-24)

Nikkhoo, M., Walter, T. R., Lundgren, P. R., & Prats-Iraola, P. (2016). Compound dislocation models (cdms) for volcano deformation analyses. *Geophysical Journal International*, ggw427. doi: https://doi.org/10.1093/gji/ggw427

Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space. *Bulletin of the seismological society of America*, 75(4), 1135–1154.

Peltier, A., Saur, S., Ballu, V., Beauducel, F., Briole, P., Chanard, K., ... van der Woerd, J. (2022). Ground deformation monitoring of the eruption offshore mayotte. *Comptes Rendus. Géoscience*, 354(S2), 1–23.

Rebischung, P., Altamimi, Z., Métivier, L., Collilieux, X., Gobron, K., & Chanard, K. (2024). Analysis of the IGS contribution to ITRF2020. *Journal of Geodesy*, 98(6), 49.

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(6), 4666–4683.

REVOSIMA. (2024a). *Monthly bulletin of the REVOSIMA, November 2024* (Tech. Rep.). IPGP, BRGM, IFREMER, CNRS. Retrieved from https://www.ipgp.fr/volcanoweb/mayotte/Bulletins-Communiques/Bulletins/2024/Revosima_bull_20241206.pdf

REVOSIMA. (2024b). *REVOSIMA: Mayotte volcanological and seismological monitoring network [Dataset]*. http://mayotte.gnss.fr/. (Accessed: 2024-04-19) doi: 10.18715/MAYOTTE.REVOSIMA

Rosen, P. A., Gurrola, E., Sacco, G. F., & Zebker, H. (2012). The insar scientific computing environment. In *Eusar 2012; 9th european conference on synthetic aperture radar* (pp. 730–733).

Rosen, P. A., Gurrola, E. M., Agram, P., Cohen, J., Lavalle, M., Riel, B. V., ... Buckley, S. M. (2018). The InSAR scientific computing environment 3.0: a flexible framework for NISAR operational and user-led science processing. In *Igarss 2018-2018 ieee international geoscience and remote sensing symposium*
(pp. 4897–4900).

République Française. (2019). *Arrêté du 5 mars 2019 portant application du décret n° 2000-1276 du 26 décembre 2000 modifié portant application de l'article 89 de la loi n° 95-115 du 4 février 1995 modifiée d'orientation pour l'aménagement et le développement durable du territoire relatif aux conditions d'exécution et de publication des levés de plans entrepris par les services publics.* (NOR : TRED1803160A, JORF n°0056 du 7 mars 2019, https://www.legifrance.gouv.fr/loda/id/JORFTEXT000038203601)

Segall, P. (2019). Magma chambers: what we can, and cannot, learn from volcano geodesy. *Philosophical Transactions of the Royal Society A*, **377**(2139), 20180158.

Shom. (2016). *MNT Bathymétrique de la façade de Mayotte (Projet Homonim).* https://dx.doi.org/10.17183/MNT_MAY100m_HOMONIM_WGS84. (Accessed: 2024-07-23)

Shom. (2024). *RONIM network. Raw sea-level heights from DZAOU tide-gauge (Shom - Département de Mayotte - Direction de la mer sud océan Indien)*. 13 rue du Chatellier 29200 BREST, France. Retrieved from https://data.shom.fr/donnees/refmar/DZAOUDZI doi: 10.17183/REFMAR#30

SONEL. (2024). *Post-processed sea-level height from DZAOU tide-gauge instrument*. Retrieved from https://www.sonel.org/spip.php?page=maregraphe&idStation=1903

Stamps, D., Kreemer, C., Fernandes, R., Rajaonarison, T., & Rambolamanana, G. (2021). Redefining east african rift system kinematics. *Geology*, **49**(2), 150–155.

Stanaway, R. (2020). *A schema for complex time-dependent transformations between geodetic reference frames* (Unpublished doctoral dissertation). UNSW Sydney.

Stanaway, R., Crook, C., Kelly, K., & Lott, R. (2024). A functional model for quantifying deformation in reference frame transformations. In *Proceedings of iugg general assembly 2023, international association of geodesy symposia*.

Tarrío, J. A., Cáceres, C., Vásquez, V., Marten, M., Inzunza, J., Isla, F., ... others (2024). Geodetic innovation in chilean mining: The evolution from static to kinematic reference frame in seismic zones. *Journal of Geodetic Science*, **14**(1), 20220173.
Tiamo, K., Rundle, J., Klein, W., Ben-Zion, Y., & McGinnis, S. (2004). Using eigenpattern analysis to constrain seasonal signals in southern california. *Computational earthquake science part I*, 1991–2003.

Verdurme, P., Gurioli, L., Chevrel, O., Médard, E., Berthod, C., Komorowski, J.-C., ... others (2024). Magma ascent and lava flow field emplacement during the 2018–2021 fani maoré deep-submarine eruption insights from lava vesicle textures. *Earth and Planetary Science Letters*, 636, 118720.

Walwer, D., Calais, E., & Ghil, M. (2016). Data-adaptive detection of transient deformation in geodetic networks. *Journal of Geophysical Research: Solid Earth*, 121(3), 2129–2152.

Walwer, D., Ghil, M., & Calais, E. (2019). Oscillatory nature of the Okmok volcano’s deformation. *Earth and Planetary Science Letters*, 506, 76–86.

Zinke, J., Reijmer, J., & Thomassin, B. (2003). Systems tracts sedimentology in the lagoon of mayotte associated with the holocene transgression. *Sedimentary Geology*, 160(1-3), 57–79.
References from the supporting information

Amoruso, A., & Crescentini, L. (2009). Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data. *Journal of Geophysical Research: Solid Earth, 114*(B2).

Amoruso, A., & Crescentini, L. (2013). Analytical models of volcanic ellipsoidal expansion sources. *Annals of Geophysics, 56*(4), 0435.

Bailey, S. (2012a). *EMPCA: Weighted Expectation Maximization Principal Component Analysis [Software]*. https://github.com/sbailey/empca. (Accessed: 2024-04-14)

Bailey, S. (2012b). Principal component analysis with noisy and/or missing data. *Publications of the Astronomical Society of the Pacific, 124*(919), 1015.

Beauducel, F., Peltier, A., Villié, A., & Suryanto, W. (2020). Mechanical imaging of a volcano plumbing system from gnss unsupervised modeling. *Geophysical Research Letters, 47*(17), e2020GL089419.

Cesca, S., Letort, J., Razafindrakoto, H. N., Heimann, S., Rivalta, E., Isken, M. P., ... others (2020). Drainage of a deep magma reservoir near Mayotte inferred from seismicity and deformation. *Nature Geoscience, 13*(1), 87-93.

Collilieux, X. (2025). *PROJ implementation of the Mayotte deformation model (2018-2022)*. Retrieved from https://doi.org/10.57932/41cf41f5-5a0c-47d0-b88d-a774e06e9012 doi: 10.57932/41cf41f5-5a0c-47d0-b88d-a774e06e9012

Crook, C. (2020). *Proposal for encoding of a deformation model* (Tech. Rep.). Land Information New Zealand. Retrieved from https://github.com/linz/deformation-model-format/blob/master/doc/JsonGeotiffDeformationModelFormat_20200501.pdf

Davis, P. M. (1983). Surface deformation associated with a dipping hydrofracture. *Journal of Geophysical Research: Solid Earth, 88*(B7), 5826-5834.

Davis, P. M. (1986). Surface deformation due to inflation of an arbitrarily oriented triaxial ellipsoidal cavity in an elastic half-space, with reference to kilauea volcano, hawaii. *Journal of Geophysical Research: Solid Earth, 91*(B7), 7429-7438.

De Zan, F., Zonno, M., & Lopez-Dekker, P. (2015). Phase inconsistencies and multiple scattering in sar interferometry. *IEEE Transactions on Geoscience and Re-
mote Sensing, 53 (12), 6608-6616.

Dupuy, S., Lelong, C., Gaetano, R., & Villers, A. (2024). *Occupation du sol - Mayotte - Carte de 2023*. CIRAD Dataverse. Retrieved from https://doi.org/10.18167/DVN1/CTNRHW doi: 10.18167/DVN1/CTNRHW

Grandin, R., & Collilieux, X. (2025). *Ground displacement at Mayotte from GNSS campaigns (2004-2023) and levelling campaigns (2006-2023)*. Retrieved from https://doi.org/10.57932/688c6c06-d553-4000-b80c-f95eb92a3b4b doi: 10.57932/688c6c06-d553-4000-b80c-f95eb92a3b4b

Le Mével, H., Gregg, P. M., & Feigl, K. L. (2016). Magma injection into a long-lived reservoir to explain geodetically measured uplift: Application to the 2007–2014 unrest episode at laguna del maule volcanic field, chile. *Journal of Geophysical Research: Solid Earth, 121*(8), 6092–6108.

Lemoine, A., Briole, P., Bertil, D., Roullé, A., Foumelis, M., Thinon, I., ... 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(1), 22–44.

Lohman, R. B., & Simons, M. (2005). Some thoughts on the use of insar data to constrain models of surface deformation: Noise structure and data downsampling. *Geochemistry, Geophysics, Geosystems*, 6(1).

Métois, M. (2022). *GPSPCA: Weighted Expectation Maximization Principal Component Analysis for GNSS time-series [Software]*. https://bitbucket.org/mmetois/geodepy. (Accessed: 2022-01-01)

Nikkhoo, M., Walter, T. R., Lundgren, P. R., & Prats-Iraola, P. (2016). Compound dislocation models (cdms) for volcano deformation analyses. *Geophysical Journal International*, ggw427.

Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space. *Bulletin of the seismological society of America*, 75(4), 1135–1154.

Peltier, A., Saur, S., Ballu, V., Beauducel, F., Briole, P., Chanard, K., ... van der Woerd, J. (2022). Ground deformation monitoring of the eruption offshore mayotte. *Comptes Rendus. Géoscience*, 354(S2), 1–23.

Ray, J., Griffiths, J., Collilieux, X., & Rebischung, P. (2013). Subseasonal gnss positioning errors. *Geophysical Research Letters*, 40(22), 5854–5860.
Rebischung, P., Altamimi, Z., Métivier, L., Collilieux, X., Gobron, K., & Chanard, K. (2024). Analysis of the igs contribution to itrf2020. *Journal of Geodesy*, 98(6), 49.

Shom. (2024). *RONIM network. Raw sea-level heights from DZAOU tide-gauge* (*Shom - Département de Mayotte - Direction de la mer sud océan Indien*). 13 rue du Chatellier 29200 BREST, France. Retrieved from https://data.shom.fr/donnees/refmar/DZAOUDZI doi: 10.17183/REFMAR#30

SONEL. (2024). *Post-processed sea-level height from DZAOU tide-gauge instrument*. Retrieved from https://www.sonel.org/spip.php?page=maregraphe&idStation=1903