Earth and Planetary Science Letters The 2007 caldera collapse of Piton de la Fournaise volcano: Source process from very-long-period seismic signals

In April 2007, Piton de la Fournaise volcano experienced its largest caldera collapse in at least 300 yr. This event consisted of a series of 48 subsidence increments characterized by very-long-period (VLP) seismic signals equivalent to M w 4.4 to 5.4. Source analysis of VLP events suggests a piston-like mechanism with a collapsing crack source corresponding to the contraction of the magma reservoir and a single force representing the collapse of the above rock column. We show that the collapse dynamics is primarily controlled by magma withdrawal from the reservoir, which was likely in a bubbly state at the time of the event. Similar to the 2018 Kilauea collapse, we observe a reduction of the time-interval between successive subsidence increments, which results from the acceleration of magma withdrawal and a progressive weakening of the ediﬁce at the beginning of the sequence. © 2019 The Author(s).


Introduction
Although caldera collapses are frequent in the evolution of basaltic systems, they have been monitored at only a few volcanoes worldwide (e.g., Kumagai et al., 2001;Michon et al., 2007;Neal et al., 2019). Such collapses usually consist of a quasiperiodic series of discrete subsidence events manifested by tilt steps and very long period (VLP) signals (e.g., Fontaine et al., 2014;Kumagai et al., 2001). While the initiation mechanisms of such phenomena remains largely enigmatic, they generally involve a lateral magma withdrawal from a subsurface reservoir through low elevation vents. Observations indicate a complex feedback process between summit collapses and the dynamics of the eruptive vent. While magma outflow from a sub-surface reservoir is usually associated with summit deflations (e.g., Epp et al., 1983), collapses can also act as a mechanism to push the magma out of the chamber. Effusive surges and short-term edifice inflations following summit collapses suggest a piston-like mechanism, in which the depressurized magma chamber is regularly repressurized by the subsiding rock column (Kumagai et al., 2001;Michon et al., 2011;Neal et al., 2019).
Caldera collapses involve different source processes like dipslip motion on ring faults and the contraction of magma reservoirs. These mechanisms are often modeled using moment tensor parameters, whose inversion on volcanoes usually yields vertical * Corresponding author.
E-mail address: zacharie.duputel@unistra.fr (Z. Duputel). compensated-linear-vector-dipoles (CLVD) representing ring faulting (e.g., Shuler et al., 2013) or tensile crack sources associated with the pressurization of magma reservoirs (Kumagai et al., 2001;Mildon et al., 2016). Such a moment tensor representation, however, ignores the effects of the detached mass that can be significant during gravity-driven caldera collapses. As far as the excitation of long period seismic waves is concerned, the descending collapse rock column induces unloading and loading of the volcanic edifice that can be modeled with a single force (e.g., Takei and Kumazawa, 1994). Because volcanoes involve mass advection processes generating forces on the Earth, the modeling of volcanic sources commonly requires the consideration of single forces in addition to moment tensor components (Nakano et al., 2003;Chouet and Matoza, 2013). A remarkable example is the 1980 Mount St. Helens eruption (USA), which combined a massive landslide represented by a horizontal single force (Kanamori and Given, 1982) and a series of volcanic explosions corresponding to the superposition of an isotropic source and a vertical force (Kanamori et al., 1984). Other examples of volcanic phenomena inducing single forces are dome collapses (e.g., Uhira et al., 1994) or magma movements in the edifice (e.g., Chouet et al., 2010).
In this study, we focus on the April 2007 collapse that affected the summit area of Piton de la Fournaise volcano located on La Réunion island (cf., Fig. 1). Although numerous collapses have occurred in the summit area since the 18th century, the 2007 event is the largest collapse reported at Piton de la Fournaise (Michon et al., 2013). It followed a particularly active period with 4 eruptions in 10 months ) that ended in a lateral low-altitude eruption with exceptionally high flow rate and a large volume of ejected lava . Early geodetic co-eruptive measurements show a significant deflation of the volcano before the collapse Fontaine et al., 2014;Froger et al., 2015), suggesting that the low-altitude eruption induced a sudden depressurization of the magma reservoir that ultimately resulted in the collapse of the rock column below the summit (Michon et al., 2007). This major failure within the edifice was accompanied by a sudden decrease of seismic velocities in the summit area (Duputel et al., 2009) and an exponentially decaying post-eruptive deflation of the central cone (Froger et al., 2015).
The surface collapse activity started on April 5 at 20:48 UTC with an earthquake of magnitude M ∼ 5 that could be observed at teleseismic distances ( Fig. 2 and Supplementary Fig. S1). This event was followed by more than 40 smaller collapses that occurred quasi-periodically until mid April 2007 (Michon et al., 2011;Fontaine et al., 2014). Each collapse was characterized by inflationary tilt steps and VLP signals (cf., supplementary Fig. S2; Fontaine et al., 2014). It was also accompanied by an intense microseismic activity below the volcano summit . The collapse resulted into a depression of ∼330 m of the Dolomieu summit crater with a total volume similar to the bulk volume of emitted lava (∼10 8 m 3 ; Michon et al., 2007;Urai et al., 2007;Staudacher et al., 2009). Although numerous studies have focused on geodetic deformation, seismicity and velocity variations associated with this activity Froger et al., 2015;Brenguier et al., 2008;Massin et al., 2011;Got et al., 2013), only a limited number of studies focused on the actual source mechanism of the VLP signals generated by collapse events (e.g., Fontaine et al., 2019). In this study, we thus combine near-field observations at the Geoscope station RER (8.5 km from the volcano summit) and surface waves observed at 25 teleseismic stations to investigate the source of collapse events that affected Piton de la Fournaise in April 2017. To mitigate any bias due to the effect of lateral heterogeneities on teleseismic surface waves, our analysis accounts for 3D heterogeneities using the spectral element approach. We then explore different source representations to constrain the different source processes involved during the April 2007 caldera collapse.

