Modeling decompression paths in a basaltic andesite magma using the nucleation and growth of plagioclase microlites

Plagioclase microlites in a magma nucleate and grow in response to melt supersaturation (Δϕplag). The resultant frozen plagioclase crystal size distribution (CSD) preserves the history of decompression pathways (dP/dt). SNGPlag is a numerical model that calculates the equilibrium composition of a decompressing magma and nucleates and grows plagioclase in response to an imposed Δϕplag. Here, we test a new version of SNGPlag calibrated for use with basaltic andesite magmas and model dP/dt for the ca. 12.6 ka Curacautín eruption of Llaima volcano, Chile. Instantaneous nucleation (Nplag) and growth (Gplag) rates of plagioclase were computed using the experimental results of Shea and Hammer (J Volcanol Geotherm Res 260:127–145, 10.1016/j.jvolgeores.2013.04.018, 2013) and used for SNGPlag modeling of basaltic andesite composition. Maximum Nplag of 6.1 × 105 cm h−1 is achieved at a Δϕplag of 44% and the maximum Gplag of 27.4 μm h−1 is achieved at a Δϕplag of 29%. Our modeled log dP/dtavg range from 2.69 ± 0.09 to 6.89 ± 0.96 MPa h−1 (1σ) with an average duration of decompression from 0.87 ± 0.25 to 16.13 ± 0.29 h assuming a starting pressure Pi of 110–150 MPa. These rates are similar to those derived from mafic decompression experiments for other explosive eruptions. Using assumptions for lithostatic pressure gradients (dP/dz), we calculate ascent rates of < 1–6 m s−1. We conducted a second set of Monte Carlo simulations using Pi of 15–30 MPa to investigate the influence of shallower decompression, resulting in log dP/dtavg from 2.86 ± 0.49 to 6.00 ± 0.86 MPa h−1. The dP/dt modeled here is two orders of magnitude lower than those calculated by Valdivia et al. (Bull Volcanol, 10.1007/s00445-021-01514-8, 2022) for the same eruption using a bubble number density meter, and suggests homogeneous nucleation raises dP/dt by orders of magnitude in the shallow conduit. Our modeling further supports the rapid-ascent hypothesis for driving highly explosive mafic eruptions.


Investigating magma ascent rates
Decompression rate affects eruption style (e.g., Eichelberger et al. 1986;Jaupart and Allegre 1991;Burgisser and Gardner 2005;La Spina et al. 2021;Bamber et al. 2020). As magmas ascend from depth, volatiles exsolve and crystals nucleate and grow in response to changes in pressure (P) and temperature (T). During rapid ascent, bubbles remain coupled to the magma resulting in explosive eruption (Eichelberger et al. 1986;Jaupart and Allegre 1991). Conversely, during slow ascent, bubbles coalesce, resulting in sufficient permeability to degas the melt and thus removes the volatile primer necessary for explosivity and results in effusive eruption (Mangan and Sisson 2000). Crystallization of microlites during Communicated by Mark S Ghiorso. decompression increases magma viscosity (Vona et al. 2011;La Spina et al. 2016;Vetere et al. 2021) and may act to either impede the ability of gas to decouple from the magma or enhance coalescence by pushing isolated vesicles together. As such, understanding the rate of magma decompression and therefore ascent rate is important for estimating eruption duration, intensity, and volcano hazards.
Several analytical and experimental methods exist for the investigation of magma decompression rate, each with their own strengths and weaknesses. Bubble and crystal textures provide a record of magma decompression or ascent path (Cashman and Marsh 1988;Blundy and Cashman 2008;Arzilli et al. 2019;Bamber et al. 2020;Marshall et al. 2022a;Valdivia et al. 2022), and thus rocks provide a valuable look into the subsurface evolution of a magma. Crystal size distributions (CSDs) of microlites can be used to approximate crystallization times when a crystal growth rate is assumed (Marsh 1988;Cashman and Marsh 1988;Murch and Cole 2019;Bamber et al. 2020;Valdivia et al. 2022). Although CSDs can be easily measured and their slopes used for interpretation of changing ascent rates, the calculations may be skewed if post-fragmentation crystallization occurs. In addition, CSDs assume constant crystal nucleation and growth rates. More sophisticated investigations involve reproducing measured microlite textures by performing magma decompression experiments (Fig. 1), during which crystal textures evolve in response to an applied perturbation in P and/or T (Geschwind and Rutherford 1995;Hammer and Rutherford 2002;Hammer 2004;Szramek et al. 2006;Castro and Dingwell 2009;Andrews and Gardner 2010;Brugger and Hammer 2010;Shea and Hammer 2013;Waters et al. 2015;Befus and Andrews 2018). Decompression experiments are effective at approximating ascent rates by producing sufficient effective undercooling (ΔT eff ) or supersaturation ( Δ ) necessary to drive crystallization, but for the most part only produce time-averaged ascent rates that do not reflect possible changes in ascent rate as a magma nears the surface. Furthermore, conducting decompression experiments can be time-consuming. Mineral breakdown reaction rims (Rutherford and Hill 1993;Browne and Gardner 2006) and compositional zoning (Waters et al. 2015) form in response to the pressure change imposed on a magma during ascent but are not always present. Melt embayments allow for diffusive modeling of elemental loss and thus ascent rates (Liu et al. 2007;Myers et al. 2016Myers et al. , 2018Barth et al. 2019). Melt inclusions and embayments are, however, not perfect storage containers. Mineral fractures may result in leakage, and diffusion modeling cannot be conducted without knowledge of initial conditions and diffusive boundary conditions. Finally, geophysical observations can be used to monitor seismicity with depth in real time and allows researchers to track magma movement during an eruption (e.g., Moran et al. 2008;Thelen et al. 2008). Most volcanoes, however, are not equipped with extensive geophysical arrays that allow precision monitoring, and geophysical observations may not distinguish between different types of subsurface volcanic activity.

