T he shape of volcanic conduits inferred from bubble size distributions

The most intense known explosive volcanic eruptions on Earth are Plinian eruptions of silicic magma. Geospeedometers indicate that Plinian magma erupts from high pressure within the magma chamber at average speeds of 0.001-1 MPa/s . At the same time, dissolved magmatic volatiles, predominantly water, nucleate about one quadrillion bubbles per cubic meter of melt, preserved as vesicles within tephra. Vesicles span several orders of magnitude in size, with power-law size distributions, and vesicularities of approximately 70±20 vol.% . To date these combined observations of relatively modest decompression rates, a reflection of the fluid dynamics of magma ascent, as well as the exceedingly large bubble number densities with power-law size distributions, a consequence of bubble nucleation and growth kinetics, have never been explained in a quantitatively coherent manner. Here we demonstrate that the integration of these observations requires that bubble nucleation commences when magma starts ascending from within the chamber and continues until the magma fragments to produce tephra, contrary to the conventional view that nucleation occurs as a burst at high supersaturation. We substantiate experimentally that nucleation in rhyolitic melt can occur continuously over prolonged time intervals. We then use integrated modeling of bubble nucleation and fluid dynamics of magma ascent to demonstrate that bubble size distributions in Plinian pyroclasts are the product of continuous nucleation throughout magma ascent, and at average decompression rates that are consistent with geospeedometers. An important outcome of our results is that the transition from magma chamber to volcanic conduit is gradual, resembling a cupola that narrows upward into a conduit.


Introduction
The abundant bubbles that form during magma ascent in the conduit power explosive volcanic eruptions. The supercritical aqueous fluid within bubbles, sometimes referred to as vapor or gas, is of lower density than the silicate melt from whence it exsolves [1,2]. By virtue of conservation of mass the volume of magma, that is silicate melt plus bubbles, increases as it rises toward the Earth's surface. This results in magma acceleration and expansion [3,4]. During explosive eruptions the associated strains and strain rates will cause the magma to fragment, which in turn is a defining process of explosive silicic eruptions [5,6,7].
Bubbles form by exsolution of dissolved magmatic volatiles as a consequence of decreasing magma pressure en route to the surface [8]. Volatile exsolution involves the simultaneous nucleation of new bubbles and diffusion from the melt into existing bubbles [9,10].Volatiles become supersaturated as magma pressure decreases, inducing bubble nucleation [8,10,11]. Diffusion of volatiles, in particular H 2 O, from the melt into existing bubbles causes bubbles to grow [9,10], while at the same time decreasing supersaturation and increasing magma viscosity [12]. A dynamical feedback between decompression, bubble nucleation and diffusion ensues [10,13]. A record of these intertwined processes is the bubble size distribution (BSD), preserved in pyroclasts as the vesicle size distribution (VSD) [14,15,16,17,18,19].
VSDs have been used to infer observationally inaccessible subsurface conditions during eruptive magma ascent [14,15,16,17,18,20,19,21,22,23]. To date there have been no quantitative models, integrating bubble nucleation and growth within a fluid dynamical model of magma ascent and decompression, which are capable of reproducing VSDs and thereby providing quantitative constraints on explosive eruptions. Moreover, exist-ing models that predict the overall number density of bubbles tend to infer magma decompression rates that are at odds with geospeedometers which yield orders of magnitude lower decompression rates [24,25]. Here we resolve both deficiencies using numerical modeling of eruptive magma ascent and vesiculation, with the latter encompassing bubble nucleation and growth. Our simulated BSDs reproduce observed VSDs at average decompression rates that are consistent with geospeedometers.

