An experimentally validated numerical model for bubble growth in magma

and non-linear dependence on the concentration of H 2 dissolved in a numerical bubblegrowth modelshavepreviously applied, continuous,in situ experimentaldata or provided as a user-friendly present a numerical bubble growth model, implemented in MATLAB, allows for arbitrary temperature and pressure pathways, and accounts for the impact of spatial variations in dissolved H 2 concen- trationonviscosityanddiffusivity.Wevalidatethemodelagainsttwosetsofexperimentaldata:(1)Newcontin-uous data for gas-volume fraction as a function of time, collected using optical dilatometry of vesiculating hydrous obsidian samples which were heated from 930 °C to 1000 °C at atmospheric pressure. This dataset captures isobaric, isothermal bubble growth under strongly disequilibrium conditions. (2) Discrete data from pub- lished decompression experiments at 825 °C and pressures from 200 MPa to ~5 MPa with decompression rates from0.1MPas − 1 to10 MPas − 1 .Theseexperimentsrepresentisothermal,decompression-drivenbubblegrowth spanning equilibrium to strongly disequilibrium conditions. The numerical model closely reproduces the exper- imental data across all conditions, providing validation against contrasting bubble growth scenarios. The validated model has a wide range of potential volcanological applications, including forward modelling of bubble growth phenomena, and inverse modelling to reconstruct pressure – temperature – time pathways from textures and volatile contents of eruptive products. 2020 The Authors. B.V. an CC


Introduction
Bubble growth drives magma ascent in the shallow crust. The rate of bubble growth is an important control on whether an eruption is effusive or explosive, because rapid growth can lead to magma fragmentation (Degruyter et al., 2012;McBirney and Murase, 1970;Sparks, 1978;Verhoogen, 1951). Accurate modelling of bubble growth processes is therefore necessary for the investigation of many physical volcanological problems. However, explicit forward-modelling of bubble growth remains challenging because it involves several non-linear and coupled physico-chemical processes (Blower et al., 2001;Liu and Zhang, 2000;Lyakhovsky et al., 1996;Navon et al., 1998;Proussevitch and Sahagian, 1998;Sparks, 1978); consequently, a numerical approach is required to model bubble growth accurately. Furthermore, the pressure and temperature changes that drive and modify bubble growth may vary in a complex way as a packet of magma moves through the volcanic system from crustal reservoir to final emplacement (Carey et al., 2013;McIntosh et al., 2014). For practical applications, therefore, a model should be flexible, accurate, and numerically stable over nonlinear temperature and pressure pathways.
Bubbles grow through equation-of-state expansion when the gas pressure exceeds the ambient pressure confining the magma. The excess pressure usually arises from a combination of decreasing confining pressure and exsolution of volatile species (principally H 2 O) into the bubble. Growth is retarded by viscous stresses arising in the shell of silicate melt that surrounds the bubble, as it flows to accommodate the growth. Importantly, both melt viscosity and the diffusivity of H 2 O

Numerical model
We adopt the mathematical framework of Blower et al. (2001), which is, in turn, built on the work of Proussevitch and Sahagian (1996). The model solves the hydrodynamics of the melt surrounding the bubble, the diffusion of H 2 O through the melt, the equation of state of H 2 O in the bubble, and the mass balance of H 2 O between the melt and the bubble. We extend the framework to include the effects of surface tension and update the component models for material properties. The numerical implementation is described in detail in Appendix A, and a downloadable, user-friendly MATLAB implementation of the model with user manual is available as an online electronic supplement.

Geometry
We adopt a 'shell model' (Blower et al., 2001;Proussevitch et al., 1993) in which the bubbly magma is conceptually divided into units composed of a spherical bubble and its associated spherical shell of melt ( Fig. 1). This approach reduces the problem to 1D spherical symmetry, which simplifies formulation and computation, but it does not capture asymmetry introduced by bubble-bubble interactions, which become increasingly important as the bubble volume fraction increases. The bubble-shell unit is used to compute the growth and resorption behaviour of a single bubble in response to pressure and temperature changes. We follow the approach of Proussevitch et al. (1993) and generalize this to the behaviour of a body of bubbly magma by assuming that all of the units are identical, hence that the evolution of physical properties of the unitporosity, bubble radius, melt-shell viscosity, etc.
is identical to that of the bulk magma. This leads to the following relationships among bubble radius R, shell radius S (where both R and S are measured from the bubble centre; Fig. 1), the gas volume fraction ϕ, and the bubble number density per unit volume of dense magma N b : where V m is the volume of the melt and V g is the volume of the gas in the shell-bubble system. We define the radial thickness of the melt shell as L = S − R (Fig. 1), such that Together Eqs. (1) & (2) capture all the geometrical relationships required in our model.

