Absence of Detectable Precursory Deformation and Velocity Variation Before the Large Dome Collapse of July 2015 at Volcán de Colima, Mexico

Introduction

Volcanic eruptions result from complex processes that include feeding of magma storage zone, magmatic intrusion, interaction with surrounding rock and hydrothermal system, or changes in the physical state, chemical and mineralogical content of reservoirs and conduits. Many of these processes produce phenomena that can be observed at the free surface of the edifice before an eruption. They are thus considered as precursors and volcanologists use them for forecasting volcanic eruptions. Because these phenomena are the basis of volcano monitoring, many volcanological research aims at detecting and interpreting them. The most widely used methods are the study of the seismic activity, the measurement of ground deformations and the analysis of gas flow and composition (Scarpa and Gasparini, 1996). Other approaches, such as magnetic and electric studies can complement the classical ones. For example perturbations of magnetic field and variations of self-potential anomalies have been observed prior to a few eruptions (Johnston, 1997; Zlotnicki et al., 2009).

Several types of seismic precursory phenomena can be observed (McNutt, 1996). They include the increase of the level and energy of seismic signals, the rise of the number of events, the migration of the seismogenic zones, variations of the focal mechanisms of volcano-tectonic (VT) events, changes in anisotropy of the seismic velocities, emergence of swarms of various types of event. By analyzing earthquakes with similar waveforms (or multiplets – Poupinet et al., 1984) or ambient noise correlation functions, variations of seismic velocity in the medium have been detected before some eruptions. In a pioneering work, Ratdomopurbo and Poupinet (1995) used seismic multiplets and a technique known as Coda Wave Interferometry (Grêt et al., 2005; Snieder, 2006), to detect a velocity increase of 1.2% several months before the 1992 Merapi eruption. At Piton de la Fournaise, La Reunion Island, Brenguier et al. (2008) demonstrated, using noise correlation functions, that the velocity decreased by about 0.05% a few weeks before several eruptions in 1999 and 2000. Since then, many studies of velocity variations at Piton de la Fournaise have been published, extending the period of analysis (Duputel et al., 2009; Clarke et al., 2013; Rivet et al., 2014), improving the technique (Clarke et al., 2013; De Plaen et al., 2016), and locating the source of perturbation (Obermann et al., 2013). Velocity variations have also been detected before some eruptions of stratovolcanoes, such as Ruapehu, New Zealand, (Mordret et al., 2010), Mt Asama, Japan, (Nagaoka et al., 2010), Miyakejima, Japan, (Anggono et al., 2012), Etna, Italy, (Cannata, 2012), Mt St Helens, USA, (Hotovec-Ellis et al., 2015) or Merapi, Indonesia, (Budi-Santoso and Lesage, 2016), and shield volcanoes, such as Kilauea, Hawaii (Donaldson et al., 2017). Sharp velocity decreases induced by large tectonic earthquakes have also been observed at some volcanoes (Nishimura et al., 2000, 2005; Battaglia et al., 2012; Brenguier et al., 2014; Lesage et al., 2014).

Deformation data have also proven useful to reveal magma accumulation inside crustal storage zones or emplacement at shallow depth before an eruptive event (Dvorak and Dzurisin, 1997; Dzurisin, 2003, 2007), such providing both long (months to years) and short (hours to days) term precursors of volcanic activity. Deformation before eruptions are observed both on basaltic volcanoes (Sturkell et al., 2006) and andesitic and rhyolitic systems (Swanson et al., 1983). In recent years, the number of volcanoes where deformation data are available have drastically increased thanks to satellite radar interferometry, which provides high spatial resolution surface displacement with a precision of a few millimeters to a few centimeters depending on the available dataset (Biggs et al., 2014; Pinel et al., 2014; Biggs and Pritchard, 2017). This large and worldwide amount of available data has enabled to statistically question the link between deformation and eruptive activity.