Continuous nucleation
The continuous power-law vesicle size distributions in pyroclasts are thought to be the result of protracted bubble nucleation throughout magma ascent [20,21]. Although such continuous nucleation is in principle consistent with nucleation theory [10,26,27], wide VSDs as observed in Plinian pyroclasts are usually not produced in nucleation experiments [28,29,30]. Moreover, because it is nearly impossible to directly observe bubble nucleation during these high pressure and temperature experiments, continuous nucleation cannot be validated through direct observation.
The objective of our experiments is to validate that bubble nucleation during decompression of hydrated rhyolitic melt is indeed continuous. The experiments consisted of three steps [31] (Supplementary Table 1; Supplementary Figure 1). First samples were hydrated at an initial pressure, P i = 160-250 MPa, and temperature, T = 850 • C, until equilibrium was achieved. Next, samples were decompressed over a time interval, 1 s ≤ t d ≤ 50 s, to their final pressure, P f = 10-100 MPa. In the final step, samples were held at the final pressure for a duration of annealing time, 0 s ≤ t post ≤ 120 s, and then were rapidly quenched. Bubble number density (BND) was subsequently measured in thin  sections of the quenched samples. The experiments were divided into six suites for each of which saturation pressure, decompression rate, and final pressure were closely similar. Within each suite samples were quenched at different t post . If nucleation was ongoing during the annealing period, then samples within a given suite should have a linear correlation between t post and bubble number density.

Modeling of magma ascent and bubble nucleation
We simulate the fluid dynamics of eruptive magma ascent together with bubble nucleation and growth during Plinian eruptions, Figure 1 (see Methods for details). We assume a cylindrical conduit whose radius may change with depth and solve the steady mass and momentum conservation equations for the change in magma velocity and pressure with depth. The latter provides the boundary condition for the integrated sub-grid bubble nucleation and growth model. We consider H 2 O the most abundant volatile phase [32] and driver of bubble nucleation [8,10] and assume heterogeneous nucleation [24,25]. Our objective is to test the hypothesis that VSDs in Plinian pyroclasts are the product of continuous nucleation during magma ascent. Each simulation is therefore constrained to maintain sufficient supersaturation for nucleation. This is achieved by varying conduit size, such that at any given depth the magma decompression rate is slightly greater than the inverse characteristic diffusion time. Conceptually this approach is similar that employed in conduit models with a near-lithostatic magma pressure [33]. Thus, aside from predicting decompression rate, bubble number density and size distribution, the model also predicts conduit radius as a function of depth.

Experimental verification of continuous bubble nucleation
Our experiments produced bubble number densities that vary over 5 orders of magnitude (10 8 -10 13 bubbles per m 3 of melt). Although supersaturation, defined as the difference between the initial and final pressures, drives nucleation, experimental BNDs are not correlated with supersaturation (Supplementary Figure  2). Instead, experiments within a suite produced BNDs that vary by up to 2 orders of magnitude and are correlated with the annealing time, t post .
Because BNDs among various suites vary significantly, a synthesis of all experiments requires non-dimensionalization of bubble number density and annealing time. t post is scaled by the nucleation time scale, τ, which represents the time that nucleation rate drops by one e-folding time. N post , which represents BND for bubbles that are nucleated during annealing time, is non-dimensionalized by the equilibrium bubble number density, N eq , the final bubble number density that a sample could have attained under the given experimental conditions.
For each experimental suite N eq and τ are estimated by minimizing the difference between predicted (equation 33) and observed bubble number densities (Supplementary Figure 3). The resultant values, shown in Supplementary Table 2, are similar for samples within a given experimental suite and predominately depend on the initial and final pressures. Non-dimensional BND is predicted to increase linearly between 0 ≤ t post ≤ τ, after which it reaches its maximum value [26]. Thus, non-dimensionalization should collapse all experiments during which nucleation was continuous onto the predicted line by equation (33). Figure 2 shows that this is the case within the expected experimental variability. Given that the maximum value of τ in our experiments was ≈1000s, our experiments substantiate that bubble nucleation in rhyolitic melt can be a continuous process for at least that long.  [33,34]. Similar results are obtained for different discharge rates and/or saturation pressures. The simulations shown in blue have no bubbles prior to eruption, whereas those in red have an initial volume fraction of 0.03 bubbles. The simulations shown in thick blue and red lines are base cases. They were obtained by adjusting conduit radius such that the decompression rate is greater than the inverse characteristic H 2 O diffusion time at all depths above where the first bubbles nucleate (for more details see equation (30) in Methods). No other parameters were adjusted to obtain the simulation results shown. The thin blue and red simulations illustrate model sensitivity to variations in predicted conduit radii. They were obtained for conduit radii that are shifted to smaller and larger values relative to the base cases.

