VolcDeGas : A program for modelling hydrogen isotope fractionation during degassing of rhyolitic melts

Magma degassing mechanisms are key determinants of explosive and e ﬀ usive eruption styles. Paired measurements of H 2 O content and hydrogen isotopic ratios (e.g. δ D) in pyroclastic and e ﬀ usive products can elucidate end-member degassing mechanisms (e.g. closed and open system) during eruption. Here we present “ VolcDeGas ”, a MatLab program that models δ D-H 2 O degassing trajectories of rhyolitic magma. Operating within an intuitive GUI, VolcDeGas calculates degassing paths based on: initial δ D (in % (cid:24) ), the H 2 O content of the melt (wt.%), degassing step size, and temperature. VolcDeGas also calculates hydrous speciation based on either empirical models or analytical data, and incorporates this e ﬀ ect in simulations. VolcDeGas solves open-, closed-, and batched-system fractionation equations, and also combines these into multi-stage (e.g. closed-to-open) degassing paths. Tests show that degassing is highly sensitive to step size, and multi-stage scenarios involving progressively declining steps (e.g. a closed- to batched-system) provide the best statistical ﬁt to pyroclastic and e ﬀ usive obsidian geochemical datasets.


Natural context: rhyolitic volcanism
Much volcanological research over the last decades has been focused on understanding the mechanisms driving silica-rich rhyolite eruptions [e.g. Alfano et al. 2011;Castro and Dingwell 2009;Dingwell 1996;Melson and Saenz 1973;Ninkovich et al. 1978;Ogburn et al. 2015;Wilson 1976]. Moderately sized rhyolite eruptions (<1 km 3 ), while rarely directly observed, often start with an explosive phase that involves fragmentation and the formation of pyroclastic columns [Castro and Dingwell 2009]. Sometime later the eruption transitions to an effusive phase whereby lava is emitted [Castro et al. 2012;Tuffen et al. 2013]. This sequence is termed an "explosive-to-effusive" eruption. Specifically how the transition between eruption styles occurs-e.g. how long it lasts and what triggers ithas long been a target of investigations [Eichelberger et al. 1986], as it bears on the way hazards are understood and eventually mitigated. Prior to the first direct observations of rhyolite activity [e.g. Castro and Dingwell 2009], it was assumed that these eruptions transitioned rapidly from explosive to effusive behaviour, this interpretation being based on the conformable superposition of obsidian lavas on co-eruptive pyroclastic fall deposits [e.g. Newman et al. 1988;Taylor et al. 1983]. Hydrous geochemical patterns (i.e. δD versus H 2 O trajectories) measured in the deposits of prehistoric eruptions were interpreted to reflect a two-staged Corresponding author: castroj@uni-mainz.de degassing sequence comprising initial closed-and later open-system degassing [e.g. Dobson et al. 1989;Eichelberger et al. 1986;Newman et al. 1988] and corresponding to the purported pyroclastic to effusive activity. Thus, the phrase "explosive-to-effusive transition" was meant to encapsulate both the physical and geochemical aspects of an eruption sequence, even though these eruptions remained unobserved. Detailed observations and geochemical measurements of the 2008 eruption of Chaitén, Chile [Castro et al. 2014] showed that rhyolite eruptions may not follow the previously accepted explosive-to-effusive paradigm. Notably, this eruption lacked a sharp explosive-to-effusive transition and was instead characterised by "hybrid" activitysimultaneous pyroclast and lava venting-during most of its duration [e.g. see also Castro et al. 2013]. The hydrous geochemical signature of Chaitén's deposits, is furthermore, consistent with a one-staged degassing mechanism that lacked a transition (to be discussed in further detail in Subsection 1.3). In this work we present VolcDeGas, a tool that can model both one-and two-staged degassing paths so as to account for the variety of measured hydrous geochemical patterns in eruption deposits.
Explosive and effusive eruption styles are influenced by the budget and behaviour of volatile components in the melt, which are dominantly H 2 O originating from the magma storage region [e.g. Castro et al. 2014;Castro and Dingwell 2009;Luhr 2001]. Due to the strong dependence of H 2 O solubility on pressure, H 2 O-rich bubbles will grow in the magma as it decompresses during ascent [Allard et al. 1991;Blank et al. 1993;

VolcDeGas
Walter & Castro, 2020Eichelberger et al. 1986Ryan et al. 2015]. Depending on many complicated trade-offs between magma porosity and permeability, and magma ascent dynamics, the magma may begin to violently fragment [Gonnermann and Manga 2007]. Experimental and analytical evidence suggests that high initial magmatic H 2 O contents lead to explosive fragmentation and lower initial magmatic H 2 O contents prompt coherent lava flows [Forte and Castro 2019]. Other variables such as temperature and decompression rate are also important for how degassing takes place [e.g. Forte and Castro 2019;Gonnermann and Manga 2007]. In our model, we address a range of H 2 O and temperature conditions on resultant degassing paths in rhyolite magma.
Hydrogen isotope ratios (e.g. the deuterium/hydrogen ratio), generally cast in terms of the variable δD, is calculated as the difference between ratio of the measured H-isotopes, D/H, in the unknown sample and a standard (e.g. Standard Mean Ocean Water: SMOW) normalised by the D/H ratio of a standard: δD "`2 H{ 1 H˘S ample´`2 H{ 1 H˘S tandard p 2 H{ 1 Hq Standardˆ1 000% (1) Such isotopic ratios are routinely measured in concert with magmatic H 2 O to synoptically track degassing processes that give rise to pyroclastic-effusive sequences [Taylor et al. 1983]. Traditionally, values of δD are plotted against the bulk-H 2 O content measured on obsidians, and the resulting plots or "degassing trends" interpreted in terms of end-member degassing mechanisms ( Figure 1).
Because deuterium fractionates preferably into volatile phase [Taylor et al. 1983], magmatic degassing results in isotopic "lightening" of the residual melt, leading to trends pyroclastic sample arrays that are characterised by subtle to sharp decreases in δD with decreasing H 2 O concentration [e.g. Castro et al. 2014;Taylor 1991]. Degassing trends will reach a minimum at the most-degassed levels, and often, lava samples' δD demarcate a strong final downturn in the degassing path ( Figure 1). This pattern has been observed at numerous rhyolite systems and has been widely interpreted as a signature "open-system" degassing [e.g. Castro et al. 2014;Eichelberger et al. 1986;Newman et al. 1988;Taylor et al. 1983].