Because the quest for eruption forerunners is of upmost importance, many scientific papers in this field present observations of precursory phenomena, investigate their interpretation and relationship with magmatic and hydrothermal processes, and discuss their use for forecasting. However, volcanic eruptions are not always preceded by precursors. There are cases where no precursory, or even co-eruptive, deformation was observed (e.g., Moran et al., 2006; Chaussard et al., 2013; Ebmeier et al., 2013; Biggs et al., 2014) or, more generally, where the pre-eruptive phenomena were too small to be detected and interpreted correctly for predictions. Such cases have been described at Popocatepetl, Mexico (Quezada-Reyes et al., 2013); Ontake, Japan (Kato et al., 2015), Soufrière Hills, Montserrat (Calder et al., 2002), and many other places. Most of the eruptive crises characterized by lack of precursors are probably not reported in scientific journals, although this information may be partly found in bulletins of volcanological observatories. This prevents from evaluating on a large scale the proportion of eruptions that are preceded, or not, by precursory phenomena.

The forecast of volcanic eruptions can follow two complementary approaches. The probabilistic one aims at estimating the probability that an eruption occurs in a given time interval either at short or at long term (Marzocchi and Bebbington, 2012). The deterministic approach tries to estimate the time of occurrence of an eruption. The most widely used method for the later approach is the material Failure Forecast Method (FFM–Voight, 1988) which gave encouraging results in a few cases (e.g., Cornelius and Voight, 1995; Kilburn and Voight, 1998; De la Cruz-Reyna and Reyes-Dávila, 2001). Nevertheless, eruption forecasting is partly an empirical task based on the knowledge of the previous volcanic activity, on the observations produced by the monitoring system, and on the experience of the volcanologists in charge. Forecasting methods still require to be improved in order the predictions to be more reliable and precise. Their performance and success rate need to be evaluated in a large variety of cases, including different types of volcanoes and various kinds and amplitudes of eruptions. In this evaluation, it is important to take into account the eruptions that were not preceded by precursors, as their occurrence can produce hazardous situations for inhabitants and visitors. Thus, published papers should also document volcanic crises characterized by lack of precursors as well as observations of forerunners that are difficult to interpret and to use for forecast. The difficulty in estimating the performance of forecasting methods is also illustrated by the fact that no more than 20% of the eruptions were anticipated by alert level changes (Winson et al., 2014). Furthermore, there is a tendency in volcano observatories to underreport intrusive episodes not leading to eruptions (Moran et al., 2011), even if some recent databases are trying to include these cases (Ebmeier et al., 2018).

In the present paper, we look for precursors of the large dome collapse of Volcán de Colima in 2015 by carrying out two approaches. First we apply techniques of seismic ambient noise correlation to estimate velocity variations in the structure by coda wave interferometry. Then we track potential deformation of the summit part of the volcanic edifice using SAR images. We show that there were no significant precursory signals in deformation and velocity variation before this major eruption, we interpret these observations from a volcanological point of view and we discuss the significance of the lack of precursor in term of hazard management.

Volcanological Setting

Volcán de Colima is located in western Mexico (19.51°N, 130.62° W, 3,860 m asl) and is currently one of the most active volcanoes in North America. It has produced at least three sub-Plinian to Plinian eruptions in 1576, 1818, and 1913, and an average of one large magnitude eruption per century (Robin et al., 1987; Luhr and Carmichael, 1990). The 1913 eruption included an opening phase producing pyroclastic flows and surges, a vent clearing phase which destroyed the summit dome, and a Plinian phase with a 23 km high column, the collapse of which generated a surge and 15 km long pyroclastic density currents (PDC). A Volcanic Explosivity Index (VEI) of 5 was calculated by Simkin and Siebert (1994) and Saucedo et al. (2010). The later authors estimated that a similar Plinian eruption would threaten more than 300, 000 inhabitants nowadays. Recent periods of moderate activity occurred in 1991, 1994, 1998–1999, 2001–2003, and 2004–2005 with alternation of dome buildings and destructions, Vulcanian eruptions, lava flows, and PDCs. The hazard associated with this activity and the observation of precursory phenomena triggered several evacuations of inhabitants of the most exposed localities (Macias, 1999).