Seismological observations and waveform simulations
We select 25 teleseismic broadband stations for which the VLP signal generated by the main collapse event is clearly visible. The resulting dataset has an azimuthal gap of 30 • and is composed of stations within 90 • of epicentral distances. Seismic waveforms are deconvolved to velocity and band-pass filtered between 50 s and 100 s. Periods shorter than 50 s are filtered out to mitigate the effect of short-scale structural heterogeneities that are not resolved in global 3D velocity models. The 100 s corner period reduces the effect of long-period noise, in particular for stations in Indian ocean southern islands and Antarctica (removing those stations would result in an azimuthal gap exceeding 130 • ). In the 50-100 s period range, the observed seismograms are dominated by fundamental mode surface waves for which a significant fraction of the energy stays near the free surface where lateral structural heterogeneities are significant at global scale (e.g., due to oceans and continents). To account for such strong lateral variations, we conduct 3-D spectral-element simulations to compute the Green's functions used for source characterization (Komatitsch and Tromp, 1999). We assume a global 3-D Earth model composed of S40RTS (Ritsema et al., 2011) and Crust2.0 (Bassin et al., 2000). In addition to 3-D variations of wave speeds and density, our simulations include the effect of ellipticity, topography, bathymetry, oceans and self-gravitation.
We complement teleseismic records with near-field observations at the Geoscope station RER, located 8.5 km away from the Piton de la Fournaise summit area. Given its proximity to the source, the RER station provides direct information on the overall duration and centroid time of the source. It is modeled using Herrmann (2013) implementation of the wavenumber integration method for the local velocity model used at the Piton de la Fournaise Volcano Observatory (OVPF; model derived from Battaglia et al., 2005). We only use the vertical component seismogram and do not incorporate horizontal channels that are significantly contaminated by tilt signals (Fontaine et al., 2014). Before inversion, the seismogram is deconvolved to velocity and band-pass filtered in the period range of 5-100 s. The resulting waveform consists of a long-period signal with a relatively simple sinusoid pulse with a duration of ∼17 s (see Supplementary Fig. S2a).

