Introduction

Eruptive cycles, evolving from dominantly effusive, dome-forming eruptions to explosive Plinian eruptions, are typical of intermediate to silicic magma systems. These cycles have different timescales and can occur between several eruptions or within single eruptive events (Cassidy et al. 2018). Eruptive style during these cycles involves persistent effusive-explosive transitions and is thought to be controlled by both pre-eruptive magma conditions (crystallinity, pressure, temperature and volatile content) and storage and conduit processes, including magma recharge and mixing, ascent, degassing, crystallization and rheological changes, with complex feedbacks between each (e.g. Villemant and Boudon 1999; Andujar and Scaillet, 2012; Martel 2012; Cassidy et al. 2018; Popa et al. 2021). Establishing links between eruptive style and the evidence for magma properties that is recorded in erupted products provides constrains on the mechanistic modelling of eruptive style and related hazards at active volcanoes.

Eruptive cycles during the late Holocene at Guagua Pichincha volcano (Ecuador) are well-documented by tephrostratigraphy and 14C dating (Barberi et al. 1992; Robin et al. 2008). Recurrent activity has been characterised by cycles of phreatic explosions and dome-forming eruptions (including Vulcanian phases), which culminate in Plinian activity (volcanic explosivity index VEI 4–5). The activity is characterised by a cycle duration of 100–200 years separated by quiescent periods of several centuries (Robin et al. 2008).

Here, we compare the two most recent magmatic eruptions of Guagua Pichincha volcano: the CE 1660 Plinian eruption and the dome-forming and Vulcanian-Subplinian activity of 1999–2001, in order to highlight potential variations in conduit processes (degassing and crystallization) that may in turn control the evolution of eruptive style at different stages of the cyclic behaviour. Previous studies of Guagua Pichincha dacites have provided a petrological analysis of the erupted products to yield estimates of eruptive temperatures and crystallization pressures (Samaniego et al. 2010). Petrophysical studies have thus far focused on the 1999–2001 period and largely on breadcrust bombs from Vulcanian phases (Wright et al. 2007). These studies do shed light on storage conditions as well as crystallization and degassing processes in the Guagua Pichincha dacitic magma. However, a comparison between products of Plinian and dome-forming/Vulcanian activity is needed because Guagua Pichincha may potentially be in the middle of an eruptive cycle evolving towards more explosive behaviour (Tadini et al. 2021).

Our study combines chemical, petrophysical and textural analyses of Guagua Pichincha dacites. We link pore properties observed in natural products to percolation models in order to improve our understanding of degassing and gas escape in the dacitic magma. We also explore the role of magma crystallinity and bulk viscosity on magma ascent rate, degassing processes and, ultimately, on eruptive style. Finally, we generalize our findings to provide a conceptual model that relates conduit processes to cyclic eruptive activity at Guagua Pichincha.

Geological setting

Guagua Pichincha is a 10-km-wide stratovolcano located on the western part of the older, composite and extinct Rucu Pichincha volcano (Robin et al. 2010). It is situated near the city of Quito, Ecuador’s capital city that has a population of ca. 3 million inhabitants (Fig. 1a). Guimarães et al. (2021) describe it as one of the most hazardous and risk-prone volcanoes in Latin America. Active monitoring is therefore continuously carried out at Guagua Pichincha by several geophysical techniques (Ramón et al. 2021). Recent volcanic activity occurred in the Cristal dome volcanic complex, which was emplaced in the Toaza amphitheatre after a major landslide 4 ka years ago (Fig. 1b; Robin et al. 2010). The Toaza edifice was emplaced in a former amphitheatre formed about 9–11 ka ago. Previous studies have recognized three major eruptive cycles over the last 2000 years, with major Plinian episodes in the first century, tenth century and in CE 1660 exhibiting VEIs of 4–5 (Robin et al. 2008, 2010). Between these major Plinian eruptions, several dome-forming, Vulcanian and Subplinian eruptions as well as episodes of phreatic eruptions and non-eruptive unrest have occurred. As such, Barberi et al. (1992) proposed that the explosive eruptions of the tenth century and CE 1660 represent the maximum eruptive events (VEI 4–5) to be expected in future eruptive episodes. Recent 1981–1998 phreatic activity and subsequent dome-forming eruptions in 1999–2001 have been interpreted as a possible initial step towards more violent explosive activity, following previous cycles of increasing explosivity towards Plinian eruption (Robin et al. 2008). For example, the last Plinian eruption of CE 1660 was preceded by what were possibly similar minor dome-forming eruptions in CE 1566, 1575 and 1582, as well as phreatic activity (Wolf, 1904).

Fig. 1
figure 1

Geological setting of Guagua Pichincha volcano and location of samples studied here. a Map of Guagua Pichincha near the city of Quito. The center of this image has a latitude and longitude of 0.203889°S and 78.511111°W, respectively. Inset on upper right corner is a Google Earth image (Imagery date: 12/14/15-newer. Google Data SIO, NOAA, U.S. Navy, NGA, GEBCO Landsat/Copernicus INEGI IBCAO). b to e correspond to Google Earth images showing the location of the main samples analysed in our study. b Overview of the outcrops where samples were collected in November 1999 and January 2020, where letters c to f correspond to the emplacement of outcrops shown in the corresponding images. c Detail of the Cristal dome complex showing the 1981 phreatic crater and the location of the current dome emplaced in late 2000 (with sampling location indicated by a red square) and remnants of the pre-existing 1660 dome. d Location of samples from fallout deposits of the three Plinian eruptions (first century, tenth century and CE 1660) near El Refugio. e Samples collected from block and ash flow deposits of the 1999 eruption westward from the vent in the Rio Cristal valley. f Pumice flow deposit of the CE 1660 Plinian eruption. North is up in (a) to (e)

The CE 1660 eruption occurred primarily in October–November 1660 with its climax on the 28 October, resulting in pumice falls and pyroclastic currents (Wolf, 1904; Barberi et al. 1992). This main explosive phase was followed by repeated ash emissions and dome emplacement in November 1660 (Barberi et al. 1992). The greatest volume of pyroclastic currents was directed westward towards uninhabited areas, whereas on the eastern side they were stopped by the topographic relief of Rucu Pichincha and flowed down the Lloa and Mindo Valleys, thereby protecting Quito (Wolf, 1904; Barberi et al. 1992; Canuti et al. 2002).

The 1999–2001 eruption was preceded by periods of phreatic activity in the nineteenth century and in the two decades preceding the eruption (1981–1998; Robin et al. 2008). Minor ash-rich phreatic explosions occurred in June 1999, with the first major magmatic eruption on the 26 September (Wright et al. 2022). Repeated episodes of dome formation and destruction occurred from late-September through mid-December 1999 and two main explosive eruptions on the 5 and 7 October, together with several additional Vulcanian eruptions (Wright et al. 2007). Dome emplacement continued in late-December followed by destruction in June 2000 (Global Volcanism Program, 2000). The current dome was emplaced in mid- to late-2000 (Mothes et al. 2002). Ash fallout threatened Quito again during May 2001 (Robin et al. 2008).

The 5 and 7 October 1999 eruptions were low intensity, pumice-rich eruptions that differed from the 1999 Vulcanian activity in their lack of dense rocks and longer durations. We refer to these eruptions as Subplinian eruptions based on their plume height of 5–10 km and dominance of pumice clasts in the deposits, although the activity was less sustained and explosivity was lower (VEI = 2) than for typical Subplinian activity.

A recent study on eruption type probability indicated that the most probable eruption type at Guagua Pichincha in the next 100 years is Vulcanian, but also proposed that Subplinian and Plinian eruptions are highly likely (Tadini et al. 2021).

Samples and methods

Bulk samples of dacites from the CE 1660 and 1999 eruptions of Guagua Pichincha were primarily collected during January 2020. Samples from the CE 1660 eruption consist of pumiceous, banded and dense dacites and originate from fallout deposits along the Lloa-Refugio road close to El Refugio (Fig. 1d) and from a pumice-flow deposit near the Rio Cristal (Fig. 1f). Samples from the 1999–2001 eruption consist of dense blocks from block and ash flow deposits collected at Rio Cristal (Fig. 1e), a dense rock sampled from the dome carapace (Fig. 1c) and pumice clasts from the 5 and 7 October Subplinian eruptions collected around the Rio Cristal valley westward of the vent. In addition, some block and ash flow samples that were analysed for density and tortuosity were collected from hot deposits ~ 9 km from the dome, near the Rio Cristal, in November 1999.