The seismic activity of Volcán de Colima is monitored since the installation of a permanent network in 1989. Reyes-Davila and De la Cruz-Reyna (2002) analyzed the behavior of the Real-time Seismic Energy Measurement (RSEM) before several moderate eruptions during a period of about 10 years. They observed that, in most cases, no clear pattern of increasing seismicity can be detected. Before some events, the seismic activity presented clear acceleration only a few hours prior the eruption onset. The precursory activity of two eruptions in 1994 and 1998 was characterized by a clear acceleration of the energy release during several days. The authors applied a version of FFM to forecast the date of the eruption either in hindsight (1994 event) or in foresight (1998 event) (De la Cruz-Reyna and Reyes-Dávila, 2001). Boué et al. (2015, 2016) proposed a Bayesian approach of FFM which uses an automatic recognition system (Benítez et al., 2007; Cortés et al., 2014) to classify and separate different types of seismic event. They processed 13 years of continuous recording of Volcán de Colima. Among 36 mild explosions, 24 were preceded by accelerations of LP events rate, 7 by a linear increase, and 3 occurred without seismic precursors. Successful forecasts were obtained for one third of the cases, while this proportion reaches 83% when some reliability criteria are fulfilled. Lesage et al. (2014) calculated and analyzed the velocity variations for the whole period 1998–2013 using data from a pair of short-period stations. They did not find any clear relationship between velocity changes and the mild to large vulcanian eruptions that occurred during this interval. They only noted that most large eruptions coincided with periods of decreasing velocity.

Limited in situ measurements of ground deformation have been performed at Volcán de Colima volcano based on tiltmeters, precise leveling, EDM and GPS campaigns (Murray and Ramírez-Ruiz, 2002; Murray and Wooller, 2002; Zobin et al., 2013). In particular, an acceleration of the summit inflation measured by EDM before the 1998 eruption was interpreted based of the FFM to predict a posteriori the time of the eruption (Murray and Ramírez-Ruiz, 2002). The application of InSAR technique to retrieve and quantify the displacement field at Volcán de Colima is challenging due to the limited coherence on the volcano slope as well as tropospheric artifacts (Pinel et al., 2011). SAR studies have thus been restricted to the 5 kilometer wide summit coherent area. The time series analysis of ASAR-ENVISAT data recorded from mid-2002 to the end of 2006 evidenced a summit subsidence reaching a rate of around 1 cm/yr and centered on the summit but enhanced on recent lava flows, which was interpreted as due to eruptive deposits load effects associated to a shallow deflating source (Pinel et al., 2011). Using high spatial and temporal SAR data acquired by TerraSAR-X satellite, Salzer et al. (2014) were able to catch a localized pre-explosive deformation induced by a transient pressurization of the shallow plumbing system before an explosion which occurred in January 2013.

After a period of quiescence since June 2011, a new phase of activity initiated in January 2013 when several moderate Vulcanian explosions destroyed the lava dome emplaced in 2007–2011. Then lava extrusion occurred with rate of 0.1–0.2 m³s⁻¹, forming a new lava dome and producing rockfalls and explosions until February-March 2014. An increase of lava extrusion and rockfall activity was observed in September 2014 with extrusion rate of 1–2 m³s⁻¹. In January-February 2015, a series of explosions produced the destruction of the dome and the end of the rockfall activity. Another lava dome was observed on 20 May 2015 and continued to grow at low rate. On 3 July, a moderate explosion occurred and was followed by a decrease of the explosive activity and an acceleration of the rockfall activity. More information on the volcanic and seismic activity of Volcán de Colima in 2013–2015 can be found in Zobin et al. (2015), Capra et al. (2016), Reyes-Dávila et al. (2016), Macorps et al. (2018), and Arámbula-Mendoza et al. (2018). On 10 July 2015, after 2 days of increased extrusion rate, a partial collapse of the dome occurred accompanied by large PDCs that reached 9.1 km on the south flank. On 11 July, 16 hours after the first events, another series of PDCs of larger size reached distances of 10.3 km from the crater. No eruptive columns were produced during the whole sequence. The total volume of material including PDCs and ashfalls was estimated to 14.2 × 10⁶ m³ (Reyes-Dávila et al., 2016). A more recent study estimated a volume of block-and-ash flow material of 7.7 ± 1 × 10⁶ m³, based on optical and field data (Macorps et al., 2018). This volume makes the July 2015 sequence the most important eruption since the Plinian event of 1913. In contrast with the common behavior of the volcano (Arámbula-Mendoza et al., 2011), the July 2015 eruption was preceded by a decrease of the rate of LP events and explosions. The main precursory phenomenon was a marked increase of the number and energy of rockfalls and PDCs that accompanied the rise of the extrusion rate (Reyes-Dávila et al., 2016; Arámbula-Mendoza et al., 2018). Although the deformations and the velocity variations associated with the previous eruptions were small or undetectable, we may expect to detect stronger forerunners in the case of the 2015 exceptional events.