Existing numerical models for magma ascent rate
To circumvent some of the disadvantages of existing experimental and analytical methods for investigating ascent rates, numerical models exist that utilize observations easily collected from rocks. Toramaru (2006) developed a magma ascent rate meter as a function of bubble number density (BND) assuming a single homogeneous nucleation event and constant decompression. Although BNDs indeed reflect changes in volatile supersaturation and decompression, extensive coalescence, multiple nucleation events, highly tortuous bubble networks (e.g., Valdivia et al. 2022), or collapsed foam textures are not representative of original BNDs and will skew ascent rate calculations. The model of Toramaru et al. (2008) uses microlite number densities (MND) to estimate ascent rates and only requires water and groundmass Si content at the point of microlite nucleation as additional inputs. But as Murch and Cole (2019) point out, the model results of Toramaru et al. (2008)  influenced by the Si content input, and an error of only 5% in Si content can result in errors in ascent rate calculations as large as 500%. In addition, both models only produce time-averaged ascent rates rather than instantaneous rates over time, and therefore do not adequately model variable ascent rates that occur in nature (e.g., Mastin and Ghiorso 2000;Moran et al. 2008;Thelen et al. 2008). More recent numerical models consider the fluid dynamics of three phase basaltic magmas and disequilibrium crystallization and vesiculation during their ascent (La Spina et al. 2016;Arzilli et al. 2019;La Spina et al. 2021;Bamber et al. 2020).

SNGPlag
Supersaturation Nucleation and Growth of Plagioclase (SNGPlag) is an iterative forward model that calculates time-dependent plagioclase crystallization, the integral of nucleation and growth, within a constant magma composition for a specified pressure-temperature-time (P-T-t) path (Andrews and Befus 2020). Comprehensive descriptions of the model can be found in Befus and Andrews (2018) and Andrews and Befus (2020) and are only summarized here. Specifically, the model tracks the numbers and sizes of plagioclase crystals within a 1 m 3 volume of magma. SNG-Plag considers nucleation and growth as functions of plagioclase supersaturation (Δϕ plag ), defined as the difference between the equilibrium volume fraction of plagioclase as determined using MELTS (Gualda et al. 2012;Ghiorso and Gualda 2015) and the modeled volume fraction. SNGPlag uses Δϕ plag rather than ΔT eff as the former can be readily determined through time whereas ΔT eff is only known at the onset of decompression. Melt decompression and/or cooling act to increase Δϕ plag . Nucleation and growth of plagioclase crystals in response to Δϕ plag drive the magma towards equilibrium, with the instantaneous nucleation and growth rates of plagioclase being functions of Δϕ plag (Befus and Andrews 2018). SNGPlag allows nucleation and growth to be path-dependent and does not assume constant nucleation and growth rates (Andrews and Befus 2020). SNGPlag can model multiple styles of decompression (e.g., linear, accelerating, paused) to investigate the style of decompression on plagioclase crystallization. Multiple decompression styles may be applied to the same simulation, such as a linear pathway that has a pause during decompression. While SNGPlag cannot provide a unique solution for natural samples, it can describe a limited range of likely decompression rates and paths (Andrews and Befus 2020).
Previous versions of SNGPlag are calibrated for felsic compositions. Here, we extend the calibration of SNGPlag to include basaltic andesite compositions using the experimental results of Shea and Hammer (2013). We then apply an inverse implementation of SNGPlag to the 12.6 ka basaltic andesite Curacautín eruption of Llaima volcano, Chile (Marshall et al. 2022a;Valdivia et al. 2022) to estimate decompression rates necessary to generate ignimbrite-forming mafic eruptions. The results and application of our modeling can be applied to similar mafic volcanic centers to investigate the conditions that result in unusually explosive mafic eruptions.

