Micro-Textural Controls on Magma Rheology and Vulcanian Explosion Cyclicity

Understanding the relationship between degassing, crystallization processes and eruption style is a central goal in volcanology, in particular how these processes modulate the magnitude and timing of cyclical Vulcanian explosions in intermediate magmas. To investigate the influence of variations in crystal micro-textures on magma rheology and eruption dynamics, we conducted high-temperature (940°C) uniaxial compression experiments at conditions simulating a shallow volcanic conduit setting on eight samples of high-crystallinity andesite with variable plagioclase microlite populations from the 2004 to 2010 Vulcanian explosions of Galeras volcano, Colombia. Experiments were conducted at different strain rates to measure the rate-dependence of apparent viscosities and assess the dominant deformation processes associated with shear. Variations in plagioclase micro-textures are associated with apparent viscosities spanning over one order of magnitude for a given strain rate. Samples with low numbers of large prismatic microlites behaved consistently with published rheological laws for crystalline dome samples, and displayed extensive micro-cracking. Samples with high numbers of small tabular microlites showed a lower apparent viscosity and were less shear-thinning. The data suggest a spectrum of rheological behavior controlled by concurrent variations in microlite number, size and shape. We use previously published micro-textural data for time-constrained samples to model the apparent viscosity of magma erupted during the 2004–2010 sequence of Vulcanian explosions and compare these results with observed SO2 fluxes. We propose that variations in magma decompression rate, which are known to produce systematic textural differences in the plagioclase microlite cargo, govern differences in magma rheology in the shallow conduit. These rheological differences are likely to affect the rate at which magma densifies as a result of outgassing, leading to magmatic plugs with a range of porosities and permeabilities. The existence of magmatic plugs with variable physical properties has important implications for the development of critical overpressure driving Vulcanian explosions, and thus for hazard assessment during volcanic crises. We suggest a new conceptual model to explain eruption style at andesitic volcanoes based on micro-textural and rheological differences between “plug-forming” and “dome-forming” magma. We advance that existing rheological laws describing the behavior of andesitic magma based on experiments on dome rocks are inappropriate for modeling large Vulcanian explosions ( ∼ 106 m3), as the magma involved in these eruptions lacks the characteristics required to form exogenous lava domes.

Understanding the relationship between degassing, crystallization processes and eruption style is a central goal in volcanology, in particular how these processes modulate the magnitude and timing of cyclical Vulcanian explosions in intermediate magmas. To investigate the influence of variations in crystal micro-textures on magma rheology and eruption dynamics, we conducted high-temperature (940°C) uniaxial compression experiments at conditions simulating a shallow volcanic conduit setting on eight samples of high-crystallinity andesite with variable plagioclase microlite populations from the 2004 to 2010 Vulcanian explosions of Galeras volcano, Colombia. Experiments were conducted at different strain rates to measure the rate-dependence of apparent viscosities and assess the dominant deformation processes associated with shear. Variations in plagioclase micro-textures are associated with apparent viscosities spanning over one order of magnitude for a given strain rate. Samples with low numbers of large prismatic microlites behaved consistently with published rheological laws for crystalline dome samples, and displayed extensive micro-cracking. Samples with high numbers of small tabular microlites showed a lower apparent viscosity and were less shear-thinning. The data suggest a spectrum of rheological behavior controlled by concurrent variations in microlite number, size and shape. We use previously published micro-textural data for time-constrained samples to model the apparent viscosity of magma erupted during the 2004-2010 sequence of Vulcanian explosions and compare these results with observed SO 2 fluxes. We propose that variations in magma decompression rate, which are known to produce systematic textural differences in the plagioclase microlite cargo, govern differences in magma rheology in the shallow conduit. These rheological differences are likely to affect the rate at which magma densifies as a result of outgassing, leading to magmatic plugs with a range of porosities and permeabilities. The existence of magmatic plugs with variable physical properties has important implications for the development of critical overpressure driving Vulcanian explosions, and thus for hazard assessment during volcanic crises. We suggest a new conceptual model to explain eruption style at andesitic volcanoes based on microtextural and rheological differences between "plug-forming" and "dome-forming" magma.

INTRODUCTION
Hazardous sequences of Vulcanian explosions are a common feature during periods of activity at arc volcanoes. These explosions are commonly associated with the development of overpressure beneath or within low-permeability, highcrystallinity magmatic plugs in the upper conduit (e.g., Sparks, 1997;Stix et al., 1997;Voight, 1999;Mueller et al., 2005;Clarke et al., 2007;Wright et al., 2007;Mueller et al., 2011a;Lavallée et al., 2015). Plugs most commonly develop in intermediate magmas (andesites and dacites) and are thought to result from decompression-induced degassing and crystallization during magma ascent from a shallow crustal reservoir (e.g., Hammer et al., 1999;Cashman and Blundy, 2000;Hammer et al., 2000;Preece et al., 2016;Bain et al., 2019a), which produces a large increase in viscosity (Sparks, 1997;Lavallée et al., 2007;Clarke, 2013), inhibiting magma ascent and rendering it prone to rupture or fragment (Lavallée and Kendrick, 2020; and references therein). Concurrently, magma permeability is reduced by outgassing and densification during compactant shear Ashwell and Kendrick et al., 2015) and fracture closure due to healing and tuffisite infills (e.g., Castro et al., 2012;Kolzenburg et al., 2012;Kendrick et al., 2016;Lamur et al., 2019), which leads to a large increase in strength (Spieler et al., 2004;Vasseur et al., 2013;Coats et al., 2018;Lamur et al., 2019). The degassing, crystallization and densification processes responsible for plug and overpressure development may be modulated by differences in magma composition, ascent rate and decompression path, petrological characteristics, rheological evolution, outgassing and shear history (e.g., Calder et al., 2015). This results in a spectrum of physico-chemical evolution pathways that promotes a variety of repose times between explosions, a range of ejected volumes and transitions from explosive activity to effusive dome-building behavior (e.g., Cassidy et al., 2018;Wallace et al., 2020).
However, the detailed mechanisms by which ascent, degassing, crystallization and densification processes may be modulated to govern this rich variety in eruptive style are complex, inter-linked, and remain incompletely understood. In intermediate magmas, crystal micro-textures of anhydrous phases such as plagioclase are known to track magma ascent processes and represent an important record of shallow sub-surface dynamics (e.g., Hammer et al., 1999;Cashman and Blundy, 2000;Hammer et al., 2000;Miwa et al., 2009;Wright et al., 2012;Preece et al., 2016;Cashman, 2020;Gaunt et al., 2020;Wallace et al., 2020). Decompression experiments have shown that the final plagioclase micro-textures are a combined result of the decompression rate, decompression path (e.g., single-step or multi-step ascent), final storage pressure, time spent at that pressure, and the composition of the exsolving volatile phase (Hammer and Rutherford, 2002;Couch et al., 2003;Brugger and Hammer, 2010a;Brugger and Hammer, 2010b;Riker et al., 2015;Cichy et al., 2017). In general, low decompression rates result in low degrees of effective undercooling (ΔT) in magmas, defined as the difference between the liquidus temperature and the subliquidus magmatic temperature (Kirkpatrick, 1981), setting a small driving force for crystallization. This results in low crystal nucleation rates and high crystal growth rates, producing microlite populations characterized by low number densities (i.e. the number of crystals per unit area or volume) and large crystal sizes, e.g., the summer 1980 explosive eruptions of Mount St Helens, United States (Cashman, 1992) and the 2006 dome-forming eruption of Merapi, Indonesia (Preece et al., 2013). Higher decompression rates result in higher ΔT, higher nucleation rates and lower growth rates, producing microlite populations characterized by higher number densities and smaller sizes, e.g., the 1991 pre-climactic eruptions of Mount Pinatubo, Indonesia (Hammer et al., 1999). Crystal morphologies change from facetted crystals at low ΔT to skeletal morphologies at higher ΔT, due to an increasing disparity between the rate of crystal growth and the rate of supply and rejection of chemical components to/from the crystal interface (Lofgren, 1974). In addition, for facetted crystals the growth rate between different crystal axes becomes more disparate with increasing ΔT (Holness, 2014), leading to more prismatic microlites at high growth rates (Hammer et al., 1999;Bain et al., 2019a). Differences in the phenocryst population resulting from variable rates of decompression and degrees of undercooling may also be expected, such as the extent of anhydrous phenocryst growth (Cashman, 2020) and the extent of breakdown of hydrous phases such as amphibole (Devine et al., 1998).
The properties of the crystal cargo of a magmatic suspension, e.g., the crystal fraction (Caricchi et al., 2007;Ishibashi and Sato, 2007;Lavallée et al., 2007), crystal aspect ratio (Mueller et al., 2011b;Klein et al., 2017) and crystal size distribution (Klein et al., 2017;Klein et al., 2018), are also known to strongly affect its rheology, by impacting the flow field due to particle-melt and particle-particle interactions. The attributes of the crystal cargo control the maximum packing fraction of the magmatic suspension (ϕ m ), which is the maximum crystal fraction at which the suspension partly locks up (Einstein, 1911;Roscoe, 1952) and adopts a different rheology encompassing deformation mechanisms such as crystal plasticity, grain sliding and microcracking Lavallée et al., 2012;Kendrick et al., 2013;Kendrick et al., 2017). Thus, even at a given crystallinity, the geometry and distribution of crystals generated during magma ascent impacts the extent of crystal interaction, exerting an important rheological control on magma behavior in the shallow conduit and during eruption. In the context of plug formation and dome-building activity at arc volcanoes, a micro-textural control on rheology (including for example, viscosity and shear-thinning behavior) can be expected to affect the rate at which magma densifies to form a lowpermeability plug (Ashwell and Kendrick et al., 2015;Bain et al., 2019a), which has important implications for outgassing behavior, overpressure development, and eruption style, size and timing. The porosity of magma may also exert an influence on its rheology, as isolated gas bubbles may act to increase or decrease viscosity depending on a competition between deforming and restoring forces (e.g., Llewellin and Manga, 2005), while connected porosity dominated by interconnected vesicles and cracks may form a pathway for gas escape and reduce shear viscosity as well as pore pressure (Kendrick et al., 2013). The original porosity of the magma has also been inferred to control the strain-dependence of viscosity evolution associated with densification (Wright and Weinberg, 2009;Ashwell and Kendrick et al., 2015). Thus, petrological and micro-structural characteristics and the rheological evolution occurring during densification and plug formation may exert an important control on the plug characteristics and associated eruptive behavior. A thorough understanding of these relationships is clearly required to aid in interpreting monitored geophysical and geochemical signals, and thus inform risk management during volcanic crises.
In this work, we investigate the effect of varying crystal microtextures on magma rheology by conducting a series of hightemperature (940°C) deformation tests on eight bomb samples ejected during Vulcanian explosions from Galeras volcano during the 2004-2010 period of activity. Experimental conditions were chosen to replicate a shallow volcanic conduit environment where intermediate magmas commonly stall, form plugs or feed dome-building eruptions. Prior to the experiments, these samples were texturally characterized to quantitatively assess the properties of plagioclase microlites, as the dominant groundmass phase. The connected porosity and permeability of the samples were measured before and after deformation to monitor the micro-structural changes that took place during the experiments. Phenocryst proportions were estimated, groundmass glass compositions were analyzed, and the viscosity of the melt phase was calculated in order to isolate the contribution of microlites to the experimentally-determined apparent magma viscosity. The experimental results were then used to calibrate an empirical viscosity model to estimate the high-temperature, strain rate-dependent viscosity of Galeras magma over the 2004-2010 eruptive period. We then compare the rheology of the Galeras magma in different phases of the eruptive period with the observed SO 2 fluxes and attempt to reconcile our multiparametric dataset into a conceptual model linking evolving eruption dynamics (i.e. changes in the volume of erupted material, the repose time prior to explosions and the presence or absence of a dome) with evolving magma rheology, modulated by changing crystal micro-textures due to magma decompression rate variations. This model sheds light on the transition from effusive dome-building behavior accompanied by occasional, small Vulcanian blasts to frequent, large Vulcanian explosions at Galeras volcano, and may have important implications for hazards and micro-textural monitoring during periods of high and low magma decompression rate at Galeras and other similar arc volcanoes.

