Constraining Spatiotemporal Characteristics of Magma Migration at Piton De La Fournaise Volcano From Pre-eruptive Seismicity.

Volcano-tectonic seismicity has been recorded for decades on various volcanoes and is linked with the magma transport and reservoir pressurization. Yet earthquakes often appear broadly distributed such that magma movement is difficult to infer from its analysis. We explore the seismicity that occurred before eruptions at Piton de la Fournaise in the last 5 years. Using template matching and relocation techniques, we produce a refined image of the summit seismicity, demonstrating that most earthquakes are located on a ring structure. However, only a portion of this structure is activated before each eruption, which provides an indication as to the direction of the future eruptive site. Furthermore, we show that the delay between the pre-eruptive swarm and the eruption onset is proportional to the distance of the eruptive fissures relative to the summit cone. This reveals that the beginning of the intrusion already bears information regarding the future eruption location.

Introduction

One goal of volcano monitoring is to provide reliable information concerning an impending eruption. With this purpose, one generally relies on the detection of anomalous signals that could reflect magma intrusions within the edifice. A major challenge is the difficulty in interpreting the different observations that are recorded during the pre‐eruptive phase. Volcano‐tectonic events are often considered to signal dynamic shear failure within the edifice. It is envisioned that a migrating magma pressure source can locally create high‐stress concentrations that ultimately result in material failure and microearthquakes. Such a relation between the propagation of a magma‐filled dike and the occurrence of earthquakes has been documented on several occasions (e.g., Battaglia et al., 2005; Einarsson & Brandsdóttir, 1980; Grandin et al., 2011; Rubin et al., 1998; Sigmundsson et al., 2015). However, it is often not easy to relate the recorded seismicity to a moving magma source. In certain circumstances, a dike can even intrude the volcanic edifice with no detectable earthquakes (e.g., Aki & Koyanagi, 1981; Wright et al., 2006). Indeed, as noted by Rubin and Gillard (1998), the occurrence of earthquakes large enough to be detected most likely corresponds to the failure of pre‐existing weak structures within the edifice in the vicinity of the propagating crack but not necessarily the exact tip of the dike. Furthermore, the final stage of the magma ascent is often accompanied by very low seismic activity only comprising low‐amplitude events that are difficult to track and locate (Eibl et al., 2017; Taisne et al., 2011; Zecevic et al., 2013). This is because the low stresses at shallow depth only permit very small earthquakes that cannot be easily detected and located. In addition, upper parts of volcanic edifices are usually highly heterogeneous, which can prevent long faults to develop. Therefore, as shallow magma chambers are usually found on basaltic volcanoes, the last stage of magma propagation between such reservoirs and the surface is difficult to track. It thus limits our ability to make predictions of the future surface eruptive activity based on seismicity recording.

In this study, we focus on the Piton de la Fournaise volcano located on La Réunion island in the Indian Ocean. This frequently active (i.e., 13 eruptions between June 2014 and July 2018; Table 1) basaltic volcano represents an excellent field laboratory that is monitored by a dense network of geodetic and seismological stations (see Figure 1; Brenguier et al., 2012; Staudacher & Peltier, 2016). Being able to predict the location of a future eruption is of great interest at the Piton de la Fournaise. Indeed, a major concern are distal flank eruptions close to inhabited areas outside of the Enclos Fouqué caldera. In addition, as the volcano is a popular tourist attraction (with over 350,000 visitors per year; Harris et al., 2017), the prediction of eruptive fissure location is important to plan the evacuation of visitors hiking into the caldera. Eruptions are generally preceded by an intense shallow seismic crisis lasting a few hours (Roult et al., 2012), which originates from a persistent structure located below the volcano summit around sea level (Lengliné et al., 2016; Sapin et al., 1996). Some shallow events are occasionally captured in the latter stage of the crisis, but they are usually difficult to locate (e.g., De Barros et al., 2013; Taisne et al., 2011). In fact, most events located before eruptions remain confined below the central cone and do not track the last stage of magma movement. Here we use a template matching technique to detect and relocate seismic events during the 13 last eruptive unrests at Piton de la Fournaise volcano. We show that a careful analysis of the pre‐eruptive seismicity provides clear indications as to the location of an oncoming eruption.

Table 1

Date, Time, Duration, and Location of the Analyzed Eruptions

Date	Time	Duration	Location
20/06/2014	21:35	19 hr	Summit area
04/02/2015	07:00	11.5 days	Summit area
17/05/2015	09:45	13.3 days	Southeast flank
31/07/2015	05:20	2.1 days	North flank
24/08/2015	14:50	58.3 days	Southwest flank
26/05/2016	04:05	26 hr	South flank
11/09/2016	04:41	6.8 days	North flank
31/01/2017	15:40	27.0 days	South flank
17/05/2017	16:10a	4.0 days	North flank
13/07/2017	20:50	45.1 days	South flank
03/04/2018	06:40	17 hr	North flank
27/04/2018	19:50	34.6 days	Southwest flank
12/07/2018	23:30	19 hr	Northwest flank