We first performed a simple visual componentry analysis of the tephra in the laboratory from the CE 1660 Plinian eruption and compared with the other Plinian eruptions of the first and tenth century (collected at the same outcrop along the Lloa-Refugio road; Fig. 1d). The tephra was divided into the fractions 2.0–2.8, 2.8–4.0, 4.0–5.6 and 5.6–8.0 mm, and clasts were separated into vesicular, dense and accidental classes based on colour, texture and shape criteria (Fig. SM1 in Online Resource 1).

Chemical analysis of the dacites

We analysed groundmass glass chemistry of a pumice clast from the CE 1660 eruption and of a dense rock from a block and ash flow of the 1999 eruption. We used a Cameca SX-100 microprobe at the Ludwig Maximilian University (LMU, Germany) with a 5-μm defocused beam at an accelerating voltage of 15 keV and a current of 5 nA in order to limit Na loss during glass measurements. For the dense rock, areas with lowest local microlite content were selected for analysis. These data were then combined with additional data measured by Wright et al. (2007) and Samaniego et al. (2010).

Petro-physical measurements

We measured the bulk density and total porosity of 149 rocks from the CE 1660 and 1999 eruptions using several techniques: geometrical measurements on cylinders, the immersion technique in water, the three-axes method and two-dimensional image analysis (Houghton and Wilson, 1989; Shea et al. 2010; Barker et al. 2012; Bernard et al. 2015; Colombier et al. 2021; Pisello, 2017). In addition, the bulk density and total porosity of 51 rocks from the first and tenth century Plinian eruptions were measured. The solid density of two pumice clasts from the CE 1660 eruption, one block and ash flow dense rock and several (> 15 clasts) pumice clasts from the 1999 eruption were obtained by helium pycnometry on powdered samples. The total porosity (\(\Phi\)) was then calculated from the bulk (\(\rho\)) and solid (\(\rho\) s) densities via the formula \(\Phi =1-\left(\frac{\rho }{\rho s}\right)\) and by image analysis in 2D (scanning electron microscope, SEM). A full description of the methodologies used for analysis of petro-physical and textural character of the samples is given in Online Resource 1.

Pore connectivity of 34 samples was measured by He-pycnometry with a Quantachrome Ultrapyc 1200e at LMU, yielding the sample specific volume including the solid phases and isolated vesicles only, and allowing us to subsequently calculate the connected porosity (Formenti and Druitt, 2003; Colombier et al. 2017). Connectivity was then retrieved by dividing this connected pore fraction by the total porosity. We also performed a floating experiment (Online Resource 1) on pumice clasts of similar size and porosities from the Plinian and Subplinian activity, in order to provide complementary information on isolated porosity and the size of vesicle throat openings (Mitchell et al. 2021).

Permeability was measured on 11 cylinders with diameters of about 20 mm and lengths of 30–50 mm, cored from the pumice lapilli (> 32 mm in diameter) and bombs, as well as on two block and ash flow samples and one dome rock (which were 10 mm in diameter and nominally 25 mm in length). All samples were dried in a vacuum oven at 40 °C for at least 48 h prior to analysis. The permeability of each sample was measured using a benchtop gas (nitrogen) permeameter at the Strasbourg Institute of Earth and Environment (see Farquharson et al. 2016; Heap and Kennedy, 2016). Permeability measurements were performed at ambient temperature and under a confining pressure of 1 MPa using the steady-state method. Samples were first left for one hour at the target confining pressure to ensure microstructural equilibration. Steady-state volumetric flow rates were then measured (using a gas flowmeter) for six different pore pressure differentials (measured using a pressure transducer; the maximum pressure differential used was 0.2 MPa). These data were used to calculate permeability using Darcy’s law and to check whether the Klinkenberg or Forchheimer corrections were needed, which were applied on a case-by-case basis (the methods of data reduction are described in Heap et al. 2017).

In addition, we obtained the electrical tortuosity of dense to vesicular blocks collected in late 1999 from valley-confined pyroclastic density current deposits, using electrical conductivity measurements on 17 cylindrical samples (25 mm in diameter and 20–40 mm in length) directly cored in the laboratory. Because the water-filled fractional porosity controls the path of the electrical flow, the cores were wrapped with Teflon ribbon to avoid core desaturation and saturated with NaCl brines for the electrical measurements (details in Jouniaux et al. 2000). The conductivity was measured on a two-electrode assemblage (Wayne-Kerr bridge), once equilibration with the saturating solution was achieved, and a frequency of 1 kHz was selected for the full salinity range of the study (0.02 to 30 g/L). Six to seven measurements were obtained with saturating solutions of increasing NaCl concentration. The measured resistivity of the cores and of the saturating NaCl solutions yields the electrical surface conductivity, Cs (following the model of Revil and Glover, 1998), which is in the range of 0.04 to 0.26 mS/m in our sample set.

Scanning electron microscope textural observations and quantitative analysis

Back-scattered electron (BSE) images were taken using a HITACHI SU 5000 Schottky FE-SEM at LMU. Phenocryst and micro-phenocryst crystallinities (Online Resource 2) were obtained via a combination of BSE and energy-dispersive X-ray spectroscopy (EDX) mapping of the thin sections, and were then corrected using 2D porosity (i.e. the crystallinities reported here represent the fraction of crystals within the groundmass). As this method is relatively time-consuming, it was applied on selected samples, namely one pumice lapillus from the CE 1660 eruption, and one dense rock and three pumice clasts from the 1999 eruption. We also determined the microlite content (on a vesicle-free basis) of the groundmass based on BSE images. The limit between microlites and micro-phenocrysts/phenocrysts was set at 30 µm, consistent with the observed textures such as change in crystal habit and lack of compositional zonation for crystals below this size. Total crystallinity (Φx-t) was calculated as the sum of micro-phenocryst/phenocryst content and groundmass microlite content over selected areas of each sample, using the equation:

$${\Phi }_{\mathrm{x}-\mathrm{t}}={\Phi }_{\mathrm{x}-\mathrm{p}}+{\left(\left({1-\Phi }_{\mathrm{x}-\mathrm{p}}\right)\right)}^{*}{\Phi }_{\mathrm{x}-\mathrm{m}}$$
(1)

where Φx-p is the phenocryst content (calculated on a vesicle-free basis), and Φx-m is the microlite content of the groundmass (excluding phenocrysts and vesicles). All crystallinities are expressed here as fractions.

Porosity was also measured in 2D on smaller lapilli by binarizing BSE images on lapilli of variable sizes (2.0–2.8 mm fraction). We measured the vesicle size distribution (VSD) and vesicle number density (Nv) of one pumice clast from the CE 1660 eruption and two pumice clasts from the 1999 eruption. From the VSD, we could obtain parameters such as polydispersivity and average vesicle diameter to use as input parameters in the percolation models (Online Resource 1). Vesicles were imaged at different magnifications in the SEM to picture the full spectrum of sizes. Only one image per magnification was used to optimize time-consuming, user dependent segmentation and separation steps. This strategy did not alter distributions or parameters of the VSD (Fig. SM5 in Online Resource 1). The measured VSDs and Nv were corrected for crystal content where observable in the BSE images. Vesicles were also decoalesced following Shea et al. (2010) where glass films could be identified or inferred, allowing us to approximate values of Nv that represent the state of the magma prior to the last stage of coalescence. Using three-dimensional numerical bubble networks of Vasseur et al. (2020), we show that our 2D procedure involving decoalescence into individual vesicles and stereological conversion yielded accurate results for Nv quantification, whereas 3D methods underestimate Nv by several orders of magnitude for highly connected pore networks such as in the pumice clasts (see Online Resource 1).

Results

