Challenges of determining frequency and magnitudes of explosive eruptions even with an unprecedented stratigraphy

Abstract

Through decades of field studies and laboratory analyses, Volcán de Colima, Mexico has one of the best known proximal eruption stratigraphies of any volcano, yet the frequency and magnitudes of previous eruptions are still poorly resolved. Hazard assessments based on models of well-known, well-mapped recent eruptions may appear to have low uncertainty, but may be biased by the nature of those events. We present a comprehensive stratigraphy of explosive eruption deposits combining new data collected as part of this study together with published and unpublished data. For the first time we have been able to model five of the best exposed and cross-correlated pre-historical Holocene explosive events at Volcán de Colima. By modelling the volumes and magnitudes of Holocene eruptions at Volcán de Colima, we are able to improve estimations of the potential range of magnitudes of future explosive eruptions, which can be incorporated into hazard assessments for nearby communities. Based on recent studies we demonstrate that these volumes may be underestimated by at least an order of magnitude, and show that even with an exceptionally well-defined stratigraphic record our understanding of the full range of explosive eruptions may still be biased.

Introduction

Volcán de Colima, situated at the western end of the Trans-Mexican Volcanic Belt (Fig. 1), is one of North America’s most active volcanoes (Cortés et al., 2010; Luhr et al., 2010). Colima is a typical arc stratovolcano with recent activity characterised by andesitic lava dome growth with numerous small, degassing eruptions, punctuated by larger Vulcanian-style explosions causing dome collapse and block-and-ash flows (Luhr and Carmichael, 1980; Medina-Martínez, 1983; Robin et al., 1987; Luhr and Carmichael, 1990b; Savov et al., 2008). In contrast, Volcán de Colima has experienced sub-Plinian and Plinian explosions throughout the Holocene. Since records began in 1519 with the arrival of the Spanish Conquistadors, there have been nine sub-Plinian or Plinian events in 1585, 1590, 1606, 1622, 1690, 1806, 1818, 1890 and 1913 (Medina-Martínez, 1983; De la Cruz-Reyna, 1993). The most violent of these occurred in 1913 and prior to that, in 1818. The 1913 and 1818 Plinian eruptions generated explosions that reportedly deposited ash in Guadalajara, 130 km to the north of Colima, and in the town of Saltillo, 725 km northeast of Colima (Waitz, 1915). Collapse of the eruption column generated pyroclastic density currents (PDCs) that travelled 15 km down ravines in the flanks of the volcano (Bretón González et al., 2002; Saucedo et al., 2005). During the 1913 eruption, the top 100 m of the current stratovolcano edifice was destroyed and left a crater 350 m deep (Waitz, 1935; Bretón González et al., 2002; Luhr, 2002). Based on the activity following the 1818 and 1913 explosions, Luhr (2002) identified an approximately 100-year cycle culminating in a major explosive Plinian eruption. Luhr (2002) suggested that Volcán de Colima is currently in the third eruption cycle.

Fig. 1

Location map of Volcán de Colima, the southernmost volcano of the Colima Volcanic Complex (CVC), in the western Trans-Mexican Volcanic Belt (TMVB) after Crummy et al. (2014). Volcanism along the western TMVB is related to the subduction of the Cocos and Rivera plates under the North American plate. Volcanoes are: Tq, Volcán Tequila; Cb, Volcán Ceboruco; Tp Volcán Tepetiltic; S, Volcán Sanganguey; LN, Las Navajas; SJ, Volcán San Juan. Cities: C, Colima; G, Guadalajara; M, Manzanillo

Tephra fallout and pyroclastic surge deposits are exposed on the adjacent but extinct Nevado de Colima stratovolcano (Luhr et al., 2010) (Fig. 2). Through detailed radiocarbon dating (181 ages) and stratigraphic correlations, Luhr et al. (2010) identified at least 25 Plinian deposits spanning the past 30,000 years, and defined a detailed stratigraphy of explosive eruptions comprising five Plinian fallout deposits erupted between ∼5000 years before present (yrs BP) and the present day. Further field campaigns carried out as part of this study have built on the work of Luhr et al. (2010). Based on field relationships, tephra geochemistry and in situ mineral chemistry (Crummy et al., 2014), we have extended the detailed explosive eruptive stratigraphy back through the Holocene and Late Pleistocene and have identified 15 distinct tephra fallout deposits that are separated by ash-rich PDC and surge deposits that erupted between ~ 13,000 years BP and the present day (Fig. 2). During this period there have also been at least four gravitational collapse events resulting in large (1 to > 10 km³) debris-avalanche deposits to the south, southeast and southwest of Volcán de Colima dated at approximately 2500, 3600, 7000 and 9600 yrs. BP (Robin et al., 1987; Luhr and Prestegaard, 1988; Luhr and Carmichael, 1990a; Siebe et al., 1992; Stoopes and Sheridan, 1992; Navarro et al., 1994; Komorowski et al., 1997; Cortés et al., 2009; Cortés et al., 2010; Norini et al., 2010; Roverato et al., 2011). The source of the pre-Historical Holocene deposits is assumed to be the summit vent of Paleofuego de Colima which suffered a sector collapse at approximately 2500 yrs. BP (Siebe et al., 1992; Navarro et al., 1994; Komorowski et al., 1997; Cortés et al., 2010). The currently active cone resides in the centre of the old Paleofuego crater (Luhr and Carmichael, 1990a; Cortés et al., 2019).