Sample Selection
The 2004-2010 period of activity at Galeras volcano ( Figure 1) comprised 17 Vulcanian explosions (VEI 2-3) and two periods of dome-building activity (in 2006 and 2008), interspersed with numerous ash and gas venting events (up to VEI 1; Vargas and Torres, 2015;Narváez Medina et al., 2017). Based on observations by the Colombian Geological Survey (Servicio Geológico Colombiano -SGC), explosions were Vulcanian in character and short-lived (typically lasting several minutes) with plume heights reaching up to 12 km following the most powerful explosions. The minimum volume of material ejected during each Vulcanian explosion, calculated from isopach maps established by the SGC, ranged from 3.5 × 10 4 to 3.56 × 10 6 m 3 and the repose time between explosions ranged between 1 and 554 days (Vargas and Torres, 2015). The minimum total volume of material produced during this period was 21 × 10 6 m 3 , with 16 × 10 6 m 3 ejected during Vulcanian explosions and 5 × 10 6 m 3 erupted in the form of domes. Dome growth was first observed on January 13, 2006 ( Figure 1D). Subsequent explosions partially destroyed this dome and renewed dome growth was noted on September 18, 2008. No pyroclastic flows or associated deposits were observed during this period. Eruptive products were limited to fallout deposits and ballistics.
We selected eight samples of high-crystallinity andesite bombs ( Figure 2A) from the study of Bain et al. (2019b), which represent juvenile blocks produced by the 2004-2010 Vulcanian explosions, collected in July 2017 with the support of the SGC. Following the explosions in 2004-2010, SGC personnel marked juvenile pyroclasts in the field using spray paint, which allowed us to identify these samples during our 2017 fieldwork. The samples were collected from the amphitheater rim, at a distance of ∼ 650-700 m from the vent ( Figure 1C). Samples were chosen in order to cover the full range in density identified in previous work on Galeras ballistics and observed in the field, and range in total porosity from 0 to 26% (Bain et al., 2019b). These bombs correspond to the "dense" and "scoriaceous" types described by Bain et al. (2019a), which were shown to have a low water content (<0.36 wt% H 2 O) in the interstitial glass, preventing them from vesiculating on the short timescales of Vulcanian eruptions (e.g., Wright et al., 2007). These bombs are therefore assumed to be representative of the properties of the magma stored in the shallow conduit prior to Vulcanian explosions at Galeras volcano during 2004-2010. Breadcrusted bombs were also generated during this sequence of explosions (i.e., the 'inflated type described by Bain et al., 2019a) but these were excluded from the present study due to the occurrence of syn-eruptive vesiculation, which overprints the physical properties of magma residing in the shallow conduit. A detailed description of bombs from this eruptive period is given in Bain et al. (2019a). Briefly, the selected samples consist of porphyritic andesite with plagioclase (20-25%; 0.5-6 mm length in section), orthopyroxene + clinopyroxene (7-10%; 0.5-5 mm) and Fe-Ti oxide phenocrysts (3-5%; 0.5-1 mm) and microphenocrysts (0.25-0.5 mm for all phases) with rhyolitic interstitial glass. Amphibole phenocrysts are rare and always show an opaque breakdown rim. Microlites (<0.25 mm) consist dominantly of plagioclase, as well as Fe-Ti oxides and pyroxene. The vesicularity of the selected samples varies widely (3-23% area) in thin sections, spanning the range of vescularity observed in these bomb types at Galeras, with vesicles ranging in size from a few tens of microns to 3-4 mm, typically with complex, branching/ polylobate shapes. A study of the porosity, permeability and porous micro-structure of these bombs can be found in Bain et al. (2019b).

High-Temperature Rheology Experiments
Parallel-plate rheometry was employed to measure the viscosity, and its strain rate dependence, of the selected samples. Cylinders 26 ± 0.22 mm in diameter were prepared from each sample ( Figure 2B), with lengths varying between 40.66 and 53.76 mm. The ends of the cylinders were ground parallel, and cores were then cleaned in an ultrasonic bath and oven-dried. Cores were deformed at high temperature in a 100 kN Instron 8862 uniaxial press in the Volcanology and Geothermal Research Laboratory at the University of Liverpool, following the instrument protocol described in Hess et al. (2007). In each experiment, a sample was heated to ∼ 940°C at a rate of 3°C/ min, using a three-zone, split cylinder furnace. Temperature was monitored via a thermocouple set in contact with the sample and was recorded at the beginning and end of each experiment (the average temperature ranged between 924 and 939°C in the different experiments). Following a dwell time of 1 h to ensure that the sample and the pistons had thermally equilibrated, the sample was deformed sequentially at three constant strain rates:  Medina et al., 2017). Red bars correspond to explosions for which time-constrained samples were collected by the SGC. Periods of dome growth, dome presence in the crater (with no continued growth) and explosions that partially or fully destroyed the domes are also indicated, based on SGC observations. Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 45 min at 10 -6 s −1 , 30 min at 10 -5 s −1 , and 15 min at 10 -4 s −1 . The applied strain rates were chosen to reflect the range of rheological conditions applicable to magma deformation in the shallow conduit, spanning the near-static ()10 -6 s −1 , below which magma rheology is largely Newtonian; e.g., Caricchi et al., 2007) to significantly non-Newtonian regime (10 -6 ∼10 -4 s −1 ; Caricchi et al., 2007) and toward the transition to brittle deformation conditions (U10 -4 s −1 ; e.g., Coats et al., 2018). The sample load and extension were recorded during the experiments. Extension was then corrected for the compliance of the apparatus using a priori correction measurements at the tested conditions in the absence of a sample (e.g., Gruber, 2018). The stress (σ, in Pa) imparted on the sample and the resulting strain (ε) during the experiment were then used to calculate the apparent shear viscosity (η a , in Pa s) at each strain rate once the sample began to deform steadily. η a was calculated using the equation of Gent (1960) for parallel-plate rheometry, adapted for use under constant strain rate (Coats et al., 2018): where F is the force (in N), h is the sample height (in m), _ ε is the strain rate (in s −1 ) and V is the initial sample volume (in m 3 ), which is assumed constant throughout deformation. Cores were subsequently cooled at a rate of 3°C/min to avoid thermal cracking ( Figure 2C shows examples of cores following rheology tests). The length and width of the cores were measured with high precision calipers before and after the rheology tests to calculate the volume of the original cylinders and post-deformation barrels, and determine changes in volume (Supplementary Material A). Table 1 provides a list of selected relevant symbols used throughout the text.