Here, we combine clast descriptions, microscopic textures, chemistry and petrology to provide a more complete picture of eruption dynamics, with particular focus on eruptions of 1660 and 1999–2001.

Macroscopic description

We define three main groups of dacite lithologies from Guagua Pichincha based on colour, shape, density, porosity and textural properties: (i) vesicular, (ii) banded and (iii) dense. Breadcrust bombs contain hybrid characteristics ranging from vesicular to dense depending on the location in the pyroclast (rim vs. core; Wright et al. 2007). Here, we classify both pumice clasts and vesicular breadcrust bomb interiors (thickly breadcrusted bombs TkB and thinly breadcrusted bombs TnB; Wright et al. 2007) as vesicular clasts. Clasts are defined as pumice when they are macroscopically homogenous and have a density lower than 1000 kg/m3. Only two clasts, with a density slightly higher than 1000 kg/m3, were also classified as pumice due to their textural similarities. Pumice clasts have varying morphologies and colour depending on the type of deposit (fallout vs. flow) and sampling location. Pumice clasts from the CE 1660 pumice flow deposit are white-beige in colour and have smooth, subrounded to rounded shapes due to abrasion and rounding of clasts during transport (e.g. Manga et al. 2011; Buckland et al. 2018; Hornby et al. 2020). Pumice clasts from the same eruption but sampled at the fallout deposit around El Refugio are subangular and have an orange-beige colour, likely as a result of alteration. Pumice clasts from the 1999 Suplinian activity are white and subangular.

Dense rocks are texturally homogenous with density values typically higher than 1800 kg/m3. Dense textures are also found in some breadcrust bomb rinds but differ from those of the dense rocks by the dominance of small, spherical isolated vesicles (Wright et al. 2007). The origin of the dense rocks observed in deposits from the CE 1660 eruption is uncertain as will be shown in more detail in the discussion. Some of these rocks likely represent partly recycled or remobilized dense particles from a previous eruptive phase (i.e. accidental lithics). Other dense clasts may have a juvenile origin with macroscopic properties that gradually transition to those of the pumice clasts. These dense rocks are grey in colour and angular to subangular in shape.

Dense rocks from the 1999–2001 eruption studied here are in contrast all juvenile in nature. These are mainly dense pyroclasts from block and ash flow deposits of the 1999–2001 eruption or rocks from the in situ lava dome in the Cristal dome complex. Dense rocks from this eruption also include the dense blocks (dense blocs DB) studied by Wright et al. (2007) that are grey in colour and angular. We also note the presence of intermediate vesiculated products in the samples collected in the block and ash flows from 1999. These samples are not classified as dense and are referred to below as block and ash flow intermediate type. Banded clasts have properties (texture, colour and shape) intermediate between vesicular clasts and dense rocks. This category includes FIB (foliated intermediate bombs) bombs described by Wright et al. (2007).

Chemistry and petrology of the Guagua Pichincha dacites

All rocks studied here are dacitic in composition and contain abundant phenocrysts and micro-phenocrysts, rhyolitic groundmass glass (± microlites) (Fig. 2; Online Resource 2), and varying vesicularities (Fig. 3). Differences are observed in whole-rock composition between the products of the CE 1660 and 1999–2001 eruptions, with the latter extending to higher SiO2 and Na2O and lower TiO2 and FeO contents (Samaniego et al. 2010) (Fig. 2), with some compositional overlap between eruptive suites. Minor differences also exist between other major and trace elements (not shown). Phenocrysts (and micro-phenocrysts) are dominantly plagioclase, amphibole, orthopyroxene, Fe-Ti oxides, minor clinopyroxene and/or olivine and scarce apatite and biotite crystals (Samaniego et al. 2010). Plagioclase and amphibole crystals are zoned in all dacites studied here, with oscillatory zoning dominating in plagioclase phenocrysts and micro-phenocrysts from pumice clasts and dense rocks, as has also been reported for breadcrust bombs (Wright et al. 2022).

Fig. 2
figure 2

Composition and phenocryst content of the Guagua Pichincha dacites. a to c Whole rock (WR) composition of the CE 1660 and 1999 products showing the evolution of TiO2, Fe2O3 and Na2O with SiO2 (Samaniego et al. 2010; Wright et al. 2007). d and e Glass composition of pumice clasts of the CE 1660 eruption and breadcrust bombs and dense rocks from the 1999 dome-forming eruption (Samaniego et al. 2010; Wright et al. 2022; this study). f The evolution of phenocryst-content (calculated vesicle-free) with porosity for the different lithologies studied here. Data were obtained by a combination of Back-scattered electron (BSE) BSE and energy-dispersive X-ray spectroscopy (EDX) mapping of thin sections (this study) and point counting (Wright et al. 2022). BCB corresponds to breadcrust bombs

Fig. 3
figure 3

Porosity and pore connectivity of pumice clasts, breadcrust bombs (BCB), dense rocks and banded clasts. a Frequency histogram showing the porosity distribution for the different types of rocks analysed here. b Pore connectivity as a function of porosity with corresponding error bars (errors are variable due to different techniques used)

Petrophysical analysis

All raw data of density, porosity, connectivity, permeability and tortuosity discussed hereafter are available in Table SM2 in Online Resource 3.

Differences in density and porosity of clasts of the two eruptions are negligible for a given lithology, and require a more extensive dataset to distinguish them statistically. Therefore, we treat those properties as equivalent for the two eruptions and highlight here differences between the lithologies generated in the two eruptions. We further include comparison with breadcrusted and dense ballistic bombs in Wright et al. (2007). Densities of pumice clasts range from 360 to 1210 kg/m3. Using the range of solid densities obtained (ρs = 2620–2670 kg/m3), this corresponds to total porosities (Φ) in the range of 0.54–0.86 (Fig. 3a). Smaller vesicular lapilli (2.0–2.8 mm fraction) have similar microtextures but substantially lower porosities compared with larger pumice clasts (0.27–0.67; Fig. SM2 in Online Resource 1), possibly due to the lack of large vesicles at small grain sizes. Dense rocks have densities in the range 1640–2630 kg/m3, corresponding to porosities of 0.01–0.38 (Fig. 3a). Banded clasts have intermediate density (1380–1800 kg/m3) and porosity values (0.32–0.48).

Pore connectivity (C) and permeability (k), plotted as functions of porosity, show distinct patterns for each lithology (Figs. 3b and 4). We observe two main patterns, one for dense rocks and one for vesicular clasts. Dense rocks (dense juveniles of the 1660 eruption and dense rocks from the 1999–2001 dome, block and ash flows and Vulcanian activity of the 1999 eruption) show a diffuse pattern with a broad range of connectivity within a low porosity window, but with only one sample at C < 0.9 (light and dark blue squares in Fig. 3b). Pumice clasts from the CE 1660 Plinian eruption show C values from 0.86 to 0.94 for a porosity range of 0.68–0.81 (orange circles in Fig. 3b). In comparison, pumice from the 1999 Subplinian eruptions show substantially lower connectivity values (0.78–0.85) for a similar porosity range (0.76–0.85) (yellow diamonds; Fig. 3b). Vesicular breadcrust bomb interiors show an increase of C (from 0.74 to 0.89) with Φ (from 0.47 to 0.71) (green crosses in Fig. 3b). Banded clasts plot at intermediate values of C and Φ (0.81–0.94 and 0.32–042) compared to dense rocks and vesicular clasts (pink triangles in Fig. 3b). The differences in micro-textures and connectivity indicated from helium pycnometry between pumice clasts from the 1999 and CE 1660 eruptions have been confirmed by a floating experiment (see Online Resource 1).

Fig. 4
figure 4

Permeability and related pore properties of pumice clasts, breadcrust bombs (BCB), dense rocks, banded clasts and intermediate vesiculated products from block and ash flows (BAF intermediate). a Permeability–porosity relationships in the dacitic rocks. Letters correspond to the Back-scattered electron (BSE) images (e) to (j). The grey area corresponds to experimental data of deGraffenried et al. 2019. b Permeability evolution with connected porosity corrected from crystal content. c Permeability-tortuosity relationship. d Permeability as a function of connectivity. The dashed line is used to guide the eye towards the separation between dense rocks and vesicular clasts. e to j BSE images showing the micro-textural variety of the samples analysed for permeability. In (j), the BSE image shows a dense band consisting of sintered ash particles (labelled “s”) in between pumiceous lapilli clasts (labelled as “lap”) in an agglomerated pumice