Fig. 2

Stratigraphic correlations across the flanks of Nevado de Colima, from Luhr et al. (2010; VF95-XX to VF02-XX) and this study (VF10-XX and VF11-XX). Section locations are given in Luhr et al. (2010) and Crummy et al. (2014). Locations are labelled following Luhr et al. (2010): VF, year, location number

The pre-Historical Holocene units have greater thicknesses, providing geological evidence of significantly larger magnitude Plinian eruptions than the 1913 eruption, on which current hazard maps are based. These deposits, and those of the 1913 and 1818 eruptions, are well studied in terms of their chemistry and petrology (e.g. Luhr and Carmichael, 1982; Robin et al., 1987; Luhr and Carmichael, 1990a; Robin et al., 1991; Luhr et al., 2006; Luhr et al., 2010; Saucedo et al., 2010; Crummy et al., 2014; Massaro et al., 2018; Crummy et al., 2019). Nevertheless, the characteristics of the eruptions that produced the pre-Historical Holocene tephra fallout deposits were unknown until this study.

Here, we present an expansion of the Luhr et al. (2010) stratigraphy of explosive eruptions from approximately 5000 to 13,000 yrs. BP, based on field campaigns in 2010 and 2011. For the first time at Colima, we present calculations of the volumes and the dispersion of tephra for five of the best-exposed and cross-correlated Holocene tephra fallout deposits sourced from the Paleofuego de Colima vent. These fallout deposits were erupted between 6000 and 4400 yrs. BP (conventional C¹⁴) (Luhr et al., 2010). This combination allows us, for the first time, to assess the completeness of the stratigraphic record, and better constrain the potential maximum magnitudes of future explosive eruptions at Volcán de Colima, which could be incorporated into hazard assessments.

Methods

During two field campaigns in 2010 and 2011, we revisited numerous sections described by Luhr et al. (2010) and collected data from eight new locations across the flanks of Nevado de Colima volcano (Fig. 2). Deposit thickness and granulometry data for tephra dispersion modelling were compiled from unpublished data collected by James Luhr, Carlos Navarro and Ivan Savov during the 1990s and 2000s, and new data collected as part of this study. The thickness of each of the five eruption deposits has been measured at 20 to 27 localities (Table 1) and granulometry samples for grain size distribution collected from 2 to 8 localities across the flanks of Nevado de Colima. Granulometry samples were weighed, and then sieved separating the sample into fractions from − 4.5 phi to 4.25 phi (22 mm to 53 μm). Each size fraction was then weighed to give the grain size distribution at each locality. The granulometry and thickness data are provided in the Additional file 1.

Table 1 Field characteristics of the pre-Historical Holocene explosive eruption deposits

We have calculated the eruption column height, total erupted mass and tephra fallout maps for individual units P, S, U, W and Y (erupted between ca. 4400 and 6000 yrs. BP; Table 1) and units W and Y combined, assuming deposition during a single eruption, using the Tephra2 dispersion model (Bonadonna et al., 2005; Connor and Connor, 2006). Individual deposit thicknesses of the five units and the combined thicknesses of units W and Y were input into the model. Tephra2 is an advection-diffusion model based on the work of Suzuki (1983) that describes diffusion, transport and sedimentation of tephra particles released from an eruption column (Connor et al., 2001; Bonadonna et al., 2005). It calculates the total mass per unit area (kg m^− 2) of tephra accumulation at individual grid locations by solving a simplified mass conservation equation. The mass conservation equation takes into account the distribution of tephra mass in the eruption column and particle settling velocity, as well as horizontal diffusion within the eruption column and atmosphere after the particle has been ejected from the plume (Connor et al., 2001; Bonadonna et al., 2005; Connor and Connor, 2006).