Porosity and Permeability Measurements
The connected porosity and permeability of each core were measured before and after the high-temperature experiments. The connected porosity was measured by gas displacement pycnometry at the University of Liverpool (initial measurements were published in Bain et al., 2019b), using an AccuPyc II 1340 helium pycnometer developed by Micromeritics. In this method, a core is sealed inside a chamber of known volume (100 cm 3 ). The instrument calculates the skeletal volume of the sample (i.e. the volume of rock plus isolated pores) by comparing the volume of gas injected to bring the chamber to a pressure of 21 psi, to the volume of the chamber calibrated at the same gas pressure. The skeletal volume returned by this method is accurate to 0.1%. The volume of pores connected to the outside of the core and allowing helium infiltration (the connected porosity, ϕ c ) is then calculated by taking the ratio of the skeletal volume to that of a perfectly solid cylinder with the core's dimensions.
Steady-state permeability (κ, in m 2 ) at effective pressures ranging from 1.91 to 2.39 MPa (confining pressure ∼6 MPa; Supplementary Material A) was measured at the University of Liverpool using a steady-state hydrostatic cell developed by Sanchez Technologies. In this setup, the confining pressure is applied using silicon oil around a Vitton jacket containing the sample assembly, while water is pumped through the sample using two servo-controlled pore pressure pumps. Darcian conditions (i.e., a laminar flow through the sample) were ensured by imposing a pressure differential (ΔP) of 0.5-1 MPa across the sample (depending on the permeability of the sample, with higher ΔP required to obtain a detectable flow in the densest samples) by decreasing the outflow pore pressure. The resulting flow rate was measured and used to calculate the Darcian permeability, as well as the Forcheimer (Whitaker, 1996) and Klinkenberg (Klinkenberg, 1941) coefficients.

Glass Composition and Melt Viscosity Estimation
Polished thin sections of each sample were prepared and groundmass glass major element compositions were measured by Electron Probe MicroAnalysis (EPMA) at the University of Edinburgh, using a Cameca SX-100 electron microprobe. Analyses were performed using an accelerating voltage of 15 kV, a beam current of 1 nA and a defocused beam of 5 μm diameter. Alkalis were measured first to avoid their loss by migration. A natural rhyolitic glass from Lipari was used for initial calibration and was analyzed on a daily basis to monitor precision and accuracy during the analyses. Totals outside the range 99-101% were discarded from the dataset. Eight analyses fitting this criterion were collected for each sample. The viscosity of the melt phase at the experimental temperature (940°C) was estimated using the results of each chemical analysis as input parameters in the viscosity calculator of Giordano et al. (2008); here, 0.36 wt% H 2 O was used as the water concentration for the groundmass glass of "dense" and "scoriaceous" bombs, which corresponds to the maximum constrained for these eruptive products from the 2004 to 2010 period at Galeras volcano (as measured through secondary ion mass spectrometry by Bain et al., 2019a), and the concentration of fluorine was simplified to be zero for all samples. The obtained melt viscosity values were then averaged for each sample. Using the melt viscosity (η m ) for each sample and the apparent viscosity (η a ) from the rheometry experiments, the relative viscosity (η r ) was then calculated at each applied strain rate, using η r η a / η m (Caricchi et al., 2007), to isolate the contribution of the crystals to the rate-dependent shear viscosity of Galeras magmas.

Textural Analysis of Plagioclase Microlites
High magnification (×1,500) Back-Scattered Electron (BSE) images of the groundmass in sample thin sections were collected using a Carl Zeiss SIGMA HD VP Field Emission Scanning Electron Microscope (SEM) at the University of Edinburgh. The BSE images were stitched using Zeiss SmartTiffV3 software to provide one continuous area of groundmass for each sample, and plagioclase microlites were outlined manually using Inkscape freeware (Inkscape Project, 2020). Polygons were then drawn on the traced images using ImageJ (Schneider et al., 2012) to select the groundmass area of interest, avoiding phenocrysts (Supplementary Material B). Whole microlites not touching the polygon sides were then counted and measured using the best-fitting ellipse tool in ImageJ. The smallest analyzed crystals were 0.1-0.9 μm in length, depending on the sample (Supplementary Material C).
Microlite areal number densities (N A , in mm −2 ) were calculated over the vesicle-free selected area. Areal crystal dimensions were measured as the short and long axes of best fit ellipses and the crystal areas of 2D intersections through the microlites were obtained. The plagioclase microlite crystallinity ϕ of the groundmass of each sample was calculated as the fraction of the selected groundmass area (on a vesicle-free basis) that is made up of plagioclase. The best-fit 3D short (S), intermediate (I) and long (L) dimensions of crystals were obtained by comparing distributions of the ratio of the 2D short/long axes of the microlites with corresponding distributions for intersections of known crystal aspect ratios using CSDslice (Morgan and Jerram, 2006). CSDCorrections software (Higgins, 2000) was used to perform stereological conversions on the 2D size distributions to obtain 3D crystal size distributions (CSDs). The ImageJ ellipse minor axis measurements were used for the calculations, as these represent the most likely intersections for crystals with prismatic morphologies (Higgins, 2000) as obtained from CSDslice. The area percentage of voids and micro-cracks were estimated with ImageJ and used as a correction for the area analyzed in CSDCorrections. The resulting three-dimensional data were binned into five size bins per decade and bins containing less than five crystals were removed (Higgins, 2000). The y axis intercept of a regression line fitted to the centers of the three smallest size bins on each CSD represents the plagioclase nuclei population density (n 0 , in mm −4 ) for each sample (Cashman, 1988), omitting any downturn of the CSD curve in small size bins resulting from incomplete compensation in the stereological calculations for the intersection probability effect for the smallest crystals (Brugger and Hammer, 2010a). The volumetric plagioclase number density (N V , in mm −3 ) was calculated as N V n 0 × (−1)/α and a characteristic crystal size (L c , in mm) was calculated as L c −1/α, where α is the slope of the regression line in the CSDs (Blundy and Cashman, 2008). The model of Klein et al. (2018), incorporating the effects of the crystal fraction and the particle aspect ratio, was used to calculate the maximum packing fraction (ϕ m ) for the groundmass of samples in this study.

Phenocryst Image Analysis
The areal percentage of phenocrysts (ϕ Ph ), including plagioclase, clinopyroxene, orthopyroxene and Fe-Ti oxides, on a vesicle-free basis, was assessed by image analysis using ImageJ on highresolution scans (4,800 dpi) of thin sections of each sample ( Figure 3). As phenocryst abundance appeared uniform in hand specimens of Galeras bombs, one representative thin section was analyzed per sample. For each sample, a large (300-557 mm 2 ) and small (59-165 mm 2 ) area of the thin section were selected for analysis (Supplementary Material C). The large/small analyzed areas varied slightly in size due to the necessity of choosing areas of the thin section with uniform thickness, as thinner areas caused variations in the color of certain phases in scanned images (e.g., groundmass glass, blue resin filling vesicles) that would impact image analysis. Each large/ small area was analyzed using two different methods to evaluate the robustness of the approach. In each method, the distinction between phenocrysts and groundmass was made on the basis of phenocrysts being distinguishable by eye in the high-resolution scans. The methods consisted of segmenting the different phases using HSB (Hue, Saturation, Brightness -Method 1) or RGB (Red, Green, Blue -Method 2) channel splitting on the selected image area followed by grayscale thresholding to select and quantify the phases of interest ( Figure 3). Air bubbles in vesicles that were erroneously segmented as one of the crystal phases were deleted, and holes in crystals were filled where appropriate before quantification.
Methods 1 and 2 yielded differences in ϕ Ph of ≤ 10% and ≤ 13% for the large and small areas of each sample, respectively. In addition, Method 1 yielded differences in ϕ Ph of 0-18% between the large and small areas, whereas Method 2 yielded differences of 1-10% for the same areas, with differences of ≤ 6% for all but one sample (Supplementary Material C). This suggests that some natural variability exists within samples and results depend to a limited extent on the choice of analyzed area within samples (e.g., whether the area is large enough to capture the largest phenocrysts), but also that Method 2 was more reproduceable than Method 1 for this suite of samples. We therefore selected the results of Method 2, and used the mean of results for the large and small areas to obtain ϕ Ph . We estimate a maximum error of ∼ 10% based on the results of Method 2 and note that our results more likely overestimate the phenocryst area due to a small amount of groundmass being erroneously attributed to the phenocryst fraction in the image analysis.