The Curacautín eruption
The Curacautín eruption occurred at ca. 12.6 ka and resulted in the deposition of the extensive Curacautín ignimbrite (Ci, Marshall et al. 2022a). The Ci is a 4.0-4.5 km 3 (dense-rock equivalent) unconsolidated basaltic andesite ignimbrite exposed radially around Llaima that flowed up to 30 km from Llaima (Marshall et al. 2022a;Naranjo and Moreno 2005), though others have mapped the Ci up to 100 km from source (Naranjo and Moreno 1991). The Ci consists of four coarse ash to fine lapilli tuff flow units (Fig. 2, Marshall et al. 2022a). Recent work by Marshall et al. (2022a) and Valdivia et al. (2022) suggests the Ci is the result of fragmentation of a rapidly ascending, non-degassed magma at a low fragmentation threshold. There is no evidence to suggest the explosivity of the Ci eruption was driven by magma-water interaction, though some evidence exists for localized phreatic activity (Marshall et al. 2022a, b).

Calibration of SNGPlag for basaltic andesite compositions
Previously published versions of SNGPlag (Befus and Andrews 2018; Andrews and Befus 2020) use nucleation and growth rates determined experimentally for the 1991 Pinatubo dacite magma with a rhyolitic melt composition. Application of SNGPlag to the Curacautín eruption necessitates acquiring plagioclase nucleation (N plag ) and growth rates (G plag ) for a basaltic andesite magma. We used the results of single step decompression experiments conducted by Shea and Hammer (2013) on the Mascota basaltic andesite. Their study includes 11 experimental runs (Table 1) with P, T, H 2 O, and compositional conditions reasonable for the Curacautín eruption (Lohmar 2008;Schindlbeck et al. 2014). Importantly, they report the plagioclase crystallinities, maximum lengths, and volumetric number densities for all runs, thereby enabling calculation of instantaneous nucleation and growth rates.

Determination of instantaneous nucleation and growth rates of plagioclase
We adapted the existing SNGPlag code written in MATLAB to find N plag and G plag that best fit the experimental observations of Shea and Hammer (2013). Briefly, we assume that the N plag and G plag have functional forms that can be described as log-normal functions of Δϕ plag ; variation of four different parameters can change the functional form to virtually any arbitrary path (Befus and Andrews 2018). We find the best fit for N plag and G plag by running SNGPlag for the known decompression experiments of Shea and Hammer (2013) across an 8-dimensional space (four dimensions for both N plag and G plag ). This results in 100,000 possible combinations of N plag and G plag . N plag and G plag were modeled using the R2 high performance computing cluster at Boise State University. The best fit N plag and G plag are those that best recover the observed results of Shea and Hammer (2013). Run parameters were taken from Shea and Hammer (2013) with each single-step run discretized into 2,500 P-T-t steps. N plag and G plag are calculated at each step as functions of Δϕ plag with the form: where x = b·exp 1 Δϕ plag , Δϕ plag = plagioclase supersaturation, and μ, σ, b, and k are fit parameters that describe the specific shape of curves that represent the mean, standard deviation, scaling with respect to Δϕ plag , and its maximum value (Befus and Andrews 2018). The input ranges and best fit calibration parameters for N plag and G plag are provided in Table 2. Values for μ, σ, and k were randomly sampled from a selected range (Table 2). For our calibration, b was set to 1. SNGPlag accounts for volume interferences ϕ int between crystals for a randomly distributed population of crystals by, where ϕ app is the apparent crystallinity, which is the sum of all crystal sizes and numbers calculated at each step divided by the system volume (1 m 3 ). From this, we obtain the equation, where plagioclase crystallinity ϕ plag is reported with overlapping crystals removed (Andrews and Befus 2020). Finally, uncertainty in N V and σ Nv is determined by, where Sn is the characteristic crystal size in a 1 mm 2 area (Andrews and Befus 2020). Optimum values for G plag and N plag were determined using least squares optimization of the calibration data (Table 2). During each step of SNGPlag, existing plagioclase grow, and new plagioclase nucleate based upon G plag , N plag , and Δϕ plag . SNGPlag produces matrices of plagioclase crystal number and size that can be binned and converted into cumulative CSDs. Because SNGPlag calculates volumetric number densities and size distributions by nucleating and growing plagioclase in a 1-m 3 model volume, we avoid the uncertainties that result from stereological conversions of 2D data. Use of a 1-m 3 model volume to effectively eliminate rounding errors and discrepancies that can occur in smaller volumes with fewer crystals.