Governing equations
A modified form of the Rayleigh-Plesset equation describes the hydrodynamics of the growth of a spherical bubble in a finite incompressible shell of liquid in the absence of inertial effects (Proussevitch et al., 1993) Proussevitch et al. (1993) and Blower et al. (2001). a) A single unit cell, comprising a spherical bubble in a concentric spherical melt-shell. The line x 1 → x 2 is a radial transect across the melt shell; other symbolic notation is defined in the main text in Section 2. b-d) Schematic representation of spatial variation in concentration of dissolved H 2 O, melt viscosity, and diffusivity of H 2 O, across the radial transect x 1 → x 2 during bubble growth.
where p g is the pressure of the gas in the bubble, p ∞ is the hydrostatic pressure acting on the outside of the liquid shell, p c = 2Γ/R is the capillary pressure at the bubble wall, for which Γ is the gas-liquid interfacial tension, η is the viscosity of the liquid in the shell, and t is time. For magmatic systems, the viscosity of the melt is a strong function of local dissolved H 2 O content; since this varies radially over time as H 2 O diffuses in and/or out of the bubble, the radial distribution of viscosity also evolves over time. Blower et al. (2001) adapted Eq. (3) to account for changes in the radial distribution of viscosity by integrating over the shell thickness in spherical coordinates. We adopt this approach, and additionally include the capillary pressure, which was neglected by Blower et al. (2001), to improve fidelity for bubbles that are sufficiently small that the capillary term is important: where R 0 and S 0 are the radius of the bubble and melt shell at t = 0, and x is the radial Lagrangian coordinate which represents the spatial coordinate in the melt shell. The use of a radial Lagrangian coordinate system simplifies the mathematical formulation of the integral because, in this system, the radial position x of any infinitesimal shell of melt does not change as the bubble grows or shrinks. The Lagrangian radial coordinate x is related to the Eulerian radial coordinate A by conservation of volume Evaluation of Eq. (4) requires the radial distribution of viscosity η(x) at time t, which depends on the spatial distribution of H 2 O concentration in the melt C(x) at time t. This is computed by solving the one-dimensional form of Fick's second law in spherical coordinates in the Lagrangian coordinate system (Blower et al., 2001;Braithwaite et al., 1999): where D is the diffusivity of H 2 O in the melt, which appears inside the outer spatial derivative because it depends on the local H 2 O concentration. The Lagrangian reference frame accounts for the advective flux of H 2 O molecules towards or away from the bubble wall as the bubble grows or shrinks (Blower et al., 2001). The radial H 2 O concentration distribution C(x) is converted to a radial viscosity distribution η(x) via a constitutive law for η(C), which is discussed later in Section 2.3.3, allowing the integral in Eq. (4) to be evaluated. The radial H 2 O concentration distribution is used to compute the mass m of gas in the bubble at time t by assuming conservation of mass of H 2 O within the bubbleshell unit where m 0 is the initial mass of gas in the bubble, ρ m is the density of the melt and C is in wt%. The mass m(t) and the bubble radius R(t) can then be converted to p g at any given pressure and temperature using an equation of state for the gas (or supercritical fluid) phase, which is discussed later in Section 2.3.2. For most situations of interests, the system of equations (Eqs. (1)- (7)) must be solved numerically to compute R(t) and ϕ(t).

Material properties
Solution of the governing equations requires models for various physical properties of the materials, as functions of pressure p g or p ∞ , temperature T, and concentration of H 2 O dissolved in the melt C. The relevant physical properties are: 1) solubility of H 2 O in the melt phase C s (p g , T); 2) diffusivity of H 2 O through the melt D(C, p ∞ , T); 3) equation of state for the gas (or supercritical fluid) phase ρ g (p g ,T); and 4) the viscosity of the melt phase η(C, T). In this section we present constituent models for these physical properties. The models are appropriate for the volatile species (H 2 O) and melt composition (rhyolite) used in the experiments that produced the data we use to validate our model (Section 3). However, the flexibility of the general mathematical framework presented above allows different constituent models to be chosen for other volatile phases and melt compositions. Effects associated with the heat of vaporization are neglected in our model. This is explored in Sahagian and Proussevitch (1996), but may be a minor effect at most magmatic conditions (Liu et al., 2005).

Solubility and diffusivity
The solubility sets the concentration boundary condition for the gasmelt interface, determining whether H 2 O will exsolve from the melt into the bubble, or resorb from the bubble back into the melt. Liu et al. (2005) present an equation for the solubility of H 2 O C s ¼ 354:94p g 0:5 þ 9:623p g −1:5223p g where p g is in MPa and T is in K. Eq. (8) is based on experiments on synthetic haplogranites and natural rhyolites at 700 to 1200°C and 0 to 500 MPa. This model is poorly constrained at low pressure and the dataset against which we validate our numerical model includes results from experiments conducted at atmospheric pressure (0.1 MPa). Therefore, we derive a low-pressure solubility model by fitting the same functional form as in Eq. (8) to C s (T) data from Fig. 5b and Table 3 of Ryan et al. (2015b), for rhyolitic melt at 0.1 MPa and 900-1100°C. Here we neglect the pressure dependence and aim to find a simple model for the temperature dependence of the solubility of water specifically for use with 0.1 MPa experiments. We find: The starting material used by Ryan et al. (2015b), and used to constrain Eq. (9), is identical to the material that we use in our 0.1 MPa experiments. We use either Eq. (8) or (9) to calculate C s , choosing the equation that is more appropriate for the relevant experimental conditions.
At low H 2 O concentration (C ≲ 2 wt%), the diffusivity of total H 2 O in rhyolite is captured by Eq. (15) of Zhang and Ni (2010): where p ∞ is in GPa and T is in Kelvin. This equation is appropriate for the continuous, in situ experiments, which are conducted at atmospheric pressure. At higher H 2 O concentrations (C N 2 wt%), diffusivity shows an exponential dependence on H 2 O concentration, hence a more sophisticated formulation for D(C, p ∞ , T) is required (e.g. Coumans et al., 2020). We use a combination of Eqs.

Equation of state
Solution of Eqs.
(1)- (7) gives the time evolution of the radius of the bubble R and the mass of H 2 O in the bubble m, from which the density of H 2 O in the bubble is calculated as where ρ g is in mol cm −3 , m is in grams, M is in g mol −1 , and R is in cm. These units are chosen to be consistent with constants used in the equation of state of Pitzer and Sterner (1994), which is used to determine the pressure in the gas (or supercritical fluid) phase where p g is in MPa, T is the temperature in Kelvin, G = 83.144 cm 3 bar K −1 mol −1 is the gas constant, and a 1 − a 10 are temperature-dependent parameters given in Pitzer and Sterner (1994).