SO 2 Fluxes
Galeras volcano has been continuously monitored by the SGC since 1988 (Narváez Medina et al., 2017). The SGC provided SO 2 flux data covering the Aug. 2004-Jan. 2010 period of activity ( Figure 1D), which included 17 Vulcanian explosions and two episodes of lava dome effusion (Narváez Medina et al., 2017;Vargas and Torres, 2015). As the catalog of data is incomplete prior to the first explosion (11-12 Aug. 2004), and as a repose Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 time (defined as the time between a Vulcanian explosion and the previous one) for this explosion is difficult to define, the time period prior to that date was discounted from this analysis. The remaining SO 2 time series (Supplementary Material D) was used to obtain a first-order estimate of the magnitude of the gas flux from the shallow conduit during the repose period prior to each explosion. The SGC used a COSPEC (Correlation Spectrometer) during traverses beneath the wind-blown gas plumes to estimate SO 2 fluxes until July 27, 2006. Thereafter, SO 2 fluxes were measured using a DOAS (Differential Optical Absorption Spectrometry) instrument ( Figure 1D), typically by scanning the plume from a single vantage point, but traverses were also occasionally undertaken. Where several flux measurements existed in the time series for a single day, the measurements for that day were averaged. The mean of the ten largest daily fluxes over the repose period prior to a Vulcanian explosion was then taken as a first-order proxy for the magnitude of the SO 2 flux during that time (Supplementary Material D). It should be noted FIGURE 4 | Results of the material properties tests on Galeras samples. (A) Stress-strain curves calculated from the high-temperature rheometry tests. Each sharp inflection corresponds to the application of a higher strain rate (_ ε 10 -6 , 10 -5 , 10 -4 s −1 ). Note that the drop in stress at the transition from 10 -5 to 10 -4 s −1 for GAL6 results from slight unloading of the sample prior to the instrument imparting the final strain rate. (B) Apparent viscosity (η a ) calculated at the three different strain rates (_ ε) imparted during the rheometry tests. Dashed lines correspond to the best fit power law relationship between strain rate and apparent viscosity for each sample. The black line corresponds to the non-Newtonian rheological law for high-crystallinity dome lavas by Lavallée et al. (2007). (C) Permeability (κ) and connected porosity (ϕ c ) of the samples measured before (circles) and after (triangles) the high-temperature experiments. Dashed lines highlight the changes in these properties due to deformation. (D) Changes in permeability (where Δκ final κ − initial κ) as a result of deformation, compared with the initial connected porosity. (E). Apparent viscosity measured at the three imparted strain rates, compared with the changes in connected porosity (where Δϕ c final ϕ c − initial ϕ c ) as a result of deformation.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 8 that the scanning DOAS instrument is installed to the west of the main crater and measurements are dependent on environmental conditions such as the presence of daylight, cloud cover, wind direction and wind speed, therefore some SO 2 fluxes may not be recorded by the instrument as a result of these limitations. In addition, during 2004In addition, during -2006 to the implementation of the DOAS instrument, SO 2 flux measurements were not recorded every day ( Figure 1D), therefore the average SO 2 flux results for 2004-2006 reflect a lower temporal resolution than for 2007-2010.

High-Temperature Rheometry Tests
The mechanical data from the rheometry tests show that each sample responded in a similar fashion to deformation ( Figure 4A). Upon the application of the first constant strain rate ( _ ε 10 -6 s −1 ), stress increased non-linearly with strain until samples reached a peak stress that characterized the shear conditions. Upon the application of higher strain rates (_ ε 10 -5 s −1 , followed by _ ε 10 -4 s −1 ), stress again accumulated non-linearly until a new steady-state was achieved. All samples reached higher stresses upon application of higher strain rates ( Figure 4A). Each sample showed a slightly different magnitude of stress-strain response to the imparted strain rates, highlighting a range of rheological behaviors under similar deformation conditions. Although the samples commonly achieved a steady-state, some samples eventually exhibited strain weakening after a period of deformation at the highest applied strain rates, characterized by decreasing shear resistance as strain accumulated (e.g., GAL4; Figure 4A). The apparent viscosity computed for each applied condition decreased linearly as a function of applied strain rate ( Figure 4B), and this characteristic shear-thinning behavior was observed for all samples. Data from all tests are provided in Supplementary Material E.

Porosity, Permeability and Volume Changes
Prior to deformation, the samples showed a wide range of permeabilities (κ 1.72 × 10 -17 -3.39 × 10 -13 m 2 ), over the measured connected porosity span (ϕ c 2-26%). Measurements of these properties following deformation show that the connected porosity and permeability of most samples increased ( Figure 4C). Only one sample (GAL16) showed a small decrease in connected porosity and two samples (GAL6 and GAL8) showed decreases in permeability ( Figure 4C). For the samples that showed a permeability increase, the change (Δκ final κ − initial κ) was greater for samples with higher initial porosity ( Figure 4D). The samples also experienced a change in connected porosity (Δϕ c final ϕ c − initial ϕ c ) as they deformed from cylinders to barrels ( Figure 4E) proportional to the change in sample volume (Supplementary Material A). The samples that displayed the highest apparent viscosities (η a 0.5-1.2 × 10 12 Pa s at a strain rate of 10 -6 s −1 ) displayed the greatest changes in connected porosity (Δϕ c 5-20%; Figure 4E), whereas samples with lower apparent viscosities (η a 0.5-3 × 10 11 Pa s at a strain rate of 10 -6 s −1 ) showed more modest changes in connected porosity (Δϕ c <5%; Figure 4E).

Plagioclase Micro-Textures
The samples exhibit a range of micro-textural parameters in 2D (N A and mean crystal area) and 3D (n 0 , N V , L c , and S/L) (Figures 5, 6, Supplementary Material C). The natural logarithm of n 0 , lnn 0 , varies between 17.6 and 21.6 and represents a useful metric to track variations in other microtextural parameters as it is proportional to crystal nucleation rate and inversely proportional to growth rate (Cashman, 1988). To a first order, changes in ln n 0 therefore track changes in ΔT in this sample set ( Figure 5). Accordingly, samples with higher n 0 feature a higher plagioclase microlite number density, as measured through N A and N V , but a lower mean crystal area and characteristic crystal size (L c ) ( Figures 6A-D). In addition, the 3D aspect ratio of plagioclase microlites from CSDslice (Morgan and Jerram, 2006), measured through S/L, ranges between 0.1 and 0.2 for samples with low n 0 (ln n 0 17-20), indicating a prismatic morphology ( Figures 6F,G). Samples with higher n 0 (ln n 0 > 20) are characterized by S/L ranging between 0.3 and 0.5, indicating a more tabular morphology. The groundmass crystallinity (ϕ), does not vary systematically with n 0 ( Figure 6E).

Comparison Between Plagioclase Micro-textures and Viscosity
A comparison of the micro-textural results with the measured apparent viscosities shows that the data display some scatter (Figure 7), as expected from natural samples with complex viscosity dependencies. However, samples characterized by higher ln n 0 (and therefore higher N A , N V and S/L, and lower mean crystal area and L c ) generally display a lower apparent viscosity at a given strain rate compared to samples with lower n 0 (Figures 7A-E,G). Samples with a lower groundmass crystallinity, ϕ, also display a lower apparent viscosity ( Figure 7F). In addition, samples with higher ln n 0 typically have a higher groundmass maximum packing fraction (ϕ m ; Figure 8A). For each sample, comparing the observed groundmass crystallinity to the theoretical groundmass maximum packing fraction, ϕ/ϕ m , demonstrates that samples with a higher apparent viscosity display a crystallinity that is closer to ϕ m ( Figure 8B), consistent with the interpretation that variations in crystal micro-textures impose constraints on the observed viscosity spectrum.

Relationships Between Phenocrysts, Melt Viscosity, Porosity and Apparent Viscosity
The samples range in phenocryst content from 44% to 61%, and the data do not show a clear correlation between the phenocryst content and the apparent viscosity measured at a given strain rate ( Figure 9A). All groundmass glasses are rhyolitic in composition, with 73.8-78.5 wt% SiO 2 and NBO/T (non-bridging oxygen to tetrahedra ratio; calculated using the method described in Mysen Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 9 and Richet, 2019) in the range 0.01-0.06 (analytical results and melt viscosity calculations are provided in Supplementary Material F). The apparent and relative viscosities, η a and η r , display a linear relationship ( Figure 9B), demonstrating that the observed differences in apparent viscosity between samples (at a given strain rate) is more likely controlled by physical attributes of the multiphase suspensions than differences in interstitial melt composition.
Comparing the initial connected porosity of the samples with the measured apparent viscosity ( Figure 9C) indicates two groups of samples with similar trends. Samples with low connected porosities (ϕ c 0-20%) displayed higher viscosities, and showed both increasing apparent viscosity (at a given strain rate) and increasing strain rate-dependence of viscosity with increasing porosity. Samples with higher connected porosities (ϕ c > 20%) displayed lower apparent viscosities ( Figure 9C) and FIGURE 5 | Representative high magnification (×1,500) BSE (back-scattered electron) images of the groundmass of the selected bomb samples from Galeras volcano, organized by ascending plagioclase microlite ln n 0 . In these images, groundmass glass and plagioclase crystals appear in similar shades of dark gray, vesicles and microcracks appear black, pyroxene crystals appear light gray and Fe-Ti oxide crystals appear white. The complete analyzed groundmass areas are provided along with plagioclase microlite tracings in Supplementary Material B.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 10 also showed increasing apparent viscosity (at a given strain rate) with increasing porosity, but without an increasing dependence of viscosity on strain rate.