Apparent Velocity Variations

Network, Data, and Method

We used data from the Volcán de Colima monitoring network which is part of the State of Colima's Seismological network (RESCO). In a first stage, it included 4 stations equipped with vertical SS-1 Ranger short-period seismometers and analog transmission. It records continuous signals since 1998. This network was completed in 2001 and 2007–2008 by 6 Guralp CMG-6TD broadband stations (Figure 1).

Figure 1. Map of the monitoring seismic network of Volcán de Colima (VC). EZV3-6 are short-period stations. Yellow lines indicates the pairs of broadband stations used to calculate the velocity variations. The red frame delineates the area displayed on the images of Figure 4.

We extracted Green functions between pairs of sites by cross-correlating ambient seismic noise recorded at the corresponding two stations (Weaver and Lobkis, 2001; Campillo and Paul, 2003). We followed a standard procedure which consists of trend and mean removal, band-pass filtering in the range [0.125–2] Hz, spectral whitening and one-bit amplitude normalization (Bensen et al., 2007; Lesage et al., 2014) and then calculating ambient noise cross-correlation functions (NCFs) for delays of ±150 s. For three Guralp stations (SOMA, WEST, INCA), we detected periods of several months where their clocks were not synchronized, with time lags of up to ±1.3 s. These periods are 12 March 2014 to 25 November 2014 for station SOMA, 8 October 2014 to 30 April 2016 for WEST, and 23 January 2015 to 30 June 2015 for INCA. They are indicated by dotted lines in Figure 2. These time drifts were corrected for stations SOMA and INCA by comparing the corresponding NCFs with those obtained when the clocks were well-synchronized (Stehly et al., 2007; Sens-Schönfelder, 2008). The time drift pattern for station WEST was more complex and corrections are not reliable enough. Thus, we did not use this station for the following analysis.

Figure 2. Apparent velocity variations estimated for pairs EFRE-INCA (red), EFRE-SOMA (black), MNGR-INCA (blue), SOMA-JUBA (green), and SOMA-MNGR (cyan). Average velocity variation for the available pairs (thick dark line). Vertical dashed green lines indicate main regional tectonic earthquakes and vertical dashed blue lines indicate main eruptions. Horizontal dotted lines at the bottom of the graphic indicate periods of bad clock synchronization for stations SOMA (red) and INCA (green). The corresponding clock delays were corrected using noise cross-correlation (see main text for details).

Daily NCFs were calculated for the vertical components of station pairs EFRE-SOMA, EFRE-INCA, MNGR-INCA, SOMA-MNGR, and SOMA-JUBA from January 2013 to April 2017. The paths between the stations belonging to these pairs go through the volcanic edifice (Figure 1). The corresponding NCFs should be thus affected by any perturbation in the volcano. Then we estimated the apparent velocity variation (AVV) by using the stretching method (Lobkis and Weaver, 2003; Sens-Schönfelder and Wegler, 2006). A NCF stacked over 2013 was used as a reference and compared to each daily NCF. The reference NCF is stretched or compressed in order to maximize the correlation coefficient between the coda of both NCFs in a delay range of [10, 80] s. Because the correlation functions obtained are asymmetrical, we used the side of the NCF which presented maximum amplitude. The corresponding stretching coefficient is equal to the negative of the apparent relative velocity variation (Lesage et al., 2014). Using this approach, we calculated time series of AVV for each station pair, as well as their average. In the Supplementary Material we compare this approach with two other methods and we show that similar results are obtained in all cases.

Results

