Monitoring ground movement at Volcán de Colima, Mexico, using Sentinel-1 data and SqueeSAR^®

Volcán de Colima (VdeC) represents the southernmost volcanic edifice of the southward younging Quaternary Colima Volcanic Complex (Luhr and Carmichael 1980), a north–south-trending chain of volcanoes and cinder cones at the southwestern edge of the Trans-Mexican Volcanic Belt (Fig. 1). Volcanism occurs at the Colima Volcanic Complex as a result of the subduction of the Cocos and Rivera plates beneath the North American plate (Garduño-Monroy et al. 1998), and rifting between the Jalisco and Michoacán blocks, forming the north–south-trending Colima graben (Allan 1986). These two processes have led to the generation of alkaline and calc-alkaline volcanism over the last 4.6 myr. The Colima Volcanic Complex comprises the active VdeC edifice itself and the now extinct Nevado de Colima and Volcán Cántaro edifices (c. 6 and 20 km north of VdeC respectively).

Major eruptive activity at VdeC has been observed to occur in cycles of c. 100 years, followed by lava dome growth and destruction, lava flows and moderate Vulcanian explosions, before culminating in a major sub-Plinian to Plinian (Volcanic Explosivity Index, VEI > 4) event terminating the eruptive cycle (Luhr 2002). Historical evidence indicates that VdeC has presented some 40 eruptive periods since written records began in 1576 (Bretón-Gonzalez et al. 2002), with at least five sub-Plinian eruptions and the last in 1913 (VEI 5; Luhr and Carmichael 1980, 1990; Robin et al. 1991; Luhr 2002).

Following the sub-Plinian eruption of 1913, VdeC underwent a c. 50 year period of intermittent dome growth within the crater. Flank effusive activity resumed in 1961 (Luhr and Carmichael 1980), and several eruptive episodes emplaced lava flows and pyroclastic deposits in 1975–76, 1981–82, 1991, 1998–99, 2001–2003, 2004–2005, 2007–11 and 2013–17 (Fig. 2; Varley 2019). The most significant recent eruption, in terms of volume, occurred in 2015, as a long-runout pyroclastic flow of c. 10 km length. This was the largest event since the 1913 cycle-ending sub-Plinian eruption and was not preceded by any detectable precursors in terms of edifice inflation or seismic velocity variations (Lesage et al. 2018). The last explosion of the most recent phase occurred on 3 February 2017, and up to June 2022 there have been only low levels of seismic activity and low gas emissions.

Several recent lava flows have been emplaced, mainly on the western, southwestern and southern flanks of VdeC (erupted in 2013, 2014, 2015 and 2016), as shown in Figure 2. Of the southern flows, the 2015 lava flow is partly overlain by a less voluminous, darker (as it appears in optical imagery; see Fig. 1c) lava flow that erupted between September and October 2016.

The VdeC lavas are almost entirely andesitic in composition, and the more recent flows in the field are present as 20–30 m thick flows, c. 75–120 m in width at the toe extending to c. 200 m in width for the main body of the flow, and of variable length. The structure of the flows typically includes a main flow body topped with levees on both sides, and occasional flow banding or ripples on the metres to tens of metres scale on the top of the flow. The main body is surrounded by a talus slope of loose andesitic boulders, blocks and gravel from the lava flow and from rockfalls. The loose material and the constant risk of rockfalls make access to the edifice for scientific observations and sample collection difficult and dangerous (Fig. 2b and c).

Owing to its highly active nature within a populated area, VdeC has been identified as a significant potential hazard. The area immediately surrounding the volcano is rural, yet c. 300 000 people inhabit the wider region (Saucedo et al. 2010), including the villages of La Yerbabuena (relocated, but still partially occupied) and La Becerrera, 8 and 12 km respectively SW of the volcano, the towns of San Marcos and Quesería 13.5 and 14.5 km to the SE, and the provincial capital, Colima City, a conurbation of c. 150 000 people 30 km to the south (National Census). Hazard maps and assessments have examined the potential hazards for these settlements and concluded that the rural villages near the ravines, such as La Becerrera and San Marcos, are most at risk from pyroclastic density currents (PDCs). Towns and cities in the wider region are likely to be affected by ash and pumice falls with low to moderate risk (Cortés and Téllez 2003; Navarro and Cortés 2003; Capra et al. 2014; Varley et al. 2019). These hazard assessments, along with active monitoring, are used in hazard planning and mitigation by authorities, including evacuations during the eruptive events in 2002 and 2015 (Cuevas-Muñiz and Gavilanes-Ruiz 2017; Rodríguez-Garcia 2019).

The volcano is actively monitored using a variety of techniques, as illustrated in Figure 3. A real-time, short-period and broadband seismometer network (RESCO, Arámbula-Mendoza et al. 2018) is employed by the University of Colima, and volcanic gas is routinely monitored using spectrometer measurements of plumes. Diffuse gas and spring water measurements have identified precursors to eruptions (e.g. Varley and Taran 2003; Varley et al. 2019). Observational monitoring from conventional cameras, stationary real-time webcams and thermal cameras are used to record data on crater and lava flow morphology as well as heat distribution from stationary locations and during overflights (Hutchison et al. 2013). Thermal remote sensing information is obtained from the MIROVA satellite system. Remote sensing studies using DInSAR have often failed to resolve millimetre-scale ground movements for large areas of the volcanic edifice (Walter et al. 2013; Salzer et al. 2014, 2017; Lesage et al. 2018), because of dense vegetation cover on the lower slopes and shifting tephra and rockfalls on the upper edifice, and because the steep-sided shape of the cone can cause radar shadow (depending on the particular SAR viewing geometry). Minor ground movements can provide invaluable early indications of imminent activity and information about the stability of the edifice, and as such the paucity in high-quality ground deformation data at VdeC is a gap that needs to be addressed.