Source mechanism of the main collapse event on April 5, 2007
As mentioned above, the April 2007 collapse can possibly involve different source processes that can be modeled at long period by a moment tensor source (e.g., slip on a ring fault, contraction of a magma reservoir) or by a single force (e.g., unloading and loading of the volcanic edifice by the detaching rock column). In the following, we confront these different source representations to the seismological observations of the main VLP event on April 5 at 20:48 UTC.

Moment tensor inversion
We employ a strategy similar to Duputel et al. (2016) where centroid moment tensor (CMT) parameters are inverted for using surface waves in a 3-D Earth model. The parameters to be determined are the 6 elements of the seismic moment tensor as well as the centroid location and source duration. Although we have tested different centroid locations, our tests show a very limited sensitivity of long period seismic waves to the source location at the and t F h are respectively the moment-tensor and force half-durations. M w is the moment magnitude, M 0 is the scalar seismic moment, F M A X is the maximum force and M C S F is the CSF value (cf., Kawakatsu, 1989). Comparison between data (black) and synthetic (red) traces are shown for representative stations. The part of the waveform used for source inversion is delimited by red dots. The blue star in right insets indicate the source location. Yellow circles show the distribution of stations used for the inversion while the red circle indicate the location of the station whose waveform is presented. The station azimuth (φ) and epicentral distance ( ) is indicated on top of each seismic trace. 100 m scale of the Piton de la Fournaise summit area. We therefore fix the centroid below the Dolomieu summit crater (21.244 • S, 55.714 • E) and only grid-search for the best point-source depth (from 0.5 km to 3.5 km). The centroid time and duration are also estimated using a simple grid-search strategy assuming an isosceles triangular moment rate function (MRF). Our preferred CMT solution for the main April 2007 collapse is presented in Fig. 2a and Fig. 3a. The optimum centroid time is 20:48:49 UTC and the source half-duration is 3 s (cf., Supplementary Fig. S3). We find a best point-source depth 2 km below the free surface (i.e. ∼500 m above sea level). This estimate is however subject to large uncertainties as waveform fits are similar if we change the source depth.
Although estimated posterior uncertainties are relatively small, results suggest a large negative correlation between isotropic and vertical compensated linear vector dipole (CLVD) components (cf., Supplementary Table S1). To explore the contribution of the isotropic component to fit the data, we conduct an inversion assuming the trace of the moment tensor to be zero (i.e., the deviatoric moment tensor). Results in the Supplementary Fig. S4 show that such source cannot fit data properly (in particular at the near-field station RER), therefore confirming that an isotropic component is necessary to explain available observations. To further assess what mechanism is requested by our observations, we estimate the minimum RMS misfit associated with every type of moment tensor source (cf., Supplementary Text S1). The resulting source-type diagram in Fig. 3b allows us to estimate what type of moment tensor source is consistent with the observed waveforms (Ford et al., 2010). Results suggest that the main collapse is consistent with a vertical closing crack, in agreement with the values and orientation of the principal axes in Fig. 3a. The source-type diagram in Fig. 3b also reveals uncertainty on the double-couple component of the source, which results from poor constraints on M rθ and M rφ , due to the small amplitudes of the associated excitation kernels (see Table S1; Kanamori and Given, 1981).
Our moment tensor model is qualitatively similar to the moment tensor solution of Fontaine et al. (2019). However, the small scalar moment (M 0 = 1.5 × 10 16 N.m), negative centroid time-delay (−4 s), very long duration (20 s) and large depth (5 km) reported by Fontaine et al. (2019) is inconsistent with both near-and farfield observations. Fig. S17 indicates that waveform predictions from such source model are significantly smaller than observations by about 1 order of magnitude, which likely results from a bug in their inversion procedure (Hejrani and Fontaine, personal communication). Beyond this scaling issue, additional differences could potentially be attributed to the use of a spherical model (PREM) and the fact that RER was not incorporated in their inversion.