Figure 2 displays the apparent velocity variations obtained for the 5 pairs of broadband stations from 2013 to 2017, together with their average. All the AVVs present short term (<month) fluctuations with amplitudes of 0.05–0.1%. The averaging of the AVVs reduces to less than 0.05% the amplitude of these fluctuations that can be considered as the noise level. Several tectonic earthquake with magnitude larger than 6 and epicentral distance smaller than 500 km from Volcán de Colima (vertical green lines in Figure 2) induced sharp velocity decreases. For example, a velocity drop of about 0.2% in average occurred on April 18, 2014 during a M = 7.1 earthquake located at 350 km from the crater. This phenomenon was reported by Lesage et al. (2014) who demonstrated that the corresponding perturbation is localized in the shallow layers of the edifice.

No clear velocity variations with amplitude larger than the noise level appeared before, during nor after the July 2015 large eruptions. A small oscillation of approximately ±0.07% can be observed during less than 2 months before the events, but it is poorly significant and its use as a precursory signal would not be reliable. A sequence of mild Vulcanian eruptions in January-February 2017 was neither associated with AVVs. In the Supplementary Material we present several arguments that support the reliability of the estimation of the velocity variations.

Summit Deformation Studied by InSAR

Data and Method

We used SAR images acquired in C-band by the European Satellite Sentinel-1A over Colima volcano since November 2014. The present study is based on 8 descending images of Track 12 (subswath 1, VV polarization, look angle over the volcano summit of 34.1°) recorded before the July 2015 eruption from November 27th 2014 until June 15th 2015 (Figure 3). Images were provided by the European Space Agency (ESA) as Single Look Complex (SLC) images and processed using the NSBAS chain (Doin et al., 2012) modified to integrate Sentinel-1 data acquired in TOPSAR mode as described in Grandin (2015). Topographic contribution was removed using the SRTM DEM at 30 m resolution. Tropospheric contributions were corrected using the ERA Interim global meteorological model provided by the European Center for Medium-Range Weather Forecast (ECMWF) as explained in Doin et al. (2009). Twenty-eight interferograms were calculated (see Figure 3 for the network) and unwrapped using the ROI_PAC branch-cut unwrapping algorithm. Unwrapped interferograms were geocoded on the 30 m resolution DEM. The phase delays of unwrapped interferograms were then inverted pixel by pixel using a least square method to solve for the cumulative phase delay through time. Pixels characterized by a large RMS (above a threshold of 0.9 rad) were discarded such that our displacement detection threshold on the remaining pixels can be estimated around 0.4 cm over the 6 month period studied.

Figure 3. Spatio-temporal distribution of Sentinel-1 data (descending track 12) used in this study. The effective spatial distance between orbits (perpendicular baseline) is represented as a function of the temporal distance (time gap) between the two satellite acquisitions. Red lines are for the interferograms computed and used for the time series inversion. Images were acquired on 23 November 2014, 17 December 2014, 10 January 2015, 11 March 2015, 4 April 2015, 28 April 2015, 22 May 2015, and 15 June 2015 at 12H50 UTC.

Results

The mean velocity of the ground surface in Line of Sight derived from the time series analysis is presented in Figure 4. The surface displacement can be retrieved from the phase of the radar signal only if the ground backscattering properties remain stable though time, insuring a good coherence. This condition is not fulfilled on the vegetated volcano slopes such that this SAR dataset only provides information in the summit area (at a distance smaller than 4 km from the volcano crater). Ash deposits resulting from explosions further limit the area with good coherence thus restricting the available information to the south-eastern part of the summit area at some distance from the dome. No information is available on the deformation of the dome itself. A subsidence signal (which appears in blue on Figure 4) is observed above and nearby the summit lava flows (Global Volcanism Program, 2015). Elsewhere no significant (superior to 6 ± 3 mm) Line of Sight displacement is evidenced over the period preceding the July 2015 eruption. Based on this observation, we can deduce that no vertical displacement larger than 1 cm can be seen over a distance of 4 km away from the summit crater. Considering a Mogi source (Mogi, 1958) localized beneath the crater, we can estimate a maximum value for the magma volume that can be stored elastically at depth below the crater without inducing detectable surface displacement. As shown in Figure 5, this value increases with the reservoir depth. We thus show that no significant volume of magma (above 3 million of m³) can have been stored shallower that 5.5 km depth during the 6 months preceding the eruption.