Key results for bubble nucleation during magma ascent
The key result is an upward narrowing conduit as a necessary requirement for continuous nucleation.
Thus, in all cases the BSDs shown in Figure 4 are model predictions that were obtained by constraining conduit radius such that the decompression rate is slightly greater than the inverse characteristic H 2 O diffusion time. Variations in conduit radius  (Table S2). All experiments fall along the trend expected for continuous nucleation, demonstrating that that bubble nucleation in rhyolitic melts is continuous across a wide range in conditions. Given the range of τ n in our experiments, nucleation was continuous for at least 1000 seconds.
primarily affect decompression rate, fragmentation depth, and bubble number density, but still yield results for all relevant parameters that are well within range of observed values for Plinian eruptions. At very low nucleation rates (J < 10 9 m −3 s −1 ) nucleation is quite sensitive to small variations in conduit radius (decompression rate), resulting in oscillations in nucleation rate. These become, however, quickly damped out as bubble number densities, N m exceed values of 10 10 bubbles per m −3 of melt. For all cases shown, the match between predicted bubble size distributions and those measured in Plinian pyroclasts is remarkable and suggests that bubble nucleation during eruptions is indeed continuous throughout magma ascent, presumably commencing shortly above the magma chamber or perhaps already during magma withdrawal from the chamber.

Simulation details
As bubbles nucleate and grow, the distance between bubbles decreases. Consequently the characteristic diffusion rate, which scales as the inverse squared distance between bubbles, increases.
To maintain supersaturation and nucleation decompression rate has to increases with distance above the magma reservoir. An upward narrowing conduit increases the average magma velocity, which increases both the hydrostatic and dynamic components of the momentum balance. At the same time, due to bubble nucleation and growth, the average water concentration within the melt decreases as the magma rises to shallower depths. Even-  tually, this results in a sufficiently large increase in viscosity, such that the pressure loss due to wall friction results in a pos-itive feedback between decompression rate, water exsolution and viscosity. In Figures 3a and 3b this manifests itself as the pronounced curvature in P m and decompression rate prior to fragmentation. the rapidly increasing decompression rate results in increasing supersaturation and, hence, nucleation rate, bubble number density, and volume fraction of bubbles (Figures 3c, 3d and 3e, respectively). Once the feedback between decompression rate, water exsolution and viscosity becomes dominant the conduit radius no longer needs to decrease in order to maintain supersaturation ( Figure 3f). Instead, It is this runaway feedback that produces the high overall bubble number density as well as the small bubble diameter tail in the size distribution. The latter can also be found in Plinian pyroclasts, where it is associated with a decrease in the power-law exponent of the distribution with decreasing bubble diameter ( Figure 4).
Upon fragmentation gas pressure inside bubbles quickly equalizes with magma pressure, as a consequence of permeable gas flow from pyroclasts into the surrounding fractures and producing a rapid increase in gas volume fraction, as the gas-pyroclast mixture expands (Figure 3e). Outgassing of the pyroclasts is simulated using a pore-pressure relaxation rate (equation (23) in Methods), based on measured pyroclast permeabilities [35]. At the same time decompression rate decreases because of the abrupt change in viscosity, which is much lower for the gas-pyroclast mixture above fragmentation than for the bubbly magma ( Figure 3b). As a consequence there is rapid water diffusion out of the melt and loss of supersaturation conditions. H 2 O continues to exsolve en route to the surface, contributing together with the expansion of the exsolved H 2 O to the continued increase in gas volume fraction (Figure 3e). The mixture of gas-pyroclasts thus accelerates and reach choked flow conditions at the conduit exit (equation (27)), as commonly assumed in eruption models [3,33,36]. We assume that bubbles do not nucleate and grow within the pyroclasts that is after fragmentation. In other words, any remaining dissolved H 2 O exsolves through diffusion into existing bubbles and then enters the gas phase that surrounds pyroclasts by permeable flow within the pyroclasts [37].