Single force inversion
To investigate possible effects due to the unloading and loading of the collapsing rock column, we also perform a centroid single force (CSF) inversion using an approach similar to Kawakatsu (1989). At long period, the effect of a detached mass can be modeled by a single time-varying force (e.g., Kanamori and Given, 1982;Ekström and Stark, 2013). To represent the effect of the descending mass acceleration and deceleration, we assume a sinusoidal source time history such that the final integral of the force equal to zero (e.g., Kawakatsu, 1989). We invert for the 3 components of the peak force vector along with the source timing and duration.
Results in Fig. 4 indicates a predominantly vertical force vector that is consistent with the downward motion of the rock column (i.e., the initial upward orientation of the force is in the opposite direction of the accelerating rock mass). In agreement with the moment tensor inversion, we obtain a centroid time at 20:48:49 UTC and a half-duration of 4 s (cf., Supplementary Fig. S3b). As far as the overall fit is concerned, the CSF solution in Fig. 2b seems to fit the data as well as the moment tensor model shown in Fig. 2a. The posterior correlation matrix presented in Supplementary Table S2 do not show significant correlations between the different components of the inverted peak force vector. If we constrain the force vector to be purely vertical, we do not notice a significant deterioration of the waveform fit (cf., Supplementary Fig. S5), which indicates that the small south-west component of the force in Fig. 4a is not required to fit the data. This is consistent with the fact that horizontal force components are associated with larger uncertainties (cf., Table S2). Force time-history. The blue dot labeled F MAX indicate the time when the force is maximum (i.e., the initial peak force). Waveform fit is presented in Fig. 2b and S18.

Joint inversion of moment tensor and force
We also conduct a joint inversion of moment tensor and force parameters. Given our limited sensitivity to the centroid depth and the force horizontal components, we assume a vertical force and fix the source 2 km below the free surface. We thus invert for the 6 elements of the moment tensor, the maximum force amplitude along with their centroid times and durations. The solution pre- and t F h denote the half-duration of the moment tensor and force source, respectively. The centroid time-shift of the moment tensor (τ MT ) and force (τ F ) are defined with respect to 20:48:43UTC. λ, δ and φ indicate respectively the eigenvalue or force, plunge and azimuth angles of the principal axes. The lower right inset in (b) indicate the force time history. The blue dot labeled F MAX indicate the time when the force is maximum (i.e., the initial peak force). Waveform fit is presented in Fig. 2c and S19. sented in Fig. 2c and Fig. 5 is obtained for a half-duration of 3 s for both single force and moment tensor sources (consistently with CMT and CSF solutions presented above).
The inverted focal mechanism in Fig. 5a is consistent with a vertically closing crack and is very similar to the one obtained in Fig. 3. The force source in Fig. 5b is consistent with a downward moving rock column as the single force model in Fig. 4. Optimum centroid times indicate that the closing crack source occurs ∼2 s before the upward force is applied. However, this delay remains small compared to the sampling rate of the seismological data (1 sample per second). The retrieved peak force amplitude is correlated with the moment tensor parameters (e.g., with the vertical CLVD component; cf., Supplementary Table S3). As in Table S1, we also notice a large negative correlation between isotropic and vertical CLVD components. Despite these tradeoffs, Fig. 2 shows that this joint model yields to a better waveform fit at the near-field station RER compared to moment tensor and single force sources presented before.