Note. Eruptive fissures are mapped in Figures 1, 2, and S3.

Only gas and no magma was associated with the eruption on 17/5/2017.

Figure 1

Relocated earthquake catalog, eruptive fissures and seismic stations. (a) Map view of the 15,598 relocated earthquakes (black circles) and volcanic fissures of eruptions between June 2014 and July 2018 (red lines). Orange triangles indicate seismic stations of the Piton de la Fournaise observatory. (b) Map view of relocated earthquakes in four different depth intervals (“elv” stands for elevation above sea level in meters). (c) Projection of relocated events onto an east‐west cross section. Topography at latitude −21.45° is shown for reference. We observe in (b) and (c) that deeper seismicity focuses on the east, while shallower events are mostly located on the west. The overall structure has an eastward dipping angle of ∼20°.

Figure 2

(a–l) Location of earthquakes and eruptive fissures. Each map presents the seismic crisis preceding an eruption (cf., date indicated on top of each subfigure). The location of eruptive fissures is indicated in red. Black circles represent the seismic events that occurred within 24 hr before the eruption. Gray circles represent the entire set of relocated seismic events. Fissure locations are based on visual inspections from field investigations with GPS tracking and photogrammetry. The same map is presented in Figure S3 for the eruption of 17/5/2017 that was only associated with gas emissions (no lava was ejected from the fissure). We also compare the azimuthal distribution of earthquakes and eruptive fissures in Figure S4.

Data and Methods

Earthquake detection is conducted using a template matching procedure following the same approach described in Lengliné et al. (2016). This detection process has been directly implemented at the Piton de la Fournaise Volcano Observatory to provide an automated detection of any seismic events. Each day, we build a template earthquake catalog from the last 1,000 events manually located in the shallow part of the volcano. Among these 1,000 template candidates, we only select events that have a depth uncertainty lower than 340 m and be located above sea level in the Enclos Fouqué caldera. We then scan the previous 24 hr of continuous seismic data to detect events that match at least one of the events in this template list. The earthquake detection is performed using the vertical components of at least three stations near the central cone (BOR, BON, SNE, FJS, RVL, or CSS). All traces are filtered between 5 and 25 Hz. For each template, we compute correlation coefficients for each time step (0.01 s) on all selected channels. We then apply a maximum filter over a duration of ±0.1 s and stack the correlation traces after applying a travel time correction. We consider all time windows where the average correlation over the last three traces is larger than 0.45 as possible detections. In the case when multiple templates lead to a common detection, we only consider the detection associated with the template that produced the highest correlation coefficient with the continuous data stream. Once detection times have been obtained, we extract the waveforms at all of the stations in the Piton de la Fournaise Volcano Observatory network for each newly detected event based on the template travel times. This results in a total set of 16,246 events between June 2014 and July 2018. This compares to the 1,635 events that are listed in the original catalog for the same time period. Almost all of these events occur during the seismic crises that precede the eruptions listed in Table 1.

In order to obtain the precise location of each event, we use a double‐difference relative location algorithm (Waldhauser, 2000). We first compute P and S wave differential travel times for all pairs of detected events on 1.28‐s‐long windows based on the maximum of the cross‐correlation function. The differential travel times of P wave are computed on vertical components, whereas the differential travel times of S waves are obtained on both horizontal components. Following Carrier and Got (2014), we apply an a priori Cauchy distribution with σ 0=0.1 s to this cross‐correlation function in order to mitigate the contamination of outliers in the data set. We discard travel time delays when the correlation coefficient multiplied by the a priori function is lower than 50%. All events are finally relocated together which results in a final set of 15,598 relocated events. We use a simple homogeneous velocity model for the relocation with a P wave velocity of 3.0 km/s and a V /V ratio of 1.73. We also account for the station elevations when computing raypaths. Events that are discarded during the relocation process correspond to earthquakes that have too few links with other events to be adequately relocated.

Results