Viscosity Model Analysis
The relationship between stress (σ) and the resulting strain rate for each sample during the high-temperature deformation  Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 12 experiments can be described by the following power law relationship (Ostwald, 1925): where k (in Pa s n ) and n are the Ostwald and non-Newtonian constants, respectively. k is sometimes referred to as the consistency and is cognate with viscosity (Mueller et al., 2011a), whereas n is a measure of the deviation from linearity of the stress-strain rate relationship (Newtonian fluids have n 1, indicating a linear increase in strain rate with increasing stress; shear-thinning fluids have n <1, indicating a power law increase in strain rate with increasing stress). As the samples all display shear-thinning behavior (High-Temperature Rheometry Tests), the best-fit power law relationship between stress and strain rate can be used to estimate n and k for each sample (Supplementary Material G). By then fitting logarithmic relationships to the variation in n and k with n 0 ( Figure 10A,B, respectively), we found that both n and k share relationships with n 0 . Thus, n 0 may be used to parameterize the apparent viscosity of the tested Galeras samples with quantified plagioclase micro-textural properties: n 0.077 ln n 0 − 0.8 and k 3.894 × 10 9 ln n 0 − 6.479 × 10 10 [4] The obtained relationship between n and n 0 has a root mean square error (RMSE) of 0.091 and the relationship between k and n 0 has an RMSE of 3.756 × 10 9 Pa s n . These relationships assume that the effects of the concurrent variations in microlite textural characteristics, crystal size distribution, porosity, melt composition and phenocrysts are taken into account by the fitting parameters, and provide a simple empirical model to estimate the strain-rate dependent flow curves of Galeras samples with known n 0 .
To model viscosity for samples from the 2004 to 2010 explosions with known eruption dates, we used ln n 0 from the micro-textural dataset of Bain et al. (2019a) and applied Eqs. 3, 4 to obtain n and k for each sample. By then applying Eq. 2 to calculate the stress at three different strain rates (_ ε 10 -6 , 10 -5 , and 10 -4 s −1 ) and assuming that η a σ/ _ ε (Caricchi et al., 2007), we then modeled the strain rate-dependent apparent viscosity of the time-constrained samples ( Figure 10C). These calculations show that samples with high n 0 (ln n 0 ∼ 21) that erupted in 2009-2010 have the lowest modeled apparent viscosities (1-3 × 10 11 Pa s over the modeled range of strain rates, 10 -6 to 10 -4 s −1 , chosen to reflect the experimental conditions, which are appropriate for magma deformation in the shallow conduit; Figure 10C). In contrast, time-constrained samples with low n 0 FIGURE 9 | Effect of the phenocryst fraction, melt viscosity and initial porosity on the apparent viscosity. Apparent viscosity (η a ) of each sample, measured at each applied strain rate (see legend), is compared with: (A) the phenocryst content (ϕ ph ) obtained from image analysis, (B) the calculated relative viscosity (η r ), and (C) the initial connected porosity (ϕ c ).
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 (ln n 0 ∼ 17) that erupted in 2006 have higher modeled apparent viscosities (2 × 10 11 -2 × 10 12 Pa s), and the predicted viscosity difference between samples with low and high n 0 is especially marked at low strain rates (e.g., 10 -6 s −1 ; Figure 10C). In addition, these trends persist when taking into account the error on the relationships in Eqs. 3, 4 ( Figure 10C and Supplementary Material G). For example, negative errors on n and positive errors on k lead to estimated viscosities that are 0.5-1 log units higher than the modeled viscosities discussed above (upper limit of gray area in Figure 10C), with the greatest error on the samples with low ln n 0 at low strain rates. We note however that this error leads to apparent viscosities of up to >10 13 Pa s, which is ∼ 1 order of magnitude greater than the highest viscosities measured in our high-temperature experiments. The error in the model toward higher viscosities is therefore unlikely to be so high, and the model likely performs well in estimating the high viscosities expected at low strain rates at Galeras. Conversely, positive errors on n and negative errors on k lead to estimated viscosities that are 0.5-1 log units lower than the modeled viscosities (lower limit of gray area in Figure 10C), with the greatest error on samples with low ln n 0 at the highest strain rates (the low viscosity error on the sample with the lowest ln n 0 was not quantifiable, see Supplementary Material G). We note that the model more poorly reproduces the lowest viscosities measured in our hightemperature experiments at the highest strain rates (i.e., estimates a lowest modeled viscosity of ∼ 10 11 Pa s compared to a lowest measured viscosity of ∼ 10 10 Pa s in our experiments). The model may therefore overestimate viscosities by up to ∼ 1 order of magnitude at the highest strain rates. Interestingly, the slight inflection in viscosity at high strain rates that appears in our modeled data at a value of ln n 0 ∼ 18.3 ( Figure 10C) is more marked and occurs at ln n 0 ∼ 19 when considering the maximum error on the lowest viscosities. This inflection could imply more complex dynamics for samples with very low ln n 0 at high strain rates and merits further investigation that is outside the scope of this study. Overall, the error in the model does not significantly impact our modeled viscosity results and tends to compound the observed variations between magma viscosity in 2004-2006 compared to 2009-2010. We also compared k and n for the samples tested in this study with the observed connected porosity changes to investigate how rheology is related to the changes that took place during deformation. The high-viscosity samples that displayed the greatest shear-thinning behavior (i.e., n ≤ 0.6) and a low consistency (k) showed the greatest increases in ϕ c (Figure 11). In contrast, the lower-viscosity samples that displayed less shear-thinning behavior (n > 0.7) showed more modest increases in ϕ c . Thus, the sample rheology impacted the evolution of connected porosity during deformation.

Comparison of Modeled Apparent Viscosities With Monitoring Data
As the rheological properties of magma control the evolution of densification and, thus, permeability in the upper conduit, we now compare the monitored changes in explosion characteristics, lava dome growth and gas emissions in 2004-2010 with the apparent viscosities modeled in the previous section. We choose to compare gas emissions rather than permeability with viscosity, as permeability is highly variable, scale dependent, and importantly, transient in nature; yet, sustained gas emissions are congruent with a magma that retains permeability and does not rapidly yield to external compressive stresses, which is dictated by viscosity (cf. Ashwell and Kendrick et al., 2015). In 2004-2006, relatively small-volume explosions (0.08-1.2 × 10 6 m 3 ; Vargas and Torres, 2015; Narváez Medina et al., 2017) occurred at Galeras, ejecting samples with the highest modeled viscosities (Figure 12). These explosions were preceded by high average SO 2 fluxes ( ∼ 7,000-9,000 tons/days), and lava dome growth was observed for a short period of time between 13 January and end of March 2006 (no lava dome was present in the crater prior to this date, yet two Vulcanian explosions occurred; Figure 12). In contrast, in 2009, frequent and comparatively higher volume explosions (0.063-3.56 × 10 6 m 3 ; Vargas and Torres, 2015;Narváez Medina et al., 2017) ejected samples with the lowest modeled viscosities (Figure 12). These explosions were preceded by considerably lower average SO 2 fluxes ( < 500 tons/day) and destroyed a pre-existing lava dome.   Figure 10C (η a , shown at near-static conditions of _ ε 10 -6 s −1 to highlight the effect of the varying properties of the crystal cargo on the macroscopic rheology without considering additional effects related to strain rate; purple stars). The "inefficient densification" phase occurred when magma with a high η a occupied the shallow conduit, gas escaped comparatively freely and dome-building occurred. A transition phase occurred when magma with a lower η a entered the shallow conduit and gas escape became more restricted. Finally, an "efficient densification" phase occurred when magma with the lowest η a promoted rapid plug formation with low SO 2 fluxes from the conduit. The sequence of Vulcanian explosions ceased when the magma supply from depth dwindled.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 In 2008, a transition period between these two end-member phases occurred, with low average SO 2 fluxes ( ∼ 1,000 tons/ day) and a single Vulcanian explosion that ejected samples with an intermediate modeled viscosity, which was followed by the extrusion of a lava dome (Figure 12). The dataset therefore suggests that the explosion periodicity scales with magma viscosity (lower frequency, higher viscosity), which is also inversely related to the erupted volume (lower frequency, smaller volume). As mentioned in SO 2 Fluxes, the SO 2 flux monitoring data from 2004 to 2006 reflect a lower temporal resolution than for 2007-2010. However, we judge that this difference does not impact our conclusions regarding outgassing patterns over time from the Galeras conduit for the following reasons. 1) Despite the difference in resolution, very large fluxes up to 18,622 tons/ day were recorded in 2004-2006, whereas the largest flux recorded in 2009 was 9,300 tons/day, with most measurements recording much lower fluxes ( Figure 1D). The conduit was therefore clearly more open for degassing in 2004-2006 compared to 2009. 2) Here we use the average of the ten largest SO 2 fluxes in the repose period prior to a Vulcanian explosion as a first-order measure of the maximum extent to which the conduit was open to gas flux. In reality permeability in the conduit is transient and cycles of decreasing SO 2 flux prior to an explosion, followed by high flux immediately following an explosion, are well-documented at Galeras (e.g., Fischer et al., 1994;Cortés and Raigosa, 1997). Our chosen measure therefore captures only the average maximum flux during a given repose period, to provide an estimate of the maximal magma permeability prior to densification and plug-formation. We note that taking the single largest flux, or the mean of the largest five or ten flux measurements, does not change the observed pattern of high emissions in 2004-2006and lower emissions in 2008.