To account for variations in wind conditions, ten years of wind data for Colima were downloaded from the National Oceanic and Atmospheric Administration (NOAA) National Centre for Environmental Prediction (NCEP) Reanalysis2 project (NOAA, 2015), from 1st January 2006 to 1st January 2016, for 12 wind levels from 4 km to 30 km a.s.l. NOAA Reanalysis wind profiles were used by Tephra2 to find the best fit to the field data. For each unit, the inversion was run for a dataset of wind over a year (2015), sampled 4 times daily (1460 files). The results were analysed and sorted according to the best fit, and those that were geologically inaccurate were discarded. Out of the remaining results, for consistency, the top 100 best fit were used to calculate the median. Tephra fallout thickness maps were generated by forward modelling the best-fit results.

Hazard curves were calculated using the output parameters from Tephra2 and the ten-year NOAA Reanalysis2 wind database to calculate the probability of mass loading exceeding a specific value at a specified location. The ten-year NOAA Reanalysis2 wind database was used to account for variations in wind conditions.

Results

Stratigraphy

The 1913 tephra fall deposit forms the topmost layer of the Holocene stratigraphy exposed on the flanks of Nevado de Colima. Luhr et al. (2010) initially described this unit as unit α, measured at 37 locations with variations in thickness from 60 cm to 6 cm at distances of 3 to 12 km from the source (the current Volcán de Colima vent). However, in some localities, unit α comprises two or three pumice-rich horizons separated by ash-rich surge deposits. Based on geochemistry and petrology, Luhr et al. (2010) identified two distinct tephra fall deposits within unit α which they interpreted as resulting from the 1818 and 1913 events. The 1818 tephra fall deposit is not preserved in all localities where the 1913 unit is exposed. Stratigraphically below unit α is a fine-grained ash deposit (unit Z) varying in thickness from 30 cm to > 2 m, measured at 24 localities (Luhr et al., 2010). Detailed charcoal sampling and subsequent radiocarbon dating revealed that this unit represents ca. 4000 years of eruptive activity at Paleofuego de Colima (Luhr et al., 2010; Table 1) including two of the gravitational collapse events at approximately 3600 and 2500 yrs. BP (Fig. 3).

Fig. 3

Uncalibrated radiocarbon C¹⁴ ages for the Holocene stratigraphy of Volcán de Colima from Luhr et al. (2010). Out of 181 dates collected, 143 are from the last 9000 years. Explosive activity at Colima appears to occur in distinct clusters (i.e. units V-Y and P-T) with gaps of 100 to 400 years. Unit Z spans ~ 4000 years of eruptive activity with numerous gaps of 100–300 years, including two sector collapse events at approximately 2500 and 3600 yrs. B.P. (Siebe et al., 1992; Navarro et al., 1994; Komorowski et al., 1997; Cortés et al., 2010)

Underlying unit Z is a series of three tephra fall deposits interbedded with ash-rich surge deposits (units U to Y; Fig. 2, Table 1). Based on the planar contacts between the units and the lack of palaeosols, Luhr et al. (2010) described these deposits as resulting from either a series of closely spaced compositionally zoned eruptions or a single event. These units range in composition from basaltic-andesite (54.2 wt.% SiO₂) in unit U to high silica andesite (60.8 wt.% SiO₂) in unit Y defining a clear magma differentiation trend from mafic to more evolved, with unit W intermediate in between the two end members (Fig. 4) (Luhr et al., 2010). Units V to Y erupted between ~ 4580 ± 50 yrs. BP and ~ 4460 ± 40 yrs. BP (Luhr et al., 2010), supporting the hypothesis that these were part of a single zoned eruption. Charcoal sampled from pyroclastic surge deposits underlying unit U gave an older age of ~ 4820 ± 50 yrs. BP (Luhr et al., 2010); however, the contacts between Unit U and the overlying unit V are planar suggesting very close succession of eruptions (Luhr et al., 2010).

Fig. 4

Whole-rock major element variation diagrams showing the compositional range of units P to Y. A = Andesite, B = Basalt, B-A = Basaltic-Andesite, TB-A = Tholeiitic Basalt-Andesite, TA = Tholeiitic Andesite. Units U, W and Y define a differentiation trend from basaltic-andesite (54.2 wt.% SiO₂) to high-Si andesite (60.8 wt.% SiO₂). Units P and S are intermediate between units U and Y