Discussion
Continuous nucleation reconciles multiple observations from Plinian eruptions. When incorporated into numerical models for eruptive magma ascent, predicted bubble number number densities and vesicularities are within the range of Plinian pyroclasts ( Figure 3). Moreover, the predicted BSDs also mirror those measured in Plinian pyroclasts (Figure 4). They span several orders of magnitude in size and the observed change in power-law exponents from d = 3-4 to d = 1-2 with decreasing vesicle size [15,19,16,14] also exists in the simulated BSDs. This change in slope of the size distribution can be attributed to the aforementioned accelerating feedback between viscous pressure loss and nucleation.
Although magma decompression rate increases continuously as magma rises toward the level of fragmentation, the predicted time-averaged decompression rates ∼10 −2 MPa s −1 are well within the 0.001 − 1 MPa s −1 range of diffusion chronometers [38]. Thus, continuous nucleation fully closes the large previous gap in decompression rate estimates between diffusion chronometers and bubble nucleation models [24,38,27].  Figure 3 (red and blue curves), together with distributions from Plinian pyroclasts (gray curves). The latter are from the 1912 Novarupta eruption in Alaska [15], the 1875 Askja eruption in Iceland [16], the 2008 Chaiten eruption in Chile [19], and the 7.7 ka Mazama eruption, Oregon [14]. Plinian BSDs have been subdivided into two overlapping power law distributions with exponents of d ≈ 3-4 and d ≈ 1-2. In the simulations the latter size fraction of the distribution represents bubbles nucleated during the rapid increase in decompression rate in the depth range where the feedback between water exsolution and viscous pressure loss accelerates up to fragmentation.
The overarching result, however, is the necessary requirement that conduits widen with depth from a few 10s of m at shallow depths to radii of ∼100s to 1000 m at H 2 O saturation depths. The magma reservoir-conduit complex thus resembles in shape a wine decanter with a wide base. Alternatively, one could view our simulated 'conduits' as the top of magma chambers that form broad cupolas which narrow upward into the developing conduit. The implication would be that bubbles start nucleating during magma withdrawal within the upper reaches of the magma chamber, with the resultant decompression rates spanning a wide range in values [39]. Motivated by this consideration we also examined the role of pre-existing bubbles prior to eruption [40]. Although pre-existing bubbles somewhat reduces the required conduit widening, they do not affect the essential character of the simulation results, including the BSDs.
We have demonstrated experimentally that bubble nucleation in rhyolitic melt can be a continuous process over time intervals of at least thousands of seconds. Building upon this insight, numerical simulations of bubble nucleation and growth, in conjunction with the fluid dynamics of eruptive magma ascent, show that vesicle size distributions in Plinian pyroclasts can be explained by continuous heterogenous bubble nucleation during magma ascent. The ensuing average magma decompression rates are consistent with diffusion-chronometry geospeedometers. Moreover, the resultant bubble size distributions are similar to vesicle size distributions measured in Plinian pyroclasts, regardless of whether there already are bubbles in the magma prior to eruption or not. A necessary condition is that conduits widen with depth and transition via a broad cupola into the upper reaches of magma reservoirs.

Methods
We model the one-dimensional fluid dynamics of steady magma flow in the volcanic conduit, coupled with the nucleation and growth of H 2 O bubbles. Simulations start at H 2 O saturation, corresponding to a pressure, P H 2 O . Assuming a constant lithostatic pressure gradient of ρ rock g, where ρ rock = 2400 kg/m 3 is rock density and g is gravity acceleration, the saturation pressure P H 2 O corresponds to a depth, Z 0 = P H 2 O /(ρ rock g). Simulations end when magma reaches the surface at Z = 0. As magma flows up the volcanic conduit magma pressure, P m , decreases due to hydrostatic and viscous pressure loss. The magma thus becomes supersaturated in H 2 O, resulting in bubble nucleation. At the same time H 2 O diffuses from the melt into existing bubbles, as they grow. The fluid dynamics of magma ascent are coupled with the nucleation and growth of bubbles through magma pressure P m , which provides a depth-varying boundary condition for the bubble nucleation and growth calculations. In the subsequent sections we will detail the methodologies for employed to model magma flow in the conduit as well as bubble nucleation and growth within the ascending magma.

