Delphine Smittarello
email: d.smittarello@phare.normalesup.org

Valérie Cayol

Virginie Pinel

Aline Peltier

Jean-Luc Froger

Valérie Ferrazzini

J-L Froger

Magma propagation at Piton de la Fournaise from joint inversion of InSAR and GNSS

est destinée au dépôt

Introduction

Dike and sill intrusions are the dominant mechanisms for basaltic magma transport within the shallowest region of the Earth's crust. They sometimes propagate tens of kilometers away from the main magma reservoir. Reaching the Earth's surface, they often lead to fissural eruptions. Near to towns or man-made structures, fissural activity and lava flows present a risk, as experienced, for instance during the May 2018 eruption of Kilauea volcano in the Lower East Rift Zone [Global Volcanism Program, 2018] or at Piton de la Fournaise (PdF) in 1977 and 1986 [START_REF] Villeneuve | Revue de la typologie des éruptions au Piton de La Fournaise, processus et risques volcaniques associés[END_REF]. At a given volcano, intrusions often vary in style, length, orientation and direction [START_REF] Chadwick | Mechanical modeling of circumferential and radial dike intrusion on Galapagos volcanoes[END_REF]Peltier et al., 2009a]. Consequently, distances between the summit and the vents, and vents' elevations, vary from one eruption to another, as observed at PdF [START_REF] Roult | A new comprehensive classification of the Piton de la Fournaise activity spanning the 1985-2010 period. Search and analysis of short-term precursors from a broad-band seismological station[END_REF]Peltier et al., 2009a]. Moreover, magma ascent can stop before reaching the ground and erupting (see [START_REF] Moran | Failed magmatic eruptions: late-stage cessation of magma ascent[END_REF] for a review). Such non-eruptive intrusions were described, for instance at Kilauea by [START_REF] Cervelli | The 12 September 1999 upper east rift zone dike intrusion at Kilauea volcano, Hawaii[END_REF], or at Natron Rift by Calais et al. [2008]. At PdF, between 1985 and 2016, at least 15% of the detected intrusions did not lead to eruptions [START_REF] Roult | A new comprehensive classification of the Piton de la Fournaise activity spanning the 1985-2010 period. Search and analysis of short-term precursors from a broad-band seismological station[END_REF]. To mitigate the hazards induced by eruptions, it is important to understand the parameters responsible for these variations. The purpose of our work is to use in situ measurement to understand the physical and mechanical parameters controlling magma propagation.

The magma propagation toward the surface is complex and might be controlled by several physical and mechanical parameters. These parameters were studied using field mapping, analog laboratory experiments and numerical studies (see [START_REF] Rivalta | A review of mechanical models of dike propagation: Schools of thought, results and future directions[END_REF] for a review). It was proven that the shape, orientation and dynamics of intrusions are strongly influenced by the host medium (rheological properties as well as local stress field), by the magma's physical properties (density and viscosity) and by the pressurization of the feeder storage zone when the intrusion remains connected with the reservoir.

Geological [START_REF] Gudmundsson | How local stresses control magma-chamber ruptures, dyke injections, and eruptions in composite volcanoes[END_REF], analog [START_REF] Rivalta | Buoyancy-driven fracture ascent: Experiments in layered gelatine[END_REF][START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF] and numerical [START_REF] Dahm | Numerical simulations of the propagation path and the arrest of fluidfilled fractures in the Earth[END_REF][START_REF] Maccaferri | A numerical model of dyke propagation in layered elastic media[END_REF] studies evidenced the role of heterogeneities in mechanical properties on dike paths, shape, opening and kinematics. Besides, pre-existing fractures have the potential to be re-used [START_REF] Delaney | Field relations between dikes and joints: emplacement processes and paleostress analysis[END_REF][START_REF] Ruch | Oblique rift opening revealed by reoccurring magma injection in central Iceland[END_REF]. More generally, it has been demonstrated that the local stress field strongly influences magma path and velocity [START_REF] Watanabe | Analog experiments on magma-filled cracks[END_REF][START_REF] Sigmundsson | Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland[END_REF][START_REF] Pinel | A two-step model for dynamical dike propagation in two dimensions: Application to the July 2001 Etna eruption[END_REF]. For instance, the stress field is locally influenced by the topography.

When a caldera collapse occurs, this stress field becomes such that dike to sill transition is favored at shallow depth. [START_REF] Corbi | How caldera collapse shapes the shallow emplacement and transfer of magma in active volcanoes[END_REF]. Melt properties are also known to influence magma propagation. In particular, magma density controls the buoyancy force available to drive vertical propagation. Magma being a complex mix of fluid, solid particles and dissolved gas, strong density variations are potentially induced by gas exsolution or crystals sedimentation during ascent. Gas accumulation at the tip may form a pocket that favors vertical migration [START_REF] Lister | Buoyancy-driven fluid fracture: similarity solutions for the horizontal and vertical propagation of fluid-filled cracks[END_REF][START_REF] Menand | A phenomenological model for precursor volcanic eruptions[END_REF][START_REF] Taisne | Magma expansion and fragmentation in a propagating dyke[END_REF]. The magma's driving force is also determined by its initial overpressure, provided by the feeder storage zone. Studying the pressure variations of the intrusion through time provides information on the efficiency of the connection with the magma feeding reservoir [START_REF] Segall | Constraints on dike propagation from continuous GPS measurements[END_REF]. Due to the lack of information on the intrusion source, models usually consider simple end members. These include either a crack disconnected from its source [START_REF] Maccaferri | A numerical model of dyke propagation in layered elastic media[END_REF], a crack fed by a constant inflow of magma [START_REF] Pinel | Magma storage and horizontal dyke injection beneath a volcanic edifice[END_REF][START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF] or by a source at constant pressure [START_REF] Pinel | The effect of edifice load on magma ascent beneath a volcano[END_REF][START_REF] Menand | A phenomenological model for precursor volcanic eruptions[END_REF]. Such approximations may not be fully realistic. More generally, due to the multiple factors that may impact the shape and dynamics of magma during propagation, more observations, providing both high spatial and temporal resolution, are required to improve models and better determine in situ parameters and boundary conditions.

The analysis of seismic signals recorded during magma propagation might provide key additional information. The maximum shear stress is expected to be close to the tip of the propagating dike, potentially inducing micro-earthquakes [START_REF] Rubin | Dike-induced earthquakes: Theoretical considerations[END_REF]].

Thus micro-earthquakes migration can be used to track magma propagation , sometimes revealing variations in the propagation velocity [START_REF] Sigmundsson | Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland[END_REF]. However, locating the very shallow and low magnitude seismicity induced by magma intrusions using P and S-waves is challenging. In addition, even though, alternative methods [Taisne et al., 2011] exist, the crack tip propagation does not always produce seismicity. The preexisting deviatoric stress may not be large enough to generate shear failure or the magnitude of induced earthquakes may be too low to be detected by certain seismic network configurations. Besides, even when present, seismic activity does not always reveal the specific location of the dike tip, but rather a volume of maximum stress change [START_REF] Dieterich | The use of earthquake rate changes as a stress meter at Kilauea volcano[END_REF]. At PdF, micro-seismicity sometimes shows that there has been deep vertical propagation [START_REF] Battaglia | Pre-eruptive migration of earthquakes at the Piton de la Fournaise volcano (Réunion Island)[END_REF] and activation of a volume beneath the summit crater [START_REF] Lengliné | Uncovering the hidden signature of a magmatic recharge at Piton de la Fournaise volcano using small earthquakes[END_REF][START_REF] Duputel | Constraining spatio-temporal characteristics of magma migration at Piton de la Fournaise volcano from pre-eruptive seismicity[END_REF], but no lateral migration.

In addition to seismic data, high rate geodetic observations have the potential to track the magma propagation, providing additional key information on the volume involved [START_REF] Einarsson | Seismological evidence for lateral magma intrusion during the July 1978 deflation of the Krafla volcano in NE-Iceland[END_REF][START_REF] Aoki | Imaging magma transport during the 1997 seismic swarm off the Izu Peninsula, Japan[END_REF][START_REF] Segall | Constraints on dike propagation from continuous GPS measurements[END_REF][START_REF] Sigmundsson | Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland[END_REF]. [START_REF] Toutain | Real time monitoring of vertical ground deformations during eruptions at Piton de la Fournaise[END_REF] used tilt records to track the maximum uplifted area during the April 18, 1990 PdF eruption. Similarly, combining tilt data and seismic background noise, [START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF] looked into the migration of the uplifted area during 9 eruptions at PdF from 2000 to 2003, revealing a recurrent pattern : fast vertical propagation (2 m.s -1 ) followed by slower lateral migration (0.2-0.8 m.s -1 ). Nevertheless, the link between the maximum uplift and the dike tip location is not straightforward. Using high-Rate GNSS measurements at Mount Etna, [START_REF] Cannavò | Real time tracking of magmatic intrusions by means of ground deformation modeling during volcanic crises[END_REF] determined a distribution of point pressure sources, but such sources are not able to account for the complex deformation pattern produced by dike intrusions and might lead to erroneous results [START_REF] Menassian | Validation of a 'displacement tomography' inversion method for modeling sheet intrusions[END_REF]]. An alternative method to determine the evolution of volume change consists of a two-step procedure. The intrusion geometry is first determined using all available geodetic data, then linear inversions are used to retrieve the opening distribution [START_REF] Aoki | Imaging magma transport during the 1997 seismic swarm off the Izu Peninsula, Japan[END_REF][START_REF] Segall | Constraints on dike propagation from continuous GPS measurements[END_REF].

However, despite their high temporal resolution, continuous GNSS and tilt data are spatially sparse. When inverted, they lead to oversimplified source geometries. On the other hand, Interferometric Synthetic Aperture Radar (InSAR) data provides information on the displacement field, in the Line of Sight (LOS) direction, which has very good spa-tial resolution (some meters), but poor temporal resolution (usually not better than 6 days).