DISCUSSION
We now summarize our interpretations and reconcile our multiparametric dataset to propose a unifying model for the observed eruption dynamics at Galeras in 2004-2010.

Relationship Between Plagioclase Micro-Textural Characteristics and Rheology
The range of plagioclase micro-textural parameters in our sample set suggests crystallization over a range of effective undercooling, ΔT. The samples cover broadly the same range in plagioclase microlite textural characteristics and show the same relationships as other analyzed dense and scoriaceous samples from the 2004 to 2010 period of activity of Galeras ( Figure 6; Bain et al., 2019a). The samples therefore adequately reflect the typical range of micro-textural characteristics observed in erupted products from this period. The parameter n 0 is proportional to crystal nucleation rate and inversely proportional to growth rate (Cashman, 1988), and therefore represents a useful metric to assess variations in ΔT in this sample set, as a proxy for undercooling and, by extension, decompression rate. As N A , N V , mean crystal area, L C and S/L vary systematically with n 0 (Figure 6), n 0 also represents a useful parameter for tracking multivariate changes in plagioclase microlite textural characteristics. The groundmass crystallinity is also expected to be an important textural control on magma rheology, but ϕ does not vary systematically with n 0 as both high growth rates characteristic of low ΔT and high nucleation rates characteristic of higher ΔT can increase ϕ. The plagioclase nuclei population density n 0 is therefore considered to be the best metric to track both the variations in undercooling driven by changes in decompression rate, and the resulting changes in the plagioclase microlite characteristics.
The samples exhibited a range of apparent viscosities at a given strain rate, which correlate with the observed variations in microtextural characteristics. Samples with micro-textures characteristic of crystallization under high degrees of ΔT (i.e., high n 0 ) exhibited lower apparent viscosities at a given strain rate than samples that crystallized under lower degrees of ΔT (Figure 7). In addition, the samples exhibited a range of viscosity relationships with strain rate. The observed behavior was compared with the empirical shear-thinning rheological model of Lavallée et al. (2007) for high-crystallinity dome lavas (black line in Figure 4B): where T is the temperature (in°C). The samples displaying higher viscosities showed a viscosity magnitude and shear-thinning behavior that agree well with the model (Figure 4B), which was developed to resolve the strain-rate dependence of apparent viscosity for four chemically contrasting dome lavas, as they displayed similar behavior and a narrow apparent viscosity range ( ∼ 0.5 log units at a given strain rate; Lavallée et al. (2007)). In our dataset however, some samples displayed viscosities that are lower than the ± 0.25 log unit described in Lavallée et al. (2007), as well as showing a weaker dependence of apparent viscosity on strain rate (lower curves in Figure 4B, high n in Figure 10A). Together, these observations suggest that the samples that crystallized under lower degrees of ΔT have a rheology consistent with material that routinely erupts effusively in the form of lava domes at andesitic arc volcanoes. In contrast, samples that crystallized under higher degrees of ΔT have a rheology that cannot be adequately described using relationships established from previous experiments on dome rocks. This indicates that chemically similar magmas evolve physically, and thus rheologically, in contrasting manners during ascent at variable rates, which should be expected to impact eruptive activity.

Material Evolution During High-Temperature Deformation
The samples can be divided into two rheological groups that responded differently to deformation. The group of samples with the highest apparent viscosities (GAL4, GAL6, GAL7) experienced the largest increases in connected porosity, i.e., through dilation, as a result of the imparted strain ( Figure 4E). The group with lower viscosities showed only modest dilation (GAL5, GAL8, GAL14, GAL18), and in one instance, a small reduction in connected porosity, i.e., compaction (GAL16). Changes in permeability as a result of deformation appear to be related to the initial porosity, with more porous samples experiencing larger increases in permeability ( Figure 4D). Samples GAL6 and GAL8 are outliers that displayed a significant permeability decrease despite recording a slight (GAL8) to moderate (GAL6) porosity increase following high-temperature deformation. As a through-going micro-fracture was identified in the highresolution scan of the GAL8 thin section, we interpret these observations as the result of pre-existing micro-fractures that enhanced the original permeability of these samples by up to two orders of magnitude compared to samples with a similar connected porosity (Bain et al., 2019b). Such pre-existing micro-fractures are likely to have healed during the hightemperature deformation (e.g., Lamur et al., 2019), which could cause the observed decrease in permeability with little impact on porosity.
Our viscosity modeling also revealed that the samples that showed the greatest dilation (increases in ϕ c ) displayed the most shear-thinning behavior ( Figure 11B). Together, these observations suggest that the samples did not deform in a purely viscous regime, and brittle behavior occurred, enhanced by the presence of rigid crystals (Kendrick et al., 2013) and vesicles that act to concentrate stresses in bubble walls (Coats et al., 2018; Figure 4D). The greater dilation in samples with higher viscosity suggests that these samples experienced more extensive micro-fracturing, and these samples were notably more friable once cooled, following the high-temperature experiments. We therefore suggest that a transition in rheological behaviorfrom deformation characterized by viscous flow with small amounts of micro-cracking and modest dilation, to highly non-Newtonian deformation characterized by extensive microcracking and large increases in connected porosityoccurs in the range n 0.6-0.72 (for the range of strain rates and strains investigated in this study). Using Eq. 4, this corresponds to a transition at a nuclei population density of ln n 0 ∼ 18.20-19.76.

Rheology of the Magmatic Plugs at Galeras Volcano
The rheological constraints on plug-and dome-forming Galeras magma obtained in this study reveal viscosity variations over more than one order of magnitude at a given strain rate. However, there is no systematic control of the phenocryst fraction on the apparent viscosity in our dataset ( Figure 9A). Viscosity differences imposed by compositional variations of the rhyolitic interstitial melt phase in these samples also cannot explain the observed viscosity span ( Figure 9B). There is a possible weak relationship between the sample porosity and apparent viscosity ( Figure 9C), suggesting that a higher initial porosity may slightly increase viscosity, possibly as a result of the creation of isolated porosity during porous network compaction.
However, monitoring the changes in the connected porosity of the samples as a result of deformation reveals two rheological groups with the following characteristics: 1) GAL4, GAL6, GAL7 high viscosity, highly shear-thinning, high dilation; 2) GAL5, GAL8, GAL14, GAL16, GAL18lower viscosity, less shearthinning, lower dilation. These groups are also related to variations in the plagioclase micro-textures: 1) GAL4, GAL6, GAL7low n 0 , N A and N V , high mean crystal area and L c , prismatic morphology (low S/L); 2) GAL5, GAL8, GAL14, GAL16, GAL18higher n 0 , N A and N V , lower mean crystal area and L c , more tabular morphology (higher S/L, with the exception of GAL8 and GAL18). The correspondence of these textural and rheological groups is striking and suggests that the observed apparent viscosity differences are predominantly related to variations in the textural characteristics of the plagioclase microlites.
Our observations therefore allow the macroscopic rheological behavior of plug-and dome-forming magma to be linked with the micro-textural characteristics of the dominant groundmass phase produced by decompression-driven degassing. Bain et al. (2019a) found that dense and scoriaceous Galeras samples that crystallized under high ΔT host high numbers of small, tabular microlites (eg lnn 0 20.85-20.95, N V 3.07-3.34 × 10 6 mm −3 , L c 3 μm, and S/L 0.45-0.48 in the Feb. 20, 2009 and Jan. 2, 2010 eruptions) resulting from rapid magma ascent rates that imparted an average decompression rate of ∼ 10 MPa/ h. Samples that crystallized under lower ΔT host lower numbers of large, prismatic microlites (e.g., lnn 0 17.39-18.17, N V 3.66-5.55 × 10 5 mm −3 , L c 7-10 μm, and S/L 0.13-0.22 in the July 12, 2006 eruption) resulting from more modest ascent rates that imparted an average decompression rate of ∼ 1 MPa/h (see further discussion of these decompression rates in A Unifying Conceptual Model for Vulcanian Explosion Dynamics at Galeras). Our results therefore suggest that samples that experienced higher decompression rates formed magma plugs with lower apparent viscosities as a direct result of those micro-textural characteristics. Conversely, samples that experienced lower decompression rates formed magma plugs with comparatively higher apparent viscosities, with a highly non-Newtonian behavior consistent with published relationships for highcrystallinity dome lavas  Figure 4B). These samples also deformed via greater proportions of microcracking, demonstrating more dilatant behavior. Systematic variations in plagioclase microlite number density, size and aspect ratio resulting from differences in decompression rate may therefore exert an important control on the rheology of intermediate composition magmas, by changing the nature and frequency of particle-melt and particle-particle interactions (Mueller et al., 2011b;Klein et al., 2017;Klein et al., 2018).