H 2 O speciation in rhyolitic melts
H 2 O in silicate melt exists as hydroxyl groups (OH´) and molecular water (H 2 O), with the relative amounts of these hydrous species being controlled by the bulk H 2 O content and temperature [Stolper 1982]. Although these two species exist, it is the molecular H 2 O that diffuses out of the melt and into bubbles [Zhang et al. 1991]. Since deuterium is preferentially fractionated into the vapour, the degree of fractionation (and path Figure 1: Qualitative trends of hydrogen isotope fractionation occurring during different degassing modes. In each case, the effect of degassing-i.e. progressive separation of H 2 O from melt-results in a decrease δD.
[A] Closed-system degassing results in a relatively subtle near-linear drop in δD.
[B] Open-system degassing (Rayleigh fractionation) results in a non-linear yet significant depletion in δD at low H 2 O contents.
[C] Batched-degassing, a hybrid mechanism comprising both closed and open-system elements, involves repetitive degassing "steps", the size of which dictate the ultimate degree of δD decrease as degassing progresses to low H 2 O contents.
of δD) depends on how the species ratios change with falling H 2 O concentration. A hydrogen isotope fractionation factor between vapour and melt in rhyolitic melts, alpha (α v-m ), accounts for this speciation dependence on bulk H 2 O content in modelling degassing trends [Newman et al. 1988]. Thus, to completely retrace the degassing history several variables need to be known or estimated: α v-m , temperature of the melt, initial bulk H 2 O content and initial δD.
Early studies performed hydrogen fractionation calculations with a constant α v-m factor [Taylor et al. 1983]; these factors were experimentally determined [Dobson et al. 1989;Newman et al. 1988;Taylor et al. 1983]. However, as the hydrous speciation of the melt (i.e. the relative proportions of OH´and H 2 O species) is dependent on both the total H 2 O content and temperature, both of which may change during an eruption, alpha factors are not likely to be constant. Accordingly, later studies used hydrous-speciation dependent α v-m factors.
In addition to syn-eruptive degassing, the uptake of externally derived meteoric water during rehydration of obsidian may alter the hydrogen isotopic signature of eruption products [e.g. DeGroat-Nelson et al. 2001;Seligman et al. 2018]. While the degree to which original magmatic isotope signatures are altered depends on both time and the isotopic composition of the external water, significant departures of δD-H 2 O from idealised magmatic degassing trends have been observed [e.g. Bindeman and Lowenstern 2016;Giachetti and Gonnermann 2013;Seligman et al. 2018]. This warrants independent control on the magmatic versus meteoric signature of the isotopic and bulk H 2 O signatures of eruption deposits.

Observations of rhyolitic eruptions
The 2008 eruption of Chaitén volcano, Chile [Alfano et al. 2011;Castro and Dingwell 2009], provided many new insights into the degassing systematics and eruption mechanisms of rhyolitic systems. This was the first opportunity to see how the actual "explosive-toeffusive" sequence came about, and given close observations [Lara 2009] it became clear that the previously assumed temporally sharp changeover of explosive-toeffusive activity [e.g. Eichelberger et al. 1986] was not the case [Castro et al. 2012]. Instead, after a period of about two weeks of purely explosive activity, a simultaneous lava and pyroclastic eruption occurred and lasted several months before solely lava effusion closed out the eruptive episode. This simultaneous, or "hybrid" eruption style [e.g. Castro et al. 2013], involved intermittent pyroclastic explosions from several vents stationed within the effusing lava itself . As observations of the Chaitén event may shed new light on how other pre-historic eruptions like Mono Craters, Big Glass Mountain, and Big Obsidian Flow (USA) may have behaved [e.g. Castro et al. 2014;Rust and Cashman 2007], here we aim to use geochemical characteristics of the eruption deposits of Chaitén and Mono Craters to model aspects of rhyolite eruption systematics. This work presents a new program that can be used to investigate how eruption styles relate to degassing mechanisms ( Figure 1) with particular emphasis on H-isotope and H 2 O content systematics in eruption deposits.