Magma flow in the conduit
We assume a vertical cylindrical conduit, whose radius can vary with depth, Z. We assume the flow is steady because the duration of magma ascent is much shorter than the duration of Plinian eruptions [41]. We assume flow is one dimensional and integrate flow properties over the cross-sectional area of the conduit. Below the level of fragmentation we define magma as the mixture of silicate melt and H 2 O bubbles. Above the level of fragmentation magma is the mixture of continuous H 2 O vapor with suspended fragments of vesicular magma, that is pyroclasts. Throughout the relative velocity between the two phases (melt and H 2 O vapor/fluid) is neglected. Below fragmentation this is justified because the buoyant rise velocity of bubbles is negligible, given the large melt viscosity [33]. Above fragmentation it is a commonly employed approximation [3,33] and one that does not affect the salient results of our simulations in any significant manner. The properties of the mixture, that is of the magma, are the volumetric average of the two phases. The flow is furthermore assumed to be isothermal, another common approximation that dose not significantly impact bubble nucleation rate [36].
With these assumptions, conservation of mass and momentum are [3,33] ∂(ρ m uA) and respectively. Here ρ is magma density, averaged over liquid and gas phases, (3) φ is the volume fraction of bubbles, u is magma ascent rate, A = πa 2 is the conduit cross sectional area, a is conduit radius, and ρ g and ρ l = 2400 kg/m 3 are gas and melt densities respectively. F fric is the frictional pressure loss estimated as ρ m u 2 f /a where f is the friction factor. Before fragmentation f = 16/Re + f 0 and after fragmentation f = f 0 . Re = 2ρ m ua/η is the Reynolds number, η is the viscosity of the mixture, and f 0 is a factor related to conduit wall's roughness and assumed to be 0.0025 [33].

By substituting equation (1) into (2) one obtains
Here and given that under constant entropy where c is the speed of sound within the magma. Equation (4) can thus be simplified to where M = u/c is the Mach number of the mixture. The sound speed before fragmentation is estimated from c 2 = K m /ρ m where K m is bulk modulus of the mixture The bulk modulus of the gas phase, K g , is calculated from equation of state [2] and the liquid is assumed to be incompressible [36]. The sound speed of the gas-pyroclast mixture after fragmentation is assumed to be equivalent to the sound speed in the gas phase [34]. Alternate models for the sound speed of the gas-pyroclast mixture [33] do not affect the outcome of the model results in terms of bubble nucleation, which occurs prior to fragmentation.

Bubble nucleation and growth
We assume H 2 O is the only volatile phase because it is most abundant and controls final bubble number density [32]. H 2 O exsolves through bubble nucleation and diffusion into already nucleated bubbles. We assume nucleation is heterogeneous and facilitated by abundant pre-existing magnetite nanolites [24,25]. We use the far field approximation to calculate H 2 O diffusion from melt into bubbles [10]. Lastly, we assume that permeable escape of H 2 O vapor from the bubbles is negligible before magma fragments, but account for it after fragmentation. The bubble nucleation and growth model is in the Langrangian frame of reference. It is integrated to the steady state equations for magma flow in the conduit, which is in Eulerian frame of reference, by considering d/dt = u∂/∂z.
The number of bubbles in volcanic systems is too high to track growth of each individual bubble. We therefore use the method of moments to calculate the evolution of the bubble population [10]. The corresponding moments of the population, µ, are defined as where R is the bubbe radius, Λ(R, t) is the bubble population per unit volume of melt within the interval of R and R+dR, and subscript k = 0-3 refers to the order of the moment. Each moment refers to a measurable characteristic quantity [10]. µ 0 is the bubble number density, that is the number of bubbles per unit volume of melt, µ 1 is the sum of bubbles radius, 4πµ 2 is the total surface area of bubbles, and 4/3πµ 3 is the total volume of bubbles. The evolution of moments through time is given by and for , k ≥ 1, by where J is nucleation rate of bubbles, G(R) is the growth rate of bubbles assumed to be equal for all bubbles and equivalent to the growth rate of a bubble with mean bubble size,R = µ 1 /µ 0 , and R c is the critical size for nucleating bubbles.