Influence of Magma Rheology on Eruption Dynamics in 2004-2010
Acquiring a number of large time-constrained samples for analyses and material properties testing is an enduring challenge in volcanology, which motivates studies crosscorrelating relationships across broader suites of samples (cf. Harnett et al., 2019;Wallace et al., 2020), such as in the present study. In the case of Galeras volcano, a limited, highly valuable suite of time-constrained samples from the SGC collection was analyzed in previous studies (Bain et al., 2019a;Bain et al., 2019b), and here we used additional samples from the same period for this rheological investigation. By estimating the apparent viscosity of time-constrained samples ( Figure 10C) based on existing micro-textural data (Bain et al., 2019a), variations in magma rheology in 2004-2010 can be linked with observations made by the SGC during the volcanic crisis.
The andesitic magma erupted at Galeras in 2004-2010 was degassed and highly crystalline (Bain et al., 2019a), consistent with extensive decompression-driven crystallization. When such magma stalls in a shallow volcanic conduit, gas bubbles are essentially immobilized in the crystal-rich slurry due to its high viscosity (Stix et al., 1997) and outgassing occurs through the creation of permeable pathways formed by connected bubbles and cracks (e.g., Westrich and Eichelberger, 1994;Edmonds et al., 2003). However, despite the overall very high viscosity of the magma erupted at Galeras in 2004-2010, our results indicate important second-order viscosity variations related to a spectrum of "low"-( ∼ 10 11 Pa s at near-static conditions of _ ε 10 − 6 s −1 ) to "high"-viscosity ( ∼ 10 12.5 Pa s at _ ε 10 − 6 s −1 ) behavior, which are predominantly caused by differences in the magma's petrological properties and are linked with variations in eruption style. In 2009-2010, the emplacement of comparatively low-viscosity magma plugs was associated with frequent, largevolume explosions preceded by low average SO 2 fluxes from the conduit ( Figure 12). As low-viscosity magma should facilitate deformation and promote compaction, we suggest that efficient densification restricted outgassing pathways, rapidly leading to the development of critically high pore pressures capable of driving Vulcanian explosions (cf. Spieler et al., 2004). An example of this behavior was the explosion that occurred on Jan. 2, 2010 ( Figure 12). This explosion was preceded by a short repose time (43 days) and a low average SO 2 flux (440 tons/day), suggesting that the conduit pressurization timescale was short due to the existence of a low-permeability plug. In contrast, the emplacement of higher viscosity magma plugs during 2004-2006 was associated with infrequent, low-volume explosions preceded by higher average SO 2 fluxes, suggesting inefficient densification and comparatively open outgassing pathways. For example, the explosion that occurred on July 12, 2006 partially destroyed the lava dome that appeared in Jan. 2006 ( Figure 12). This explosion was preceded by a long repose time (230 days) and a high average SO 2 flux (6,690 tons/day), suggesting that magma densification was sluggish, magma permeability was high and outgassing pathways were open, leading to long timescales for pressurization of the upper conduit.