Modelling theory and degassing path calculations
Hydrogen-isotope fractionation patterns in rhyolitic magmatic systems have been explained in a relatively small number of ways ( Figure 1): 1) by a closed system; a simple conservative mass balance between melt and vapour; 2), by an open system involving Rayleigh fractionation, and 3) as a batched-system [Taylor 1991] comprising combinations of closed-and open-system degassing. Modelling these isotopic and degassing patterns assumes that the outgassed part of the volatile system was in equilibrium with the remaining parts of the system (i.e. the quenched silicate glass). Analysis of eruption deposits, therefore, track the remaining mass of H 2 O and H-isotopic composition of the glassy pyroclasts rather than the once-present vapour. Addition-ally, in all model systems, the step size-defined as a fixed amount or the percentage of bulk water loss at each degassing interval-is an initial condition and is required for calculating each point in the degassing history. In the following sections we describe each system in detail including their respective H-isotopic fractionation equations.
1.4.1 Closed-system degassing Closed-system degassing takes place when volatile components begin to exsolve into vapour bubbles [e.g. Anderson and Fink 1989;DeGroat 1997] and may continue to shallow (hundreds of meters) levels in a feeder conduit, potentially driving explosive fragmentation [e.g. Eichelberger et al. 1986;Newman et al. 1988]. In this scenario, melt and fluid and/or vapour stay in contact and it is assumed that compositional equilibrium is maintained between these phases. As no material is added or released from the magmatic system, these conditions comprise a mass-conservative system, and consequently the isotopic composition of hydrogen in the melt is a strong function of that in the bubbles. This "equilibrium buffered state" whereby the ambient pressure controls the amount of H 2 O exsolutionand therefore the hydrogen fractionation degree-is not to be confused with buffered states within an opensystem, whereby the gas-rock ratio is so high that equilibrium is maintained between melt and the permeating gas [e.g. Rust et al. 2004]. Closed-system conditions can however give way to open-system degassing when a magma reaches the critical vesicularity to stimulate effective melt-gas separation via permeable gas flow [Eichelberger et al. 1986].
Isotopic fractionation in the closed-system is modelled by the mass balance equation [Taylor et al. 1983]: where δD f is the final δD value (in % ) in the melt, δD i is the initial δD value of the melt, F is the fraction of H 2 O remaining dissolved in the melt, calculated as final H 2 O at the end of the process divided by the initial H 2 O content in the melt at the start of the process (H 2 O f /H 2 O i ), and α v-m is the bulk H-isotope fractionation factor between vapour and melt at specific temperature and pressure conditions. Here, a "process" is defined as a cycle or a discrete degassing "step" which releases a fixed aliquot (e.g. a percent of the total H 2 O budget) of vapour from the melt. Furthermore, in terms of modelling the closed system, the initial δD and initial bulk H 2 O content are thought to be constant through the process. The isotopic shift in a perfectly closed system follows a nearly linear trend in wt.% H 2 O-δD space and is relatively subtle in terms of the amount of fractionation that occurs with degassing. In other words, a large drop in bulk H 2 O content (several wt.%) results in only slight (few tens of % ) decrease in δD.

Open-system degassing
During open-system degassing, volatiles are continually exsolved from melt into bubbles and, due to high magma porosity and permeability, these volatiles are immediately removed from the bubbly melt through pathways (e.g. magmatic foam, cracks, or the brecciated conduit walls; Eichelberger et al. [1986]; Gonnermann and Manga [2003]). This highly effective degassing mechanism drives very strong isotopic fractionation, as equilibration between dissolved and exsolved H 2 O is prevented by the constant removal of gas. Open system isotope fractionation is described by the Rayleigh fractionation formula [Taylor 1991]: where δD f is the final δD value (in % ) in the melt at the end of the process, δD i is the initial δD value of the melt at the beginning of the process, F is the fraction of H 2 O remaining dissolved in the melt, calculated as the final H 2 O at the end of the process divided by the initial H 2 O content in the melt at the start of the process , and α v-m is the bulk H-isotope fractionation factor between vapour and melt at specific temperature and pressure conditions. Here, in contrast to the closed-system case, the initial δD (δD i ) as well as initial bulk H 2 O content (H 2 O i ) are not constant because H 2 O vapour gets continually removed and consequently, the initial δD and H 2 O must be accordingly adjusted to lower values after each degassing step. In addition, degassing step size has an influence on the calculated δD value in the open system case because each successive "new" degassing aliquot must derive from an increasingly depleted hydrous reservoir. Therefore, calculations (in the open system) result in quite large isotopic fractionation (several tens of % ) at low H 2 O content. Thus, H 2 O-poor and strongly D-depleted lavas have been interpreted to be indicators of open-system degassing [e.g. Taylor et al. 1983].