Modeling conditions
Modeling the Ci CSDs using the best fit N plag and G plag rates requires realistic or plausible values for P i , P f , T, dP/dt, and volume fraction phenocrysts. Schindlbeck et al. (2014) calculated Ci crystallization temperatures of ~ 1,110 ± 45 °C using the olivine-and clinopyroxene-liquid thermobarometer of Putirka (2008), water content of 1.4 ± 0.32% using the plagioclase hygrometer of Lange et al. (2009), and storage pressures between 400 and 600 MPa corresponding to depths of up to 18 km, though work by Lohmar (2008) suggests that crystallization occurred at ≤ 7 km. Marshall et al. (2022a) measured phenocryst content of Curacautín pyroclasts from < 1 to ~ 3.5%, and Lohmar (2008) measured up to 7% phenocrysts. Valdivia et al. (2022) estimated dP/dt for the Ci from 0.36 to 2.6 MPa s −1 using the bubble number density decompression rate meter of Toramaru (2006). Finally, experiments by Arzilli et al. (2019) found the conditions required for basaltic magmas to erupt as high explosivity events are temperatures < 1100 °C, syn-eruptive crystal content ≥ 30%, and melt viscosities of 10 5 Pa s. Our modeling consisted of 100,000 simulations with initial and final conditions selected in a random Monte Carlo scheme from a range of defined inputs (Table 3). We conducted experiments with P i between 110 and 150 MPa due to the limited range of experimental pressure data used for SNGPlag calibration. We did not assume a 400-600 MPa P i because those pressures calculated by Schindlbeck et al. (2014) are total pressure whereas those of SNGPlag are the partial pressure of water (e.g., P H2O ). Starting phenocryst content was 5 vol.%. P f was set to 10-40 MPa. All simulations were run at T = 950-1050 °C; T i and T f were allowed to vary independently. We used average dP/dt of 1-1000 MPa h −1 (0.0003-0.3 MPa s −1 ). 40% of runs were linear decompressions, 30% accelerating, and 30% were twostep decompressions, whereby there was a pause following initial linear decompression and subsequent post-pause decompression was either linear or accelerating. A subset of experiments was declared to "fragment" at a pressure P frag of 20-80 MPa during the simulations; these runs had dP/ dt of 1-20 MPa h −1 prior to fragmentation and increased to 30-400 MPa h −1 following fragmentation. Runs that fragmented experienced cooling ΔT frag of up to 60 °C, the upper bound suggested by Mastin and Ghiorso (2001) for adiabatic cooling of an erupting mixture of gas and ash. 10 9 -10 13 6.0677 × 10 9 10 -10 -10 -5 2.2003 × 10 -8

Comparison of natural and modeled CSDs
Cumulative CSDs of natural samples describe the number of plagioclase crystals that are larger than each size bin. Using counting statistics, we can convert that size relationship into an uncertainty bound (σ CSD ) at each size, CSD = √ n bin , where n bin is the number of microlite counts per size bin. The upper and lower bounds then define an envelope for natural CSDs (Fig. 3). Therefore, with higher n bin , our uncertainty becomes smaller. Our modeled CSDs therefore have an effective uncertainty of zero as the billions of crystals compose each bin. This is not to say the modeling here is perfect, but rather that uncertainty is orders of magnitude greater in measurements of the natural samples.

Modeling limitations
The experiments of Shea and Hammer (2013) were mostly quenched at higher pressures, with only two experiments decompressed to P f of 22 and 10 MPa and ΔT > 113 °C (Table 1). Those two experiments produced the highest plagioclase crystallinities of 34.8% and 46.1%, respectively. No experiments have been conducted at conditions where the melt viscosity should be highest. As such, our N plag and G plag for very high Δϕ plag are extrapolated, although we note that any decompression path other than single-step will have some crystallization prior to reaching lower P, and therefore have a lower Δϕ plag than a single step run initially has at the same pressure. SNGPlag does not consider any unique conduit geometries or eruption style (e.g., dike geometry, ring faulting during eruption) that may impact late decompression or ascent dynamics. Shearing along conduit margins is not considered in SNGPlag but has been shown to impact crystallization (Vetere et al. 2021). Nucleation delay inversely correlates with ΔT and may suppress crystallization up to ~ 100 + h at low ΔT (Rusiecka et al. 2020), which could impact CSDs. SNGPlag does not explicitly consider nucleation delay in the crystallization calculations, but low growth rates acting upon small numbers will result in very low numbers of detectable crystals at early times or low ΔT. Finally, the only volatile species considered in our modeling is H 2 O, although the presence of CO 2 or another volatile species should only affect the crystallization of plagioclase insofar as it reduces the partial pressure of H 2 O.

Instantaneous nucleation and growth rates of plagioclase
Instantaneous N plag and G plag curves have similar geometries (Fig. 4). The maximum N plag of 6.1 × 10 5 cm −3 h −1 is reached at Δϕ plag = 44 vol.%. The maximum G plag of 27.4 µm h −1 is reached at Δϕ plag = 29 vol.%. There is very little N plag activity at Δϕ plag < 10%, but the G plag of these early crystals is quite high. N plag and G plag beyond maximum Δϕ plag are extrapolated and may not be representative of nature.

