22 years of satellite imagery reveal a major destabilization structure at Piton de la Fournaise

Flank destabilization of volcanoes is a major hazard because it can trigger tsunamis and large earthquakes, accounting for a significant fraction of volcanic fatalities. On oceanic islands, debris avalanche deposits are widespread, and multiple geological, geophysical, physical and analog modeling studies indicate that magmatic activity can trigger flank destabilization. Magma commonly follows preferential intrusion paths or rift zones, controlled by the edifice stress field and structural weaknesses, and expressed at the surface by aligned cones and fissures. Two principal magma-induced destabilization mechanisms have been proposed: (1) dyke intrusions in near-vertical rift zones coupled with slip on low-angle décollement or detachment faults (documented at Kilauea, Mount Etna, and Cumbre Vieja), and (2) sill intrusions undergoing coeval opening and shear along a pre-existing detachment fault (proposed for Piton des Neiges). In the latter case, magma guided by a prior fault may intrude non-orthogonally to the minimum principal stress, allowing coeval opening and slip; repeated sill emplacement can also promote hydrothermal alteration and creep that weaken the flank. At Piton de la Fournaise (Réunion Island), collapse of the eastern flank is a major hazard. Submarine debris avalanche deposits up to ~100 km³ extend tens of kilometers offshore. In 2007, a major eruption produced up to 1.4 m eastward and 0.35 m upward motion of the eastern flank, and continuous slip at ~1.4 cm/yr has been observed since, with accelerations linked to magmatic activity. A shallow detachment-like fault beneath the eastern flank has been proposed, with modeled depths of ~0.5–1.5 km. These observations raise the question of the link between current magmatic activity and potential catastrophic collapse. Here, leveraging the exceptional activity and geodetic monitoring of Piton de la Fournaise since 1998, the study images in 3D the shallow plumbing system and structures accommodating flank slip, determining the geometry of 57 magmatic intrusions (1998–2020) via state-of-the-art inverse modeling, and relates them to known rift zones.

Prior work has established that ocean-island volcano flanks commonly fail and that magmatic processes can trigger such instability. At Kilauea, dyke intrusions within vertical rift zones push the flank seaward and slip is accommodated on a deep décollement at the edifice–seafloor interface, occasionally producing large earthquakes (e.g., M7.2 in 1975; M6.9 in 2018). At Etna and Cumbre Vieja, shallow detachment-type faults accommodate dyke-induced flank motion. A contrasting mechanism is inferred at Piton des Neiges, where field studies identified a detachment intruded by stacked sills; magnetic fabrics and structural relationships suggest that sill injections can induce rapid co-intrusive slip along the fault. Because magma can be guided by pre-existing structures, intrusions may not be oriented perpendicular to the least principal stress and can undergo coeval opening and shear. Over longer timescales, heat and fluids from cooling sills can promote hydrothermal alteration and creep, weakening the detachment. At Piton de la Fournaise, bathymetry documents large submarine collapse deposits, and geodetic observations since 2007 reveal ongoing, accelerating flank slip correlated with intrusive activity. Previous inversions of 2007 co-/post-eruptive deformation suggested a shallow detachment sub-parallel to topography at ~0.5–1.5 km depth. However, prior geodetic analyses focused on discrete dykes and lacked a comprehensive 3D characterization of continuous failure structures and their relationship to intrusion pathways.