The entire set of relocated earthquakes is presented in Figure 1. We observe that almost all events focus on a ring structure located below the summit of the volcano between a depth of 200 m and 1 km above sea level. This structure was already evidenced by Lengliné et al. (2016) using data from the 2014 and 2015 eruptions. It delineates a persistent feature of the volcanic edifice with a dip of 20° toward the east (Figure 1). Its location is considered to correspond to a zone of weakness or stress concentration at the roof of the magma reservoir (located around sea level according to previous studies; Lénat et al., 2012; Peltier et al., 2005). Our updated relocated catalog seems to show more details than previous studies (e.g., Lengliné et al., 2016; Sapin et al., 1996): we notice secondary seismicity alignments inside the main ring structure such as a N45° lineation that is clearly visible to the northwest of the Dolomieu crater in Figure 1b. This structure certainly corresponds to a ring fault system hosting repetitive collapses of the volcano summit (e.g., Staudacher, 2010). Alignments of seismicity inside the main ring structure might outline secondary faults or weak structures associated with the different collapsed areas that affected Dolomieu and Bory summit craters in the past (Michon et al., 2013). Earthquake swarms that occur on the ring structure before each eruption seem to have distinct locations. Indeed, the correlation matrix highlighting the similarity of the recorded signals shows that events are mostly correlated with earthquakes belonging to the same pre‐eruptive swarm (cf., Figure S1 in the supporting information). One consequence of this poor correlation among seismic swarms is that most earthquakes are detected shortly before or after their template events (cf., Figure S2).

To analyze the pre‐eruptive activity, we focus on events that occurred within 24 hr before each of the 13 eruptions that occurred between 2014 and 2018 (see Table 1). Eruptions at Piton de la Fournaise are typically preceded by a seismic crisis a few hours before the eruption onset. We report in Figure 2 the location of earthquakes preceding each eruption along with the corresponding eruptive fissures. For all eruptions, we readily observe that only a portion of the ring structure is activated during the pre‐eruptive seismic crisis. Furthermore, we clearly see that the activated portion of the ring provides a good indication of the future eruptive fissure location. More specifically, the azimuth of the eruptive site relative to the summit area can be extrapolated from the location of earthquakes on the main ring structure. We can see in Figures 2, S3, and S4 that we have a perfect match for all 13 analyzed eruptions. It suggests that, although individual events do not allow us to track the magma propagation all the way up to the surface, there is a clear link between earthquake locations on the ring structure and the existence of a moving magma body.

To investigate how the pre‐eruptive seismicity is connected to the mass magma transport within the edifice, we measure the timing of the seismic crises relative to the eruption onsets. Pre‐eruptive seismic crises are conspicuous on the seismograms presented in Figure 3: eruptions are systematically preceded by a significant increase of seismicity for a duration that varies significantly between each eruption (from a few tens of minutes to more than 8 hr). We can also note that the earthquake rate drops before the onset of the eruptions. As mentioned before, this suggests that the ultimate propagation of the magma up to the eruptive site is mostly aseismic, which makes it difficult to track using earthquakes. Some events can still be detected during this phase (see Figure 2), but they are difficult to locate because they are usually visible only at a few stations (e.g., De Barros et al., 2013; Taisne et al., 2011). We can however investigate if there is any relationship between the timing of the seismic swarms and the location of the eruptive sites. Figure 3 shows the pre‐eruptive time delay of the seismic swarms as a function of their distance to the eruption sites. Pre‐eruptive time delays correspond to the time difference between the beginning of the seismic crisis and the eruption onset. They are computed using a simple approach, considering that the seismic crisis starts when the earthquakes rate reaches a threshold of 150 events per hour. The distance to the eruption site is defined as the shortest distance between earthquakes occurring at the beginning of the seismic crisis and the eruptive fissures. There is a clear correlation between the time delay and the distance of the seismic swarms to the eruption sites. This simply indicates that an eruption occurring at a large distance from the magma reservoir will also occur at a large delay after the onset of the seismic crisis. This is consistent with previous studies showing that longer seismic crisis are usually associated with eruptions at larger distances from the summit (Aki & Ferrazzini, 2000; Peltier et al., 2005; Roult et al., 2012). A rough linear fit on the observed scaling allows us to constrain a migration speed for the magma of the order of 0.1 m/s.

Figure 3

Time evolution of pre‐eruptive seismic crises. (a) Pre‐eruptive swarm time delay as a function of distance between earthquakes and eruptive fissures. The pre‐eruptive time corresponds to the time difference between the swarm onset and the eruption onset. (b–d) Examples of continuous seismograms on the vertical component recorded before three eruptions at station BOR, located on the southwest rim of the Bory crater (cf., Figure 1). The continuous red line refers to the seismic crisis onset. The dashed red line at time 0 indicates the start of the eruption.

Discussion