Model results
The large parameter space over which we modeled the Curacautín eruption includes many runs that are physically unrealistic; we applied filters to remove those results. Our filters identified runs that begin and/or end > 10 °C above the plagioclase liquidus and removes them. This reduced the number of model runs from 100,000 to 13,283 (Table 4). Because our decompression rates vary in an exponential fashion, it is not appropriate to compare them in linear space, so we report our average decompression rates as log 2 values. For example, three decompression rates of 1, 9, and 80 MPa h −1 would yield a linear average rate of 30 MPa h −1 , but a more representative average rate is obtained in log space and yields 9 MPa h −1 .
Average Unit 1 dP/dt are 53-93 MPa h −1 for L1, 46-89 MPa h −1 for L4, and 87-95 MPa h −1 for L6 (Figs. 5,6). Average Unit 2 dP/dt are 62-93 MPa h −1 (L8). Average Unit 3 dP/dt are the slowest at 6-55 MPa h −1 (L10). Conversely, average Unit 4 dP/dt are the fastest at 104-141 MPa h −1 (L18 , Table 4). Unit 1 average durations of decompression t avg are between 1.40 and 4.08 h for L1, Fig. 3 Example of how uncertainty is shown on our crystal size distribution (CSD) model runs (Supplemental Materials). The blue line is the natural cumulative CSD and the pink lines are the 2σ error bounds calculated for each bin. Notice how 2σ decreases with smaller microlite sizes. This is a result of the higher number of microlites counted in the natural samples at these size ranges. The increase in 2σ near the y-intercept (gray field) results from a relative decrease in the number of smallest crystals counted in 2D measurements of the natural sample (Fig. 4) (Valdivia et al. 2022) 2.40-4.69 h for L4, and 1.79-1.96 h for L6. Unit 2 t avg are between 1.69 and 2.02 h (L8). Unit 3 t avg are between 3.56 and 16.13 h (L10). Unit 4 t avg are between 0.87 and 0.96 h (L18, Figs. 5, 6, Table 4).

Discussion
Plagioclase nucleation and growth rates N plag and G plag curves (Fig. 4) for the basaltic andesite Curacautín magma have similar shapes, but very different magnitudes in comparison to those determined for the 1991 Pinatubo dacite (Befus and Andrews 2018). The Curacautín magma reaches a maximum N plag = 6.1 × 10 5 cm −3 h −1 at Δϕ plag = 44 vol.% which is an order of magnitude lower than the Pinatubo dacite at the same Δϕ plag (Fig. 4A). Conversely, the maximum Curacautín G plag of 27.4 µm h −1 is reached at Δϕ plag = 29 vol.%, whereas the 1991 Pinatubo G plag for the same Δϕ plag is 6.0 µm h −1 and does not reach 27.4 µm h −1 until Δϕ plag ≲ 52 vol.% (Fig. 4B). Indeed, G plag is more than an order of magnitude higher in the mafic composition for Δϕ plag ≲ 5%. Our modeled N plag and G plag suggest that although plagioclase nucleates more than an order of magnitude slower in basaltic andesites than in dacites at similar Δϕ plag the growth rate G plag in the mafic composition is generally an order of magnitude faster. Significantly, the difference in volumetric growth rate is ~ 1000 times greater in the basaltic andesite (the linear growth rate G plag raised to the third power). That is, a smaller number of crystals are able to more rapidly grow and thus reduce Δϕ plag in the mafic magma as compared to more silicic magmas. This explains the predominance of acicular plagioclase microlites commonly observed in the pyroclasts of mafic explosive eruptions (Constantini et al. 2010;Arzilli et al. 2019;Bamber et al. 2020;Rowe et al. 2021;Marshall et al. 2022a).

