An experimentally-validated numerical model of diﬀusion and speciation of water in rhyolitic silicate melt

The diﬀusion of water through silicate melts is a key process in volcanic systems. Diﬀusion controls the growth of the bubbles that drive volcanic eruptions and determines the evolution of the spatial distribution of dissolved water during and after magma mingling, crystal growth, fracturing and fragmentation, and welding of pyroclasts. Accurate models for water diﬀusion are therefore essential for forward modelling of eruptive behaviour, and for inverse modelling to reconstruct eruptive and post-eruptive history from the spatial distribution of water in eruptive products. Existing models do not include the kinetics of the homogeneous species reaction that interconverts molecular (H 2 O m ) and hydroxyl (OH) water; reaction kinetics are important because ﬁnal species distribution depends on cooling history. Here we develop a ﬂexible 1D numerical model for diﬀusion and speciation of water in silicate melts. We validate the model against FTIR transects of the spatial distribution of molecular, hydroxyl, and total water across diﬀusion-couple experiments of haplogranite composition, run at 800–1200 (cid:1) C and 5 kbar. We adopt a stepwise approach to analysing and modelling the data. First, we use the analytical Sauer-Freise method to determine the eﬀective diﬀusivity of total water D H 2 O t as a function of dissolved water concentration C H 2 O t and temperature T for each experiment and ﬁnd that the dependence of D H 2 O t on C H 2 O t is linear for C H 2 O t K 1 : 8 wt.% and exponential for C H 2 O t J 1 : 8 wt.%. Second, we develop a 1D numerical forward model, using the method of lines, to determine a piece-wise function for D H 2 O t C H 2 O t ; T ð Þ that is globally-minimized against the entire experimental dataset. Third, we extend this numerical model to account for speciation of water and determine globally-minimized functions for diﬀusivity of molecular water D H 2 O m C H 2 O t ; T ð Þ and the equilibrium constant K for the speciation reaction. Our approach includes three key novelties: (1) functions for dif-fusivities of H 2 O t and H 2 O m , and the speciation reaction, are minimized simultaneously against a large experimental dataset, covering a wide range of water concentration (0 : 25 (cid:2) C H 2 O t (cid:2) 7 wt.%) and temperature (800 (cid:3) C (cid:2) T (cid:2) 1200 (cid:3) C), such that the resulting functions are both mutually-consistent and broadly applicable; (2) the minimization allows rigorous and robust analysis of uncertainties such that the accuracy of the functions is quantiﬁed; (3) the model can be straightforwardly used to determine functions for diﬀusivity and speciation for other melt compositions pending suitable diﬀusion-couple experiments. The modelling approach is suitable for both forward and inverse modelling of diﬀusion processes in silicate melts; the model is available as a MATLAB script from the electronic supplementary material.


INTRODUCTION
The diffusion of water through silicate melts is one of the most important processes in volcanic systems. Volcanic eruptions are driven by the buoyancy arising from bubbles that grow when volatiles that are dissolved in the magma diffuse into them (Blower et al., 2001;McBirney and Murase, 1970;Sparks, 1978). Water is the most important volatile component because it is usually the most abundant, and because the dissolved total water concentration C H 2 Ot strongly affects the viscosity of the melt (Hess and Dingwell, 1996). As bubbles grow, water exsolves at the bubble-melt interface, setting up a gradient in its concentration, which drives the diffusion of water through the melt. Diffusion of water also plays a central role in shaping a wide range of other physical and chemical processes in the volcanic system, including welding of pyroclasts in the conduit (Gardner et al., 2018;Wadsworth et al., 2014) and in ignimbrites (Sparks et al., 1999), outgassing of fractured magma (Castro et al., 2012) and pyroclasts (Hort and Gardner, 2000), and bubble resorption during pressurization or cooling (McIntosh et al., 2014;Watkins et al., 2012). Consequently, considerable experimental efforts have been spent in quantifying water diffusivity in silicate melts, and it is known to depend on bulk melt composition, temperature, and C H 2 Ot Wang et al., 2009;Zhang et al., 2017;Zhang and Ni, 2010;Zhang et al., 1991).
An accurate, quantitative model for diffusion is required for both forward and inverse modelling of volcanic processes. Forward models that include diffusion are routinely used to explore a wide range of volcanic and magmatic processes (e.g., Proussevitch and Sahagian, 2005). Inverse modelling, in which the spatial distribution of dissolved water in eruptive products is used to infer their eruptive history, is an emerging area (Ferguson et al., 2016;Humphreys et al., 2008;McIntosh et al., 2014;Watkins et al., 2009), which has the potential to transform our understanding of pre-and syn-eruptive processes. Diffusivity models are often built on data collected from diffusion-couple experiments, in which pieces of glass with different dissolved water contents are placed in contact and held at elevated temperature and pressure, during which time the step profile in water concentration is smoothed by diffusion. Suites of samples quenched after different dwell times are analysed and the spatial distribution of water used to reconstruct diffusivity (Fanara et al., 2013;Nowak and Behrens, 1997;Zhang and Ni, 2010). However, these studies rarely propagate uncertainties through to the resulting diffusivity models. This is particularly problematic for inverse models, because model uncertainty affects the quality and uniqueness of the interpretations.
Reconstructing eruptive histories through inverse modelling is particularly effective when water speciation is accounted for (McIntosh et al., 2014), so accurate models for speciation are also required. Water exists as at least two species when dissolved in silicate melts: molecular water (H 2 O m ) and hydroxyl groups (OH). Molecular water is the diffusing species, whilst OH is essentially immobile, since it is bound to the silicate network (Zhang et al., 1991). Water species interconvert with bridging oxygens O according to an equilibrium speciation reaction where the equilibrium position depends mainly on temperature and dissolved water content. Consequently, analysis of the spatial distributions of the two species, and their relative concentrations, supports the deconvolution of the temperature and pressure histories of a sample (McIntosh et al., 2014). We develop the first 1D numerical model for water diffusion and speciation reaction (including kinetics) in silicate melts, which we validate against data from diffusion-couple experiments. The numerical model includes component models for diffusivity and speciation of water as functions of temperature and dissolved water concentration. We use the numerical model to minimize simultaneously the diffusivity and speciation models from the same experimental data, and carefully propagate experimental uncertainties.
We anticipate three categories of users of the work presented in this article. For convenience, we guide those users to the following sections: (1) Those seeking a validated model for diffusivity of total water with well-constrained confidence intervals are directed to Section 4, and Eqs. (13) and (14) and Table 2.
(2) Those seeking validated and internally-consistent models for diffusivity of molecular water and the equilibrium constant of the speciation reaction, with well-constrained confidence intervals, are directed to Section 5, and Eqs. (23) and (25), and Table 4. (3) Those seeking a flexible numerical model that can be used to forward model diffusion and speciation, to inverse model the spatial distributions of water species in natural or experimental samples, and to design experiments and analyse the resulting samples, are directed to section 5, and to the online electronic supplement, where a downloadable MATLAB implementation is available.