Melt viscosity
The melt viscosity η is computed as a function of temperature and local H 2 O concentration using the empirical parameterization of Giordano et al. (2008). This parameterization has been calibrated for a range of melt compositions common in terrestrial volcanic rocksmafic to silicic, subalkaline to peralkaline, and metaluminous to peraluminousagainst experiments spanning a temperature range of 525-1705°C and H 2 O concentrations from 0 to 8 wt%. The model is based on the Vogel-Fulcher-Tamman equation for non-Arrhenian temperature dependent viscosity: where η is in Pa s, and J 1 , J 2 , and J 3 are parameters that depend on melt composition, computed following Giordano et al. (2008).

Numerical methods
A numerical solution has been developed, in the MATLAB programming environment, to solve the governing equations presented in Section 2.2 and compute the required material properties via the auxiliary equations presented in Section 2.3. The principal output is the bubble radius as a function of time R(t), which is computed by solving Eq. (4) for dR/dt and integrating. The solution of Eq. (4) requires that the spatial distribution of dissolved H 2 O be known at each timestep; this is determined through solution of Fick's second law (Eq. (6)) via the method of lines. This approach involves approximating the spatial derivatives using finite differences such that the partial differential equation is transformed into a coupled system of ordinary differential equations in time, which we solve using MATLAB's in-built solver ODE15s (Shampine and Reichelt, 1997). Solutions are computed on a radial line through the melt shell, discretized into nodes that are spaced logarithmically to give the highest spatial resolution near the bubble wall, where the concentration gradient is highest.

Initial conditions and boundary conditions
The numerical model requires initial values for the bubble radius R and shell thickness S, the pressure in the bubble p g and at the shell's outer wall p ∞ , the isothermal system temperature T, the mass of H 2 O in the bubble m 0 , and the spatial distribution of H 2 O concentration in the shell C(x). In principle, these quantities can be set to any value that is within the range of validity for the models used to compute material properties. In practice, we initialize the model to equilibrium conditions, such that: (1) the temperature, pressure, bubble radius, and mass of H 2 O in the bubble satisfy the equation of state (Eqs. (11) and (12)); (2) the pressure in the bubble balances the sum of the ambient pressure and capillary pressure, i.e. p g = p ∞ + p c ; (3) the H 2 O concentration in the shell is uniform (i.e. the same for all x). Once the model is initialized, p ∞ and T may be varied arbitrarily as functions of time to drive changes in the modelled system. Initial bubble radius can be set by either defining an initial gas volume fraction ϕ and computing the radius from bubble number density N b , or by setting the radius directly. The model does not include the bubble formation process and cannot be initialized with a bubble radius (or gas volume fraction) of zero; hence, when modelling growth immediately following bubble formation, bubble radius must be initialized to a small but non-zero value. In practice, when we wish to model bubble growth from formation, we set the initial gas volume fraction to the arbitrary value of ϕ 0 = 10 −6 , which is chosen to be small enough that subsequent growth is independent of the choice of ϕ 0 .
Boundary conditions at the inner and outer edges of the melt shell are required. At the inner edge (i.e. the bubble-melt interface) the concentration of dissolved H 2 O is set to the solubility value, determined for current values of p g and T, via Eq. (8) or (9). At the outer edge, a zeroflux boundary condition is imposed by setting the spatial gradient of the concentration of dissolved H 2 O to zero. A further condition is that the total mass of H 2 O in the bubble-shell unit (i.e. both dissolved and as a free phase) is conserved throughout the simulation.

Availability of code
Details of the implementation of the numerical model, including pseudo-code, are provided in Appendix A. A user-friendly version of the code, implemented in MATLAB, is available to download in the Supplementary information. A brief user guide is available as an electronic supplement.

Experimental methods
We validate the model presented in Section 2 against two contrasting and complementary suites of experimental data. The first dataset, generated specifically for this study, captures the continuous evolution of the gas volume fraction ϕ of a vesiculating rhyolite melt during heating at 1 atm. This approach, which is chosen principally because the resulting data are well suited to model validation, also has direct application to the natural scenario of thermal vesiculation of magma (Lavallée et al., 2015). Continuous data provide a strong test of the bubble growth model because: (1) they contain a much larger number of individual ϕ(t) datapoints than conventional discontinuous datasets from high-pressure, high-temperature (HPHT) experiments; (2) a single sample can be tracked through all stages of bubble growth; and (3) the dataset is not complicated by poorly-constrained resorption processes during quench (McIntosh et al., 2014).
The second dataset is drawn from published HPHT decompression experiments (Burgisser and Gardner, 2004;Hamada et al., 2010;Mourtada-Bonnefoi and Laporte, 2004). While these experiments are discontinuous in nature, they complement the continuous experiments because: 1) they include significant changes in pressure, replicating conditions that are typical of natural volcanic systems; and 2) they extend the range of methodologies and compositions considered. Fig. 2 shows, schematically, the different pressure and temperature pathways represented by the two suites of experimental data.