Based on detailed stratigraphic correlations using a combination of field mapping, petrology, geochemistry (Crummy et al., 2014) and the radiocarbon dates of Luhr et al. (2010), we have re-interpreted the stratigraphy below unit U, from approximately 5000 to 13,000 yrs. BP (Fig. 2; Table 1). Underlying unit U is a complex sequence of pyroclastic surge and fall deposits which Luhr et al. (2010) defined as units P to T, produced by a closely spaced series of sub-Plinian to Plinian eruptions akin to units U to Y. Through our interpretation of the stratigraphy, it has become clear that these units are more complex than previously thought, with numerous (> 10) surge and fall deposits erupted over an 800-year period between approximately 5200 and 6000 yrs. BP (Figs. 2 & 3).

Below unit U is unit T which comprises a pumice-rich fall deposit bound by ash-rich horizons. The contacts between these horizons are gradational, and all three are not always present. Below this sequence is a dark brown-orange basaltic-andesite (55.7–58.5 wt.% SiO₂) scoria fall deposit (unit S). Grain size variations within this unit suggest it is not a single fall deposit but a series of two or three pulsatory fall deposits interbedded with surge deposits. In some exposures there are clear variations in grain size with a coarse (8 cm scoria) layer underlain by a finer grained (up to 2 cm clasts) fall deposit (section VF95–09). Interbedded with these are distinct ash-rich lenses, which in some localities are identified as well-defined pyroclastic surges (section VF10–07). In other locations, unit S appears as a single fall deposit with little internal variation (section VF11–08). Dating of charcoal found within the binding pyroclastic surge horizons yield ages of 5430 ± 50 to 5500 ± 60 yrs. BP for unit S (Luhr et al., 2010; Table 1). Unit Q consists of two tephra fall deposits separated by an organic-rich ash layer (section VF95–09). In some localities, there is no ash layer preserved (section VF10–01); while in others only one fall deposit is observed (section VF97–03). Charcoal samples within this unit yield ages of 5870 ± 80 and 5940 ± 110 yrs. BP (Luhr et al., 2010). Below unit Q is a basaltic-andesite (56.1–59.4 wt.% SiO₂) pumice fall deposit varying in thickness from 6 to 42 cm across the northern CVC. Radiocarbon ages of 5980 ± 50 and 6150 ± 40 from binding ash horizons gives an approximate age of 6000 yrs. BP for unit P (Luhr et al., 2010). Underlying unit P is a distinctive massive ash layer varying in thickness from 30 to > 100 cm (Fig. 2) with pumice clasts and charcoal fragments giving it a distinct speckled appearance (unit O). Eleven radiocarbon dates from unit O give an age range from 6200 ± 40 to 6950 ± 50 yrs. BP (Luhr et al., 2010). Because of the similarity in appearance and the broad age range within this unit we have interpreted it as representing a series of pyroclastic surges and ash-rich fall deposits, similar to unit Z, albeit deposited over a shorter timescale.

Correlation of units below unit O has proved difficult due to limited exposure of older deposits and/or lack of dateable charcoal. However, we have identified three chemically distinct marker horizons (units N, F and D; Crummy et al., 2014), which we have used to stratigraphically correlate some of the interbedded units (Table 1, Fig. 2). Although numerous older tephra fall deposits have been observed, due to limited or no exposures and lack of appropriate charcoal samples for dating, we have been unable to correlate below unit D (12,460 ± 60 and 13,350 ± 130 yrs. BP) with any confidence.

Reconstruction of Holocene explosive eruptions

Eruption column height and volume

Estimated eruption column heights for units P to Y have large uncertainties with ranges from 10 to 30 km above sea level (a.s.l). Median values show little variation between units ranging from 22.5 km a.s.l. in unit W to 25 km a.s.l in unit S (Fig. 5a). The total erupted mass is better constrained as the thickness and grain size distribution of the deposit are directly proportional to erupted mass (Scollo et al., 2008). Unit Y has the highest erupted mass, with a median of 3.6 × 10¹¹ kg and a maximum of 7.7 × 10¹¹ kg (out of 100 simulations; see methods section), while unit U has the lowest with a median erupted mass of 1.3 × 10¹¹ kg and a maximum of 4.1 × 10¹¹ kg (Fig. 5b). These values equate to average erupted volumes of 0.36 and 0.13 km³ for units Y and U, respectively, assuming an average basaltic-andesite deposit density of 1000 kg/m³. The median total erupted masses for units W, S and P are 1.4 × 10¹¹, 1.8 × 10¹¹ and 1.7 × 10¹¹ kg, respectively. Volume estimates calculated using exponential thinning of the deposit after Pyle (1989) are equivalent to the minimum volume estimated by Tephra2 (Fig. 5b), but have high uncertainty due to the need to interpolate smooth isopach contour lines through sparse data.

Fig. 5