BACKGROUND
The solubility and diffusivity of water have been widely investigated in rhyolitic melt compositions because of its significant impact on the physical properties of the melt. Thorough reviews of the behaviour of water in silicate melts can be found in Zhang and Ni (2010) and . Here we first briefly describe how water is dissolved in rhyolite melts and transported through the silicate network, then summarize modelling methodologies to date and introduce our approach.

Speciation and kinetics
Infrared spectra of hydrous volcanic glasses (Stolper, 1982a,b) indicate that they contain at least two dissolved water species: molecular (H 2 O m ) and hydroxyl (OH). This can be explained by reaction of molecular water with bridg-ing oxygens (O ) within the silicate network to produce hydroxyl groups, according to the speciation reaction (Eq. (1)). The equilibrium constant of the speciation reaction (Eq. (1)) is defined as: where X OH , X H 2 Om , and X O are, respectively, the mole fractions of hydroxyl water, molecular water, and bridging oxygens normalized to a single oxygen basis. Both the equilibrium and kinetics of the speciation reaction are strongly temperature-dependent. Comparison of hightemperature, in situ measurements Behrens, 1995, 2001;Shen and Keppler, 1995) and roomtemperature measurements of rapidly quenched experiments (Ihinger et al., 1999;Zhang, 1994;Zhang et al., 1997bZhang et al., , 1995Zhang et al., 2003) shows that the hightemperature equilibrium speciation can be preserved only if quench occurs rapidly from relatively low temperatures near the glass transition (e.g., <800 K) . Quenching from higher temperatures results in observed speciation equivalent to speciation at temperatures near the glass transition. The equilibrium constant K for the speciation reaction follows an Arrhenius dependence on temperature T and can be described using a symmetric ternary regular solution model (Hui et al., 2008;Ihinger et al., 1999;Silver and Stolper, 1989;Silver et al., 1990;Zhang et al., 1991)  Silver and Stolper (1989) for a full description of the physical meanings of d 1 Á Á Á d 4 . Quantification of the speciation kinetics allows measurements of water speciation to be used to calculate cooling rates (Zhang et al., 1997b. The forward reaction rate k f (rate of consuming H 2 O m to produce OH) is related to the backward reaction rate k b (rate of consuming OH to produce H 2 O m ) via the equilibrium constant (Zhang et al., 1997b) leading to equations for the rate of change of each species (Zhang et al., 1997b) These previous studies of speciation and kinetics enable us to model the speciation reaction at arbitrary total water concentrations, temperature and cooling rates in silicic melts. In principle the same approach can be applied to any melt composition where appropriate experimental data are available.

Diffusion
Early work (e.g., Delaney and Karsten, 1981;Shaw, 1974) identified that the apparent diffusivity of total water H 2 O t is dependent on its concentration. Coupling of diffusion and speciation through Fourier-transform infra-red spectroscopy (FTIR) analysis of diffusion-couple experiments  allows a better understanding of the transport mechanism of water through rhyolite, and shows that the diffusivity of hydroxyl water D OH is much smaller than the diffusivity of molecular water D H 2 Om in rhyolitic melts even when OH is the dominant species (Zhang et al., 1991). This implies that water is transported by rapid diffusion of H 2 O m and undergoes subsequent equilibration with the aluminosilicate network to form OH, via Eq. (1). Overall these data indicate that: (1) at low water concentrations ( K 2 wt.%) the diffusivity of H 2 O t and H 2 O m is approximately proportional to the concentration of total water (Nowak and Behrens, 1997, Fig. 6b therein;; (2) at higher H 2 O t concentrations, H 2 O t and H 2 O m diffusivities increase exponentially with the concentration of total water (Ni and Zhang, 2008;Nowak and Behrens, 1997;Wang et al., 2009); and (3) the apparent diffusivity of H 2 O t can be expressed as a function of the diffusivity of H 2 O m and the speciation if local equilibrium between H 2 O m and OH is maintained (Wang et al., 2009;, leading in first approximation to the expression at time t, assuming that D OH is negligible Diffusivities for hydrous species have also been determined for dacitic (Ni et al., 2009), andesitic (Ni et al., 2013), basaltic (Zhang et al., 2017), and phonolitic/trachytic compositions (Fanara et al., 2013). The diffusivity of hydroxyl water appears to be more significant in mafic compositions, especially at low H 2 O t concentrations (Zhang et al., 2017). The non-negligible D OH in mafic compositions potentially restricts application of Eq. (6) to more silicic compositions.

Previous modelling efforts
Due to the concentration dependence of diffusivity, a straightforward error-function fit is not applicable. There have been two approaches to calculating concentration dependent water diffusivities: (1) Boltzmann-Matano analysis, including the Sauer-Freise method (Fanara et al., 2013;Nowak and Behrens, 1997;Sauer and Freise, 1962); and (2) forward numerical modelling . The Boltzmann-Matano approach determines diffusivity graphically from the shape of the inferred relationship between diffusivity and C H 2 Ot , while forward numerical approaches model diffusion using finite difference or finite element methods, iteratively adjusting input parameters to obtain a good fit to measured concentration profiles. The Boltzmann-Matano approach has the advantage of not requiring prior knowledge of the functional relationship to describe diffusivity, however, this approach can only handle isothermal and isobaric data, and cannot account for speciation, and errors increase as concentration gradients become small. Conversely, numerical methods are more computationally expensive and require prior knowledge of the functional form; however, they can handle speciation and arbitrary pressure-tempera ture-time pathways.
To exploit the strengths of each of these methods, we developed a stepwise approach that applies both methods to a single experimental suite. We first use the Sauer-Freise method (Sauer and Freise, 1962) to determine the functional form of water diffusivity and its dependence on concentration and temperature. Once a functional form is determined, we then use a forward numerical approach to determine more precisely the diffusion relationship. A unique aspect of this study is the addition of the speciation reaction kinetics by incorporating Eqs. (2)-(6) into the numerical framework.

Samples
We use the experimental diffusion-couple samples from Nowak and Behrens (1997). Diffusion-couple experiments are preferred over dehydration experiments for modelling purposes because the boundary conditions are simple (no flux), the initial water concentrations are known, and diffusion profiles are easily modelled with a 1D geometry. Full details of sample synthesis and experimental conditions are given in Nowak and Behrens (1997). Briefly, the starting materials are haplogranitic glass close to the ternary Albite-Orthoclase-Quartz (AOQ) minimum composition (Ab 38 -Or 34 Qz 28 ) at 5 kbar (Holtz et al., 1992). The haplogranitic composition is a widely used synthetic analogue for granite eutectic melts and is, therefore, a suitable analogue for associated rhyolites. Crystal-free glass pieces with different, homogeneous water contents were synthesised and their water concentration determined by Karl-Fischer titration. Glass blocks with differing water contents were placed in contact and held at a range of temperatures (800-1200°C ) at a pressure of 5 kbar for various dwell-times (30-720 min) in an internally heated pressure vessel (IHPV). Run conditions are given in Table 1. Samples were heated linearly to experimental temperature at $100°C/min and were cooled by switching off the furnace, giving a maximum cooling rate of $200°C/min; experimental pressure was maintained throughout quench. Temperature change during heating and cooling was recorded, and we use the data to derive model temperature curves that we use for numerical modelling of the experiments (details in Appendix A). After quenching, $500 lm thick slabs of the experimental charge were cut perpendicular to the diffusion interface, doubly polished, and mounted over a 2 mm wide by 20 mm long slit in a glass slide in preparation for analysis by Fourier Transform Infrared Spectroscopy (FTIR). We used samples from experiments that had dwell times P30 min. Nowak and Behrens (1997) also performed experiments at pressures below 5 kbar, but these were prone to exsolution, and the resulting bubbles complicate FTIR analysis. Consequently, we focus exclusively on the products of experiments at 5 kbar.

FTIR analysis
For this study we use two analytical FTIR datasets from the same experimental samples. The data analysis that forms the core of this study -the minimization of model coefficients against water species -was performed using data collected as part of the Nowak and Behrens (1997) study. They determined concentrations of H 2 O m and OH as functions of position along the diffusion axis (which we define as the x-axis) of the diffusion-couple samples (see Nowak and Behrens (1997) for full details of the analytical Table 1 Diffusion-couple experimental samples and conditions -full details in Nowak and Behrens (1997). Water concentrations are for initial values of high-end/low-end of the diffusion couple.  Nowak and Behrens (1997). † New speciation data acquired for error analysis purposes. methodology; an example dataset is presented in Fig. 2 of that work). We present these data for each sample in the online supplementary material. We also performed new analysis of a subset of the same samples (see online supplementary material), explicitly for the purpose of determining and propagating analytical uncertainties, especially those arising from uncertainty on measured absorbance peak height; Nowak and Behrens (1997) quote an error of 0.003 absorbance units.
We analysed this subset of diffusion-couple experiments using a ThermoFisher Nicolet iN10 spectrometer at the University of Bristol, UK. Individual analytical points were measured along the x-axis with a spacing of 20 lm and along the y-axis with a spacing of 150 lm, using a motorized stage. Spectra were collected using an 80 lm aperture and KBr beamsplitter with a wavelength range of 675 cm À1 to 6000 cm À1 and 3.85 cm À1 spectral resolution. Each analysis consisted of 128 scans, which took $25 seconds in total to acquire. The absorption bands at 5230 cm À1 and 4520 cm À1 were used to obtain peak heights for H 2 O m and OH respectively (Stolper, 1982a,b). The sample thickness was recorded using a micrometer by Nowak and Behrens (1997) with an associated error of ±2 lm. Concentrations were calculated using the Beer-Lambert law: where C is the concentration of species i (H 2 O m or OH) in wt.% H 2 O, W is the molar mass of H 2 O (18.02 g mol À1 ), A is the baseline-corrected height of the species absorption peak (arbitrary units), q is the density of the sample (g l À1 ), which depends on the concentration of water (see Appendix B), d is the thickness of the sample (cm), and e is the molar absorption coefficient of the species (l mol À1 cm À1 ). We use linear molar absorption coefficients determined experimentally for AOQ glass : For consistency with previous workers (Behrens et al., 1996;Nowak and Behrens, 1997) we applied a linear baseline correction by measuring peak height from a tangent to the minima on either side of the H 2 O m band at 5230 cm À1 , which we extrapolate to the OH peak at 4520 cm À1 . We determined the mean and standard deviation of A for the 5230 and 4520 cm À1 peak heights by averaging the spectra collected along the y-axis at each x-position. Methodologies used for propagation of analytical uncertainties are presented in Appendix B.

DIFFUSIVITY OF TOTAL WATER
Some studies of water diffusion neglect the speciation of water in silicate melts to approximate the process as diffusion of total water (e.g., Fanara et al., 2013;Nowak and Behrens, 1997). This greatly simplifies calculations at the expense of losing information that is contained in the relative abundances of the molecular and hydroxyl water spe-cies. Experimental studies have demonstrated that the diffusivity of total water D H 2 Ot is dependent on water concentration and temperature (Fanara et al., 2013;Nowak and Behrens, 1997;Wang et al., 2009;. The 1D diffusion equation for concentration-dependent diffusivity is given by Fick's second law: where t is time in seconds, x is the spatial position along the diffusion axis in metres, and the diffusivity has units m 2 s À1 . We adopt two complementary approaches to determining the diffusivity of total water as a function of concentration and temperature from the diffusion-couple data: the analytical Sauer-Freise method; and direct numerical solution of Eq. (9) using the method of lines. The Sauer-Freise method (Sauer and Freise, 1962) allows the diffusivity at a given water concentration to be calculated without prior knowledge of the functional form of D H 2 Ot C H 2 Ot ð Þ, hence without making assumptions about the physical mechanism of diffusion. A limitation of the method is that it has low accuracy when concentration varies slowly with position, as is the case close to the ends of a diffusion-couple profile. By contrast, numerical methods are more accurate over the entire concentration range, but a functional form for D H 2 Ot C H 2 Ot ð Þ must be assumed. Our approach is therefore to use the Sauer-Freise method to determine a functional form that is then implemented and parameterized numerically.

Total water diffusivity via the Sauer-Freise method
The Sauer-Freise method (Fanara et al., 2013;Nowak and Behrens, 1997;Sauer and Freise, 1962) allows D H 2 Ot C H 2 Ot ð Þ to be determined from data derived from a diffusion-couple experiment, in which C H 2 Ot x ð Þ is a step function at t ¼ 0, via integration of Eq. (9). First, C H 2 Ot x ð Þ data from the diffusion-couple experiment are normalized so that concentration ranges between 0 and 1: A smooth function is then fitted to the normalized water concentration data C Ã H 2 Ot x ð Þ. A Monte Carlo approach is adopted in which fitting is repeated against multiple synthetic datasets in order to account for analytical uncertainty on the raw data; this is the most effective way to robustly propagate analytical uncertainties into the modelling fits, and is facilitated by the short run-time of the computational approach. Synthetic datasets are computed by sampling from a normal distribution using the uncertainty as a function of water concentration as determined by propagation of error through the Beer-Lambert law (Appendix B). We tested polynomial and rational functions of various orders and found that 5th order polynomials give a good fit whilst minimizing high-frequency oscillations, consistent with Fanara et al. (2013) and Nowak and Behrens (1997). The integrated form of Eq. (9) is: where À ½ x indicates that the quantity in brackets is evaluated at position x, and t is the effective duration of the experiment, which is longer than the dwell time because it includes the higher temperature parts of the ramp and quench, see Nowak and Behrens (1997) for details. The D H 2 Ot x ð Þ values computed from the dataset via Eq. (11) are then transformed to D H 2 Ot C H 2 Ot ð Þ via the polynomial function, and Eq. (10). Results for a range of diffusioncouple experiments conducted at different temperatures are shown in Fig. 1.
As expected, diffusivity increases with dissolved water concentration (Zhang and Ni, 2010). The D H 2 Ot C H 2 Ot ð Þ curves show exponential behaviour over the majority of their range (i.e. linear dependence of ln D H 2 Ot on C H 2 Ot , as shown in Fig. 1b), but the low water concentration end of each curve shows a marked deviation from exponential dependence, approximating a linear dependence of D H 2 Ot on C H 2 Ot (most apparent in the inset to Fig. 1a). This is consistent with previous work that has proposed a linear dependence of D H 2 Ot on C H 2 Ot at low water concentrations, and an exponential dependence at higher water concentrations, as determined using numerical methods Wang et al., 2009;. While our results using the Sauer-Freise method support the observations of these previous numerical studies, there are two complications. Firstly, the steepening of the ln D H 2 Ot C H 2 Ot ð Þ curves towards the low water end is seen for all datasets, including those for which the 'low' water concentration end of the diffusion couple is rel-atively high -note the red and dark blue curves in Fig. 1b. Secondly, most datasets also show a slight steepening of the ln D H 2 Ot C H 2 Ot ð Þcurve at the high water concentration end. It is known that the Sauer-Freise method is most accurate when dC H 2 Ot =dx is steep, and is less accurate when it is shallow , and that diffusion-couple experiments tend to have shallow gradients in water concentration near their lowest and highest water concentrations. We therefore hypothesise that the non-exponential behaviour could be an artefact of the Sauer-Freise method. Alternatively, if the linear dependence of D H 2 Ot on C H 2 Ot at low water concentrations is real, the Sauer-Freise method may not be well-suited to quantifying that effect. Later, in Section 4.2, we test these hypotheses using the numerical forward model.
The temperature dependence of diffusivity is shown for three concentrations of total water in Fig. 2. The D H 2 Ot T ð Þ curves follow an Arrhenius relationship of the form where c and m are respectively the y-intercept and gradient of the curves, and T is the temperature in Kelvin. Previous studies have also demonstrated that water diffusivity has a modest dependence on pressure (Ni and Zhang, 2008;Nowak and Behrens, 1997). However, the experimental suite to which we have access does not include enough runs at different pressures (all but three were run at 5 kbar) to permit a robust characterization of pressure-dependent behaviour.
The dependence on water concentration can be included in Eq. (12) by making c and/or m functions of C H 2 Ot . Based on inspection of the shape of the curves in Fig. 1, and on the previous work that identifies a linear-to-exponential Fig. 1. Diffusivity of total water D H2Ot ð Þon linear scale (a) and log scale (b) as a function of concentration of dissolved water C H2Ot ð Þ, derived using the Sauer-Freise method (Section 4.1, Eq. (11)) for diffusion-couple experiments conducted at different temperatures. We chose a representative experiment at each temperature that was conducted using the greatest difference in water concentration. Dark lines show best fit values, lighter-coloured bands indicate the uncertainty (2r) based on propagation of analytical uncertainties using a Monte Carlo approach (Appendix B). We determined the uncertainty as a function of water concentration by re-analysing water speciation in selected samples from Nowak and Behrens (1997) and then applied this to the full experimental dataset. Wang et al., 2009;, we propose a piecewise functional form with a linear dependence on concentration at low water concentration, and an exponential dependence at high water concentration: where a 1 Á Á Á a 6 are empirical coefficients, and C Brk H 2 Ot is the water concentration at which the functional form changes from linear to exponential. We can ensure that the piecewise function is continuous and smooth by requiring the value and gradient of D H 2 Ot at C Brk H 2 Ot to be the same for both the linear and exponential functions. This imposes the following constraints which reduce the number of independent parameters in Eq.
(13) to five: a 1 , a 3 , a 4 , a 6 , and C Brk H 2 Ot . Fitting of Eq. (13) (subject to the constraints in Eq. (14)) to the Sauer-Freise data provides an 'initial guess' for the numerical forward model described in the next section.

Total water diffusivity via numerical forward modelling
Numerical forward modelling of the water distribution in the diffusion-couple samples involves solution of Fick's second law for the case of concentration dependent diffusivity (Eq. (9)). This, in turn, requires a function for D H 2 Ot C H 2 Ot ; T ð Þ to be assumed. Based on the analysis of the Sauer-Freise data, we use the piecewise model presented in Eqs. (13) and (14). We also test an exponential relationship with first order polynomials where b 1 Á Á Á b 4 are empirical coefficients. Forward modelling allows us to test and minimize the D H 2 Ot C H 2 Ot ; T ð Þ models directly against C H 2 Ot x ð Þ data from the diffusioncouple experiments, hence this method is not susceptible to any potential artefacts associated with the Sauer-Freise method, identified in Section 4.1.
Similar to Goudarzi et al. (2016), we solve Eq. (9) using the method of lines in which the spatial derivative is discretized using the finite difference method while the time derivative is maintained as continuous. This results in a set of coupled ordinary differential equations (in time) which are solved using the numerical solver ODE15s (Shampine and Reichelt, 1997) in MATLAB. Details of the numerical approach are given in Appendix C. The forward numerical model is coupled with a least squares minimization algorithm to determine the values of the parameters in the D H 2 Ot C H 2 Ot ; T ð Þ equations (piecewise and exponential, separately) that give the best fit to the data. We use the 'Levenberg-Marquardt' least squares algorithm within the MATLAB optimization and minimization toolkit. We determined an initial guess for the diffusion coefficient parameter values (a 1 , a 3 , a 4 , a 6 , C Brk H 2 Ot , b 1 , b 2 , b 3 , b 4 ) using the Sauer-Freise data (Appendix D). We then calibrate the model to C H 2 Ot x ð Þ experimental data from all the 14 diffusion couples simultaneously, such that derived parameter values represent a best-fit to the whole dataset. The initial spatial distribution of water is taken from Nowak and Behrens (1997) and the evolution of concentration at each spatial node with time is computed; an example is shown in Appendix C. The concentration of water as a function of position at the final time is then compared to the observed data for fitting purposes; an example fit is shown in Fig. 3. The fitting yields best-fit and 5-95% confidence intervals on the values of the parameters in the diffusion models, which are presented in Tables 2 and 3. The misfit between the output of the numerical forward model, using the globally-minimized best-fit parameter values for the piecewise and exponential models, and data from all 14 diffusion-couple experiments, is plotted in Fig. 4. The misfit between the data and the models of Nowak and Behrens (1997), Zhang and Ni (2010), and Fanara et al. (2013) are also shown. The results show that the piecewise model gives a better fit to the data than the exponential model, particularly at low water concentrations. This suggests that the deviation from exponential dependence of diffusivity on water concentration at low concentration is real and not simply an artefact of the Sauer-Freise method. The breakpoint -i.e. the concentration below which the exponential model no longer accurately captures the data -is found to be C Brk H 2 Ot ¼ 1:8 wt. %. Both the piecewise and exponential models are more precisely defined at higher water concentrations (C H 2 Ot > 4 wt.%) than at lower water concentrations (C H 2 Ot < 4 wt.%), as demonstrated by the width of the 1 sigma error bounds. The piecewise function is slightly more  Fig. 1b at 0.5, 2, 3.5, and 5.5 wt.% total water concentration; colours indicate temperature as for Fig. 1. accurate than the Nowak and Behrens (1997) model; this may, in part, be a result of the greater number of fitting parameters in the piecewise function. It is also slightly more accurate than the models of Zhang and Ni (2010) and Fanara et al. (2013), but we note that these were developed using different experimental datasets.

DIFFUSIVITY AND SPECIATION OF WATER
Models for the diffusion of total water are a simplification because they neglect the role played by speciation of water. As presented in Section 2, water in silicate melts exists as two species -molecular water and hydroxyl water -which interconvert via the speciation reaction (Eq. (1)). Both the equilibrium constant and the reaction rate are temperature dependent, hence information about the thermal history of a sample is contained in the relative abundances of the water species, for samples quenched significantly above the glass transition temperature T g (e.g., Nowak and Behrens, 2001). In previous work, water speciation was determined independently and then applied to diffusion experiments in order to develop models of molecular water diffusivity (Zhang, 1999). By developing a model that includes the kinetics of the speciation reaction we can now account for the change in speciation during quench enabling us to minimize for K and D H 2 Ot simultaneously, from the same dataset. In this section, we modify the numerical forward model to capture both diffusion and speciation.
Since hydroxyl water is bound to the silicate structure of the melt and diffuses much more slowly than molecular water for rhyolitic compositions (i.e. D OH ( D H 2 Om ), we can neglect diffusion of hydroxyl water, and treat molecular water as the only diffusing water species Zhang et al., 1991). For mafic compositions, D OH will need to be included in the numerical formulation. The rate of change of concentration of molecular water at any point is therefore given by the combined effect of a diffusive component (transport of molecular water) and a speciation component (conversion of molecular water into hydroxyl species, and vice versa), whereas the rate of change of concentration of hydroxyl water depends only on the speciation reaction.
The change in molecular water concentration owing to the diffusion of molecular water is captured by Fick's second law: @X H 2 Om @t The change in molecular water concentration owing to species interconversion can be derived from Eqs. (2), (4), (5) and (E.2) (in Appendix E, which describes the conversion between concentration and mole fraction) Fig. 3. Numerical forward modelling of the spatial distribution of dissolved water in a typical diffusion-couple experiment using a piecewise (Eq. (13) and (14)) and an exponential model (Eq. (15)) for D H2Ot C H2Ot ; T ð Þ . Model parameters are minimized globally against the dataset of 14 experiments from Nowak and Behrens (1997); parameter values are given in Tables 2 and 3. The light shaded region represents the 2 standard deviation error (2r) as a function of concentration as determined in Appendix B. We have chosen the most conservative error estimates from Appendix B and applied them to the Nowak and Behrens (1997) data for illustrative purposes. Table 2 Best-fit parameters for Eqs. (13) and (14), based on global minimization of a numerical forward model to C H2Ot x ð Þ data from the 14 diffusioncouple experiments.
Piecewise model parameters (Eqs. (13) and (14)  The net rate of change of molecular water concentration is therefore given by: and, from the stoichiometry of Eq. (1), the rate of change of hydroxyl water concentration is given by: In order to solve Eqs. (18) and (19), we need equations for the forward reaction rate k f , the diffusivity of molecular water D H 2 Om , and the equilibrium constant K.
The kinetics of the speciation reaction have been determined empirically via controlled cooling-rate experiments (Zhang et al., 1997b. A model for the forward reaction rate was proposed by Zhang et al. (1997b); in Appen-dix F we extend the validity of that model to a wider range of cooling rates (0.00017-94 K/s) and total water concentrations in the range 0:5 C H 2 Ot 7:7 wt.%, using data from . This yields: where k f has units of s À1 and errors represent the 95% confidence interval of the model fit.
The diffusivity of molecular water is related to the diffusivity of total water for the case where species equilibrium (Eq. (1)) is reached at every point along a concentration profile (Wang et al., 2009) Manipulation of Eqs (2) and (E.2) gives which can be recast as a quadratic equation in X H 2 Om , and solved to give X H 2 Om X H 2 Ot ; K ð Þ . This, in turn, allows the derivative term in Eq. (21) to be determined, yielding an expression for the diffusivity of molecular water: It remains to find an expression for the equilibrium constant -we do this by fitting a numerical forward model to the experimental diffusion couple data Nowak and Behrens (1997). The equilibrium constant follows an Arrhenius dependence on temperature (Nowak and Behrens, 2001;Zhang and Ni, 2010). We consider two functional forms, an ideal mixing model and a non-ideal (regular solution) mixing model, in which the equilibrium constant is a function of the mole fractions of the dissolved water species (Hui et al., 2008;Ihinger et al., 1999;Silver and Stolper, 1989;Silver et al., 1990;Zhang et al., 1991) where the coefficients c 1 , c 2 , and d 1 Á Á Á d 4 must be determined from experimental data. Previous work on rhyolites has indicated that water mixes ideally with silicate melt for C H 2 Ot < 2:5 wt.%, but that a regular solution model is appropriate for higher concentrations (Ihinger et al., 1999).

Numerical forward modelling of diffusion and speciation
The numerical forward model is an extension of the one used in Section 4.2, and details of the implementation are given in Appendix C. We solve Eqs. (18) and (19) to determine the time-evolution of the distribution of molecular and hydroxyl water in the diffusion-couple samples. The model uses Eq. (20) for the forward reaction rate k f , and Eq. (23) for the diffusivity of molecular water D H 2 Om . This, in turn, uses the piecewise model (Eqs. (13) and (14)) for the diffusivity of total water D H 2 Ot , with coefficients minimized against C H 2 Ot using the numerical forward model for total water (Table 2), as described in Section 4.2. The equilibrium constant K is determined using either the ideal mixing model (Eq. (24)) or the regular solution model (Eq. (25)). The only unknowns in the numerical forward model are the coefficients of Eqs. (24) or (25), and these are determined by globally minimizing to the suite of C H 2 Om x ð Þ and C OH x ð Þ data from the diffusion couples. Note that, since all the constituent models include temperature as a variable, the model captures the diffusion-reaction behaviour through the higher temperature parts of the heating ramp and quench process, as well as through the hightemperature experimental dwell.

Ideal mixing
The numerical forward model was used to model profiles for C H 2 Om x ð Þ and C OH x ð Þ for the diffusion-couple samples assuming an ideal mixing model for the equilibrium constant (Eq. (24)); example results are shown in Fig. 5. Note that the final speciation is not representative of the high temperature equilibrium but is 'frozen in' during quench Zhang and Ni, 2010). We sweep through different values of the coefficients c 1 and c 2 and compute the least squares error simultaneously across the entire suite of diffusion couples for which speciation data are available (11 of 14 experiments).
The parameter sweep demonstrates significant covariation between c 1 and c 2 , such that no meaningful global minimum can be found. Rather, there is a 'low-error trench' (Fig. 6) describing a continuum of pairs of coefficient values that give a similarly good fit to the data. Previously suggested values of c 1 and c 2 from studies that determined the equilibrium constant from FTIR measurements of spe- Fig. 5. Numerical forward modelling of the spatial distribution of molecular and hydroxyl water in the diffusion-couple experiments assuming ideal mixing (Eq. (24)) and the piecewise diffusion equation. Fits to two representative datasets are presented. The minimized coefficients of Eq. (24) are c 1 ¼ 26:6 and c 2 ¼ À4:179, which are similar to those presented in Nowak and Behrens (2001). The light shaded regions represent the 2 standard deviation analytical error (2r) error as a function of concentration as determined in Appendix B. cies concentrations in quenched samples (Behrens and Nowak, 2003;Zhang et al., 1997a) and in situ at magmatic temperatures (Nowak and Behrens, 2001) fall within this low-error trench (Fig. 6). This demonstrates that: (1) the numerical forward modelling approach that we adopt yields a K T ð Þ relationship that is consistent with previous work; and (2) published models of the form of Eq. (24) perform similarly well for the current dataset. Fig. 5 shows that the numerical forward model with ideal mixing tends to underestimate C OH at intermediate C H 2 Ot . This suggests that a simple ideal mixing model is insufficient to describe the data. We therefore focus on the 'regular solution' as the simplest extension of ideal mixing, where the entropy of mixing is identical to the ideal mixing solution, but all configurations of the three species (H 2 O m , OH, O ) are not energetically equivalent (Hui et al., 2008;Ihinger et al., 1999;Silver and Stolper, 1989).

Non-ideal mixing
As for the ideal mixing model, we modelled C H 2 Om x ð Þ and C OH x ð Þ for the diffusion-couple samples, using the regular solution model to describe the equilibrium constant (Eq. (25)). The coefficients of Eq. (25) were globally minimized against the data using a non-linear fitting algorithm in MATLAB. We use a fixed value for d 2 of À2.87 K (Ihinger et al., 1999) to reduce the number of free parameters from 4 to 3 and thus minimize the other coefficients, which are presented in Table 4. We chose this parameter because 1000Rd 2 represents the standard state enthalpy change for the reaction among species when H 2 O t approaches 0 (Ihinger et al., 1999), where R is the gas con- Fig. 6. Sum of squared errors between the ideal mixing model and data as a function of the coefficients of Eq. (24), c 1 and c 2 . There is considerable covariance between the coefficients, such that a meaningful global minimum cannot be determined. Previous estimates of the coefficients fall within a 'low-error trench' of coefficient values: Z1997 = Zhang et al. (1997a), NB2001, 2003= Nowak and Behrens (2001 and Behrens and Nowak (2003). Table 4 Globally-minimized best-fit parameters for Eq. (25) (non-ideal mixing/regular solution model), based on fitting the numerical forward model to diffusion-couple data.  stant. Representative results are presented in Fig. 7 and show a significant improvement in the fit to C OH x ð Þ data. We also determined the values for the equilibrium constant at 'apparent equilibrium' K ae ð Þ using the 'frozen-in' speciation (Eq. (2)) from the full suite of samples for which speciation data are available, and compared these with the values computed numerically using the ideal mixing model and the regular solution (Fig. 8). The temperature of 'apparent equilibrium' T ae ð Þ is calculated using the geospeedometers of Zhang et al. (1997b) and  (Appendix F). These results show a significant improvement in fit for the regular solution model, suggesting that non-ideal mixing is occurring (although we note that this study was not specifically designed to investigate melt structure or species interactions).

IMPLICATIONS AND CONCLUSIONS
We have developed a multi-tiered approach to investigate speciation and diffusion of water in a synthetic rhyolite analogue (AOQ) using a combination of diffusion-couple experiments and computational analysis. The results include: (1) a diffusion model for total water D H 2 Ot C H 2 Ot ð Þ that explicitly accounts for the change in behaviour from a linear to exponential dependence as concentration increases; (2) a modelling framework that uses both total water diffusivity and the speciation kinetics to give an internally consistent minimization for the equilibrium coefficient K; and (3) a flexible numerical model that can handle both speciation and diffusion of water for arbitrary P-T-t pathways. In the following sections we discuss the implications of these results.

The piece-wise model yields improved accuracy
Because the dependence of D H 2 Ot on water concentration is linear at low C H 2 Ot  but becomes exponential at higher C H 2 Ot , we use a piecewise function with both linear and exponential components (Eq. (13) and (14)) to model the experimental suite of diffusion-couple samples of Nowak and Behrens (1997). A comparison between the piecewise diffusion model and the exponential model of Nowak and Behrens (1997) show a good agreement in calculated diffusivities (Fig. 9), particularly at higher water concentration, but at lower water concentration and higher temperature the models deviate strongly. For example, at 1200°C and $0.25 wt.% the difference in diffusivity between the two models is approximately one order of magnitude. For the range of temperature, water concentration and experimental conditions (Table 1) investigated in the Nowak and Behrens (1997) dataset, our new piecewise equation provides a small increase in accuracy as demonstrated by the direct comparison of misfits between models and data (Fig. 4). The improved accuracy is most marked at low water concentrations, presumably because the original experimental dataset is heavily weighted towards elevated water concentrations (e.g., >1 wt.%) which therefore dominate the fit of the exponential model adopted by Nowak and Behrens (1997). This highlights the importance of including explicitly the linear behaviour for applications modelling water diffusion in volcanic scenarios where water concentrations are low, such as Fig. 8. The temperature of apparent equilibrium T ae is defined as the calculated 'equilibrium' temperature of the quenched speciation. Other definitions include the closure temperature, the fictive temperature, or the glass transition temperature of the speciation reaction (Zhang, 1994). T ae can be thought of as the apparent temperature at which the final speciation is 'frozen in' or quenched. We computed the temperature of apparent equilibrium T ae for the whole suite of samples for which speciation data are available by re-arranging Eqs. (F.2) and (F.3) in Appendix F, using our expression for k f (Eq. (F.5)), and the expression for quench rate as a function of temperature (Appendix A, Fig. A.1a). The equilibrium constant at apparent equilibrium K ae is computed using Eq.
(2): for each of the calculated T ae from the speciation data we compute a model K ae value using the Ideal mixing and Regular solutions along with the speciation data. The purpose of this figure is to demonstrate that the observed final speciation can be modelled well by the regular solution, supporting previous interpretations of non-ideal mixing (Hui et al., 2008;Ihinger et al., 1999). Note that here, the quench history (dT =dt) is the same for all samples (Appendix A), so T ae is effectively a proxy for total water concentration. The scatter observed at higher T ae (low 1000=T ae ) could be due to the increased relative error in the species as total water, hence also 1000=T ae decreases. Fig. 9. Plot of total water, D H2Ot as a function of C H2Ot and temperature for the piecewise diffusion Eq. (this study, solid lines), for comparison with the exponential model of Nowak and Behrens (1997), and the diffusion Eq. for rhyolite (see Eq. (14) in Zhang and Ni, 2010). in the upper conduit, following emplacement of obsidian flows, or within welded ignimbrite deposits.
The piecewise function explicitly includes a breakpoint where the behaviour changes from linear to exponential. For the haplogranitic system investigated we find a best fit value of the breakpoint at $1.8 wt.% C H 2 Ot , which is consistent with observations from previous studies (Ni and Zhang, 2008;Wang et al., 2009;. A key question for future studies is whether this breakpoint marks a significant structural change that affects the diffusion behaviour in silicate melts, as previously suggested (e.g., Behrens and Nowak, 1997;Persikov et al., 2014) and if so, whether the location of the breakpoint is dependent on the melt composition.

Constraints on equilibrium constant, K
A key innovation of this study is the inclusion of both speciation and kinetics into the numerical framework. The addition of speciation to diffusion models enables diffusioncouple experiments to be used to estimate the equilibrium constant K given a known T-t history (including the quench). We demonstrated this by minimizing for the functional form of the equilibrium constant K for the well-studied haplogranitic system. Our approach is ready for application to other systems (e.g. dacite, andesite and peralkaline melts) when data are available, and has the advantage that the equilibrium constant for speciation can be determined on the same sample and at the same experimental conditions as the diffusion coefficients. This internal consistency will help to minimise uncertainties. Although this study only examined samples at a fixed pressure, this approach can be easily applied to experiments at a variety of pressures. One caveat is that for more depolymerised compositions (andesite and basalt) diffusion of OH becomes non-negligible (Zhang and Ni, 2010;Zhang et al., 1991) and may need to be included in the model. Ideally, experiments would also be performed at different cooling rates to investigate whether the forward reaction rate k f could be obtained by appropriate minimization of the model. This would create a complete and internally consistent set of models from a single dataset rather than building consecutive models from different sets of experiments and compositions, as is current practice.

Limitations
The approach and physical/chemical parameters that we present here have important limitations based upon the experimental data and functional forms used. The experimental data used for determining the rate of forward reaction k f ranged from 400°C to 1100°C, C H 2 Ot concentrations from 0.5 to 7.7 wt.%, and cooling rates from 0.00017 to 94°C/s (Zhang et al., 1997b. The experimental data used for determining D H 2 Ot ranged from 800°C to 1200°C and 0.01 to 7.75 wt.% H 2 O t at 5 kbar. Therefore, we have extrapolated the diffusivity function to lower temperatures during quench. However, we note that the diffusivity is very slow at these temperatures, and the time spent below the calibrated temperature window is short, so our extrapolation has a minimal effect on the results. Where we have validated the model, it can be confidently applied. Extrapolation outside of this range may be reasonable depending on the conditions for example, the diffusion functional form used in this study has been applied to lower temperatures (400-550°C) in other studies (Zhang et al., 1991). The model presented here considers disequilibrium processes via the forward reaction rate and uses the complete quench history, but we do not consider the diffusivity of OH, which can be important for compositions other than rhyolite (Zhang and Ni, 2010;Zhang et al., 1991). Low temperature hydration of obsidian may occur by different mechanisms to high-temperature water diffusion (e.g. Anovitz et al., 2008), so use of our model in these environments is not recommended. Additionally, our diffusivity equation is only valid for pressures of 5 kbar. However, the numerical framework and approach provided here provides a means to expand the calibration parameter space to develop a broader suite of internally consistent models for volcanological applications.

Volcanological applications of the new model
We anticipate that the model will be applied to two categories of study: (1) inverse modelling of experimental data -as in the current manuscript -in order to determine diffusivity, equilibrium constant, and speciation kinetics from experimental data; (2) forward modelling of diffusion and speciation in natural volcanic and magmatic scenarios. Inverse modelling, in the most general case, requires only that the experimental conditions and material composition are accurately known. It is hoped that the availability of the model will encourage researchers to design and execute experimental campaigns to quantify diffusion and speciation for a range of natural melt compositions. Forward modelling using the full speciation and diffusion model requires that the following component models for the material of interest are known: equilibrium coefficient K X; T ð Þ, where X represents mole fractions of the relevant water species; forward reaction rate k f X ; T ð Þ; and either the diffusivity of molecular water D H 2 Om C H 2 Ot ; T ð Þ , or the diffusivity of total water D H 2 Ot C H 2 Ot ; T ð Þ . Note that, in the general case, each of these component models may also be a function of pressure.
The most obvious practical volcanological application of our new numerical framework is the use of precise forward numerical modelling of speciation and diffusion to invert measurements on natural samples to determine their eruptive P-T-t pathways. During magma ascent, decreasing water solubility causes water to diffuse into bubbles at a rate that depends on the evolving water concentration (and temperature for non-isothermal ascent). During and after eruption, cooling causes the solubility of water in the melt to increase (resulting in diffusion back into the melt, e.g., McIntosh et al., 2014) and simultaneously alters the equilibrium speciation towards more H 2 O m -rich compositions (e.g. Nowak and Behrens, 2001). Therefore, coupling the diffusion and speciation model presented here with a bubble growth model (e.g., Blower et al., 2001;Proussevitch and Sahagian, 1998;Proussevitch et al., 1993), would enable a robust investigation of the P-T-t his-tory of magma during the final stages of eruption and emplacement. This would provide a new approach to quantitative understanding of syn-eruptive magma decompression and cooling that could be used to reconstruct conduit processes during volcanic eruptions, including volcanological phenomena such as transitions in eruptive style, with application to silicic volcanism around the world.
FUNDING JPC, EWL, and MCSH acknowledge support from the UK Natural Environment Research Council (NERC) via grant NE/N002954/1.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  A1. (a) The cooling rate (dT =dt) versus the difference between temperature (T ) and ambient temperature (T a ) for runs with initial temperature T 0 ¼ 1200°C (orange circles) and T 0 ¼ 800°C (green circles). Linear regression of the data gives cooling constant h ¼ À0:357 s À1 . (b) The temperature versus time data during cooling from 1200°C and 800°C collected during the experimental runs of Nowak and Behrens (1997). The black line represents the regressed cooling model (Eq. (A.1)) which defines the T t ð Þ cooling function for use in the numerical model. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) This is a quadratic equation in C H 2 Ot and is solved analytically using the quadratic formula. Uncertainties on the quantities in Eq. (B.2) are propagated through to determine the uncertainty on C H 2 Ot using a Monte Carlo approach, described later in this appendix. The value of q, and the associated uncertainty, is computed for each C H 2 Om , C OH data pair via Eq. (B.1), Eq. (7) is then used to compute the values of C H 2 Om and C OH and their associated uncertainties for each measurement point. Uncertainties are propagated through Eqs. (B.1) and (7) according to standard error formulae: where dS is the uncertainty on quantity S, da is the uncertainty on quantity a, etc. These formulae are predicated on the assumption that uncertainties are normally distributed and are uncorrelated. We note that the uncertainty on concentration varies with absolute concentration: uncertainty at low water concentration (C H 2 Om K 1 wt.%) is typically 10-15% relative, but <5% relative at higher water concentrations (Fig. B1). Of the six samples that were re-analysed (6 of 14), the most conservative estimate of uncertainty comes from sample AOQD013, which covers a wide range in C H 2 Ot . The fractional uncertainty (2r C i ð Þ=C i , where C i is the concentration of species i and r is its uncertainty) for each datapoint in the Nowak and Behrens (1997) dataset was assumed to be the same as the fractional uncertainty determined from sample AOQD013 at the same species concentration, interpolated as necessary.

Monte Carlo method for uncertainty propagation
We use a Monte Carlo approach to propagate the analytical uncertainties on the quantities in Eq. (B.2) through to determine the uncertainty on C H 2 Ot . Note that we adopt this approach because the quadratic formula contains the same term more than once, violating the assumption of uncorrelated errors in the standard error formulae (Eq. (B.3)). For each analytical point, we generate 10,000 different synthetic values of the relevant measured quantities, drawn from a normal distribution about the mean measured value, with standard deviation taken as the uncertainty: e.g., 10,000 different values of the absorption Fig. B1. Uncertainty of the measured species concentrations and how they relate to the absolute value of the concentration. Note that relative uncertainty rises as the absolute concentration decreases. coefficient e H 2 Om are drawn from a normal distribution with mean 1.79 and standard deviation 0.02, etc. The value of C H 2 Ot is computed for each of the 10,000 synthetic sets of values, and the uncertainty on C H 2 Ot is taken as the standard deviation of the resulting distribution.
We adopt a similar approach to propagating the uncertainties in the diffusion-couple data through subsequent analysis. Each time a function or model is to be fitted to experimental data -i.e., C H 2 Ot x ð Þ, C H 2 Om x ð Þ, and C OH x ð Þ -we produce between 2000 and 10,000 synthetic datasets in which the concentration at any x-position is drawn from a normal distribution with mean equal to the measured value at that point, and standard deviation given by the associated uncertainty. The function or model is then fitted to each synthetic dataset in turn, resulting in a distribution of best-fit or minimized parameters. The distribution is then used to determine the mean and standard deviation on the best fit values or minimized parameters. Quantities that are derived in this way are presented as the mean value, and the standard deviation is taken as the associated uncertainty.

APPENDIX C. NUMERICAL FORWARD MODEL
The numerical forward model solves Fick's second law (Eq. (9) of the main text) using the method of lines, via MATLAB's ODE solver ODE15s. Eq. (9) is repeated here for convenience: The spatial dimension x is discretized into a number n of equally-spaced blocks along the length of the diffusion couple, which is composed of a 'wet' glass part of length L w and a 'dry' glass part of length L d . The finite difference grid is block-centred, such that the ODEs are solved at the midpoint of each of the blocks -the 'nodes' -which have position À L w ðC:2Þ where i ¼ 1 Á Á Á n, such that x ¼ 0 marks the boundary between the wet and dry glasses. The concentration of total water at each node is represented by a value in wt. % in a 1 Â n matrix: C H 2 Ot x i ð Þ. For nodes with x i < 0, the initial concentration of dissolved water is set to the 'wet' glass value; for nodes with x i > 0, it is set to the 'dry' glass value. The outer boundaries have a Neumann boundary condition, such that the diffusive flux of water is zero (i.e. D H 2 Ot @C H 2 Ot =@x ¼ 0 at x 1 and x n ).
The function that is called by the solver ODE15s is represented by the following pseudo-code: 1. Determine the temperature at the current timestep (Appendix A). 2. Calculate the diffusivity of total water D H 2 Ot at each node using Eqs. (13) and (14) (piecewise model) or Eq. (15) (exponential model), based on current temperature, and the concentration of total water at each node. 3. Determine spatial gradient in water concentration @C H 2 Ot =@x at each node using the diff() function on the concentration array C H 2 Ot x i ð Þ. 4. Determine the diffusive flux at each node J ¼ ÀD H 2 Ot @C H 2 Ot =@x; set J ¼ 0 at x 1 and x n . 5. Determine the gradient of the diffusive flux @J =@x at each node using the diff() function on the flux array J x i ð Þ; this is equal to the À@C H 2 Ot =@t at each node, according to Fick's second law.
In this way, ODE15s solves @C H 2 Ot =@t for each node throughout the duration of the diffusion-couple experiment. The final output is the concentration matrix C H 2 Ot x i ð Þ at the end of the experiment, but note that the spatial distribution of water is computed throughout the experimental runtime; an example is shown in Fig. C1, in which the evolution of the spatial distribution of total water over time