Continuous isobaric experiments
These experiments use a tholeiitic rhyolite obsidian from Hrafntinnuhryggur, Krafla (Iceland) as starting material (see glass composition for sampling site 'AO' given in Tuffen and Castro, 2009). Ryan et al. (2015b) used glass from the same outcrop location to collect the data that we use to calibrate the solubility model for atmospheric pressure (Eq. (9)). They measured the initial H 2 O concentration of the glass to be C i = 0.114 ± 0.013 wt% using Fourier transform infrared spectroscopy (FTIR).
Continuous ϕ(t) data were collected during heating of the obsidian using a Hesse Instruments EM-201 Optical Dilatometer, following Wadsworth et al. (2016). Obsidian cores of 3 to 5 mm diameter and height are loaded into a mini-tube furnace and back-lit with a halogen lamp such that the silhouette is projected towards a CCD camera capturing images continuously at 1 Hz. Images are stored for analysis when the temperature recorded by the thermocouple under the sample plate changes by N1°C, or when the area of the image occupied by the sample has changed by more than 1% since the last stored image. The samples were heated at a constant rate to a dwell temperature in the range 930 − 1000°C. Experimental conditions are given in Table 1. The evolving gas volume fraction of the sample during the experiment is reconstructed from the imagery by applying a Canny edge detection algorithm to find the outer edge of the silhouette, then applying the solid of revolution technique (Wadsworth et al., 2016) to find the volume V(t). We correct for melt/glass expansion with changes in temperature by correcting the data to background expansion observed prior to bubble nucleation. The gas volume fraction is calculated as where V 0 is the initial, bubble-free volume of the sample. All sample dimensions were measured after experiments to confirm the volume calculations from image analysis; the overall uncertainty was found to be of the same order as the uncertainty derived from the pixel resolution, which is~0.6% relative.

Decompression experiments
We compile published data from decompression experiments performed by Burgisser and Gardner (2004), Mourtada-Bonnefoi and Laporte (2004), and Hamada et al. (2010). A summary of relevant experimental details is given below, and full details can be found in the original works. Burgisser and Gardner (2004) conducted experiments on natural rhyolite glass, with 75.6 wt% SiO 2 , from the Panum Crater dome (CA, USA). Experimental charges were held with excess H 2 O at 825°C and 150 MPa in an externally heated pressure vessel for 5 days to ensure uniform saturation. Samples first underwent a nucleation step, in which the pressure was dropped from 150 MPa to 100 MPa 'instantaneously', and held for 10 to 15 min to produce a unimodal population of 'seed' bubbles at equilibrium and with a bubble number density that is repeatable within a narrow interval ( Table 2). Samples were then decompressed at different controlled rates, to different final pressures p f using a stepwise method to approximate linear decompression. Upon reaching final pressure, the samples were quenched rapidly, cooling below the glass transition in around 3 to 10 s. Relevant experimental parameters are summarized in Table 2.
Mourtada-Bonnefoi and Laporte (2004) and Hamada et al. (2010) conducted experiments on natural rhyolite glass, with 77.4 wt% SiO 2 , from Güney Dagi, Turkey. Both studies followed similar experimental procedures, which differed from Burgisser and Gardner (2004). Hydrated glass samples were made with 6.83 to 7.05 wt% H 2 O (Mourtada-Bonnefoi and Laporte, 2004) and 6.6 wt% H 2 O (Hamada et al., 2010), which acted as the starting material for decompression experiments. For both studies, the samples were loaded into an externally heated pressure vessel and first pressurized to~250 MPa and then heated to either 700°C (Hamada et al., 2010) or 800°C (Hamada et al., 2010;Mourtada-Bonnefoi and Laporte, 2004) iñ 60 to 90 min. These temperatures and pressure conditions, at which the samples were H 2 O undersaturated, were chosen to ensure that the samples were bubble-free before decompression, in contrast with the experiments of Burgisser and Gardner (2004). Experimental charges were decompressed at different rates to a variable final pressure, then quenched isobarically. In this study, we use data only from experiments in the Mourtada-Bonnefoi and Laporte (2004) and Hamada et al. (2010) suites that were decompressed at 1 MPa s −1 in order to provide continuity between the datasets collected at 700°C and 800°C. Relevant experimental parameters are summarized in Table 2.  (Burgisser and Gardner, 2004;Hamada et al., 2010;Mourtada-Bonnefoi and Laporte, 2004). For each plot, i, ii, and iii represent times for which schematic representations of the bubble-shell unit are shown underneath. In a): bubble growth (i, ii, iii) is driven by the increase in temperature, which reduces H 2 O solubility. In b): bubble growth (i, ii) is driven by decompression, while bubble shrinkage (iii) is driven by decreasing temperature, which causes volume decrease by equation of state and increased solubility, which drives resorption of H 2 O into the melt shell. All these processes are captured by the numerical model (Section 2).

Modelling experimental data
The bubble growth model does not capture the bubble nucleation process; consequently, certain assumptions are required to initialize the model for application to the experimental data. For experiments that begin undersaturated, we assume an initial gas volume fraction of ϕ 0 = 10 −6 (a value sufficiently small that the subsequent growth is not dependent on this choice of ϕ 0 ) at H 2 O saturation, where the bubble radius is computed from the chosen initial N b using Eq. (1). For the range of N b used in these experiments initial radii are from 0.5 to 5 μm. It is important that the initial bubble size is larger than the critical bubble , so that surface tension does not cause the bubble to shrink rapidly to a singularity. For some experiments there are material quantities (for instance bubble number density or initial dissolved H 2 O concentration) that are not known sufficiently accurately to allow the initial conditions to be defined a priori. These quantities are determined by minimization to the experimental data. The specific steps required to apply the numerical model to the various experiments are presented below.