Seismicity is not always a straightforward indicator of magma propagation. Although volcano‐tectonic earthquakes can sometimes be interpreted as a direct marker of dike pathways (e.g., Rubin et al., 1998; Sigmundsson et al., 2015), it can also correspond to distributed damage caused by pressure variations in the magma reservoir (Carrier et al., 2015) or alternatively as the result of pre‐existing suitably oriented structures already near to failure (Rubin & Gillard, 1998). At Piton de la Fournaise, almost all pre‐eruptive seismicity is focused on the same ∼20° eastward dipping ring‐shaped cluster below the summit and above a shallow magma reservoir located around sea level (see Figure 1; Peltier et al., 2009). This cluster certainly represents an area of weakness outlining the piston structure associated with the repetitive formation of pit craters and caldera collapses in the volcano summit (Lénat et al., 2012; Michon et al., 2013). This interpretation seems consistent with the existence of a shallow low‐density fractured column beneath the central cone that has been inferred from gravity measurements (Gailler et al., 2009). Such a weakened zone that is prone to failure can facilitate the propagation of magma intrusions. These intrusions can in turn maintain damages within the edifice due to the associated stress changes and induced hydrothermal alteration (Heap et al., 2015; Mordensky et al., 2018; Pola et al., 2012; Wyering et al., 2014).

Although pre‐eruptive swarms occur consistently on the main ring structure, our results show a clear correlation between earthquake locations and the azimuth of eruptive sites relative to the central cone (Figure 2). In addition, the delay between the seismic crisis and the eruption onset correlates relatively well with the distance of the eruption fissure relative to the magma chamber (Figure 3). To gain further insights into this pre‐eruptive activity, we also investigate the evolution of earthquake locations during each eruptive unrest. Interestingly, we notice a systematic migration of the earthquake activity during the seismic crises (Figures 4 and S5). These seismicity migrations correspond primarily to an east to west propagation at a speed of the order of 0.1 m/s. Westward propagations are usually associated with an upward seismicity migration at a similar speed (see right subplots in Figure 4). This earthquake migration speed is consistent with the overall magma migration velocity inferred in section 3. Although the horizontal direction of earthquake migrations is sometimes consistent with eruptive fissure locations (e.g., Figure 4b), the evolution of seismicity is often not in agreement with the final location of the eruption site (e.g., Figure 4a where eruptive fissures are located at the southeast of a westward migrating swarm).

Figure 4

Pre‐eruptive earthquake migrations. (a) Migration observed before the 17/5/2015 eruption. (b) Migration observed before the eruption on 27/4/2018. Earthquakes are color coded as a function of time before eruption in the map showed on the left. Subfigures on the right depict earthquake coordinates as a function of time. Black arrows in (b) indicate a shallow cluster of seismicity that occurred close to eruptive fissures about 40 min before the eruption on 27/4/2018. Similar maps are presented in Figure S3 of the supporting information for the 11 other eruptions analyzed in this study.

A simple interpretation of these results is that the observed seismicity responds to blade‐like dike intrusions propagating vertically from the roof of the magma reservoir. We can suppose that only the weakened part of the edifice located close to the intruding dike will be activated seismically. The upward seismicity migration from east to west could then be attributed to the eastward dip of the ring structure. Indeed, because of this structure geometry, an upward migrating stress source will cause the failure of events located more and more to the west. Once this intruding dike has inflated vertically, its lateral propagation direction is mostly controlled by its location relative to the volcano summit (i.e., the stress field created by the topography) and the existence of preferential intrusion paths (e.g., rift zones; Michon et al., 2015). This interpretation is consistent with the analysis of geodetic time series that evidenced two phases of dike propagation: first, a vertical migration when the magma leaves the reservoir and then a lateral migration in the direction of the eruptive site (Peltier et al., 2005, 2007).

Our findings suggest that a careful analysis of the seismicity during the onset of a magma intrusion provides direct information on the future eruption location. More specifically, at a given time after the beginning of a seismic crisis, we can predict the azimuth and minimum distance of the eruptive site relative to the summit area. After a seismically active phase at the beginning of the intrusion, the magma propagates almost aseismically up to the surface. Although seismicity is very low after the seismic crisis, a few events can still be detected during this second phase. This is illustrated in Figure 4b where a shallow seismicity cluster is activated close to the final eruptive site 40 min before the eruption (cf., black arrow in Figure 4b). Our tests show that including more stations away from the central cone does not allow to increase the number of detections at shallow depth. This is mainly due to the limited number of templates located out of the ring structure. Shallow events are difficult to locate because they are usually observed only at a few stations with small signal‐to‐noise ratio (hence, they cannot be used as templates in our detection and relocation procedure). The development of a dedicated approach tracking the location of such shallow microearthquakes is of major interest at Piton de la Fournaise, where the prediction of the location of an eruption is important to plan rescue and evacuation operations. Some approaches have been proposed to get approximate locations based on time‐averaged seismic amplitudes or amplitude ratios (De Barros et al., 2013; Taisne et al., 2011). Beyond the implications for monitoring purposes, the development of new approaches enabling accurate localization of individual induced microearthquakes is important to improve our understanding of the dynamics of magma propagation.

Supporting Information S1

Click here for additional data file.