3D Diffusion of Water in Melt Inclusion‐Bearing Olivine Phenocrysts

Olivine‐hosted melt inclusions are an important archive of pre‐eruptive processes such as magma storage, mixing and subsequent ascent through the crust. However, this record can be modified by post‐entrapment diffusion of H+ through the olivine lattice. Existing studies often use spherical or 1D models to track melt inclusion dehydration that fail to account for complexities in geometry, diffusive anisotropy and sectioning effects. Here we develop a finite element 3D multiphase diffusion model for the dehydration of olivine‐hosted melt inclusions that includes natural crystal geometries and multiple melt inclusions. We use our 3D model to test the reliability of simplified analytical and numerical models (1D and 2D) using magma ascent conditions from the 1977 eruption of Seguam volcano, Alaska. We find that 1D models underestimate melt inclusion water loss, typically by ∼30%, and thus underestimate magma decompression rates, by up to a factor of 5, when compared to the 3D models. An anisotropic analytical solution that we present performs well and recovers decompression rates within a factor of 2, in the situations in which it is valid. 3D models that include multiple melt inclusions show that inclusions can shield each other and reduce the amount of water loss upon ascent. This shielding effect depends on decompression rate, melt inclusion size, and crystallographic direction. Our modeling approach shows that factors such as 3D crystal geometry and melt inclusion configuration can play an important role in constraining accurate decompression rates and recovering water contents in natural magmatic systems.