Continuous, isobaric experiments
We apply our model to the in-situ bubble growth experiments using a parameter sweeping approach in which we vary the initial H 2 O concentration C i and bubble number density N b ; the ranges were chosen based on the experiments of Ryan et al. (2015b). We set the initial temperature as the temperature at which the melt becomes saturated for the assumed C i , as determined by the solubility law (Eq. (9)). The thermal history T(t) follows the experimental pathway from the point of saturation. Pressure p ∞ = 0.1 MPa throughout. Model outputs are compared with the experimental data to investigate goodness-of-fit and each model run is calibrated to the relevant ϕ(t) dataset using the root mean square error (RMSE) between data and model. In addition to quantifying the overall goodness-of-fit, the calibration provides best-fit values and uncertainties for poorly constrained input parameters and quantifies the model sensitivity to input parameters.

Decompression experiments
The Burgisser and Gardner (2004) experiments are modelled from the beginning of the decompression step and samples are assumed to start in equilibrium with the conditions at the end of the nucleation step. The initial gas volume fraction ϕ i and H 2 O concentration C i are calculated using the solubility law (Eq. (8)) and the conditions of the nucleation step (pressure drop from 150 to 100 MPa at 825°C) under the assumption that all H 2 O exsolved from the melt goes into the bubbles. Model runs are performed by parameter sweeping over bubble number density, which is varied over the experimentally observed range of 0.4 × 10 12 to 10 × 10 12 m −3 ( Table 2), and over pressure pathways that bracket the experimental values: decompression rates of 10, 1, and 0.1 MPa s −1 to final pressures p f in the range 95 to 10 MPa. When p f is reached the models follow an isobaric linear quench pathway (Fig. 2b) of 100 K/s, consistent with the experiments.
The Mourtada-Bonnefoi and Laporte (2004) and Hamada et al. (2010) experiments are modelled together, since they adopt the same experimental procedure. We use the initial H 2 O concentration C i reported in the original studies, which ranges from 6.6 to 7.05 wt%. The initial gas volume fraction ϕ i was 10 −6 and bubble number density N b had a range of roughly 0.5 × 10 12 to 10 × 10 12 m −3 . The initial pressure is set at a value consistent with the pressure at which bubble nucleation occurs as extracted from Hamada et al. (2010) Fig. 4 (p i = 52 MPa at T = 700°C, and p i = 85 MPa at T = 800°C). Model runs are performed at a decompression rate of 1 MPa s −1 . When p f is reached, the models follow an isobaric linear quench pathway of q quench = 100 K/s (Fig. 2b) consistent with the experiments. Both suites of decompression experiments include samples with H 2 O contents that are higher than the stated range of validity of the diffusivity model presented in Eq. (10), so we use instead the diffusivity model encapsulated in Eqs. (1), (7a), (13), and (14) of Zhang and Ni (2010). This model is more complex than Eq. (10) but is valid over the appropriate range of H 2 O concentrations.
Each model run yields a complete ϕ(t) pathway for the sample from the start of the experiment to the end of the quench. For each set of modelling runs that bracket the conditions for an experimental suite, the minimum and maximum ϕ at each p f is computed from among the model runs. These limits define a model range for comparison with the experimental data.

Modelling of continuous, isobaric experiments
The samples heated at atmospheric pressure show a characteristic sigmoidal ϕ(t) growth curve (Fig. 3). Sample growthwhich corresponds to bubble growthis slow initially, then accelerates as temperature increases and melt viscosity drops, before slowing again as the sample approaches equilibrium saturation at atmospheric pressure and the experimental dwell temperature. Results are qualitatively similar to those of Ryan et al. (2015b) except that their experiments do not show the slow and accelerating early growth phases, presumably because their samples are loaded directly into a furnace pre-heated to the dwell temperature, such that the samples heat rapidly and most of the growth occurs at the dwell temperature. The slow and accelerating growth phases indicate that the high viscosity of the melt at lower temperatures plays an important role in limiting bubble growth rate in our experiments.
The best-fit model ϕ(t) data show close agreement with the experimental data (Fig. 3a, sample K_2_975), faithfully capturing all stages of accelerating and decelerating growth. The companion plot (Fig. 3b) shows the dependence of goodness-of-fit on the values of the two swept parameters: bubble number density N b , and initial H 2 O concentration C i . For this dataset there is a clear minimum of root-meansquare-error (RMSE) at N b = 2.11 × 10 10 m −3 and C i = 0.1106 wt%. Labelled model ϕ(t) curves show how the model deviates from the data as the bubble number density and initial H 2 O concentration are varied, illustrating the sensitivity of the goodness-of-fit to those parameters (Fig. 3c). This analysis shows that the initial H 2 O concentration controls the equilibrium bubble volume fraction, and the bubble number density controls the rate of approach to equilibrium. Higher initial H 2 O concentration leads to higher equilibrium bubble volume fraction at constant pressure simply because there are more exsolved moles of H 2 O at equilibrium. Higher bubble number density leads to faster bubble growth because the distance over which H 2 O must diffuse through the melt to exsolve into the bubble is shorter, and because the shell of viscous melt that provides resistance to growth is thinner. The model can reproduce growth curves closely for all experiments in which the final bubble volume fraction ϕ ≲ 0.4. However, the quality of constraint on the best-fit parameters depends on the nature of the experimental data, as shown in Fig. 4. For example, the data from sample K_3_975 (Fig. 4a) are very similar to those for sample K_2_975 (Fig. 3a), except that the equilibrium bubble volume fraction plateau is not as well defined. This results in a distinctive covariation of the swept parameters, manifesting as a trench-like RMSE minimum in the goodness-of-fit plot (Fig. 4b). This is more pronounced for sample K_4_1000 (Fig. 4c, d) which does not reach the equilibrium plateau during the experimental run time, and still more pronounced for sample K_6_1050 (Fig. 4e, f) which does not clearly reach the inflection between accelerating and decelerating growth. The result is that equilibrium gas volume fraction is poorly constrained, leading to poor constraint on initial H 2 O concentration. This, in turn, facilitates covariation of bubble number density: for instance, a lower initial H 2 O concentration can be compensated by a higher bubble number density to give a similar ϕ(t) curve over the data range. Where there is strong co-variation, if an independent constraint on one of the variables (C i or N b ) is available, then minimization of the model can be used to determine the value of the other variable. Fig. 4g, h shows a sample that approaches ϕ~0.25. The data come close to the equilibrium plateau, so the minimum is well defined, but because the final gas volume fraction is low the constraint on the initial H 2 O concentration is rather poor. This manifests as a wider RMSE trench than is found for samples that approach a higher equilibrium gas volume fraction.
For samples that reach a final bubble volume fraction ϕ ≳ 0.4, the model is not able to capture accurately the whole of the sample growth curve. Fig. 5 presents data for sample K_1_1050, which approaches an equilibrium gas volume fraction ϕ = 0.64. This equilibrium value is used to compute the initial H 2 O concentration C i = 0.114 wt% (via the equation of state and Eqs. (1) and (8)). A suite of model runs that sweep over bubble number densities are shown in Fig. 5. The data closely track the model for a small range of N b for ϕ ≲ 0.25, but deviate to successively lower N b for ϕ ≳ 0.25, indicating that the model overestimates bubble growth rate for higher ϕ. This implies that some of the model assumptions are violated at higher gas volume fractions.  Table 1 for experimental conditions) which reaches an equilibrium gas volume fraction, ϕ ≈ 0.4. a) Observed data and the best fit model. b) Goodness-of-fit between the model and the data, represented as root mean squared error (RMSE), as a function of C i and N b . The 'bullseye' pattern indicates a clear global minimum in C i and N b , which are therefore well constrained. c) Model outputs for four different combinations of C i and N b (see legend) demonstrate the effect of each parameter on growth rate and equilibrium gas volume fraction.

