Satellite radar data reveal short-term pre-explosive displacements and a complex conduit system at Volcán de Colima, Mexico

Introduction

The conduit systems feeding dome-building eruptions can be highly dynamic, with irregular changes in the location and style of magma extrusion (Bernstein et al., 2013). At some volcanoes, the timing and dimension of extrusive activity may be forecasted based on precursory signals such as earthquakes, tremors and long periodic seismicity (Kilburn and Voight, 1998; Neuberg, 2000; Bean et al., 2014), as well as ground displacements on different scales (Surono et al., 2012; Di Traglia et al., 2013) that can be readily identified and have been successfully used for timely evacuation. A large variety of conduit processes may generate precursory signals, including gas accumulation beneath a volcanic dome (Johnson et al., 2008), episodic slip on the conduit margins (Anderson et al., 2010), changes in a dynamic plumbing system (Kahl et al., 2011) or the accumulation of magmatic material at shallow levels (Ratdomopurbo et al., 2013).

However, identifying which periods of unrest will lead to an eruption still remains a challenge (Sparks et al., 2012). Improving our understanding of the processes generating short-term pre-eruptive geophysical signals and ultimately avoiding false alarms require close observations of the dome feeding system and positioning of instrumentation in proximity of the vent. Such observations are rare at explosive volcanoes due to the vulnerability of near summit stations. Therefore, detailed observations of processes in a volcano feeding system are mostly obtained at non-explosive volcanoes, or where complete seismic and geodetic monitoring was achieved. A constrained path the magma propagates through could only be described in a few of those cases, such as Stromboli (Chouet et al., 2008; Di Traglia et al., 2013), Mount Etna (Kahl et al., 2011), Piton de la Fournaise (Massin et al., 2011) or Hawaii (Chouet and Dawson, 2011). At explosive dome-building volcanoes, the location of the conduit is more difficult to constrain. In some cases, we can observe summit deformation caused by the pre-eruptive intrusion of magma to very shallow levels of the volcano (Dzurisin et al., 2008; Saepuloh et al., 2013). Since the newly injected material remains there, this deformation is static and lacks co-explosive subsidence (Ratdomopurbo et al., 2013). However, short term pre-eruptive deformation originating from conduit processes of smaller magnitude are rarely detected with geodetic instrumentation, yet may be related to a considerable proportion of pre-eruptive seismic signals.

We present a model of the uppermost conduit system of Volcán de Colima (Figure 1), a dome building volcano located in Western Mexico, based on geodetic data acquired shortly prior to an explosion. The eruptive history indicates that Volcán de Colima may experience Plinian eruptions with VEI (volcanic explosivity index) >4 approximately every 100 years, the last of which occurred in 1913 (Bretón et al., 2002). Between 2007 and June 2011, the activity at Volcán de Colima was characterized by periods of dome extrusion, as well as pyroclastic flows and explosions. After the last explosive event in June 2011, Colima entered a period of quiescence, until a new period of eruptive activity was initiated by a vulcanian explosion on January 6th, 2013, preceded only by 3 days of increased seismicity (Figure 2). Smaller explosions occurred on January 11th and 13th and were preceded by smaller increases in seismicity. The explosions were followed by relatively low seismicity rates. The pre-explosive seismic activity was observed at the nearest seismic station (EZV4), located at 1.5 km from the summit (Figures 1, 2). A similar type of seismicity has previously been observed at Volcán de Colima (Arámbula-Mendoza et al., 2011) and is possibly related to fluid movement and rock fracturing in the conduit (Bean et al., 2014).

FIGURE 1

Figure 1. Location map and overflight photograph of Volcán de Colima. The steep-sloped volcanic cone is approximately 4–6 km wide, and completely covered by both the ascending and descending TerraSAR-X satellite data (blue and red boxes in map view). The yellow star marks the seismic station EZV4 described in this paper.

FIGURE 2

Figure 2. Seismicity increase in January 2013 and TerraSAR-X interferograms. (A) Time plot with cumulative number of seismic events recorded at EZV4 (blue), timing of TerraSAR-X satellite acquisitions (red arrows) and timing of significant explosions (yellow stars). The interferograms cover the quiescent period, the pre-explosive period and the co-explosive period. (B,C) The re-wrapped pre-explosive interferograms (December 26th–January 6th) show strong deformation affecting the summit area and are used in our modeling study. The descending data (B) were acquired on track number 128, the ascending data (C) on track 121. Each color cycle corresponds to 1.55 cm of line-of sight (LOS) displacement. The arrows indicate the horizontal projection of the LOS and the radar azimuth direction. The seismic station EZV4 lies outside the deforming region.

We use a satellite radar dataset at unprecedented resolution to measure the surface displacements prior to the first explosion at thousands of points on the volcano summit and crater region. The timing of the satellite imagery precedes the 6th January 2013 explosion by hours only (Figure 2). We exploit the radar data interferometrically (InSAR), a method now used routinely to analyze large scale deformation at many volcanoes (Fournier et al., 2010; Ebmeier et al., 2013), and allowing the identification of the depths and geometries of volcanic magma chambers (Dzurisin, 2007). Recent studies have also revealed the importance of both timeliness and spatial resolution when identifying processes associated with basaltic eruptions (Bagnardi et al., 2013; Richter et al., 2013). However, this is the first time that transient plumbing of a dome feeding system was detected as a pre-explosive inflation episode using InSAR. This allows the modeling of the conduit associated with an explosion, linking the deformation to the pattern of pre-explosive seismicity and drawing conclusions on the processes generating these signals.

Methods

Data Processing and Generation of Displacement Maps