Batched-system degassing
In batched-system or "multi-step open/closed" degassing [Taylor 1991] batches of exsolved gases are held resident in the magma for periods of time (e.g. in bubbles or cracks; Castro et al. [2012]), and then get episodically released from the system in outgassing pulses. These successive gas-release events could, in principle, correspond to real explosions or blasts observed during activity [Castro et al. 2014;Schipper et al. 2013].
To simulate a batched-system, both open-, and closed-system degassing elements are required; collectively these two essential steps constitute a hybrid degassing cycle. That is, the initial exsolution stage takes place in a finite closed system step (using the appropriate mass-balance Equation 2) with gas staying trapped in contact with melt. This closed-system is then dis-rupted during a subsequent opening of this system, which then removes that aliquot of exsolved gas. The hybrid degassing cycle then repeats-all the while simultaneously accounting for the updated and incrementally lower H 2 O budget and δD values-until the point that the system eventually becomes highly depleted in deuterium at low total H 2 O contents. As Castro et al. [2014] showed, this gas release mechanism best explains large H-isotopic fractionation they measured on a suite of pyroclastic and effusive products from the 2008 rhyolite eruption of Chaitén volcano [Castro and Dingwell 2009].
To simulate batched degassing, the formula of the closed system is used in a series of degassing increments, but the initial and δD i get continually updated at the beginning of each step. As will be shown, the degree of fractionation during batched process is heavily influenced by the step size. Furthermore, this system can be simulated with variable step sizes (e.g. by implementing a percentage rather than fixed amount of water loss at each step), to mimic the observation of variably sized explosions during hybrid activity [e.g. Schipper et al. 2013].