Modelling of decompression experiments
Each decompression experiment yields data only for the final state of the sample at the end of the experimental run. However, application of the numerical model allows us to reconstruct the evolution of the gas volume fraction throughout the experiment. Fig. 6 shows two modelled ϕ(t) curves for decompression experiment sample ABG27 (see Table 2 for experimental conditions); the equilibrium growth curve is also plotted. The two model runs are identical except that: 1) they assume different values for the bubble number density; and 2) one assumes instant quench, equivalent to stopping the model before cooling, while the other run includes a realistic quench rate, and so captures bubble  Table 1 for experimental conditions). a, c,e,g) Observed data and the best fit model. b,d,f,h) Goodness-of-fit between the model and the data, represented as root mean squared error (RMSE), as a function of C i and N b . Best fit parameters for K06_975 are C i = 0.1106, log 10 N b = 10.326, for K07_1000 are C i = 0.1100, log 10 N b = 10.177, for K11_1050 are C H 2 O = 0.1065, log 10 N b = 10.544, and for K07_1000 are C H 2 O = 0.1035, log 10 N b = 10.544. The equilibrium gas volume fraction is not as well defined by the observed data for these samples it is in sample K05_975 (Fig. 3); this manifests as a best-fit trench of covariation, rather than a 'bullseye' pattern, indicating that the best fit parameters are not as well constrained. For sample K11_1050 (e, f) the data stop well short of the equilibrium gas volume fraction and a similarly good fit is obtained for pairs of values anywhere along the trench of covariation. Sample K10_1050 (g, h) has a low equilibrium gas volume fraction compared with other samples, which results in a broader trench.  Table 1 for experimental conditions) which reaches a relatively high equilibrium gas volume fraction of ϕ ≈ 0.64. The solid lines represent distinct model runs at different bubble number densities. The initial H 2 O concentration for the model runs was estimated from the equilibrium gas volume fraction (see text for details). Note that at ϕ ≳ 0.4, no single model curve captures the data, which cross model curves towards lower bubble number densities. Fig. 6. Application of numerical model to experiment ABG27 (Table 2). Decompression rate is 1 MPa s −1 from initial pressure 100 MPa to final pressure 60 MPa. Black line represents the equilibrium growth curve for the initial H 2 O concentration and dwell temperature. Model runs for two different bubble number densities are shown, bracketing the range of N b observed in the experimental samples (Table 2): 'max. run' with N b = 10 × 10 12 m −3 is stopped at p f , mimicking an instantaneous quench; 'min. run' with N b = 0.4 × 10 12 m −3 includes quench at 100 K/s at p f . The difference between the max. and min. runs for the range of experimental conditions across the suite of samples in Table 2 defines the range of ϕ shown in the shaded bands in Fig. 7.
shrinkage and resorption during isobaric cooling. In both cases, the curves show accelerating growth during isothermal decompression. The run that includes modelling of the quench process shows decelerating shrinkage during isobaric cooling. The shrinkage is driven by a combination of the equation-of-state decrease in gas volume as temperature drops, and by resorption of H 2 O from the bubble, back into the melt, as the solubility increases with decreasing temperature (Eq. (8)). Analytical evidence for thermally-driven resorption in the same suite of samples is presented by McIntosh et al. (2014). Owing to the resorption during quench, the peak gas volume fraction attained by the sample is substantially higher than the final, measured value. The two curves in Fig. 6 show the effect of varying bubble number density and quench cooling rate on the ϕ(t) evolution of a sample. Higher bubble number density leads to faster growth, shifting the modelled ϕ(t) towards the equilibrium curve. As before, growth is faster at higher bubble number density because the distance over which H 2 O must diffuse through the melt to exsolve into the bubble is shorter, and because the shell of viscous melt that provides resistance to growth is thinner. Therefore, the gas volume fraction reached at the final pressure is higher, and nearer the equilibrium value, for higher bubble number density. In general, a faster quench allows less time for the sample to resorb and shrinkthe instant quench run can be thought of as the limiting case of fast coolinghence the final gas volume fraction is higher for a faster quench.
We conducted model runs of the sort shown in Fig. 6 for the range of experimental conditions across the suite of samples in Table 2 (i.e. samples from Burgisser and Gardner (2004), Mourtada-Bonnefoi and Laporte (2004), and Hamada et al. (2010)). Values of bubble number density and cooling rate during quench (instant quench and linear quench at 100 K/s) are bracketed to produce a maximum and minimum curve for final gas volume fraction as a function of final pressure, for three different values of decompression rate (Burgisser and Gardner samples) and for two different temperatures (Mourtada-Bonnefoi and Laporte, and Hamada et al. samples). Model output curves and corresponding experimental data are plotted in Fig. 7. The agreement between model and data is good, with all experimental datapoints falling within the relevant range of model results. Both disequilibrium bubble growth and thermal resorption during quench tend to produce final samples with lower gas volume fractions than would be found under equilibrium conditions. For example, at a decompression rate of 1 MPa s −1 and a final pressure of 60 MPa, accounting for disequilibrium and resorption could decrease final ϕ from the equilibrium value of 0.25, to as little as 0.13 (Fig. 6). The model can be used to deconvolve the effects of disequilibrium and quench resorption, allowing the pre-quench gas volume fraction to be reconstructed for the experimental samples.