Decompression rates
Natural plagioclase CSDs for the Ci are concave upward at the finest size bins (Valdivia et al. 2022). Valdivia et al. (2022) divided Ci CSDs into two segments based on linear regression fitting. Using experimentally derived growth rates of 10 -4 mm s −1 (Arzilli et al. 2019), 2 × 10 -5 mm s −1 (Arzilli et al. 2015), 10 -6 mm s −1 (Shea and Hammer 2013), and 10 -7 mm s −1 (Arzilli et al. 2015), they calculated timescales of crystallization from 2 s to 1.2 h for the smallest size fraction of plagioclase microlites in CSDs, and 8 s to 5.0 h for the largest size fraction. Here, we use cumulative natural CSDs for fitting to our modeled CSDs to remove downturns at the smallest size fractions observed by Valdivia et al. (2022).
Using the 1% population of isolated Ci vesicles, Valdivia et al. (2022) calculated average dP/dt for the Ci magma of 0.84-1.95 MPa s −1 for Unit 1, 0.36 MPa s −1 for Unit 2, 2.60 MPa s −1 for Unit 3, and 0.55 MPa s −1 for Unit 4 using the BND meter of Toramaru (2006), with a minimum average dP/dt for the Curacautín eruption of 1.4 MPa s −1 . Our average modeled dP/dt rates (0.18 × 10 -2 -3.9 × 10 -2 M Pa s −1 ) are approximately two orders of magnitude slower  Table 4 Summary of 150,000 SNGPlag results. Three natural crystal size distributions (CSDs) for each sample were modeled. For each model run, SNGPlag generates a series of fits to the natural CDSs, denoted below as CSD fit. Fit f1 corresponds to the model runs that best fit the natural CSDs bins within 2σ (out of 31 total bins). Fit f2 is the second best fit and is determined by removing one bin from the total bins that fit. than the rates calculated by Valdivia et al. (2022, Fig . 7, Table 4). The bubble textures investigated by Valdivia et al. (2022) represent two distinct phases of Curacautín magma evolution. The highly tortuous vesicle network of > 99% pore volume is indicative of relatively slow ascent (e.g., Marshall et al. 2022b), whereas the small, isolated vesicles likely formed during an episode of homogeneous nucleation very late in ascent or syn-eruptively at low pressures in the shallow subsurface (Mangan and Sisson 2000) where dP/ dt are greatest. Conversely, our average dP/dt modeled with SNGPlag represent pressures from 10 to 150 MPa where rates of decompression begin slow and increase over time (Supplemental Material). Together, our work and that of Valdivia et al. (2022), suggests that decompression (and therefore ascent) rates increase by orders of magnitude at the shallowest conduit depths (Fig. 7) or during fragmentation. Interestingly, the decompression rates modeled here are in agreement with those modeled by La Spina et al. (2021) for highly explosive eruptions, and the decompression rates calculated by Valdivia et al. (2022) are consistent with rates associated with high fountaining (La Spina et al. 2021) and may further support the boiling over eruption model proposed in Marshall et al. (2022a). The low water content of the Ci melt (1.1 ± 0.32%; Schindlbeck et al. 2014) suggests storage at shallower depths, or water undersaturation. We conducted a second, smaller set of experiments (n = 50,000) at P i = 15-30 MPa and P f = 3-10 MPa to investigate crystallization over a shorter decompression window from shallower depths (Table 4). Average dP/dt for Unit 1 are 24-59 MPa h −1 (L1), 34-46 MPa h −1 (L4), and 46-61 MPa h −1 (L6). Unit 2 dP/dt are 66-75 MPa h −1 . Unit 3 average dP/dt are 8-61 MPa h −1 . Finally, average Unit 4 dP/dt are 74-80 MPa h −1 . These rates tend to be slower than those modeled for deeper chamber conditions but are generally within the same order of magnitude (Table 4). Because Schindlbeck et al. (2014) estimated a chamber depth of ~ 18 km for the Curacautín magma, the dP/dt calculated with P i up to 150 MPa are likely a more reasonable approximation of Curacautín decompression (Fig. 7).
The dP/dt modeled here for the Curacautín magma are similar to dP/dt calculated or estimated for other mafic eruptions using decompression experiments and diffusion modeling, but are orders of magnitude lower than mafic dP/ dt calculated from bubble textures (Fig. 8). Homogeneous bubble nucleation events create densely packed networks of bubbles at very shallow depths where rates of dP/dt are highest (Mangan and Sisson 2000), and thus dP/dt determined from bubble textures may only reflect very shallow ascent conditions and not be representative of conditions from deeper in the conduit. Conversely, our modeling here reflects ascent rates integrated over the entire conduit and not just the shallowest depths and likely records more of the

Magma ascent rates
Decompression rates do not have the same relationship to ascent rate at all volcanoes. This results from differences in lithostatic or magmastatic pressure gradients at different volcanoes, which is impacted by factors such as crustal thickness, country rock compositions and densities, conduit geometry, and elevation. In addition, particular decompression speedometers may be sensitive to the partial pressure of a particular volatile species, not total pressure (P total ); SNGPlag is sensitive to P H2O , which is less than P total when the system is water undersaturated or saturated with a mixed volatile phase. Here, we consider two simplified scenarios to derive first order estimates of magma ascent rate from our modeled decompression rates, and then compare those rates with a calculated lithostatic pressure gradient (dP/dz) for the crust beneath Llaima.
Our first estimate assumes that P H2O = P total and that there is no other volatile species in our system. This of course is an oversimplification as there would be some amount of P CO2 present as well as others volatile species in minor concentrations. If we also assume that a dP/dz = 90 MPa per every 4 km is reasonable for a mix of mafic lavas and granitic plutons (Naranjo and Moreno 2005), then we obtain average Unit 1 ascent rates for the Ci of 0.66 ± 0.67-1.13 ± 0.78 m s −1 for L1, 0.66 ± 0.58-1.10 ± 0.86 m s −1 for L4, and 1.07 ± 0.80-1.17 ± 0.80 m s −1 for L6. Our Unit 2 (L8) average ascent rates are 0.77 ± 0.37-1.14 ± 0.81 m s −1 . Unit 3 (L10) average ascent rates are the slowest at 0.08 ± 0.01-0.68 ± 0.58 m s −1 . Conversely, Unit 4 (L18) average ascent rates are the fastest at 1.28 ± 0.09-1.74 ± 0.84 m s −1 . Due to our assumptions and simplifications, these rates should be considered a minimum (Fig. 9).
Our second calculation combines our modeling parameter space with a chamber depth estimate of 18 km (Schindlbeck et al. 2014). If we assume the Ci magma is water undersaturated, then we can expect the magma resided at a greater depth prior to decompression. Using a maximum P i during SNGPlag modeling of 120 MPa, we obtain an effective dP/dz in P H2O of 60 MPa per every 9 km. Using these new assumptions, our ascent rates for the Ci magma increase. Average ascent rates for Unit 1 become 2.23 ± 2.27-3.89 ± 2.04 m s −1 (L1), 1.90 ± 1.97-3.72 ± 2.90 m s −1 (L4), and 3.61 ± 2.68-3.96 ± 2.71 m s −1 (L6). Unit 2 average ascent rates are 2.58 ± 1.24-3.86 ± 2.73 m s −1 . Unit 3 average ascent rates are 0.27 ± 0.02-2.30 ± 1.97 m s −1 . Finally, Unit 4 average ascent rates are 4.31 ± 0.30-5.86 ± 2.85 m s −1 . Because this second set of ascent rates assumes the same decompression rates as our first scenario but over a greater depth, they should be considered maximum estimates (Fig. 9).