Introduction
Water plays a significant role in the formation and physico-chemical evolution of magma.The addition of water can induce melting in the sub-arc mantle wedge by suppressing the solidus (Grove et al., 2006;Till et al., 2012), and is ultimately responsible for volcanism in subduction zone settings (Grove et al., 2012).Water can act as a network modifier in silicate melts meaning it can exert an important control on the rheological properties of magma such as viscosity and yield strength (Dingwell, 1996;Gonnermann & Manga, 2007;Russell et al., 2002), which in turn may exert first order controls on the storage depths of arc magmas (Rasmussen et al., 2022).The exsolution of water and other volatiles into a vapor phase can also affect magma compressibility (McCormick-Kilbride et al., 2016), eruptibility (Stock et al., 2016), and ultimately eruptive style (Cassidy et al., 2016).Furthermore, the addition of water to crystallizing magma can have a large influence on mineral stability and subsequently chemical differentiation (Gaetani et al., 1993); a notable example being the differences in the tholeiitic and calc-alkaline trends (Zimmer et al., 2010).
Measuring reliable water contents in magmas is therefore of great importance for understanding volcanic systems, but has proven challenging.As magma ascends to the surface, the solubility of water in the melt decreases meaning it exsolves into a vapor phase once kinetic and physical barriers to bubble nucleation are overcome (Gonnermann & Manga, 2007).Degassing can further drive the buoyant ascent of magma and lead to fragmentation and explosive eruptions; however, it means initial water contents cannot be recovered from the erupted matrix glass.Petrologists have subsequently turned to olivine-hosted melt inclusions in order to measure preeruptive volatile contents in basaltic magmas to understand primary mantle melt compositions (Sobolev & Chaussidon, 1996) and storage conditions at the time of entrapment (Ruth et al., 2018;Wieser et al., 2021).This approach has been facilitated by advances in microanalytical techniques such as secondary ion mass spectrometry (SIMS) and Fourier transform infrared spectroscopy (FTIR), but assumes that post-entrapment modification such as crystallisation (Steele-Macinnis et al., 2011), decrepitation (Maclennan, 2017), or diffusive loss through the crystal host (Bucholz et al., 2013;Gaetani et al., 2012) are minimal or can be corrected for.The water content of olivine-hosted melt inclusions can be dominantly modified by diffusive loss or gain of hydrogen (as H + ) through the olivine lattice (Barth & Plank, 2021;Barth et al., 2019;Gaetani et al., 2012;Hartley et al., 2015).The ability of melt inclusions to retain their water content depends on a combination of geometrical factors (e.g., crystal size, melt inclusion size and melt inclusion position) and physico-chemical parameters (e.g., magma decompression rate, temperature, degassing style, diffusivity, and olivine-melt partitioning behavior).The diffusion of H + , which can involve multiple mechanisms (Barth et al., 2019;Ferriss et al., 2018) and inter-site reaction (Jollands et al., 2019), is typically very rapid meaning most melt inclusions rarely preserve their initial water contents even for ascent timescales on the order of hours (Barth & Plank, 2021;Gaetani et al., 2012).Consequently, initial water contents typically have to be reconstructed using relationships between water and incompatible major elements (e.g., K) or trace elements (e.g., Ce) across the melt inclusion population (Barth et al., 2019;Hartley et al., 2015).
While rapid diffusive loss of water through the olivine lattice hinders initial water content estimation, it does offer considerable promise as a chronometer to track decompression rates during final magma ascent (Barth et al., 2019;Le Voyer et al., 2014;Newcombe, Plank, Barth, et al., 2020;Peslier & Luhr, 2006).Popular olivine geospeedometers (e.g., Fe-Mg interdiffusion, Ni, Ca) used to time pre-eruptive magma mixing and ascent (Couperthwaite et al., 2020;Kahl et al., 2023;Mutch, Maclennan, Shorttle, et al., 2019;Rae et al., 2016;Rasmussen et al., 2018) typically respond to chemical changes induced at depth and do not have the temporal resolution with current microanalytical techniques to track processes over timescales of minutes to hours (Bradshaw & Kent, 2017).Fluid-mobile monovalent cations, such as hydrogen and lithium, can produce resolvable diffusion profiles over minutes to hours at magmatic temperatures (Lynn et al., 2018;Newcombe, Plank, Barth, et al., 2020).Their composition in the melt responds directly to the degassing process meaning diffusion gradients can provide high resolution temporal information about final magma ascent (Lynn et al., 2018;Neukampf et al., 2021).These chronometers typically provide an average estimate of magma decompression from vapor saturation, but can be combined with other complementary methods such as volatile diffusion in melt embayments (deGraffenried & Shea, 2021;Ferguson et al., 2016;Moussallam et al., 2019;Myers et al., 2019), crystal size distributions (Cashman & Blundy, 2000;Gurioli et al., 2005), mineral growth or dissolution (Neave & Maclennan, 2020), or bubble number densities (Hajimirza et al., 2021;Myers et al., 2021;Shea et al., 2010) to obtain a more complete picture of a magma's ascent history.Extracting this kind of information about magma ascent opens up the exciting possibility to link magma decompression rates to physico-chemical parameters and the predictions of conduit models (Barth et al., 2019).In fact, correlations between magma decompression rate and eruption style have recently been established, with basaltic eruptions with higher mass eruption rates being associated with faster decompression rates (Barth et al., 2019).
Diffusive loss of hydrogen through olivine has subsequently become a popular chronometer and has been used in two main ways: measuring and modeling water loss from a population of olivine hosted melt inclusions (Barth et al., 2019) or measuring water gradients directly in olivine (Newcombe, Plank, Barth, et al., 2020).The former offers the advantage that the water content of melt inclusions can be precisely measured by SIMS, whilst the latter can be used in olivine crystals that do not contain melt inclusions.Both of these modeling approaches have used assumptions either about crystal geometry or diffusion in multiple phases.3D spherical approximations that assume infinite diffusivity in the melt inclusion and isotropic diffusion in the olivine host have been commonly used (Chen et al., 2011(Chen et al., , 2013;;Cottrell et al., 2002;Gaetani et al., 2012).Diffusion of water in 1D versus 3D in nominally anhydrous minerals has been tested by Thoraval and Demouchy (2014).They compared 1D plane models to 3D spherical and prismatic morphologies for different H + diffusion mechanisms in olivine.They concluded that 1D models can be accurate if the diffusive anisotropy is above an order of magnitude.1D approximations along the [100] axis subsequently became favored due to the high anisotropy measured for H + in Febearing olivines with diffusion along the [100] potentially being more than 10 times faster than the [010] or [001] axes (Barth et al., 2019).The melt inclusion approach of Barth et al. (2019) assumes that diffusive loss dominantly occurs along the [100] axis, diffusivity in the melt inclusion is infinite and constant, and that the crystal is symmetric about the center of the melt inclusion (reflective boundary condition in the center of the melt inclusion).Likewise, Newcombe, Plank, Barth, et al. (2020) apply a 1D approximation along [100], but do not formally incorporate any complexities associated with melt inclusions.In addition, both modeling approaches fail to account for the interaction of multiple melt inclusions in a single crystal and how that could influence water loss.
In this study, we present a new multiphase finite element model for the concurrent diffusion of water in the melt inclusion and the surrounding host olivine.This approach facilitates the use of complex 3D geometries associated with idealized olivine crystal morphologies, and can simulate diffusive loss from multiple melt inclusions akin to crystals observed in natural volcanic systems.We first present analyses of the water content of melt inclusionbearing phenocrysts from the 1977 eruption of Seguam volcano, Aleutians.The well-characterized decompression history of this eruption based on multiple petrological methods makes it an ideal candidate to test our modeling approach.We then perform 2D model inversions on these phenocrysts to guide our selection of input parameters for our 3D models.Finally, we use our 3D models to explore uncertainties with assumptions based on geometry and dimensionality.We test uncertainties associated with 1D and 2D approximations and analytical solutions during magma ascent with a degassing boundary condition, and we test how multiple melt inclusions interact during water loss.We then discuss modeling caveats and provide recommendations for future modeling practice.

1977 Eruption of Seguam, Aleutian Arc
Application of geospeedometers based on the dehydration of melt inclusion and olivine macrocryst populations has the greatest potential in water-rich arc magmas, where the water contents of melt inclusions and gradients in the host crystals can be easily resolved using current microanalytical methods.In this study, we shall focus on erupted products from the March 1977 eruption of Seguam volcano to provide a framework for modeling water loss in olivine crystals from arc systems.
Seguam Island is an ∼200 km 3 volcanic complex in the central Aleutian Island arc, with a measured eruption chronology dating back to 318 k.y. in the Pleistocene (Jicha & Singer, 2006).The volcano has been historically active with eight eruptions over the past 200 years: notably two basaltic and basaltic-andesite fissure eruptions occurring in 1977and 1992-1993(Jicha & Singer, 2006;Miller et al., 1998).The 1977 eruption took place between the 6th and 8th of March, and was associated with 8 lava fountains erupting from a ∼1 km long radial rift on Pyre Peak (Miller et al., 1998).The lava fountains reached heights of up to 90 m, culminating in a maximum volcanic explosivity index estimate of 1 for the eruption.The estimated volume of erupted basaltic material, including lava flows and pyroclastic deposits, is ∼0.06 km 3 based on stratigraphic analysis of present day deposits (Jicha & Singer, 2006).The crystal cargo contains olivine, clinopyroxene and plagioclase (Lloyd et al., 2016;Zimmer et al., 2010).The olivine crystals are melt inclusion rich (Figure 1).Attempts to characterize the pre-eruptive water content of the 1977 Seguam magma were made by Zimmer et al. (2010), who found that melt inclusion suites in forsteritic olivines (Fo 80-85 ) from basaltic tephra have been modified the least by degassing and preserve water contents up to 3.7 wt%.Plagioclase and clinopyroxene melt inclusions from more evolved samples have much lower water contents (<2.5 wt%), which have been attributed to degassing (Zimmer et al., 2010).Comparison of the water content of clinopyroxene-hosted melt inclusions in lava flow and tephra samples show diffusive loss of water in the lava flow in less than 1 hr (Lloyd et al., 2016).Syneruptive decompression rates have recently been estimated by Newcombe, Plank, Barth, et al. (2020) using a combination of volatile diffusion in olivine-hosted melt embayments and modeling of water concentration gradients across melt inclusion free olivine phenocrysts.Newcombe, Plank, Barth, et al. (2020) found that decompression rate estimates for the 1977 Seguam eruption from both of these methods were fairly consistent and ranged from 0.02 to 0.23 MPa s 1 .

Measurement of Water Concentration Gradients in Olivine Phenocrysts by SIMS
Concentrations of volatiles (H 2 O, CO 2 , S, Cl, F) and P were characterized using the Cameca IMS 6f ion microprobe at Carnegie Earth and Planets Laboratory, following previously developed analytical protocols (Hauri, 2002;Koga et al., 2003;Kumamoto et al., 2017;Mosenfelder et al., 2011;Newcombe, Plank, Barth, et al., 2020).The analytical protocol is described in detail by Newcombe, Plank, Barth, et al. (2020) and is summarized here.Samples were gold coated and placed into the sample exchange chamber of the SIMS one to 3 days prior to the beginning of the analytical session.A Cs + primary beam with a current of ∼15-20 nA was tuned to achieve an approximate beam diameter of ∼20 μm, and charge compensation was provided by an electron flood gun.The sample was pre-sputtered for 120 s with the rastered primary beam in order to remove surface contamination, and a field aperture was used to ensure the collection of ions from the central ∼10 μm of the measurement area.Negatively charged 12 C , 16 O 1 H , 19 F , 30 Si , 31 P , 32 S , and 35 Cl ions were detected using an electron multiplier.Five cycles of data were collected at each point.Electrical glitches occasionally produced "empty" cycles in which no counts were recorded; these cycles were deleted prior to data reduction.The mass resolving power was sufficient to resolve 16 O 1 H from 17 O.Data reduction was approached following a protocol similar to that described by Kumamoto et al. (2017).Basaltic glasses and grains of olivine and orthopyroxene with well-characterized volatile contents were used as standards (see Data Sets S1 and S2).The calibrations of water concentration in olivine used in this study are based on absolute water concentrations of olivine and orthopyroxene standards determined by Bell et al. (2003).Counts of volatile and 31 P ions were divided by counts of 30 Si in order to account for instrumental drift and multiplied by SiO 2 /50 (where SiO 2 is the independently measured silica concentration of the sample in wt%).Drift in the H 2 O analyses was further monitored and corrected for using frequent analyses of "Herasil 102" silica glass and basaltic glass ALV-519-4-1.The detection limit (Long & Winefordner, 1983) was calculated using replicate analyses of "Suprasil 3002" silica glass, and was found to be <4 ppm H 2 O (this estimate is an upper bound, given the nonzero water concentration of Suprasil glass and the fact that the quality of the vacuum improved over the course of each session).

Measurement of Major Element Compositions by Electron Microprobe
Olivines were analyzed for major, minor, and trace elements using a Cameca SX100 electron microprobe at the American Museum of Natural History (AMNH) and a JEOL JXA 8900R electron microprobe at the University of Maryland (UMD).The AMNH data were previously published by Newcombe, Plank, Barth, et al. (2020), whilst the UMD profiles are shown the in Supplementary Data.American Museum of Natural History analyses were conducted using a 15 kV accelerating potential, a 10 nA beam current and a focused beam (with nominal diameter ∼1 μm).UMD analyses were conducted using a 20 kV accelerating potential, a 20 nA beam current and a focused 1 μm beam.On-peak counting times at AMNH varied between 20s (Mg, Si, Ca, Mn, Al, Fe), and 40s (Ti and P).At UMD, on-peak counting times varied between 20s (Mg, Si, Fe) and 40s (Ca, Ni). Background counting times were set to 50% of the on-peak counting times.Analyses were corrected for inter-run calibration offsets using factors determined by replicate analyses of San Carlos olivine.Replicate analyses of San Carlos olivine yielded average 2 RSDs of <3% for SiO 2 , MgO, and FeO.

Numerical Modeling of Water Diffusion in Melt Inclusion-Bearing Olivine Phenocrysts During Syn-Eruptive Magma Decompression
In this section we first review important parameters required by our models including: the olivine-melt partition coefficient and diffusion coefficients of H species in the olivine and the melt.We then go on to describe our 2D inversion and 3D modeling approaches.

Approach to Modeling the Olivine-Melt Partition Coefficient
Hydrogen dominantly dissolves as OH groups in olivine, however it can dissolve as both molecular H 2 O and OH species in the melt (Demouchy & Mackwell, 2006;Stolper, 1982).This speciation can subsequently affect partitioning behavior and be associated with different partition coefficients.For simplicity, and to maintain consistency with previous studies (Barth et al., 2019;Gaetani et al., 2012;Newcombe, Plank, Barth, et al., 2020), we have used a constant partition coefficient based on the H 2 O equivalents in the olivine and melt: where C H 2 O Ol and C H 2 O MI are the bulk H 2 O equivalent concentrations in the olivine and melt inclusion respectively.To avoid confusion with the diffusion coefficient, D, we opted to use K to represent the partition coefficient, however we appreciate that this has traditionally been used to represent equilibrium constants and exchange coefficients.The olivine-melt partition coefficient for bulk water can range from 0.0005 (Le Voyer et al., 2014;Newcombe, Plank, Barth, et al., 2020) to 0.003 (Hauri et al., 2006;Koga et al., 2003;Towbin et al., 2023).

Approach to Modeling Diffusion in the Host Olivine
There have been many experimental studies that have tried to characterize the diffusion mechanisms of H + through the olivine lattice during either hydration or dehydration (Barth et al., 2019(Barth et al., , 2023;;Jollands et al., 2019;Kohlstedt & Mackwell, 1998;Mackwell & Kohlstedt, 1990;Padrón-Navarta et al., 2014;Peslier et al., 2015).Multiple diffusion mechanisms have been observed, which explains the six orders of magnitude variation in experimentally derived diffusion coefficients at magmatic temperatures (∼1,000 °C).The fastest mechanism measured in pure forsterite, called proton-polaron exchange, is associated with a flux of H + being chargebalanced by a flux of electrons from Fe 2+ to Fe 3+ (Kohlstedt & Mackwell, 1998;Mackwell & Kohlstedt, 1990).This redox process progressively reduces Fe 3+ .A slower mechanism associated with M-site vacancies has also been measured where a vacant M-site charge balanced by 2H + exchanges with either Fe 2+ or Mg 2+ parallel to the main H + gradient.The slowest mechanism is associated with T-site vacancies, where a tetrahedral vacancy charge balanced by 4H + exchanges with Si 4+ (Jollands et al., 2019;Padrón-Navarta et al., 2014).Recently, it has been shown that H + may migrate via a coupled reaction-diffusion process which could involve exchange between different sites and subsequent diffusion via the different mechanisms (Jollands et al., 2019).Despite this complexity, dehydration experiments on natural olivines can describe bulk loss of H + with simple Arrhenius relationships (Barth et al., 2019;Ferriss et al., 2018;Wallace et al., 2021).Dehydration experiments conducted by Barth et al. (2019Barth et al. ( , 2023) ) on Fe-bearing olivines from Cerro Negro (Fo 79-81 ) and Etna (Fo 90 ) show significantly higher diffusivities for bulk H + than those observed in pure synthetic forsterite.They found that diffusion was considerably faster along the [100] direction, which is consistent with the proton-polaron mechanism.Barth et al. (2019) attribute their higher diffusivities to the higher Fe content of their olivines.
Given the wide range in recorded diffusivities, it is therefore important to select the appropriate diffusion mechanism and diffusion coefficient for the system of interest.In this study, we are interested in modeling dehydration in Fe-bearing olivines analogous to the crystal cargo observed at Seguam (Fo 80-85 ).We have opted to use the diffusion coefficient of Barth et al. (2019), which was derived from experimental dehydration of ∼Fo 80 olivines: This diffusion coefficient was calibrated using 2 experiments run at 800 and 1,000 °C, and subsequently requires a small extrapolation up to temperatures appropriate for Seguam magma (∼1,070 °C).Barth et al. (2019) show that extrapolation of their diffusion coefficient to 1,100 ± 12 °C can lead to calculated diffusivities that range from 10 9.87 to 10 9.67 m 2 s 1 .This would represent an upper limit of the extrapolated range.To account for the covariance in uncertainty structure of the parameters contributing to the diffusion coefficient (Mutch et al., 2021), we generated a covariance matrix using the error propagation of Barth and Plank (2021).Given the lack of experimental data available for dehydration in Fe-bearing olivines, this covariance is just a basic estimate, however it implicitly accounts for the uncertainties associated with extrapolation.

Approach to Modeling Diffusion in Melt Inclusions
Diffusivity of water in the melt inclusion was calculated using the diffusion coefficient of Ni and Zhang (2018).
Here the diffusivity of total water is dependent on temperature, pressure, and melt composition.Most importantly, there is a dependence of diffusivity on water content, which adds extra complexity that affects the efficiency of the computational model.This is because the diffusion equation now takes a non-linear form and has to be solved using iterative solvers at each time step.The expression of Ni and Zhang (2018) involves a lot of input parameters, which can also slow down the models.For 3D models with a high number of mesh points we used an empirical polynomial expression derived from solutions to the Ni and Zhang (2018) This parameterization makes computation of D MI much faster given the smaller number of variables, however it can only account for the effect of water content and is only applicable to melt compositions similar to Seguam melt inclusion Seg1-MI1 (Newcombe, Plank, Zhang, et al., 2020) at a temperature of 1,070 °C.In this study we assume isothermal ascent (Newcombe, Plank, Zhang, et al., 2020) and that the major element composition of the melt inclusions do not change during ascent.Seguam melt inclusions suffered very little post-entrapment melting or crystallization during ascent: Rasmussen et al. (2020) calculate the extent of post-entrapment crystallization to be 1 to 4 wt% (where negative values indicate post-entrapment melting).Compositional changes in response to such small quantities of post entrapment crystallisation or melting are unlikely to significantly affect total water diffusivity.For models with fewer mesh points and which were run at different sets of temperatures (e.g., 2D inversions, see below) we used the original parameterization of Ni and Zhang (2018).

2D Inversions of Natural Melt Inclusion-Bearing Olivine Phenocrysts From Seguam
There are still some uncertainties associated with parameters that control water loss from olivine crystals and their associated melt inclusions.The most important parameters that have first order experimental approximations, but have yet to be widely applied to natural samples, are the diffusive anisotropy and the olivine-melt partition coefficient for water.In order to select parameters that are representative of the 1977 Seguam eruption for our 3D modeling, we ran 2D inversions on the selected crystals which were measured by SIMS.Bulk water measurements were made in the melt inclusion and in profiles close to the [100] and [001] direction.We generated 2D finite element meshes of the crystals using Gmsh (Geuzaine & Remacle, 2017) and crystal outlines generated in ImageJ (Schindelin et al., 2012).The melt inclusions and the host olivine crystal were flagged as different domains in the mesh.We used electron backscatter diffraction (EBSD) measurements of the crystals to calculate the angles between the [100] and [001] axes of the mesh with the main crystallographic axes, which was then used to adjust the diffusion coefficients in olivine down these directions.
The 2D inversions were conducted using DFENS (Mutch, Maclennan, Holland, & Buisman, 2019;Mutch et al., 2021), which employs a Monte Carlo Bayesian method using Nested Sampling (Buchner et al., 2014;Feroz et al., 2009).This approach allows for parameter estimation as well as a means to characterize the associated uncertainties.A description of the numerical methods are the same as those used in the 3D modeling and are discussed in Section 3.3.5.We inverted for 7 different parameters including: decompression rate, H + in olivine diffusive anisotropy, the olivine-melt partition coefficient, the initial water content, temperature, and the underlying parameters that are in the diffusion coefficient of H + in olivine (e.g., the exponent of D 0 and the activation energy).We used log uniform prior distributions for decompression rate (0.01-0.5 MPa s 1 ), diffusive anisotropy (1-100) and the olivine-melt partition coefficient (0.0001-0.002) to get better representation of smaller values over multiple orders of magnitude.The decompression rate prior encapsulates the range of decompression rates already measured in the 1977 Seguam eruption (Newcombe, Plank, Barth, et al., 2020).The anisotropy and partition coefficient priors cover the range of values observed from natural and experimental studies (Barth & Plank, 2021;Hauri et al., 2006;Le Voyer et al., 2014;Newcombe, Plank, Barth, et al., 2020;Towbin et al., 2023).A uniform prior was used for initial water content (measured melt inclusion content to 6 wt%), which captures the variability observed in the melt inclusion population and also accounts for considerable water loss upon ascent.A Gaussian prior was used for temperature (1,070 ± 20 °C) using an estimate based on melt thermometry and the 1σ uncertainties conducted on Seguam matrix glass using the Sugawara ( 2000) thermometer (Newcombe, Plank, Zhang, et al., 2020).Multivariate Gaussian priors were used for log 10 D 0 ( 5.0 ± 0.4) and E a (125 ± 10 kJ mol 1 ), which were obtained from the experimental results of Barth et al. (2019) and error propagation by Barth and Plank (2021).We used our derived covariance matrix to account for the correlated nature of the uncertainties in the diffusion coefficient parameters in the same way as Mutch et al. (2021).
In each simulation we generated a H 2 O solubility curve which was used as the degassing boundary condition during decompression.These were calculated using a Python version of VolatileCalc (Newman & Lowenstern, 2002;Rasmussen et al., 2020) using the major element content of Seguam tephra glass, the initial water content generated from the prior and a CO 2 content of 900 ppm (Moore et al., 2015).We used the basalt version of VolatileCalc (Newman & Lowenstern, 2002) and set the SiO 2 content of the Seguam melt to the maximum accepted value of 49 wt%, given that the Seguam tephra glass slightly exceeds this value (50.55 wt%).We assumed that there was closed system and equilibrium degassing with 2 wt% excess CO 2 at the start of each run to fit with the observations made by Newcombe, Plank, Barth, et al. (2020) and Rasmussen et al. (2020).Each model had 5,000 steps, and the degassing curve terminated at the pressure that corresponded to 0.3 wt% water, which is the water content of glass measured at the edge of Seguam embayments (Newcombe, Plank, Barth, et al., 2020).To ensure that the inversions were completed in a timely manner, each model realization was run with 100 time steps, with the size of the time step being controlled by the corresponding decompression rate.Mesh resolution was maintained at 10-15 μm in the olivine and 7-10 μm in the melt inclusion.

3D Numerical Modeling of Water Diffusion in Melt Inclusion-Bearing Olivine Phenocrysts
Modeling diffusive loss of water from olivine-hosted melt inclusions has typically involved using 1D analytical solutions of the diffusion equation, either using spherical approximations (Chen et al., 2013;Cottrell et al., 2002), or assuming diffusion along a 1D plane where diffusion along the [100] axis is dominant (Barth et al., 2019).3D finite difference models have also been used, in which the geometry of the olivine crystal was approximated as a sphere, cuboid or prism (Le Voyer et al., 2014;Thoraval & Demouchy, 2014).Spherical diffusion models fail to properly account for diffusive anisotropy, whilst 1D numerical attempts (Barth et al., 2019;Newcombe, Plank, Zhang, et al., 2020) may not fully encapsulate the complete 3D geometry of natural or idealized olivine crystals (Shea et al., 2015;Welsch et al., 2012).Furthermore, previous modeling attempts assumed that diffusivity of water in the melt inclusion was infinitely faster than diffusion in the olivine, meaning any changes in diffusivity in the melt inclusion due to changes in composition or intensive parameters were not fully accounted for.
Geochemistry, Geophysics, Geosystems 10.1029/2023GC011365 Solving partial differential equations using finite elements has started to gain prominence in diffusion chronometry applications in igneous petrology (Mutch, Maclennan, Holland, & Buisman, 2019;Mutch, Maclennan, Shorttle, et al., 2019;Mutch et al., 2021).This is because finite elements offer an efficient way of solving the diffusion equation for complex geometries, such as those observed in natural crystals.In this study, we use FEniCS (Alnaes et al., 2015), a free open-source Python-based software that can exploit parallel processing, to develop a new 3D finite element diffusion model that can track the diffusive loss of water from olivine-hosted melt inclusions.A full description of the numerical implementation of FEniCS to diffusion modeling in igneous systems can be found in Mutch et al. (2021).The present model includes two separate domains meaning it can be used to track changes in water concentration across the melt inclusion and the host olivine crystal.Diffusion across two domains was modeled as follows: where C, D and ρ denote the concentration of water, the bulk diffusivity of water, and the density of the melt inclusion (Ω MI ) and olivine (Ω Ol ) domains as marked by the respective subscripts.
The two domains are separated by a boundary Γ.On this boundary we have the conditions.
where the first of these conditions represents the partitioning of an element with partition coefficient K, and the second represents conservation of mass across the interface.
We can then introduce the following rescaling: With this scaling the system, we can represent the problem with a single diffusion equation in both domains, with standard continuity relationships across the interface Γ, Geochemistry, Geophysics, Geosystems 10.1029/2023GC011365 We can then solve Equation 11 by standard methods with a spatially varying ρ and D, without needing to treat the interface in a special manner.This kind of approach, in which there is a mass balance at the interface is also applied in multiphase exchange thermometry and chronometry in igneous and metamorphic systems (Müller et al., 2013).Like previous studies (Cottrell et al., 2002), we assumed that the density ratio of olivine to melt was equal to 1.2.We verified our model against 1D analytical solutions from Zhang (2009), which is shown in the Supplementary Information S1.
Linear continuous Galerkin finite elements were used to represent concentration, whilst a discontinuous Galerkin function space (DG0) was used for the spatially variable ρ and D fields.A fixed Dirichlet boundary condition was maintained at the exterior boundary of the crystal and was modified based on the composition of olivine in equilibrium with the exterior melt as calculated by the basalt version of VolatileCalc (Newman & Lowenstern, 2002;Rasmussen et al., 2020).Across both domains, olivine equilibrium compositions were calculated.
The composition in the melt inclusion was then calculated using the partition coefficient at the end of each time step.The diffusion equation was solved at each time step using an iterative Newton solver.A Crank-Nicholson (second order) time-stepping scheme was used in the models.The number of time steps was typically set to 300, with the size of the time step being adjusted according to decompression rate.An efficient algebraic multigrid preconditioned conjugate gradient solver was used to solve the resulting systems of linear equations in 3D (Balay et al., 2023).
Finite element meshes were generated using Pygmsh and Gmsh version 3.0.0(Geuzaine & Remacle, 2017).An idealized morphology of the olivine crystal structure was based on Welsch et al. (2012).We refined the mesh at the boundary between the olivine and the melt inclusion.Meshes typically had 400,000 vertices.This was to balance spatial resolution with computational efficiency.The accuracy of finite element analysis depends on the mesh resolution and sampling across the physical domain.We tested for convergence by running multiple models at different spatial resolutions (see Supporting Information S1).For the mesh resolutions that we employed in our models (4-12 μm), the model converged within 0.05 wt% of the true solution.Once the mesh was defined, we then labeled the separate domains and boundaries accordingly.

Seguam Crystal Cargo
Seguam olivines show prominent zoning in bulk water content, which decreases away from the large central melt inclusions toward the crystal edge (Figure 2).Water concentrations are typically 13-24 ppm next to the melt inclusion and are 4-6 ppm next to the crystal edge.Analyses measured directly next to the melt inclusions may be elevated in water content due to contamination from the melt inclusion itself.The zoning profiles are most prominent along [100] (the a-axis) and gradually decrease from the melt inclusion to the crystal edge.
Profiles measured along [001] (the c-axis) show sinusoidal profile shapes with steep gradients close to the melt inclusion and crystal edge and a plateau in between.This is most prominently shown in Seg1-MI1.Melt inclusion water compositions range from 3.5-3.84wt%.Zoning in forsterite content is relatively minor with crystal cores of Fo 83 and thin normal zones extending to Fo 81-82 at the crystal edge.Seg4-MI1 is the exception, with reverse zoning from the crystal center (Fo 82 ) to a mantle of Fo 83 which in turn is surround by a normally zoned outer rim.

Seguam 2D Inversion Results
We performed 2D inversions on three melt inclusion-bearing phenocrysts from Seguam (Seg1-MI1, Seg4-MI1, Seg13-MI1).The models that occupy the maximum likelihood space fit the data well (Figure 2), with the general exception of points adjacent to the melt inclusions.These analyses may have been contaminated by the melt inclusion giving them their elevated water concentration, or they may be the product of a partitioning mechanism that is not incorporated in the current model framework (Lynn & Warren, 2021).
Each crystal generally required 10,000 to 40,000 Monte Carlo realizations to reach convergence.Trade-offs in the posterior distributions of the seven inverted parameters for one of the crystals (Seg13-MI1) can be observed in Figure 3. Strong trade-offs can be observed between decompression rate, temperature and the diffusion coefficient parameters.Minor trade-offs between the decompression rate, the partition coefficient and initial water content have also been observed, with weak negative correlations between initial water content and the partition Geochemistry, Geophysics, Geosystems coefficient, and decompression rate and the partition coefficient.Given all of these parameters play a major role in controlling water loss from the melt inclusion, changes in one parameter can come at the expense of the others.For example, higher temperatures are associated with faster diffusivities and thus require a higher decompression rate to match any potential water loss.Conversely, correlations with the partition coefficient likely reflect the model's ability to fit the SIMS profile data.
Inversion results for all three crystals are shown in Figure 4. Posterior distributions for temperature and the diffusion parameters (D 0 and E a ) typically agree with their Gaussian priors, and thus reflect their initial uncertainties.Decompression rate, anisotropy and partition coefficient posteriors fall within their respective priors and show the greatest inter-crystal variability which indicates they are being controlled by the data.
Estimated linear decompression rates show some agreement with the range estimated by Newcombe, Plank, Barth, et al. (2020) by volatile diffusion in embayments (0.02-0.13 MPa s 1 ) and water gradients in olivine phenocrysts (0.04-0.23 MPa s 1 ).Seg13-MI1 falls within this range with a median decompression rate of 0.07 MPa s 1 and with a 1σ uncertainty of 0.01 MPa s 1 .Seg1-MI1 and Seg4-MI1 have faster decompression rates than those measured by Newcombe, Plank, Barth, et al. (2020) et al. (2019).These anisotropy distributions have long tails that extend to higher anisotropies.These tails are likely artifacts that reflect the model's lack of sensitivity to high anisotropies given the spatial resolution of the SIMS profiles.Higher resolution analyses will be needed to get a more accurate assessment of the diffusive anisotropy.Finally, initial water contents generally converged toward the melt inclusion value, which reside at the lower end of the priors and is below the estimated maximum for Seguam (4.2 wt%).This suggests that the melt inclusions lost nearly no water during magma ascent.The models, however, may be preferentially converging to lower initial water content values given the large uncertainties in other free-floating parameters such as the partition coefficient, diffusion coefficient and anisotropy.If these parameters can be associated with tighter priors, then initial water may be more accurately constrained by these type of inversions.

Variability in Inverted Parameters at the 1977 Eruption of Seguam
Of the three olivine crystals from the 1977 eruption of Seguam that we conducted 2D inversions, we have identified two distinct types based on inverted decompression rates and geochemical parameters associated with diffusion, notably diffusive anisotropy and the partition coefficient.Seg1-MI1 and Seg4-MI1 show consistent decompression rates, which are a factor of 5 greater than the rate determined by Seg13-MI1.Seg13-MI1 and Seg1-MI1 show consistent partition coefficients and to a lesser degree diffusive anisotropy, which differ to those estimated for Seg4-MI1.In the case of anisotropy, our inversions suggest that H + diffusion in Seg4-MI1 is almost isotropic, which conflicts with recent experimental observations (Barth et al., 2019;Ferriss et al., 2018).However, other low anisotropy examples have also been found in the Seguam crystal cargo by Newcombe, Plank, Barth, et al. (2020).The discrepancies in inverted parameters that we observe across the Seguam crystal cargo could be due to several different reasons.They could relate to the underlying major and minor element chemistry of the olivines.Forsterite profiles measured across each crystal parallel to the SIMS profiles show that Seg1-MI1 and Seg13-MI1 have relatively flat cores between Fo 83 and Fo 84 and a simple outer rim, whilst Seg4-MI1 contains an inner core of Fo 82 surrounded by a reversely zoned mantle (Fo 83 ) and then an outer rim (Figure 2).These different zoning patterns not only suggest that these crystals may have undergone slightly different magma ascent histories, but also they may have different defect populations that could affect the diffusion of H + (Jollands et al., 2019).This highlights the limitation of applying a single Arrhenius relationship that has been calibrated on a limited range of compositions to chemically variable natural crystal cargoes.These discrepancies could then be associated with caveats in our diffusion modeling methods which are discussed in more detail below.
Variations in decompression rate observed across the modeled crystal population could also be the result of limitations in the diffusion modeling, or could reflect natural physical processes.Magma may undergo complex ascent histories in which velocity may change across the width of the conduit, conduit flow regime may change, turbulence might be induced by rough conduit walls, or the magma may stall at different points in the conduit (Gonnermann & Manga, 2007;Myers et al., 2019).Water in olivine crystal chronometry is sensitive to changes in the local environment, meaning individual crystals may be recording different aspects of these complexities.For example, entrained crystals that ascended in different parts of the conduit may record different decompression rates.Alternatively, crystals erupted at different points during a single eruption may record changes in ascent histories.Due to the large uncertainties associated with individual decompression rates of each crystal, these individual histories may be difficult to disentangle.Improving the uncertainties on the diffusion coefficients by further understanding the associated diffusion mechanisms could help to resolve this issue.

Monte Carlo Simulation of Magma Ascent
Crystal geometry and diffusive anisotropy have been shown to play an important role in diffusion modeling of major (e.g., Fe-Mg interdiffusion) and minor elements (Ni, Mn, Ca) in olivine (Mourey et al., 2023;Shea et al., 2015).Given the very high diffusivity of H + in olivine, which can be up to eight orders of magnitude faster than Fe-Mg interdiffusion (Dohmen & Chakraborty, 2007), multi-dimensional and geometrical considerations may become even more important for water loss or gain during magmatic processes.Recent work has argued that 1D models of diffusive loss of water from melt inclusions and the host olivine along the [100] direction are the best approximations given the very high diffusive anisotropy observed in dehydration experiments of Fe-bearing olivines (Barth & Plank, 2021;Barth et al., 2019;Ferriss et al., 2018;Thoraval & Demouchy, 2014).The modeling framework proposed by Barth et al. (2019) for water loss from olivine-hosted melt inclusions is a 1D model from the center of the melt inclusion to the nearest crystal edge along the [100] crystallographic axis.This model applies a reflective Neumann boundary condition at the center of the melt inclusion and a fixed Dirichlet boundary condition at the edge of the crystal.The model therefore assumes that crystal is symmetrical about the melt inclusion (i.e., the melt inclusion is at the center of the crystal) and that the flux out of the crystal is dominated by the shortest [100] direction along the ascent pathway.Here we applied a Monte Carlo approach with the aim to understand the role of 3D crystal geometry and melt inclusion configuration on estimates of water loss from melt inclusions during magma ascent.This approach was then used to test the efficacy of 1D and 2D numerical models as well as analytical solutions at constraining melt inclusion water loss and thus retrieving decompression rates.

Model Set Up
We based our models on physical and chemical observations of the 1977 Seguam eruption.Table 1 shows all of the parameters that were used in the Monte Carlo models, as well as their associated ranges.Parameter distributions are shown in Figure 5.A full description of all of the selected parameters can also be found in Supporting Information S1.Melt inclusion position along the main crystallographic axes show non-uniform distributions, particularly along the [010] and [001] directions.This is because we applied conditions to prevent melt inclusions from being positioned too close the crystal edge.We set the maximum distance along each crystallographic axes to be crystal axis length minus the melt inclusion radius plus a tolerance.The tolerance was one quarter of the melt inclusion radius.We then implemented an algorithm to ensure that the mesh generation failed if the melt inclusion intersected any crystal edge.The smaller length of the [010] axis means that most meshes were successfully generated if the melt inclusion was positioned toward the center of the [010] axis.This was particularly true if the melt inclusion radius was large with respect to the [010] axis length.Olivine crystals are typically tapered toward the ends of the [001] and [100] axes.This means there is a higher chance for larger melt inclusions to intersect crystal edges, and that meshes in which the melt inclusions were closer to the center of the crystal pass the edge intersection criteria and are thus successfully generated.
In total, 1015 model configurations were run with each using the three different anisotropies.Of these, 463 model sets were run using the Seguam-derived partition coefficient (0.000459), 265 model sets were run using a partition coefficient of 0.001, and 287 model sets were run with a randomly selected partition coefficient.In each Monte Carlo simulation a 3D forward model was run and the composition at the center of the melt inclusion was tracked through time (Figure 6).The water content at the end of the decompression path was recorded.The mesh resolution was determined by both the crystal size and melt inclusion radius in order to balance accuracy and computational efficiency.The mesh point spacing was determined to be 60 times smaller than the length of the [001] axis from the center of the crystal.This resulted in mesh resolutions of 4-12.5 μm, which is comparable to other 3D modeling efforts that have explored H + diffusion in olivine (Lynn & Warren, 2021).The mesh resolution in the melt inclusion was determined to be one twelfth of the melt inclusion radius, which resulted in mesh resolutions of 0.8-8.3μm.The melt inclusion was surrounded by a mesh of the same resolution, which extended to a distance based on the melt inclusion radius.This then transitioned to the mesh resolution of the host olivine over the lengthscale of the melt inclusion radius.In each model a constant number of 300 time steps were used, which resulted in individual time steps varying from 2.1 s to 3.8 hr depending on the decompression rate and starting pressure.
We then ran equivalent 1D and 2D models using the same physico-chemical parameters in each 3D Monte Carlo model.The same mesh resolutions and time steps were used to maintain consistency with the 3D models.We use two types of 1D model termed here symmetric and asymmetric models.For the symmetric models, we generated a mesh that extends from the center of the melt inclusion to the closest edge of the crystal along the [100] direction.
In this instance it is assumed that the crystal along the [100] axis is completely symmetric around the center of the melt inclusion, which is implemented using a Neumann or no-flux boundary condition at the center of the melt inclusion.This is to simulate the approach recently adopted by Barth et al. (2019).For the asymmetric models (Figure 6), we generated a mesh that extends from one crystal face to the other along the [100] axis, and which goes through the center of the melt inclusion.For models in which the melt inclusion is not in the center of the Note.Additional information and justification for each parameter is also included.L [001] is the length of the [001] axis from the center of the crystal.
Geochemistry, Geophysics, Geosystems 10.1029/2023GC011365 crystal, this creates an asymmetric 1D arrangement along [100].In this model type, Dirichlet boundary conditions were applied at the edges of the crystal.We also generated two types of 2D model, the first uses a 2D section of the crystal that includes the center of the melt inclusion using the (010) plane, whilst the second uses the (001) plane (Figure 6).In these 2D models, the corresponding anisotropy was used.
The 1D and 2D models were run under the same conditions as the 3D models, and the water content at the center of the melt inclusion at the end of the decompression pathway was recorded for comparison.We also ran Nested Sampling Bayesian inversions in order to extract magma decompression rates.In this case we used the final melt inclusion water content of the 3D models for each anisotropy to fit the 1D and 2D models.The inversions were conducted using PyMultinest (Buchner et al., 2014;Feroz et al., 2009) with the only invertible parameter being decompression rate.All other physico-chemical parameters were fixed to the values used in the 3D forward models.We used a log uniform prior for decompression rate which encapsulated a much wider range than was used in the Monte Carlo forward models (0.00001-0.5 MPa s 1 ).The decompression rate that maximized the log likelihood function was then used as the inverted decompression rate.

Comparing 1D and 3D Numerical Models
In general, the 1D models underestimate the amount of water loss when compared to the 3D models (Figure 7).The differences in the 3D versus 1D models are lowest at low or high percentages of water loss.The former are likely to be associated with high decompression rates or large melt inclusions meaning there is insufficient time for discernible water loss between the two types of models.The latter are most likely associated with smaller inclusion sizes and low decompression rates where the melt inclusion is approaching equilibrium.This discrepancy between model types is greatest when the 1D models are compared to 3D models conducted at low anisotropy (3.0), but diminishes for 3D models conducted at higher anisotropy (15.6 and 35).For example, the 1D models typically under predict melt inclusion water loss by up to 40%-50% when compared to 3D models that use an anisotropy of 3.0, whereas this typically falls within 30% for the 15.6 and 35 anisotropy models.This effect can simply be explained by the extra amount of water loss along the extra dimensions in the crystal.The effect would be greatest at lower anisotropies as faster diffusion along these extra dimensions would not be properly accounted for in the 1D models.
A similar effect can be observed in inverted decompression rates, where 1D models underestimate decompression rates relative to the 3D models (Figure 8).This is because the 1D models lose less water than their 3D counterparts, therefore lower decompression rates are required by these models to match the final water content of the 3D models.At high diffusive anisotropies (15.6 and 35) nearly all of the inverted decompression rates from both sets of 1D models falls within a factor of 10 of the 3D model decompression rates.For the 1D symmetric models 88%-94% of inversion results were within a factor of 5 of the 3D decompression rates, whilst 78%-90% of the 1D asymmetric models fell within a similar range.Only 40%-56% of 1D symmetric models and 14%-26% of 1D asymmetric models gave inverted decompression rates within a factor of 2 of the 3D decompression rates.At low diffusive anisotropies (3.0), the 1D models generally performed more poorly with 47% of 1D symmetric and 27% of 1D asymmetric decompression rates falling within a factor of 5 of the 3D results.
The shortest distance along [100] from melt inclusion edge to crystal edge (MI [100] ) normalized to [100] axis length (L [100] ) has one of the largest effects on 1D inverted decompression rates relative to constrained 3D value (defined as a ratio, R 1D/3D ).For melt inclusions situated close to the center of the crystal (MI [100] /L [100] is close to 1), the R 1D/3D of both asymmetric and symmetric models is similar (Figure 9).For models with low anisotropy (3.0), R 1D/3D is approximately 0.1-0.2.At higher anisotropies (15.6 and 35), R 1D/3D is 0.3-0.5.In this geometrical configuration, the symmetric and asymmetric models are almost identical.As the melt inclusion is moved closer to the crystal edge (MI [100] /L [100] approaches 0), the R 1D/3D of the asymmetric models approaches 0.75 for the low anisotropy configuration and approaches 1.0 for the high anisotropy configurations.This is because the diffusive flux along the short dimension of [100] becomes dominant and diffusive loss from other directions becomes less significant.In this geometrical configuration, the 1D approximation can create a more accurate result.
The 1D symmetric models actually overpredict decompression rate as the melt inclusion moves closer to the crystal edge along [100] (i.e., the R 1D/3D approaches values above 1.0).This effect is likely due to the trade-off between the water loss associated with extra dimensions and the artificial increase in the flux out of the melt inclusion in symmetric models by assuming that the crystal is symmetrical about the melt inclusion.If the melt inclusion is positioned off-center and close to the edge of the crystal (i.e., the 1D profile is highly asymmetrical), the symmetric model will significantly underestimate the distance from the center of the melt inclusion to the furthest crystal edge.This means that more water will be lost over the same time period because the length scale of diffusion is shorter in the symmetric model.Melt inclusion size relative to the crystal size also affects 1D models ability to replicate the 3D decompression rates (Figure 9).Smaller melt inclusions tend to underpredict decompression rates to a greater degree than larger inclusions when compared to the 3D models.This is the case for both 1D symmetric and asymmetric models, and is because smaller melt inclusions are more sensitive to water loss from the other dimensions.Furthermore, larger melt inclusions are more likely to have a shorter distance between the edge of the inclusion and the edge of the crystal along the [100] axis, which would increase the validity of a 1D approximation.To try and account for geometrical effects, such as inclusion location and size, on 1D inverted decompression rates we fitted empirical models to the distributions shown in Figure 9 using the form: where the coefficients (a 1 -a 4 ) for each parameter and anisotropy are shown in Table 2, MI [100] is the location of the melt inclusion along [100], r MI is the radius of the melt inclusion, and L [100] is the length of the [100] axis from the center of the olivine to the edge.This empirical fit generates a scaling factor (R 1D/3D ) that could potentially correct decompression rates calculated by the 1D symmetric (Barth et al., 2019) and asymmetric methods.This can be used as a practical tool to generate more realistic decompression rates and may help to reduce the amount of scatter in decompression rates estimated by 1D symmetric and asymmetric melt inclusion models.

Comparing 2D and 3D Numerical Models
Two dimensional models also consistently underestimate water loss (Figure 7) and decompression rates (Figure 8) compared to 3D models.Given the computationally intensive nature of the 2D inversions, only 49 sets were completed.Overall, the 2D models produce predictions that are closer to the 3D models and are associated with tighter distributions than the corresponding 1D models.For example, nearly all of the 2D models return water loss values within 20% of the 3D values regardless of diffusive anisotropy, decompression rate and melt inclusion configuration.Meanwhile, 60%-70% of the 2D inverted decompression rates are within a factor of 2 of the 3D decompression rates, albeit with a much smaller sample size.This is because the 2D models provide a more accurate representation of crystal and melt inclusion geometry than 1D models, and can also account for water loss in two dimensions rather than one.The ability of 2D models to faithfully capture water loss from a melt inclusion therefore depends on the flux of water along the additional unconstrained dimension.This will largely depend on the location of the section in which the 2D model is conducted, and its distance from additional edges of the crystal in the third dimension.It is therefore more crucial for the 2D models to capture the water loss along the shorter dimensions and those with higher diffusivities.In a typical idealized olivine crystal, this would be the (001) sections.
The choice of 2D section may also depend on what other elements are being measured and modeled in the same crystal.Preferred sections for Fe-Mg interdiffusion need to include the [001] axis, as diffusion is 6 times faster in this direction (Dohmen & Chakraborty, 2007;Shea et al., 2015).The (001) section would therefore be non-ideal, which conflicts with the desired 2D section for H + modeling.In this instance, a compromise would need to be made if Fe-Mg and H + need to be modeled in the same crystal.The best solution would be to use the (010) section, as this contains both the [001] axis for Fe-Mg interdiffusion and the [100] axis for H + diffusion.The benefits of having Fe-Mg well constrained would probably outweigh the loss of accuracy of H + given that the diffusive anisotropy of Fe-Mg is much lower and that 2D H + can still be conducted on the (010) plane.

Comparing Anisotropic Analytical Solution and 3D Numerical Models
Analytical solutions offer the most efficient means for calculating diffusive water loss from olivine-hosted melt inclusions (Bucholz et al., 2013;Gaetani et al., 2012).These methods are, however limited by simplifying assumptions about geometry and diffusive anisotropy.For example, many studies have used the analytical solution of Qin et al. (1992), which assumes isotropic diffusion in a spherical host mineral.Given the potentially highly anisotropic nature of H + diffusion in olivine, analytical solutions that capture this behavior may be more Note.This transforms decompression rates estimated using 1D symmetric and asymmetric models into an equivalent 3D estimate.
appropriate.Analytical solutions for the diffusive equilibration of an elliptical inclusion in an unbounded host mineral can be adapted to solve diffusive equilibration of a spherical melt inclusion in an anisotropic host.The equilibration timescale, τ, can be expressed as: where ρ Ol and ρ MI are the densities of olivine and the melt inclusion respectively.D [100] , D [010] and D [001] are the olivine diffusivities along the [100], [010] and [001] directions.R F is a Carlson completely symmetric elliptic integral of the first kind, and r MI is the melt inclusion radius.The derivation of this solution is available in Supporting Information S1.This analytical solution assumes that the melt inclusion is sufficiently small that it sees the crystal as an effectively infinite domain, and that diffusion in the inclusion is rapid, such that the diffusivity of H + in the crystal sets the rate of equilibration.We can use the estimate of equilibration time scale to approximate the behavior of water in the melt inclusion when the boundary conditions change over time.If C i is the melt inclusion water concentration and C 0 is the concentration of water in the surrounding melt, the equilibration on a time scale τ can be represented by the simple ordinary differential equation: If C 0 (t) were constant, this would be exponential decay.The solution of the above ordinary differential equation can be written in integral form as and the integral can be approximated by the trapezoidal rule if C 0 (t) is given, as would be the case for a water solubility model.
In terms of melt inclusion water loss, the analytical solution presented in Equation 15 also underestimates water loss relative to the 3D numerical models.It performs relatively well with over 90% of the water loss estimates from the analytical solution falling within 10% of the 3D models regardless of the degree of equilibration and the diffusive anisotropy (Figure 8).These underestimations also translate into decompression rates for the same reasons as the 1D and 2D numerical models as discussed above.Over 90% of inverted decompression rates fall within a factor of 5 of the 3D values, whilst 50%-80% fall within a factor of 2. The performance of the analytical solution appears to improve with decreasing anisotropy.Like the 1D numerical models, we can then compare the ratio of inverted decompression rates from the analytical models and original 3D decompression rates (R A/3D ) to parameters related to melt inclusion size and configuration (Figure 9).Melt inclusion position along the [100] axis and melt inclusion radius play a major role in determining the accuracy of the inverted decompression rate from the analytical solution.The reason the analytical solution underestimates water loss is likely related to the assumption that the olivine crystal is an infinite domain, which requires the melt inclusion to be sufficiently small not to feel far-field effects.If the inclusion is large relative to the crystal, and it is situated close to the edge along the [100] axis, then the analytical model will underestimate water loss.
The anisotropic analytical solution and 1D numerical models could offer complementary ways to obtain magma decompression rates for different olivine-melt inclusion configurations that retain some fidelity to the 3D numerical solution but are also computationally more efficient.If the melt inclusion is close to the center of the crystal and is small relative to the host crystal then the analytical solution (Equation 15) could be used to extract decompression rates within a factor of 2. Conversely, if the melt inclusion is situated close to the crystal edge along the [100] axis, then the 1D approximations (with our empirical correction) could offer an accurate solution.2D inversions could offer a solution that could be applied over a much wider range of geometrical configurations, provided the effects of diffusion along the third dimension can be mitigated.However, 2D models are more computationally intensive.All of the model configurations and subsequent recommendations discussed thus far are valid for crystals that contain a single melt inclusion in a crystal with an idealized olivine morphology.These may need to be modified if the crystal contains multiple inclusions or has a more complex crystal shape.

The Role of Multiple Melt Inclusions in Controlling Water Loss or Gain
Most attempts to model water loss from melt inclusions assume that the olivine contains a single melt inclusion located in the center of the crystal (Bucholz et al., 2013;Chen et al., 2013;Qin et al., 1992).This has primarily been done to simplify the modeling (e.g., radial analytical or numerical solutions), and could reliably be applied to natural systems in some instances (Barth et al., 2019;Chen et al., 2013;Gaetani et al., 2012).Many natural olivine crystals, however, can contain multiple melt inclusions with different sizes and configurations (Wallace et al., 2021).This depends on the growth history of the olivine, which may undergo periods of crystal growth, resorption and regrowth (Wallace et al., 2021).Rapid dendritic growth induced by high degrees of undercooling to form skeletal or hopper crystals followed by slow maturation of crystal faces can cause melt inclusion entrapment (Faure et al., 2003;Mourey & Shea, 2019;Welsch et al., 2014).Alternatively, melt inclusions can form as a result of irregular growth due to the attachment of small phases (e.g., spinel) to the surface of the crystal, or from the sealing off of melt embayments formed by partial resorption (Roedder, 1979).Combinations of these processes mean that some crystals can contain hundreds of melt and fluid inclusions (Métrich & Wallace, 2008;Wallace et al., 2021).
Given the diversity of melt inclusion configurations observed in the crystal cargoes of single eruptions (Wallace et al., 2021), it is therefore important to be able to understand water loss or gain from olivine crystals that contain multiple melt inclusions.This may be critical for getting a representative view of the crystallisation and decompression history of a given magma or volcanic system.The 3D multiphase finite element model that we have developed as part of this study allows us to incorporate multiple melt inclusions in order to understand their interaction during magma ascent.Here, we ran additional forward models to compare water loss from olivine crystals that contained a single central spherical melt inclusion (Figure 10) with models in which this central inclusion is surrounded by 18 other spherical melt inclusions (Figure 11).Three variations of the models were run where we varied the size of central melt inclusion (radius of 25, 50, and 75 μm).The size of the surrounding melt inclusions were fixed with a radius of 50 μm.The crystal size (695 μm along the [001] axis from the center of the crystal), the olivine-melt partition coefficient (0.000459), initial water content (4.2 wt%), final matrix glass water content (0.3 wt%) were fixed to representative values for the 1977 eruption of Seguam.Decompression rate was varied from 0.001 to 0.8 MPa s 1 .To assess the role of diffusive anisotropy in olivine on melt inclusion interaction, we ran additional models that compare water loss from a single 25 μm single central melt inclusion with models in which there is an additional 75 μm inclusion positioned along [100], [010], or [001] (Figure 12).All of these model configurations were run with different anisotropies (3.0, 15.6, 35.0) and decompression rates.
Our models show that the addition of multiple melt inclusions could have a large effect on melt inclusion water loss during magma decompression.Isolated melt inclusions (i.e., a single melt inclusion in a crystal) can lose more water during magma ascent than melt inclusions that are adjacent to or are surrounded by neighboring inclusions (Figure 13).This difference can range from 0.5 wt% to over 1.5 wt% depending on melt inclusion size and configuration.3D model results show that the "haloes" of elevated water content surrounding the melt inclusions overlap and interact (Figures 11 and 12).Our results demonstrate that large melt inclusions are able to buffer water loss from adjacent small melt inclusions, and clusters of melt inclusions are able to "shield" each other from the effects of diffusive water loss in response to degassing of the host magma.Furthermore as the water content of the melt inclusions drop (particularly for low decompression rates), the diffusivity of water in the melt could potentially decrease to below that of the olivine along the fast [100] direction.Hence the melt inclusions could act as a slower pathway for water to escape.
The effect that multiple melt inclusions have on water loss appears to also depend on both decompression rate and melt inclusion size (Figure 13).This shows a near-Gaussian distribution with a peak at decompression rates around 0.001-0.01MPa s 1 , and likely relates to the amount of time in which diffusion is allowed to operate.At high decompression rates, there is insufficient time for diffusion to remove water from the central melt inclusion in either model scenario (i.e., the diffusion fronts have not approached the center of the crystal).Meanwhile at low decompression rates, there is plenty of time for diffusion to occur meaning the central melt inclusion is close to Geochemistry, Geophysics, Geosystems 10.1029/2023GC011365 MUTCH ET AL.
equilibrium for both model scenarios upon eruption.For melt inclusion shielding to reach full effect, the magma needs to be ascending at a rate in which diffusion has sufficient time to operate, but also for disequilibrium to maintained.This would occupy the speedometer domain described by Barth and Plank (2021).
The size of the shielding effect is sensitive to melt inclusion size (Figure 13).The difference in water loss from the central melt inclusion is greater if it has a smaller radius.This is because smaller melt inclusions are more sensitive to water loss.The shielding effect is more pronounced in larger central melt inclusion models (50-75 μm radius) that were run using low anisotropies (3.0).This is because in these models the crystal loses more water because diffusivity along the [010] and [001] is faster meaning the larger melt inclusions can be more sensitive to water loss.Given the greater utility of smaller melt inclusions in ascent rate speedometry, knowing about potential shielding effects from nearby melt inclusions is even more important.Our models may be extreme examples, particularly with regards to melt inclusion sizes, but they do highlight the importance of considering multiple melt inclusions for magmatic applications.For hygrometry, melt inclusions that are surrounded by other melt inclusions, particularly along the [100] direction, may provide more faithful estimates of initial water content provided the crystal has not undergone complete equilibration.For speedometry, it may be crucial to incorporate multiple melt inclusions in modeling if they are situated in the same plane along the fast direction.

Caveats, Complexities and Recommendations
The diffusion models in this study are simplified and based on a series of approximations and assumptions.This approach may be sufficient for understanding how melt inclusion interaction and geometry may impact water loss or gain during simple ascent histories, however it may not completely capture complexities that are observed in natural systems.Here we list some potential caveats and complexities associated with model assumptions and natural systems, but also pragmatic recommendations based on the insights of our modeling and parameter exploration.

Non-Diffusive Modification of Melt Inclusion Water Content
We assume that the water content of the melt inclusion is solely controlled by diffusive loss through the host olivine during decompression.Furthermore, the melt inclusion maintains a fixed volume.Our models predict that there should be water zoning within the melt inclusion that occurs in response to this process.In reality, there are additional processes that can modify melt inclusion water content and volume upon ascent including: • Vapor bubble growth.Water may partition into a growing vapor bubble during the final stages of decompression, meaning it may be unaccounted for if the bubble composition is not measured.Water gradients may develop around vapor bubbles and likely mark the complexities associated with water diffusing into the host crystal and the vapor bubble, which may be further compounded by the bubble moving around in the inclusion.• Post-entrapment crystallisation.As the magma ascends to the surface, changes in P-T conditions may affect phase stability and cause crystallisation or dissolution on the melt inclusion walls, which may in turn cause the size of the inclusion to change and for the concentration of water in the melt inclusion to increase or decrease.This process will depend on the magma composition and the ascent history, and could potentially be reconstructed by looking at cooling histories of melt inclusions (Newcombe et al., 2014;Newcombe, Plank, Zhang, et al., 2020).Thermal history modeling applied to MgO profiles measured in Seguam olivine hosted melt inclusions suggests they did not undergo any net cooling upon ascent and that any growth on the inclusion walls was likely quench crystallisation (Newcombe, Plank, Zhang, et al., 2020).

Diffusion Mechanisms
The choice of the H-in-olivine diffusion coefficient can also be associated with significant caveats.These include: • Extrapolation of the diffusion coefficient beyond experimental conditions can introduce further uncertainty given that the exact dependence of bulk H + diffusivity on olivine composition has yet to be properly quantified (Barth & Plank, 2021).Iron content has been shown to play a major role (Barth et al., 2019(Barth et al., , 2023)), but the concentration of trace elements and water itself could be important (Jollands et al., 2019;Tollan et al., 2018).
Lithium has shown the potential to couple with trace element zoning patterns, most notably phosphorus (Lynn et al., 2020); H + could show similar behavior.• Consideration of intersite reaction and the changing availability of defects (Ferriss et al., 2018;Jollands et al., 2019).Reaction rates and site availability for defects can change over time, meaning diffusivity could change as dehydration progresses (Ferriss et al., 2018).• Uncertainty associated with diffusive anisotropy.Experimental observations (Barth et al., 2019;Ferriss et al., 2018), and some of our 2D inversions indicate high anisotropies (10-40) whilst melt inclusion-free crystals from Seguam do not appear to show evidence of high anisotropy, with estimated values of 2-3 (Newcombe, Plank, Barth, et al., 2020).It is unclear what could cause this discrepancy.Barth and Plank (2021) suggest that the conditions under which dehydration experiments are conducted may not fully replicate the processes that occur in nature.Alternatively, limitations in analytical precision, drift and spatial resolution (particularly along the slow directions), along with trace element coupling may generate profiles with apparent high or low anisotropy.

Partitioning Behavior
The partition coefficient of water between the olivine and the melt can exert a large control on water loss and the shape of diffusion of bulk H + profiles in olivine.We have assumed a constant partition coefficient based on measured profiles in natural Seguam crystals, however additional behavior may need to be considered including: • Different incorporation mechanisms have been suggested for H + in olivine that may dominate under different P-T conditions and mineral compositions (Barth & Plank, 2021;Danyushevsky et al., 2002;Gaetani et al., 2012;Portnyagin et al., 2019).If this were indeed the case, the partition coefficient would be continually changing along the decompression pathway.• We assumed a constant partition coefficient that is based on H 2 O equivalents in the olivine and melt, which does not consider the effects of changing speciation in the melt.Hydrogen in olivine dissolves as OH groups and hydrogen in melt dissolves as both molecular H 2 O and OH (Demouchy & Mackwell, 2006;Stolper, 1982).Different partition coefficients for different species may need to be considered if different species become dominant along the ascent pathway.• Non-instantaneous and non-equilibrium partitioning.Our models assume that equilibrium between the crystal and the external melt is obtained instantaneously when the boundary condition is being changed in response to degassing.Whether equilibrium partitioning is attained instantaneously will depend on the size of the concentration jumps in each step, which will be controlled by ascent rate or style of degassing.This issue can be dealt with by introducing a kinetic parameter that can account for fractional attainment of equilibrium (Dohmen & Chakraborty, 2003).

Style of Magma Ascent and Degassing
Assumptions have been made in the modeling about the style of magma ascent and degassing which may not be appropriate to some natural systems.We assumed a linear decompression rate and equilibrium closed-system degassing as an external boundary condition.Additional complexities include: • Magma decompression is expected to be highly non-linear (Barth et al., 2019;Gonnermann & Manga, 2007;Hajimirza et al., 2021;Su & Huber, 2017).Exsolution of water and subsequent nucleation and growth of bubbles at shallow depth increases the buoyancy of the magma and drives acceleration (Gonnermann & Manga, 2007).Further changes in magma viscosity due to crystallisation and water loss, in addition to friction along the conduit wall may need to be considered.Models that assume a constant decompression rate can underestimate the total ascent time because more time has to be spent at depth in order to compensate for the reduced amount of time at shallow depths where most water loss takes place (Barth & Plank, 2021;Barth et al., 2019;Hajimirza et al., 2021).The degree of non-linearity would need to be known a priori, potentially from the predictions of conduit models, in order to constrain faithful ascent timescales.• The assumption of equilibrium degassing may not be valid in some systems.Equilibrium degassing may be acceptable for basaltic systems, where bubble nucleation is not necessarily inhibited by the low viscosities.For systems with higher viscosities, such as andesites and rhyolites, consideration of disequilibrium degassing may be more important.deGraffenried and Shea (2021) showed that this is the case for modeling decompression rates from melt embayments in quartz crystals in rhyolitic systems.For simple decompression histories in basaltic systems, the assumption of closed system degassing is also reasonable because the growing vapor bubbles are unlikely to segregate from the ascending melt unless they reach a critical radius (Vergniolle & Jaupart, 1986).• Non-isothermal decompression and magma stalling upon ascent.Individual crystal populations may undergo different ascent histories and thermal conditions.Magmas may stall in the crust, and undergo partial or even full re-equilibration of water in the crystal and melt inclusion record.These stalling periods may be associated with cooling and crystallisation and can allow sufficient time for bubbles to segregate from the melt.This can exert changing boundary conditions on the crystal, and can also affect the rate of diffusive exchange.Furthermore, complex flow patterns within the conduit itself may cause magma to sink which may even result in rehydration (Wei et al., 2022).

Morphology of Natural Crystals
For our 3D models, we have assumed that the melt inclusions are spherical and that the host olivine has an idealized morphology based on the forms of Welsch et al. (2014).Even though these approximations are better than previous modeling efforts ( Le Voyer et al., 2014;Newcombe, Plank, Barth, et al., 2020), natural melt inclusions and olivine crystals can take different morphologies, which will depend on the mechanisms of entrapment and crystal growth (Wallace et al., 2021).
• Olivine crystals can take a range of forms from blocky, to skeletal and hopper crystals, which depends on the degree of undercooling and growth rates (Mourey & Shea, 2019).• Melt inclusions can also take ellipsoidal forms, which could depend on the minimization of surface tension (Wallace et al., 2021).• Our models also assume that crystal and melt inclusion morphology remain constant during ascent.This may not necessarily be the case due to crystal growth, resorption, attrition and fracture.

Recommendations
Despite the complexities associated with modeling diffusion of H + in natural systems, simplified models can still provide powerful insights into the dynamics of magma ascent, particularly if large crystal and melt inclusion populations are considered.Our modeling approach is relatively flexible and can be adapted once we gain more information about natural systems.Here we provide some recommendations based on our modeling insights that can help the community develop more comprehensive modeling strategies.
• The anisotropic analytical solution we provide in Equation 15 can provide a fairly accurate solution relative to the 3D model provided its conditions are met.For instance, the melt inclusion needs to be small relative to the host and in the center of the crystal.This method is the most computationally efficient and could be applied to large populations of measured melt inclusions quickly.• One dimensional numerical models also provide a computationally efficient solution; however, these are the least accurate when compared to the 3D model.Equation 14 can be used to correct 1D decompression rates, and is not limited by the location of the melt inclusion along [100], or its size.Corrected 1D decompression from melt inclusions close to crystal edges could therefore be used to complement decompression rate estimates using the analytical solution.• Two dimensional models may offer a compromise for providing accurate solutions at a wider range of crystal shape and melt inclusion configurations.This approach can also account for multiple melt inclusions if they are in the same 2D plane.However, 2D inversions will be more computationally intensive.Two dimensional modeling on the (010) plane may offer the best compromise for timescale estimation for Fe-Mg interdiffusion and H + diffusion as it contains the fast direction for both of these systems.• When it comes to constraining initial water contents, melt inclusions that are surrounded by other inclusions (particularly large inclusions along [100]) may be the best to target, as the melt inclusion shielding effect may limit water loss from these inclusions.

Conclusions
Measuring and modeling water loss from olivine hosted melt inclusions play a crucial role in understanding the evolution of magmatic systems from the mantle to the surface.This includes estimating the water content of primary mantle melts, estimating magma storage depths in combination with other volatile species (e.g., CO 2 ), and understanding and estimating timescales of processes associated with final magma ascent.Most attempts at modeling water loss from olivine-hosted melt inclusions apply simplifying assumptions associated with crystal and melt inclusion geometrical configurations, such as spherical or 1D models.
We applied 2D inversion models to measured water profiles in olivine crystals from the 1977 eruption of Seguam.In this model, the covariance between D 0 and E a was considered in the error propagation.This approach generates median decompression rates of 0.07-0.36± 0.04 MPa s 1 , which is consistent with previous estimates, and likely represents natural variability in final magma ascent rates.
Guided by the results of these 2D inversions, we have developed a new finite element 3D multiphase diffusion model that accounts for water loss from both olivine and melt inclusions.This model uses an idealized olivine crystal morphology and can include multiple melt inclusions.We used this model to assess the accuracy of simplified model configurations including 1D models with symmetric and asymmetric crystal geometries around the melt inclusion, 2D models sectioned through (010) and (001), and an anisotropic analytical solution that can be applied in an infinite domain.We find that the 1D, 2D and analytical models generally underestimate the amount of water loss from melt inclusions compared to the 3D models.These models also underestimate decompression rates compared to the 3D models.This is because additional water is lost along the additional dimensions that are not accounted for, or in the case of the analytical solution, the melt inclusion geometry falls outside the limits in which the solution is valid.The amount in which these models underestimate water loss and decompression rate depend on the position of the melt inclusion along the [100] axes, inclusion size and the sectional cut of the crystal.If the melt inclusion is located close to the crystal edge along [100], water loss and decompression rate predictions from 1D models closely match those from the 3D models, because overestimations in diffusive flux out of the melt inclusion counter balance the extra water loss from the unaccounted dimensions.We developed an empirical correction that aligns predictions from 1D models closer to those from the 3D models.If the inclusion is small relative to the host and is situated away from the edge of the crystal, the anisotropic analytical solution that we present provides an accurate and computationally efficient way of determining water loss and decompression rate.
We also developed 3D models that include multiple melt inclusions in order to replicate some crystals found in nature.We find that melt inclusions shield each other and reduce the amount of water loss when compared to crystals that only contain a single melt inclusion.This shielding effect depends on decompression rate, the size of the melt inclusions, and is enhanced if the melt inclusions are situated along the [100] direction.Crystals with dense populations of melt inclusions may therefore provide a more reliable estimate of magmatic water contents prior to ascent driven degassing.
Our modeling approach has shown that factors such 3D crystal geometry and the configuration of melt inclusion populations can play an important role in constraining accurate decompression rates, and for recovering water contents in natural magmatic systems.

Figure 1 .
Figure 1.Transmitted light (a, c, e) and reflected light (b, d, f) photomicrographs of polished Seguam olivine crystals used in this study.(a, b) Shows sample Seg1-MI1, (c, d) shows Seg4-MI1, and (e, f) shows Seg13-MI1.The location of measured secondary ion mass spectrometry profiles adjacent to the exposed melt inclusion are shown in the reflected light images.Profiles measured close to the [100] crystallographic direction are marked in red, whilst profiles measured close to [001] are shown in blue.Scale bars are shown next to the reflected light images.

Figure 2 .
Figure 2. Measured water concentrations, major element contents and best fit Bayesian inversion models for Seg13-MI1 (a), Seg1-MI1 (b), and Seg4-MI1 (c).The first column shows color maps and contours of 2D best-fit models for each crystal.The location of secondary ion mass spectrometry (SIMS) spots are shown in red.The second and third columns show SIMS water profile data measured along [100] and [001] respectively, with the best fit curves from the 2D model shown in red.Error bars on SIMS data are shown and may be smaller than data symbols.Insets show forsterite profiles (X Fo , mol fraction) measured in the same locations.Melt inclusions are represented by gray regions.

Figure 3 .
Figure 3. Density plots showing the posterior distributions of the 2D Bayesian inversion calculations performed on sample Seg13-MI1.Inverted parameters include log decompression rate (log 10 dp/dt), the anisotropy of H + diffusion in olivine (D [100] /D [001] ), the olivine-melt partition coefficient for water (K), initial water content (H 2 O i ), temperature (T ), log 10 D 0 and the activation energy (E a ) for the diffusion of H + in olivine.Prominent trade-offs between decompression rate, temperature, and the partition coefficient can be observed.

Figure 4 .
Figure 4. Stacked histograms showing the posterior distributions from the 2D Bayesian inversion calculations for Seguam olivine samples Seg1-MI1 (dark blue), Seg13-MI1 (light blue) and Seg4-MI1 (slate blue).Inverted parameters shown include: (a), the diffusive anisotropy of H + in olivine (D [100] /D [001] ); (b), log decompression rate (log 10 dp/dt); (c), the partition coefficient (K); (d) temperature (T ); (e), log 10 D 0 for H + diffusion in olivine; (f), the activation energy (E a ) for H + diffusion in olivine; and (g), the initial water content (H 2 O i ).The red curve in panel (a) is a log normal distribution fitted to the total anisotropy distributions for Seg1-MI1 and Seg13-MI1.The red dashed red line show the modal value selected for further modeling.The red dashed lines shown in panel (b) show the upper and lower bounds of decompression rates estimated by Newcombe, Plank, Barth, et al. (2020).The red dashed line in panel (c) shows the median K estimated from the distributions of Seg1-MI1 and Seg13-MI1, which was then used in further modeling.The red curves in panels (d-f) show the Gaussian prior distributions with mean values marked by the dashed lines.These mean values were used in later models.

Figure 5 .
Figure 5. Distributions of the physical and geometrical parameters that were varied during the 3D Monte Carlo modeling for all of the partition coefficients that we used (K = 0.000459, 0.001, and a random subset between 0.0004 and 0.002).(a) Shows the log of magma decompression rate (log 10 dp/dt MPa s 1 ).(b) Shows crystal size which has been parameterized as the length of the [001] axes from the center of the crystal (μm).(c) Shows the melt inclusion size which has been parameterized as melt inclusion radius (μm).(d-f) Show the location of the melt inclusion along the [100], [010], and [001] directions respectively (MI).These position values have been normalized to the length of the corresponding axes (L).Values of MI/L close to 0 are close to the center, whilst values close to 1 are close to the crystal edge.(g)Shows the distribution of water solubility curves that correspond to different initial water contents in the melt and starting pressures.These curves represent the equilibrium degassing pathways that were used in each model realization, and thus relate to the exterior boundary condition.Solubility curves were calculated using VolatileCalc(Newman & Lowenstern, 2002).

Figure 6 .
Figure 6.Example model configurations used in the Monte Carlo simulations.(a, f) Show the initial 3D models, which have been sectioned along (010) and (001) through the center of the melt inclusion.(b, g) Show 2D models based on (010) sections through the center of the melt inclusion.(c, h) Show 2D models based on (001) sections through the center of the melt inclusion.(d, i) Show 1D models along the [100] axis that incorporate both edges of the crystal (termed asymmetric models here).(e, j) Show 1D models with reflective boundary at the center of melt inclusion and only incorporates the shortest direction along the [100] axis (termed symmetric models).Melt inclusions are represented by gray regions.

Figure 7 .
Figure 7. Results of the Monte Carlo models showing melt inclusion water loss.Comparisons are made between the 3D models and the 1D asymmetric models (a), 1D symmetric models (b), 2D models sectioned along (010) (c), 2D models sectioned along (001) (d), and the anisotropic analytical solution presented in Equation 15 (e).Each column shows model results using different olivine diffusive anisotropies (D [100] /D [001] ).1:1 lines are shown in black with ±20% envelopes shown with dashed lines.Points have been color-coded based on log 10 decompression rate of the 3D models.

Figure 8 .
Figure 8.Comparison of 3D model decompression rates to inverted decompression rates obtained from 1D asymmetric models (a), 1D symmetric models (b), 2D models sectioned along (010) (c), 2D models sectioned along (001) (d), and the anisotropic analytical solution presented in Equation 15 (e).Each column shows model results using different olivine diffusive anisotropies (D [100] /D [001] ).1:1 lines are shown in black with dashed lines marking points within a factor of 5. Points have been colorcoded based on melt inclusion radius.

Figure 9 .
Figure 9. Geometrical controls on modeled decompression rates.Each panel shows how the position of the melt inclusion along the [100] axis affects the ratio of simplified model to 3D modeled inversion decompression rates (dp/dt R model /3D ).For 1D numerical models, R model/3D is the same as R 1D/3D in the main text.For the analytical solution, R model/3D is the same as R A/3D in the main text.MI [100] is the distance from the nearest crystal edge to the edge of the melt inclusion along the [100] axis.L [100] is the distance between the crystal edge and the center of the crystal along [100].The melt inclusion is close to the edge when MI [100] /L [100] is close to 0, and is close to the center of the crystal when MI [100] /L [100] is close to 1. Point have been colored based on the size of the melt inclusion radius (r MI ) relative to the length of the [100] axis (L [100] ).The black line shows where 3D and 1D model inversions are equal (dp/dt R model/3D = 1).The red dashed lines show empirical model fits through the 1D model data.(a) Shows results for 1D asymmetric models, (b) shows results for 1D symmetric models, and (c) shows results from the anisotropic analytical solution.Each column shows model results using different olivine diffusive anisotropies (D [100] /D [001] ).

Figure 10 .
Figure 10.Configurations of 3D models that contain a single central melt inclusion.2D sections through the center of the crystal are shown perpendicular to the main crystallographic directions.Models with 25, 50, and 75 μm radius central inclusion are shown in panels (a), (b), and (c).All melt inclusions are represented by gray regions.Results of models with a decompression rate 0.05 MPa s 1 and a diffusive anisotropy of 15.6 are shown.

Figure 11 .
Figure 11.Configurations of 3D models that contain multiple melt inclusions.2D sections through the center of the crystal are shown perpendicular to the main crystallographic directions.Models with multiple melt inclusions with a central melt inclusion with 25 μm, 50 and 75 μm radius are shown in panels (a), (b), and (c).The surrounding melt inclusions in these models have 50 μm radius.All melt inclusions are represented by gray regions.Results of models with a decompression rate of 0.05 MPa s 1 and a diffusive anisotropy of 15.6 are shown.

Figure 12 .
Figure 12.Model configurations of 3D models in which there are 2 melt inclusions: a 25 μm radius central inclusion and a 75 μm radius inclusion situated along one of the crystallographic axes.The large melt inclusion was placed along the [100] axis (a), the [010] axis (b), and the [001] axis (c) 250 μm from the central inclusion.All melt inclusions are represented by gray regions.Each column shows a 2D section through a 3D model.Results of models with a decompression rate of 0.05 MPa s 1 and a diffusive anisotropy of 15.6 are shown.

Figure 13 .
Figure 13.The role of multiple melt inclusions in controlling water loss during magma ascent.(a) Comparison of water loss from a central melt inclusion in models with a single melt inclusion and models in which a central melt inclusion is surrounded by melt inclusions with 50 μm radius.See Figures 10 and 11 for configurations.Each point shows the difference in water loss between the single melt inclusion models and multiple melt inclusion models at different decompression rates.The lines are skewed Gaussian peak fits to the data.The radius of the central melt inclusion was varied (r = 25, 50, and 75 μm), which is shown by the different shades of blue.Different panels show different model anisotropies (3, 15.6, 35) (b) Comparison of water loss from single melt inclusion models (radius of 25 μm) with models in which a 75 μm radius melt inclusion has been positioned along one of the principle crystallographic axes ([100], [010], and [001]), which corresponds to the different shades of blue.See Figure 12 for configurations.Different panels show different model anisotropies (3, 15.6, 35).

Table 1
Parameters Used in the 3D Monte Carlo Modeling With Their Modeled Ranges and Distribution Types Towbin et al. (2023)al value (0.0009 ± 0.0003) fromTowbin et al. (2023).c Random subset spanning full range of natural and experimental measurements.Uniform Size of the crystal based on the length of the [001] axis from the center of the crystal.

Table 2
Coefficient Values for the Empirical Correction Shown in Equation14