Boxplots showing the range (25th to 75th percentile), median, and upper and lower extremes of (a) the eruption column height (km a.s.l.) and (b) the total erupted mass (kg). Eruption column heights were difficult to constrain due to the proximity of the deposits to the vent, as different column heights can produce the observed deposit. Total erupted mass was easier to constrain as the thickness and grain size distribution of the deposit is dependent on the erupted mass. Diamonds show erupted mass estimates (converted from volume) using the Pyle (1989) method

Results from simulating a single eruption that deposited tephra fall units Y and W, gave an average eruption column height of 23 km a.s.l., with a maximum of 30 km a.s.l., and an average erupted mass of 8.3 × 10¹¹ kg with a maximum of 4.4 × 10¹² kg. Considering a typical arc basalt-andesite density of 1000 kg/m³, these equate to volumes of 0.83 km³ and 4.4 km³.

Tephra dispersion maps

Tephra dispersion maps based on our new field and geochemical cross correlations reveal that the main axis of dispersion is to the northeast for units P to Y (Fig. 6). This is consistent with modelling results of the 1913 deposit and eyewitness accounts and reports of ash fall from the 1913 eruption (Bretón González et al., 2002; Saucedo et al., 2010; Bonasia et al., 2011; Connor et al., 2019). The tephra dispersion maps for units Y, S and P reveal tephra dispersal to the north-northeast, with the town of Ciudad Guzman (23 km to the northeast of Volcán de Colima) experiencing ash fall thicknesses of 5 to 10 cm. However, unit S is not as widely dispersed as unit P or the larger in volume unit Y. The isopach map for unit W reveals tephra dispersion to the northeast and only 1 to 2 cm thick ash fall in Ciudad Guzman. Unit U has a similar erupted mass to units W and S of 1.48 × 10¹¹ kg (~ 0.1 km ³), but the dispersal map reveals a thicker ash fall deposit in Ciudad Guzman of 10 cm (Fig. 6). The tephra dispersion map for a single unit W-Y eruption reveals tephra accumulation in Ciudad Guzman of 10–20 cm (Fig. 6). The model input parameter ranges and best-fit values used to derive the tephra dispersion maps in Fig. 6 are summarised in Table 2.

Fig. 6

Tephra dispersion maps of the Holocene units Y to P generated using the parameters in Table 2. Dispersal is to the northeast over Ciudad Guzman, which has estimated ash thicknesses of 1–10 cm for individual units Y to P and 10–20 cm for the combined unit W-Y, making it susceptible to ash fall from potential future Plinian eruptions

Table 2 Model input parameter ranges and best-fit eruption parameters based on geological data for each of the Holocene units

We have calculated hazard curves at various locations along the axis of dispersion from the town of Tuxpan, 25 km east-northeast, to Guadalajara International Airport, 120 km north-northeast of Volcán de Colima for eruptions of Volcanic Explosivity Index (VEI; Newhall and Self, 1982) 3, 4, 5, and 6 (Fig. 7). According to the modelling results of Connor et al. (2019), using just the geological record, the 1913 explosion was a VEI 3 with an eruption mass range of 1 × 10¹⁰ to 1 × 10¹¹ kg. Units P to Y were VEI 4 events (1 × 10¹¹ to 1 × 10¹² kg), while the combined unit W-Y eruption was a VEI 5 event (1 × 10¹² to 1 × 10¹³ kg). The probability for tephra accumulation exceeding 100 kg/m² (threshold for roof collapse depending on roof type, among other factors; Jenkins et al., 2015) for a VEI 3 explosion in Ciudad Guzman is less than 1%. For a VEI 4 event, this increases to 6%. For VEI 5 and 6 explosions, the probability of tephra accumulation exceeding 100 kg/m² is approximately 30% and 65%, respectively. Assuming an average basaltic-andesite density of 1000 kg/m³, this is equivalent to a tephra thickness of approximately 10 cm. In Ciudad Guzman, there is a 90% probability of tephra accumulation from a VEI 3 explosion from Volcán de Colima (Fig. 7b), while at Guadalajara International Airport, 120 km north-northeast, there is a 20% probability (Fig. 7e).

Fig. 7

Hazard curves showing the probability of mass loading of tephra exceeding a specified amount for VEI 3, 4, 5 and 6 eruptions for various towns along the principal axis of dispersion: a) Tuxpan, 25 km east-northeast of Volcán de Colima; b) Ciudad Guzman, 23 km northeast; c) Sayula, 40 km north-northeast; d) Zacoalco de Torres, 79 km north-northeast; and e) Guadalajara International Airport, 120 km north-northeast of Volcán de Colima

Discussion

Completeness of the explosive record