Fig. 7
Curacautín ignimbrite (Ci) decompression rates (dP/dt) modeled using SNGPlag plotted with respect to Ci stratigraphy (m) (Marshall et al. 2022a) along with the dP/dt calculated by Valdivia et al. (2022) from x-ray computed microtomography 3D renderings and using the bubble number density rate meter of Toramaru (2006). Sample names are provided in red and associated units are plotted along the right y-axis. SNGPlag curves are provided for all three crystal size distribution fits (see explanation in Table 3). dP/dt results plotted are those from the 100,000 model run dataset (Table 4)

Difficulty of fitting smallest CSD microlites
SNGPlag struggles to fit the smallest crystal sizes in the observed plagioclase CSDs. This may be a result of the tighter 2σ bounds at smaller sizes because the number of crystals exceeding those sizes is large, therefore decreasing the uncertainty allowed in the model fits (Fig. 3).
Alternatively, poor fits at small microlite sizes may result from the range of P values reported in the experiments of Shea and Hammer (2013) that we used for calibration of SNGPlag N plag and G plag rates (Table 1). The lowest P f used for calibration are 10 and 22 MPa, but these were only two out of the eleven experiments, whereas the other nine were conducted to 42 ≤ P f ≤ 100 MPa (Shea and Hammer 2013).

Fig. 8
The range of decompression rates (dP/dt) for mafic magmas estimated using different methods. Blue = decompression experiments. Green = diffusion modeling. Red = bubble number density (BND). Black = SNGPlag. SNGplag modeling has the most overlap with decompression experiments and diffusion modeling. The similarity of our modeled dP/dt to decompression experiments is likely due to the way SNGPlag is calibrated using Shea and Hammer (2013) data. dP/dt calculated using BND data are consistently orders of magnitude higher. This may be a function of using bubbles from homogeneous nucleation events which occur at very shallow depths and reflect moments of very high dP/dt (Mangan and Sisson 2000) Fig. 9 Curacautín magma ascent rates (m s −1 ) versus decompression rate in both log 2 dP/dt and dP/dt (MPa h −1 ). Polynomial fits to our minimum and maximum end-member estimates for lithostatic pressure gradient (dP/dz) and that of Schindlbeck et al. (2014) Fig. 6 Because N plag and G plag are not linear with respect to Δϕ plag (Fig. 4), they should be higher in experiments conducted at very low P, corresponding to high Δϕ plag , although this extrapolation does not account for much high viscosities potentially retarding the nucleation and growth rates. Our modeled P f and calibration P f stop at 10 MPa, but natural plagioclase textures likely continue to record shallower conduit conditions. In this scenario, we would expect crystallization of a higher number of smaller plagioclase microlites, which may have produced the densely crystalline Ci pyroclasts (Marshall et al. 2022a, b;Valdivia et al. 2022).