α v-m factor and speciation of H 2 O in rhyolitic melts
Regardless of the degassing system at hand (closed, open, or batched), variations of the alpha fractionation factor that occur during progressive H 2 O loss from the melt must be calculated [e.g. Newman et al. 1988]. The bulk H 2 O hydrogen fractionation factor-hereafter referred to as "alpha factor"-between vapour and melt (α v-m ) in rhyolitic melts, is defined as: where all variables above are defined as before.
The influence of the hydrous species ratio on the α v-m factor was first examined by Newman et al. [1988]. They determined that since H 2 O is present in the glasses [Stolper 1982] as both hydroxyl groups (OH´) and molecular H 2 O, the α v-m factor should also change when the proportions of those two species change with decreasing bulk H 2 O content. Similarly, Holloway and Blank [1994] postulated that as the relative proportions of dissolved hydrous species change when the total amount of H 2 O is declining, the proportion of OH´(which has a large D/H fractionation with water vapour) increases as degassing progresses. Consequently the bulk fractionation factor between the melt and vapour will increase as magma degasses. Moreover, the greatest change in in isotopic composition of the melt resulting from degassing will occur at low total dissolved H 2 O content. Experimental calibrations of the effects of temperature and total H 2 O content on hydrous speciation are available in empirical studies 3(1): 155 -168. doi: 1 .3 9 9/vol. 3. 1.155168 [Stolper 1982;1989].
In order to calculate the α v-m factor while accounting for its dependency on hydrous speciation, [Newman et al. 1988] used polynomial curve fits to plots of "total wt.% H 2 O" versus "wt.% as molecular H 2 O", which were in turn based on their own measurements of hydrous speciation in obsidians. Alternatively, another way of determining the hydrous species in the melt and α v-m , is the "Sil/Hol" model, which was originally developed by Stolper [1989] and later modified by Silver et al. [1990]. This "regular solution model" is implemented in a program by Holloway and Blank [1994] that numerically determines permissible values of speciation, output as "total wt.% H 2 O" versus " wt.% as molecular H 2 O" plots using empirical parameters.
To keep the variables needed to calculate degassing paths in VolcDeGas to a minimum, the following equations were used, derived from the ideal mixing model of Stolper [1989]. This model considers that the thermodynamic properties of the solution are analogous to those of a mixture of ideal gases. Firstly, the proportions of hydroxyl groups and the molecular H 2 O are governed by the equilibrium speciation equation with associated constant, K: where square brackets refer to the mole fraction of hydrous species in the melt on a single oxygen basis and O denotes anhydrous oxygen. The K factor is dependent on temperature and has been experimentally determined [e.g. Ihinger et al. 1999;Zhang 1999], whereas the temperature (T ) dependent K models of Zhang [1999] and Ihinger et al. [1999] produce nearly the same results. The model of Zhang [1999] has been utilised in this work: lnpKq " 1.89´ˆ3 120 T˙ (  6) where T is in kelvin.
VolcDeGas models the hydrous silicate melt as an ideal mixture of H 2 O molecules, hydroxyl groups and oxygen atoms, as proposed by Silver and Stolper [1985]. The following formula is used to calculate the mole fraction of OH´(XOH m ) [Silver and Stolper 1985]: where X bulk is the mole fraction of total H 2 O in the melt and can be calculated by [Silver and Stolper 1985]: where W is the molecular weight of anhydrous silicate melt on a single oxygen basis (" 30.04 g mol´1, based on the molecular weight of SiO 2 , 60.08 g mol´1), calculated by dividing the molecular weight of the substance by the number of oxygen atoms in the substance. Each hydrous species (OH´, H 2 O) possesses its own alpha fractionation factor and is calculated independently, for example, the alpha factor for molecular H 2 O (α v-m H 2 O ) is determined as such: Here δD represent values of molecular water in the melt phase. "Vapour" refers to the gas phase. The hydroxyl group fractionation factor also employs a similar formula.
Once the alpha factors of each hydrous species are determined, VolcDeGas calculates the bulk fractionation factor with the following equation [Rust et al. 2004;Taylor 1991]: where X is the atomic fraction of hydrogen incorporated as molecular H 2 O dissolved in the melt and (1-X) is the atomic fraction of hydrogen as OH´in the melt. However, if H 2 O speciation data is available (e.g. from Fourier-transform infrared spectroscopy-FTIRmeasurements), the following formula can be applied: where H 2 O m,f rac is the weight fraction of molecular H 2 O in relation to the total H 2 O content and OH g,f rac is the weight fraction of the hydroxyl groups. In this case, the speciation-specific alpha factors, α pv-mqH 2 O and α pv-mqOH are fixed values (0.9896 and 1.0415, respectively), after Rust et al. [2004].

Program functionality and verification
The program VolcDeGas was written in MatLab and uses the formulae described in the previous sections to calculate entire degassing histories (i.e. spanning a range of high-to low-H 2 O content). Degassing paths may be modelled as either one-stage or two-stage, and

VolcDeGas
Walter & Castro, 2020 may comprise open-, closed-, and batched system segments, as well as combinations thereof. The program is ideally suited for rhyolitic eruptions, whose glassy eruptive products have been routinely measured to yield δD versus bulk H 2 O content relations, therefore providing a natural degassing trend to simulate with VolcDeGas.
VolcDeGas is operated by way of a user-friendly GUI in which the required input parameters (e.g. temperature, initial δD, initial bulk H 2 O content, step size, degassing mechanism transition point, speciation data), and plot layouts can be easily set (Supplementary Material, Figures A1-A2). Each parameter's input location in the GUI is indicated by a series of clearly defined labels that appear by hovering the mouse over them, allowing the user to understand which functions the buttons/parameters carry out. Once these parameters are loaded, VolcDeGas automatically calculates and displays the best-fit degassing paths to natural, user-input δD-H 2 O data (e.g. based on lowest root-mean-square error (RMSE), and highest coefficient of determination (R 2 ); Supplementary Material). Owing to the underlying Matlab code, which is thoroughly described in the dialog window in the program interface, the user can perform literally thousands of calculations per minute, a feat not possible by employing standard spreadsheet and graphing programs [e.g. Castro et al. 2014]. An in-depth description of the code can be found in the code comments and the program manual (Supplementary Information).
In order to test the modelling capabilities of VolcDeGas, we used published bulk H 2 O and H-isotopic data from the 1340 C.E. Mono Craters, CA eruption as a basis for model simulations [Dobson et al. 1989;Newman et al. 1988]. In addition, we also applied VolcDeGas to data from the 2008 Chaitén Eruption; those tests are discussed at the end of this section.
Our approach with respect to the Mono Crater's eruption was to compare our model fits to natural data with the degassing paths fit by Dobson et al. [1989].  Figure 3C).  Figure 2C is reproduced from DeGroat [1997] to demonstrate the differences between alpha factor values computed with a modified version of the "Sil/Hol" Program [Holloway and Blank 1994;Silver et al. 1990] and using natural speciation data of Newman et al. [1988]. According to Figure 2, the primary difference between alpha factor simulations ( Figure 2A) and literature-derived data fits (Figure 2B-C) is the curved versus linear projections of alpha factor. Despite this difference, the alpha values are relatively similar, even though the offsets grow larger as the H 2 O increases. As we show below, despite some alpha factor inequalities, models will result in just a few per mil difference in the final calculations of δD.

Modelling pyroclastic H 2 O-δD patterns in the
Mono Craters 1340 C.E. eruption Figure 3A and 3B show the results of VolcDeGas calculations of H 2 O-δD trends for the Mono Craters system comprising an inferred closed-to open-system degassing history [Dobson et al. 1989] and using the previously determined alpha factors (e.g. Figure 2A, B). Here, the same parameters (e.g. starting H 2 O and δD values) as previously described are applied, with the exception of temperature, which is set to 800°C. The results demonstrate that the computations of the degassing history that utilise natural speciation data are the best match to the natural data ( Figure 3B) and those analogous calculations from literature ( Figure 3 [ Dobson et al. 1989]). A lower quality fit to data arises from using the purely simulated hydrous speciation (Figure 3A), illustrating that it is preferable to use available natural speciation data for the simulations. Such variations of δD between Figure 3A and Figure 3B, C only amount to 3.3 % in the closed-system, and up to 8 % in the open-system case (Figure 3). In the next section, we demonstrate how VolcDeGas can be used to calculate permissible degassing paths and resolve best-fit magmatic degassing parameters in another natural rhyolitic system, Chaitén volcano (Chile).

Modelling pyroclastic H 2 O-δD patterns in the Chaitén 2008 eruption
In this section, we demonstrate a case study involving a rhyolitic volcanic eruption that was directly observed and therefore provides essential information about the style of explosive and effusive activity that was responsible for the samples that have been measured for their  [Barnes et al. 2014], respectively. In [C] the linear fits to speciation values of Newman et al. [1988] and of the Silver/Holloway Program [Holloway and Blank 1994;Silver et al. 1990] are compared.  [Dobson et al. 1989;Newman et al. 1988].
[A] and [B] show the computed values by VolcDeGas. The blue line displays the VolcDeGas simulated closed-, to open-system degassing graph, the red dots the sample data and the black cross shows the transition where the degassing system changes.
[C] presents a published model fits to data from the literature [Dobson et al. 1989], with red dots being the sample data [Newman et al. 1988]. Numbers on two plots in [C] indicate the alpha factor values used by Dobson et al. [1989].  Figure 2A), and in [B] on measured speciation values from Barnes et al. [2014] (see also Figure 2B).