Interferometric Synthetic Aperture Radar (InSAR) is a method used to measure ground displacements in the line of sight (LOS) of the satellite that occurred between two SAR acquisitions. We consider SAR images on ascending and descending tracks acquired every 11 days by the German radar satellite TerraSAR-X (TSX) in spotlight mode. This type of acquisition yields a spatial resolution of over 2 m, which is unprecedented at Volcán de Colima and most other volcanoes worldwide. We focus our study on images taken on December 15th, December 26th, January 6th, the day of the first explosion, and January 17th, following another explosion (Figure 2). The acquisitions on January 6th fall into the period of increasing seismicity and precede the explosion by 19 and 7 h on the ascending and descending track, respectively. This imagery therefore allows us to analyze the deformation during the quiescent period during which no seismicity was detected, during the pre-explosive seismicity increase, as well as the time frame covering the explosions. Our data extends further back into the quiescent period in 2012, for simplification we only show the last interferogram.

In order to quantify the surface displacements, we used the processing chain implemented in the ROI_PAC software (Rosen et al., 2004). The Digital Elevation Model (DEM) used for correction of topographic effects was kindly provided by the Universidad Nacional Autonoma de Mexico, based on airborne LIDAR data acquired in 2005, with a horizontal resolution of 5 m. We superimposed a photogrammetric DEM of the volcano summit from December 2011 on the LIDAR data, in order for our DEM to contain the dome built in the years 2007 to 2011 (James and Varley, 2012).

Source Modeling

We use the Boundary Element Method to simulate the observed surface displacement and herewith to quantify parameters such as the pressurization and the geometry of the source. This modeling method is based on the analytical solution of a triangular dislocation (TD) element in an elastic full-space (Yoffe, 1960). The TD elements allow us to discretize any curved surface, such as topography and the geometry of complex sources without any discontinuity (Kuriyama and Mizuta, 1993; Maerten et al., 2005). In our model, we simulate the pressurization of the source by defining a traction boundary condition at the source TD elements. The topography on the other hand is simulated as a traction-free surface, i.e., the traction at the centroids of all surface TD elements is zero (Walter et al., 2005). By computation of the stress influence coefficient matrix of our model and considering the imposed boundary conditions we are able to solve for the unknown slip and opening components at the TDs (Jeyakumaran et al., 1992), which allows the calculation of the displacements at any point in the model. We have chosen this modeling method because it considers topographic effects (Cayol and Cornet, 1998) and also the interaction of multiple sources of pressure.

Through a combination of forward modeling and optimizations we minimize the misfit between the modeled displacements and the data. Forward models first aim at investigating the effect of the source type on the modeled surface displacements and confine the solution space. We use these search parameter ranges to find the solution with minimum misfit by performing a non-linear optimization using a genetic algorithm (Haupt and Haupt, 2004). In our optimization we use a population size of 40–50, mutation rate of 0.2, selection rate of 0.5 and a maximum iteration of 50, which lead to stable convergence of the results. We ran up to 5 optimizations for each model setup, each run taking approximately 20 h on a standard PC.

To prepare the data set for the model optimization we reduce the number of data points in the InSAR displacement maps from several million pixels to 6000 data points using a Quadtree subsampling approach (Jónsson et al., 2002). Quadtree subsampling is an irregular subsampling approach assigning higher point densities in areas of larger displacement gradients. In our subsampled data, the summit region of strong deformation is represented by a higher number of points than the far-field; and by weighting the points equally, we assign a more relevant role to the summit data during the optimization. This irregular subsampling is essential, since the deforming region only covers a fraction of the full dataset. It also allows us to limit the impact of far-field data, which are not relevant for our model, on the optimization results.

The topographic mesh consisted of 3250 TD elements and was constructed based on the same DEM as used for the InSAR data processing (Persson and Strang, 2004). We used a mesh spanning a radial distance of 1 km to the summit crater, and reaching a resolution of 5 m in the center, corresponding to the summit of the volcano and region of highest displacements.

We chose a Young's modulus of 1 GPa and a Poissons Ratio of 0.25 for our model. The low Young's modulus reflects the fragmented, uncompacted lava and ash matrix, which commonly composes the summit structure at dome-building volcanoes. Similar values have been constrained at Montserrat (Voight et al., 1999) and at Merapi Volcano (Beauducel et al., 2000).

Results

InSAR

The pre-explosive ascending and descending interferograms (December 26th–January 6th) both reveal significant deformation that is focused around the upper edifice of Volcán de Colima (Figure 2). Coherent interferometric phase is found at pixels almost throughout the volcano edifice, except for some parts of the dome summit and of the flank areas. Loss of coherence is possibly due to the high displacement gradients and partial snow cover near the summit, and vegetation on the lower flanks. The seismic station EZV4 lies outside the strongly deforming area.

The unwrapped interferograms covering the quiescent period (December 15th to December 26th) show that no major deformation took place during the 11 to 22 days prior to the January 6th explosion, which is in strong contrast with the pre-explosive interferograms described above (Figure 3). This change in detectable surface displacements is confirmed in both ascending and descending tracks. The descending interferogram shows a butterfly pattern, with oppositely directed displacements approaching 6 cm in the satellites line of sight (LOS) on either flank of the volcano. The displacements in the ascending interferogram are, however, more asymmetric, with strong movements also approaching 6 cm in LOS on the western flank, yet little LOS displacements are observed on the eastern flank. In both datasets, the deformation is confined to the summit area. The displacement patterns on the ascending and descending interferograms are found to be consistent, suggesting a complex source of deformation likely to be related to the explosion that removed a part of the dome only few hours later.

FIGURE 3

Figure 3. InSAR data shows deformation shortly prior to the January 2013 eruption. Top panels: no significant displacements are observed at the summit during the quiescent period (15th December–26th December). Bottom panels: pre-explosive deformation (26th December–6th January, 11—1 days prior to the explosion). The arrows indicate the satellite LOS and azimuth direction.