Interpreting the Curacautín eruption
Rapid magma ascent rates are often invoked to explain mafic Plinian and ignimbrite-forming eruptions (Sable et al. 2006(Sable et al. , 2009Vinkler et al. 2012;Arzilli et al. 2019;Bamber et al. 2020;Marshall et al. 2022a;Valdivia et al. 2022). Supersaturation resulting from rapid ascent drives plagioclase nucleation and crystallization. Our modeling here reveals that N plag in the basaltic andesite Ci is considerably lower than N plag in dacites, but maximum G plag of 7.6 × 10 -7 cm s −1 is volumetrically up to 1000× greater than dacite G plag at the same Δϕ plag . Our G plag is one order of magnitude lower than the ~ 3-5 × 10 -6 cm s −1 measured by Vetere et al. (2021) during basaltic andesite viscosity experiments. Those authors argue for the importance of shear rate being considered in models of magmatic and volcanic processes, which is not something considered in this version of SNGPlag (Table 2). Indeed, shear rate and its impact on viscosity would impact our G plag and may help explain conduit processes proposed by Marshall et al. (2022b).
Our modeling here suggests that rapid dP/dt produced the plagioclase microlite textures observed in Ci pyroclasts (Table 3; Marshall et al. 2022a, b;Valdivia et al. 2022). Such extensive crystallization would have increased the magma viscosity to the point that vesicles would begin to distort and wrap around the nucleating and rapidly growing acicular plagioclase. This explains the highly tortuous 99% interconnectivity vesicle population textures identified by Valdivia et al. (2022). Highly tortuous vesicle networks inhibit outgassing, which in turn enhances the overpressure necessary for brittle fragmentation.
The three sets of magma ascent rates we estimated here using different dP/dz reasonable for the South Central Volcanic Zone of Chile offer a first-order look into the ascent rates that drove the Curacautín eruption (Fig. 9). Minimum ascent rates of 0.1-1.7 m s −1 using a dP/dz of 22.5 MPa km −1 are similar to the ascent rates of 0.1-2.0 m s −1 we estimated using the dP/dz of Schindlbeck et al. (2014). Conversely, a dP/dz of 20 MPa per every 3 km yields ascent rates up to 3× faster (Fig. 9).
Unit 1 ascent rates are variable between 0.6 and 1.3 m s −1 and increase slightly to 0.8-1.3 m s −1 in Unit 2. Unit 3 ascent rates drop by an order of magnitude to 0.1-0.8 m s −1 and suggests modulation of the magma flux during the Curacautín eruption. Unit 4 has the fastest magma ascent rate of 1.3-2.0 m s −1 and represents the final pulse of the Ci eruption. Valdivia et al. (2022) calculated vesicle overpressures necessary to fragment the Ci magma between 3.8 and 5.1 MPa. Such a low fragmentation threshold combined with the rapid dP/dt calculated here implies a limited decompression history prior to climatic fragmentation. Because the Ci was produced during a single eruptive event (Marshall et al. 2022a), changes in magma ascent rate did not likely result from changes in shallow magma storage or magma recharge, but rather changes in vesiculation or conduit/vent geometry during eruption. Discriminating between those different parameters is beyond the scope of the current version of SNGPlag.
The rapid G plag calculated in our modeling would generate acicular plagioclase morphologies that produce highly tortuous vesicle networks that inhibit outgassing. Following fragmentation, decompression and ascent rates of the gaspyroclast mixture are orders of magnitude greater than the original bulk magma and suggest there is little time between fragmentation and eruption. In the case of the Ci, the time period between fragmentation and eruption likely generated the highly crystalline groundmass of l < 10 μm plagioclase microlites that overprints sutures between fused domains of heterogeneous vesicle textures discussed in Marshall et al. (2022b). These results help elucidate the still poorly understand conduit processes that impact how mafic magmas can erupt as large, explosive events.

Conclusions
Plagioclase nucleation and growth rates, N plag and G plag , respectively, differ substantially between mafic and felsic magmas. Those differences can affect eruption style. Modeled maximum N plag for the 12.6 ka basaltic andesite Curacautín eruption are orders of magnitude lower than those for the 1991 Pinatubo dacite (Fig. 4); G plag , however, is up to 10× greater in mafic magmas than felsic magmas, resulting in volumetric growth rates ~ 1000× greater in mafic magmas than felsic ones. This result explains the predominately acicular nature of plagioclase microlites in the products of mafic explosive eruptions attributed to rapid ascent rates.
The dP/dt modeled here for the Ci eruption using SNGPlag are between 10 -3 and 10 -1 MPa s −1 and are similar to dP/dt measured experimentally for similar compositions and known eruption styles (e.g., Arzilli et al. 2019;La Spina et al. 2021;Bamber et al. 2020). We were able to fit the majority of CSD bins to the natural samples. Unlike decompression experiments which must follow some particular decompression pathway(s) (Fig. 1), our modeling applies instantaneous N plag and G plag to produce thousands of possible decompression pathways to derive the most likely decompression scenario, and thus reflect the total decompression path of the Ci magma. Our modeled dP/dt are ~ 2 orders of magnitude lower than those calculated by Valdivia et al. (2022) for the same eruption. This difference reflects time-integrated rates recording most of magma decompression and ascent presented here, whereas those of Valdivia et al. (2022) were calculated using the BND meter of Toramaru (2006) on a homogenous nucleation event in the shallow conduit. Importantly, these two sets of dP/dt reveal that decompression (and therefore magma ascent) of the Curacautín magma increased by orders of magnitude following the onset of fragmentation and record the explosive nature of the eruption. In addition, such a dramatic change in ascent rate would have similar impacts on Δϕ plag (Fig. 4), resulting in the crystallization of the l < 10 μm population of unbroken plagioclase microlites identified by Marshall et al. (2022b) and may explain the rapid crystallization times Valdivia et al. (2022) calculated from plagioclase CSDs.
Future work is necessary to fully describe the effects of decompression on crystallization and eruption processes described here. Integrating a viscosity component into SNG-Plag would allow us to investigate viscosity's role on ascent dynamics, which has profound impacts on degassing and crystallization and may help explain the textures reported in Marshall et al. (2022b). Plagioclase is not the only crystal phase in many mafic eruptions, and future modeling should consider additional crystal phases such as pyroxenes and olivine in addition to plagioclase. Finally, future decompression experiments conducted to very low P i (and therefore higher melt viscosity) would enhance the calibration parameter space of SNGPlag and allow for the investigation of crystallization at the shallowest depths of conduits where microlites are likely to crystallize most extensively.