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

The geometry of the volcanic conduit is a main parameter controlling the dynamics and the style of volcanic eruptions and their precursors, but also one of the main unknowns. Pre-eruptive signals that originate in the upper conduit region include seismicity and deformation of different types and scales. However, the locality of the source of these signals and thus the conduit geometry often remain unconstrained at steep sloped and explosive volcanoes due to the sparse instrumental coverage in the summit region and difficult access. Here we infer the shallow conduit system geometry of Volcan de Colima, Mexico, based on ground displacements detected in high resolution satellite radar data up to seven hours prior to an explosion in January 2013. We use Boundary Element Method modeling to reproduce the data synthetically and constrain the parameters of the deformation source, in combination with an analysis of photographs of the summit. We favour a two-source model, indicative of distinct regions of pressurization at very shallow levels. The location of the upper pressurization source coincides with that of post-explosive extrusion; we therefore attribute the displacements to transient (elastic) pre-explosive pressurization of the conduit system. Our results highlight the geometrical complexity of shallow conduit systems at explosive volcanoes and its effect on the distribution of pre-eruptive deformation signals. An apparent absence of such signals at many explosive volcanoes may relate to its small temporal and spatial extent, partly controlled by upper conduit structures. Modern satellite radar instruments allow observations at high spatial and temporal resolution that may be the key for detecting and improving our understanding of the generation of precursors at explosive volcanoes.


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 preeruptive 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 domebuilding 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).
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 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. 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.

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 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.
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).

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.
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.

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 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.

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. 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.
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 twosource 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 preexplosive 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

Frontiers in Earth Science | Volcanology
June 2014 | Volume 2 | Article 12 | 6 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.

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.
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,

www.frontiersin.org
June 2014 | Volume 2 | Article 12 | 7 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 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. 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 domebuilding 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.
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 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. 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 offcentered 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.