The co-explosive interferograms (January 6th–January 17th) are strongly decorrelated in contrast to the pre-explosive interferograms due to the deposition of ash and other material. We can however still identify fringes locally in the summit area in both satellite look directions (Figure 4). The spacing of the fringes in the co-explosive interferogram is found to be almost identical with that of the “pre-explosive” fringes, but opposite in sign. Also, the asymmetry of the pre-explosive deformation is visible in the co-explosive interferograms: the descending interferogram exhibits two lobes of displacements with opposite signs, while only a single displacement lobe is visible in the ascending interferogram. The co-explosive displacements therefore appear to be the evidence for the reversing of the process producing the pre-explosive displacements. The deformation is hence considered elastic, and no significant amount of magma has intruded and arrested within the volcano summit.

FIGURE 4

Figure 4. InSAR data prior to and covering the January explosion. The co-explosive deformation in LOS (in color) is shown overlying the pre-explosive inflation fringes (grayscale) in the descending and ascending interferograms (A and B, respectively). The shaded relief of the summit dome can be seen in the center, not covered by data.

Source Modeling

Our modeling efforts are focused on explaining the pre-explosive displacements. We first attempted to model these by using single pressurized sources of various shapes, i.e., volumetric bodies and planar (sill- or dike-like) sources (Dzurisin, 2007). The better fit to the data was achieved using pressurized cylindrical and ellipsoidal sources, the inversion results of which are shown in Figure 5. We assess the ability of the models to explain the observations by investigating the patterns of the simulated surface displacements, residuals as well as plots showing the data and modeled displacements along an E-W profile crossing the summit.

FIGURE 5