Discussion
The results and analysis presented in Section 4 validate the numerical model that we develop in Section 2 over a wide range of experimental conditions. Modelling of the continuous, isobaric experiments (Sections 3.1 and 4.1) constitutes a high-resolution test of the model under strongly disequilibrium conditions in which bubble growth is driven by supersaturation of the melt at constant pressure and varying temperature. Modelling of the decompression experiments (Sections 3.2 and 4.2) tests the model under conditions that are more typical of volcanic and magmatic systems, in which isothermal bubble growth is driven by decreasing pressure, generating varying degrees of disequilibrium. Together, these tests demonstrate that the model can be applied successfully over arbitrary pressure-temperature-time pathways.

Scope and validity
The continuous data demonstrate that the model closely reproduces the ϕ(t) evolution of samples with ϕ ≲ 0. 4 (Figs. 3,4,and 5). For samples that reach equilibrium at a gas volume fraction below this threshold (Figs. 3 and 4) the model captures the whole of the bubble growth, including its approach to equilibrium, with a high degree of fidelity. For the sample that reaches equilibrium at gas volume fraction above this threshold (Fig. 5) the model reproduces the early part of the bubble growth curve (i.e. when ϕ ≲ 0.4), but overestimates bubble growth rate for ϕ ≳ 0.4. This mismatch most likely arises because the model does not explicitly account for bubble-bubble interactions. These interactions become increasingly important as gas volume fraction increases and bubbles begin to impinge on one another, breaking the spherical symmetry that is assumed by the model (Fig. 1). Bubble interaction can lead to coalescence at high gas volume fraction (Blower et al., 2001;Giachetti et al., 2019;Mueller et al., 2005), effectively reducing the bubble number density, leading to slower bubble growth. We note, however, that the model accurately predicts a final gas volume fraction as high as ϕ~0.8 for the decompression experiments (Fig. 7).
The best-fit residual between model and data for each of the samples presented in Figs. 3-5 is plotted in Fig. 8. For all samples, regardless of final gas volume fraction, the residual between model and data is less a) b) Fig. 7. Application of the numerical bubble growth model to decompression experiment datasets (Table 2). a) Experimental data of Burgisser and Gardner (2004) and model outputs for decompression rates of 10, 1 and 0.1 MPa s −1 at 825°C. Black diamond represents a sample that has only undergone the initial decompression (150 MPa to 100 MPa) to nucleate bubbles. Solid line represents the equilibrium growth profile. b) Experimental data of Mourtada-Bonnefoi and Laporte (2004) and Hamada et al. (2010) and model outputs for decompression rate of 1 MPa s −1 at 700 and 800°C. Solid lines represent the equilibrium growth profile for the two temperatures. The coloured/shaded bands (a, b) represent the range in modelled values of final ϕ, resulting from bracketing values of bubble number density and quench rate (Fig. 6). For both a) and b) the range of bubble number densities are N b = 0.4 × 10 12 m −3 to N b = 10 × 10 12 m −3 demarking the minimum and maximum edges of the shaded areas.
than 5% for ϕ ≲ 0.3, and for samples with a final gas volume fraction ϕ ≲ 0.4, the residual is within 2% for the entire bubble growth history. Fig. 8 also shows that there is a slight systematic tendency for the model to underestimate growth rate at low gas volume fraction, and overestimate growth rate at high gas volume fraction. The most likely explanation is that the samples do not contain a population of identical monodisperse bubbles, as is implicitly assumed in the modelling, but a population of bubbles that nucleated at different times, with different initial radius, and with non-uniform spatial distribution. Consequently, the growth of the whole sample is a convolution of the growth of a population of bubbles that are growing at different rates, to different final (local) gas volume fractions. The convolution of many growth curves will tend to produce a global growth curve that is shallower than the growth curve for a best-fit 'average' bubble, which would manifest as a residual with the shape shown in Fig. 8. The implication is that more accurate modelling of the growth of a heterogeneous sample could be achieved by convolving individual growth curves for bubbles with different characteristics, representative of the sample. This could be straightforwardly implemented with the existing numerical code through simple post-processing of cohorts of model runs.
While the model has been validated against samples of rhyolitic composition, it can be applied to silicate melts of arbitrary composition if component models for the relevant material properties (solubility, diffusivity, and viscosity: see Section 2.3) are available. The MATLAB code was designed with ease-of-use in mind, and different models for material properties can be straightforwardly implemented as functions that are called by the main routine. Currently, the model is limited to systems in which H 2 O is the only significant volatile species but including other gases would be straightforward where appropriate models for the equation of state, solubility, and diffusivity are available.
A more significant limitation is that the model assumes perfect closed-system degassingi.e. that a bubble remains in, and only interacts with, the aliquot of melt in which it nucleated. The model is not, therefore, appropriate for systems in which gas is lost from a bubble by escape through a permeable foam, or in which bubbles ascend through the melt with an appreciable slip velocity.

