Volcanic activity can induce flank failure, sometimes generating large earthquakes and tsunamis. However, the failure structures have never been fully characterized and the failure mechanism is still debated. Magmatic activity is a possible trigger, either through fault slip, which might be induced by dyke intrusions, or through sill intrusions, which might be undergoing coeval normal displacements and slip. At the Piton de la Fournaise volcano, satellite imagery combined with inverse modeling highlights the pathways of 57 magmatic intrusions that took place between 1998 and 2020. We show that a major arcuate dyke intrusion zone is connected at depth to a sill intrusion zone, which becomes a fault zone towards the sea, forming a spoon-shaped structure. Some sills are affected by coeval normal displacement and seaward slip. Overall, the structure is characterized by a continuum of displacement from no slip, to sheared sills and finally pure slip. Repeated intrusions into this spoon-shaped structure could trigger catastrophic collapses.
Flank destabilization of volcanoes is a major hazard as it can trigger tsunamis and large earthquakes, which, when combined, account for 24% of volcanic fatalities1. Large flank destabilizations are common on oceanic islands, as demonstrated by bathymetric studies which show ubiquitous debris avalanche deposits2. Evidence from geological3,4, geophysical observations5,6,7,8,9, and analysis using physical and analog models10,11,12,13,14 indicate that flank destabilization at these volcanoes may be triggered by magmatic activity.
This activity often occurs along preferential intrusion paths, generally referred to as “rift zones”15, which are, in turn, controlled by the stress field of the edifice10,16 and by structural weaknesses17,18. Rift zones can be fed by magma either sporadically or continuously, and vary in size from a kilometer to tens of kilometers long. Despite these differences, their surface expression is characterized by a high density of pyroclastic cones and eruptive fissures aligned along one or several preferential directions along which magma may be guided away from a central vertical conduit.
To date, two mechanisms of magma-induced flank destabilization have been proposed. The first mechanism (Fig. 1a), corresponding to subvertical intrusions (dykes) in rift zones coupled to a low angle fault5,6,7,10,11,12,19, was proposed for Kilauea (Hawaii, U.S.A), Mount Etna (Italy) and Cumbre Vieja (La Palma, Canary Islands) volcanoes. At Kilauea, geodetic data reveal that intrusions within vertical rift zones push the southern flank of the volcano seaward5,7, episodically leading to flank failure. Flank slip is accommodated by a deep décollement-type fault, located at the interface between the sea floor and the edifice, sometimes generating large earthquakes, such as the M7.2 in 197520 and the M6.9 in 201821. At Etna and Cumbre Vieja, the faults accommodating dyke intrusions are detachment-type faults, located at a shallow level in the edifice6,19. The second mechanism (Fig. 1b), corresponding to subhorizontal intrusions (sills) undergoing coeval opening and slip3,22, was proposed for Piton des Neiges, the extinct neighbor volcano of Piton de la Fournaise (Réunion Island). Identification of a detachment-type fault intruded by a pile of sills, separating an extinct gabbroic magma chamber from debris avalanche deposits, led to the suggestion that sill injections into the fault could induce rapid co-intrusive slip3, as shown by the seaward oriented magnetic fabric observed in the sills4. Because a previous structure acts as a guide for the intrusions, magma is not necessarily emplaced perpendicular to the minimum principal stress, in which case the intrusion may undergo coeval opening and slip23. If the previous structure is a fault and the fault extends beyond the intrusion zone22, the intrusion of magma may induce failure of the whole fault, possibly leading to flank collapse. In the longer term, heat from the cooling sills in the detachment promotes syn-deformation hydrothermal alteration which weakens the surrounding rocks, leading to low-temperature creep3,24.
At Piton de la Fournaise, collapse of the eastern flank represents one of the major hazards. Bathymetric studies have shown that debris avalanche deposits of up to 100 km3 are present on the submarine flanks of the volcano at distance up to 80 km from the shoreline25. Moreover, it has recently been demonstrated that flank slip is still ongoing: in March–May 2007, the largest eruption of the last hundred years took place on the lower eastern flank of the volcano. Synthetic aperture radar interferometry (InSAR) showed that the eruption was associated with an eastward slip of the eastern flank of the volcano of up to 1.4 m, with 0.35 m of upward motion26 (see Fig. 2). Since then, InSAR and global navigation satellite system (GNSS) data have shown that there is a continuous slip at a rate of 1.4 cm/yr9,27,28, which accelerates with magmatic activity29,30,31. The presence of a fault under the eastern flank (Grandes Pentes area in Fig. 2) was proposed to explain the co- and post-eruptive deformation recorded during the 2007 eruptive crisis23,32,33. Inverse modeling of post-eruptive deformation32 points to a fault plane sub-parallel to the topography with a depth ranging from ≈500 to ≈1500 m depending on the parameters considered in the inversion (Fig. 3). These observations raise the question of a link between the present magmatic activity and the possibility of a catastrophic flank collapse at the volcano. In addition, Piton de la Fournaise is one of the most active volcano in the world with an average of one eruption every five months and one of the best monitored with continuous displacement monitoring by both InSAR and GNSS since 1998.
Here, we take advantage of the exceptionally high activity of Piton de la Fournaise and the large amount of displacement data to image the shallow plumbing system in 3D, as well as the structures which accommodate flank slip. The geometry of 57 magmatic intrusions that occurred between 1998 and 2020 are determined relying on state-of-the-art inverse modeling in a fully 3D framework, allowing for curved intrusions and stress changes. We link the imaged intrusions to the rift zones previously recognized at the surface. A major spoon-shaped active structure, accommodating coeval opening and slip is highlighted. We compare this structure to other volcanoes and discuss the implication for the stability of the volcano.
Piton de la Fournaise plumbing system
Piton de la Fournaise is located in the south-east of Réunion Island in the Indian Ocean (Fig. 2a). The volcano is buttressed by the older volcanic edifice of Piton des Neiges, located to the north-west. Piton de la Fournaise is built on an intrusive gabbroic complex possibly originating from an older edifice (Les Alizés), which has been inferred as lying to the east of the volcano from drill-hole34 and gravity studies35 (Figs. 2a and 3).
Located inside a U-shaped collapse structure (Enclos Fouqué), the current volcanic center has formed a 400 m-high central cone with a 1 km-wide and 300 m-deep summit crater (Dolomieu) (see Fig. 2b). Several rift zones (Fig. 2c), radiating from the central cone, have been identified from locations and orientations of superficial morpho-structural markers, such as eruptive fissures or cinder cones. The main rift zone has an overall NE-SE arcuate shape passing through the central cone36,37,38. It extends beyond the walls of the Enclos Fouqué, both to the north and south, and reaches down to the shore. Additional, less active, intrusive directions have also been identified: one, N120∘, running from the southeast edge of the Dolomieu crater to the base of the central cone18,37, an East Rift Zone striking N60∘36,38, a N300-325∘ direction, and a volcanic alignment toward the south-south-west37.
The deeper part of the plumbing system has been interpreted mainly through seismicity (Fig. 3). A sub-vertical conduit, extending between −5 km and sea level, and located 500–1000 m west of the summit was outlined by seismic event migrations during the 1998 and 2014 eruptive crises39,40. Its location was confirmed by P-wave velocity tomography41 and the modeling of the deformation recorded at the beginning of the 2014 crisis42. Around sea level, the presence of a main shallow reservoir is postulated from a discontinuity in the location of micro-earthquakes39 and the P-wave velocity structure41, as well as from pre-eruptive deformation8,42. Pre-eruptive seismic swarms also occur before eruptions between sea level and 1 km above sea level (asl)40,43,44. However, no event migration has been identified above 1 km asl, preventing to precisely image the shallow plumbing system. The shallower part is known through deformation. The analysis of InSAR and GNSS data associated with 16 previous eruptions32,33,45,46,47,48,49,50 highlights shallow discrete dykes, but no continuous structure. Our study provides the first extensive image of Piton de la Fournaise shallow plumbing system.
Inverse modeling of intrusions
Our inverse modeling approach has two characteristic features (see “Methods”). The first is that intrusions can be curved. Indeed, curved intrusions are sometimes required to account for the strong asymmetrical displacement50,51,52,53, characteristic of most Piton de la Fournaise eruptions. This asymmetry is best explained by curved intrusions corresponding to sill-to-dyke transitions48,49,50. Intrusions are meshed by triangular elements, making it possible to create such complex geometries. The second feature is that stress changes are determined on fractures, as opposed to most geodetic studies which determine displacement on fractures20,54. Here, we chose to determine overpressure and shear stress changes (normal and tangential to the fracture elements, respectively), as these quantities are more physical and more informative than displacements. This study is the first that systematically searched for possible curvatures and allowed coeval shear and normal stress changes, when it improved data fit.
Between March 1998 and December 2020, 54 intrusions leading to eruptions and 15 failed eruptions (i.e., intrusions that did not reach the surface, but created a seismic crisis and deformation) occurred. All were captured by either InSAR or GNSS, or both (see “Methods”). Among these intrusions, we analyzed 57 events and discarded 12 events which did not produce deformation detectable by InSAR (>2 cm). Of the 57 intrusion models considered for our analysis, 16 had already been determined and 41 are brand new computations. Among the already determined models, 8—based on InSAR data—were not recomputed (March 1998, July and September 1999, February 2000 and June 200045,48, March–May 200732,33, May 201649 and July 201750), while 8 others—solely based on GNSS data—were recomputed based on InSAR data (August and September 200346, May and August 2004, February, October, November and December 200547). Indeed, the GNSS network is located around the summit and misses displacements associated with eruptions beyond the summit area49.
Of these models (see Supplementary Figs. S6–S62 for detailed model results), 8 have large uncertainties because of low signal-to-noise ratios of displacement measurements, which come from low displacement amplitudes (<5 cm). These models mainly correspond to small volume (<0.3 Mm3) summit intrusions in the Dolomieu crater55. Because of their small volume, these intrusions are not expected to provide any relevant information about the plumbing system at the scale of the edifice. The high uncertainties of the deformation data and the negligible associated volumes led us to discard them in our analysis.
All inverted models show a consistent spatial distribution, with swarms of intrusions close to each other, in agreement with the similar displacement patterns observed by InSAR at the volcano49,50. With the exception of one model (February 2007 at the summit), all models delineate six preferential intrusion zones radiating from the central cone. Five coincide with the rift zones previously identified from surface features18,36,37,38, while the sixth, corresponding to a sill intrusion zone at depth, is identified for the first time by our 3D comprehensive study (Fig. 4).
All intrusions are rooted beneath the Dolomieu crater at 0.5–1.5 km asl, consistent with the analysis of continuous tilt and GNSS data which show that magma starts migrating vertically in the conduit beneath the summit before being transmitted to the vent either vertically for summit eruptions or laterally for rift zone eruptions and sill intrusions49,56. Our study confirms that the Dolomieu crater is a turning point for magma.
A major spoon-shaped intrusion zone accommodating seaward slip
The path with the highest number of intrusions (22 events, represented in red in the Fig. 4), connects at the surface to the already recognized major NE-SE rift zone shown on Fig. 236,38. Since this rift zone accounts for the largest intruded volume (33.1 × 106 m3, totaling 45% of the intruded volume), it confirms its role as the main preferential intrusion zone. The lack of eruptions outside the Enclos Fouqué since 1998 prevents to image the 3D shape of this rift zone outside the caldera wall. Intrusions in the summit area are planar and vertical, while those away from the summit to the NE and SE curve horizontally and vertically. Together, they form a spoon-shaped structure with a horizontal curvature, as shown by the range of strikes from N30∘ to NS on the north flank and NS to N120∘ on the south flank (Fig. 4a), and a vertical curvature, as shown by the transition from gently dipping to highly dipping intrusion with increasing elevation (Fig. 4b). This structure is consistent with a joint analysis of continuous GNSS and InSAR data which determined that an intrusion may start as a sill beneath the summit, propagate laterally before becoming vertical to reach the vent49.
A sub-horizontal preferential intrusion zone is shown through 7 sill models beneath the eastern flank (shown in orange on Fig. 4). This intrusion zone can also be considered as major as it corresponds to the second largest intruded volume after the main rift zone (25.9 × 106 m3), amounting to 35% of the intruded volume. All 7 intrusions have a low eastward dip (0–20∘), following the topography. At depth, they connect to the main NE-SE intrusion zone, extending the spoon-shaped structure to the east (Fig. 5).
Both the NE-SE and the sill intrusion zones form a single major structure, as indicated by the overlap of 95% confidence intervals of neighbor best-fit models and by the similar volumes that the NE-SE and the sill intrusion zones accommodate (45 and 35% of the intruded volume, respectively), which confirms a continuity of dilations. Together, both intrusion zones amount to 80% of the injected volume, while constituting only 50% of the number of intrusions. 3D animations highlight this structure in the Supplementary Movies S1, S2. The overall shape of the structure was estimated by quadratic polynomial regression considering the mesh points of (1) best-fit models and of (2) models randomly generated within the 95% confidence intervals. The surface determined is shown in black in Fig. 5.
This main spoon-shaped structure is characterized by vertical dips in the rift zone, progressively decreasing towards the east to reach ≈15∘ beneath the east flank, where the structure follows the topography. Dips of this structure are in agreement with field observations. Indeed, at Piton des Neiges57 and several other volcanoes worldwide (Koolau, Waianae, Tutulia, Tenerife, Stromboli, Piton de la Fournaise23), intrusion outcrops show preferential dips that correlate with emplacement depths (Fig. 6): sub-vertical intrusions are mainly found in the shallow part of edifices (<1 km), while low dip intrusions (10–40∘) are found at greater depths (>1 km). This distribution is close to that inferred from our models (Fig. 6), indicating that our curved models are highly plausible. Moreover, eruptions that occurred between 1998 and 2020 show vent locations and orientations similar to those of the 1932–1992 period (see Fig. 2b), which indicates that the results determined for the recent period probably also apply to the last 100 years.
A striking result is the continuously increasing slip towards the east. Indeed, the four westernmost sill intrusions show pure opening, while the three easternmost ones (January 2004 eruption, September 2020 failed eruption, and October 2019 eruption from west to east) require additional shear stress to the overpressure in order to explain the meter scale lateral displacement measured for each (see individual models in Supplementary Figs. S21, S31, S58, and S61). January 2004 and September 2020 events have shear stress changes of the same order as the overpressures. These sheared sills confirm the field observations3,4 of sheared intrusions accommodating flank slip. The easternmost intrusion (October 2019) is a shallow plane, predominantly undergoing slip, corresponding to a shear stress change one order of magnitude greater than the overpressure (see Supplementary Table S1). Thus, this intrusion behaves almost like a fault. East of the seven sills, the source of March–May 2007 flank displacement was only affected by shear stress, suggesting that this fracture was a fault32,33.
Beneath the Dolomieu crater, the spoon-shaped structure follows the top of the pre-eruptive seismic swarm (Fig. 7). This swarm is consistent with a critically stressed volume beneath this structure and a mobile, less stressed, volume above. However, further east, the structure does not coincide with the seismicity of the eastern flank, which follows a 45∘ dip as shown by locations from the Observatoire Volcanologique du Piton de la Fournaise from the Institut de Physique du Globe de Paris (OVPF-IPGP) (in blue on Figs. 3 and 7). Therefore, this seismicity is unlikely to be directly related to slip of the eastern flank or to the propagation of intrusions. However, the temporal correlation of east flank seismicity with intrusive events in the upper part of the edifice above sea level (as in September 2020,31), suggests that this seismicity is triggered by stress transferred by magma intrusions.
Secondary intrusion zones
Four secondary intrusion zones are highlighted. Two are located to the east of the summit. One, following a N60∘ direction, is inferred from six intrusive events (in blue on Fig. 4). This intrusion zone has been alluded to in previous work. It could correspond to the East Rift Zone36,37,38, which has a similar strike (N50-70∘), but is located closer to the central cone (Fig. 2). Dykes are planar and sub-vertical (see the stereographic projection in Fig. 4). They are also shallow (rooting between 1400 and 1800 m asl) and short, making the total volume of intruded magma, the smallest (1.33 × 106 m3, 2% of the intruded volume) of all the intrusion zones. The other intrusion zone, following a N120∘ direction, is delineated by 10 intrusive events (Fig. 4 in magenta). It has been described previously18,37, but it was assumed to be limited to the summit cone. Our models show that it extends up to the Grandes Pentes. This magma pathway accommodates 12% of the intruded magma volume (8.5 × 106 m3), which makes it the third-largest intrusion zone. Dips are mainly sub-vertical (Fig. 4), sometimes with curvature at depth. Two intrusions are planar and have 60∘ dips, resulting in 60∘ being the most frequent dip in Fig. 4. The bottom of this structure is shallow, following a mean slope of ≈10∘, sub-parallel to the topography (Fig. 4).
West of the summit, two secondary intrusion zones are indicated by curved eastward-dipping intrusions, which root beneath the summit before becoming vertical toward the vent. A N210∘ intrusion zone (Fig. 4 in cyan) is delineated by four eruptions. This path has been described previously18. It appears to be an extension of a volcanic alignment on the south-western flank of the volcano37 (SSW on Fig. 2). It is the fourth largest in terms of intruded volume (3.33 × 106 m3), but nevertheless only corresponds to 4% of the total intruded volume. Finally, a N300∘ intrusion zone (Fig. 4 in yellow) is indicated by two intrusions, representing 2% of intruded volume (1.41 × 106 m3).
It is likely that these secondary intrusion zones have been induced by the main spoon-shaped structure rather than inherited from crustal fractures. Indeed, the N120∘ rift zone is shallow with a lower limit following that of the NE-SE preferential intrusion zone (1000 m asl beneath the Dolomieu crater and 500 m asl beneath the Grandes Pentes). The shallow depth of this rift zone confirms, as previously suggested58, that it is probably induced by stress build-up in the main rift zone rather than by an inherited crustal fracture. The N60∘ intrusion zone also consists of shallow dykes. Its strike has the same angle from the SE branch of the main rift zone as the N120∘ intrusion zone from the NE branch of the main rift zone (Fig. 7), indicating that it might also result from stress build-up from the main rift zone. The N210∘ and the N300∘ intrusion zones could be related to the southern and northern extensions of the NE and SE branches of the main rift zone, respectively. Indeed, it is known that fluid-filled fractures tend to propagate along the same direction59. These rift zones lie in the continuation of the NE and SE branches of the main rift zone. The repeated intrusions in this rift zone might have increased the tensile stress along the main rift zone directions, despite it crossing the summit.
Combining InSAR and inverse modeling, we obtained an exceptionally high level of detail on a continuum of fractures, dips, and displacements (Figs. 5 and 7). The upper part of the imaged structure, which channels dykes, could result from tectonic and topographic loading, combined with buoyancy forces, as suggested for Sierra Negra volcano, Galapagos Islands60. However, its lower part is probably inherited. Indeed, slip is repeatedly induced along the sill part of the structure, indicating that magma is not intruded perpendicular to the least principal stress, but is instead guided by a pre-existing structure22,23. Moreover, the base of the structure is correlated with several known features in the edifice (Fig. 7). Its lower limit connects eastward to the base of the Piton de la Fournaise edifice61 (around sea level), which is known to guide magma, as indicated by a horizontal shift in seismic events recorded during the 1998 magma migration39. It could also correspond to ductile hyaloclastite layers34 or to the roof of the hypovolcanic complex of Les Alizés34,35, whose depths are in the same range.
The overall spoon-shaped geometry (Fig. 5) is suggestive of a rotational landslide, which could act as a destabilization surface. However, this structure differs from classical landslides, in that it accommodates magma. At the head, the structure undergoes pure opening pushing the eastern flank seaward. The associated deformation steepens the edifice slopes52,58 and transfers shear stress, which favors slip on the sub-horizontal surface. Further down, as the structure becomes horizontal, slip coeval with opening takes place, ultimately leading to pure fault slip at the toe of the structure (Fig. 7b, c). The easternmost part of the structure is probably locked as indicated by the abrupt transition from slip to a lack of slip towards the east.
The eastern flank steady displacement9,28,62 takes place in the same area as the co-intrusive flank slip (Fig. 7a). This displacement shows accelerations after intrusive events29,30,31. Thus, the slip surface is affected by rapid co-intrusive displacement, and slower steady inter-eruptive displacement. Rapid co-intrusive slip takes place during sill emplacement, such as in January 2004, or fault activation, such as in 2007. In contrast, slower steady displacement may arise from low temperature creep induced by sill emplacement which heats up the host rocks and induces fluid circulation, promoting hydrothermal alteration and weakening of the slip surface3,24. Creep on these surface is able to accelerate in response to co-intrusive stress perturbation. To the difference of Kilauea, where slip on a decollement (a reverse fault) is driven by continuous magma intrusions or the creep of dense cumulates in the rift zones5,63, the steady flank displacement takes place on a detachment (a normal fault) and thus may be driven by gravity, as assumed at Etna6 or Cumbre Vieja19. Both the steady and the co-intrusive slip release stresses within the rifts zones and probably participate to shape the rift zone as an arcuate NE-SE rift zone13. Thus, we propose a spoon-shaped intrusive structure characterized by a positive feedback: slip at the base of the structure is induced in response to dyke opening in NE-SE rift zone, which in turn favors the opening of new dykes in the arcuate NE-SE intrusion zone.
Because of friction, the fault to the east of the sill zone locks the system, accumulates stress, and sporadically relaxes them, either in slow slip events, as in 200732,33, or suddenly, possibly triggering large earthquakes and tsunamis, as proposed at several ocean islands2,5,10,11,12,13,14.
The mechanism revealed by our study is in agreement with that inferred at Piton des Neiges where sills intrude a preexisting fault and induce slip (Fig. 1b;3,4,22,64). To the difference of these studies, we observe that fault slip may be induced even in absence of sill intrusion. Thus, the flank slip mechanism we propose is hybrid between models proposed for shield volcanoes: (1) dyke intrusions coupled with slip on low angle faults, whether decollement or detachments5,6,10 (Fig. 1a); and (2) sheared sills intruded into faults3,22 (Fig. 1b).
From the geometry of the structure, we estimate the average length of the failure surface of a potential collapse to be about 8.5 ± 2.5 km for a width of 10 ± 2.5 km, while the mean slope of the destabilization surface from head to toe is about 15∘. The estimated volume above the surface is 45 km3 (ranging between 14 and 96 km3 in the 95% confidence interval of the polynomial regression). This volume lies within the range of massive flank collapse scars observed on other oceanic islands. At Fogo (Cape Verde), the destabilization surface of the last collapse event is estimated to be ≈10 km wide and ≈12 km long for a volume of 110 km3 65, and, at Tenerife (Canary Islands), the failure surface of the Icod landslide is estimated to be 15–20 km by 20–25 km66 for an estimated deposit volume of around 80 km3 67. The structure we have imaged concerns Piton de la Fournaise, but it highlights the importance of comprehensive long-term studies of intrusive activity. At other volcanoes, similar geodetic analysis could reveal that displacements continuum—ranging from pure opening to coeval opening and slip, and to pure slip—induced by the intrusion of magma into spoon-shaped structures also control flank movements.
InSAR and GNSS data
The 57 intrusions with significant deformation (>2 cm) that took place between 1998 and December 2020 at Piton de la Fournaise were imaged by InSAR or by GNSS. InSAR data were acquired by the Indian Ocean InSAR Observatory (OI2), a component of the French National Service of Volcanological Observation (INSU/CNRS)53, which carries out continuous InSAR monitoring on Piton de la Fournaise using data from nine satellite missions (see Supplementary Fig. S5,53), while GNSS data were acquired through reiteration campaigns and by the permanent network of OVPF-IPGP. OI2 imaged 55 of the total of 57 eruptions and failed eruptions with at least one usable ascending or descending interferogram. For two eruptions (November 2002 and May 2003), for which InSAR data are not available, data from the GNSS campaigns were used.
We carried out inverse modeling in order to characterize intrusions. Inverting consists of finding the intrusion that minimizes the difference between the observed and modeled displacements, in the least square sense. Model computations are based on a 3D Mixed Boundary Element Method68, which assumes a linearly elastic, homogeneous, and isotropic medium with a Young’s modulus of 5 GPa and a Poisson’s ratio of 0.25, following previous estimations for the volcano45,69. Realistic topographies are taken into account. At Piton de la Fournaise, failure to take topography into account can induce errors in the estimation of an intrusion depth and volume45 by 30 and 20%, respectively. Our models take fracture curvature into account, which can vary along the strike and dip directions. The along dip curvature accounts for the clear asymmetry of surface displacement often observed for magma intrusion at Piton de la Fournaise50,53. Previously, this asymmetry has been explained by (1) low eastward dip of dykes as in an elastic homogeneous medium45,48,70, or (2) complex mechanical behaviors, such as elastic heterogeneities51 or elasto-plasticity52. The complex mechanical behaviors were inferred from qualitative 2D data comparisons. In contrast, the homogeneous elastic models are 3D and inversions are used to explain the data quantitatively. Using these models, it was determined that along-dip curvature provides a better compromise between the data fit and the model simplicity than planar sources as indicated by the lower Akaike Information Criterion48,50. Thus, we approximate intrusion geometries by quadrangular surfaces which could be curved, with curvatures along the strike or dip directions. These fractures could be linked to the surface topography or not, depending on whether the intrusion led to the eruption of magma or resulted in a failed eruption. An intrusion geometry is defined by a set of geometrical parameters (see Supplementary Figs. S1–S3). Our models assume that fractures deform in response to constant overpressure, and constant shear stress changes, resulting in two extra parameters. The overpressure corresponds to the difference between the magma pressure and the normal stress exerted by the host rock on the fracture surface, while the shear stress change corresponds to the shear stress exerted by the host rock. We considered constant stress changes, as all parts of magma intrusions are assumed to be hydraulically connected and depth-dependent gradients can not be constrained on the scale of the intrusions at this volcano48. Using constant stress changes on fractures leads to smoothly varying displacements32,68, similarly to observations of natural fractures71. Our models differ from those traditionally used to invert geodetic data which determine fracture displacement, rather than stress7,20. Our rationale is twofold. First, the stress boundary conditions on fractures are physically more realistic than displacement boundary conditions71. Second, stress boundary conditions provide more information. For instance, the host rock stress may be inferred from the determined stress change72.
Inversions are conducted in two stages: a search and an appraisal stage. The search stage explores the parameter space by combining forward model computation and a neighborhood inversion algorithm45,73. The appraisal stage74 is based on the Bayesian inference, and uses models computed in the search stage. A mean model, one- and two-dimensional marginal probability density functions (PPD) are computed. The one-dimensional PPDs provide the confidence intervals, while the two-dimensional PPDs provide trade-offs between parameters. More detailed explanations on the iterative procedure, subsampling, and appraisal are given in the Supplementary.
Statistical representation of the models
When inverting ground displacement, given the modeling assumptions, the result is not a unique best-fit source but a family of possible geometries and locations that equally account for the data within their uncertainties. For each inversion, we statistically define the preferential locations and the uncertainties of intrusions, based on the triangular element of fracture models. We randomly generated meshes in the 95% confidence interval of the best-fit model obtained from the appraisal stage. The family of mesh elements is used to graphically represent (1) the spatial distribution of the element centroids and their posterior probability density, which reflects the level of confidence of the presence of the intrusion at a given point in space, and (2) the dip and strike of the element normals on a stereographic projection (the poles), in a similar way to that used in structural geology. This procedure is shown in diagrams in the Supplementary Fig. S4. Pole distribution is computed using a Schmidt Distribution (which counts points within a counting circle comprising 1% of the total area of the hemisphere)75.
Interferograms used in this study are available for registered users on the CASOAR link (https://wwwobs.univ-bpclermont.fr/casoar). GNSS data used were acquired through reiteration campaigns conducted by OVPF-IPGP and are available on request.
Inverse models were conducted using the DefVolc software, which is available online on the OPGC website (http://opgc.fr/defvolc).
Auker, M. R., Sparks, R. S. J., Siebert, L., Crosweller, H. S. & Ewert, J. A statistical analysis of the global historical volcanic fatalities record. J. Appl. Volcanol. 2, 1–24 (2013).
Mitchell, N. C., Masson, D. G., Watts, A. B., Gee, M. J. & Urgeles, R. The morphology of the submarine flanks of volcanic ocean islands: A comparative study of the Canary and Hawaiian hotspot islands. J. Volcanol. Geotherm. Res. 115, 83–107 (2002).
Famin, V. & Michon, L. Volcano destabilization by magma injections in a detachment. Geology 38, 219–222 (2010).
Berthod, C. et al. Evidence of sheared sills related to flank destabilization in a basaltic volcano. Tectonophysics 674, 195–209 (2016).
Cayol, V., Dieterich, J. H., Okamura, A. T. & Miklius, A. High magma storage rates before the 1983 eruption of Kilauea, Hawaii. Science 288, 2343–2346 (2000).
Puglisi, G. et al. Dynamics of Mount Etna before, during, and after the July–August 2001 eruption inferred from GPS and differential synthetic aperture radar interferometry data. J. Geophys. Res.: Solid Earth 113, 1–20 (2008).
Montgomery-Brown, E. K. et al. Spatiotemporal evolution of dike opening and décollement slip at Kilauea Volcano, Hawai’i. J. Geophys. Res.: Solid Earth 116, 1–14 (2011).
Peltier, A. et al. Deep fluid transfer evidenced by surface deformation during the 2014–2015 unrest at Piton de la Fournaise volcano. J. Volcanol. Geotherm. Res. 321, 140–148 (2016).
Chen, Y. et al. Long-term ground displacement observations using InSAR and GNSS at Piton de la Fournaise volcano between 2009 and 2014. Remote Sens. Environ. 194, 230–247 (2017).
Dieterich, J. H. Growth and Persistence of Hawaiian Volcanic Rift Zones. J. Geophys. Res. 93, 4258–4270 (1988).
Iverson, R. M. Can magma-injection and groundwater forces cause massive landslides on Hawaiian volcanoes? J. Volcanol. Geothermal Res. 66, 295–308 (1995).
Elsworth, D. & Voight, B. Evaluation of volcano flank instability triggered by dyke intrusion. Volcano Instability Earth Other Planets 110, 45–53 (1996).
Walter, T. R. & Troll, V. R. Experiments on rift zone evolution in unstable volcanic edifices. J. Volcanol. Geotherm. Res. 127, 107–120 (2003).
Walter, T. R. et al. Rift zone reorganization through flank instability in ocean island volcanoes: An example from Tenerife, Canary Islands. Bull. Volcanol. 67, 281–291 (2005).
MacDonald, G. A. Volcanoes (Prentice Hall, 1972).
Walter, T. R., Klügel, A. & Münn, S. Gravitational spreading and formation of new rift zones on overlapping volcanoes. Terra Nova 18, 26–33 (2006).
Valentine, G. A. & Krogh, K. E. Emplacement of shallow dikes and sills beneath a small basaltic volcanic center—The role of pre-existing structure (Paiute Ridge, southern Nevada, USA). Earth Planet. Sci. Lett. 246, 217–230 (2006).
Michon, L., Saint-Ange, F., Bachelery, P., Villeneuve, N. & Staudacher, T. Role of the structural inheritance of the oceanic lithosphere in the magmato-tectonic evolution of Piton de la Fournaise volcano (La Réunion Island). J. Geophys. Res.: Solid Earth 112, 1–21 (2007).
González, P. J., Tiampo, K. F., Camacho, A. G. & Fernández, J. Shallow flank deformation at Cumbre Vieja volcano (Canary Islands): Implications on the stability of steep-sided volcano flanks at oceanic islands. Earth Planet. Sci. Lett. 297, 545–557 (2010).
Owen, S. E. & Bürgmann, R. An increment of volcano collapse: Kinematics of the 1975 Kalapana, Hawaii, earthquake. J. Volcanol. Geotherm. Res. 150, 163–185 (2006).
Lin, J.-T., Aslam, K. S., Thomas, A. M. & Melgar, D. Overlapping regions of coseismic and transient slow slip on the Hawaiian décollement. Earth Planet. Sci. Lett. 544, 116353 (2020).
Chaput, M., Pinel, V., Famin, V., Michon, L. & Froger, J.-L. Cointrusive shear displacement by sill intrusion in a detachment: A numerical approach. Geophys. Res. Lett. 41, 3307–3314 (2014).
Cayol, V. et al. Sheared sheet intrusions as mechanism for lateral flank displacement on basaltic volcanoes: Applications to Réunion Island volcanoes. J. Geophys. Res.: Solid Earth 119, 4001–4016 (2014).
Famin, V. et al. Localization of magma injections, hydrothermal alteration, and deformation in a volcanic detachment (Piton des Neiges, La Réunion). J. Geodyn. 101, 155–169 (2016).
Oehler, J. F., Lénat, J. F. & Labazuy, P. Growth and collapse of the Reunion Island volcanoes. Bull. Volcanol. 70, 717–742 (2008).
Froger, J. L. et al. Time-dependent displacements during and after the April 2007 eruption of Piton de la Fournaise, revealed by interferometric data. J. Volcanol. Geotherm. Res. 296, 55–68 (2015).
Peltier, A. et al. Long-term mass transfer at Piton de la Fournaise volcano evidenced by strain distribution derived from GNSS network. J. Geophys. Res.: Solid Earth 120, 3099–3117 (2015).
Poland, M. P., Peltier, A., Bonforte, A. & Puglisi, G. The spectrum of persistent volcanic flank instability: A review and proposed framework based on Kilauea, Piton de la Fournaise, and Etna. J. Volcanol. Geotherm. Res. 339, 63–80 (2017).
OVPF. Bulletin mensuel de l’Observatoire Volcanologique du Piton de la Fournaise. Tech. Rep. http://www.ipgp.fr/fr/ovpf/bulletin-mensuel-mardi-1-septembre-2020 (2020).
OVPF. Bulletin mensuel de l’Observatoire Volcanologique du Piton de la Fournaise. Tech. Rep. http://www.ipgp.fr/fr/ovpf/bulletin-mensuel-mardi -1-decembre-2020 (2020).
OVPF. Bulletin mensuel de l’Observatoire Volcanologique du Piton de la Fournaise. Tech. Rep. http://www.ipgp.fr/fr/ovpf/bulletin-mensuel-lundi-5-octobre-2020 (2020).
Tridon, M., Cayol, V., Froger, J. L., Augier, A. & Bachèlery, P. Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement. J. Geophys. Res.: Solid Earth 121, 7846–7866 (2016).
Cayol, V., Tridon, M., Froger, J., Augier, A. & Bachèlery, P. Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement. In IAVCEI 2017 Scientific assembly, Forstering Integrative Studies of Volcanoes, Abstr. 1209 p.180 (2017).
Rançon, J. P., Lerebour, P. & Augé, T. The Grand Brule exploration drilling: New data on the deep framework of the Piton de la Fournaise volcano. Part 1: Lithostratigraphic units and volcanostructural implications. J. Volcanol. Geotherm. Res. 36, 113–127 (1989).
Gailler, L. S. et al. Gravity structure of Piton de la Fournaise volcano and inferred mass transfer during the 2007 crisis. J. Volcanol. Geotherm. Res. 184, 31–48 (2009).
Bachèlery, P. Le Piton de la Fournaise (Ile de la Reunion). Etude volcanologique, structurale et petrologique. Piton de la Fournaise, Reunion; volcanologic, structural and petrographic study. Ph.D. thesis (1981).
Bonali, F. L., Corazzato, C. & Tibaldi, A. Identifying rift zones on volcanoes: An example from La Réunion island, Indian Ocean. Bull. Volcanol. 73, 347–366 (2011).
Michon, L., Ferrazzini, V., Di Muro, A., Villeneuve, N. & Famin, V. Rift zones and magma plumbing system of Piton de la Fournaise volcano: How do they differ from Hawaii and Etna? J. Volcanol. Geotherm. Res. 303, 112–129 (2015).
Battaglia, J., Ferrazzini, V., Staudacher, T., Aki, K. & Cheminée, J. L. Pre-eruptive migration of earthquakes at the Piton de la Fournaise volcano (Réunion Island). Geophys. J. Int. 161, 549–558 (2005).
Lengliné, O., Duputel, Z. & Ferrazzini, V. Uncovering the hidden signature of a magmatic recharge at Piton de la Fournaise volcano using small earthquakes. Geophys. Res. Lett. 43, 4255–4262 (2016).
Prôno, E., Battaglia, J., Monteiller, V., Got, J. L. & Ferrazzini, V. P-wave velocity structure of Piton de la Fournaise volcano deduced from seismic data recorded between 1996 and 1999. J. Volcanol. Geotherm. Res. 184, 49–62 (2009).
Beauducel, F. et al. WebObs : The volcano observatories missing link between research and real-time monitoring. Front. Earth Sci. 8, 1–22 (2020).
Lénat, J. F., Bachèlery, P. & Peltier, A. The interplay between collapse structures, hydrothermal systems, and magma intrusions: The case of the central area of Piton de la Fournaise volcano. Bull. Volcanol. 74, 407–421 (2012).
Peltier, A. et al. Changes in the long-term geophysical eruptive precursors at Piton de la Fournaise: Implications for the response management. Front. Earth Sci. 6, 1–10 (2018).
Fukushima, Y., Cayol, V. & Durand, P. Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise. J. Geophys. Res. B: Solid Earth 110, 1–15 (2005).
Peltier, A., Staudacher, T. & Bachèlery, P. Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data. J. Geophys. Res.: Solid Earth 112, 1–16 (2007).
Peltier, A. et al. Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion hotspot) inferred from deformation and geochemical data. Earth Planet. Sci. Lett. 270, 180–188 (2008).
Fukushima, Y., Cayol, V., Durand, P. & Massonnet, D. Evolution of magma conduits during the 1998–2000 eruptions of Piton de la Fournaise volcano, Réunion Island. J. Geophys. Res.: Solid Earth 115, B10204 (2010).
Smittarello, D. et al. Magma propagation at Piton de la Fournaise from joint inversion of InSAR and GNSS. J. Geophys. Res.: Solid Earth 124, 1361–1387 (2019).
Dumont, Q., Cayol, V. & Froger, J.-L. Mitigating bias in inversion of InSAR data resulting from radar viewing geometries. Geophys. J. Int. 227, 483–495 (2021).
Letourneur, L., Peltier, A., Staudacher, T. & Gudmundsson, A. The effects of rock heterogeneities on dyke paths and asymmetric ground deformation: The example of Piton de la Fournaise (Réunion Island). J. Volcanol. Geotherm. Res. 173, 289–302 (2008).
Got, J. L., Peltier, A., Staudacher, T., Kowalski, P. & Boissier, P. Edifice strength and magma transfer modulation at Piton de la Fournaise volcano. J. Geophys. Res.: Solid Earth 118, 5040–5057 (2013).
Richter, N. & Froger, J. L. The role of Interferometric synthetic aperture radar in detecting, mapping, monitoring, and modelling the volcanic activity of Piton de la Fournaise, La Reunion: A review. Remote Sensing 12, 1019 (2020).
Montgomery-Brown, E. K., Segall, P. & Miklius, A. Kilauea slow slip events: Identification, source inversions, and relation to seismicity. J. Geophys. Res. 114, 1–20 (2009).
Peltier, A., Staudacher, T. & Bachèlery, P. New behaviour of the Piton de La Fournaise volcano feeding system (La Réunion Island) deduced from GPS data: Influence of the 2007 Dolomieu caldera collapse. J. Volcanol. Geotherm. Res. 192, 48–56 (2010).
Peltier, A., Ferrazzini, V., Staudacher, T. & Bachèlery, P. Imaging the dynamics of dyke propagation prior to the 2000–2003 flank eruptions at Piton de la Fournaise, Reunion Island. Geophys. Res. Lett. 32, 1–5 (2005).
Chaput, M., Famin, V. & Michon, L. Sheet intrusions and deformation of Piton des Neiges, and their implication for the volcano-tectonics of La Réunion. Tectonophysics 717, 531–546 (2017).
Michon, L. et al. Edifice growth, deformation and rift zone development in basaltic setting: Insights from Piton de la Fournaise shield volcano (Réunion Island). J. Volcanol. Geotherm. Res. 184, 14–30 (2009).
Maccaferri, F., Bonafede, M. & Rivalta, E. A quantitative study of the mechanisms governing dike propagation, dike arrest, and sill formation. J. Volcanol. Geotherm. Res. 208, 39–50 (2011).
Davis, T., Bagnardi, M., Lundgren, P. & Rivalta, E. Extreme curvature of shallow magma pathways controlled by competing stresses: Insights from the 2018 Sierra Negra eruption. Geophys. Res. Lett. 48, 1–10 (2021).
Gailler, L. S. & Lénat, J. F. Internal architecture of La Réunion (Indian Ocean) inferred from geophysical data. J. Volcanol. Geotherm. Res. 221–222, 83–98 (2012).
Peltier, A. et al. Long-term mass transfer at Piton de la Fournaise volcano evidenced by strain distribution derived from GNSS network. J. Geophys. Res.: Solid Earth 120, 1874–1889 (2015).
Denlinger, R. P. & Morgan, J. K. in Instability of Hawaiian Volcanoes 148–177 (U.S. Geological Survey Professional Paper 1801, 2014).
Chaput, M., Famin, V. & Michon, L. Deformation of basaltic shield volcanoes under cointrusive stress permutations. J. Geophys. Res.-Solid Earth 119, 5384–5397 (2014).
Martínez-moreno, F. J. et al. Investigating collapse structures in oceanic islands using magnetotelluric surveys: The case of Fogo Island in Cape Verde. J. Volcanol. Geotherm. Res. 357, 152–162 (2018).
Boulesteix, T., Hildenbrand, A., Gillot, P.-y & Soler, V. Geomorphology Eruptive response of oceanic islands to giant landslides: New insights from the geomorphologic evolution of the Teide - Pico Viejo volcanic complex (Tenerife, Canary). Geomorphology 138, 61–73 (2012).
Ablay, G. & Hürlimann, M. Evolution of the north flank of Tenerife by recurrent giant landslides. J. Volcanol. Geotherm. Res. 103, 135–159 (2000).
Cayol, V. & Cornet, F. H. 3D mixed boundary elements for elastostatic deformation field analysis. Int. J. Rock. Mech. Min. 34, 275–287 (1997).
Cayol, V. & Cornet, F. H. Three-dimensional modeling of the 1983–1984 eruption at Piton de la Fournaise Volcano, Réunion Island. J. Geophys. Res. 103, 18025 (1998).
Froger, J. L. et al. The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry. Geophys. Res. Lett. 31, 1–5 (2004).
Zeller, S. S. & Pollard, D. D. Boundary conditions for rock fracture analysis using the boundary element method. J. Geophys. Res. 97, 1991–1997 (1992).
Wauthier, C., Cayol, V., Kervyn, F. & D’Oreye, N. Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an InSAR analysis. J. Geophys. Res.: Solid Earth 117, B05411 (2012).
Sambridge, M. Geophysical inversion with a neighbourhood algorithm—1. Searching a parameter space. Geophys. J. Int. 138, 479–494 (1999).
Sambridge, M. Geophysical inversion with a neighbourhood algorithm—2. Appraising the ensemble. Geophys. J. Int. 138, 727–746 (1999).
Vollmer, F. W. C program for automatic contouring of spherical orientation data using a modified Kamb method. Comput. Geosci. 21, 31–49 (1995).
Chaput, M., Famin, V. & Michon, L. Deformation of basaltic shield volcanoes under cointrusive stress permutations. J. Geophys. Res.: Solid Earth 119, 274–301 (2014).
Duputel, Z., Lengliné, O. & Ferrazzini, V. Constraining spatiotemporal characteristics of magma migration at Piton De La Fournaise volcano from pre-eruptive seismicity. Geophys. Res. Lett. 46, 119–127 (2019).
Lénat, J. F., Bachèlery, P. & Merle, O. Anatomy of Piton de la Fournaise volcano (La Réunion, Indian Ocean). Bull. Volcanol. 74, 1945–1961 (2012).
We thank the Sentinel-1 ESA team, especially P. Potin and Y.-L. Desnos for having made possible routine Sentinel-1 stripmap acquisitions on La Réunion. We are thankful to Philippe Durand (CNES) to have been a strong advocate of these stripmap acquisitions. This work was supported by the TelluS program of CNRS/INSU. This work is a contribution to the EUROVOLC project, under the EU Horizon 2020 Research and Innovation Action, grant No. 731070. This is a Laboratory of Excellence Clervolc contribution number 518. InSAR data were acquired within the framework of the Indian Ocean InSAR Observatory (OI2/OPGC/SNOV/INSU). Inversions have been performed on the supercomputer facilities of the Mésocentre Clermont-Auvergne. We are thankful to Raphaël Paris, Denis Andrault and Jean-François Moyen for proofreading this manuscript. We thank Frances van Wyk de Vries for proofreading the English of this manuscript. We thank Eleonora Rivalta who also provided useful comments.
The authors declare no competing interests.
Peer review information
Nature Communications thanks Vincent Famin, Emily Montgomery-Brown, and the anonymous reviewer(s) for their contribution to the peer review of this work.
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Dumont, Q., Cayol, V., Froger, JL. et al. 22 years of satellite imagery reveal a major destabilization structure at Piton de la Fournaise. Nat Commun 13, 2649 (2022). https://doi.org/10.1038/s41467-022-30109-w