As with connectivity, we see increases in permeability as a function of increasing porosity for vesicular, intermediate and dense textural types (Fig. 4a). In vesicular clasts, k increases from 10−14 to 10−11 m2 as porosity increases from 0.47 to 0.77 (green crosses in Fig. 4a). The pumice samples from the CE 1660 eruption all have a porosity of ~ 0.75, but a permeability range that varies by about an order of magnitude, from ~ 5 × 10−14 to ~ 5 × 10−13 m2 (orange circles in Fig. 4a). Dense rocks increase in permeability from 10−17 to 10−12 m2 over a lower porosity window (0.01–0.3) (light and dark blue squares in Fig. 4a). Again, the banded clasts form an intermediate pattern between the vesicular clasts and dense rocks (pink triangles in Fig. 4a).

The electrical tortuosity data on the dense rocks are compared to data for breadcrust bombs (and one banded clast) obtained by Wright et al. (2009) (Fig. 4c). Tortuosity is lower (6–14) for the dense rocks than for the breadcrust bombs (14–53), and no clear trend is observed for the dense rocks (blue squares; Fig. 4c). In contrast, breadcrust bombs show a clear decrease of tortuosity with increasing porosity (green crosses; Fig. 4c). One banded clast and the block and ash flow rocks with an intermediate porosity plot at intermediate tortuosity (8–20) (pink triangle and red squares in Fig. 4c).

2D Micro-textural analysis

Based on microtextures, we subdivide the above textural categories into (Fig. 5): (i-1) pumice clasts, (i-2) vesicular breadcrust bombs (type TnB and TkB; Wright et al. 2007) characterised by rim-to-core textural variations due to expansion and chilling; (ii) banded clasts with intermediate porosities, including FIB bombs with pervasive flow-banding (Wright et al. 2007); (iii-1) type 1 dense clasts with absence of diktytaxitic porosity (void space in fully crystalline groundmass), rare silica-phase crystals and important shearing in vesicle clusters; (iii-2) type 2 dense with abundant diktytaxitic porosity largely filled by a silica phase (presumably in part cristobalite; Kushnir et al. 2016); and (iii-3) type 3 dense with bands of diktytaxitic porosity (with ubiquitous silica phase) similar to type 2 that alternate with denser bands dominated by cracks. Banding in type 3 dense rocks consists of anastomosing shear bands with crack-like pore textures and broken crystals similar to those in the shear partitioning experiments of Laumonier et al. (2011), and therefore these dense rocks were not classified as belonging to the banded population that shows both pumice-like and dense-like textural components.

Fig. 5
figure 5

Back-scattered electron (BSE) images at different magnifications showing the different categories of dacites. Breadcrust bombs previously described by Wright et al. (2007) are not included in this figure. For a given category, the BSE image magnification increases from left to right. a to c Pumice clast (Sample 98c-1660–38) from the CE 1660 eruption with large voids associated with broken crystals (a) and abundant small, isolated vesicles in between larger ones (red dashed line in c). d to f Banded clasts from the CE 1660 eruption with alternation of zones with high porosity and degree of coalescence and areas dominated with isolated vesicles. g to i Type 1 dense rocks from the 1999 eruption with abundant sheared vesicles organised in clusters. j to l Type 2 dense rocks from the 1999 eruption with diktytaxitic porosity and a silica phase that appears in the interstices between microlites and as large fractured patches in macropores (l). m to o bands of diktytaxitic porosity that alternate with denser bands dominated by cracks in a type 3 dense rock

We observed all three types of dense rocks in deposits from the 1999 block and ash flow deposits. The dome carapace/talus sample studied here corresponds to a type 1 dense rock. Bombs also consist of dense varieties, as described in Wright et al. (2007), and type 1 in nature. Figure 4 illustrates how such variable textures are related to permeability data.

We also note the presence of agglomerated pumice lapilli clasts that contain bands of different porosity and micro-textures and show evidence of ash sintering between pumice clasts (Fig. 4j; cf. Giachetti et al. 2021). This feature was not observed macroscopically. Although this aspect is beyond the scope of our current study, we note that the pumice clast with the lowest permeability corresponds to such agglomerated pumice (Fig. 4).

Crystallinity and crystal textural features in the dacites

Broken phenocrysts are ubiquitous in the pumice clasts (not observed in the dense block and ash flow and dome rocks), with different degrees of breakage ranging from simple cracking to boudinage textures (cf. Giachetti et al. 2010) to completely shattered crystals (Fig. 5a, b). We note that the degree of crystal breakage appears overall lower in the pumice clasts of the 1999 Subplinian activity compared to those from the CE 1660 eruption.

Crystal content measured combining BSE images and EDX mapping yielded a phenocryst/micro-phenocryst crystallinity of ~ 33 vol.% and 20–24 vol.% in pumice clasts from the CE 1660 and 1999 eruptions, respectively, and ~ 48 vol.% in the dense rocks (Fig. 2f). Groundmass microlite content (< 30 µm) is highly variable in the dacites and depends mostly on sample type (Fig. 6). No microlites were observed in pumice clasts from the CE 1660 eruption, where the interstitial films between groundmass vesicles are completely glassy (Fig. 6a). In contrast, rare microlites were observed in the pumice clasts from the 1999 explosive eruptions although the microlite content could not be quantified due to low amount and heterogenous distribution. Microlite contents of 12–23 vol.% have been measured in another sample set (Wright et al. 2022) for the pumice from the 1999 eruption. The microlite content ranges from 9 to 20 vol.% in breadcrust bombs (Wright et al. 2022). High microlite content (mostly plagioclase) was observed for all of the dense rocks studied here. We could not quantify precisely the exact microlite content for the dense rocks at the sample scale because of strong local heterogeneities in micro-textures. We could however quantify a microlite content in the groundmass of approximately 10 vol.% locally. Considering that some areas are fully crystalline, we estimate an average groundmass crystallinity of about 50 vol.%. Plagioclase microlite shape evolves from dominantly equant and tabular in type 1 dense rocks (Fig. 6d) to hopper, skeletal and elongate morphologies in types 2 and 3 dense rocks (Fig. 6e to i).

Fig. 6
figure 6

The variety of groundmass crystallinity in the dacitic rocks. a Microlite-poor groundmass in pumice clasts. b Microlite-bearing bands (mb) that alternate with microlite-free bands (mf) in a banded clast, with an overall low microlite content and spherical vesicles. c Microlite-rich area in the same banded clast as in (b) with tabular shapes for microlites. Note the abundance of oxide microlites in the groundmass. d Groundmass of a type 1 dense rock with dominantly equant and tabular shapes for plagioclase microlites, with rare swallowtail and elongate shapes also observed. e Type 3 dense rocks where hopper, skeletal and elongate morphologies dominate for plagioclase and cohabit with equant oxide microlites. f to i Groundmass of type 2 and 3 dense rocks with similar plagioclase and oxide microlites as in (e) but with an additional silica phase. The amount of silica increases from (f) to (i) and evolves from interstitial silica in between plagioclases to large patches surrounded by a rim of vesicles

In some dense rocks, oxide nanolites are present but are heterogeneously distributed (Cáceres et al. 2020). Dense rocks (types 2 and 3) also commonly contain a silica phase, either in the interstices between plagioclase crystals in the groundmass (Fig. 6e to h), as large fractured patches in macropores (similar to “fish-scale” cracking of cristobalite; Damby et al. 2014) (Fig. 5l) or surrounded by a rim of small vesicles (Fig. 6h, i).

Some banded clasts are characterised by abrupt variations in microlite content and groundmass characteristics, where microlite-poor areas coinciding with spherical vesicles transition to microlite-rich zones with highly deformed vesicles (Fig. 4i, Fig. 5d). We also note that the groundmass in dense areas of banded clasts differ from those of dense rocks, with more abundant oxide microlites in dense bands of at least some banded clasts (Fig. 6c).