Figure 5. Boundary element modeling results. Three different model setups are shown: a cylindrical source (A), an ellipsoidal source (B) and a combination of two ellipsoids (C). For each, we present the best-fitting source geometry and its parameters. In addition, we provide maps showing the modeled displacements and residuals in both look directions for each setup. Under each map are profiles (P-P') of the data (red) and the modeled displacements (blue). These profiles comprise a projection of all values within a 40 m-wide strip centered on the line of profile. The marked regions in the profiles (a–i) are referred to in the main text.

The overall results for a single ellipsoid or cylinder are quite similar in the modeled displacement range, patterns and fit to the data, however the source parameters are slightly different in terms of size and depth. The distribution of the residuals suggests that the cylindrical source provides a slightly better fit to the data; however, neither source can explain the characteristic asymmetric patterns observed in the displacements. The E-W profiles in particular reveal the inconsistent shapes of the two curves: while the observations show an exponential pattern on both the eastern and western segments of the descending data, the modeled displacement curves are bell-shaped and often divergent when approaching the summit (Figure 5, profiles P-P', labels a-f). The lack of displacements in the ascending data up to the immediate summit proximity could not be reproduced by any of the single-source models (Figure 5, profiles P-P', labels c, f), without significantly compromising the fit to the data elsewhere. The displacement pattern could therefore not universally be explained by any of the single sources investigated. We attributed this to the real source geometry being more complex in shape than our simple single-source models. Furthermore, the explosion left an almost circular shallow crater in the summit dome (see Supplementary Figure 2), which indicates pressurization at levels which are more shallow than those suggested by our single-source optimization results. Therefore, we approximate a more complex source geometry by using a combination of two sources.

In order to stabilize our results in the two-source model optimization, we narrow down the solution space by first modeling the data on the volcano flanks (excluding a radius of 350 m about the summit), using an ellipsoidal source, through an inversion (Supplementary Figure 1). The residuals for this source were in the mm-range everywhere outside the 350 m radius, however, as expected, very high within. We kept this source model and added a second, shallow triaxial ellipsoid to simulate the displacements within the near field of the summit. In the following optimization we constrain the search ranges for the deep source parameters based on the result of the previous inversion. The limits for the geometrical parameters for the shallow ellipsoid were given by the summit topography, which it may not intersect. We also found that the E-W profiles of the displacements are very sensitive to the position, E-W extension and pressurization of the shallow source; hence these parameters appear well constrained in our setup. Including the data within the 350 m radius again, we initialize the inversion with two sources, the result of which is shown in Figure 5C.

The best-fitting two-source model is made up by a deep, vertically elongated source located 150–300 m depth below the summit, and a shallow, horizontally elongated source, located just a few tens of meters below the summit and inside the dome. The deep source has a near circular horizontal aspect ratio, and is located within the volcanic cone at a centered position with respect to the main edifice. The distinct shallow source, on the other hand, is significantly offset to the West, extended in North-South direction and almost perfectly centered beneath the active dome. Although location and geometry of these two sources are distinct, a similar amount of pressurization (22 and 25 MPa for the shallow and deep sources, respectively) was estimated in our optimization. The similar degree of pressurization of the two sources suggests a physical connection between them.

The residual displacements are reduced overall in the two-source model, and distributed more equally over the ascending and descending datasets. However, from the RMS values alone, the improvement of the data fit is not significant enough to justify the two-source model over the single-source models. We still favor the two-source model since it can better reproduce the complexities in the observed displacement gradients, such as the increase toward the summit and the existence of both long and short wavelength signals. The modeled displacement patterns of the two-source model follow the observations more closely, in particular the lack of displacements in the eastern ascending data (Figure 5i), as well as the exponential trend in the other segments (Figures 5g,h). Also, short wavelength changes in the observed displacements can now be reproduced (Figure 5i).

Discussion

Our results indicate that the summit area of Volcán de Colima was first at rest, and then unexpectedly underwent strong deformation just during a few days prior to the January 6th explosion. This we attribute to shallow sources of pressurization located within the shallow volcanic conduit. Source complexities are found concerning the location and orientation of the two sources, determined from a lateral shift from the centralized deeper source to an off-centered shallow source, and from a change of a vertically extended radial symmetric geometry at depth to a horizontally extended, fracture-like geometry near the surface. The pre-explosive pressure increase lies in the order of 22–25 MPa, values that are realistic for plumbing systems of this setting (Voight et al., 2010; Lavallée et al., 2012). During the period covering the explosion, the summit area subsided again, and did so by an amount and in a spatial pattern similar to the pre-explosive inflation. Our study therefore shows, for the first time, on a very small spatial and temporal scale, the transient pre-explosive pressure increase in the uppermost conduit system of Volcán de Colima.

Observations of precursory deformation with InSAR are relatively rare at explosive volcanoes, which has been attributed to the characteristics of the magma reservoir (Chaussard and Amelung, 2012; Ebmeier et al., 2013), a lack of pressurization (Chaussard et al., 2013), the steep topography (Pinel et al., 2011) or temporal aliasing (Fournier et al., 2010). The short temporal baseline and high resolution of our satellite imagery, in combination with a high quality DEM, reduces errors typically associated when using standard resolution InSAR in this type of terrain (Pinel et al., 2011). In the case of Colima, the timing of the satellite acquisitions was critical for detecting the deformation.

Pre-explosive pressurization previously observed at the conduit of Montserrat (Voight et al., 2010) may resemble the situation observed at Volcán de Colima, as in both cases pressurization has occurred within the edifice, very close to the actual eruption center, but the case presented here stands out for two reasons: firstly, because of the short period during which the pressurization developed, and secondly, because the high spatial resolution of the displacement data in the summit area allowed resolving the summit deformation and modeling the underlying sources. The first point implies that the conduit system feeding dome explosions at Volcán de Colima was temporally plugged, leading to the January 2013 volcano explosion; the second point allows us to draw conclusions on the geometry of the conduit system and its effect on the generated displacements.

Before discussing the broader implications of these findings, we address the limitations of our study.

Limitations

Data limitations

The number of seismic events at Volcán de Colima began increasing exponentially 3 days prior to the eruption. Similarly, the rate of pre-eruptive deformation can increase toward the explosion (Sigmundsson et al., 2010), therefore, the 12 h time difference lying between the ascending and the descending slave images may have caused a certain amount of additional displacements in the descending interferogram. We cannot account for such a possible mismatch, as any model based solely on one data set would be poorly defined. However, due to the similar spatial extent and maximum amounts of displacements in both look directions we believe that the effect may be neglected in this study. Other instrumentation, such as GPS, tiltmeters or EDM, could have provided a higher time resolution and hence allowed the consideration of the 12-h gap between our TerraSAR-X datasets. However, given the complexity of the deformation and the positioning of other instrumentation at greater distance from the summit, any such point-wise measurements would not have allowed constraining the two-source model as we have done here.

Model limitations

The main limitations of our modeling are the assumptions of elastic deformation and a homogenous medium. Material heterogeneities, such vertical and lateral changes in the elastic parameters, are known to strongly affect the source evolution, e.g., by arresting the ascent of magmatic material when it enters a stiffer medium (Maccaferri et al., 2011) but also the displacement pattern at the surface (Manconi et al., 2007). At Volcán de Colima, changes in the material properties may be found for instance at the contact between the dome and the pre-existing crater, but are also produced by the thermal, structural and compositional changes relating to the conduit and the volcanic cone itself. This may lead to distinct zones of pressurization, such as in the dome interior. Whether the models' complexity is an artifact of the homogeneous material assumption, or whether the distinct two zones of pressurization result from mechanical contrasts in the dome and edifice itself, cannot be assessed using the available surface deformation data alone, and necessitates additional information on the mechanical heterogeneities.

The near summit region may also experience non-elastic deformation. The very near field of the summit dome in particular (<50 m away from the dome) may have responded in a brittle fashion, which would explain the high displacement residuals in this region. However, brittle deformation did not have a strong influence on the pre-explosive displacements at Volcán de Colima. This is evidenced by, firstly, the continuity of the signal in the InSAR data suggesting no faulting near the surface. Secondly, the explosion itself reflects the failure of the dome material occurring after the observed deformation. Thirdly, the co-explosive interferograms (Figure 4) show displacements in opposite sign and similar spatial wavelength to the pre-explosive deformation up to very close to the summit, the process leading to the displacements therefore needs to be reversible. We conclude that our elastic model assumption is well justified.

The pressure changes suggested by our models are on the order of 20 MPa, which might appear rather high if compared to the rock mass strength values of below 10 MPa as determined at similar settings (Sparks, 1997; Lavallée et al., 2012). However, the required pressurization in geodetic models is dependent on the values used to describe the material properties in the model, which, at Volcán de Colima, are difficult to constrain. Also, in this particular case, the presence of a solidified plug and a cooled dome may lead to strength values which are considerably larger than 20 MPa (Schultz, 1995). Additionally, shear along the conduit walls may also contribute to the deformation during dome-building eruptions (Beauducel et al., 2000; Green et al., 2006). Mechanical testing of the material properties is therefore a prerequisite to further improve models and understand the nature of precursory deformation.

Due to the long computation time, our optimization does not allow to constrain true confidence bounds of our model results based on multiple optimizations. However, we investigated the development of the individual model parameters during the inversion (see Supplementary Materials) and analyzed the range in which the parameters converge. While this does not provide model uncertainties, it allows estimating the stability of the model results and how well the parameters are constrained. We conclude that the main features of the source model on which we base our interpretations are robust.

Implications

In order to corroborate our modeling results and put them into relation with the volcano's activity, we analyze a set of photographs taken during an airplane overflight on January 10th, only 4 days after the explosion (Figure 6). We geocoded the photograph showing the western summit and the January 6th explosion crater using distinct features visible both in the photograph and in the high-resolution DEM of the volcano summit. The extrusion of the uppermost part of the conduit material can be identified within the explosion crater. This location almost perfectly agrees with the position of the shallow pressurization source, which lies just beneath the dome. The figure also provides the projected location and approximate extension of the deep spheroid, positioned central to the edifice and offset with respect to the dome, which has important implications for the plumbing system and magma pathway.

FIGURE 6

Figure 6. Geocoded photograph of the summit taken on January 10th during an airplane overflight. The circular markers show the approximate projected locations of the upper source (in orange) and the lower source (in green), as well as their approximate extents. The white arrow points to the location of the new extrusion, located at the center of the crater left by the January 6th explosion. This extrusion was removed by the following explosion on January 11th.

The good spatial agreement found between the locations of the shallow source and renewed dome extrusion as seen in aerial photographs suggests that the region was first blocked and pressurized, therefore detectable as a precursory signal by InSAR data, and began serving as a feeding system within less than 4 days following the explosion.

Pressurized regions along a volcanic conduit system may develop at geometric bends or changes in the conduit (Hautmann et al., 2009). It has been proposed that the upper conduit at Volcán de Colima may curve away from the center due to a strong plug of degassed material located beneath the summit crater (Lavallée et al., 2012). In this view, the feeder is centralized under the edifice, but bending westward toward the active dome close to the surface. Our InSAR data and model results agree with this conceptual model derived from independent data by Lavallée et al. (2012), and adds additional constraints on the localization of this postulated curvature. We postulate that our deeper source, centered at a depth of around 220 m below the summit, is the result of a geometrical change produced by the plug. The shallow source in turn reflects the transition into the summit dome.

Possibly, a larger scale structural trend may also play a role in defining the geometry of the sources. Volcán de Colima is located in an extensional regional stress field directed N-S, and in an E-W trending volcanic graben, leading to gravitational spreading of the volcanic complex toward the south (Norini et al., 2010). It has been suggested that these trends may also be reflected in dome structures (James and Varley, 2012). However, the axial symmetry of the deeper source suggests that the shape of the edifice plays a stronger role than regional trends in controlling the geometry of the shallow conduit system at Volcán de Colima. The N-S extension of the shallow source is more likely to be controlled by the shape and local stress field of the dome, overflowing to the West (Walter et al., 2013), which in turn is a consequence of the off-centered position of the vent. The N-S structures visible in the 2007–2011 Colima dome are therefore likely to be dominantly controlled by the local rather than the regional stress field.

We therefore describe the current system feeding dome-building eruptions at Colima by a conduit-like source within the main edifice, overlain (with a lateral offset) by a horizontally extended source. Two shallow but distinguished zones of pressure buildup have also been localized at other volcanoes, one located within or just beneath a dome with a visco-elastic dome cap (Johnson et al., 2008), and one located few hundred meters below, within the edifice (Holland et al., 2011). Our results suggest that a combination of processes was taking place at Colima prior to the January 6th explosion. A physical connection, such as a fracture network, could have developed between the shallow and the deeper source, accompanied by pre-explosive seismicity, until the pressure overcame the strength of the dome cap. The explosions released the pressure built up in the conduit, leading to co-explosive summit subsidence (Figure 7). The prolonged low rates of seismic events following the explosions suggest that the pressure release was not followed by immediate re-pressurization. Possibly, the increase in seismicity visible prior to the January 14 explosion suggests another pressurization cycle, smaller in magnitude, took place then. However, due to temporal aliasing, this deformation would not be visible in our interferograms.

FIGURE 7

Figure 7. Conceptual model. Transient pressurization at distinct locations of the conduit at Volcán de Colima. The sources are located in the uppermost region of the edifice. The depth of the lower source reflects the bifurcation of the conduit and presence of a plug, the upper source is located at the transition between the crater floor and the volcanic dome.

Our study differs to other investigations of pre-eruptive inflation (Fournier et al., 2010; Lu et al., 2010; Chaussard and Amelung, 2012) in that there are no significant amounts of new deeply sourced material detected in our InSAR data, but the bulk of the deformation was caused at very shallow levels. We do not attribute the deformation to the inflation of a deep magma reservoir, but to transient pressurization of the dome feeding system, possibly only lasting a few days. The ultimate trigger for this pressurization however still remains unclear. Camera observations show significant snowfall beginning around 4 days prior to the explosion, associated with increased fumarolic activity, suggesting that precipitation may have percolated into the edifice leading to a pore pressure buildup at depth.

The time span over which displacements were observed as well as the proximity of the deformation sources to the summit underline the need for a configuration of the monitoring instruments that allows the detection of such short term pre-explosive changes. The spatial pattern of the pre-explosive surface displacements at Volcán de Colima was mainly controlled by the internal structure of the summit, with the size and position of the conduit, the presence of a plug and the shape of the edifice all playing a role. Therefore, these factors may also affect the distribution of pre-explosive signals generated by conduit pressurization at other volcanoes, i.e., the shallow pressurization of a geometrically simple conduit system may lead to signals smaller in magnitude and spatial amplitude, and hence even more difficult to detect. This type of precursors may be detected and early warning measures implemented only when monitoring systems operate at a sufficiently high spatial and temporal resolution.

Conclusions

The January 6th 2013 Colima explosion is a rare example where satellite radar interferometry allowed the quantification of activity changes only few hours prior to an explosive volcanic eruption. The geometry of the reactivated conduit and the location of renewed activity can be constrained using high-resolution InSAR data. Ground displacement on the volcano flanks are best explained by a cigar-shaped source centered within the main edifice and subject to pressurization. The displacements at the volcano summit, in turn, can only be explained by an additional source located inside the dome.

The two-source model is supported by independent observations, such as the characteristics and location of the explosion crater as evidenced by photographs of the summit. It is also in agreement with the conceptual model presented by Lavallée et al. (2012), based on geological data. Our results support the existence of a plug beneath the crater floor at Volcán de Colima, implying the presence of a curved ascent path leading to off-centered extrusion sites at the surface. We associate the lower source with the depth of the plug, which we now speculate to reach a depth of around 200–250 m. The geometry of the upper source appears to be controlled by the local stress field of the dome, while the geometry of the lower source may be an approximation of more complex flow patterns, or imposed by the presence of the plug. These conduit complexities are directly reflected in the distribution of the observed displacements, suggesting that the spatial dimension of pre-explosive signals at other volcanoes may also depend on the properties of the conduit system.

The spatiotemporal characteristics of the signals shown here suggest that the occurrence of explosive eruptions lacking deformation precursors may be due to the poor spatial and/or temporal resolution of geodetic monitoring data in the summit area. In our study, the timing of the satellite acquisitions was crucial for detecting the deformation. The transient nature of the process generating them suggests that temporal aliasing, rather than a lack of pressurization, limits our ability to detect pre-explosive deformation by InSAR.

We therefore believe future satellite missions with a high revisiting frequency, as well as the installation of near-summit monitoring stations, will reveal more details about the upper conduit system at other volcanoes. The combined analysis of short-term pre-explosive seismic and geodetic data of high resolution will improve our understanding of the generation of these signals their relevance for eruption forecasts.

Author Contributions

Jacqueline T. Salzer processed the radar images, analyzed the modeling results, performed the comparison with the seismic and photographic data, prepared the figures and wrote the manuscript. Mehdi Nikkhoo developed the Boundary Element Modeling code and modeled the displacement data. Henriette Sudhaus contributed to the analysis and modeling of the InSAR data and interpretation of the results. Thomas R. Walter supervised the study, adding relevant scientific content throughout its development. Raúl Arámbula and Gabriel Reyes-Dávila provided and analyzed the seismic records. Mauricio Bretón provided camera data as well as background information on the January 6th explosion.

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

This work is supported by the Initiative and Networking Fund of the Helmholtz Association, in the framework of the “Helmholtz Alliance: Remote Sensing and Earth System Dynamics,” GFZ Potsdam and the TOPOMOD training project in the Marie Curie Actions-Initial Training Networks (ITN). The TerraSAR-X images were provided by the German Aerospace Center (DLR proposal ID 1505). We thank Lucia Capra from the Universidad Nacional Autónoma de Mexico and Michael James from Lancaster University for providing the DEMs used in this study. The photographs were provided by Protección Civil Jalisco. The authors wish to thank Jose Manuel Alvarez Nieves for technical assistance and company in the field at Volcán de Colima, and Eoghan Holohan and Michele Pantaleo for their comments and constructive suggestions. Figures 3, 5, and Supplementary Figures 1, 4 were produced with the aid of GMT (Wessel et al., 2013).

Supplementary Material

The Supplementary Material for this article can be found online at: http://www.frontiersin.org/journal/10.3389/feart.2014.00012/abstract

References

Anderson, K., Lisowski, M., and Segall, P. (2010). Cyclic ground tilt associated with the 2004–2008 eruption of Mount St. Helens. J. Geophy. Res. 115, B11201. doi: 10.1029/2009jb007102

CrossRef Full Text

Arámbula-Mendoza, R., Lesage, P., Valdés-González, C., Varley, N. R., Reyes-Dávila, G., and Navarro, C. (2011). Seismic activity that accompanied the effusive and explosive eruptions during the 2004–2005 period at Volcán de Colima, Mexico. J. Volcanol. Geotherm. Res. 205, 30–46. doi: 10.1016/j.jvolgeores.2011.02.009

CrossRef Full Text

Bagnardi, M., Amelung, F., and Poland, M. P. (2013). A new model for the growth of basaltic shields based on deformation of Fernandina volcano, Galápagos Islands. Earth Planet. Sci. Lett. 377–378, 358–366. doi: 10.1016/j.epsl.2013.07.016

CrossRef Full Text

Bean, C. J., De Barros, L., Lokmer, I., Metaxian, J.-P., O' Brien, G., and Murphy, S. (2014). Long-period seismicity in the shallow volcanic edifice formed from slow-rupture earthquakes. Nat. Geosci. 7, 71–75. doi: 10.1038/ngeo2027

CrossRef Full Text

Beauducel, F., Cornet, F. H., Suhanto, E., Duquesnoy, T., and Kasser, M. (2000). Constraints on magma flux from displacements data at Merapi volcano, Java, Indonesia. J. Geophy. Res. 105, 8193–8203. doi: 10.1029/1999jb900368

CrossRef Full Text

Bernstein, M., Pavez, A., Varley, N., and Calder, E. S. (2013). Rhyolite lava dome growth styles at Chaitén Volcano, Chile (2008-2009): interpretation of thermal imagery. Andean Geol. 40, 295–309. doi: 10.5027/andgeoV40n2-a07

CrossRef Full Text

Bretón, M., Ramirez, J. J., and Navarro, C. (2002). Summary of the historical eruptive activity of Volcán De Colima, Mexico 1519–2000. J. Volcanol. Geotherm. Res. 117, 21–46. doi: 10.1016/S0377-0273(02)00233-0

CrossRef Full Text

Cayol, V., and Cornet, F. H. (1998). Effects of topography on the interpretation of the deformation field of prominent volcanoes—Application to Etna. Geophys. Res. Lett. 25, 1979–1982.

Chaussard, E., and Amelung, F. (2012). Precursory inflation of shallow magma reservoirs at west Sunda volcanoes detected by InSAR. Geophys. Res. Lett. 39, L21311. doi: 10.1029/2012gl053817

CrossRef Full Text

Chaussard, E., Amelung, F., and Aoki, Y. (2013). Characterization of open and closed volcanic systems in Indonesia and Mexico using InSAR time series. J. Geophy. Res. 118, 3957–3969. doi: 10.1002/jgrb.50288

CrossRef Full Text

Chouet, B., and Dawson, P. (2011). Shallow conduit system at Kilauea Volcano, Hawaii, revealed by seismic signals associated with degassing bursts. J. Geophy. Res. 116, B12317. doi: 10.1029/2011jb008677

CrossRef Full Text

Chouet, B., Dawson, P., and Martini, M. (2008). Shallow-conduit dynamics at Stromboli Volcano, Italy, imaged from waveform inversions. Geol. Soc. (Lond.) 307, 57–84. doi: 10.1144/sp307.5

CrossRef Full Text

Di Traglia, F., Del Ventisette, C., Rosi, M., Mugnai, F., Intrieri, E., Moretti, S., et al. (2013). Ground-based InSAR reveals conduit pressurization pulses at Stromboli volcano. Terra Nova 25, 192–198. doi: 10.1111/ter.12020

CrossRef Full Text

Dzurisin, D. (2007). Volcano Deformation: Geodetic Monitoring Techniques. Berlin: Springer.

Dzurisin, D., Lisowski, M., Poland, M. P., Sherrod, D. R., and Lahusen, R. G. (2008). Constraints and Conundrums Resulting from Ground-Deformation Measurements Made During the 2004-2005 Dome-Building Eruption of Mount St. Helens. Washington, DC: US Geological Survey Professional Paper.

Ebmeier, S. K., Biggs, J., Mather, T. A., and Amelung, F. (2013). On the lack of InSAR observations of magmatic deformation at Central American volcanoes. J. Geophy. Res. 118, 2571–2585. doi: 10.1002/jgrb.50195

CrossRef Full Text

Fournier, T. J., Pritchard, M. E., and Riddick, S. N. (2010). Duration, magnitude, and frequency of subaerial volcano deformation events: new results from Latin America using InSAR and a global synthesis. Geochem. Geophy. Geosys. 11, Q01003. doi: 10.1029/2009gc002558

CrossRef Full Text

Green, D. N., Neuberg, J., and Cayol, V. (2006). Shear stress along the conduit wall as a plausible source of tilt at Soufrière Hills volcano, Montserrat. Geophys. Res. Lett. 33, L10306. doi: 10.1029/2006gl025890

CrossRef Full Text

Haupt, R. L., and Haupt, S. E. (2004). Practical Genetic Algorithms. New York, NY: John Wiley & Sons.

Hautmann, S., Gottsmann, J., Sparks, R. S. J., Costa, A., Melnik, O., and Voight, B. (2009). Modelling ground deformation caused by oscillating overpressure in a dyke conduit at Soufrière Hills Volcano, Montserrat. Tectonophysics 471, 87–95. doi: 10.1016/j.tecto.2008.10.021

CrossRef Full Text

Holland, A. S. P., Watson, I. M., Phillips, J. C., Caricchi, L., and Dalton, M. P. (2011). Degassing processes during lava dome growth: insights from Santiaguito lava dome, Guatemala. J. Volcanol. Geotherm. Res. 202, 153–166. doi: 10.1016/j.jvolgeores.2011.02.004

CrossRef Full Text

James, M. R., and Varley, N. (2012). Identification of structural controls in an active lava dome with high resolution DEMs: Volcán de Colima, Mexico. Geophys. Res. Lett. 39, L22303. doi: 10.1029/2012gl054245

CrossRef Full Text

Jeyakumaran, M., Rudnicki, J. W., and Keer, L. M. (1992). Modeling slip zones with triangular dislocation elements. Bull. Seismol. Soc. Am. 82, 2153–2169.

Johnson, J. B., Lees, J. M., Gerst, A., Sahagian, D., and Varley, N. (2008). Long-period earthquakes and co-eruptive dome inflation seen with particle image velocimetry. Nature 456, 377–381. doi: 10.1038/nature07429

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jónsson, S., Zebker, H., Segall, P., and Amelung, F. (2002). Fault slip distribution of the 1999 mW 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. Bull. Seismol. Soc. Am. 92, 1377–1389. doi: 10.1785/012167990922

CrossRef Full Text

Kahl, M., Chakraborty, S., Costa, F., and Pompilio, M. (2011). Dynamic plumbing system beneath volcanoes revealed by kinetic modeling, and the connection to monitoring data: An example from Mt. Etna. Earth Planet. Sci. Lett. 308, 11–22. doi: 10.1016/j.epsl.2011.05.008

CrossRef Full Text

Kilburn, C. R. J., and Voight, B. (1998). Slow rock fracture as eruption precursor at Soufriere Hills Volcano, Montserrat. Geophys. Res. Lett. 25, 3665–3668. doi: 10.1029/98gl01609

CrossRef Full Text

Kuriyama, K., and Mizuta, Y. (1993). Three-dimensional elastic analysis by the displacement discontinuity method with boundary division into triangular leaf elements. Int. J. Rock Mech. Mining Sci. Geomechanics Abstracts 30, 111–123. doi: 10.1016/0148-9062(93)90704-H

CrossRef Full Text

Lavallée, Y., Varley, N., Alatorre-Ibargüengoitia, M., Hess, K.-U., Kueppers, U., Mueller, S., et al. (2012). Magmatic architecture of dome-building eruptions at Volcán de Colima, Mexico. Bull. Volcanol. 74, 249–260. doi: 10.1007/s00445-011-0518-4

CrossRef Full Text

Lu, Z., Dzurisin, D., Biggs, J., Wicks, C., and Mcnutt, S. (2010). Ground surface deformation patterns, magma supply, and magma storage at Okmok volcano, Alaska, from InSAR analysis: 1. Intereruption deformation, 1997–2008. J. Geophy. Res. 115, B00B02. doi: 10.1029/2009jb006969

CrossRef Full Text

Maccaferri, F., Bonafede, M., and Rivalta, E. (2011). A quantitative study of the mechanisms governing dike propagation, dike arrest and sill formation. J. Volcanol. Geotherm. Res. 208, 39–50. doi: 10.1016/j.jvolgeores.2011.09.001

CrossRef Full Text

Maerten, F., Resor, P., Pollard, D., and Maerten, L. (2005). Inverting for slip on three-dimensional fault surfaces using angular dislocations. Bull. Seismol. Soc. Am. 95, 1654–1665. doi: 10.1785/0120030181

CrossRef Full Text

Manconi, A., Walter, T. R., and Amelung, F. (2007). Effects of mechanical layering on volcano deformation. Geophy. J. Int. 170, 952–958. doi: 10.1111/j.1365-246X.2007.03449.x

CrossRef Full Text

Massin, F., Ferrazzini, V., Bachèlery, P., Nercessian, A., Duputel, Z., and Staudacher, T. (2011). Structures and evolution of the plumbing system of Piton de la Fournaise volcano inferred from clustering of 2007 eruptive cycle seismicity. J. Volcanol. Geotherm. Res. 202, 96–106. doi: 10.1016/j.jvolgeores.2011.01.008

CrossRef Full Text

Neuberg, J. (2000). Characteristics and causes of shallow seismicity in andesite volcanoes. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 358, 1533–1546. doi: 10.1098/rsta.2000.0602

CrossRef Full Text

Norini, G., Capra, L., Groppelli, G., Agliardi, F., Pola, A., and Cortes, A. (2010). Structural architecture of the Colima Volcanic Complex. J. Geophy. Res. 115, B12209. doi: 10.1029/2010jb007649

CrossRef Full Text

Persson, P., and Strang, G. (2004). A simple mesh generator in MATLAB. SIAM Rev. 46, 329–345. doi: 10.1137/S0036144503429121

CrossRef Full Text

Pinel, V., Hooper, A., De La Cruz-Reyna, S., Reyes-Davila, G., Doin, M. P., and Bascou, P. (2011). The challenging retrieval of the displacement field from InSAR data for andesitic stratovolcanoes: case study of Popocatepetl and Colima Volcano, Mexico. J. Volcanol. Geotherm. Res. 200, 49–61. doi: 10.1016/j.jvolgeores.2010.12.002

CrossRef Full Text

Ratdomopurbo, A., Beauducel, F., Subandriyo, J., Agung Nandaka, I. G. M., Newhall, C. G., Suharna Sayudi, D. S., et al. (2013). Overview of the 2006 eruption of Mt. Merapi. J. Volcanol. Geotherm. Res. 261, 87–97. doi: 10.1016/j.jvolgeores.2013.03.019

CrossRef Full Text

Richter, N., Poland, M. P., and Lundgren, P. R. (2013). TerraSAR-X interferometry reveals small-scale deformation associated with the summit eruption of Kīlauea Volcano, Hawai‘i. Geophys. Res. Lett. 40, 1279–1283. doi: 10.1002/grl.50286

CrossRef Full Text

Rosen, P. A., Hensley, S., Peltzer, G., and Simons, M. (2004). Updated repeat orbit interferometry package released. Eos Trans. Am. Geophys. Union 85, 47. doi: 10.1029/2004eo050004

CrossRef Full Text

Saepuloh, A., Urai, M., Aisyah, N., Sunarta, W. C., Subandriyo, and Jousset, P. (2013). Interpretation of ground surface changes prior to the 2010 large eruption of Merapi volcano using ALOS/PALSAR, ASTER TIR and gas emission data. J. Volcanol. Geotherm. Res. 261, 130–143. doi: 10.1016/j.jvolgeores.2013.05.001

CrossRef Full Text

Schultz, R. A. (1995). Limits on strength and deformation properties of jointed basaltic rock masses. Rock Mech. Rock Eng. 28, 1–15. doi: 10.1007/bf01024770

CrossRef Full Text

Sigmundsson, F., Hreinsdottir, S., Hooper, A., Arnadottir, T., Pedersen, R., Roberts, M. J., et al. (2010). Intrusion triggering of the 2010 Eyjafjallajokull explosive eruption. Nature 468, 426–430. doi: 10.1038/nature09558

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sparks, R., Biggs, J., and Neuberg, J. (2012). Monitoring volcanoes. Science 335, 1310–1311. doi: 10.1126/science.1219485

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sparks, R. S. J. (1997). Causes and consequences of pressurisation in lava dome eruptions. Earth Planet. Sci. Lett. 150, 177–189. doi: 10.1016/S0012-821X(97)00109-X

CrossRef Full Text

Surono, J. P., Pallister, J., Boichu, M., Buongiorno, M. F., Budisantoso, A., Costa, F., et al. (2012). The 2010 explosive eruption of Java's Merapi volcano—A ‘100-year’ event. J. Volcanol. Geotherm. Res. 241–242, 121–135. doi: 10.1016/j.jvolgeores.2012.06.018

CrossRef Full Text

Voight, B., Hidayat, D., Sacks, S., Linde, A., Chardot, L., Clarke, A., et al. (2010). Unique strainmeter observations of Vulcanian explosions, Soufrière Hills Volcano, Montserrat, July 2003. Geophys. Res. Lett. 37, L00E18. doi: 10.1029/2010gl042551

CrossRef Full Text

Voight, B., Sparks, R. S. J., Miller, A. D., Stewart, R. C., Hoblitt, R. P., Clarke, A., et al. (1999). Magma flow instability and cyclic activity at soufriere hills volcano, montserrat, british west indies. Science 283, 1138–1142. doi: 10.1126/science.283.5405.1138

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Walter, T. R., Acocella, V., Neri, M., and Amelung, F. (2005). Feedback processes between magmatic events and flank movement at Mount Etna (Italy) during the 2002–2003 eruption. J. Geophy. Res. 110, B10205. doi: 10.1029/2005jb003688

CrossRef Full Text

Walter, T. R., Legrand, D., Granados, H. D., Reyes, G., and Arámbula, R. (2013). Volcanic eruption monitoring by thermal image correlation: pixel offsets show episodic dome growth of the colima volcano. J. Geophys. Res. 118, 1408–1419. doi: 10.1002/jgrb.50066

CrossRef Full Text

Wessel, P., Smith, W. H. F., Scharroo, R., Luis, J., and Wobbe, F. (2013). Generic mapping tools: improved version released. Eos Trans. Am. Geophy. Union 94, 409–410. doi: 10.1002/2013eo450001

CrossRef Full Text

Yoffe, E. H. (1960). The angular dislocation. Philos. Mag. 5, 161–175.