Thus, InSAR data can be used to image complex displacement patterns. Addressing this complexity with simple analytical models is not satisfactory and requires the combination of numerical models with inversions to reveal the complexities of geometric sources [START_REF] Trasatti | Numerical inversion of deformation caused by pressure sources: application to Mount Etna (Italy)[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF]. Because temporal resolution of InSAR is low (of the order of days) while the intrusion process is generally short (some tens of minutes to hours at PdF), the opportunity to capture an ongoing eruption is rare. One interesting example concerns an intrusion at Fernandina (Galapagos), where [START_REF] Bagnardi | A new model for the growth of basaltic shields based on deformation of Fernandina volcano, Galápagos Islands[END_REF] used an InSAR image acquired during the propagation phase to reveal the intrusion geometry at an early stage. They showed that the intrusion orientation had switched from a sill to a radial dike. In another rare case in May 2016, such an acquisition by Sentinel-1 at PdF also imaged the displacement field produced by propagating magma.

In order to overcome limitations specific to point measurements and InSAR, a combination of these two types of measurements have been proposed to track the opening of dikes [START_REF] Hamling | Geodetic observations of the ongoing Dabbahu rifting episode: new dyke intrusions in 2006 and 2007[END_REF][START_REF] Fukuda | Mixed linear-non-linear inversion of crustal deformation data: Bayesian inference of model, weighting and regularization parameters[END_REF][START_REF] Auriac | Iceland rising: Solid Earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling[END_REF][START_REF] Sigmundsson | Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland[END_REF]. The proposed methods are derived from the Okada formulation [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF] and ignore topography, and medium heterogeneities, sometimes providing erroneous source determinations [START_REF] Masterlark | Transient volcano deformation sources imaged with interferometric synthetic aperture radar: application to Seguam Island, Alaska[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]. This led [START_REF] Currenti | Toward a realistic deformation model of the 2008 magmatic intrusion at Etna from combined DInSAR and GPS observations[END_REF] to propose fully 3D heterogeneous elastic models to determine the distribution of openings on a planar vertical dike during the 2008 intrusion at Mount Etna.

However, inverting surface displacement allows either overpressure or opening distributions on fractures to be determined. It is more meaningful physically to determine the overpressure for the following reasons. Pressure boundary conditions are inherently in agreement with the physics of intrusions, where hydraulic connectivity results in homogeneous overpressure. Homogeneous overpressure leads to smoothly varying openings with no need for smoothness-constraints to regularize the inversion [START_REF] Zeller | Boundary conditions for rock fracture analysis using the boundary element method[END_REF].

Lastly, fewer parameters are required for overpressure determination than opening determinations, making these models more likely from a statistical point of view [START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF].

In this study, we propose an integrated approach which complements the high spatial resolution of InSAR with the high temporal resolution of GNSS data in order to (i) determine the complex geometry of an intrusion located beneath a realistic topography, (ii) resolve pressure changes as the intrusion propagates. We apply this strategy to the May 26, 2016 eruption at PdF, a highly active volcano characterized by a dense monitoring network (Fig. 1) and routine InSAR monitoring, providing one of the best spatio-temporal ground deformation data coverage available for a volcano. Taking advantage of the SAR image that was acquired during the magma propagation, we discuss the importance of an integrated approach to solve trade-offs between inverted parameters. Focusing on the lateral propagation phase, we determine the time evolution of the pressure source. We then discuss the characteristics of the intrusion propagation path, geometry and dynamics and their implication in terms of physical processes.

Geological and historical background

Piton de la Fournaise volcano

PdF (La Réunion Island, France) is a basaltic intra-plate volcano reaching 2632 m above sea level (a.s.l.) where fissure eruptions resulting from dike intrusions frequently occur. The Central Cone of PdF (400 m high, 1500 m radius)is located within a major Ushaped structure called the Enclos Fouqué caldera (Fig. 1). Based on the distance of the eruptive fissures from the summit, eruptions are classified as summit, proximal (when they occur on the flanks close to the summit), or distal (for distances greater than 4 km from the summit). Since 1972, the percentage of eruptions classified as each type are 21%, 68%, and 11%, respectively [Peltier et al., 2009a].

Since 1979, the Observatoire Volcanologique du Piton de la Fournaise managed by the Institut de Physique du Globe de Paris (OVPF-IPGP) has monitored the volcano to follow seismicity (40 seismometers), ground deformation (via 24 GNSS stations, 9 tiltmeters, 3 extensometers) and gas emission. The network of permanent stations is now one of the densest in the world, with distances between stations of less than 1 km in the summit area, and less than 3 km in the proximal area, (Fig. 1C) making the volcano one of the best monitored and allowing to anticipate [START_REF] Peltier | Changes in the long-term geophysical eruptive precursors at Piton de la Fournaise: Implications for the response management[END_REF] and follow eruptions with precision. In addition to in situ networks, the Indian Ocean InSAR Observatory Service (OI 2 ), carries out systematic acquisitions of InSAR data. Intermittent InSAR monitoring was performed within this framework at the volcano from 1998 , and this monitoring became systematic in 2003.

In April 2007, a distal eruption at low elevation drained the shallow magma reservoir, leading to the most voluminous lava flow in historical records, and collapse of the summit caldera associated with a 1 km × 0.8 km wide and 330 m deep crater [START_REF] Michon | April 2007 collapse of Piton de la Fournaise: a new example of caldera formation[END_REF]Peltier et al., 2009b]. A large seaward displacement (up to 1.4 m) of the east flank, evidenced by InSAR, has been associated with this major eruption [START_REF] Froger | Timedependent displacements during and after the April 2007 eruption of Piton de la Fournaise, revealed by interferometric data[END_REF][START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF]. After a series of small summit eruptions and non-eruptive intrusions in 2008-2009, the volcano remained quiet for 3.5 years. But from 2014 to September 2018, activity picked up, with 14 further eruptions which alternated between north and south flanks location (Fig. 1C) [START_REF] Peltier | Deep fluid transfer evidenced by surface deformation during the 2014-2015 unrest at Piton de la Fournaise volcano[END_REF][START_REF] Peltier | Changes in the long-term geophysical eruptive precursors at Piton de la Fournaise: Implications for the response management[END_REF][START_REF] Coppola | Shallow system rejuvenation and magma discharge trends at Piton de la Fournaise volcano (La Réunion Island)[END_REF][START_REF] Gurioli | Integrating field, textural, and geochemical monitoring to track eruption triggers and dynamics: a case study from Piton de la Fournaise[END_REF].

The May 2016 Eruption

On May 25, 2016 at 19:40 (all times are in Coordinated Universal Time (UTC)), after ten days of slight inflation and low-level seismic activity (∼5-10 Volcano-Tectonic earthquakes (VT) per days), a seismic crisis began, characterized by ∼2000 VT, with small magnitude M<1. Eight hours and 25 minutes later, on May 26 at 04:05, a 100 m long eruptive fissure opened. This proximal eruption occurred along the Southeast Rift Zone, 2800 m southeast of the summit cone, at an elevation of 1890 m. A small volume of lava, estimated around 0.5 Mm 3 from satellites and photogrammetry measurements was emplaced (Fig. 1C). This eruption lasted for 27 hours, and only a few direct observations were available due to its location and unfavorable weather conditions. Witnesses reported moderate lava fountain activity a few hours after the eruption onset.

We focus on the May 2016 eruption for two main reasons. First, a Sentinel-1 SAR image was acquired on May 26 at 01:45 UTC, i.e. during the seismic crisis that preceded the eruption onset. By combining this image with the previous and the following Sentinel-1 acquisitions for PdF, we were able to compute two interferograms. The first interferogram (hereafter refered to as S1 D1 for Sentinel-1 Descending 1) imaged the displacement related to the earliest part of the magma propagation. The second interferogram (hereafter refered to as S1 D2 for Sentinel-1 Descending 2) imaged the displacement related to the shallowest part of the magma propagation to the surface and to the opening of the eruptive fissure. Such data provide unique constraints that may provide insight into the magma propagation. Second, considering the location of the vent, the seismic crisis preceding the eruption was notably longer than usual. Based on a compilation of distal and proximal eruptions that occurred at PdF from 1985 to 2016, seismic crisis generally last from 20 minutes to 18 hours (Fig. 2A). The longer the seismic crises, the further it is from the summit and the lower in elevation the fissure location [START_REF] Aki | Seismic monitoring and modeling of an active volcano for prediction[END_REF][START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF]. A linear trend can be plotted which corresponds to a mean velocity for lateral propagation of 0.19 m.s -1 , (Fig. 2A). For the May 2016 eruption, considering the summit-vent distance (2800 m), the seismic crisis was expected to last ∼3h30min, which is significantly shorter than the observed 8h25min duration. Another central question concerns the reason why the magma stalled for 5 hours before erupting. If this delay was due to a slower than usual propagating intrusion then, according to [START_REF] Corbi | Understanding the link between circumferential dikes and eruptive fissures around calderas based on numerical and analog models[END_REF] it should have been more sensitive to the stress field. We postulated that this would make the study of this intrusion particularly informative.

Available observations for the May 2016 Eruption

Final ground surface displacement distribution

Thanks to the routine InSAR monitoring at PdF, the May 2016 eruption was imaged by acquisitions from the Cosmo-Skymed (CSK, X-band) and the Sentinel-1 (S1, C-band) satellites in Stripmap mode along ascending and descending tracks (Figs. 3 and4, Table S1). The computed interferograms are high spatial resolution (5 m by 5 m pixels).

The main processing steps are described in the Supporting Informations Text S1 [START_REF] Cnes | Philosophy and Instructions for the Use of the DIAPASON Interferometry Software System Developed at CNES[END_REF][START_REF] Massonnet | Radar interferometry and its application to changes in the Earth's surface[END_REF][START_REF] Chen | Network approaches to two-dimensional phase unwrapping: intractability and two new algorithms[END_REF].