Below fragmentation
We use classical nucleation theory to estimate nucleation rate of stable bubble nuclei as a function of supersaturation pressure. Nucleation rate of stable bubbles are [11] and bubbles are stable if they are larger than a critical size, R c , given by Here T is the absolute temperature, k B is the Boltzmann constant, γ is the surface tension of bubble nuclei, P n is the pressure inside a bubble nucleus, and P m is pressure in the mixture. p n is related to the saturation pressure of volatiles, p sat , through [42] f (p n , T )p n = f (p sat , T )p sat e Ω(pm−psat)/k B T , where f (p, T ) is the fugacity coefficient, and Ω is the volume of volatile molecules. W is the change in free energy as a result of nucleation of a bubble and is given by where α is the heterogeneous nucleation factor and depends on the contact angle, θ, between bubble nuclei and pre-existing crystals. It is defined as and we assume nucleation is facilitated by magnetite crystals with θ = 145 • [43,44].
The pre-exponential factor, J 0 , in equation (12) is defined as where n 0 is the concentration of volatiles molecules in the melt, D is the diffusion coefficient, a 0 is the average distance between volatiles molecules.
After nucleation bubble nuclei grow by H 2 O diffusion, because the concentration of H 2 O at the bubble-melt interface, C R , is lower that the concentration in the surrounding melt, C m . The H 2 O flux into bubbles is given by where D is diffusion coefficient, r is the distance from bubble's center, R is bubble radius, and c is the water concentration in the surrounding melt given by ∂c ∂t + dR dt At low supersaturation pressure the left hand side in equation (19) can be neglected [45]and the concentration gradient at the melt-vapor interface can be approximated as The mass of H 2 O inside bubbles, m g , will increase due to diffusion at a rate dm g dt = (4πµ 2 )ρ m q, whereas the bubble growth rate is given by [46] G(R) =R 4η where the inertia terms are neglected [10]. Here η is the viscosity of melt, p g is the pressure the H 2 O fluid (vapor) inside bubbles, estimated using the equation of state for H 2 O.

Above fragmentation
After fragmentation we assume bubbles that bubbles no longer nucleate or grow [37]. In other words, J = 0 and G(R) = 0. H 2 O, however, continues to exsolve into bubbles with a rate given by equation (18). Exsolved H 2 O escapes by permeable flow from within the pyroclasts into the gas phase surrounding the pyroclasts. Consequently, the gas pressure within pyroclasts decreases at a rate of where τ k is the characteristic time scale for permeable outgassing estimated from Darcy's law as Here l ≈ 10cm is the characteristics length scale [37], k ≈ 10 −12 m 2 is permeability [35], η g = 10 −5 Pa.s is viscosity of the gas phase, and β is the compressibility of the gas phase, and φ is porosity of bubbles in the pyroclasts.
Lastly, the concentration of dissolved H 2 O and thus the saturation pressure, both below and above fragmentation, decrease as a result of the diffusion of water into bubbles. The resultant conservation of H 2 O requires that where ρ l is the melt density, assumed to be constant throughout magma decompression, and m c is the mass of a bubble nuclei estimated from equation of state.