Data: 57 magmatic intrusions (54 eruptions and 15 failed eruptions; 12 events with <2 cm deformation discarded) between March 1998 and December 2020 were imaged by InSAR (OI², nine satellite missions) and/or GNSS (OVPF-IPGP campaigns and permanent network). Two events (Nov 2002, May 2003) used GNSS only. Previously published geodetic models were incorporated where appropriate; several GNSS-only models were recomputed using InSAR. Inverse modeling: A fully 3D Mixed Boundary Element Method was used, assuming a linearly elastic, homogeneous, isotropic medium (Young’s modulus 5 GPa; Poisson’s ratio 0.25) with realistic topography. Fractures (intrusions or faults) were represented by quadrangular surfaces meshed with triangular elements that can be curved along strike and/or dip, allowing sill–dyke transitions and capturing asymmetric surface deformation. Unlike conventional displacement-based inversions, the method solves for constant boundary stress conditions on fractures: a uniform overpressure (magma pressure minus normal host-rock stress) and a uniform shear stress change (tangential to the fracture). These stress boundary conditions produce smooth, physically realistic displacement fields and can better inform host-rock stress. Curvature and allowance for coeval shear and normal stress changes were included when they improved data fit (as indicated by statistical criteria such as Akaike Information Criterion). Inversion workflow: A two-stage process was used. (1) Search stage: parameter space exploration via forward computations combined with a Neighborhood Algorithm to find models that minimize least-squares misfit between observed and modeled displacements. (2) Appraisal stage: Bayesian inference on the ensemble from the search stage to compute mean models, one- and two-dimensional marginal posterior probability density functions (PPDs), confidence intervals, and parameter trade-offs. Uncertainty and statistical representation: For each event, random meshes were generated within the 95% confidence interval of the best-fit model. The ensemble of triangular elements was used to map spatial probability density of intrusion locations, and to compute stereographic projections (poles) of element orientations (Schmidt distribution). This yielded probabilistic 3D geometry, dip/strike distributions, and confidence envelopes for each intrusion zone. Event selection and quality: Eight small-volume (<0.3 Mm³) summit intrusions with low deformation signals (<5 cm) had large uncertainties and were excluded from synthesis. GNSS-only data near the summit were recognized as insufficient to capture distal deformation; where possible these were updated with InSAR-based inversions. Supporting analyses: The 3D intrusion models were compared with seismicity (pre-eruptive swarms and eastern flank events), known rift-zone orientations, and topographic/structural constraints. Polynomial regression surfaces (quadratic) were fit to mesh points of best-fit and 95% CI models to estimate the overall geometry of the primary structure and to quantify its dip variation with elevation and its potential failure-surface extent and volume.

- Comprehensive 3D imaging of 57 intrusions (1998–2020) reveals a major arcuate NE–SE dyke intrusion zone that connects at depth to an eastward, low-dip sill intrusion zone, together forming a spoon-shaped structure from summit to offshore. - Volume partitioning by intrusion zones: NE–SE rift zone: 22 events; 33.1×10^6 m³ (~45% of total intruded volume). Sill zone beneath eastern flank: 7 events; 25.9×10^6 m³ (~35%). N120° zone: 10 events; 8.5×10^6 m³ (~12%). N210° zone: 4 events; 3.33×10^6 m³ (~4%). N300° zone: 2 events; 1.41×10^6 m³ (~2%). N60° zone: 6 events; 1.33×10^6 m³ (~2%). Together, the NE–SE and sill zones account for ~80% of intruded volume but only ~50% of intrusion count. - Geometry: Near the summit, intrusions are planar and sub-vertical; away from the summit they curve laterally and vertically, transitioning to low-dip sills (~0–20° eastward) beneath the eastern flank, with the overall structure following topography and dipping ~15° at its base. Dip–depth distributions from models agree with global field datasets: shallow intrusions (<1 km) tend to be sub-vertical, deeper intrusions (>1 km) have lower dips (10–40°). - Displacement continuum and increasing slip eastward: Western sills show pure opening; three more easterly intrusions (Jan 2004, Oct 2019, Sep 2020) require significant shear stress in addition to overpressure, indicating coeval opening and slip (sheared sills). The easternmost (Oct 2019) behaves predominantly as a fault (shear stress ≫ overpressure). Beyond the sills, the March–May 2007 flank-slip source is modeled as pure shear (fault), indicating a transition from sills to fault slip toward the sea. - Structural continuity and confidence: Overlapping 95% confidence intervals of neighboring models and similar accommodated volumes in the NE–SE and sill zones confirm a single connected spoon-shaped structure. Polynomial-regressed surfaces fit to model ensembles depict its spoon-shaped geometry in map and cross-section. - Seismicity relation: Beneath the summit, the structure overlies the pre-eruptive seismic swarm, consistent with a critically stressed volume below and a mobile volume above. Eastern flank seismicity dips ~45° and is temporally correlated with intrusive events above sea level, likely reflecting stress transfer rather than direct slip or intrusion propagation. - Secondary intrusion zones: Four additional, shallower rift/intrusion zones (N60°, N120°, N210°, N300°) are delineated; they collectively accommodate a minor share of volume (each 2–12%). Their geometry and shallow depths suggest induction by stress build-up from the main spoon-shaped structure rather than inheritance from deep crustal faults. - Hazard metrics: From the structure geometry, the potential failure surface associated with a future collapse is estimated to be ~8.5±2.5 km long and ~10±2.5 km wide, with mean slope ~15°. Estimated unstable volume above the surface is ~45 km³ (95% CI: 14–96 km³), comparable to volumes of large ocean-island flank collapses. - Kinematics and forcing: The eastern flank exhibits both rapid co-intrusive slip (during sill emplacement or fault activation) and slower, steady inter-eruptive slip (~1.4 cm/yr reported previously), with accelerations following intrusive events. Co-intrusive dyke opening in the arcuate NE–SE zone pushes the flank seaward and transfers shear onto the low-dip sill/fault surface, producing a continuum from opening to sheared sills to pure fault slip.