VLP precursors and smaller collapse events
Although smaller collapses that followed the April 5 event are not visible at teleseismic distances, the associated VLP events are well recorded at the RER Geoscope station (e.g., Fontaine et al., 2014;Fontaine et al., 2019). To obtain a catalog of these smaller VLP events, we use the ObsPy implementation of the STA/LTA algorithm on the vertical component of RER between April 5 and April 14, 2007. We use a short time average (STA) window of 20 s and a long time average of 400 s with thresholds of 3.4 and 3.0 for on and off trigger times. This detection procedure results into a catalog of 48 collapse events (including the main collapse) that are listed in Supplementary Table S4. As shown in Fig. 7, the main collapse event was also preceded by several VLP events. This is consistent with Fontaine et al. (2019) that reported 4 VLP events accompanied by tilt steps within 20 h before the main collapse.
Although the small amplitude of these events complicates any automated STA/LTA detection procedure, we can visually identify at least 7 VLP precursors on Fig. 7 that seem coincident with tilt changes reported by Fontaine et al. (2014) (see Tables S9-S11).
As for the main collapse event, the smaller collapses manifest as relatively simple long-period waveforms. Supplementary Fig.  S2a shows that waveforms generated by the successive collapses have a similar shape, indicating that these events correspond to similar source mechanisms. We therefore fix the source geometry according to CMT and CSF models obtained for the main event and investigate smaller collapses by only inverting the source depth, timing, duration and size (i.e., seismic moment or force magnitude). Results are summarized in Supplementary Tables S4-S6. The magnitudes of smaller events ranges from 4.4 to 5.0 and their half-duration from 3 s to 7 s. Waveform fits presented in the Supplementary Figs. S6-S10 show that joint moment tensor and force models usually perform better than single moment tensor or force sources. This is not surprising given that the joint model incorporates more parameters than the other source models.
As shown in Fig. S2b, VLP precursors are also relatively similar to the main collapse event, even if they appear to have a longer duration. As described before, we investigate the source of these VLP precursors assuming a similar source geometry as the main event. Results in Tables S9-S11 and Fig. S11 suggest that the joint moment tensor and force source is more appropriate to fit the observations. In this case, the observed differences with the main collapse event is attributed to a larger time-difference between the CMT and CSF sources. However, this remains quite uncertain as the analysis is based on waveforms recorded at a single station.