Model simulation
The parameters in the governing system of equations are either specified or calculated from existing formulations: H 2 O solubility [47], diffusion coefficient [48], equation of state [2], fugacity coefficient [2], surface tension [27], melt viscosity [12], and the molecular volume of H 2 O [1]. For a given simulation we integrated equations (7), (10), (11) for k = 1 through 3, as well as equations (23) and (25) using the ode15s function of MATLAB ® . The boundary condition at the initial depth, z 0 , are p m (z 0 ) = ρ r gz 0 , C m (z 0 ) = C eq (p m (z 0 )), Here ρ r = 2400 kg/m 3 is the rock density, C eq are the equilibrium water concentration, and N 0 and R 0 are the number density and radius of pre-existing bubbles. The boundary conditions at the surface are We assume magma fragments at a critical porosity, φ cr . For each simulation run, we vary φ cr to meet the boundary condition at the surface. We estimate bubble size distribution for each simulation by post-processing the simulation results. We discretize z into multiple bins and estimate number of bubbles nucleated at each bin as and the final size of bubbles nucleated at each bin as [36] The objective of our model is to find conditions at which bubble nucleation is continuous from nucleation onset until magma fragmentation. Nucleation is driven by supersaturation pressure, p sat − p m and is thus controlled by the competition between decompression rate, d p m /dt and d p sat /dt. The latter is proportional to the diffusion rate and is estimated from equation (25). To maintain sufficient supersaturation for nucleation, we assume that the decompression rate at any given depth is greater than diffusion rate, that is where λ > 1 is a constant. This approach requires that d p sat /dt > 0, that is sufficient number of bubbles are nucleated, to be able to estimate decompression rate. We use equation (30) at minimum bubble number density of 1 mm −3 . From decompression rate we estimate conduit cross sectional area, dA/dz from equation (7) analogous to the approach in Mastin and Ghiorso [33]. After finding an initial conduit radius we used the MATLAB ® pchip function to smooth the obtained function and assess the sensitivity of model predictions to variations in conduit radius.

Nucleation time scale in experiments
Bubble nucleation is driven by supersaturation pressure, ∆P ss = P sat −P, defined as the difference between pressure at which H 2 O would be saturated and the sample's pressure, P. ∆P ss increases as P decreases.
Considering that diffusion is ineffective during decompression [27], P sat is expected to remain the same as the initial pressure, P i , throughout decompression. ∆P ss thus reaches the maximum potential supersaturation pressure when sample pressure reaches the final pressure, ∆P ss = ∆P max , Supplementary Figure 1. Consequently, the nucleation rate in the experiments is expected to increases as sample pressure decreases toward P f , because ∆P ss ≈ ∆P and nucleation rate scales as exp(−1/∆P 2 ss ) [11]. Nucleation rate reaches a maximum value at P f and nucleation continues at this rate until the sample is quenched, unless at some point diffusion becomes non-negligible. Thus, the majority of bubbles are expected to nucleate once ∆P has reached its maximum ∆P max = P i − P f . If (and once) enough bubbles have nucleated to decrease the characteristic diffusion time to the point where water diffusion into existing bubbles becomes non-negligible, then ∆P ss < ∆P max and nucleation rate decreases.
We used the nucleation model described in Yamada et al. [26] to quantify nucleation rate during samples hold time from observed bubble number densities. The model was obtained by analytical solution of bubble nucleation and growth formulations described in Toramaru [10]. Nucleation rate, J(t), as a function of time is given by Here t is time at the final pressure, J s is the steady nucleation rate and τ is the nucleation time scale, which is the e-folding time of nucleation. τ n scales inversely with J s as [26] τ = k τ × J −2/5 s .
The coefficient k τ depends on melt properties, in particular the diffusion coefficient. Here we assume it is constant because the variability of k τ across the conditions of our experiments is relatively weak [26].
Integration of J(t) through time yields the bubble number density that are nucleated during annealing time as a function of time,N post (t), given by N post = N m + N 0 is related to observed bubble number densities where N 0 is the bubble number density at time t = 0. In our experiments N 0 accounts for bubbles that were nucleated during decompression and is obtained for each suite from N m in samples with t post = 0. Furthermore, Γ(x) is the gamma function and Γ(s, x) is the normalized lower incomplete gamma function, Both Γ(x) and Γ(s, x) arise out of the integration of J(t) (Equation 31).
At the nucleation time scale, t τ n , nucleation rate decreases and bubble number density eventually reaches the equilibrium value of Lastly, the predicted non-dimensional bubble number density is estimated from equations (33) and (36) N post (t) N eq = Γ 2 5 , (t/τ) 5/2 .