The 3D geodetic inversions demonstrate that Piton de la Fournaise is governed by a single, connected spoon-shaped destabilization surface that couples vertical dyke intrusion at the head to low-dip sill intrusion and, ultimately, to fault slip toward the toe. This structure provides a mechanistic bridge between the two classically proposed models: dyke-induced slip on low-angle faults (as at Kilauea/Etna) and sill-induced coeval opening and slip along a detachment (as at Piton des Neiges). The observed displacement continuum—from pure opening near the summit, through sheared sills beneath the east flank, to pure slip at the outer toe—directly addresses how magmatic activity can progressively load and activate a destabilization surface. The upper (dyke) part of the structure may be shaped by tectonic/topographic loading and buoyancy, whereas the lower sill/fault portion appears inherited, guided by pre-existing weaknesses (e.g., base of the edifice near sea level, ductile hyaloclastite layers, or the roof of the Les Alizés hypovolcanic complex). Dyke opening in the arcuate NE–SE rift zone drives seaward motion that increases shear on the low-dip surface, establishing a positive feedback: intrusions in the head promote slip at the base, which in turn may help maintain the arcuate rift geometry. The easternmost fault section appears locked most of the time, accumulating stress and releasing it via slow slip (e.g., 2007) or potentially in larger, hazardous earthquakes and tsunamis. The findings reconcile geodetic observations of co- and inter-eruptive deformation with seismological constraints, while explaining the spatial distribution of secondary rift zones as stress-induced features. They also emphasize that long-term intrusive activity and repeated stress perturbations can prime the flank for large-scale failure. The kinematic and geometric framework likely extends to other ocean-island volcanoes where similar spoon-shaped intrusion–fault systems could control flank movements.

This study provides the first long-term, comprehensive 3D image of Piton de la Fournaise’s shallow plumbing system and associated destabilization surface, using 22 years of satellite geodesy and advanced inverse modeling. A major arcuate NE–SE dyke zone connects to a low-dip sill zone and evolves into a fault toward the sea, forming a spoon-shaped structure that accommodates a displacement continuum from pure opening to pure slip. The structure channels ~80% of intruded magma volume, explains observed co-intrusive and steady flank motions, and yields quantitative hazard metrics (failure-surface extent and potential unstable volume on the order of tens of cubic kilometers). The results reveal a hybrid destabilization mechanism combining dyke-induced loading with sill-guided coeval slip and faulting, and highlight a positive feedback between rift-zone opening and basal slip. Future work should include: (i) extending similar long-term geodetic inversions to other ocean-island volcanoes to test the ubiquity of spoon-shaped destabilization structures; (ii) integrating inelastic and heterogeneous rheologies to refine stress and slip estimates; (iii) coupling geodetic inversions with seismic, petrological, and hydrothermal observations to assess temporal evolution of weakening and creep; and (iv) improving offshore imaging to constrain the distal fault geometry and its locking–failure behavior.

- Modeling assumptions: elastic, homogeneous, isotropic medium; constant overpressure and shear stress on fractures; depth-dependent stress gradients are not constrained at intrusion scale. These simplifications may bias depth and volume estimates (topography omission can induce ~30% depth and ~20% volume errors; here, realistic topography is included, but material heterogeneity and inelastic effects are not). - Data coverage: No eruptions outside Enclos Fouqué since 1998 limit imaging of the rift zone outside the caldera walls. Two events rely on GNSS only, and GNSS network geometry around the summit misses distal deformation; several small summit events (<5 cm signal; <0.3 Mm³) were excluded due to high uncertainty. - Seismic resolution: Lack of event migration above ~1 km asl hinders precise imaging of the very shallow plumbing; eastern flank seismicity at ~45° dip is only indirectly linked via stress transfer. - Inversion non-uniqueness: Results represent families of acceptable models within 95% confidence intervals; eight inversions have large uncertainties due to low signal-to-noise ratios. - Generalizability: While the intrusion orientations over 1998–2020 match earlier decades, extrapolation to deeper structures or to other volcanoes requires caution and further site-specific validation.