Figure 4. Mean velocity in the Line of Sight (LOS) obtained by time series analysis of the descending track images acquired between (A) November 2014 and May 2015 (7 images used) and (B) November 2014 and June 2015 (8 images used). The velocity is shown only for pixels where the phase information remains coherent through time and having a RMS smaller than 0.4 cm. Coherence is more reduced on interferograms calculated using the image acquired on 15 June 2015, which explains the difference in surface coverage on the two times series presented. The hatched square is the reference area where the displacement is set to zero. Cyan area is the surface affected by layover where SAR data cannot provide information due to the geometry of acquisition. Circles indicates the distance from the volcano summit. Dashed red lines correspond to the location of the lava flows causing the measured subsidences. The hillshade background was realized thanks to airborn LIDAR DEM acquired in 2004 by the Colima volcano observatory.

Figure 5. Maximum volume of magma stored as a function of the magma reservoir depth. Estimation is made considering a Mogi source centered below the crater and using the fact that no vertical displacement larger than 1 cm is observed from the crater toward a lateral distance of 4 km. Volume estimations of deposits are represented for comparison as horizontal black lines (plain line for Reyes-Dávila et al., 2016, dashed line for Macorps et al., 2018).

Discussion and Conclusion

Variations of seismic velocities have been observed as a precursory phenomenon of volcanic eruptions relatively recently thanks to the development of the techniques of continuous recording, of seismic multiplet analysis and of ambient noise cross-correlation. It is thus important to evaluate the frequency of its occurrence before eruptions and to investigate its origins and relationships with other observations. Several physical processes have been proposed to explain the velocity variations in volcanoes. Changes in the ground water level, related with precipitation, modify the pore pressure and can induce detectable velocity changes. This process is the source of seasonal effects that can be corrected when sufficient observations are available (Sens-Schönfelder and Wegler, 2006; Hotovec-Ellis et al., 2015; Rivet et al., 2015). Strong topographical changes related to caldera formation at Miyakejima and Piton de la Fournaise volcanoes were also accompanied by velocity increases or decreases (Duputel et al., 2009; Anggono et al., 2012; Clarke et al., 2013). Sharp temporary decreases in velocity have been observed in relation with strong ground shaking due to the passing of the seismic waves generated by large tectonic earthquakes (Battaglia et al., 2012; Brenguier et al., 2014; Lesage et al., 2014; this study). This phenomenon has been associated with the presence of highly pressurized hydrothermal or magmatic fluids at depth (Brenguier et al., 2014) or with mechanical softening and nonlinear elastic behavior of granular material in the shallow layers of volcanoes induced by the ground shaking (Lesage et al., 2014). The velocity variations observed before numerous eruptions of Piton de la Fournaise volcano have been related to overpressure induced by magma intrusion through a model of dilatancy and empirical laws linking perturbations of porosity and volume with changes in shear velocity (Brenguier et al., 2008). During the last days before the large 2010 eruption of Merapi volcano, rapid fluctuations of velocity have been interpreted as the consequence of the progressive fracturation and healing of the plug due to pulses of magma intrusion (Budi-Santoso and Lesage, 2016). However, in most studied cases, the velocity variations are interpreted as the result of the dependency of physical parameters of rocks to stress (Birch, 1960; Ratdomopurbo and Poupinet, 1995; Cannata, 2012; Hotovec-Ellis et al., 2015; Donaldson et al., 2017; Lamb et al., 2017). When an increasing effective pressure is applied to volcanic rocks, which are porous and pervasively microcracked, the most complient cracks and pores close yielding the elastic modulus and seismic velocities to increase (Vinciguerra et al., 2005; Stanchits et al., 2006; Nara et al., 2011; Heap et al., 2014). This process is involved in the increase of velocities with depth in volcanic structure (e.g., Lesage et al., 2018 for a review). When deviatoric stresses are applied, cracks normal to the axis of the maximum principal stress close, while those parallel to it remain almost unaffected (Nur, 1971). At high deviatoric stress level, a new population of cracks appears in the direction parallel to the maximum stress axis. This damaging process induces strong velocity decrease (Lamb et al., 2017). Thus there are complex relationships between changes of the stress field and the field of velocity variations.