The radiocarbon ages presented by Luhr et al. (2010) show almost continual explosive activity at Paleofuego de Colima throughout the Holocene. However, analysis of the data reveals distinct clusters of activity and gaps in the time line of eruptions of a few hundred years (Fig. 3). Unit O represents over 700 years of activity, deposited between 6200 ± 40 to 6950 ± 50 yrs. BP. Within this unit are multiple gaps of ~ 100 years. Overlying unit O are units P to T, a series of at least 10 pyroclastic surges and tephra fall deposits erupted over 800 years. Following the deposition of these units, is a period of ~ 520 years with only three radiocarbon ages. Units V-Y erupted during an approximately 120-year period as a single eruption following the deposition of unit U. Luhr et al. (2010) suggested that unit U deposited as part of the unit V-Y eruption series; however, there is a ~ 200-year gap between units U and V. Following the units U-Y series of eruptions, there is a gap of ~ 420 years. During this period there is a single charcoal date of 4240 ± 110 yrs. BP followed by a 280 year gap with no radiocarbon ages. Unit Z represents approximately 4000 years of activity at Paleofuego de Colima during which, the edifice suffered at least two large gravitational collapse events dated at approximately 2500 and 3600 yrs. BP (Siebe et al., 1992; Komorowski et al., 1997; Cortés et al., 2009; Cortés et al., 2010). Within unit Z there are no clear horizons relating to these events. We have been unable to correlate any of the gravitational collapse events and resultant debris avalanche deposits with the Holocene explosive stratigraphy exposed on the flanks of Nevado de Colima. The gravitational collapses may have destroyed evidence of previous eruptions, or they may have resulted in a change in behaviour of the volcanic activity.

Gaps in the stratigraphic record may represent a lack of preservation of explosive fall deposits. The geological record is biased towards larger magnitude eruptions, which are more likely to be preserved (e.g. Kiyosugi et al., 2015). Unit Z represents ~ 4000 years of ash accumulation with no pumice or scoria fall deposits. Luhr et al. (2010) related this to preservation, but this raises the question of why unit Z is so unique within the Volcán de Colima explosive stratigraphy. Unit O is a massive ash rich unit, similar in appearance to unit Z, but only represents ~ 700 years of activity. We suggest that unit Z may represent a change in behaviour of the volcanic activity that could be related to the gravitational collapse events that occurred ~ 2500 and ~ 3600 yrs. BP, within the time frame of unit Z. Further work should be concentrated on the unit Z ash deposits in order to fully investigate these potential preservation issues or changes in behaviour of the explosive activity at Volcán de Colima resulting from the sector collapse events.

Since the eruption of unit Y at approximately 4400 yrs. BP, there are no preserved tephra fall deposits until those of the 1818 and 1913 eruptions. However, based on historical observations we know that there have been at least nine sub-Plinian or Plinian events and numerous smaller scale Vulcanian explosions since 1519 (Bretón González et al., 2002). The 1818 and 1913 explosions are described as the most violent of the historical events and as such left traceable deposits (unit α). It may well be that there has been another switch in activity during historical times, with Volcán de Colima eruptions becoming more violent.

Wind direction at Colima displays a strong seasonal variation, with dominantly easterly winds during the summer months (June to September) and dominantly westerly and south-westerly winds in the winter months (November to March; Fig. 8). During April, May and October the winds are variable, related to the seasonal change in wind direction. Assuming similar wind conditions to present day, the tephra dispersion maps for units P to Y are consistent with eruptions during the winter months (November to March); however, it is unlikely that all Plinian eruptions at Colima occurred during November to March. Paleofuego de Colima has experienced numerous sector collapse events resulting in large debris avalanche deposits to the south. The youngest collapse has been dated at ~ 2500 yrs. BP (Siebe et al., 1992; Navarro et al., 1994; Komorowski et al., 1997; Cortés et al., 2010); therefore, it is possible that more explosive sub-Plinian and Plinian eruptions have occurred in the past that have been destroyed to the south and are not preserved in the tephra record to the north and northeast of Colima. During our field campaigns, we discovered a white pumice fallout deposit with clasts up to 15 cm, 12 km southeast of the volcano. Well preserved and abundant charcoal from this deposit was dated at ~ 7600 yrs. BP (Savov, unpublished) but we can find no correlation with the rest of the Holocene tephra deposits exposed to the north and northeast of Colima. This and other tephra deposits from highly explosive eruptions may have occurred during the spring, summer and autumn months, depositing tephra to the southeast, southwest and west. Future work should be concentrated in these areas with the aim to establish (geographic) links between tephra deposits.

Fig. 8