VolcDeGas
Walter &Castro, 2020 H-isotopic andH 2 O contents [e.g. Castro et al. 2014]. Our main goal is to illustrate how, using the Chaitén dataset, a user would be able to produce a best-fit degassing trend to measured data, which in turn serves to demonstrate the effectiveness and versatility of the program. It is not the goal however to interpret the volcanological underpinnings of the dataset, as was done previously by Castro et al. [2014]. Degassing simulations were performed using FTIRbased bulk H 2 O and hydrous speciation data of pyroclastic obsidians from the 2008 Chaitén eruption [Castro et al. 2014;2012;Lara 2009, Figure 4A]. We first used VolcDeGas to test for effects of external hydration (e.g. meteoric water sourced), then computed alpha factors based on species-specific alpha factors (α pv-mqH 2 O and α pv-mqOH ) from Rust et al. [2004] and the non-hydrated speciation data. VolcDeGas makes speciation extrapolations via polynomial fitting (Figure 4).
Testing for hydration involved first generating hydrous speciation graphs with the program, which are polynomial curve fits to natural speciation and bulk H 2 O content data (Figure 4). The fitted natural speciation data, which are derived from Castro et al. [2012], allowed for extrapolation of OH/H 2 O values, the trends of which comprise a clear envelope corresponding to a high-temperature "magmatic" array ( Figure 4A). Hydration is recognised in samples that contain too much molecular H 2 O with respect to the curves, which is a hallmark of late-stage alteration of hydrothermal fluids ( Figure 4A). Data points that fall off the speciation fits (by~0.3 wt.%) are suspect and filtered out as they are not suitable for the best-fit alpha factor calculations. This is important because altered speciation values may lead to alpha factor values under 1. In such cases, VolcDeGas will return an error, as the degassing trend would lead to more highly enriched δD with falling bulk H 2 O content, instead of the depleted δD values typical of magmatic degassing. Figure 5A and 5B show degassing paths calculated for variable step sizes in open-and batched-degassing systems. As can be seen in Figure 5A, although the step size in the open-system is changed substantially, (from 0.5 wt.% down to 0.001 wt.% steps) the influence on the δD value is relatively negligible across all degassing intervals. These step size changes result in δD differences of only 3 % at 0.5 wt.% bulk H 2 O and of 4.9 % at 0.001 wt.% H 2 O. Furthermore, the open-system plot resolution for bigger step sizes is lower as less points are calculated. It consequently appears as if there is a strong fractionation in δD over the last 0.5 wt.% H 2 O. In reality, this is only an optical artifact as the space between the steps are connected (interpolated) by a line which has no calculated intermediate values.
In contrast to the open system, the batched-degassing system ( Figure 5B) is highly sensitive to step size variations. In particular, δD differences at variable step sizes sum up to over 7 % at 0.5 wt.% bulk H 2 O content and reach~170 % at 0.001 wt.% H 2 O. These ex-treme fractionation differences occurring at the lowest of H 2 O contents (0.001 wt.%) are meant to demonstrate the theoretical maximum disparity in the step-size effect. In natural magmatic systems, melts rarely degas to values lower than about 0.05 wt.%, however, even at these values small modelled steps may still result in large δD offsets (~70-100 % ; Figure 5B). This observation highlights the importance of step size during the end-stages of degassing. Figure 5B as well as other calculations included in the supplementary information support the idea that simulated batched magma degassing can mimic either open-, or predominantly closed-system behaviour by mere selection of the appropriate step sizes. A large step size contributes to a closed-system trend (i.e. less overall δD depletion), whereas a Thus, VolcDeGas runs of the batched system, in which a percentage of H 2 O is lost at each step, may best approximate the natural ssystem [e.g. Castro et al. 2014]. To better illustrate this point, Figure 6 shows a batched system calculation with an initial δD value of´61 % , and a variable step size, implemented by setting an H 2 O loss of 40 % at each step (see Supplementary Material for additional conditions). Goodness-of-fit for this simulation is reasonable, as evidenced by the adjusted coefficient of determination (R 2 adj ) and R 2 , as well as the RMSE (Figure 6).
In Figure 7, we investigate the possibility of a twostage history for the Chaitén eruption, one involving early closed-system degassing followed by batch degassing (both at a presumed eruption T of 800˝C; [Castro and Dingwell 2009]). This scenario may for example mimic degassing across an entire eruptive sequence, from high-energy, pulsating Plinian activity (e.g. purely closed-system dynamics; [e.g. Newman et al. 1988]) to later hybrid activity involving both rhyolite pyroclast and lava eruption (e.g. batcheddegassing; [Castro et al. 2014;). Of the two-stage best-fit models ( Figure 7A, B), the closed-, to batchedsystem with a small, fixed step size (0.001 wt.% H 2 O) exhibits slightly better results ( Figure 7A), as reflected by calculated goodness-of-fit (R 2 and RMSE). However, both batched-system models exhibit similar fit statistics, suggesting that either is a good choice to model the Chaitén sequence ( Figure 7B).
Another factor that affects goodness-of-fit and therefore must be taken in to account in order to model a two-stage system, is the transition point, defined as the H 2 O content at which the closed-system gives way to the batched-system. Transitions can also be generally defined as changeovers between any other combination of degassing styles (e.g. closed-to open-system degassing). In the previously modelled closed-to-batched system with the constant step size ( Figure 7A), the bestfit transition point is at 0.5 wt.% H 2 O, whilst in the system with variable step size ( Figure 7B), the transition occurred at 0.6 wt.% H 2 O. These transition points may seem very similar but the differences are significant and exemplify the sensitivity of VolcDeGas' fitting process, which entails an iterative determination of the "most likely" transition point through statistical minimisation of RMSE values (Supplementary Material). This objective way of determining a transition point could offer an advantage over previous manual transition point placement methods [e.g. Dobson et al. 1989].
In summary, two-stage degassing systems (Figure 7) fit the sample data significantly better than the previously analysed one-stage model ( Figure 6). In particular, the relevant RMSE-value is around 0.5 lower, whereas the R 2 is~0.3 higher. As a result, twostage closed-, to batched-system degassing paths can be viewed as the best-fit for the available data points. In the following section, batched systems with constantand with variable-step sizes are discussed and compared based on available literature data as well as observations from the eruptions themselves, to identify the most probable system of the two.

Discussion
The VolcDeGas program is a useful tool for investigating feasible degassing histories of rhyolitic magmatic systems, provided the appropriate D/H isotopic and bulk H 2 O content data are available. Because the program integrates all of the required parameters in degassing calculations-e.g. real and extrapolated hydrous speciation, alpha factors, temperature, and degassing step size (both fixed and variable)-the user can investigate the influence of these variables on one-, and two-stage degassing systems in an objective manner. The program is furthermore not a "black box" insofar as it provides an intuitive GUI for entering the data, setting key parameters and visualising output in a variety of multiaxis plots. A view of the code itself (see VolcDeGas.m) shows a detailed description of how the program functions. As a result, after a short familiarisation with this program, the user can quickly calculate the different kinds of degassing systems, while being able to comprehend how the program executes these calculations.
Of the many initial conditions that influence isotopic fractionation calculations, the partitioning of deuterium and hydrogen between vapour and melt, encapsulated in the alpha factor, imparts a strong influence on degassing trends. This is because the alpha factor is a function of hydrous speciation in the melt, which in turn varies as degassing progresses. That VolcDeGas is able to calculate and update the alpha factor according to real or modelled speciation values during a simulation is a major advantage that requires relatively little user intervention, however, the user can also use this functionality to spot potential sources of error in the simulation that stem from natural processes (e.g. hydration). For example, after loading hydrous speciation data, the user can inspect the data and associated VolcDeGas speciation curve fits identify anomalous values (Figure 4). Altered samples (e.g. those with high molecular H 2 O) will change the speciation graph extrapolations significantly, and in turn calculated alpha factors may be less than 1. VolcDeGas will flag these instances and output an error message and stop the simulation, thereby precluding counterintuitive, positively sloping degassing trends. Finally, as some degassing paths yielding very good fits to natural data were calculated utilising alpha factors ( Figure 3B) that are based on natural hydrous speciation values, we recommend wherever possible, to use available natural speciation-based calculation of alpha factor, instead of the temperature-dependent ideal mixing model species calculations.
Test models show that in addition to alpha factor, the step size-the amount of H 2 O lost during Figure 6: Modeled, best-fit, single staged batchedsystem trend with variable step size for the 2008 Chaitén eruption products. Simulations performed at an assumed eruption temperature of 800˝C [Castro and Dingwell 2009]. Data from Castro et al. [2014] are shown as red dots. Here each step invokes initial closed-and later open-system release of 40 % of the H 2 O budget. The fit is particularly good at low H 2 O contents, in materials that have undergone the most degassing.
each interval-is also hugely important in the form of degassing trends and yielding good fits to natural data (e.g. Figure 5). Batched-degassing simulations exhibit particularly strong sensitivity to this variable ( Figure 5B). In the most extreme case, δD differences between moderate (0.5 wt.%) versus small steps (0.001 wt.%) amount to more than 100 % variation in the small window of degassing from about 0.5 wt.% to 0.1 wt.% (e.g. Figure 5B). This clearly shows that in the low-H 2 O content interval, δD in a magma responds sensitively to each progressive degassing step, with the implication that extreme deuterium-depletion in lavas could portend many repetitive small steps rather than large pulses in the waning stages of an eruption. Indeed, VolcDeGas models of the Chaitén pyroclastic and lava data set support this conclusion (Figures 6 and 7B), and exemplify how progressively decreasing step sizes may provide an excellent fit to natural products. The Chaitén simulation, shown in Figure 6, employs a degassing step that removes 40 % of the remaining H 2 O at each interval, which is tantamount to having variable, diminishing step sizes. The result is a relatively smooth trend and good fit for a one-stage system to natural data [see also e.g. Castro et al. 2014], which provides good evidence for a repetitive degassing process having been operant in the natural eruption.
In addition to single-stage (e.g. closed, open, or batched) degassing trajectories, VolcDeGas is capable of calculating two-step degassing paths with good precision, and with a relatively small number of parame- Step sizes differ from [A] a small fixed step size 0.001 wt.%, to [B] a variable one comprising 15 % H 2 O loss at each step. Despite these differences, model fits to data are virtually identical. In addition, the goodness-of-fit of these two-stage systems is better than observed in the singlestage batched system with variable step size ( Figure 6). All simulations were performed at a presumed eruption temperature of 800˝C [Castro and Dingwell 2009]. ters required as initial inputs. Especially relevant is the capacity for VolcDeGas to find the best "transition point" whereby the degassing mechanism changes from one mode to another (e.g. purely closed to open). Traditional closed-to-open, two-stage degassing interpretations of natural data [e.g. Newman et al. 1988;Rust et al. 2004] are based in part on evidence in the H 2 O contents of melt inclusions, which may contain several wt.% H 2 O, and pyroclastic obsidians, which can be significantly degassed (down to a few tenths of a wt.% H 2 O), together suggesting that initial Plinian activity must be fueled by vigorous closed-system degassing involving fragmentation as the main mech-anism to liberate volatiles. Eventually a "transition point" is reached, marking an abrupt change in the degassing mode (i.e. to an open system). The modelled trends should therefore be "kinked". This is indeed well depicted by VolcDeGas' closed-to-open system degassing simulation of the D/H-H 2 O data for the Mono Craters (Figure 3). These calculations yield plots that are very good matches to the δD and H 2 O values in Mono Craters' pyroclasts and lavas, and even reproduce the transition points and fits to these data from the published literature [Dobson et al. 1989;Newman et al. 1988]. This good correspondence verifies this program's ability to reproduce and confirm previously published model trends.
Two-step degassing simulations may also comprise closed-to batched-system degassing sequences, as was done for example for Chaitén D/H-H 2 O data (Figure 7). In this simulation ( Figure 7) the initial closedsystem degassing trend can be representative of the first explosive eruption(s) that must have liberated large quantities of H 2 O. Following this, the progression gives way to a batched system, demarcated by a kink in the graph and corresponding to the best-fit transition point (Figure 7). The batched system with variable step size ( Figure 7B) is able to explain the degassing data after the closed-system phase and could represent hybrid explosive-effusive activity characterised by declining explosive energy with time, which in turn releases fewer volatiles from each successive batched degassing step.

Summary and conclusions
In summary, the VolcDeGas program is effective at modelling degassing paths as well as being a tool for fitting and interpreting natural data that can lead to new insights into magmatic processes. Our understanding of rhyolite eruptions can in turn be more quantitatively based, ultimately fostering a better connection between observations of activity and measured geochemical degassing proxies (e.g. H 2 O and D/H). Establishing how step size variations modelled in VolcDeGas may scale to actual volcanic eruption processes-in particular the repetitive blast events that characterise hybrid activity [Castro et al. 2014;Saubin et al. 2016]-is one goal that may improve coupling natural field observations to analysis of eruption products. For example, it is currently unknown how big the individual blasts are and to what extent each blast "bleeds" off vapour pressure from the underlying magma body. Therefore, future field and analytical studies could profitably focus on obtaining estimates of degassing amounts during repetitive volcanic blasts. Linking these data with VolcDeGas could pave the way to understanding how degassing proceeds in a temporal context and deliver the root causes of different eruption styles.