Total crystallinities of pumice clasts and breadcrust bombs also differ. Our calculations yield Φx-t = 33 vol.% for the pumice clasts from the CE 1660 eruption (in which no microlites were observed), 20–38 vol.% for the pumice from the 1999 eruption using a range of microlite contents from 0 (this study) to 23 vol.%, and 49–59 vol.% for the interiors of breadcrust bombs (Wright et al. 2022). For the percolation model, we use the average phenocryst content of the pumice clasts of the 1999 eruption (23 vol.%) because the microlite content measured in our study is negligible. For the dense rocks, we use the average microlite content of 50 vol.% discussed previously, which gives Φx-t ~ 70 vol.%.

Vesicle size distribution and vesicle number density

Vesicle size distributions and vesicle number density (Nv) data reveal further differences between pumice from the CE 1660 and 1999 eruptions (Fig. 7). We first note that the pumice from the CE 1660 eruption has a vesicle size distribution (VSD) that extends to smaller vesicle sizes, which is reflected in the high vesicle number density for this pumice (Nv ~ 108 mm−3). The distribution is bimodal with a population of large vesicles that may also correspond to voids associated to phenocrysts. Pumice clasts from the 1999 eruption lack the smallest vesicles, and thus have lower Nv values of 106–107 mm−3. Vesicle size distributions are also bimodal for these pumice.

Fig. 7
figure 7

Vesicle size distribution and vesicle number density of pumice clasts of CE 1660 and 1999. a to c Vesicle size distribution for each clast. d Cumulative vesicle number density distribution for the three pumice clasts. Note the similarity in the vesicle size distribution (VSD) and vesicle number density (Nv) distributions between the two 1999 pumice clasts

Discussion

Here, we use textural and compositional characteristics of clasts from studied eruptions of Guagua Pichincha volcano to better understand magma storage, crystallization, ascent, degassing and gas migration.

Storage conditions, magma ascent and rheological evolution of the dacitic magma

Previous studies yielded estimates of the pre-eruptive physical conditions at Guagua Pichincha volcano. Samaniego et al. (2010) calculated a pre-eruptive temperature of 814 ± 25 °C based on an amphibole-plagioclase geothermometer and a pre-eruptive pressure of 125 ± 68 MPa based on amphibole barometry. This estimate was corroborated by seismic data, indicating a reservoir located around 5 km below the summit (Garcia-Aristizabal et al. 2007), within the saturation depth range calculated based on H2O contents of 2.7–4.6 wt.% measured in glass inclusions within plagioclase and pyroxene phenocrysts in breadcrust bombs and pumice clasts from the 1999–2001 eruption (Wright et al. 2022). Future studies could investigate whether glass inclusion volatile contents (and degassing paths) differ between Plinian and dome-forming eruptions, as they do at eruptions elsewhere around the world (Villemant and Boudon, 1999; Andujar and Scaillet, 2012; Popa et al. 2021).

All dacites analysed here are characterised by a relatively high phenocryst content (range of 20–51 vol.%; Fig. 2f). Such a high phenocryst content can have several effects on conduit processes. First, it may substantially increase bulk magma viscosity and subsequently reduce its ascent rate, especially in the case of dense rocks with > 40 vol% phenocrysts. Such dense magma may come close to a rheological eruptibility limit commonly suggested to be at > 50 vol.% phenocrysts (e.g. Marsh, 1981; Takeuchi, 2011; Popa et al. 2021), although some degree of phenocryst growth likely also occurs within the conduit if ascent is sufficiently slow. We note that the range of water and crystal contents discussed above overlap the conditions at the effusive-explosive transition discussed by Popa et al. (2021). In addition to this rheological effect of crystals, such high phenocryst contents may strongly influence degassing processes during magma ascent.

Effect of crystals on degassing in the dacitic magma

In addition to their rheological effect discussed above, crystals strongly affect degassing processes such as bubble nucleation and coalescence and thereby influence permeability development. The effect of crystals is complex and includes several processes that together may reduce the percolation threshold and enhance outgassing at low porosity and possibly deep in the conduit, including bubble deformation (Shields et al. 2014; Oppenheimer et al. 2015; Lindoo et al. 2017), channelling and gas migration (Oppenheimer et al. 2015; Parmigiani et al. 2017; Collombet et al. 2021), brittle-viscous coalescence (Colombier et al. 2020), porosity redistribution via strain-localization (Laumonier et al. 2011; deGraffenried et al. 2019), heterogenous bubble nucleation and an increase in the bubble number density (Hurwitz and Navon 1994; Burgisser et al. 2017, 2020; Cáceres et al. 2020), and simply enhanced bubble coalescence (Blower, 2001; Cáceres et al. 2022).

Despite this complex and variable effect of crystals on pore connectivity and permeability, we explore here a purely geometric influence of crystals. We first consider the relationship between the melt and bulk porosity for different groundmass crystallinities using the following equations:

$${\Phi }_{\mathrm{m}}={\Phi }_{\mathrm{t}}/\left(1+\left({\Phi }_{\mathrm{x}}({\Phi }_{\mathrm{t}}-1\right)\right)$$
(2)
$${\Phi }_{\mathrm{t}}=\left({\Phi }_{\mathrm{m}}\left(1-{\Phi }_{\mathrm{x}}\right)\right)/\left(1-\left({\Phi }_{\mathrm{m}}{\Phi }_{\mathrm{x}}\right)\right)$$
(3)

where Φm is the melt porosity (excluding crystals, that is the volume of pores divided by volume of pores plus glass phase), Φt is the total porosity (volume of pores divided by the total sample volume) and Φx is the melt crystallinity (vesicle-free). This simple formulation first shows that the difference between the total and melt porosities increases as the melt crystallinity increases (Fig. 8a). Second, we can model the effect of melt crystallinity on connectivity-porosity relationships by correcting both the total and connected porosities for crystal content using Eq. 3. Using this method, we assume that connectivity is only affected by the geometric addition of crystals which is, as discussed before, a simplification of natural systems. We show the influence of crystals on connectivity-porosity relationships for three different percolation models, those of Burgisser et al. (2017), Giachetti et al. (2019), and Vasseur et al. (2020). These models depend on different parameters such as the polydispersivity (parameter S for Vasseur et al. 2020 and σ/d for Burgisser et al. 2017) of the vesicle size distribution, the average bubble diameter (d; Burgisser et al. 2017) and the effective barrier to coalescence (ε; Blower, 2001; Giachetti et al. 2019) (see Online Resource 1 for a more detailed description of these models and their implications).

Fig. 8
figure 8

Effect of crystallinity on melt porosity, porosity–connectivity relationships and the percolation threshold. a Relationship between total and melt porosity when corrected for crystal content (vesicle-free). b to f Geometrical effect of crystallinity on porosity-connectivity trends for the models of Burgisser et al. (2017), Giachetti et al. (2019), and Vasseur et al. (2020) for different input parameters of the models. For each case, the increase in crystal content (number labelled on each curve in vol%) causes a decrease in the percolation threshold. ε is the effective barrier to coalescence used in the Giachetti et al. (2019) model, d and σ are the average (or mean) bubble diameter and polydispersivity from the Burgisser et al. (2017) model, respectively, and S is the polydispersivity parameter from the Vasseur et al. (2020) model

None of these percolation models includes the effect of crystals. We incorporate crystallinity for each of the modelled C–Φ relationships using Eq. 3. We observe that the addition of crystals systematically reduces the percolation threshold Φc (here expressed as a total porosity) (Fig. 8). At a given crystal content, the degree of shift for Φc towards lower porosities depends on the initial values of Φc and therefore on the different model parameters (σ/d, ε, and S). In the following discussion, we will use such values derived from vesicle size distributions of the Guagua Pichincha dacites, as well as measured crystallinities on these dacites, as inputs in the percolation models to compare natural and modelled connectivity–porosity relationships.

Degassing and gas percolation in the vesicular dacitic magma

Pore connectivity and percolation models

Vesicular rocks (pumice clasts and breadcrust bombs) were erupted during explosive events at Guagua Pichincha, including the CE 1660 Plinian eruption and the Subplinian episodes of the 5 and 7 October 1999.

Connectivity variation with porosity in vesicular volcanic rocks provides information about the evolution of bubble connectivity (and the percolation threshold) during vesiculation, and in some rocks as modified by bubble coalescence, deformation, collapse and subsequent outgassing, and the development of micro-fractures (Farquharson et al. 2015; Colombier et al. 2017, 2021; Gonnermann et al. 2017; Bernard and Bouvet de Maisonneuve 2020; Vona et al. 2020; Mitchell et al. 2021). The range of connectivity in explosive products, and inversely the range of isolated porosity (Fig. 3b), reflect the variable contributions of different bubble populations as confirmed by textural observations. BSE images indicate areas with ubiquitous, small, isolated vesicles (~ 1–20 µm) with relatively thick glass films in the interstices between larger pores (Fig. 5c). These large pores can themselves be distinguished between (i) vesicles showing ubiquitous coalescence features with thinner walls (Fig. 5a to c), and (ii) even larger (500 µm to several mm) inter-crystalline void/vesicle space related to crystal breakage (Fig. 5a). Such qualitative observations are confirmed by the VSD data that reveal bimodal size distributions (Fig. 7). We combine the isolated porosity fraction from helium pycnometry with the VSD data to estimate the sizes of isolated vesicles by assuming that they correspond systematically to the smallest vesicles. This yields a critical vesicle diameter threshold of ~ 5 µm that separates isolated and connected vesicles for the analysed pumice from the CE 1660 eruption, and a range of 20–25 µm for the two pumice clasts from the 1999 eruption.

We now compare our measured connectivity relationships for pumiceous samples (those least affected by bubble collapse and outgassing from our sample suite) with estimates from the three percolation models discussed above, corrected for crystallinity. First, the model of Vasseur et al. (2020) was applied using values of vesicle polydispersivity from the VSD data (Fig. 9a). Note that polydispersivity (S) values closer to 0 indicate more polydisperse samples; a value of 1 indicates monodispersivity. S values were 0.19 and 0.17–0.18 for pumice from the CE 1660 and 1999 eruptions, respectively. Average total crystal contents were 33 and 23 vol.% for pumice from the CE 1660 and 1999 eruptions, respectively. We corrected the modelled values of total and connected porosity using the above values of S (described in Colombier et al. (2021) and as based on the Vasseur et al. (2020) model) and measured values of crystallinity using Eq. 3. The modelled values come close to the measured values of C–Φ for pumice clasts and indicate higher connectivity values for the pumice from the CE 1660 eruption (Fig. 9a) at equivalent porosities. As such, this model allows us to include the effect of polydispersivity on the evolution of connectivity in vesicular Guagua Pichincha dacites when corrected for crystal content and is consistent with relative differences between eruptive samples. However, differences between model results for the two eruptive suites are small and do not fit data as well as models below, perhaps because this model assumes spherical and fully overlapping bubbles, which is more adequate to model percolation in static, low viscosity magmas than in the highly viscous rhyolitic melts of Guagua Pichincha dacites studied here. The modelled C–Φ curves indicate a percolation threshold Φc of 0.44 and 0.49 for the pumice clasts from the 1660 and 1999 eruptions, respectively. These values of Φc should be considered with caution because of the previously mentioned model simplifications.

Fig. 9
figure 9

Comparison of porosity-connectivity trends obtained from the different percolation models with natural data. a Comparison with the Vasseur et al. (2020) model. b Burgisser et al. (2017) model. c Giachetti et al. (2019) model. Each modelled curve has been corrected for crystal content given by textural analysis. d Compilation of the modelled curves for pumice clasts of CE 1660 and 1999, and breadcrust bombs of 1999. e Modelled Φ–C curves in the full connectivity range with onset of connectivity at the percolation threshold. For the Vasseur et al. (2020) model, we have non-zero values of connectivity at porosities below the percolation threshold due to finite domain-size effect (see Colombier et al. 2021). The percolation threshold for this model is hence shown with a yellow line at a connectivity of ~ 0.03

The model of Burgisser et al. (2017) also allows us to model the evolution of connectivity with porosity using bubble diameter and polydispersivity from the VSD data (Fig. 9b) and by correcting for crystal content as formerly discussed. This model also fits the natural pumice data relatively well and predicts clear connectivity differences between pumice from the CE 1660 and 1999 eruptions on a C–Φ plot based on input parameters taken from their respective VSDs. Pumice clasts from the 1999 eruption have C–Φ model relationships that imply a higher Φc ~ 0.80 than for pumice clasts from the CE 1660 eruption (Φc ~ 0.58). These are substantially higher values of Φc compared to the Vasseur et al. (2020) model (Φc = 0.44–0.49).

In contrast to these two first models, the model of Giachetti et al. (2019) has input parameters that cannot be derived from the VSD data. We instead modified values of the parameter ε (which varies from 0 to 1, where higher values indicate higher barriers to coalescence) to best fit the natural samples. Giachetti et al. (2019) assessed ε values of 0.20–0.35 for Plinian pumice, a narrower range of 0.25–0.32 yielded a better fit to our data for the pumice clasts from the CE 1660 eruption (Fig. 9c). A higher range of ε = 0.34–0.40 was used to fit the trends of the pumice from the 1999 eruption (Fig. 9c). This model indicates a percolation threshold of 0.48–0.58 and 0.66–0.73 for the pumice clasts from the CE 1660 and 1999 eruptions, respectively. The upper bound values are in good agreement with those obtained using the model of Burgisser et al. (2017). Higher ε for the pumice clasts from the 1999 eruption may be due to a higher melt viscosity and greater resistance to coalescence, which tends to shift Φc to higher values (Blower, 2001). This seems to confirm the previous observation that vesicles are interconnected at larger sizes in the pumice clasts from the 1999 eruption compared to those from the CE 1660 eruption. We used the total range of ε estimated in pumice (0.25–0.40) for the breadcrust bombs and also corrected for the total crystal content, which is higher in the breadcrust bombs (0.49–0.59) than in pumice clasts. We observe that the modelled values fit perfectly the range of C–Φ values for the breadcrust bombs, and suggest a percolation threshold range of 0.37–0.56, on average lower than for the pumice clasts due to the higher crystal content (Fig. 9c and d).

In summary, the combination of these models indicates a lower percolation threshold and consequently higher connectivity values at a given porosity for the pumice clasts from the CE 1660 eruption compared to those from the 1999 eruption, consistent with helium pycnometry data (Fig. 9d). We emphasize that crystals may have exerted additional complex influences on gas percolation than those accounted for in our simple geometric considerations. Textural observations in pumice clasts indicate that on average vesicle walls/films are thinner, vesicles are larger and coalescence is more pervasive in the vicinity of phenocrysts and micro-phenocrysts, as reported for other natural phenocryst-rich dacites at Pinatubo volcano in the Philippines (Hammer et al. 1999) and in experimental phenocryst-bearing foamed obsidians (Cáceres et al. 2022). This implies either that (i) heterogenous nucleation initiated early around these phenocrysts in the conduit (Gardner and Denis, 2004), allowing for more time for growth and coalescence or (ii) nucleation occurred in a single and/or continuous event but growth and coalescence were favoured around crystal phases (Cáceres et al. 2022). Irrespective of these two hypotheses, crystals likely promote early development of system-spanning coalescence in the conduit of Guagua Pichincha and possibly deep outgassing via a range of processes (Collombet et al. 2021).

The evolution of permeability in the vesicular magma

Permeability–porosity relationships (Fig. 4) provide complementary insights on percolation and outgassing processes. Permeability depends on porosity but also on the characteristics of the void space (nature of pores, aperture size, tortuosity, connectivity, anisotropy, presence of fractures) and the comparison to connectivity is therefore not straightforward (Le Pennec et al. 2001; Wright et al. 2009; Farquharson et al. 2015, 2016; Kushnir et al. 2016; Colombier et al. 2017). Permeability of vesicular clasts is on average higher than that of dense rocks (Fig. 4a), even for similar connectivity values (Fig. 4d), which implies that connectivity alone is not responsible for these differences. Vesicular breadcrust interiors are characterised by higher average tortuosity (Wright et al. 2007) values than dense rocks, with a decrease of tortuosity with increasing porosity (therefore during expansion of the breadcrust bomb) (Fig. 4c). A lower tortuosity may explain why two dense rocks have a higher permeability than some vesicular clasts, even if at lower values of connected porosity (Fig. 4b). The expansion pattern observed in breadcrust bombs of a decrease of tortuosity with increasing porosity, indicates that the pumice clasts also have a low tortuosity.

Permeability values of vesicular clasts studied here straddle experimental data obtained by decompression-induced vesiculation of crystal-bearing magmas with 20–40 vol% phenocrysts (Fig. 4a; deGraffenried et al. 2019). The authors showed that the presence of at least 20 vol.% phenocrysts promotes bubble coalescence around phenocrysts, which subsequently acts to increase permeability and reduce the percolation threshold to 45–50 vol.% (total porosity). These observations are consistent with the above inference that connectivity increases and percolation threshold decreases with crystallization, although absolute values of modelled percolation threshold above are higher than those of deGraffenried et al. (2019) due to incorporation of effects of polydispersivity and viscous resistance to coalescence (Burgisser et al. 2017; Giachetti et al. 2019).

DeGraffenried et al. (2019) also compared the timescales of permeable gas escape and magma ascent for vesicular magmas at shallow depth in the conduit (< 25 MPa) and define a critical threshold at a permeability of ~ 10−13 m2. For permeabilities below this threshold, the timescale for outgassing is substantially higher than magma ascent, and therefore no outgassing (nor bubble collapse and porosity reduction) occurs in vesicular magmas. This was likely the case for breadcrust bombs at Guagua Pichincha that were emplaced at such shallow conditions and were impermeable at the time of fragmentation because the rapidly quenched rims consist of isolated vesicles only and porosities close to zero (Wright et al. 2007). Pumice clasts were likely fragmented at a depth greater than 25 MPa and they have permeabilities in the range of 10−14–10−13 m2 (Fig. 4a). However, these pumice clasts measured for permeability correspond to coarse lapilli (> 32 mm in diameter) or bombs and therefore likely experienced some post-fragmentation vesiculation (Mitchell et al. 2018; Trafton and Giachetti, 2021). In addition, final porosities observed in vesicular clasts generally do not reflect pre-explosive (lower) porosities of magma in the conduit (Burgisser et al. 2019; Collombet et al. 2021). Taken together, this indicates that the vesicular magma has lower porosity, connectivity and permeabilities (lower than 10−14 m2) at the time of fragmentation as compared to that preserved in pyroclasts. We can therefore consider this vesicular magma as largely impermeable at the time of fragmentation, a status that promoted bubble overpressure, fragmentation, and explosive activity. Deeper in the conduit, the timescale for magma ascent is much greater, and therefore deep permeable gas escape via bubble coalescence and channelling is plausible as discussed previously (Parmigiani et al. 2017; Collombet et al. 2021).

Degassing, percolation and outgassing in the dense dacitic magma

Dense rocks at Guagua Pichincha were mostly produced during dome-forming eruptions although they are also encountered in the deposits of Plinian eruptions. Their emplacement and associated eruptive style will be considered here as effusive (emplacement of lava domes), although they can also have a hybrid history between this initial effusive emplacement and subsequent explosive events, as in the case of dome-destructing Vulcanian events or dome-collapse block-and-ash pyroclastic currents (e.g. Colombier et al. 2017; Cassidy et al. 2018; Popa et al. 2021).

The textural properties of the dense rocks are more complex than those of the vesicular clasts, and porosity in this lithology comprises vesicles with variable degrees of shearing, diktytaxitic porosity, large macropores associated with phenocrysts and micro-cracks (Figs. 3, 4 and 5). These different void space geometries confer variable throat sizes, tortuosities, connectivities and consequently permeability on the magma. In addition, the nature of this void space changes substantially as a function of porosity and the type of dense rock.

The dense rocks are characterised by overall high pore connectivity (except for one sample) at a low porosity window compared to vesicular clasts, which indicates a low percolation threshold, possibly due to a combination of shearing, high crystal content and an abundance of cracks (Colombier et al. 2017, 2020; Lamur et al. 2017) or pore collapse following vesiculation (Colombier et al. 2017). The higher estimated total crystallinity (~ 70 vol.%) coupled to the geometric effect of crystals on the C–Φ trends (Fig. 8) implies a decrease of the percolation threshold. Crystals also tend to increase connectivity and reduce Φc via several more complex processes such as channelling, enhanced deformation and/or brittle-viscous coalescence (e.g. Oppenheimer et al. 2015; Parmigiani et al. 2017; Colombier et al. 2020).

Shearing likely also contributed to reduction in percolation threshold. Shearing has been shown to reduce the percolation threshold during experimental vesiculation of rhyolitic melts with complex implications for eruptive style (Okumura et al. 2013). Textural evidence for shearing has been observed in the type 1 dense rocks (Fig. 4g; Fig. 5gi). In addition, shear-localization in crystal-rich magma can cause porosity redistribution and enhance magma outgassing at very low porosity (Laumonier et al. 2011). Such porosity redistribution may produce banding in which high, diktytaxitic porosity domains alternate with low porosity domains, such as observed in the type 3 dense rocks (Fig. 5mo). Lastly, meso- and micro-scale cracks also tend to promote an increase of connectivity, permeability and a reduction of the percolation threshold (Heap et al. 2014; Heap and Kennedy, 2016; Colombier et al. 2017; Kushnir et al. 2017; Lamur et al. 2017), recalling that post-emplacement thermal cracking can be substantial.

Despite these variable pore properties and processes, the continuous increase of permeability with porosity indicates a transition between crack-dominated, diktytaxitic to vesicle-dominated pore networks in this population of dense rocks (Figs. 4 and 5), also proposed by Le Pennec et al. (2001) and Kushnir et al. (2016) for Merapi volcano (Indonesia). The permeability–porosity relationship also indicates a low percolation threshold achieved by the range of processes discussed above. As tortuosity is rather low in the dense rocks, we propose that the range of permeability is the result of variations in void space aperture size, decreasing from vesicle clusters to diktytaxitic pores to cracks.

Evidence for textural mingling at Guagua Pichincha and the origin of banded clasts

Banded textures were commonly observed at both macroscopic and microscopic scales in the different lithologies from both the 1660 and 1999 eruptions. We interpret this banding as a textural mingling with (i) local variations in groundmass crystallinity that subsequently affect vesicle shape (Fig. 4i) and (ii) variations in the degree of coalescence and vesicle connectivity (Fig. 5e, f). Areas with a lower microlite content are commonly characterised by spherical and/or smooth vesicle shapes, which alternate with microlite-rich areas and more irregular vesicles (Fig. 4i). Evidence for such textural mingling between dense, microlite-rich areas and more pumice-like textures also comes from petrophysical measurements. The measured values of porosity, connectivity, permeability and tortuosity indicate trends that are intermediate between those of the vesicular clasts (Breadcrust bombs and pumice) and the dense rocks (Figs. 3 and 4). Such textural mingling near the surface has been previously proposed for dome-forming eruptions and effusive-explosive transitions, possibly due to differences in ascent rate between the dense stalling magma and vesicular, volatile-rich magma underneath (Cassidy et al. 2015). Rosi et al. (2004) explained similar textural mingling at Quilotoa volcano (Ecuador) by shearing and re-heating near the conduit walls. This may explain the higher oxide microlite content observed in banded clasts. Despite this apparent textural mingling, we do not rule out an initial chemical mixing, but this aspect is beyond the scope of our study and could be further investigated.

Implications for eruptive cycles at Guagua Pichincha

From the previous discussion, we can clearly show that there were different processes and conditions prevailing in the conduit, manifest as different eruptive styles during the CE 1660 and 1999–2001 eruptions and although not specifically discussed here, some of which may be revealed through compositional differences between eruptive suites. The difference between explosive Plinian and mostly effusive, dome-forming eruptions is usually inferred to be driven by variations in magma ascent rates (Rutherford, 2008; Cassidy et al. 2018). Faster ascent rate in the case of explosive eruptions may reduce gas escape and time for microlite crystallization during ascent in the conduit, promoting fragmentation by overpressure and an acceleration in explosive activity (Cassidy et al. 2018). Although we have no direct constraint on absolute values for ascent rate, we will use textural observations to provide some elements of discussion.

We here link the insights gained on degassing processes occurring in the conduit for the different lithologies to eruptive style during eruptive cycles at Guagua Pichincha. We categorize eruptions in two endmembers: (i) highly explosive pumice-rich Plinian eruptions and (ii) dominantly effusive, dense-dominated, dome-forming eruptions. We treat the October 1999 Subplinian eruptions as a pumice-rich eruptive style with processes and products that differ from those of Plinian and dome-forming activity.

We propose a schematic eruptive model that links storage conditions, ascent rate, conduit processes and ultimately eruptive style at Guagua Pichincha volcano (Fig. 10). Although this model is schematic, it allows us to describe the parameters controlling century-long eruptive cycles at Guagua Pichincha.

Fig. 10
figure 10

Schematic model of conduit processes explaining the differences in eruptive style at Guagua Pichincha volcano and the link with magma properties at storage conditions. Φxi: initial crystal content (phenocryst content). Φmicrolites: Microlite content of the groundmass. Φc: modelled percolation threshold. ε: effective barrier to coalescence fitted from the Giachetti et al. (2019) model. Nv: Vesicle number density. dP/dt: decompression rate. BCB: breadcrust bomb

Plinian eruptions

The deposits from the CE 1660 eruption are largely dominated by pumice clasts, with a minor amount of dense rocks, banded clasts and breadcrust bombs. A juvenile origin for these dense pyroclasts is questionable despite a rather gradual textural transition in density towards the pumice population. The presence of dense clasts in the CE 1660 Plinian deposits indicates that (i) a plug/dome formed at the onset of the eruption, (ii) densification occurred along the conduit margins during ascent or (iii) pre-existing domes were destroyed and excavated by the CE 1660 event. Plinian phases are commonly associated with the destruction of pre-existing lava domes (Maeno et al. 2019), which was likely the case for the CE 1660 eruption of Guagua Pichincha. Dense pyroclasts have been observed in all Plinian deposits (Fig. SM1 in Online Resource 1) and they are texturally indistinguishable to those of the 1999–2001 dome-forming phases. Despite all these considerations, a juvenile origin for these dense rocks cannot be guaranteed for the Plinian case. Future studies could investigate potential chemical markers that may allow us to link these dense clasts to some specific eruption cycle. Regardless of the origin of these dense clasts, we here discuss the eruptive processes during Plinian eruptions in the light of pumice clasts as these are the dominant type of lithology for these eruptions.

Pumice clasts of Plinian eruptions are generally characterised by a paucity of microlites, and connectivity and permeability data that indicate a high percolation threshold at ~ 0.58–0.68, but lower than for the 1999 pumice-rich explosive events, likely due to a lower viscosity yielding less resistance to coalescence. Plinian pumice clasts are also characterised by a high vesicle number density and ubiquitous broken crystals, likely reflecting high decompression rates in the dacitic magma during Plinian eruptions. Hence, we propose that ascent rates were high during Plinian eruptions, limiting the potential for outgassing and promoting explosive fragmentation. Although the present study is more focused on the CE 1660 eruption for the Plinian case, we may as a first approximation generalize our findings to the first and tenth century Plinian eruptions as textural observations and connectivity data indicated overall similar patterns.

The October 1999 pumice-rich eruptions

Pumice-bearing explosions were also produced during Subplinian activity on the 5 and 7 of October 1999. Pumice from the October 1999 explosive eruptions are characterised by different micro-textures, vesicle size distributions, a higher isolated vesicle fraction and higher inferred ε values than CE 1660 pumice, which indicate percolation in the presence of a more viscous melt than was the case for the CE 1660 pumice. This observation, coupled with the presence of more abundant microlites, a lower degree of crystal breakage (Miwa and Geshi, 2012) and a lower vesicle number density (Toramaru, 2006) in these pumice clasts, would be consistent with a model of higher initial viscosity and slower decompression rates.

Effusive dome-forming eruptions

The eruptive phases of the 1999–2001 eruption, except the two October 1999 pumice-rich eruptions, produced dominantly dense dacitic rocks. The dense rocks from the 1999–2001 eruptive period originate from dome-forming phases that were punctuated by frequent dome collapse and block and ash flow generation and Vulcanian eruptive pulses (Wright et al. 2007).

Dense rocks have a broad range of connectivity and permeability at low porosity that indicate a low percolation threshold. They also have highly variable textures ranging from crack-dominated, diktytaxitic to porous networks with clusters of coalesced, sheared vesicles. These features are typical of crystal-rich intermediate magmas with rhyolitic melt composition (Kushnir et al. 2016). Outgassing in the dome is favoured by the high crystal contents via a range of processes. In addition to the geometric influence of crystals on percolation, we identified the following effects that promote outgassing and a low percolation threshold in the dense rocks: (i) cracking (Kushnir et al. 2017; Lamur et al. 2017), (ii) shearing likely along conduit walls and in the dome carapace (Okumura et al. 2013), (iii) porosity redistribution via shear-localization (Laumonier et al. 2011) and (iv) enhanced coalescence around phenocrysts (Cáceres et al. 2022). All these processes are enhanced in crystal-rich magmas.

These dense magmas were also characterised by substantially higher phenocryst content (Fig. 2) than other lithologies, causing a higher bulk viscosity. They were therefore likely emplaced with lower ascent rates compared to the dacitic magma feeding Plinian to Subplinian eruptions. These lower ascent rates and decompression at shallow levels of the conduit with longer time available for crystallization caused a higher microlite content, which would have further increased the bulk magma viscosity. A combination of slower ascent and high phenocryst content may favour outgassing even in the deepest parts of the magma column (Collombet et al. 2021) thereby reducing the explosive potential and promoting dome-forming activity.

Conclusion

This study describes conduit processes of the CE 1660 Plinian eruption and the 1999–2001 dome-forming to Subplinian eruption of Guagua Pichincha volcano, Ecuador. We found that the difference between explosive pumice-rich Plinian or Subplinian phases and dome-forming activity was associated with an increase of phenocryst content for effusive phases. This led to a rheological change with an increase in viscosity and consequently slower magma ascent in the case of dome-forming eruptions. In addition, the higher crystallinity influenced magma degassing by decreasing the percolation threshold promoting permeable outgassing and effusive activity. On the other hand, we infer that magma during pumice-rich eruptions ascended faster and remained mostly impermeable until fragmentation due to a high percolation threshold. Connectivity–porosity relationships and textural analysis further revealed differences between pumice clasts erupted during Plinian and Subplinian activity, with a faster magma ascent indicated for the former case as evidenced by the higher vesicle number density, ubiquitous broken crystals, and percolation models that indicate a lower resistance to coalescence.

In this study, we have made progress in understanding the role of conduit processes on shifts in eruptive style and cyclicity at Guagua Pichincha volcano. More broadly, our findings provide insights into vesiculation and percolation in crystal-rich, intermediate magmas and implications for eruptive transitions at silicic volcanoes worldwide. Remaining questions include: What controlled the rheological variations at the onset of the different eruptive phases? What was the influence of small, but clear bulk compositional differences between products of the CE 1660 and 1999–2001 eruptions? What was the importance of magma recharge and mixing on these variations? More work to link triggering processes in the storage region to magma ascent and conduit processes studied here would be useful. Future work could therefore focus on a more detailed petrological analysis to better constrain differences in recharge and mixing for a range of eruptive style at Guagua Pichincha. Volatile contents of the pumice clasts from the CE 1660 eruption could also be further investigated in order to compare with the 1999 products. Together, these studies would allow us to better assess eruptive processes and related hazards at Guagua Pichincha volcano.