A variety of methods have been employed in the detection and monitoring of ground disturbance at VdeC since the 1990s. Early levelling measurements indicated persistent subsidence of c. 90 mm a⁻¹ at the uppermost volcanic dome (Murray and Wooller 2002). Global positioning system (GPS) measurements have been available from the wider region since the 1990s, although selective availability was in operation before 2000, and therefore measurements of ground motion were limited in accuracy. Murray and Wooller's measurements indicated a total subsidence of c. 23 mm a⁻¹ over the c. 30 years period between 1982 and 2002, caused by gravitational spreading, cooling and compaction. Other studies have shown small-scale continual subsidence, hypothesized as being caused by pressure changes within a sub-volcanic magma chamber between 1982 and 1999, as well as periodic and localized, larger movements, as much as c. 4 m horizontal shortening between 1997 and 1999, with acceleration just before eruption in 1998 (Murray and Wooller 2002).

Regular measurements between 1997 and 2000, using electronic distance measurements (EDM) between established instrument base stations and installed reflectors, showed a maximum distance change of 50 mm day⁻¹ through inflation and deflation (Ramírez-Ruiz et al. 2002). More recently, localized dome growth of 0.3 m over the period of a few minutes has also been demonstrated using sub-pixel offset tracking from very high-resolution optical imagery in 2011 (Walter et al. 2013).

Early monitoring using Interferometric Synthetic Aperture Radar (InSAR) showed no measurable deformation between 2002 and 2010 (Chaussard et al. 2013) but only single interferograms were used, which would have contained a variety of atmospheric effects, and such interferograms are sensitive to only metre-scale changes at best. Later studies using Differential InSAR (DInSAR) and high-resolution X-band SAR (TerraSAR-X, TSX) data indicated c. 50 mm of pre-eruptive deformation in 2013 (Walter et al. 2013). High-resolution TSX Spotlight SAR data have also been used to reveal dome subsidence and small-scale lateral movements; this work also indicated that the strongest subsidence was associated with higher temperatures (Salzer et al. 2014, 2017).

Later work, using Sentinel-1 (STL-1) SAR images, acquired between 2014 and 2016, and Short BASeline (SBAS) processing of data from both ascending and descending orbits revealed subsiding lava flows but also large areas of coherence loss, as described by Lesage et al. (2018) and Cararra et al. (2019). Those researchers noted that lava flow motion contributes significantly to the recorded displacements several months after its eruption and emplacement. They also mapped the deposits and estimated their volumes using digital elevation models (DEMs) and optical images. Cararra et al. (2019) estimated a mean subsidence rate of <129 ± 9 mm a⁻¹ (with a maximum of 182 ± 9 mm a⁻¹), and an average of 26 ± 11 mm a⁻¹ horizontal motion (with a maximum of 63 ± 11 mm a⁻¹ on the steepest slopes) on the 2014 El Zarco lava flow on the southwestern flank of VdeC.