Wind rose diagrams for Colima from 2015. Wind data were downloaded from NOAA Reanalysis2 (NOAA, 2015). The coloured bars indicate the frequency of counts (%) in the direction the wind is blowing, and the colours correspond to wind speeds. There is strong seasonal variation in wind direction between the summer months and winter months, with dominantly westerly and south-westerly winds in November to April, and dominantly easterly winds in June to September. May and October are variable reflecting the seasonal change in wind direction

Sources of uncertainty around modelling tephra fall deposits

As with all modelling, there are uncertainties around the results related to the deposits themselves and the numerical modelling, which we have tried to address. The tephra fall deposits exposed on the flanks of Nevado de Colima represent only the coarse fractions (> 250 μm, with median phi of 0 to − 1.6), which is consistent with their relatively proximal location in respect to the vent (4 to 12 km; Fig. 2). Our modelling results could therefore be considered minimum estimates, as they may not take into account the finer fractions, which are dispersed away from the volcano (e.g. Bonadonna and Houghton, 2005; Cashman and Rust, 2016; Costa et al., 2016). Rose and Durant (2009) estimate that tephra generated during silicic Plinian eruptions that produce abundant PDCs contain large proportions (30 to > 50%) of very fine ash (< 30 μm), for example the 1980 Mount St Helens eruption. The tephra deposits preserved on the flanks of Nevado de Colima could therefore represent < 50% of the total erupted deposit. In contrast, it may be that the deposits resulted from coarse-grained eruptions, in which case the total grain size distribution estimated by the inversion modelling are true reflections of the eruption and the modelled erupted masses are not significantly underestimated.

The five individual Holocene tephra fallout deposits modelled here, units P to Y, have minimum estimated erupted masses varying from 1.3 × 10¹¹ kg to 3.6 × 10¹¹ kg (0.13 to 0.36 km³). Connor et al. (2019) modelled the historical 1913 eruption deposit using the same method applied here. Their results yielded a minimum erupted mass of 4–7 × 10¹⁰ kg. However, Saucedo et al. (2010) and Bonasia et al. (2011) modelled a much larger erupted mass of 2–6 × 10¹¹ kg for the 1913 eruption based on historical accounts of ash fall and field measurements. Connor et al. (2019) carried out a comparison between modelling using only field measurements and using field measurements combined with historical observations. Their results revealed an order of magnitude difference, with the modelling incorporating historical observations yielding a larger erupted mass of 1–5 × 10¹¹ kg, comparative to the results of Saucedo et al. (2010) and Bonasia et al. (2011). If this relationship holds for older units, it suggests that our modelling of the proximal tephra fallout deposits could underestimate the erupted mass by an order of magnitude suggesting that the pre-historical Holocene eruptions were much larger magnitude than we initially thought.

Conversely, these units could represent a number of distinct, but individually small eruptions, built up over a period of potentially hundreds of years. Unit α, as described by Luhr et al. (2010), can appear as a single fall deposit in some localities when in fact it is a combination of the 1818 and 1913 eruptions. Luhr et al. (2010) distinguished between these eruptions based on geochemistry. Through detailed whole-rock and mineral geochemical and petrological analysis (Crummy et al., 2014), we have not seen any such distinctions within the units, therefore we believe these deposits represent single explosive events.

Another source of uncertainty concerns the relationship between the tephra layers preserved in the stratigraphy compared to the original fallout (Engwell et al., 2013). Unconsolidated tephra layers are highly susceptible to surface erosion by wind and rain (Collins et al., 1983). Slope angle and vegetation cover can also have a significant impact on tephra preservation (Cutler et al., 2016; Dugmore et al., 2018). We have measured the tephra unit thicknesses at up to 30 locations across an area of 130 km², reducing the uncertainty related to variability in the preserved deposit thickness and uncertainty related to field measurements (Engwell et al., 2013).

There are also sources of uncertainty related to the modelling. Many combinations of input parameters can match the observed deposit thickness and grain size distribution (Connor and Connor, 2006). The model therefore uses the downhill simplex algorithm (Nedler and Meade, 1965; Connor and Connor, 2006), a non-linear inversion model, to search for the best-fit eruption parameter set that minimises the error between the measured tephra thickness and grain size distribution, and the model outputs (Connor and Connor, 2006). Simulations are run in parallel on multiple processors which allows rapid implementation of the physical model and a fully probabilistic analysis of tephra fall hazard (Bonadonna et al., 2005). One such example of the modelling uncertainty is with the eruption column height estimates, which show a wide variation in values for all of the units. The proximity and clustering of the samples close to the vent means that the grain size distribution reflects only the coarse fractions, and does not reflect the very fine ash particles that are carried tens or hundreds of kilometres from the vent. As a result, there is no unique solution for eruption column height for the measured data: a variety of different column heights may reproduce the eruption deposit based on currently available data. This uncertainty is consistent with results from sensitivity analysis of Tephra2 carried out by Scollo et al. (2008), and other studies (e.g. Cobeñas et al., 2012).