This is supported by the coeval observations of ground deformations and velocity changes (Clarke et al., 2013; Rivet et al., 2014; Donaldson et al., 2017; Hirose et al., 2017). For example, Rivet et al. (2014) showed that the velocity decreases during long-term inflations of Piton de la Fournaise and increases during deflations of the edifice. Moreover a strong velocity decrease was observed at the time as a large movement of the East flank of this volcano before the major 2007 caldera collapse (Clarke et al., 2013). However, even in a simple elastic half-space, a source of increasing pressure produces both extensional strain in the region above the source and compressional strain in the surrounding volume (Pinel and Jaupart, 2003; Budi-Santoso et al., 2013; Got et al., 2013; Donaldson et al., 2017). Therefore, the relationship between pressure evolution in the magmatic system, ground deformations, and velocity variations may be relatively complex.

On the other hand, the amplitude of the velocity changes before eruptions are very small. In many studies, the published values, generally obtained by averaging the values estimated using many station pairs, are of the order of a few tenths of percent (Brenguier et al., 2008, 2016; Mordret et al., 2010; Donaldson et al., 2017). Only in a few cases, the relative velocity variation calculated using single station pairs or using seismic multiplets reaches about one percent (Ratdomopurbo and Poupinet, 1995; Nagaoka et al., 2010; Anggono et al., 2012; Hotovec-Ellis et al., 2015; Budi-Santoso and Lesage, 2016). Thus, the ability in detecting such small variations depends on the ratio of the amplitude of the signal of interest and that of the spurious fluctuations due to non-volcanic phenomena and to the non-stationarity of the sources of seismic noise (Stehly et al., 2007). It depends also on the relationship between the duration of the windows used to calculate the correlation functions of seismic noise and the characteristic time of the processes producing velocity variations.

At Volcán de Colima, although the 2015 dome collapse was the most important event in term of volume of emplaced material since the Plinian eruption of 1913, no clear signals of deformation and velocity variation could be detected before and during the eruptions, besides a good level of reliability. The main precursory phenomenon was the increase of rockfall activity that was interpreted as a consequence of accelerated extrusion rate and that was probably accompanied by strong dome modifications. Unfortunately, cloudy conditions due to rainy season prevented visual and photographic observations from quantifying the dome evolution.

In the study period 2013–2017, the amplitude of the common velocity fluctuations is about 0.05% for characteristic times less than a few months. Some sharp decreases of velocity followed by slow recovery are related with tectonic earthquakes, especially that of April 18, 2014. Other variations, such as that occurred on mid-September 2013, are not related to any known phenomenon. In May-June 2015, a sequence of increase-decrease-increase of the velocity coincides with the extrusion of a dome that was first observed on May 20 (Arámbula-Mendoza et al., 2018). However, the velocity recovered the value corresponding to its general trend by the first days of July. It is thus not possible to consider this May-June sequence as a direct precursor of the July dome collapse. The Green functions extracted by noise correlation are predominantly surface waves (Shapiro and Campillo, 2004). In the frequency range used to estimate velocity variations (0.125–2 Hz), the corresponding sensitivity kernels indicate that they are sensitive to velocity perturbations of the medium up to 2–3 km below the surface, i.e., at depth were possible sources of deformation are expected (Salzer et al., 2014).