A UNIFYING CONCEPTUAL MODEL FOR VULCANIAN EXPLOSION DYNAMICS AT GALERAS
The multi-parametric findings presented in this study advocate for explicit relationships between magma decompression rate, its physical and thus rheological evolution, outgassing, eruption style and explosion recurrence rate. Based on a comparison of Galeras crystal micro-textures (N A and ϕ) with the results of decompression experiments on hydrous rhyodacite by Brugger and Hammer (2010b), Bain et al. (2019a) estimated that average decompression rates evolved from ∼ 1 MPa/h at the onset of the eruptive phase to ∼ 10 MPa/h toward the end of 2009, though decompression rates at Galeras are generally uncertain due to scarce data. Although a more precise quantification of decompression rates is not yet available, we can reasonably consider that magma was decompressed at a higher average rate in 2009 than in 2004-2006, and we use the estimations of Bain et al. (2019a) to represent "low" ( ∼ 1 MPa/h) and "high" ( ∼ 10 MPa/h) decompression rates at Galeras. These decompression rates are low to moderate compared to decompression rates reported at other andesitic volcanoes straddling the effusive-explosive transition (Cassidy et al., 2018).
Magma plugs emplaced in the shallow conduit (<0.5 km) of Galeras following high average rates of decompression ( ∼ 10 MPa/h) feature a groundmass with large numbers of small, tabular plagioclase crystals (Bain et al., 2019a; Figure 13A). Despite the abundance of microlites in the groundmass, crystal-crystal interactions are comparatively low in these magma plugs, setting a high potential maximum packing fraction in the magmatic suspension and producing a comparatively low apparent viscosity (modeled at ∼ 10 11 Pa s at near-static conditions of _ ε 10 − 6 s −1 ). This low viscosity is compounded by the inherent moderately shear-thinning behavior of such magmas at the high strain rates expected from the inferred rapid magma ascent rates. This situation should promote rapid densification of the magma at shallow pressures, thus efficiently restricting permeable pathways, consistent with the low SO 2 fluxes monitored in 2009-2010 at Galeras. This is consistent with the study of Bain et al. (2019b), where samples with plagioclase micro-textures typical of high decompression rates were found to be associated with porous micro-textures characteristic of extensive densification and a lowpermeability magma plug ( ∼ 10 -16 m 2 ). Under conditions of high effective undercooling resulting from such high decompression rates, crystallization is likely to continue within and below this low-permeability plug, driving further magma degassing but with a restricted ability to outgas. The low permeability of this plug material is likely to promote rapid development of the pore overpressure needed to trigger Vulcanian explosions (e.g., Mueller et al., 2008;Lavallée and Kendrick, 2020). Moreover, such dense plug material is likely to confer a high tensile strength to the magma (Hornby et al., 2019), requiring high pore overpressures to develop through gas accumulation below the plug to bring it to failure. This efficient "plug-forming" regime is therefore characterized by short repose times between explosions, large stress drops upon failure of strong dense magma, and high volumes of ejected material ( Figure 13A). Consequently, during periods of high magma decompression and ascent rates, Vulcanian explosions are likely to become frequent (repose times of tens of days) and large ( ∼ 10 6 m 3 ). These explosions are likely to be preceded by low measured SO 2 fluxes ( ∼ 1,000 tons/day). In this regime, hazards associated with Vulcanian blasts, namely ballistic ejection and ash plumes, should be expected. As these explosions are likely to be large in magnitude, they may result in extensive destabilization of more vesicular magma residing deeper in the conduit and empty the conduit to significant depths, potentially feeding large eruption columns and column-collapse pyroclastic flows at the larger end of the explosivity spectrum. This latter scenario has not occurred at Galeras since the 1988 reactivation (there is evidence of pyroclastic flows in historic activity, most recently in 1936, but the mechanism of generation of these flows is unclear; Banks et al., 1997), however the pumice flows related to Vulcanian explosions at Soufrière Hills Volcano (Druitt et al., 2002) may represent an example of this deeper magma destabilization.
In contrast, magma plugs emplaced following low average rates of decompression ( ∼ 1 MPa/h) feature a groundmass with low numbers of large, prismatic crystals (Bain et al., 2019a; Figure 13B). The phenocryst fraction is also likely to be higher in this magma, due to high crystal growth rates (Cashman, 2020) (this was not observed in our dataset, however we note that sample GAL5 features the highest microlite number density and the lowest phenocryst content, supporting the expectation that the phenocryst fraction will be lower in more rapidly decompressed samples). As the total crystallinity (microlites and phenocrysts) is likely to be high relative to the maximum packing fraction (at least compared to the plug materials formed during rapid ascent), slowly decompressed magma is likely to exhibit a comparatively high apparent viscosity (modeled at ∼10 12 Pa s at near-static conditions of 10 -6 s −1 ). Crystal-crystal interactions are also comparatively high in these magma plugs, which therefore have a lower ϕ m (we estimate that the plagioclase crystal sizes and shapes produced under such decompression conditions at Galeras reduced ϕ m by ∼5%; Figure 8A). The viscosity contrast between magmas with these textural attributes and those encountered in plug materials having undergone rapid decompression is on the scale of an order of magnitude, implying that deformation timescales will be approximately ten times slower. Thus, on eruptive timescales, this high viscosity should result in inefficient magma densification in the shallow conduit, accompanied by a higher probability of deformation through micro-cracking due to the dilatant behavior identified in such samples. High-viscosity magmas are more likely to experience brittle failure at a given strain rate (e.g., Dingwell and Webb, 1989;Lavallée et al., 2008;Coats et al., 2018;Wadsworth et al., 2018) and in shallow magmatic environments, where large pore pressures can develop and external stresses are low (i.e., magmastatic and lithostatic loads), brittle behavior is most commonly dilatant, causing the creation of permeable pathways (e.g., Kendrick et al., 2013;FIGURE 13 | Summary diagrams showing the relationships between average magma decompression rate, petrological and rheological evolution, eruptive behavior and timescales. The sketch demonstrates how (A) "efficient densification" and (B) "inefficient densification" end-member rheological regimes are dictated by magma decompression rate and nucleation/crystallization kinetics (illustrated by the insets of plagioclase texture (in black)). This petrological evolution impacts magma rheology as (A) high crystal number densities (and other accompanying micro-textural attributes) favor comparatively lower viscosities, promoting densification and a low permeability that inhibits outgassing, causing explosions at an increasing recurrence rate. In contrast (B), slow decompression and ascent favor the crystallization of fewer larger crystals that promote magmas with high permeability and viscosity, conducive to protracted dome emplacement and extensive outgassing, which alleviate the periodicity of explosive events.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 611320 Lavallée et al., 2013;Lavallée and Kendrick, 2020). In this regime, sluggish densification and increased micro-crack density should therefore promote gas fluxing through the shallow conduit, consistent with the high SO 2 fluxes monitored in 2004-2006 at Galeras. This is consistent with the results of Bain et al. (2019b), where samples with crystal micro-textures typical of low decompression rates were found to be associated with porous micro-textures that had not undergone significant densification and were characteristic of a comparatively high-permeability plug ( ∼ 10 -12 m 2 ). Under the low degrees of ΔT imparted by low decompression rates, crystallization can continue within and below the high-permeability plug, driving further magma degassing, with abundant outgassing pathways that result in longer time periods for the necessary overpressure to build up and trigger a Vulcanian explosion. The high porosity of the magma plug indicates a comparatively lower tensile strength (Hornby et al., 2019), requiring lower pore overpressures to trigger the rupture of the plug. This regime of activity ( Figure 13B) is therefore characterized by comparatively long repose times between explosions (hundreds of days), small stress drops due to the poor capacity to build pressure, and hence smaller ejected volumes ( ∼ 10 5 m 3 ), preceded by high measured SO 2 fluxes ( ∼ 10,000 tons/day). Furthermore, the low strain rates associated with the inferred low ascent rates should promote a high apparent viscosity response in the magma characteristic of slowly-emplaced plugs ( Figure 4B). As the rate-dependence of viscosity of these magmas is stronger (i.e., a lower n than for plugs formed by rapid decompression, indicating highly shear-thinning behavior), strain localization near conduit margins should be enhanced, increasing the likelihood of plug flow and a switch to effusive behavior promoting exogenous dome growth (e.g., Hale and Wadge, 2008). This may explain why rheological laws developed for high-crystallinity dome lavas  also appear well-suited to describing the behavior of magma plugs emplaced following low magma decompression rates. Magma deformation in these high strain zones is also more likely to undergo crystal interactions and brittle fracturing Lavallée et al., 2008;Kendrick et al., 2017;Wallace et al., 2019), further enhancing outgassing and delaying the development of critical overpressure (e.g., Kendrick et al., 2013;Lavallée et al., 2013). The observed changes in connected porosity ( Figures 11A,B) suggest that the most efficient outgassing conditions that might favor dome extrusion due to extensive micro-cracking occur when magma with a low nuclei population density (i.e., ln(n 0 ) below the identified transition at 18.20-19.76) is emplaced in the upper conduit, which agrees well with the range of n 0 for samples produced in explosions associated with the presence of lava domes at Galeras volcano (Bain et al., 2019a). This conceptual model ought to be tested against larger datasets, and expanded across different magma compositions (i.e., across the spectrum of plug-and dome-forming magmas from high-crystallinity basaltic andesite, to andesite and dacite) before it can be applied in a predictive fashion. However, we surmise that if this model is more broadly applicable, andesitic lava domes erupted during periods of cyclical Vulcanian explosions should generally be associated with small volume explosions and long repose times, except in cases where explosions are driven by sudden decompression resulting from lava dome collapse (e.g., Voight and Elsworth, 2000) or when high rates of secondary mineral precipitation reduce the lava dome permeability (e.g., Horwell et al., 2013;Heap et al., 2019).
Periods of low average magma decompression rate are therefore likely to be associated with lava dome extrusion and small-volume Vulcanian explosions, if the permeability of the dome and magma in the shallow conduit is insufficient to regulate pore pressures developing in deeper magma. In this "domeforming" regime, Vulcanian blasts are likely to affect a comparatively restricted area around the vent due to their small magnitude, however hazards associated with dome collapse should be anticipated, namely block-and-ash flows and explosions driven by sudden decompression of the conduit. We note that dome collapse has not historically occurred at Galeras due to the growth occurring within the restricted crater of the summit cone.
In nature, we expect a spectrum of behavior ranging between the defined end-member regimes to occur, such as in 2008 at Galeras, as a result of variations in magma decompression rates and styles. Our study does not address the causes of the variations in magma ascent rates, however these could result from buoyancy differences, perhaps as a result of tapping a magma reservoir that was stratified in terms of temperature or volatile content, or from variations in magma chamber overpressure over time (e.g., Melnik and Sparks, 1999). Such differences in temperature and volatile content would undoubtedly influence the magma rheology, (e.g., Giordano et al., 2008), especially as magma evolves physically (via crystallization, vesiculation and vesicle compaction) and chemically (via degassing, and crystallization) upon ascent under evolving P-T-strain conditions (impacted by ascent, shear heating, latent heat of crystallization and outgassing, etc.). Here, we do not attempt to resolve all of these processes but assess the rheological properties of magma emplaced in the shallow conduit. As we experimentally demonstrate, petrological variations drive differences in apparent viscosity of Galeras magmas of over an order of magnitude for a given strain rate. However, other processes could also impact this rheological behavior. Magma ascending at a slower rate would in principle undergo less viscous heating and have more time to lose heat to wall rocks, whereas magma ascending more rapidly may be expected to retain more heat when it reaches the shallow conduit. Slight temperature differences associated with these different ascent pathways may be expected to act as a positive feedback effect, with higher temperatures contributing to lower the viscosity of the "low"-viscosity, rapidly decompressed andesitic magma and lower temperatures contributing to increasing the viscosity of slowly decompressed magma. Furthermore, in terms of the dense magma plugs considered here, the groundmass glass of all samples from 2004 to 2010 is highly degassed (<0.36 wt% H 2 O; Sample Selection). However, Bain et al. (2019b) found evidence of a slightly greater extent of degassing (i.e., a lower H 2 O concentration of ∼0.1 wt%) in porous samples that were interpreted to have undergone more sluggish densification, as a result of pathways for degassing remaining open for longer. Conversely, dense samples having undergone extensive densification displayed a slightly higher H 2 O concentration in the groundmass glass (∼0.3 wt%) due to the rapid disruption of degassing pathways during efficient densification. This slight difference in melt volatile content should act to increase the magma viscosity in high-viscosity (slowly decompressed) plugs and decrease the magma viscosity in low-viscosity (rapidly decompressed) plugs. The positive feedbacks related to temperature and volatiles are therefore likely to compound the differences in eruptive style described above.
This conceptual model ties together observations from monitoring data and the physical, textural and rheological characteristics of erupted products at Galeras volcano, revealing a link between magma decompression rate, degassing and crystallization processes, magma rheology, plug densification, outgassing, overpressure development, explosion recurrence rate and eruption style. As such, the monitoring of crystal micro-textures in ash during and between explosions in volcanic systems where these relationships have been constrained could provide an additional line of evidence for interpreting volcanic behavior (e.g., Wright et al., 2012;Wallace et al., 2020), for eruption forecasting and for hazard mitigation during prolonged volcanic crises, when combined with geophysical (e.g., seismic and ground deformation data) and geochemical (e.g., SO 2 flux data) monitoring. If crystal micro-textures in ash produced in gas-and-ash venting events at persistently active arc volcanoes (e.g., Wright et al., 2012;Hornby et al., 2018;Murch and Cole, 2019) could be monitored in near-real time (e.g., Gaunt et al., 2016), this could provide valuable information concerning magma ascent rates, rheology, permeability, fragmentation mechanisms and the range of hazards to be anticipated during eruptions.

CONCLUSION
Micro-textural, physical and rheological investigations of ballistic bomb samples from the 2004-2010 Vulcanian explosions of Galeras volcano suggest that high magma decompression rates promote efficient densification due to a comparatively low apparent viscosity arising from the micro-textural results of crystallization under high degrees of undercooling. In turn, efficient densification fosters the development of a dense, lowpermeability magma plug impeding outgassing. This "plugforming" regime promotes frequent, large explosions preceded by relatively low average SO 2 fluxes. Conversely, low decompression rates lead to inefficient densification due to a higher apparent viscosity and stronger non-Newtonian behavior arising from crystallization under low degrees of undercooling. This "dome-forming" regime is characterized by high average SO 2 fluxes and dome extrusion accompanied by small-volume explosions with long repose times. We suggest that this integrated textural-physical-rheological conceptual model should be rigorously tested against other datasets and time-constrained eruptive products at well-monitored volcanoes, as it may improve near real-time assessment of hazards during volcanic unrest, which commonly feature variable discharge, decompression and ascent rates. We advance that existing monitoring techniques ought to be routinely supplemented by systematic sample collection and detailed micro-textural analyses of juvenile ash particles to aid in the interpretation of volcanic behavior during protracted periods of unrest at andesitic arc volcanoes.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
AB, EC, JC, and GC conceptualised the research project. EC, JC, and GC supervised the research. AB, JK, AL, and YL performed the physical and material properties tests. DM and RT provided monitoring data from Galeras volcano. AB performed the microtextural analysis and wrote the article. All authors provided comments on the article.

ACKNOWLEDGMENTS
AB thanks Viviana Burbano, Joao Lages, Bertilda Botina and Carlos Estrada for fieldwork assistance at Galeras. The authors thank Mike Hall, John Craven and Chris Hayward for technical assistance at the University of Edinburgh. We also thank the reviewers for their thoughtful comments that helped to improve this manuscript, as well as Fabio Arzilli and Valerio Acocella for editorial assistance. This paper is dedicated to Viviana Burbano, who is sadly missed.