Descending interferograms (CSK D and S1D) show a single lobe centered to the southeast of the Central Cone, with maximum Line Of Sight (LOS) displacements of about 40 cm and 35 cm toward the Sentinel-1 and CSK satellites, respectively (Fig. 4).

Ascending interferograms (CSK A and S1A) show patterns with two lobes. The western lobe has maximum LOS displacements of 17 cm and 12 cm towards the Sentinel-1 and CSK satellites, respectively, while the eastern lobe has displacements of up to 10 cm and 13 cm away from the Sentinel-1 and CSK satellites, respectively. The small differences between Sentinel-1 and CSK interferograms are well explained by the differences in LOS (see Table S1).

Based on deformation patterns (https://wwwobs.univ-bpclermont.fr/casoar), the May 2016 eruption shares strong similarities with 5 eruptions which have taken place since 1998 (Fig. S1), that in June 2000, October 2010, May 2015, January and July 2017. For the June 2000 eruption, vertical ascent beneath the Dolomieu crater followed by lateral migration of the intrusion to the surface was determined from tilt and seismic data by [START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF]. Ascending and descending InSAR data which captured this eruption were inverted by [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF], showing that the intrusion responsible for the displacement was a curved dike inclined to the northeast. All the others occurred after the major April 2007 collapse and no study of their deformation sources has been published yet.

Ten continuous GNSS stations (cGNSS) of the OVPF network (Fig. 1) recorded displacements due to the May 2016 magma propagation. The daily position of each station is calculated using the GAMIT/GLOBK software package [START_REF] Herring | Introduction to Gamit/Globk[END_REF], providing solutions with standard deviations of 0.5 cm and 1 cm for horizontal and vertical components, respectively. GNSS data processing is described in the Supporting Information Text S2 [START_REF] Herring | Introduction to Gamit/Globk[END_REF][START_REF] Nikolaidis | Seismic wave observations with the Global Positioning System[END_REF]. Displacements between August 28, 2015 and May 30, 2016 have also been measured in static-rapid mode at 80 locations of the campaign GNSS network, on the central cone. Comparing cGNSS data with data from this campaign GNSS network (Fig. 5), we find that it provides less precise solutions. For example, horizontal and vertical uncertainties are around 1-2 cm and 4-5cm respectively. However, it improves the spatial resolution of GNSS measurements. Records from daily GNSS measurements show that the displacements recorded during the inter-eruptive period were low with respect to the co-eruptive displacements (3-5 cm versus tens of cm) and can be neglected. Consequently, we assume that the deformation measured with the campaign GNSS and detected by InSAR is mainly due to the May 25 magma propagation.

The cGNSS stations which recorded the largest displacement are located at the southeast border of the Dolomieu crater (DERG 21.2±2.2 cm, DSRG 18.5±2.3 cm, SNEG 13.5±2.2 cm) or at the eastern base of the cone (FERG, 21.2±3.3 cm) (see Table S2). Because the campaign GNSS network has more stations on the volcano than the continuous network, and due to its different spatial distribution, it recorded larger displacements, of up to 24.5±5.5 cm to the east of Dolomieu summit crater.

Comparing InSAR data with the campaign GNSS and the cGNSS data, both projected along the LOS, we assess the consistency between these data. On Figure . 4, this consistency is illustrated by the similarity between the background color, representing In-SAR displacement, and colors filling the squares and circles representing the GNSS continuous and campaign measurements, respectively, with the difference being less than the standard deviation on measurements. Such consistency was also established by [START_REF] Peltier | Assessing the reliability and consistency of InSAR and GNSS data for retrieving 3D-displacement rapid changes, the example of the 2015 Piton de la Fournaise eruptions[END_REF] for the 4 eruptions of 2015 at PdF. However, even though both GNSS networks recorded significant deformation, the stations are far from the area of maximum deformation measured by InSAR. The campaign GNSS network is centered on the Central Cone while the deformation is the greatest in the southeast area. Therefore, we postulate that since GNSS and InSAR data are consistent and campaign GNSS have large uncertainties, GNSS data do not provide any useful constraints on the final location of the intrusion (see section 5.1). In the following, to simplify, we do not include GNSS data when modeling the final geometry of the intrusion.

Temporal observations

From the temporal evolution of the tilt, seismicity and Real-time Seismic Amplitude Measurement (RSAM) ratio between DSM (summit) and FOR (southeast flank) stations, the unusually long May 2016 pre-eruptive seismic crisis can be divided to four phases (Fig 2B). It started at 19:40 UTC (phase Ia) with an intense seismic swarm, mostly affecting the summit as shown by RSAM ratio increase, during which about 430 microearthquakes occurred. We managed to locate wenty three of these events which were all beneath the summit. Thirty minutes later (i.e. at 20:10, phase Ib), tilt inside the Enclos Fouqué accelerated (Fig. 2B), showing an inflation of the summit crater. From 20:40 (phase II), tilt direction at a station on the southeastern border of the Dolomieu crater (DSOi) abruptly changed while seismicity rates decreased up to 22:00 UTC with, 2 out of 370 micro-earthquakes located. Seismicity then stabilized (phase III) at a lower, but still above baseline, level of approximately 100 earthquakes per hour for about 4 hours, until the onset of a second intense seismicity swarm from 02:15 to 03:15 UTC on May 26 (phase IVa). The decrease of the seismicity rate (phase IVb) coincided with a strong decrease in RSAM ratio between DSM and FOR stations, in agreement with the location in seismicity on the southeast flank (9 out of 600 events recorded were located). The eruption started at 04:05. During phases III and IV, from 22:00 to the eruption onset, tilt signal variations decreased rapidly and no significant signal was recorded.

A seismic swarm such as observed in phase I, followed by a decrease of seismic activity and tilt rotation as seen in phase II is often observed at PdF [START_REF] Roult | A new comprehensive classification of the Piton de la Fournaise activity spanning the 1985-2010 period. Search and analysis of short-term precursors from a broad-band seismological station[END_REF].

The classic interpretation is that magma firstly propagates vertically beneath the summit and then laterally up to the eruption onset or the end of the intrusion [START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF]. However, phase III and especially phase IV are unusual for this volcano. Another uncommon point is the occurrence of many Long Period events (LP) recorded during phases II (31 events) and III (18 events). Such events at PdF are generally associated with lateral propagation [START_REF] Aki | Seismic monitoring and modeling of an active volcano for prediction[END_REF]. Tilt data provide valuable insights into the very beginning of the propagation (phases I and II), but because the tiltmeters are centered on the summit, this dataset is then insensitive to changes that may have occurred during lateral propagation. By chance, an additional Sentinel-1 image along a descending track was acquired during phase III at 01:45 on 26 May, allowing two interferograms to be produced covering the beginning (S1 D1) and the end (S1 D2) of the propagation (Figs. 3, 4E-F, S2 and Table S1). Comparing S1 D1, which recorded the first part of the intrusion propagation, with the Sentinel-1 descending (S1 D) interferogram covering the whole propagation (Fig. 4D andE andS2), shows that most displacement had already occurred before 01:45 : the maximum LOS displacement of S1 D1 (38 cm) achieved almost the maximum displacement of S1 D (40 cm). S1 D2, which recorded the final part of the propagation (Fig. 4F andS2) only shows a small asymmetric displacement reaching ∼20 cm toward the satellite, close to the eruptive fissure, and a large wavelength phase shift which corresponds to a maximum of ∼5 cm toward the satellite east of the summit cone between May 26 at 01:45 and June 7.

However, such a displacement is not shown by the cGNSS measurements, for which only 1 cm of eastward displacement was recorded, along with negligible northward displacement and subsidence. This suggests that at least part of the phase shift revealed by InSAR could be of atmospheric origin. cGNSS data are also processed with the TRACK software providing epoch by epoch solutions (i.e. one position every 30s). The GITg station, located 4 km from the summit, is used as a reference. In order to increase the signal to noise ratio, we corrected the timeseries from multi-paths [START_REF] Nikolaidis | Seismic wave observations with the Global Positioning System[END_REF] and performed a Principle Component Analysis (PCA). The standard deviation on the post-processing cGNSS data was computed for horizontal (1-2 cm) and vertical (3-4 cm) components (see details in the Supporting Information Text S2). Both co-eruptive displacements measured with daily solutions and post-processed cGNSS, taking the reference on May 25, at 17:00 UTC, are in very good agreement. Time evolution of the cGNSS displacements (Figs. 5 andS3) confirms that most of the deformation occurred at the beginning of the crisis. This timing is consistent with the four phases previously described based on the time evolution of seismicity, RSAM and tilt. From 20:10 (phase Ib), summit stations (SNEg, DERg, DSRg) moved in a centrifugal pattern. From 20:25, stations further to the east (FERg) and south (FOAg) started moving northeastward and southwestward, respectively. At that time, displacements were mainly horizontal, with amplitudes of the order of their standard deviation. This result is in agreement with the classic pattern for the volcano where tilt and seismic data are the first to record unrest. From 20:40 (phase II), all summit displacement vectors turned to the northwest then stopped between 21:00 and 21:15, while the stations at the base of the cone, FERg and FOAg, were still moving slowly to the east and south, respectively. Coincident with this change in direction, the stations started moving vertically. DSRg reached an amplitude of 12 cm, greater than the horizontal displacements. No significant deformation was recorded after 22:00 (phases III and IV) by the cGNSS network. This is consistent with InSAR data, which showed that the final part of the eruption (S1 D2) took place in an area lacking cGNSS monitoring (Fig. 4F andS2). Despite its density, the permanent network is too sparse to comprehensively sample the displacement field produced during the May 2016 eruption. This observation supports the idea that an integrated approach is required to fully interpret the ground deformation time series in term of magma propagation dynamics.

Inverse Modeling Methods and Strategies

In order to determine the source of ground deformation measured by InSAR and GNSS, we use a 3D-Mixed Boundary Element forward model [START_REF] Cayol | Analyse élastostatique tridimensionnelle du champ de déformations des édifices volcaniques par éléments frontières mixtes[END_REF]Cornet, 1997, 1998a] combined with a neighborhood inversion algorithm [Sambridge, 1999a;[START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]. We first invert the 2 descending (CSK D and S1 D) and 2 ascending (CSK A and S1 A) interferograms covering the whole eruptive crisis (see Figure 3 and Table S1) to determine the final shape of the magma intrusion (Fig. 6A). In most cases at PdF, once the eruptive dike is emplaced, even if an eruption is ongoing, no significant deformation is observed, probably indicating that the conduit used to supply lava remains the same [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF]. Because the cGNSS time series evidence no significant deformation except during the 8 hours preceding the onset of the eruption, we consider that deformation for the whole eruption is the same as deformation resulting from the propagation of the intrusion to the surface. We then track lateral magma propagation by inverting displacements recorded by the cGNSS network for progressively longer durations, such that the determined area corresponds to the cumulative pressurized part of the intrusion.

Our approach, hereafter called the Projected Disk method, consists of using the mesh determined by the inversion of the four interferograms spanning the whole eruptive crisis as an a priori for the source location and shape, inverting for a circular area projected on this mesh (Fig. 6B).

Forward modeling

Ground displacement for linear elastic homogeneous media is computed using a 3-D Mixed Boundary Element Method [START_REF] Cayol | 3D mixed boundary elements for elastostatic deformation field analysis[END_REF]. This method combines a direct method with a displacement discontinuity method to allow for fast and precise computations. The method considers tensile cracks and shear fractures as well as realistic surface topographies, as neglecting topography was shown to induce errors in source depth and volume characterization [Cayol and Cornet, 1998b;[START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Masterlark | Magma intrusion and deformation predictions: Sensitivities to the Mogi assumptions[END_REF]. The edifice is assumed to have a Young's Modulus of 5 GPa and a Poisson's ratio of 0.25, as proposed by [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF] for PdF. Boundaries are discretized with triangular elements. The topographic mesh is constructed from a Digital Elevation Model that includes the bathymetry. The mesh extension is about five times greater than the deformed area in order to limit edge effects to a few percent [START_REF] Cayol | Analyse élastostatique tridimensionnelle du champ de déformations des édifices volcaniques par éléments frontières mixtes[END_REF]. To provide the best compromise between computation time and accuracy, the mesh is dense close to the eruptive fissure (50 m between nodes) and progressively coarser further away (Fig. S4).

To determine the geometry of the displacement source, only 9 geometrical parameters (Table 1 and Fig. 6A) are required. Because of the reduced number of parameters required compared to kinematic inversions for instance, this method is statistically more significant [START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF]. The location of the eruptive fissure is known from field observations and used as an input to the model. Because the May 2016 eruption has one eruptive fissure, we search for a quadrangle-shaped intrusion linked to this fissure by a single echelon. This quadrangle could be curved along the strike and the dip direction.

Preliminary models show that, as for the June 2000 eruption [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF], such a curvature significantly increases the data model agreement. To follow the propagation of the deformation source before it reaches the surface, we use the Projected Disk method which only requires 3 geometrical parameters as an a priori on the geometry is provided of the curved quadrangular fracture, previously determined from InSAR data inversion (Table 1 and Fig. 6B). We define the mean plane for the a priori curved quadrangular fracture. On this plane a circular source is defined which is projected to the curved quadrangular fracture. Parameters for the inversion of this source are the radius and location of the source before projection. This method allows us to determine a source consistent with both InSAR and cGNSS data , reducing the number of inverted parameters and fixing the depth by the use of an a priori mesh. Because the disk can be partially projected beyond the mesh, we interpret the area change and not the radius change.

Prescribed boundary conditions are pressures changes. The ground surface is assumed to be stress-free and intrusions are assumed to open as a response to the difference between the magma pressure and the normal stress exerted by the host rock to their surface. We assume these stress changes are homogeneous. Linearly varying overpressure could also be assumed, but a previous study [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF] showed that overpressure gradients could not be constrained at PdF.

Non linear Inversions

Following, [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF] and [START_REF] Wauthier | Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an InSAR analysis[END_REF][START_REF] Wauthier | Nyamulagira's magma plumbing system inferred from 15 years of InSAR[END_REF], we use a Monte Carlo Neighborhood Algorithm to invert the deformation data. The inversion process involves determining the forward models which best fit the data. The fit between modeled and observed data is expressed through the misfit function, χ 2 , calculated in the least square sense for each forward model as:

χ 2 = (d o -G(m)) T C d -1 (d o -G(m)) (1)
where d o represents the vector of observed displacements, G(m) is the vector of modeled displacements such that m is the set of model parameters and G represents the model prediction. C d is the covariance matrix on the data. It accounts for correlated noise related to the atmospheric contribution to InSAR data, and for the variance of cGNSS data, as well as for the modeling errors. For explanations about construction of the cGNSS and InSAR data vector, as well as the covariance matrix, refer to Supporting Information Text S3 [START_REF] Tarantola | Inverse problem theory: Method for data fitting and model parameter estimation[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Welstead | Fractal and wavelet image compression techniques[END_REF][START_REF] Jónsson | Fault slip distribution of the[END_REF][START_REF] Sudhaus | Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: application to the June 2000 Kleifarvatn earthquake, Iceland[END_REF][START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF] and Figure S5.

We also calculate the % explained data for the i th interferogram (%E d i ) as:

%E d i = 1 - u obs i -u mod i T u obs i -u mod i u obs i T u obs i * 100 (2)
where u i obs and u i mod are the data and model vectors, respectively, corresponding to all pixels in the masked interferograms. This value is independent of the covariance matrix, the subsampling and the dataset used to run the inversion.

The inversion consists of two stages: a search and an appraisal stage. The search stage [Sambridge, 1999a] aims to find an ensemble of models that preferentially sample the good data-fitting regions of the parameter space. The appraisal stage [START_REF] Sambridge | Geophysical inversion with a neighbourhood algorithm -II. Appraising the ensemble[END_REF] concerns the estimation of model uncertainties. Following the Bayesian inference, the population of models calculated during the first stage is resampled with a Monte Carlo integration procedure, enabling the calculation of the posterior probability density function (PPD) without any new calculation of the forward problem (see details in the Supporting Information Text S4 [Sambridge, 1999a[START_REF] Sambridge | Exploring multidimensional landscapes without a map[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Wauthier | Magma sources involved in the 2002 Nyiragongo eruption, as inferred from an InSAR analysis[END_REF]).

Strategy for the inversion of the GNSS time series

The intermediate S1 D1 interferogram provides good spatial-resolution information on the displacement with a good spatial resolution, while the intrusion was underway. We used this data independent of the GNSS time series to validate our Projected Disk method. We invert the shape and location of the intrusion using either GNSS, InSAR (S1 D1) or both datasets together. We find that introducing a priori on the source location is needed to be able to invert cGNSS data and that the Projected Disk method leads to consistent results between InSAR and cGNSS data even when one datatype is omitted (Table 2). We next perform independent inversions considering progressively longer durations of the cGNSS recording starting at t 0 =17:00 UTC. The first time step in the inversion corresponds to the moment when the cGNSS signal surpasses the standard deviation, i.e. at 20:40. Length of the time steps is adjusted to the displacement change, using smaller time steps when larger changes are observed. We invert sources for 29 successive time intervals, with 5 min time steps between 20:40 and 21:30, then 10 min time steps up to 22:30 and finally 30 min time steps until the eruption onset at 04:05. As the intrusion proceeds, changes in the GNSS signal become smaller (Fig. 5). One more calculation at 01:45 was made for comparison with InSAR data (S1 D1) which lead to 30 time steps computed.

Results

Final geometry of the intrusion derived from four interferograms

Inverting the 4 InSAR data spanning the whole eruption, we find that the best-fit final geometry for the intrusion is a low dipping intrusion at a mean depth of 800±60 m a.s.l., slightly dipping to the northeast (∼8 °). We refer to this intrusion as a sill for two reasons: i) the dip is of the order of the topographic slope above the intrusion; ii) in a similar basaltic shield volcano context at Fernandina, Galapagos, such a low dipping intrusion was also referred to as a sill. This sill is rooted to the southeast of the Dolomieu crater, and covers a quadrangular area of about 2700 m × 600 m. It reaches the surface turning into a sub-vertical 880±16 m high dike (Tables 2 andS3, Fig. 7) feeding the eruptive fissure. A plot of the opening shows that the northern corner of the sill did not really open (Fig. 7A) and that the intrusion shape is probably less angular at depth than predicted by the model. Tests of joint inversions conducted using GNSS data in addition to

InSAR, even if their weight is increased to be the same as InSAR data, do not modify the obtained model or the data fit (see Supporting information S5 [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF][START_REF] Duputel | Accounting for prediction uncertainty when inferring subsurface fault slip[END_REF][START_REF] Auriac | Iceland rising: Solid Earth response to ice retreat inferred from satellite radar interferometry and visocelastic modeling[END_REF][START_REF] Yabuki | Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip[END_REF][START_REF] Akaike | Seasonal adjustment by a Bayesian modeling[END_REF][START_REF] Fukuda | Mixed linear-non-linear inversion of crustal deformation data: Bayesian inference of model, weighting and regularization parameters[END_REF][START_REF] Sigmundsson | Segmented lateral dyke growth in a rifting event at Bárðarbunga volcanic system, Iceland[END_REF]) confirming that the addition of GNSS data is not necessary. Our best-fit model well explains the 4 interferograms (Table 2), with a better fit (95%) for the descending tracks than for the ascending ones (82%) and an underestimation of the eastward displacement east of the summit cone (Fig. 8). The appraisal (Fig. 7B) does not evidence any trade-off between parameters. We estimate the volume of this intrusion to be 2.5 Mm 3 with an average opening of 0.5 m (Table 2). Consequently, 2.5 Mm 3 of the magma was stalled at depth, while 0.5 Mm 3 of lava was emitted. Neglecting volume changes due to compressibility and degassing, an approximation of the total volume supplied to the volcano leads to 3 Mm 3 . We estimate that 80% of the magma remains trapped at depth, which is unusual at this volcano, where most of the magma is usually emitted. Inversions of InSAR displacements for 5 eruptions from 1998 to 2000, and inversions of campaign GNSS data for 6 eruptions between 2004 and 2006, [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion hotspot) inferred from deformation and geochemical data[END_REF] led to estimate that 17% and 15 %, respectively, of the magma supplied remained to contribute to the intrusive growth of the edifice.

Temporal evolution

Initial upward propagation from tilt and seismicity

From 20:10 to 20:40, vertical cGNSS displacements were four times smaller than their standard deviation, and only recorded at the 4 summit stations. As a consequence, the source can not be constrained from these data. However, tilt and seismicity [START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF][START_REF] Battaglia | Pre-eruptive migration of earthquakes at the Piton de la Fournaise volcano (Réunion Island)[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF], which usually record unrest earlier than the cGNSS data, suggest that PdF intrusions first propagate upward below the summit crater then laterally. In coherence with these authors, we will assume that the first minutes of the crisis correspond to an upward dike propagation (phase Ib). Assuming the reservoir is around sea level [START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF][Peltier et al., , 2009a;;Lengliné et al., 2016, e.g.] and the intrusion stops propagating vertically when it reaches the lowest point in the identified sill (600 m a.s.l.), we can compute a vertical propagation velocity of 0.3 m.s -1 . As noted before for five eruptions that took place between 1998 and 2000 [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF], InSAR data show no evidence of the vertical conduit connecting the reservoir and the modeled intrusion. This can be explained either by the conduit being too narrow to be distinguished from the signal of the intruded sill, or by the closure of this conduit as magma propagated upward, consistent with the propagation of a buoyant magma.

Beginning of the propagation from the first intermediate InSAR data : a validation of the Projected Disk method

We take advantage of the S1 D1 interferogram computed using a Sentinel-1 image acquired during the propagation, 2h20 before the eruption onset, to validate our Projected

Disk method and to gain information on the position of the intruded magma at this stage.

We perform 3 inversions, first using only cGNSS data, second using only the S1 D1 In-SAR data and last inverting both datasets jointly. We postulate that the robustness of the method is reflected by its ability to retrieve data omitted in the inversions. Thus we also compute the % of explained data on the data (InSAR or cGNSS) omitted in the inversion.

We find that the Projected Disk method works well for explaining data (either InSAR or cGNSS) omitted in the inversion (Table 2). For all dataset used, the horizontal location of the fracture is well retrieved and parameters are close (Fig. 9). When only cGNSS data are used, the source overpressure and volume are overestimated as 60% and 35%, respectively, and the source area is underestimated by 10% compared to the joint inversion of cGNSS and InSAR data. This is probably related to the spatial distribution of cGNSS data with respect to the intrusion. Inverting only InSAR data leads to equally low misfits for cGNSS or InSAR data as for jointly inverted cGNSS and InSAR data (Table 2), confirming that when InSAR data are available, adding the cGNSS data does not improve the fit (Supporting information Text S5). Finally, joint inversion of these two datasets leads to the shortest confidence intervals (Table . S4). However, using only the cGNSS data, the fit of the omitted InSAR data remains acceptable (77% Table 2). This, together with the ability of the method to determine the model parameters (Fig. S6), indicates that this method is a satisfactory way of tracking magma from cGNSS data.

The inversion shows that, at 01:45, 77% of the area had already opened and that this opening mainly concerned the sill part of the intrusion. This sill had an average opening of 1 m, while the dike part of the intrusion remained closed. We estimate the volume of magma in the intrusion to be 2.6 Mm 3 (Table 2) which is consistent but slightly larger than the volume of magma (2.5 Mm 3 ) for the whole intrusion (Table 2).

Insights into the magma propagation dynamics from cGNSS data

At 20:40, the intrusion started as a sill with a 0.5 Mm 2 area, located southeast of the Dolomieu Crater, at 800 m a.s.l. (Figs 10, S7 and movie S1 in the Supporting Information). The source area slowly increased until 21:10. At 21:15, a sudden southward propagation, with an area increase of up to 3.8 Mm 2 was associated with a sudden decrease in pressure. At 21:30, the intrusion had reached its maximal extension and the area became stable. We estimate the mean velocity of the propagation front as being 0.6 m.s -1 between 20:40 and 22:00, with a peak of 2 m.s -1 around 21:15. From 21:30, while the area remained unchanged, pressure and volume slowly increased up to midnight. No significant evolution was evidenced after that time. In particular, neither the opening of the vertical fracture nor the opening of the eruptive fissure were recorded on the cGNSS network as stations are too far from the most deformed area. However, the interferogram S1 D2 covering the end of the eruption reveals that deformation was focused close to the location where the fissure later reached the surface.

Differentiating volume with respect to time provides insight into the magma inflow evolution (Fig. 10B). After a brief increase from 400 ±100 m 3 .s -1 to 900 ±100 m 3 .s -1 from 20:40 to 21:15, the flow rate strongly decreased up to 22:00 and stabilized at 30 ±10 m 3 .s -1 , without any significant variations until the eruption onset at 04:05.

End of the propagation from the second intermediate InSAR data

Inverting the S1 D2 interferogram spanning the end of the propagation with the Projected Disk method, the best model is a 880 m high subvertical dike with an average opening of 0.4 m (Fig. 11A, Table S5), which feeds the eruptive fissure. We determined a dike volume of 0.3 Mm 3 (Table 2). The source covering the beginning (blue sill) and the end (red dike) of the eruption cover complementary areas (Fig. 11) and are consistent with the seismicity that could be located (Fig. 7A). Consistency of both sources is also shown by the fact that the interferograms covering the whole eruption are well explained by the sum of the models (Fig. 11, last row) corresponding to the first part (Fig. 11, first column) and the second part (Fig. 11, last column) of the eruption derived from the S1D1 and S1D2 InSAR data, respectively. Our best model (Fig. 11B and Table 2) explains 83% of the InSAR data and is also consistent with a model derived using a quadrangular source (Supporting information Text S6).

We assume that the decrease in seismicity after the second major peak at 03:15 (Fig. 2B, phase IVb) coincided with the onset of vertical propagation, for the following reasons: (1) the second peak occurred at a time intermediate between the sill propagation and the final dike propagation as indicated by the intermediate SAR acquisition, (2) at PdF earthquake swarms are often associated with the onset of propagation, (3) laboratory hydraulic fracturing experiments show that the onset of acoustic emission (the equivalent of our in situ micro-earthquakes) is an indication of the initiation of fracturing and that acoustic emissions drop as soon as the fracture starts propagating [START_REF] Zoback | Laboratory hydraulic fracturing experiments in intact and pre-fractured rock[END_REF].

Based on this assumption, we determine a vertical propagation velocity of 0.3 m.s -1 corresponding to 880 m of propagation in 50 min.

Discussion

Modeling limitations

Forward modeling used in this study relies on linear elasticity, and a homogeneous and isotropic medium. Solution of the inverse problem solved here is non unique, such that one displacement field could be explained by combinations of different source geometries or by numerous mechanical properties, for instance, inclusion of anisotropy or plasticity. Because no inclined dikes are visible on the walls of the PdF summit crater, [START_REF] Got | Edifice strength and magma transfer modulation at Piton de la Fournaise volcano[END_REF] assumed that the observed asymmetric displacements could be attributed to vertical dikes emplaced in a plastic medium. However their models are only 2D. Also, most intrusions found in the plumbing system of the dormant neighbor Piton des Neiges volcano are inclined outward, which is consistent with a brittle edifice [START_REF] Chaput | Sheet intrusions and deformation of Piton des Neiges, and their implication for the volcano-tectonics of La Réunion[END_REF] and dikes determined from the inversion of Piton de la Fournaise intrusions [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion hotspot) inferred from deformation and geochemical data[END_REF]. Thus, an elastic and homogeneous medium explains, at least, the first order displacement field. However, we acknowledge that, at the dike tip, plastic behavior probably occurs due to stress singularities [START_REF] Jaeger | Fundamentals of rock mechanics[END_REF]. Moreover, for the May 2016 eruption, residuals east of summit cone have a different sign for the descending and ascending interferograms (Figs. 7B and11B). This is consistent with an eastward residual displacement of the flank that our model is not able to explain. Considering the short time-scale studied and the fact that an elastic framework seems sufficient to explain most of the observations with a minimum number of parameters, this hypothesis is the most likely from a statistical point of view. However, we also note the need for models able to properly take into account more complex mechanical properties.

Integration of InSAR and cGNSS data

Despite the fact that PdF has one of the densest cGNSS networks in the world, this study shows some limitations of using only this cGNSS network when attempting to track proximal or distal intrusions. It must be noted that large areas are unmonitored to the north and southeast of the Central Cone where, in the past decades, numerous fissures opened and lava flows were emplaced [START_REF] Roult | A new comprehensive classification of the Piton de la Fournaise activity spanning the 1985-2010 period. Search and analysis of short-term precursors from a broad-band seismological station[END_REF]. Because of this configuration, cGNSS measurements during the May 2016 proximal eruption missed the maximum deformation and provided little information on the deformation gradient. Another limitation of the summit coverage of the cGNSS network is that it is often blind to shallow distant processes, such as the final opening of the eruptive fissures. While, deformation associated with this opening is large (tens of centimeters) in amplitude but it only affects a small area (Fig. 4F). In order to combine the advantages of the high temporal resolution of cGNSS data and the good spatial resolution of InSAR data, we developed an approach in which the temporal evolution of the source is tracked by inverting for a circular pressure source projected on the a-priori source determined from InSAR data covering the whole eruption (Projected Disk method). We also tried another method which was free of such a priori. This method consisted of inverting the migrating soure as a planar ellipse (see Supporting Information Text S7, Fig. S8 and Table S6). Both methods were tested for their ability to determine the source characteristics from cGNSS data alone. Taking advantage of the Sentinel 1 image acquired when the intrusion was underway, we compared models from inversion of cGNSS data only with models from joint inversion of InSAR and cGNSS for the same time step. This validation test revealed that the Ellipse method was unable to determine sources consistent with InSAR data from cGNSS data only (Figs. S9 and S10, Table S7 andS8). However, despite the cGNSS network configuration not being optimal for this proximal eruption, using the Projected Disk method, the inversion of cGNSS data led to source parameters consistent with inversions on InSAR data and joint datasets and the fit of InSAR data omitted in the inversion is acceptable (77%) (Figs. 9 and S9). The mains limitation concerns the source pressure (60%) and volume (35%) with respect to those estimated using both cGNSS and InSAR data.

Monitoring with high spatial and temporal resolution analogue experiments of dikes and cone sheet propagation, [START_REF] Guldstrand | Experimental Constraints on Forecasting the Location of Volcanic Eruptions from Pre-eruptive Surface Deforma-tion[END_REF] proposed a method which used the detection of the maximum uplift to predict the vent location. However, even at some of the best monitored volcanoes, this approach is far from being applicable : (i) Propagation duration is usually too short to be detected by InSAR monitoring. Moreover, the processing of InSAR data is not instantaneous. (ii) cGNSS networks would need to be very dense. Here, with one of the densest networks worldwide, part of the propagation of proximal or distal eruptions is still undetected by the network. (iii) the maximum uplift does not always provide information on the vent location. The May 2016 eruption demonstrates that fissure eruptions sometimes start as subhorizontal sills. In this case, the maximum uplift is roughly central to the source rather than being located at the intersection between the fracture and the ground.

Origin of the step-wise propagation

The May 2016 intrusion is characterized by two changes in direction (vertical to horizontal, then horizontal to vertical), as well as a step-wise horizontal propagation, with one acceleration (21:15) and a pause before the final vertical propagation. Changes in the propagation velocity can be attributed to different propagation directions (vertical / lateral / vertical). Other changes, occurring while magma is propagating laterally along a planar intrusion, probably have different origins. We investigate here the origin of the dike to sill rotation, whether observed changes in lateral velocity are real and, if so, their origin, the possibility that the intrusion was fed by a single batch of magma quickly disconnected from the reservoir, and finally the potential role of gas exsolution and east flank slip.

Required conditions for dike to sill rotation

A sudden change from vertical to lateral propagation is commonly observed at PdF [START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF]. In the May 2016 case, changes in intrusion directions are an interesting feature. Here we investigate whether this dike to sill rotation was a response to a rotation of the minimum principal stress [START_REF] Menand | Dyke propagation and sill formation in a compressive tectonic environment[END_REF][START_REF] Maccaferri | A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation[END_REF] or whether its geometry is influenced by a preexisting discontinuity.

Sills can alternate with dikes, due to stress rotation induced by successive intrusions.

If stresses in the edifice are close to isotropic, as assumed for the dormant neighbor volcano Piton des Neiges [Chaput et al., 2014a] or for Fernandina in the Galapagos [START_REF] Chadwick | Mechanical modeling of circumferential and radial dike intrusion on Galapagos volcanoes[END_REF], dike intrusions would provide enough stress change to induce a rotation of the minimum principal stress, further inducing sill intrusions. A relaxation mechanism such as slip on the east flank would then be required to start initiate dike intrusions again. While relaxation of the east flank was demonstrated in 2007 [START_REF] Brenguier | Strain accommodation by slow slip and dyking in a youthful continental rift, East Africa[END_REF][START_REF] Froger | Timedependent displacements during and after the April 2007 eruption of Piton de la Fournaise, revealed by interferometric data[END_REF] and confirmed to be ongoing [START_REF] Peltier | Long-term mass transfer at Piton de la Fournaise volcano evidenced by strain distribution derived from GNSS network[END_REF][START_REF] Chen | Long-term ground displacement observations using InSAR and GNSS at Piton de la Fournaise volcano between 2009 and 2014[END_REF], systematic studies of eruption sequences [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion hotspot) inferred from deformation and geochemical data[END_REF] show no clear evidence of an alternation of dikes and sills at a given location, though these studies covered only limited periods of time.

The rotation could also be favored by the unloading induced by successive caldera collapses [START_REF] Corbi | How caldera collapse shapes the shallow emplacement and transfer of magma in active volcanoes[END_REF]. Among the 6 eruptions we identified as displaying similar displacement patterns (see Fig. S1), both modeled intrusions (June 2000 and May 2016) [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF] present similar shapes and it is likely that all these eruptions share the same source shape, of a sill turning into a dike on its western edge. Only one occurred before the last major collapse which affected the summit area in 2007 [START_REF] Staudacher | The April 2007 eruption and the Dolomieu crater collapse, two major events at Piton de la Fournaise (La Réunion Island, Indian Ocean)[END_REF] by deepening the crater by 330 m. The other five took place after this event. Thus, it is likely that the 2007 caldera collapse created a stress field more favorable to sill intrusions.

Another possibility is that dikes turn into sills, when they encounter a more rigid material, as shown by analog experiments [START_REF] Kavanagh | An experimental investigation of sill formation and propagation in layered elastic media[END_REF] and numerical experiments [START_REF] Maccaferri | A numerical model of dyke propagation in layered elastic media[END_REF], or because they encounter a preexisting discontinuity [START_REF] Delaney | Field relations between dikes and joints: emplacement processes and paleostress analysis[END_REF][START_REF] Maccaferri | A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation[END_REF]. Because the modeled structure has a slight northeast dip, its geometry could correspond to the upper limit of a flank sliding structure repeatedly activated during magma injection. Because the intrusion geometry is aligned with large scarps in the caldera rim, this structure perhaps corresponds to the limit of the Grand Brûlé depression (Fig. 1) [START_REF] Bachèlery | Evolution volcano-structurale du Piton de la Fournaise depuis 0[END_REF][START_REF] Merle | Hybrid collapse mechanism at Piton de la Fournaise volcano, Reunion Island, Indian Ocean[END_REF],

and may connects further east to the flank detachment associated with the 2007 eruption [START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF].

Step-wise lateral propagation

The mean velocity estimated for the lateral propagation is 0.6 m.s -1 , which lies in the range of previous estimates of 0.2-0.8 m.s -1 made for PdF, obtained from tilt and seismicity data [START_REF] Peltier | Imaging the dynamics of dyke propagation prior to the 2000-2003 flank eruptions at Piton de La Fournaise, Reunion Island[END_REF]. However, both the location of the propagating fracture front (Fig. 10A) and the temporal evolution of the area and pressure, (Fig. 10B) reveal sudden variations in lateral propagation velocity, corresponding to an instantaneous velocity of 2 m.s -1 between 21:10 and 21:15, then a pause from 22:00 up to the final vertical propagation that started at 3:15.

One of the questions is whether the 21:15 velocity jump is real whether it is an artifact induced by the limited number of stations to the southeast of the volcano. In order to test this hypothesis, we performed synthetic simulations where we simulated the cGNSS displacements induced by a source spreading and/or migrating at a constant rate (Details are provided in the Supporting Information Text S8, Figs S11 andS12). We determined that displacement changes similar to those observed could be obtained provided the source was both extending and migrating at the same time. Thus, considering the uncertainties on cGNSS measurements, we cannot completely exclude that a continuous evolution of the area and pressure induced the observed cGNSS data, and that the interpretation of a velocity jump in the magma propagation is an artifact.

However, the number of earthquakes provides a strong constraint, and can be used as a proxy to the propagation history. Laboratory experiments of hydraulic fracturing show that propagation of a hydraulic fracture proceeds in several stages [START_REF] Zoback | Laboratory hydraulic fracturing experiments in intact and pre-fractured rock[END_REF].

The fracture first initiates when the tensile strength is overcome (the initiation stage), and it takes some time and a higher pressure for a macroscopic fracture to form and for fluids to start propagating in the fracture (the breakdown stage). At this stage, pressure drops.

These different stages are associated with different rates of micro-fracturing as recorded by acoustic emissions. The initiation stage corresponds to an increase in micro-fracturing, and the breakdown stage corresponds to a decrease in seismicity rate. Using rate and state friction [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF], it is possible to relate the micro-fracturing to the stress history.

Thus, during injection, prior to the onset of fracture propagation, stress and seismicity rate increase. When the fracture starts propagating, pressure and host rock stress decrease, and seismicity returns to a background rate. Note that, after a stress decrease, seismicity does not immediately stop, as the time to return to a background rate is an inverse function of the stress rate [START_REF] Toda | Evidence from the AD 2000 Izu islands earthquake swarm that stressing rate governs seismicity[END_REF]. During the May 2016 PdF magma propagation, seismicity showed two peaks during the lateral magma propagation phase (Fig. 2B) at 21:15 and at 21:45. The first peak correlates with a rapid front advance (Fig. 10A), corresponding to an area increase and a pressure decrease (Fig. 10B). The second peak does not correspond to any area or pressure change, but this lack of evidence might be related to the cGNSS network configuration, and the lack of to stations southeast of the volcano. Because of the simultaneous timing of seismicity rate and propagation steps inferred from cGNSS, we can conclude that the 21:15 propagation step is not an artifact induced by the cGNSS network configuration. It is also likely that there was another propagation step at 21:45.

The residuals observed in the InSAR measurement (Figs. 7B and11B) suggest an eastward displacement of the east flank that the model is not able to explain. Strain acceleration evidenced by surface displacements and seismicity rates during the intrusion initiation indicates that damage and strain softening might occur in the edifice [START_REF] Carrier | A damage model for volcanic edifices: Implications for edifice strength, magma pressure, and eruptive processes[END_REF]. Similarly, as proposed by [START_REF] Got | Edifice strength and magma transfer modulation at Piton de la Fournaise volcano[END_REF] for 2007, and frequently suspected during eruptive crises at PdF, stress accumulation may cause the eastern flank to exceed a plastic threshold, leading to eastward slip which may further facilitate lateral migration of an intruded sill. It is likely that the stress decrease due to eastward flank slip led to a burst in magma inflow (Fig. 10B) which could explain the jump in lateral propagation observed at 21:15.

Limited magma supply and rapid disconnection from the reservoir

An unusual feature of the magma propagation during the May 2016 eruption is the 5 hours pause (phase III) observed from 22:00 (Fig. 10), leading to an exceptionally long duration of the intrusion phase of 8h25min. Generally, to reach that distance from the summit, it should have only lasted 3h30 (Fig. 2). The small volume lava flow (0.5 Mm 3 ) could have been fed with no or little input of new magma from the reservoir after the sill injection. A rough estimation of the volume of magma involved in the whole eruption (3 Mm 3 ) corresponds to the sum of the magma volume estimated from the inversion of interferograms covering the whole eruption (2.5 Mm 3 , Table 2) and the erupted volume (0.5 Mm 3 ). Since we ignore the effect of magma compressibility, while summing estimated volume of a source at depth and volume of lava flows, it is likely that this total volume is slightly overestimated. However, this volume is only slightly (0.4 Mm 3 ) larger than the volume emplaced in the sill (2.6 Mm 3 , Table 2) prior to the final dike propagation and eruption. Moreover, the overpressure for the whole eruption (1.5 MPa, Table 2) is 0.7 MPa smaller than the overpressure (2.2 MPa, Table 2) estimated for the sill part of the eruption, indicating that the eruption opened the system and released some pressure. Lastly, the lava flow volume emitted only represents 20% of the magma supplied to the volcano during this eruption, whereas it usually represents 85% [START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion hotspot) inferred from deformation and geochemical data[END_REF][START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF], suggesting that the volume supplied for this eruption was small and mostly trapped at depth.

Because the magma inflow rate was low, and the pressure increase from 22:00 was small (Fig. 10), the intrusion might have lacked the overpressure required to overcome the barriers to propagation and thus been more prone to stop than an intrusion with a larger overpressure [START_REF] Maccaferri | A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation[END_REF]. The low magma supply rate favors a scenario with a single small batch of magma disconnected from its feeding reservoir. This bottom closure could be due to a small reservoir rapidly depressurized by the intrusion propagation.

This intrusion could have stopped before reaching the surface either due to the surrounding stress field [START_REF] Maccaferri | A numerical model of dyke propagation in layered elastic media[END_REF] or to cooling [START_REF] Taisne | Effect of solidification on a propagating dike[END_REF], leading to a failed eruption. However, propagation restarted vertically, without any evidence of a new burst in influx of magma.

Final dike propagation : gas exsolution or East flank slip ?

The final dike is radial and rooted at the western edge of the deformation pattern where the thickness of rocks above the sill is at its minimum (880 m below the surface).

As a consequence, the dike initiated where the lithostatic stress was minimum, and where it was easiest to intrude. Two factors might have triggered the final intrusion either independently or in conjunction. One internal factor is related to the role of gas, while the other external factor is related to the role of the east flank slip (Fig. 12).

As intrusions propagate and reach higher elevations, gas exsolves, migrates and accumulates at the intrusion tip. In basaltic magmas, accumulated buoyant bubbles can even separate from the surrounding melt [START_REF] Menand | A phenomenological model for precursor volcanic eruptions[END_REF][START_REF] Maimon | The propagation of a dyke driven by gas-saturated magma[END_REF]. [START_REF] Menand | Gas segregation in dykes and sills[END_REF] demonstrated that, because of cooling and solidification, gas segregates more efficiently in sills than dikes, and that gas segregation is more efficient for low supply rates. Here, the conditions are met for efficient gas segregation : the intrusion is a sill and the supply rate is low. Moreover, the dike part of the intrusion is located at the highest elevation of the sill, where gas would tend to accumulate. Although, there was little visual observation for this eruption, lava fountaining activity was reported, characteristic of a high gas content [START_REF] Coppola | Lava discharge rate and effusive pattern at Piton de la Fournaise from MODIS data[END_REF]. Laboratory observations show that the part of fractures filled with gas have less lateral extension than the liquid filled part of fractures [START_REF] Menand | A phenomenological model for precursor volcanic eruptions[END_REF]. Here, the intrusion geometry corresponds to a large sill, connected to a narrow dike. This narrowing pathway is consistent with upward fracture propagation triggered by gas segregation. Another argument in favor of the role of gas, is that 49 Long Period (LP) events were recorded during phases II and III of the seismic crisis. According to [START_REF] Zecevic | Investigating the source characteristics of long-period (LP) seismic events recorded on Piton de la Fournaise volcano, La Réunion[END_REF], LP events are not common at PdF. When they do occur, they are generally associated with lateral propagation which leads to proximal or distal eruptions [START_REF] Aki | Seismic monitoring and modeling of an active volcano for prediction[END_REF]. The source of LP events at PdF is not well understood, but the resonance of gas filled fractures is an efficient way of generating them [START_REF] Chouet | Long-period volcano seismicity: its source and use in eruption forecasting[END_REF].

Sill intrusion is expected to favor flank slip [START_REF] Famin | Volcano destabilization by magma injections in a detachment[END_REF][START_REF] Got | Edifice strength and magma transfer modulation at Piton de la Fournaise volcano[END_REF]Chaput et al., 2014b]. It relaxes stress perpendicular to the rift and may allow the dike to intrude along the rift. Since 2007 the east flank is known to slip at a steady rate of 0.02 m.yr -1 [START_REF] Chen | Long-term ground displacement observations using InSAR and GNSS at Piton de la Fournaise volcano between 2009 and 2014[END_REF], and transient slip sometimes occurs as a response to intrusions [START_REF] Brenguier | Strain accommodation by slow slip and dyking in a youthful continental rift, East Africa[END_REF][START_REF] Peltier | Long-term mass transfer at Piton de la Fournaise volcano evidenced by strain distribution derived from GNSS network[END_REF][START_REF] Froger | Timedependent displacements during and after the April 2007 eruption of Piton de la Fournaise, revealed by interferometric data[END_REF]. Three independent data suggest that a limited transient east flank slip carried on after the initial slip at 21:15: records on FERG station (up to 1 cm of eastward displacement during the week after the eruption), second intermediate InSAR S1 D2, and residuals from the inversion of the four interferograms spanning the whole eruption. We assume that such a transient slip would favor gas exsolution leading in turn to reactivation of the propagation vertically.

Conclusion

In this study, we combined InSAR and continuous GNSS (cGNSS) data to gain insights into the dynamics of a magmatic intrusion at Piton de la Fournaise. An InSAR image was acquired as the May 2016 eruptive crisis was underway, allowing us to test a specific approach dedicated to integrating InSAR and cGNSS data. We use the deformation source determined from the inversion of InSAR data covering the whole eruption as an a priori for tracking a circular pressurized area from cGNSS data. This method estimates the source depth intrinsically, while this depth is wrongly estimated when only the cGNSS data is used without any a priori derived from InSAR.

The best modeled intrusion corresponds to a 2700 m long sill, dipping to the northeast and turning into a 880 m high dike located on the western limit of the sill, along the south east rift zone. This intrusion has a similar sill to dike shape as the one determined for the June 2000 event. Because InSAR patterns are similar for the June 2000, October a complex 3D geometry and potentially controlled by multiple factors. Assessing the relative contributions of the proposed processes i.e. stress inherited from the 2007 caldera collapse, the east flank instability and the gas exsolution, on the control of the magma propagation dynamics will be the subject of further work. parameters for a disk projected on an initial mesh (Projected Disk method). (See Table 1 for the description of all parameters.) S3).

Figure 1 .

 1 Figure 1. Piton de la Fournaise (PdF) setting. (A) Location of La Réunion Island. (B) Relief of the PdF volcano. Black lines refer to the axe of the 3 rift zones. (C) Locations of the lava flows emplaced from June 2014 to May 2016 and a portion of the OVPF monitoring network. Black and grey diamonds represent the locations of the permanent GNSS stations. The black circled diamonds refer to the GNSS stations used in this study, the grey diamonds are stations not used here. Red dots and blue triangles represent the location of the OVPF tiltmeters and seismometers, respectively. 'Dol' refers to the Dolomieu crater. Coordinates in kilometers (WGS84, UTM 40S).

Figure 2 .

 2 Figure 2. (A) Distance of the highest point of eruptive fissures from the center of the Dolomieu crater as a function of the seismic crisis durations. Blue squares and green triangles represent summit eruptions and flank eruptions between 1985 and 2010 (from Roult et al. [2012]), respectively. Orange diamonds represent eruptions from 2014 and 2015. The purple triangle represents the May 2016 eruption. The green line represents the linear trend. (B) Top Ratio of RSAM amplitude between DSM and FOR stations (see Fig. 1 for location). Middle NS (black dashes) and EW (blue line) tilt and temperature (orange line) variations at DSOi station (see Fig. 1 for location). Bottom Histogram of the number of volcano-tectonic (VT) events per 5 min periods during the May 2016 seismic crisis (left axis). The red curve represents the cumulated number of VT events(right axis). The blue curve represents the cumulated number of long period events (LP, left axis).Colored areas with numbers I, II, III and IV refer to the different stages of the intrusion propagation described in 3.2. The black arrow represents the acquisition time of the Sentinel image used to compute interferograms S1 D1 and S1 D2 (Fig.3). Times are UTC.

Figure 3 . 1 Figure 4 .

 314 Figure 3. (Top) Time spanned by each interferogram used in this study. The vertical red area represents the eruption. (Bottom) Zoom, colored areas corresponding to numbers I, II, III and IV represent the different phases of the intrusion propagation as revealed by the seismic and tilt data (see Fig. 2 and paragraph 3.2 for a description). Phase II also corresponds to the period when significant deformation was recorded on the cGNSS stations. The vertical dashed line at 01:45 marks the time of the intermediate SAR acquisition by Sentinel-1

Figure 5 .

 5 Figure 5. GNSS data for the May 2016 eruption. (A) Map view of the horizontal displacements from cGNSS data and campaign GNSS data. Colored circles represent the evolution of the horizontal displacement with time for the 10 cGNSS stations. Black arrows represent the cumulated displacement of the 10 cGNSS stations between May 25 at 17:00 and May 26 at 04:05. Standard deviations for horizontal and vertical components are 1-2 cm and 3-4 cm, respectively. Grey arrows represent displacements from the campaign GNSS network between August 31, 2015 and May 27, 2016. Ellipses represent their 95% confidence intervals. The lava flow is shown in purple. Coordinates are in kilometers (WGS84, UTM 40S). (B) Cumulated displacement (in meters)is plotted as a function of time for the 10 cGNSS stations which recorded the intrusion. Crosses represent the sampling dates used in the time step inversion. Black arrows mark the eruption onset.Colored areas with numbers I, II, III and IV refer to the different stages of the intrusion propagation described in section 3.2. Times are UTC.

Figure 6 .

 6 Figure 6. Representation of the parameters defining the fractures used as deformation sources. (A) 9 parameters for a curved quadrangular-shaped intrusion linked to the eruptive fissure by a single echelon. (B) 3

Figure 7 .

 7 Figure 7. A. Geometry and opening of the source (best-fit model) determined by inversion of the 2 descending and 2 ascending interferograms spanning the whole eruption. Black and red stars represent hypocenters of located micro-earthquakes occurring on May 25 and 26, respectively. B. One-dimensional (diagonals)and two-dimensional (off-diagonals) marginal posterior probability density (PPD) functions plotted with the maximum PPD (blue) and mean (red) models. Note that the PPD are zoomed in on the best fitting parameter range which is much smaller than the explored intervals (see TableS3).

Figure 8 .

 8 Figure 8. Unwrapped displacements along the LOS in meters, from the interferograms (top), model (middle) and residuals (bottom). The black segment indicates the location of the eruptive fissure. Black arrows represents the LOS direction. Coordinates are UTM.

Figure 9 .

 9 Figure 9. Projected Disk method results. A.Comparison between the intrusion determined from InSAR data covering the whole eruption (grey) and the intrusion determined for the first part of the eruption using the Projected Disk method with cGNSS data, InSAR S1 D1 data or both datasets (blue) B. Best (blue squares) and mean (orange diamonds) model parameters for the Projected Disk method using either cGNSS, InSAR S1 D1 data and both datasets. Error bars represent the 95% confidence intervals.

Figure 10 .

 10 Figure 10. A. Location of the overpressure source determined from the temporal inversion (colors) for 14 time steps between 20:40 and 4:00 out of the 30 time step inverted. The pressurized area is assumed to be a disk projected on the a priori mesh (in grey) determined from the inversion of the four InSAR data covering the whole eruption. B. Results of the temporal inversion. Time evolution of the overpressure, the area and the volume of the source inverted for each time step. The best-fit model (blue line), the mean model (orange dashs), and the lower and upper limits of the 95% confidence interval (black dots), respectively, are shown. Circles shown at 01:45 represent results from joint inversion of InSAR and GNSS data. The flux is obtained by deriving the best-fit volume curve. Colored areas with numbers I, II, III and IV refer to the different stages of the intrusion propagation described in 3.2. Times are UTC time.

Figure 11 .

 11 Figure 11. A. Deformation sources determined from the inversion of Sentinel-1 interferograms covering the beginning (S1 D1, blue) and the end (S1 D2, red) of the propagation using the Projected Disk method. Inset represents a zoom on of the intrusion. The grey mesh is the a priori curved quadrangular source determined from the inversion of the four InSAR data covering the whole eruption. B. Comparison of data, modeled displacements and residuals for the best-fitting inverted sources of both parts of the eruption (before and after 01:45). The black segment indicates the location of the eruptive fissure. The first and last columns concern the inversions at the beginning (blue) and end (red) of the propagation. The four columns in the middle concern the whole intrusion. The first row is data, the second row is the model displacements determined for the beginning of the eruption (blue source in A). The third row compares the residuals between the interferograms covering the whole eruption (second to fourth columns) to the interferogram corresponding to the final part of the eruption (last column). The fourth row presents the model displacements obtained for the end of the propagation (red source in A). The bottom row is the remaining residuals. Arrows represent the LOS.

Figure 12 .

 12 Figure 12. Summary diagram illustrating the scenario proposed for the May 2016 intrusion propagation. Phase I : Vertical ascent first evidenced by seismicity increase (a) then in addition by tilt (b). Phase II : Lateral propagation (left) potentially favored by flank slip right). Phase III : Pause with eventual gas accumulation (left)and/or flank creeping (right). Phase IV : Final dike propagation marked by seismicity increase (a) then decrease (b) before the eruption onset.

Acknowledgments

The PhD fellowship of Delphine Smittarello is funded by the French Ministry for Higher Education and Research. This research was supported by the INSU-CNRS ALEAS program. Part of this work was also supported by the ANR (Agence Nationale de la Recherche) through the SLIDEVOLC project (contract ANR-16-CE04-004-01). This research was partly supported by the French Government Laboratory of Excellence initiative ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Fund. The seismic and GNSS data used in this paper were collected by Observatoire Volcanologique du Piton de la Fournaise / Institut de Physique du Globe de Paris (OVPF/IPGP). We thank the team that manage the OVPF, install, maintain and monitor stations. Raw seismic and GNSS data are available on the website: http://volobsis.ipgp.fr. Processed data are available by contacting the OVPF. We thank the European Space Agency for providing the Sentinel-1 data and the Italian Space Agency for providing the Cosmo SkyMed data. All interferograms used in this study were computed with Diapason ((c) TRE ALTAMIRA) by the Indian Ocean InSAR Observatory Service OI 2 , part of the National Service for Volcanological Observations (SNOV) of the French National Research Council Institute for Earth Sciences and Astronomy (CNRS/INSU). The interferograms are available on the web interface of the OI 2 platform, by following the link https://wwwobs.univ-bpclermont.fr/casoar. Seismic catalog, GNSS time series, all interferograms used in this study and the meshes considered in the boundary element modeling (topography and final intrusion shape) can be downloaded via the link http://doi.org/10.5281/zenodo.2548686 The matlab inversion codes and a users manual can be provided on demand by Yo Fukushima (fukushima@irides.tohoku.ac.jp) This paper has benefited from interesting and helpful discussions with A. Walsperdorf, R. Grandin, M. Métois for GNSS post-processing; B. Valette about inversion and bayesian inference, P. Lesage and F. Brenguier about seismic data, Alexis Hrysiewicz about In-SAR. We also thank T. Shreve and F. van Wyk de Vries for english reviewing of this manuscript. Mike Poland, Scott Henderson and an anonymous reviewer who helped to improve this manuscript. This is Laboratory of Excellence ClerVolc contribution number 328 and IPGP contribution number 4003.

Confidential manuscript submitted to JGR-Solid Earth 2010, May 2015, May 2016, January and July 2017 eruptions, these eruptions probably correspond to the same intrusion geometry. Because most of these intrusions occurred after the 2007 caldera collapse, it is possible that they were guided by the stress field resulting from the collapse. Alternatively, their geometry, as well as their location, suggest that they might follow a pre-existing discontinuity, which might correspond to the upper limit of the Grand-Brûlé depression.

The May 2016 intrusion is characterized by two changes in direction from vertical to lateral, then from lateral to vertical. From the seismicity, it is likely that the reservoir started fracturing at 19:40 (Fig. 12, phase Ia) and from the tilt, that the dike started propagating upward at 21:10 (phase Ib). Then, at 20:40, the intrusion turned laterally to the southeast as a sill (phase II). Inverting cGNSS data at 30 successive time steps, we found that lateral magma propagation was unsteady. Magma slowly propagated until 21:15, then the sill area suddenly increased to reach its maximum dimension at 22:00. The velocity jump is confirmed by a peak in seismicity and could be associated with east flank sliding.

Pressure and volume kept on increasing until midnight, while the intrusion was paused.

After 5 hours of stand-by, the increase in rate of micro-earthquakes marked fracture initiation (phase IVa). This was followed by a decrease in the rate of micro-earthquakes which marked the onset of upward propagation at 03:15 (phase IVb), leading to a fissural eruption at 04:05. From volume budgets and the estimated magma flux, it is likely that the eruption was fed by little or no input of new magma after the initial supply of a single batch of magma. This could explain why the intrusion lasted longer than expected considering the distance between the eruptive fissure and the summit. This intrusion could have been non-eruptive if it had not been for the additional factors favoring the ultimate The proposed method relies on the integration of InSAR and cGNSS data. Since it requires InSAR data covering an intrusion, so far, it can not be used in a predictive way.

However, it is a powerful way to track magma propagation fronts, and pressure changes after an eruption. It provides a valuable picture of an unsteady process characterized by