VdeC is included in the list of over 180 volcanoes that are routinely monitored using InSAR and STL-1 data by the Centre for Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET). The regularly produced interferograms from their analysis are made publicly available (https://comet.nerc.ac.uk/volcanoes/colima). They indicate no significant activity, or rather, that any movement there is too small to be detected by DInSAR, and thus VdeC is not currently (as of June 2022) categorized as imminently hazardous.

Gaining quantitative measurements of volcanic eruptive volumes or cooling and contraction rates is considered challenging and, for field data collection, potentially hazardous. Optical imaging is often prohibited by bad weather and clouds of volcanic ash or gas, and extraction of elevation changes, in general, cannot be done accurately enough using photogrammetric methods, although James and Varley (2012) had some success in demonstrating dome height changes in 2010–11 using such methods. Imaging using SAR, on the other hand, allows all-weather operation, and interferometric methods (InSAR) allow precise measurements of topography and ground movement, over variable and potentially long time-frames, and so provide a compelling approach to furthering the understanding of volcanic eruptive processes and the hazards they pose (Wicks et al. 2006; Zhong and Dzurisin 2014). It has been thought that the use of InSAR to provide such information for andesitic stratovolcanoes was not possible because of their steep topography, rapidly fluctuating covering of tephra and, in many cases, densely vegetated slopes, all of which cause low coherence and potentially noisy results.

There have been a number of investigations into post-emplacement lava subsidence at other volcanoes, such as Etna (Briole et al. 1997; Stevens et al. 2001) and Piton de la Fournaise (Chen et al. 2018b; Richter and Froger 2020), all of which demonstrate that lava compaction continues for several years after emplacement and that the detection of persistent lateral movement of flows is an important hazard indicator. Briole et al. (1997) observed that subsidence also extends beyond the extents of the lavas, at rates of at least 12 mm a⁻¹, and that some of the observed subsidence is due to relaxation of the substrate as well as compaction. Pinel et al. (2011) conducted a study of ground motion at VdeC between 2002 and 2006 using Envisat ASAR data. They found no evidence of magmatic injection but detected a small amount of subsidence (c. 10 mm a⁻¹) around the summit and a maximum of c. 40 mm a⁻¹ subsidence around the 1998 lava flows. They also demonstrated some summit deflation and stated that the steep-sided volcano geometry induced strong atmospheric effects in the results.

Obtaining line of sight (LoS) point measurements, from both ascending (east-looking) and descending (west-looking) orbital passes, allows the horizontal and vertical components of ground motion to be geometrically decomposed, provided that the ground motion vectors are at some angle from the orbital track directions (preferably but not necessarily perpendicular to it). Ability to detect LoS ground motion can be severely reduced when the motion direction is dominantly north–south (i.e. nearly parallel to the orbital track direction), but having observations from both orbits can help to constrain the east–west component of movements that have NE–SW and NW–SE orientations. In this case, the lava flows of interest have erupted on the western and southeastern flanks and their movement vectors have been captured well by both ascending and descending orbital passes. Chen et al. (2018a) also performed 2D decomposition of TerraSAR-X and Tandem-X SAR observations from both orbits to derive vertical and horizontal motions at Piton de la Fournaise. They were also able to derive lava thickness using multiple DEMs (with vertical accuracies of a few metres) and observed that although some areas showed strong correlation between thickness and subsidence rate, other areas showed far more complex relationships. They concluded that other factors were of influence in those areas, such as pre-existing structural weaknesses, and that such areas should be the focus of attention for flank-instability hazard assessment. Ebmeier et al. (2012) were able to estimate flow thicknesses at Santiaguito, Guatemala over a 9 year period from DInSAR inversion, with an accuracy of a few metres. Cararra et al. (2019) used high-resolution TerraSAR-X InSAR and Pleiades DEM to estimate average lava flow thickness in the west and SW of VdeC at 17.3 ± 1.5 m and 19.3 ± 1.5 m, with maxima of 51.5 ± 1.5 m and 65.4 ± 1.5 m, respectively.

Long time-series analysis using InSAR can only be done with multi-interferogram persistent scatterer or short-baseline methods (Ferretti et al. 2001, 2014; Hooper 2008). Persistent scatter InSAR (PSI) is excellent for measuring and monitoring small-scale ground motions over long time periods but relies on sufficient stable scatterers over the study time period and is therefore ideal in urban environments where man-made objects provide a dense network of objects that are considered not to change position greatly over time. Performing PSI is very challenging in rural environments where small random changes in the natural scatterers can cause a total loss of coherence between observations. Volcán de Colima therefore presents an extremely challenging target for persistent scatter (PS) techniques because, although the topmost edifice of the volcano is devoid of vegetation and is composed of blocky lavas, there are frequent eruptions and deposition of small amounts of tephra on the surface. The Short BASeline (SBAS) technique allows spatial averaging to alleviate the problem of spatial decoherence, allowing regions to be tracked over longer time series, but it cannot provide the spatial precision that is desirable here. SqueeSAR® provides an advanced, hybrid technique that combines the advantages of both techniques by allowing the characterization and monitoring of both persistent and distributed scatterers throughout the time series. Distributed scatterers can be rocky outcrops, patches of stable soils and regolith, and areas of mature and weathered land surface. This processing method developed and patented by TRE ALTAMIRA (Ferretti et al. 2011), has been developed for monitoring in a variety of situations including volcano monitoring, and has been used here to tackle the problem of small-scale movements and changeable, non-urban surface geometry.

There is a long global archive of optical imagery from the Landsat satellite series, which offers a powerful ability to map and discriminate past lava flow patterns and extents. There are also now many sources of very high-resolution (VHR) optical imagery (offering between 1 and 4 m spatial resolution) allowing very detailed mapping of surface features. However, such data cannot be acquired during poor weather and they cannot provide precise quantitative information about surface change. Routine imaging using SAR, along with DInSAR processing, allows precise and accurate measurement of ground surface change on a centimetre scale and is complementary to optical imagery but cannot provide measurement over long time periods because of temporal decoherence. Such long-term interferometric analysis can feasibly extend back to 1991 with ERS and Envisat SAR data but routine imaging using STL began only in 2015. For long-term monitoring of small-scale ground movements, either PSI or Short Baseline multi-interferogram techniques are needed, in conjunction with multi-temporal stacks of SAR images (>20); such techniques allow measurement of ground deformation with accuracies on a scale of millimetres. Long time-series and bidirectional acquisitions (from both east-looking ascending and west-looking descending orbital observations) allow decomposition and retrieval of both 1D (Line of Sight, LoS) and 2D (vertical and horizontal) ground motions using multi-interferogram techniques. Long-time series, multi-interferogram analysis, multi-sensor SAR temporal image stacks, have been used to demonstrate complex vertical and horizontal spatiotemporal patterns of ground motion at Piton de la Fournaise to highlight the complex hazard potential there (Richter and Froger 2020).

This work aims to demonstrate the effectiveness of multi-temporal DInSAR, using SqueeSAR®, for deriving precise ground deformation measurements in the very challenging environment of an active stratovolcano. The work's objectives are to derive and interrogate vertical and horizontal measurements of ground elevation change over time, to contribute to better understanding of post-emplacement behaviour, cooling and contraction rates, and to demonstrate the potential of the technique, not only for volcanic hazard monitoring but also for revealing small-scale surface expression of a dynamic volcanic system.

The European Space Agency (ESA) Copernicus programme satellites, STL-1A and -1B, capture SAR images every 12 days, with 6 day separation between the two. The satellite image footprint covers some 3600 km², and the two sensors capture C-band SAR images (wavelength 5.3 cm), with incidence angle of 29.1–46.0°, 250 km swath (interferometric wide mode) with three sub-swaths, and with c. 30 m spatial resolution routinely across the globe. The datasets are freely available for download and thus provide a very important and cost-effective resource for hazard monitoring.

This study uses a multi-temporal stack of 147 STL-1A and -1B SAR scenes (from both ascending and descending orbits), acquired between May 2017 and February 2019, that covers a region around VdeC. The SAR data were processed using the SqueeSAR® processing chain at TRE ALTAMIRA. Multi-temporal InSAR measurements are always relative and made with reference to a stable reference point located somewhere within the processed region; in this case the reference point lies near the city of Guzman (19.965°N, 103.464°W) some 24 km from VdeC.

Given the near-polar orbit and right-looking illumination geometry of STL SAR, one might be forgiven for expecting the ascending orbit LoS data to provide the best observations of recent lava flows as these are mainly on the western flank of VdeC, which faces the radar antenna on its east-looking orbital pass. However, the steep-sided summit geometry and the significant foreshortening and layover effects mean that data from the descending orbit (west-looking) produce more coherent measurement points across the western flank (a phenomenon also observed by Kubanek et al. 2015). Although there are fewer points on the western flank in the ascending orbit data, the morphology of the subsidence patterns across the flows is more clearly defined than in the descending orbital dataset.

We also expect that atmospheric effects are significant at VdeC, given the high topographic relief and sizeable area, as confirmed by Pinel et al. (2011). The SqueeSAR® processing chain incorporates the estimation, interpolation and filtering of these spatially correlated and topographic phenomena via its atmospheric phase screen (APS), which is performed after initial candidate point selection and phase unwrapping and before final persistent scatterer/distributed scatterer (PS/DS) identification; this significantly improves the quality of the final PS/DS point selection. Estimation of the APS also relies on the setting of a common reference point, which should be in a stable area and where the signal-to-noise ratio is high. At this position, the unwrapped phase contains residual topography and changeable atmospheric effects, which are then separated using a stochastic model. The stable reference point chosen in this dataset is located on the eastern side of the city of Guzman to the NE of VdeC. Guzman is located on the margin of the active rift and regularly experiences significant fault-related ground subsidence; thus, the reference point is located on the stable eastern footwall block, which is not moving significantly (shown in Fig. 1b).

The NASA Shuttle Radar Topographic Mission (SRTM) global DEM v3 1 arc-second product (gridded at 30 m), which is freely available for download, has been used for this work. From it, a shaded relief display has been generated to help visualize the map and morphology of VdeC eruptive products and to geolocate the SqueeSAR® results.

Observational cameras are sited at distances around the volcano and provide regular footage of the summit to catch fumarolic and eruptive activity, and to monitor barrancas (ravines) for lahars. The University of Colima publishes weekly bulletins on the current activity for the volcano (freely available at https://portal.ucol.mx/cueiv/; bulletins for the period in this study are available upon request). These also include webcam images where relevant (i.e. when tremors are recorded that are accompanied by fumarolic activity at the volcano's summit).

The unique SqueeSAR® processing approach allows the exploitation of both distributed and persistent scatterers (DS and PS) and can detect non-linear deformation. SqueeSAR® combines the benefits of PSI and SBAS approaches. In addition to strong PS point scatterers, DS are identified through a spatial adaptive filter to find statistically homogeneous adjacent pixels across a pixel neighbourhood (Fig. 4a). The DS processing employed in the SqueeSAR® algorithm is different from other SBAS algorithms because information from all possible interferograms in the time series as well as the resulting coherence matrices are used and because the DS phase values are retrieved before any phase unwrapping is performed (Ferretti 2014).

One of the main challenges in using multi-interferogram techniques is understanding the accuracy of the measurements. Because all InSAR measurements are made relative to a reference point, the absolute accuracy is dependent on the reliability of the reference point. Because the results of differential interferometry show relative change rather than absolute values, the reliability of measurement points is best assessed by the precision of their measurements. Each measurement point in a SqueeSAR® dataset has individual precision values that are determined by the measurement point's scattering properties. The mean precision of an entire dataset is influenced by many factors; for example, the number of images in the data stack. The precision of the measurement point's geolocation, and the precision of the measurement point's displacement time series both need to be considered. The error on a measurement point (MP) c. <1 km from the reference point, derived from a stack of >30 SAR scenes, for a 2 year period, processed using SqueeSAR®, is widely stated as having a standard deviation of <1 mm a⁻¹; a single displacement value may have a precision of ±5 mm (Ferretti 2014). A full description of the processing approach of the SqueeSAR® method is not within the scope of this paper and more detailed explanations have been given by Ferretti et al. (2001, 2011).

The SqueeSAR® processing of these data has generated ground movement point data of exceptional quality for such a challenging environment. Processing of the SAR scenes from ascending and descending stacks respectively provides excellent coverage of PS and DS points measuring 1D LoS ground displacement values (velocity) on the unvegetated upper slopes of the volcano. Variable point densities, as high as 8500 km⁻², were achieved across the barren upper edifice of VdeC (see Fig. 5), and the coherence threshold of 0.6 has meant that no PS or DS points were retrieved across the densely vegetated lower slopes.

Use of single orbit data provides 1D LoS measurements of the ground. Use of data from both ascending and descending SAR orbits means that the 2D vertical and horizontal components of movement can be decomposed geometrically. The decomposition requires the same target to be visible from both ascending and descending orbits, which is not always achievable. The geometries of the LoS observations from opposing orbits are shown in Figure 4b. The LoS observations describe downward ground motion, away from the sensor, as negative values (subsidence) and upward motion as positive values (heave). If the two observed movements are negative in both the ascending and descending observations, the net movement in vertical direction can be derived trigonometrically, as being downward. If motion is positive in ascending and negative in descending observations, the derived net motion will be westward, and so on.

Because the SAR acquisition dates and PS and DS point positions are different in the results of ascending and descending image stack processing, the displacement measurements must be interpolated and resampled to produce a ‘pseudo-PS’ (and -DS) displacement map, so that the 2D vertical and horizontal components of displacement can be calculated. For SqueeSAR® processing of Sentinel-1 data the spatial sampling is c. 80 m (100 m in this study) and temporal sampling is 12 days (Bischoff et al. 2020). The spatial and temporal detail of the 2D data are thus reduced slightly but the separation of vertical and horizontal movement components is vital to better understand the nature and likely causes of movements.

Analysis of both 1D (single orbit) LoS and 2D vertical displacement results show a number of very clear patterns across the edifice of VdeC. The cumulative ground displacement over the 2 year period is illustrated as a series of time-slice maps of ground measurement point values that have been interpolated using a simple inverse distance weighted average (IDW) method, shown in Figure 6. Calculation of the relative total change over the 2 year period shows that almost every measured location (within the exposed edifice) has subsided by at least 5 mm, and that there are no areas of significant heave. The most noticeable ground movements are located over the most recent lava flows, which erupted between 2014 and 2016; these all show steady and pronounced subsidence between 2017 and 2019, several years after emplacement (Fig. 7). The subsidence curves of these flows demonstrate very well the clear correlation between flow dimensions and subsidence rate. The 2014 El Zarco and Playon lava flows are of similar size and volume (El Zarco slightly larger) and their subsidence rates are very similar, almost linear, with Playon reaching −60 mm and El Zarco reaching −70 mm over 2 years. The mapped subsidence patterns of these two flows, shown in Figure 6, show some spatial variations, with the areas of peak subsidence being situated in the upper and centre parts of the Playon flow and near the flow-front of the El Zarco flow; differences that may be a function of localized pre-flow topography. The Crater N flow is very short and shows the least subsidence (maximum −60 mm), and the 2015 Montegrande flow is the thickest (of these recent flows) and hence has the greatest potential to subside (maximum −120 mm). Subsidence profiles along the lengths of these four flows, approximately parallel to flow direction, are shown in Figure 8. The profiles reveal the variability of subsidence both between and along the flows. These spatiotemporal variations hint at thickness variations, caused by infilling pre-flow topography, as the main control on subsidence variation along the flow.

The most prominent area of subsidence lies within the extents of the 2015 and 2016 Montegrande lava flows on the southern flank of the edifice, both of which flowed southwards from the crater and overlie one another, the volume of the 2016 flow being much smaller than that of the 2015 one. The higher subsidence rates in the centre of the 2016 lava flow may be explained by overlapping cooling–contraction patterns of the two flows, which overlie one another, although the area of 2015 flow that is exposed beyond the southernmost extent of the 2016 flow appears not to have subsided significantly at all, suggesting that the main area of subsidence is occurring entirely within the later 2016 flow. The subsidence could be explained by gravity consolidation of the loose material of the 2015 PDC, under the weight of the lava flow, but this PDC deposit was largely erosional on the steeper slopes, where the lava flow was subsequently emplaced; thus there is very little material beneath and so this really cannot explain the contraction occurring in the 2015 and 2016 lava flows.

There are several other less obvious but notable features. The interior of the crater can also be seen to subside steadily, by more than 40 mm a⁻¹, beginning on the western side but gradually affecting the whole crater area. The area around and between the four most recent flows is also subsiding, less so than the flows, but by more than 40 mm a⁻¹ in some places, and this could be caused by the loading effect of the adjacent more recent lavas. There is also a localized (isolated) area of subsidence, located lower down on the western flank, spatially between the 2013 and 2014 flows, which is subsiding by >20 mm a⁻¹. The reason for this is unclear but it may represent some older deposits (see also Fig. 2a), which may still be cooling and contracting or compacting. Similarly, a small area in the SE of the edifice that is also subsiding by as much as 40 mm a⁻¹ is coincident with an area of lavas erupted in the mid-1970s. Examination of the coherence values for points in these areas (Fig. 4a and b) indicates that although point densities are lower in this area, and there are gaps in coverage caused by the shadowing effect of the summit, and there is no reduction in the reliability of the measurements, coherence values are mainly above 0.8. A further area of minor subsidence (<40 mm a⁻¹) lies to the NW of the crater and coincides with a lava flow erupted in 2004.

The patterns of vertical and horizontal movements are shown in map form in Figure 9 and as two average time series plots in Figure 10. Figure 9a shows the regular distribution of the interpolated points derived by the decomposition of vertical and horizontal components of ground motion. Figure 9b shows the magnitude of the vertical movements and indicates that thermal subsidence is constrained mainly but not entirely to the area of the four most recent lava flows, and that the upper part of the edifice (above 2200 m) is subsiding by at least 5 mm a⁻¹.

The horizontal component of movement is shown in Figure 9c and d. The interpolated point locations are presented as directional symbols in Figure 9c (triangles pointing left and right to indicate western and eastern movement respectively), which show clearly that there are very few point locations where there is no horizontal movement at all. Some areas are moving eastwards, and these are largely constrained to the lower western and eastern slopes. Figure 9d shows that the entire upper western flank of the VdeC edifice (above c. 2700 m elevation) is moving laterally (westwards) by nearly 20 mm a⁻¹ and that the movement is not confined to the subsiding lava flow extents but also to the areas between them. This is a clear indication of a net westward subsidence of the western flank. The relative contributions of horizontal and vertical movements equate to a slope of c. 70° and this appears consistent across the western flank of VdeC. Such a high angle would be greater than any reasonable angle of repose for the talus material on the western side of the cone, and thus suggests that the movements seen here are caused by in situ subsidence of the cone rather than by any distributed rockfalls.

The time series of individual PS or DS points around the summit and crater of VdeC reveal several different and interesting patterns, as illustrated in Figure 11. The small area to the NW of the crater, including part of the crater rim, displays fairly steady subsidence of c. 50 mm over the 2 year study period (>20 mm a⁻¹) but there are several short-term upward perturbations to that nearly linear trend; that is, short periods of heave lasting a few weeks to a couple of months; these are illustrated in Figure 10a and c. At first these short perturbations may appear to be noise or random fluctuations, but closer examination reveals that they have consistent character and timing, and they can be observed over a significant area. The two most prominent examples show a gradual rise in elevation, of nearly 10 mm, and a sharp fall back to the dominant subsiding trend. The first begins on 31 July 2017, peaks around 24 August and subsides sharply between 28 and 30 August 2017. The second begins on 23 October 2017, peaks around 22 November and subsides suddenly between 25 November and 4 December 2017. Further, very short-lived heave ‘spikes’ occur between 11 and 21 May 2018, and again in July 2018 (see Fig. 12). During this period the degassing from the volcano was generally low (from below the detection level of about 0.1 kg s⁻¹ to 4.2 kg s⁻¹) with some pulsing observed. The highest fluxes were recorded during February to April 2017, with a small increase in October and December 2017 and February 2018. Flux measurements are generally carried out for a few hours on a weekly basis using UV spectroscopy, which was insufficient to determine any correlation between the deformation and gas pulses. The SOMA station, which is situated 1.7 km from the crater on the NW flank of the volcano, recorded generally low seismicity but an increase was detected toward the end of August 2017. The real-time seismic energy measurement (RSEM) and cumulative RSEM for the 8 month period between May and December 2017 (Fig. 13) shows that seismicity gradually increased toward the end of August 2017. On 28 November 2017, a large-amplitude high-frequency volcanic earthquake was recorded, which is represented by the highest RSEM value of that month. There were further M_w 8.2 and M_w 7.1 earthquakes in southern and central Mexico, and these correspond to the maximum RSEM values in September 2017. The weekly seismic bulletin No. 45, dated 24 November 2017, also states that some ‘27 high frequency events, 8 long-term events, 1 h of tremors and 3 landslides’ were observed and recorded at that time.

These same patterns can be observed in 1D LoS measurements of both ascending and descending orbits, although their magnitude is higher in the ascending data. Examination of the 2D vertical derivation reveals that these short movements have a clear vertical component of heave and collapse but have also a horizontal component; in this area, the ground moves westward (outwards) at the same time as it heaves and then moves eastward (inwards) as it collapses downward (see Fig. 13). The temporal correlation of subsidence, heave and lateral motion and their coincidence with evidence of seismic and pulsing fumarolic activity, lead us to postulate that these short-term perturbations may be caused by gas build-up and subsequent degassing events, as others have suggested elsewhere (Walter et al. 2013), or by fluid migration (Petrillo et al. 2019). However, during phases with no significant explosivity or effusion, and low gas fluxes (such as 2017–19), despite the possibility of complex variations in permeability and gas flow (Kolzenburg et al. 2019), it is difficult to attribute the deformation to variations in fluid flow. It is also possible that such periodic changes could be caused by high-frequency small movements of magma (Chen et al. 2018a). Concentration of shear stress can be produced by magma motion, and modelling has shown that even low levels can produce edifice tilt and deformation (Marsden et al. 2019), which may be detectable at surface. Ascent of magma through the conduit, around the ductile–brittle transition depth, can also trigger seismicity (also shown by Varley et al. 2010) and can cause both upward and outward movement of the edifice.

On the eastern flank of the summit, the time series plots show a more gradual but repeated pattern of subsidence and heave, and over a much larger area (see Fig. 11b and c). These movements are observable only in the LoS measurements from the ascending orbit. There are very few descending orbit LoS measurement points here because of a loss of coherence over the time series, partly because of the shadowing effect of the summit but possibly also from shifting materials on the ground. This area experienced subsidence of c. 10 mm a⁻¹ (between May and August 2017), followed by heave and recovery back up to the same elevation as the start (by mid-March 2018), then subsidence of nearly 20 mm a⁻¹ (by September 2018) and finally recovery again back up to the starting elevation by February 2019. The reason for this longer term apparent cyclicity is unclear at this stage but may be related to gas build-up and release or to localized magma injection or movement within the shallow magma plumbing system and thermal subsidence. Infiltration, after strong seasonal rains, has the potential to influence these subsidence patterns through increasing surface movement of deposited pyroclastic materials. Examination of highly variable rainfall patterns across the period of this study reveals no apparent correlation with ground motions seen in this study, the majority of which are strongly linear. The smaller scale perturbations on the NW summit (Fig. 11a and c) occur over too long a period to be specific weather events or atmospheric effects (occurring over several days to several weeks), as the latter are usually spatially broad and temporally variable, and over too short a time period to be controlled by seasonal rains.

These ground movements are clearly complex in four dimensions and too complex to reveal using 2D time series illustrations alone. Time slices of measurement point ground movement values have been produced from the ascending orbit LoS data (to allow visualization of patterns on the NW crater and eastern summit flank). These slices of ground measurement values were then interpolated using a simple IDW method to produce time series maps of the summit region. These are shown in Figure 14, to allow comparison of the spatiotemporal relationships between the ground motion patterns and the lava flow boundaries and the crater extent. It can be seen that, in addition to the steady subsidence of the Montegrande and Crater North lava flows, several areas move upwards and downwards repeatedly; these movements are not always related to lava flow boundaries and their causes are unclear.

From the SqueeSAR® 1D LoS results, measures of ground movement acceleration have also been generated and are illustrated in Figure 15. These reveal whether negative displacement rates are slowing down (positive acceleration) or are speeding up (negative acceleration). Black continuous outlines in Figure 15 indicate the extents of the four recent lava flows, and white dashed outlines (numbered 1–4) represent the areas used to derive the average acceleration time series plots in Figure 16.

The results for VdeC are complex, they differ slightly between ascending and descending orbits and there are several patterns of note. Ground movement around the upper part of the edifice, around the crater (Fig. 15, white outlines 1 and 3), which has been shown to be subsiding (i.e. has negative velocity), also has negative acceleration (red) and thus the subsidence appears to be speeding up; this appears most noticeably on the western side of the summit and in the ascending orbit data. Conversely, the acceleration patterns in both ascending and descending orbit measurements indicate that the 2015 and 2016 Montegrande lava flows (Fig. 15, outline 2), which have also been shown to be subsiding significantly, have positive acceleration (blue), indicating that the subsidence is slowing down, as is to be expected some 4–5 years after emplacement. A second smaller area of positive acceleration and slowing subsidence is visible in the descending orbit data (Fig. 15b, outline 4). Comparison of the positions of areas 1 (speeding up) and 4 (slowing down) with the lava flow extents reveals that these patches of differing acceleration, although in close proximity, are not in exactly the same location. The patterns differ between ascending and descending orbit data because of slight differences in slope geometry with respect to the incidence angle. The reasons behind these spatially localized variations in ground motion behaviour remain unclear.

The western summit (Fig. 16, 1, green dots) lies on the uppermost part of the 2014 El Playon lava flow and shows an acceleration curve that is steady and becomes downward sloping (convex upward) after August 2018, indicating that subsidence began very gradually and then started to speed up. The Montegrande flows (2, blue dots) show an acceleration curve that is steady but gently concave upward, indicating that subsidence is gradually slowing down. The area on the eastern summit (3, red dots) shows a complex pattern of repeated periods of slowing and speeding up, which are temporally correlated with the periods of heave and subsidence shown in Figure 11c, but overall and relative to the starting date the acceleration is negative. The fourth area of interest on the western flank (4, cyan dots), which appears only in the descending orbital acceleration data (Fig. 15b), shows steep and steady acceleration at first and appears to be slowing slightly toward the end of the time series.

Although these data reveal complex, precise and detailed information about the subtle and dynamic centimetre-scale surface changes of this active volcano, the usefulness of the work is limited by the relatively short period of the time series. A longer time series, of perhaps 10–15 years, is desirable but there is no single SAR sensor that can currently provide such a dataset on which to perform multi-temporal DInSAR analysis at comparable resolutions and with the same wavelength. ERS-1 and -2 SAR sensors provide coverage from 1992 to 2000, and Envisat ASAR from 2002 to 2010, whereas Sentinel-1 became available only at the end of 2014. There are other SAR sensors that provide higher spatial resolution and shorter wavelength but extracting reliable millimetre-scale measurements over a long time period from different sensors with different orbital geometries, although achievable, is problematic and relies on regular coverage by such sensors throughout the entire time period. These small-scale measurements of ground motion would be even more compelling if they could be combined with conventional ground levelling measurements at some fixed reference points to provide independent validation. Availability of an accurate, relatively high-resolution DTM (with vertical accuracy of a few metres) that predates the recent eruptions would have allowed the estimation of flow thickness, as well as correlation with likely cooling and contraction rates and with the evidence of thermal subsidence provided here.

Although the lava flows of the most recent eruptive episode (2013–17) show dramatic evidence of cooling, contraction and post-emplacement flow (Fig. 17), evidence for longer-term persistent movement can still be discerned for lavas and pyroclastic flows emplaced in 2004 (16 years following emplacement) and potentially even the 1975–76 flows (nearly 40 years following emplacement). This demonstrates that although these lava flows may have cooled and solidified at the surface, they may remain plastic and liable to creep on the large scale for decades following emplacement. This might have important implications for volcano flank stability. Minor rockfalls occur frequently across zones of the western and southwestern flanks of the VdeC edifice (Mueller et al. 2013) and individual block scatterers may shift, topple or fall. If ‘rockfall’ is a continuous process it might appear, over time, as though the surface was moving, when there are simply many very small rockfalls over an area. Such rockfalls are most prevalent on the steepest slopes near the crater. Therefore, if the apparent subsidence signal could be explained by rockfall over an area, we would expect to see a higher subsidence rate near the crater, and a lower rate or no apparent subsidence at the foot of the flow, and perhaps an apparent heave around edges of the flow (as rockfalls accumulate material on the edges and show an average rise there). This is not what we observe in these data, and gradual subsidence provides a better explanation of the patterns across the edifice of VdeC. If blocks were to move across a substantial area, by more than a couple of centimetres (half the wavelength), the coherence would be lost between SAR observations, whereas the coherence has been shown to be very high throughout the time series. Periodic and irregularly spatiotemporally distributed rockfalls are more likely to be the cause of the small-scale variability of the measurement signal or ‘noise’ that we do see on individual measurement points. In general, and especially when averaged over an area, these data show patterns of very steady, persistent subsidence, or of periodic heave and subsidence, over the time period of the study.

Deformation monitoring is a key tool to monitor volcanic activity, from small changes caused by minor explosions and variations of fluid migration within the plumbing system, to broader deformation of the edifice during eruptive periods. The results here paint a picture of a very dynamic environment for an apparently quiescent period for VdeC. The short-term heave and subsidence of the crater rim during this period may indicate that gases or fluids preferentially find pathways to escape at the rims of the crater. Field evidence of currently active gas pathways is provided by the abundance of fumaroles on the northwestern, northeastern and southeastern edges of the crater, in contrast to their relative absence in the centre. This pattern was also postulated by Lavallée et al. (2012) during the dome-forming eruptions of the 1998–2011 eruptive phase, in which a central solidified plug in the crater forced gases to escape around the crater margins. Such degassing patterns have not decreased or changed greatly since the 1990s. The nature and origin of the longer amplitude ground motions detected around the crater, and on some flanks, are uncertain because of a lack of clear activity at the surface, but it can be speculated that these movements may relate to injection of gases or to the movement of potentially fresh magma within the shallow plumbing system. Such suggestions may be confirmed with further interpretation of more detailed seismic data and petrological work when the volcano next erupts.

One of the most important findings of this work is that the western part of the edifice of the volcano appears to be both subsiding and moving westwards by upwards of 20 mm a⁻¹. This is potentially significant for VdeC as flank collapse is a common hazard at Colima Volcanic Complex stratovolcanoes (Robin et al. 1990), and the current edifice is thought to have undergone partial collapses twice in the last c. 3.5 kyr (Cortés et al. 2019) as well as being at risk of block and ash flow inundation of nearby areas (Saucedo et al. 2005; Sulpizio et al. 2010). Although more regular volcanic activity can significantly affect the smaller villages surrounding VdeC, the larger settlements such as Colima City are built upon former debris avalanche deposits from flank collapses and such collapses could pose a significant risk to these larger settlements to the south. Even a partial collapse, similar to the collapse of the dome and partial breach of the crater rim during the July 2015 eruption, can create long-runout pyroclastic flows, which would be a significant threat to settlements such as La Yerbabuena. The finding of subsidence across the entire western flank of the volcanic edifice merits further analysis to assess the potential threat of flank collapse, and necessary hazard planning and mitigation, should lava and pyroclastic flows preferentially flow towards the west.

We conclude that, in general, insufficient InSAR analysis has been conducted at VdeC over the last decade, particularly when the volcano was actually active. The analysis presented here for the period after the last eruptions (May 2017–February 2019) provides some tantalizing evidence of small-scale, short time-frame fluctuations of VdeC, rather as though the volcano is ‘breathing’ up and down. To fully demonstrate such activity here, a longer time series would be needed, which spans the entire period of activity from at least 2011 (the end of the previous phase) to the present; such a time series would fill the gaps between existing rather focused studies and will form one of the future developments of this work. Combination of different monitoring tools (SqueeSAR®, seismic monitoring, gas monitoring), along with a time-constrained understanding of the magmatic activity from petrological studies, will provide an integrated and powerful means to understand volcanic hazards at VdeC.

The small-scale, incremental ground movements demonstrated here are measurable only using DInSAR time series analysis, such as via SqueeSAR®, PSInSAR or SBAS. These movements would be undetectable in single interferograms (and barely detectable in differential interferograms) between SAR observations separated by a few days, such as are published by COMET in their routine volcano monitoring service. Similarly, interferograms made between SAR observations separated by longer time-periods (e.g. several months or years) would probably show significant or total coherence loss because of shifting ground materials and/or tephra deposition. Interferograms have been published by COMET that show beautiful fringe patterns representing a few centimetres of motion, between 3 December 2019 and 26 January 2020. These movements appear to affect the entire volcano (Vdec and Nevado de Colima) systematically, but within these large-scale patterns the small-scale lava subsidence patterns on VdeC are obscured.

These SqueeSAR® results, based on routinely acquired and freely available Sentinel-1 SAR, and detailed analysis presented here, have revealed a series of spatially and temporally complex patterns of behaviour at one of the world's most active volcanoes. The SqueeSAR® processing has generated data of unprecedented quantity and quality in this highly challenging rural environment where persistent scatterers are in short supply. This short 2 year time series shows that, although VdeC is not currently classified as in eruption, there are ground movements occurring at spatial, temporal and displacement scales that are measurable effectively using multi-temporal DInSAR. These movements reveal that the recently erupted lavas are still cooling and contracting, that the summit and crater are undergoing periodic vertical and horizontal adjustments, caused either by gas build-up and release or by magma motion within the plumbing system, and that the entire edifice of the volcano is subsiding very slightly.

Our results show steady thermal subsidence of recent lava flows (more than 40 mm a⁻¹), which is now slowing. We can also see westward and vertical movement of the entire western flank and the recent lavas there, by c. 20 mm a⁻¹. There is periodic minor and localized uplift and occurrence of subsidence events, which may be correlated with degassing and collapse at the crater, but which may also be contributed to by magma motion in the conduit. The significance of short-term patterns is not well understood, but this work provides evidence of subsurface activity that is detected by webcams. Petrological analysis and real-time in situ data are clearly needed to assist interpretation. In view of the apparent switching of recent eruptive activity towards the western flank, there are implications for geohazard warnings regarding potential sector or flank collapse. The last major sector collapse occurred around 3600 yr BP (Cortés et al. 2010, 2019) and it has been suggested that even a minor flank collapse could threaten local populations; therefore further work to assess the causes of these fluctuating ground movements is needed to improve regional geohazard assessment.

This work demonstrates the value of retrospective analysis of any active but poorly monitored volcano in a remote area, without use of in situ instrumentation. This approach also has the potential to allow the linkage of surface evidence of change to deeper magmatic processes, but this needs further work or application to a different volcano, where considerably more geodesic instrumentation and data exist.

We are extremely grateful to TRE ALTAMIRA for processing the Sentinel-1 SAR time series using the SqueeSAR® engine and making them available for this analysis. The provision of RSEM data and description by D. V. Bracamontes, at the Universidad de Colima in Mexico, is gratefully acknowledged.

PJM: conceptualization (lead), data curation (lead), formal analysis (lead), investigation (lead), methodology (lead), project administration (lead), visualization (lead), writing – original draft (lead), writing – review & editing (lead); CAB: conceptualization (supporting), data curation (supporting), formal analysis (equal), methodology (supporting), writing – review & editing (supporting); GH: formal analysis (supporting), methodology (supporting), validation (supporting), writing – review & editing (supporting); CMP: conceptualization (supporting), methodology (supporting), validation (supporting), writing – review & editing (supporting); NRV: writing – review & editing (supporting); GN: data curation (supporting), investigation (supporting); AF: writing – review & editing (supporting)

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

The data that support the findings of this study are available from TRE-ALTAMIRA but restrictions apply to the availability of these data, which were used under licence for the current study and so are not publicly available. Data may be made available from the authors upon reasonable request and with permission of TRE-ALTAMIRA.