Emplacement of magma inside the crust is generally expected to induce surface displacements (Dzurisin, 2007; Tibaldi, 2015). A volume of the order of the one emitted during the July 2015 event (14 million of m³) would require to have been emplaced either at some lateral distance from the crater (more than 5 km laterally away) or below the crater at more than 8.5 km depth in order to remain undetected in our dataset. If we instead consider the volume estimation for Block and Ash Flow deposits from Macorps et al. (2018), which is expected to be closer to the DRE volume, the minimum depth would rather be 7.5 km. This threshold depth is derived neglecting the magma compressibility, which may be large for bubble-rich magma and thus may reduce the surface displacement produced by magma emplacement (e.g., Johnson et al., 2000; Rivalta and Segall, 2008). Accounting for the effect of magma compressibility could increase the amount of magma possibly emplaced at shallow depth without significant surface displacement, such that the potential storage zone could be shallower but even so it would have to remain below 5.5 km depth in order to be consistent with the SAR dataset. Many eruptions occur without any detected surface deformation neither in the pre- or co-eruptive phase. For instance, Ebmeier et al. (2013) clearly evidenced a statistically significant lack of deformation for active volcanoes of the Central American Arc. The reason evoked for this absence of detectable surface deformation are eruptions fed directly by rapid magma ascent from deep magma storage zones, a diffuse and extended shallow storage system made of several vertically elongated cracks or large magma compressibility due to high volatile content. Chaussard et al. (2013) also described an absence of deformation during the eruptive activity at several volcanoes among which Volcán de Colima. They explained this behavior by an open system where the presence of a permanent conduit allows magma to rise toward the surface without pressurizing the reservoir.

The seismic activity of Volcán de Colima in 2015 is mainly composed of LP events, small high-frequency events only detected by the closest station to the summit, volcanic tremor, small explosions and numerous rockfalls (Reyes-Dávila et al., 2016; Arámbula-Mendoza et al., 2018). Almost no VT events are detected in this edifice. In the days preceding the dome collapse, a clear decrease of the number of explosions and LP events was observed as well as an accelerated rate of rockfalls generated by instabilities of the front of several lava flows. Thus, the main geophysical observations, including the seismic activity, the deformations and the velocity variations, are all consistent with an open magmatic system in which no marked pressurization occurred at shallow level before the July 2015 eruptions.

At Soufriere Hills volcano, Montserrat, in 1996–1998, most large dome collapses occurred during periods of high extrusion rate, intense hybrid seismic events activity and cyclic deformations (Calder et al., 2002). However some of the large collapses occurred while no magma was extruding and they were not preceded by seismic activity. These events were interpreted as structural failures of steep crater walls (Calder et al., 2002). The two dome collapses of July 2015 at Volcán de Colima appear to be intermediate cases, as they occurred during an episode of high rate of extrusion and were preceded by a decrease of seismic activity. These dome collapses probably resulted from a mechanical instability of the crater walls triggered by the strong magma flow through an open conduit. This type of eruptive event, which is not preceded by usually observed precursors, is thus difficult to forecast with classical monitoring methods. However, it would be important to identify this kind of situation in the future in order to manage better the corresponding hazards. The integration of new observations and analysis methods to the monitoring system may also help detecting forthcoming eruptions. For example, pixel offsets tracking methods applied to optical images acquired at small distance from the dome might bring useful information regarding the dome growth rate (Salzer et al., 2016) while high resolution SAR images can provide information on the dome deformation in quiescent periods (Salzer et al., 2017). Indeed, a hypothetical observation of both summit deformation and velocity variations, with amplitude larger than the usual fluctuations (i.e., >0.2%) and not related with strong tectonic earthquakes, could indicate a possible impending eruption and should be taken into account by the warning system of Volcán de Colima.

Author Contributions

PL coordinated this study, carried out the calculation of the seismic noise correlation functions and of the velocity variations and wrote the main parts of the manuscript. AC processed the Sentinel-1 data under the supervision of VP. Together they interpreted Insar results and wrote the sections about deformations. RA-M provided the seismic data used in this study and information on the eruptive activity of the volcano. All the authors revised the whole text and figures.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the European Space Agency for providing Sentinel-1 data. This study was supported by CNES through the TOSCA project MerapiSAR. We are grateful to Gabriel Reyes-Dávila, Carlos Ariel Ramírez, Alejandro Velasquez Martinez, and Miguel González Amezcua, from RESCO, for their support in keeping the network working. The monitoring network of the Volcán de Colima is partially maintained by the CONACYT project Atención a Problemas Nacionales 2015-916. We thank the reviewers and editors for their suggestions that helped improving the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2018.00093/full#supplementary-material

References