Magnitudes of pre-historical Holocene explosive events

The tephra dispersal modelling results reveal minimum erupted volumes of 0.12 to 4.4 km³ for the five best exposed and cross-correlated pre-Historical Holocene explosive eruptions. The principal axis of dispersion was to the northeast and north-northeast over the town of Ciudad Guzman, with minimum ash fall thicknesses of up to 10 cm from individual unit P to Y explosions in Ciudad Guzman, and 10–20 cm for the combined unit W-Y deposits assuming deposition during a single eruption.

If our results for the erupted mass, which are modelled based only on the geological record, are underestimated by an order of magnitude, then consequently ash fall thicknesses could have been much greater. Hazard curves generated for various locations along the axis of dispersion reveal that the tephra accumulation can vary dramatically with eruption magnitude. By incorporating the historical observations of the 1913 eruption into their model, Connor et al. (2019) increased the estimated volume of the 1913 eruption by an order of magnitude from VEI 3 to VEI 4. If this relationship holds for the pre-Historical Holocene units modelled here, these events could have been VEI 5 events, which has significant implications for tephra accumulation. In Ciudad Guzman, this increases the probability from approximately 6% to 30% for 100 kg/m² tephra load (10 cm thickness; Fig. 7b). This demonstrates that the tephra fall hazard assessment based on a reasonable worst-case scenario can underestimate the maximum hazard magnitude if only the geological record is considered. We recommend that probabilistic hazard assessments for tephra fallout consider a broad range of explosive eruption scenarios, including eruptions smaller and larger than those identified in the geologic record and modelled. Preservation biases the geologic record towards larger eruptions. Our results suggest that lack of preservation of distal deposits in the geologic record may bias models of eruptive volumes to smaller values.

Conclusions and implications for disaster risk management

Comprehensive stratigraphic correlations identified through field mapping, radiocarbon dating and petrological and geochemical fingerprinting by Luhr et al. (2010), Crummy et al. (2014) and in this study have resulted in possibly one of the best stratigraphic records of explosive eruptions at an arc stratovolcano throughout the Holocene and into the Late Pleistocene.

Inverse modelling based on geological data for five of the best exposed and correlated pre-Historical Holocene tephra fall deposits reveal minimum erupted volumes of 0.12 to 4.4 km³. Results from a recent study by Connor et al. (2019) on the 1913 event reveal that simulations run based on only geological data can underestimate the erupted volume by at least an order of magnitude. Simulations were run based on thickness and grain size distribution data collected from proximal locations, therefore representing only the coarse fractions. Consequently, our modelling results may be underestimated. By modelling tephra dispersion from events in the stratigraphic record, this work shows that the maximum magnitude event at Colima could be a VEI 5. By using the hazard curves for this newly estimated maximum magnitude event, rather than the maximum magnitude event recorded by either the historical or geological record alone, it is apparent that the expected tephra accumulation for a worst-case scenario event increases by at least an order of magnitude.

This demonstrates that hazard assessments based on historical activity or the geological record alone can underestimate the maximum magnitude of potential future eruptions. Knowledge of the volcano from historical records, and detailed tephra-stratigraphic work and modelling often sparse data are necessary to address these potential biases. However, we have demonstrated that even with detailed stratigraphy and inverse modelling presented here, our understanding of the full potential range of explosive events remains limited.

Volcán de Colima is a frequently erupting volcano with many communities living nearby. Over 1.7 million people live within a 100 km radius of Volcán de Colima, and approximately 500,000 people live within a 30 km radius (INEGI, 2010). Large magnitude eruptions would have significant impacts on the local and regional population, infrastructure and economy. Ash fall thicknesses of a few millimetres to tens of centimetres can result in transport problems and damage to critical infrastructure (e.g. Wilson et al., 2012; Jenkins et al., 2015; Hayes et al., 2017). The potential impacts from a future Plinian eruption at Volcán de Colima are therefore significant on a local scale, but could, with the Manzanillo-Guadalajara highway passing within 20 km to the east of the volcano, also have considerable impacts on regional and national scale. Ash from the 1913 Plinian eruption reportedly fell in Guadalajara, which now hosts one of Mexico’s busiest international airports. A large magnitude explosive event at Volcán de Colima could therefore have far-reaching impacts.