Discussion and conclusion
In this study, we investigate the caldera collapse that affected the summit of Piton de la Fournaise in 2007. The surface collapse activity started on April 5 with a M w ∼ 5.4 VLP event visible at teleseismic distances. This main event was then followed by 47 smaller collapses with similar source mechanisms. Using a moment tensor parameterization, the source of the main collapse is best represented by a vertically closing crack (cf., Fig. 3). Our results show that ring-faulting is not the dominant source mechanism as a CLVD source is unable to fit the observed VLP signals ( Fig. S4; Shuler et al., 2013). Using a single force representation, the inversion yields an initially upwards vertical force that is consistent with the downward displacement of the rock column above the magma reservoir. This single force model fits the observations as well as a moment tensor source even if it relies on fewer model parameters (cf., Fig. 2a-b). Model selection using log model evidence (LME), Akaike (AIC) and Bayesian (BIC) information criteria thus favor the simpler force model over a moment tensor description (see Supplementary Text S2 and Table S7). The joint inversion of force and moment tensor parameters suggests the combination of an initially upward vertical force and a closing crack. Despite its larger number of parameters, LME, AIC and BIC criteria favor this joint model over previous source parameterization because it is more consistent with observations (especially at the near-field station RER, see Table S7). The combination of a moment tensor and a force also results in better waveform fits and lower AIC and BIC values for smaller VLP events that followed the main collapse (see Supplementary Fig. S6-S10 and Table S8). Although the joint model seems to be the most appropriate according to LME, AIC and BIC criteria, it is important to note that this solution is derived from noisy data and is affected by tradeoffs between source parameters. In the following the different source representations will be taken into account to evaluate the dependance of estimated quantities on the assumed parameterization. A simple interpretation of the results presented above is that the collapsing crack source corresponds to the vertical contraction of the magma reservoir while the vertical force represents the resulting collapse of the overlying fractured rock column. This interpretation seems consistent with our joint source models showing that the closing crack source usually occurs before the vertical force is applied to the edifice (cf., Fig. 4 and Table S6). The combination of a collapsing crack and a single force has been already inferred to explain long-period signals observed on volcanoes (e.g., Nakano et al., 2003) in particular during caldera collapses (e.g., Kikuchi et al., 2001). Moreover, previous seismological and geodetic studies have shown that closing crack sources can correspond to magma chamber deflations during a caldera collapse (e.g., Riel et al., 2015;Mildon et al., 2016). The use of a single force has also been suggested to represent a subsiding caldera block (Kumagai et al., 2001).
Assuming a simple piston cylinder to represent the subsiding caldera block as illustrated in Fig. 6a, we can use the different source models discussed above to estimate the downward displacement of the collapsing rock column (cf., Supplementary Text S3). We assume a piston diameter of 450 m consistently with Urai et al. (2007), Michon et al. (2011) and Fontaine et al. (2019). The piston extends vertically from the roof of the magma reservoir (at an elevation of ∼300 m; Peltier et al., 2008;Duputel et al., 2019) to ∼1700 m above sea level at the upper limit of CLVD and normal faulting events observed during the collapse sequence . For each model, we compute a lower-bound estimate of downward displacement for a density ρ = 2.7 g cm −3 (taken from the velocity model used at OVPF) and an upper-bound using ρ = 1.6 gcm −3 (according to the gravity study of Gailler et al., 2009). Moment tensor crack estimates use Lamé parameters λ = 6-11 GPa and μ = 8 − 13 GPa corresponding to the aforementioned values of ρ along with P and S wave speeds derived from local velocity models (Battaglia et al., 2005;Mordret et al., 2015). The results shown in Fig. 6b and Fig. S12 indicate a total piston displacement ranging from z ∼ 190 m to z ∼ 635 m depending on the source model and the value of ρ.
The 450 m estimate proposed by Michon et al. (2011) fall in the middle of this range. Interestingly, our joint source model produces consistent piston displacement estimates whether they are derived from moment tensor or force parameters. In this case, we estimate a total displacement z ∼ 330 m assuming ρ = 1.6 g cm −3 , in agreement with direct observations reporting a depression depth of ∼330 m in the Dolomieu crater (Michon et al., 2007;Urai et al., 2007). The single force model is associated with larger estimates (i.e., z ∼ 370-635 m) that are also consistent with field observations if we assume a denser piston cylinder. The moment tensor estimates z ∼ 220-370 m are also consistent with the observed final depth of the caldera. As shown in Fig. S12, the evolution of the piston displacement is globally consistent with Fontaine et al. (2019) that estimated a final depth of 340 m. For simplicity, we will now assume ρ = 1.6 g cm −3 for the moment tensor and joint model and ρ = 2.7 g cm −3 for the single force model.
The above estimates of vertical displacement can be used to evaluate the volume change of the magma reservoir after each collapse. Dividing this volumetric change by the time-delay T  in Tables S9-S11 of the online supplementary). Events F1-F2, F3, F4 and F5 correspond respectively to events E2, E3, E4 and E5 reported by Fontaine et al. (2019). Events F6-F7 are coincident with a tilt step that occurred a few minutes before the main collapse (e.g., Fontaine et al., 2014). since the previous collapse, we obtain a volumetric change rate α = z S/T (where z is the downward displacement increment and S is the cross-sectional area at the base of the piston). This volumetric change rate shown in Fig. 6b can be compared directly to the magma outflow rate. It has been showed that the high-frequency ground motion velocity amplitude near the eruption site has similar evolution to the emission rate during the April 2007 eruption Coppola et al., 2009;Michon et al., 2011). Using the approach of Hibert et al. (2015), we therefore estimate the magma outflow rate from continuous records at station TKR, located close to eruptive fissure (see Supplementary Text S4). The estimated lava extrusion rate shown in black in Fig. 6b is consistent with volumetric change rates estimated from our VLP catalog. This confirms that the 2007 Piton de la Fournaise collapse is primarily controlled by lateral magma withdrawal from the shallow reservoir. Our magma outflow estimates are lower than preliminary volumetric change rates obtained empirically by Michon et al. (2011) from seismic acceleration at the RER station reaching up to 1000 m 3 s −1 . On the other hand, our average emission rate of 61 m 3 s −1 is roughly consistent with field observations of Staudacher et al. (2009) that reported an overall mean lava outflow of 52 m 3 s −1 (estimated from the volume of the lava field and the duration of the eruption) and peak rate larger than 200 m 3 s −1 on April 6 (a lower bound estimate due to hidden flow within lava tubes). Such emission rates are exceptionally large compared to classical emission rates of Piton de la Fournaise (Hibert et al., 2015) and are of the same order as effusion rates observed on other volcanoes when the caldera collapse influenced the emission rate (e.g., Gudmundsson et al., 2016).
As reported for other caldera collapses (e.g., Shuler et al., 2013), the duration of individual collapse events are significantly longer than what is typically observed for earthquakes of similar sizes. Our duration estimates ranges from 6 s to 14 s, while standard scaling laws for earthquakes (e.g., Duputel et al., 2013) predict durations ranging from 0.6 to 2 s for M w = 4 and M w = 5 events, respectively (cf., Supplementary Tables S4-S6). In addition, the source duration of VLP events does not seem to scale with their size. This clearly suggests that the underlying physical mechanism is not the same as earthquakes that are commonly observed on tectonic fault systems. Using a simple piston spring-block model (Kumagai et al., 2001), the duration of collapse events is controlled by the geometry of the subsiding caldera block but also by the properties of the magma reservoir. In this model, the weight of the collapsing piston is balanced by friction on the surrounding ring fault and by pressure in the magma reservoir. When the magma reservoir depressurization exceeds static friction, the piston moves down into the chamber which suddenly increases its internal pressure and generate VLP signals. Assuming that the inter-event timedelay T is large compared to the duration of individual collapses, we can write (see Supplementary Text S5): where t h is the piston source half-duration, V 0 is the initial volume of the magma reservoir and κ is the bulk modulus of the magma.
Using geodetic estimates of the magma chamber volume (Peltier et al., 2008), we can interpret our estimates of t h for each collapse to evaluate κ. Despite the large uncertainties, results shown in Fig.  S13 typically ranges from κ = 10 8 Pa to κ = 10 10 Pa, which indicates that the magma in the reservoir was in a bubbly state during the caldera collapse (Huppert and Woods, 2002). The increase of source durations observed from April 6 to April 7 in Tables S4-S6 suggests a reduction of the magma bulk modulus implying that the magma becomes more bubbly. This might correspond to the exsolution of gas following the depressurization of the chamber when a significant increase in the magma outflow rate is estimated (cf., Fig. 6b). Such interpretation seems consistent with field observations indicating a dramatic increase in the height of lava fountains to more than 200 m simultaneous with the reduction of κ and the increase of magma extrusion rate . The magma bulk modulus might also be impacted by a possible input of deep magma (as suggested by Fontaine et al., 2019). Such variation of the magma bulk modulus due to the influx of gas-rich magma has been proposed during the 2000 Miyake-jima collapse by Michon et al. (2011). An interesting aspect of the 2007 Piton de la Fournaise collapse is the evolution of the time-interval T between successive collapse increments shown in Fig. 6c. As reported previously (e.g., Michon et al., 2007), T gradually decreases from about 2 h at the beginning of the collapse to ∼30 min on April 6 at 12:00 GMT and then increases again until the end of the collapse. As detailed in text S5, the time-interval in our piston spring-block model can be expressed as: where μ is the difference between static (μ S ) and dynamic (μ D ) friction coefficients and F n is the effective normal force exerted on the ring-fault delimiting the piston cylinder. Although the decrease of κ due to the exsolution of gas will tend to increase the time-interval between collapse events, the estimated variations of κ seems to have a moderate effect on T (as also inferred by Michon et al., 2011). On the other hand, we see in Fig. S14 that the variations of T can be partly explained by the increase of the magma outflow rate α. As noted by Stix and Kobayashi (2008) and Michon et al. (2011), the variations of time-interval are also controlled by changes in frictional resistance (i.e., μ F n in eq. (2)). Fig. 6d shows the evolution of F = μ F n estimated from collapse event durations and vertical displacements (see eq. (23) of the Supplementary Text S5). Fig. S15 shows that similar estimates of μ F n can be obtained if we use the time-interval T and magma outflow α from seismic amplitudes (see eq. (25) of the Supplementary Text S5). These results suggest that the initial evolution of T can be explained by a decrease of the friction parameter μ F n at the beginning of the collapse. Consistently with Michon et al. (2011), our results also indicate that larger values of T at the end of the sequence corresponds to an increase of μ F n during the last collapse events (see Fig. S15).
As proposed by Michon et al. (2011), such a transitory decrease of frictional resistance might correspond to the reduction of the effective normal force F n caused by upward migrations of hydrothermal fluids. Such weakening could also be caused by frictional melting or the activation of magma pockets intruding in the ring-fault system. Recent observations of exhumed caldera faults suggest that frictional melting play a critical role in caldera collapses (e.g., Han et al., 2019). Moreover, seismicity and field observations at Piton de la Fournaise show that the ring-fault system is a preferential path of magma intrusions (e.g., Staudacher, 2010;Duputel et al., 2019). Both interpretations seem plausible given satellite observations showing a ring-shaped thermal anomaly in the Dolomieu crater during summit subsidence and the existence of small lava flows emitted from the edge of the collapsed piston on April 6 (Urai et al., 2007). As mentioned above, the surface collapse activity was preceded by a series of VLP precursors visible in Fig. 7. As suggested by (Fontaine et al., 2019), these events might correspond to deep collapses at depth before the onset of the surface subsidence. Such precursory VLP activity has also been observed before the 2000 Miyake-jima caldera formation (Kobayashi et al., 2012). These precursors have a longer duration that the VLP signal of the main collapse event (cf., Fig. S2b). Assuming the same source geometry, this could be explained by a larger centroid time difference between the vertical contraction of the magma reservoir and the resulting collapse of the overlying fractured column. Such larger time difference could potentially result from the fact that most of the piston remains locked during the early stage of the sequence. After the main collapse, the piston falls down more easily while the frictional resistance is reduced as shown in Fig. 6d. To a smaller extent, we notice some differences between the main collapse and the following smaller collapses in Fig. S2a. There are also discrepancies in VLP signals emitted at the end of the collapse sequence (e.g., Fig. S8-S10). These discrepancies can be explained by differences in focal depth and centroid time (cf., Table S4-S6) that can possibly be linked with the weakening of the edifice. However, such interpretations remain uncertain given our limited sensitivity to depth and the availability of a single broadband station to characterize small VLP events before and after the main collapse.
The 2007 Piton de la Fournaise summit subsidence outlines the hazard caused by caldera collapses on basaltic volcanoes. In this context, a parameter of primary importance is the critical volume fraction of evacuated magma to trigger the collapse. According to the volume change derived after the main collapse, this critical volume fraction is only a few percent (between 1% and 4% of the shallow magma reservoir volume). This corresponds to critical underpressures in the magma chamber ranging from -5 to -70 MPa, which is in rough agreement with geodetic estimates (ranging from −5 to −10 MPa; Peltier et al., 2009;Got et al., 2013). Previous estimates by (Michon et al., 2011) yielded a larger critical volume fraction of 8.2% corresponding to larger critical underpressures (from −45 MPa to −150 MPa). It is interesting to notice that the 2007 Piton de la Fournaise collapse has characteristics that are similar to other basaltic caldera forming eruptions. The 1968Fernandina, 2000Miyakejima and 2018 Kilauea collapses are characterized by quasi-periodic series of incremental subsidence events manifested by tilt steps and VLP signals similar to what is observed at Piton de la Fournaise. These events are also driven by magma withdrawal from shallow reservoirs through low-elevation eruptive vents. For both Fernandina and Kilauea, the time interval T between successive collapses follows an evolution that is quite similar to what is observed at Piton de la Fournaise (Stix and Kobayashi, 2008;Michon et al., 2011;Neal et al., 2019). Despite these similarities, there are also significant differences. For example, the total duration of the whole caldera collapse ranges from about 2 days at Piton de la Fournaise to more than 3 months at Kilauea. Some of these discrepancies can possibly be explained by differences in the summit caldera block geometry, reservoir volume, edifice strength, magma properties and lava extrusion rates. The relative importance of these different parameters are still poorly constrained given the challenging circumstances in which such collapses have been observed. New near-field seismogeodetic observations of the recent Kilauea summit subsidence therefore appear promising to better understand the mechanisms at the origin of caldera collapses.