Applications
The principal application of the numerical model is expected to be in forward modelling of bubble growth and resorption in magma in response to changes in pressure and/or temperature. Previous modelling in this area has focussed primarily on the growth of bubbles in response to decompression during ascent of magma that feeds a volcanic eruption (e.g., Liu and Zhang, 2000). The validation of our model for strongly disequilibrium conditions, as well as for conditions under which temperature change is driving bubble growth and resorption, demonstrates that the model can also find application in scenarios that involve more complex pressure-temperature-time pathways. Natural scenarios, of this sort, to which the model could be confidently applied include: the thermal vesiculation of magma undergoing viscous or frictional heating under shear (Lavallée et al., 2015); repeated cycles of decompression, quenching, and heating in pyroclastic material that accretes to conduit walls (Gardner et al., 2018); cyclic pressurization-depressurization of magma subjected to seismic shaking (Manga and Brodsky, 2006); post-eruptive expansion of bubbly pyroclasts (Giachetti et al., 2010;Kaminski and Jaupart, 1997); and post-emplacement cooling-driven resorption in pyroclasts (Watkins et al., 2017).
The model could also be applied in the inverse sense, to reconstruct pressure-temperature-pathways from analytical measurements of volatile distributions around vesicles preserved in volcanic glass, such as those observed by (McIntosh et al., 2014). Similarly, the model could be applied to determine and correct for the effects of resorption during quench to reconstruct bubble sizes and gas volume fractions in natural samples before quench.
Potential applications are not limited to natural scenarios but also include modelling of HPHT experiments. For example, forward modelling could be used to identify suitable experimental conditions to investigate interesting phenomena, or to produce samples with particular textural characteristics; inverse modelling could be used to correct for the effects of resorption during quench.

Conclusions
We describe and validate an easy-to-use MATLAB implementation of a numerical model for the growth and resorption of bubbles in magma undergoing arbitrary pressure-temperature-time pathways. Validation against new, continuous in situ data for natural samples vesiculating at atmospheric pressure demonstrates that the model accurately captures bubble growth for ϕ ≲ 0.4, but indicates that bubble-bubble interactions, which are not captured by this, or any other bubble growth model described in the volcanological literature, become important at higher gas volume fractions. Validation against decompression experiments from the literature demonstrates that the model accurately captures decompressive growth and thermally-driven resorption of bubbles for experiments at widely differing degrees of disequilibrium, for conditions relevant to the syn-eruptive ascent of magma. The model is available via the online supplementary materials.
The modelling approach and implementation is flexible and extendable, opening avenues for future research. Cohorts of model runs for bubbles with differing initial conditions (timing of nucleation, size of nucleus, size of melt shell) could be convolved to capture more accurately the bulk growth of heterogeneous magmas. Constituent models for material properties, such as solubility and diffusivity of H 2 O, and viscosity of the melt, can be changed by the user to facilitate application to a range of magma compositions. Additional volatile phases could be added straightforwardly where suitable solubility and diffusivity models are available. For example, the addition of CO 2 would facilitate application to basalts, where significant CO 2 concentrations are common. The addition of volatile species with solubilities and diffusivities that differ from those of H 2 O (CO 2 , SO 2 , Cl, F, etc.) would allow more accurate and unique inversion of measured gas compositions to reconstruct Fig. 8. Residuals between the gas volume fractions of the best fit models ϕ model and the in situ experimental data ϕ obs  from this study, showing close agreement over the full bubble growth history for samples with final gas volume fraction ϕ ≲ 0.4. There is a slight tendency for the model to underestimate growth (+ve residual) at low ϕ and overestimate (−ve residual) at higher ϕ. The residual is more pronounced for the sample with the highest final gas volume fraction (K_1_1050), see main text for discussion. degassing depths. Similarly, the representation of H 2 O as two interconverting species (molecular and hydroxyl water) using, for example, the model of Coumans et al. (2020), would enable volatile distributions around vesicles to be used to determine decompression pathways with greater fidelity.
10. Compute the viscosity of the melt at each node η(x) using the Giordano et al. (2008) model (Section 2.3.3). 11. Compute the integrated viscosity across the melt shell using trapezoidal numerical integration.
13. Return a vector with to the ODE15s solver with À∂J=∂t and dR=dt as the last value. In this way, ODE15s simultaneously solves ∂C H 2 Ot =∂t (following Eq. (A.4)) for each node and dR=dt throughout the duration of the simulation. The final output is a matrix with C H2Ot (x i , t) and R(t) for any time of potential interest. Please see the downloadable scripts and user manual for more information.