Linking fluid dynamics and olivine crystal scale zoning during simulated magma intrusion

The compositional zoning styles of natural crystals produced during magma intrusion can be used to investigate the structure of magmatic plumbing systems and its relation to expressions of volcanic unrest (seismic, deformation, volatiles). However, magma intrusion is a progressive, dynamic process and yields non-monotonic heterogeneities in physio-chemical variables such as complex spatial variations in temperature and liquid composition with time. Such changes in variables are difficult to incorporate in models of crystal zoning in natural systems. Here we take another approach by integrating the results of a numerical multiphase simulation of melt arrival in an olivine-rich reservoir with models of chemical re-equilibration of olivine. We evaluate the diversity of chemical zoning styles and the inferred time scales using Fe–Mg diffusion in olivine for a limited range of system geometries and time-composition-temperature values. Although our models are still a large simplification of the processes that may occur in natural systems we find several time-dependent and systematic relations between variables that can be used to better interpret natural data. The proportions of zoned and unzoned crystals, the zoning length scales, and the calculated diffusion times from the crystals, vary with time and the initial position of the crystal in the reservoir. These relationships can be used, for example, to better constrain the plumbing structure and dynamics of mafic eruptions from monogenetic volcanoes by detailed studies of changes in the zoning of the crystal cargo with eruptive sequence. Moreover, we also find that the time scales obtained from modeling of crystals at a single temperature and boundary condition tend to be shorter (> about 25%) than the residence time, which could also be tested in natural studies by combining crystal time scale records with monitoring datasets.

2018). However, magmatic systems and the crystal cargo are heterogeneous at a variety of scales, and magma additions into pre-existing reservoirs likely involves a wide range of physical and chemical changes that are difficult to uniquely account for in models of crystal zoning (Bergantz et al. 2015;Cheng et al. 2017;Costa et al. 2010;Macdonald et al. 2008;Nakamura 1995;Probst et al. 2018;Wallace and Bergantz 2005).Although it is possible to construct relatively sophisticated models of compositional changes at the boundary of crystal and temperature with time (e.g. Costa et al. 2008;de Maisonneuve et al. 2016), it is unclear how to extend that to the complexities of natural events.
Modeling of crystal-fluid interactions is not straightforward, as the presence of crystals greatly affects the mechanics of a magma even at crystal concentrations of as little as 25 vol. %, as the system can form crystal-contact networks (Philpotts et al. 1999), yielding fragile networks by the transmission of stress by crystal-crystal contacts called force chains Sun et al. 2010;Sandnes et al. 2011;Estep and Dufek 2012). Thus, simple suspension rheology is inadequate to resolve crystal transport and reaction at the crystal scale. Investigating the mechanics involved in an open-system intrusion requires an approach that resolves both the micro (crystal-scale) physical interactions and the meso-scale circulation. Bergantz et al. (2015) used a discrete element method-computational fluid dynamics approach that includes crystal-scale frictional, collisional, translational, and buoyant forces, as well as viscous particle-fluid coupling under basaltic conditions. The particle-based simulation records the physical interactions and circulation of all crystals and creates a times series for each individual crystal as a record of changes in temperature and composition. Linking the results of such fluid-dynamic simulations with the processes of chemical re-equilibration within the crystal allow us to test the effects of changing boundary and temperature conditions into the time scales that are obtained from natural crystals ).
An obstacle in the application of crystal zoning studies to magmatic systems is that sample sizes are very small relative to total crystal contents, and sampling itself can be very incomplete. There is a need for an interpretive framework that will move the sparse crystal chemical data to the mesoto-macro scale. Hence the goal of this work is to be a proofof-concept for an approach that combines an intentionally simplified simulation of an open-system event with measures of zoning that are widely applied to complex natural systems. The objective is not to produce a highly resolved simulation of a specific natural system, but instead to explore how the most important degrees of freedom present in a generic, simple system may produce and inform the observed diversity in a naturally occurring crystal cargo. A related objective is to illuminate the statistics of establishing ergodicity in any study that employs crystal zoning. Each magmatic system is the result of a unique process of construction and open-system events, differing intensive variables and local kinetics. Hence we focused on, and demonstrated various measures in, the crystal cargo related to the system-wide behavior and global concepts at the meso-to-macro scale in the crystal cargo, which can't be otherwise attainable under current sampling protocols. Our intention is not to provide a detailed description of a particular geological episode but instead demonstrate various statistical approaches that are needful in the development of an interpretive framework.
Here we focus on the interaction and reequilibration of olivine in basaltic melts and Illustrate should be illustrate the relationships between the time scales of the different processes, the effects on crystal zoning patterns, and show how these can be recognized and accounted for by large numbers of crystals to create an ergodic population and down-sampling strategy. We also discuss the potential use of our results for the interpretation of natural crystals in basaltic eruption sequences and the monitoring signals that precede them.

Methods: fluid dynamic simulations and crystal scale modeling
Numerical simulations of magma intrusion in reservoir collaborators (Bergantz et al. 2015, 2017;Carrara et al. 2019;Schleicher and Bergantz 2017;Schleicher et al. 2016) have used a discrete element method to perform computational fluid dynamic models of a crystal-free basaltic liquid intruding at the base of a crystal-rich mush of olivine. In their simulation magma is intruded from below and through a crystal-rich pile. They found that most magma interactions occurred within the so-called "mixing bowl", which is a wedge-shaped area of the crystal mush right above where the intrusion occurs. The geometry of the mixing bowl is set by the failure along initially frictional contacts of the settled crystals, producing what is known as melt-present faults, widely observed in the plutonic environment (e.g. Bergantz et al. 2015). As most of the circulation is within the mixing bowl, this is where thermo-chemical effects of melt on the crystals will be maximal. The simulations incorporate a range of realistic parameters for the fluid and crystal dynamics (details in Appendix of Bergantz et al. (2015;) and track the position and a parametrized boundary temperature of hundreds of thousands of crystals. A first order finding of these simulations is that crystals may experience a wide range of environments over time and thus a wide range of crystal zoning patterns can be produced even during a single open-system event.
Here we use the results from Schleicher and Bergantz (2017) to evaluate the importance of the changes in environment that individual crystals may experience during a situation that tries to simulate magma mixing/mingling. We use these results to make some inferences and relate them to studies of natural crystals. The numerical simulations of Schleicher and Bergantz (2017) used the following range of physical and chemical parameters: a domain size of about 2.6 × 1.3 m, an initial thickness of the crystal pile of 0.8 m (= H o ), and a olivine crystal diameter of 4 mm, there was only one crystal size in these simulations. The Reynolds is the magma input velocity, d is the length scale of the entrance of magma input, ρ f and μ f , are the fluid density and viscosity, respectively)] varies between about 0.1 and 1. The composition and temperature of the melt was parameterized with a scalar whose value was equal to zero for the resident melt, and one for the intruding melt, with values in between reflecting mixing proportions between resident and new intruding melt. The time (t) of the simulations was 100 s but was normalized with t* = (U o /H o ) × t, where U o is as above, and H o is the vertical thickness of the crystal mush. The simulations were run for t* values of 0-2.75.
To combine the results of these simulations with the olivine zoning patterns modeled here time needs to be rescaled: the 100 s that were simulated are not representative of the times that are typically recorded from olivine zoning in basaltic magmas, which range from a few days to a few months (e.g. Albert et al. 2016;Lynn et al. 2017;Ruth et al. 2018;Pankhurst et al. 2018). Thus, we chose a time of 0-100 days as the residence time ( Fig. 1) that is representative of those from natural crystals, and we focus on the changes of the system at 10, 40, 70 and 100 days. We note that in the fluid dynamic simulations the intrusion begins at t = 0 and continues for 60 days, and then stops. This means that the duration of the intrusion is relatively long compared to the longest time residence time of 100 days. To conserve the t* of the simulations requires increasing the thickness of the initial crystal pile by about × 1000 (so instead of about 0.8 m it would be 800 m), and decreasing the input magma velocity by × 100. Such changes would increase the Reynolds number by × 10, but this can be compensated by an increase of the width of the incoming magma by the same amount (see Online supplementary Appendix 1 and Table 1s). Thus, although our total times are quite different from those of the fluid dynamics simulations, the scaling of the system to larger sizes still preserves the main relations of the simulations and thus are still a good proxy for the dynamics during melt and crystal interactions.
For the models of olivine re-equilibration, we tracked the changes of the crystal position and local scalar value at the boundary of all the 23,822 crystals from the "mixing bowl". Thus, all the active crystals are sampled so it is a perfectly ergodic sample of the process. We used the scalar value as a proxy for the temperature (T) that the crystal experienced over time. The initial condition is a scalar value of zero corresponding to T = 1150 °C and is uniform throughout the crystal pile, and a value of one corresponds to T = 1200 °C, which is that of the intruding magma. Mixing of resident and new melt was quantified by a linear relation between the scalar and temperature values. We used a basaltic andesite composition and MELTS (Ghiorso and Sack 1995) to derive a relation between the forsterite content (Fo as chemical indicator of the Fe/Mg) and the temperature (detailed in Online supplementary Appendix 1; the same as used in Costa et al. (2008)): Fig. 1 Illustration of the different time scales we employed in this study. Delay time is the time interval between when the start of the intrusion and when the crystal first records it as a change at the boundary. This interval varies from crystal to crystal as they stay in different parts of magma reservoir. We call forward time that which is obtained from diffusion models that use the boundary and temperature evolution of the fluid-dynamic simulations. We also calculate an inverse time from the crystal zoning via diffusion models as if they were natural crystals. Finally, the total time from intrusion start to eruption is the residence time where, T ranges from 1150 to 1200 °C, at an oxygen fugacity buffered at the N-NiO reaction and 1 bar. The resident magma was thus made of olivine crystals that were unzoned and of Fo = 0.74 throughout. The intruding magma had the same basaltic composition but was at 1200 °C and is crystal free, but would be in equilibrium with about Fo = 0.82.

Modeling of olivine reaction progress and chemical zoning
For the re-equilibration of olivine we considered that the compositions of the crystals were modified only by diffusion from the rim towards the interior, and thus the driving forces for the chemical zoning were the changes of melt and temperature at the boundary with time. Schleicher and Bergantz (2017) used a parametrized solution for the convective crystal dissolution (Zhang 2008) and obtained very high dissolution rates that would basically dissolve all the crystal in the mush within several hours. However, the parametrized model is valid for an infinite melt with a constant composition, and thus yields higher dissolution rates than for the case of a crystal mush with a liquid composition that is changing with time due to dissolution, growth, and diffusion. Several analytical solutions exist for dissolution and crystal growth on the convective regime (e.g. Zhang 2008), but these are also for semi-infinite liquid condition. Currently there is no widely accepted model or parameterizations to properly couple the effects of magma dynamics with crystal growth and dissolution. Thus, in addition to a simplified magma chamber our model calculations also simplify the details of the crystal-melt boundary interactions, as crystal growth, dissolution, and diffusion are occurring concurrently, and these processes will have an effect on the temperature, liquid compositions, and thus the details of the crystal and liquid dynamics. However, studies of natural olivine crystals can typically distinguish between the effect of crystal growth and diffusion, and quite a few of them have shown that the Fe-Mg zoning of the olivine can be better explained by diffusion (using arguments of diffusion anisotropy and modeling of multiple elements; e.g. Costa and Dungan (2005); Costa et al. (2008); Kahl et al. (2015); Shea et al. (2015)) than by crystal growth. Despite the limited range of parameter values that we have been able to explore, and the simplification of the processes at the crystal-melt boundary, our approach is more sophisticated than previous models as it captures the complexities of changing boundary conditions with time and the effects it may have in crystal zoning that is dominated by diffusion.
For modeling the zoning patterns, we used Fick's second law of diffusion and the diffusion coefficient of Fo (D Fo ) including compositional dependence and the relation from (Dohmen and Chakraborty, 2007a;Dohmen and Chakraborty, 2007b) parallel to c axis (in µm 2 s −1 ): With fO 2 at the Ni-NiO buffer, D o = 6.17 × 10 -10 , E = 201 kJmol −1 , T in o C. We performed two types of modeling: (1) diffusion in the crystal using the temperature and boundary composition changes with time as indicated by the fluid dynamic simulations, which give us the "forward time" (Fig. 1), and (2) diffusion models of the same profiles as one would do if they were natural crystals, which we call "inverse time" (Fig. 1). In these models the temperature and boundary evolution with time are unknown, and we used two approaches as a proxy for the temperature: (2a) using a mean temperature of 1175 °C for all crystals, and (2b) using the maximum Fo measured as a proxy for temperature so that each crystal may be modelled at a different temperature according to its profile. For the initial profiles in both cases we either used a homogenous composition throughout the profile for simple reverse profiles, or step-wise changes in Fo for complexly zoned crystals. These are standard conditions used in many studies (e.g. Costa et al. 2008). In addition to the forward and inverse times, we use the term "residence time" to indicate the time since the beginning of melt intrusion and the quenching of the system (e.g. due to eruption), and the "delay time" is the time since the beginning of intrusion and the moment that the crystals first record a change at their rim ( Fig. 1). In the models we tracked the compositional changes and calculated the diffusion profiles and times for 23,822 crystals. The models were run in the supercomputer cluster (Komodo) of Nanyang Technological University with a total of 70 nodes installed with Centos 6.5 OS. The nodes have 24 cores intel processors (Intel(R) Xeon(R) CPU E5-2609 0 @ 2.40 GHz) and 128 GB memory each.

Results
Our model of the coupling of the fluid dynamic simulations with diffusion in the crystal emphasizes two main controls on the nature of the crystal cargo. One is that there is a time delay between the beginning that new liquid enters the reservoir and the time the crystal rims record it, and this delay changes with time and space, and should also be dependent on the duration of the intrusion. The other is that crystal rims can experience multiple compositional changes at their margins over time, and thus • 10 3(0.9−Fo) • exp −E 8.31 × 10 −3 × (T + 273) • 10 12 the observed zoning after quenching does not necessarily reflect the complex history of the crystal. While both of these are to be expected, here we show how their impacts can be independently quantified. To link the simulations and diffusion modeling to applications for natural crystals, we calculate the forward and inverse times from the crystal zoning, and we evaluate their different times using two different temperature proxies and zoning styles. The different topics of the results are presented in sequence below.

Delay between the beginning of intrusion and when the crystals first start recording it
One explicit result of the numerical simulations of magma intrusion by (Bergantz et al. 2015;Schleicher and Bergantz. 2017) is that different crystals encounter the new melt (as changes in the boundary composition) at different times. This means that there is a delay between the time of the beginning of the intrusion and the moment each crystal records it with a change at the boundary which we call delay time. We tracked the position and changes at the boundary of the 23,822 crystals from different positions in the 'mixing bowl' and calculated the time delay since the beginning of the intrusion, and how the response with the residence time (Fig. 2). We find that the magnitude of the delay depends on the residence time. For example, after a residence time of 10 days, most crystals (about 80%) have not yet recorded the intrusion (Fig. 2), meaning that if we were to quench the system (e.g. during an eruption) 80% crystals would show no zoning. For a longer residence time of 40 days, many more crystals have started to record the intrusion, but there still is about 15% of them which have not (Fig. 2). After more than about 70 days of residence all crystals have recorded at least one change at a boundary (Fig. 2). However, the delay time varies depending on the original location of the crystal when the intrusion first started: crystals proximal to the new intrusion (from the bottom of the mixing bowl in this case) show less delay than those at further away (at the top of the mixing bowl); this can be expected as the simulations we used are for a melt intrusion from the bottom of the crystal pile. Moreover, the delay probably also depends on the duration of the intrusion. Such delays in changes at the boundary imply that the crystals would only start diffusing significantly after the first intrusion event, and thus they will record shorter times than the real intrusion time. It also implies that quenching of the system at a given time could yield a significant number of unzoned crystals. These findings are of general value and applies for different configurations of melt intrusion and to models that incorporate growth or dissolution, because the crystals will only start reacting once their environment has been perturbed.

Olivine zoning patterns
The coupled modeling of fluid dynamic simulations and diffusion of Fo in olivine crystals produces a significant variety of zoning patterns (Fig. 3). Most crystals have simple reverse zoning patterns, with a higher and variable Fo at their rims (Fo 75 to about Fo 82 ) than at their interiors. Others are more complexly zoned, with a reverse intermediate zone of higher Fo (Fo 75 to Fo 77 ) than their interiors, before they are normally zoned with lower Fo at their rims. The crystals that have not recorded the intrusion will not be zoned, and thus have homogenous Fo profiles (so-called unzoned patterns).
As the residence time increases the proportion of zoned crystals increases, from about 22% after 10 days, to about 83% after 40 days (Fig. 3), and the Fo profiles are longer and smoother as crystals have more time to reequilibrate. After about 70 days all crystals are zoned, and about 70-80% of them show simple reverse zoning patterns. Moreover, the compositional variability (mean and standard deviation) of the crystal rims change over time (Fig. 3). This may reflect the progressive mixing of the liquids, although it is difficult to interpret in detail for natural systems because we can't account uniquely for crystal growth, dissolution, or the changes in the liquid composition due to diffusion with the crystals.
Although the simple zoning patterns of some crystals may lead to the inference that the conditions at their margins have been relatively constant over time, this is actually not the case, and they may have undergone several compositional reversals at their rims before the last one that produced the normal zoning (see animation of boundary layer changes over time in Online supplementary Appendix 2). Complexly zoned crystals have also recorded a range of compositions over time and this is more completely expressed in the profiles, although still underestimates the variety of changes they have actually undergone. The types of zoning patterns we simulated are also found in basaltic andesitic rocks from magmas that have undergone mafic replenishment such as the Llaima (Chile) and Etna (Italy) volcanoes (e.g. Kahl et al. 2015;Ruth et al. 2018).
The crystal zoning profiles of Fo we simulated have an associated forward time that is significantly shorter than the residence times (Fig. 4). This difference is mainly because different crystals feel the intrusion at different times. The relative error of the forward times (= 100 × forward time/ of olivine at 70 days shows that all crystals are zoned (d) zoning patterns of olivine at 100 days are similar to those at 70 days, but there appears to be a higher proportion of simple vs complex zoned crystals. We note that these are time integrated profiles but that the crystal may have experienced a wide range of compositions at the boundary over time as shown in the animation in the Online supplementary Appendix 2 [forward time − residence time]) is close to 100% for most crystals (about 80%) for a short residence time of 10 days, but as the residence time increases the error of forward times decreases, with about 50% of crystals having an error of < 30% for a residence time of 100 days. However, the spread of calculated forward times for a given residence time may also be influenced by the total duration of the intrusion.

Forward and inverse models, and proxies for the temperature history to calculate inverse times from crystal zoning
Another aspect that is important to address are the time scales that can be obtained from modeling the crystals if we did not have a priori knowledge of their temperature and boundary compositions changes. This is the common situation for most studies of natural systems, and we call these inverse times ( Fig. 1; see Methods section). We obtained this from modeling the profiles produced by coupling the simulations and diffusion in the crystal shown in the previous sections but with a constant temperature and boundary composition. The temperature history is very important, as diffusion rates and thus calculated times are very temperature dependent (e.g. Costa et al. 2008). For systems that oscillate around a mean, it is possible to use a mean temperature to get meaningful results without needing to incorporate the full complexity (e.g. Costa et al. 2008;Lasaga and Jiang 1995). For systems with a protracted history of cooling or heating such as for metamorphic rocks, the concept of characteristic temperature was successfully demonstrated (Chakraborty 2006;Ganguly 2002).To illustrate the role of temperature we performed inverse diffusion models using two cases: one with the mean temperature of the two end-members (1175 °C), and another where we used the maximum Fo content recorded in the crystal as a proxy for the temperature. These are common approaches to approximate the relevant temperature of the process when studying natural crystals (e.g. Costa et al. 2008).
One consequence from the discussion above is that the forward times are typically shorter than the residence time due to the time delay of the crystals to begin to record the intrusion at their boundary ( Figs. 1 and 2). Such an effect also applies to the calculated inverse times, and comparison of forward and inverse times provides a direct expression of the typical simplification employed when modeling crystals with constant temperature and boundary conditions when they have experienced many. We have found that inverse times can be longer or shorter than the forward times (Fig. 5), and that there is no major difference between the models that use the mean or the maximum Fo as a proxy for temperature (Fig. 5). The relative differences (Fig. 5) between the forward and inverse times show a positively skewed distribution when the residence time is < 40 days, but they become a normal distribution with time (e.g., > 70 days; Fig. 5). A broad conclusion from this analysis is that the inverse and forward times tend to overlap with increasing residence time, and that whether a mean temperature or the maximum Fo content is used as a proxy for temperature does not give very different results.
Another aspect that needs to be quantified is the relationship between the inverse time as could be derived from studies of natural crystals, and the residence time (Fig. 6). This gives an indication of the influence of the delay and the uncertainties on the temperature and composition history on the derived time from modeling the crystal zoning patterns. We find that for short residence times (e.g. 10 days) the inverse times are much shorter than the residence times, but as residence time increases to 70 and 100 days the inverse times are normally distributed, with a mean that is about 25% lower than the residence time (Fig. 6). The exact values of these relationships are difficult to extrapolate to other conditions of time and diffusivity but give an idea of the general trends and dynamics of the system.

Effect of zoning profile type on calculated times
We also explored whether some types of zoning patterns better record the residence times than others. This is an indirect way to identify the effect of changes at the boundary with time, as the complexity of the zoning patterns is a record (although imperfect) of the variety of environments  in panels a, c, e, g). Inverse times were calculated using a mean temperature of 1175 °C (symbols in blue), and the maximum Fo of each crystal as a proxy for temperature (symbols in orange). The times of all diffusion models are in general shorter than the residence time (e.g. symbols are inside the box defined by the residence time in panels a, c, e, g) because of the delay time (e.g. Figs. 1 and  2). Panels b, d, f, and h show the relative error between the inverse and forward times using the two temperature models: I-F Error (%) = 100 × [(inverse time − forward time)/forward time]. See text for discussion the crystals have experienced. As noted above (Fig. 3) we found two types of zoning profiles, those that are simple and reverse in Fo, and those that are complex which show reverse followed by normal Fo zoning. We found that the crystals with complex zoning patterns give "mean" inverse and forward times that are in general closer to the residence time than those with reverse zoning profiles (Fig. 7). In addition, the mean inverse times of crystals with complex profiles, and using the maximum Fo as proxy for temperature give are the closest to the residence times (Fig. 7).

Relationship between crystal zoning types, zoning length and times
The proportions of the various crystal zoning patterns, the diffusion distance, and to some extent the variability of rim compositions are related, and systematically change with residence time (Fig. 8). With increasing residence time, the proportion of zoned crystals increases, and the proportion of simple reversed zoning dominates over the complex zoning. With increasing time, the length scale of the crystal zoning increases, and the variability of the rim compositions first increases and then stabilizes (Fig. 8). The exact values of the crystal proportions, diffusion lengths and other changes with time should vary from system to system, but their interrelationship can be used to better understand the processes and time scales obtained from modeling of olivine in natural conditions. For example, characterization of the zoning patterns of large numbers of crystals including documentation of those that are not zoned, carries information about the time evolution of the system which can be combined with the time scales from modeling each profile at an individual level as is commonly done. These relationships can have a diagnostic role for determining whether a magma has been subject to multiple distinct intrusion events and the residence time.

Estimating the number of crystals required to capture the distributions of calculated times
A significant open question in geochemistry is the extent to which any given sample set is ergodic. That is, is the sample set 'complete' with respect to both the natural endmembers and the distributions between them. In the discussion above we have been using a very large number of crystals to characterize the different times and dynamics of the system, indeed all the available crystals. However, to be able to compare our results to natural systems it is important to consider the effect of the sample size on the distributions of inverse times (Fig. 9), as it is currently not the normal practice to study more than a few tens to a hundred crystals from natural systems. The shapes of the distributions of inverse times vary with residence time Fig. 6 Comparison of inverse times (here shown only those calculated with the temperature inferred from the maximum Fo) and residence times. Because of the delay and the changes of temperature and boundary composition that crystals are experiencing the inverses can be shorter or longer than the residence times (inset). This translates to relative errors of the inverse time (= 100 × inverse time/ [inverse time -residence time]) that are very large for short residence time (e.g. black line for residence time of 10 days) but decrease and become more normally distributed with residence time (e.g. red line for residence time of 100 days) Fig. 7 The effect of the zoning profile style on the various diffusion times and comparison with the residence times (R = 10, 40, 70, 100). The mean inverse time of reverse zoning patterns (circle) and complex zoning patterns (square) with 1σ error are calculated with the mean temperature (1175 °C),and using the maximum Fo observed as a proxy for temperature. The inverse diffusion time obtained from crystals with complex zoning are closer to the residence time than that of simple zoning crystals, with the models using the maximum Fo as proxy for temperature giving the closest to the residence time. For clarity the various times associated to a given residence time (R = 10, 40, 70, 100) have been connected by a line and position within the mixing bowl. We find little difference down-sampling from tens of thousands or thousands of crystals to several hundred crystals (> 500). Depending on the original distribution, about between 30 and 100 crystals are necessary to be able recover the original distribution (Fig. 9). Although these are relatively high numbers of data to be obtained and profiles to model, they are not unrealistic, especially for the Fo zoning in olivine which can be calibrated with greyscale of electron backscattered images and a few quantitative analyses per Diffusion distance (i.e. length of the profile with a chemical gradient in micrometers 2 ) of all zoned crystals increases with the residence time. c Standard deviation of Fo composition at the rim of 1000 crystals also varies with time. Although the exact values of the changes should vary with the details of each simulation scenario the observation that they are related could be used for internal consistency of data obtained from natural systems Fig. 9 Examples of the effect of the number of profiles used on the obtained distributions for a residence time of 100 days. Panels a to d show the effect of down sampling on the inverse times starting from all the crystals of the mixing bowl (23,822) to only 30. Because of the distribution of times is very wide with a very large standard deviation, the effect of decreasing the number of crystals to < a few hundred is very important and it is not possible to properly recover original distribution. Panels e to h show the distribution the inverse times of crystals from only the 25% upper part of the system (schematically noted in panel e. These have a narrower and better defined normal distribution and thus down-sampling to much smaller numbers of crystals (e.g. 30) still gives distributions that are close to the original one section (e.g. Pankhurst et al. 2018). Similar numbers of crystals have been proposed to be necessary for considering of multidimensional effects on one dimensional profiles and random sampling of two dimensional sections to retrieve accurate time scale information (Shea et al. 2015).

Discussion
Our aim is to exemplify the controls on, and the complexity of, the crystal cargo in a highly simplified basalt reservoir containing a crystal mush, undergoing an open-system event. We did this by employing a multiphase numerical model of discrete element-fluid dynamics where the physical interactions are resolved at the crystal-scale. The temperature timeseries of every crystal was then subject to a diffusion model as a post-processing step. Hence the reaction progress itself was not incorporated in the physical model. This produces negligible error in the dynamics as the circulation and corresponding crystal response will not be much affected by small changes in liquid density or changes in crystal size or density. Our goals are to demonstrate the origins of fundamental features commonly observed in the crystal record, which can be expected in a wide variety of magmatic systems. This is a step-forward in illuminating how distinct populations arise, are dispersed, and the challenges in their interpretation.
Our approach of modeling the crystal response as diffusion is that it is more applicable to natural systems for which it can be shown that most of the zoning in the olivine is due diffusion, which is the case for quite a few studies (e.g. Albert et al. 2015;Kahl et al. 2017;Ruth et al. 2018;Saunders et al. 2012). From the times we calculated, the inverse times are those that could apply to natural datasets where olivine zoning is the primary record of magma intrusion. Below we discuss how our results can be applied to better interpret the olivine record of eruptions sequences, and the relation to other observations that could be used to gather complementary information on volcanic processes and the residence time.

Effect of the eruption on sampling of the plumbing system and crystal diffusion times
The fluid dynamic simulations and crystal diffusion modeling show that even a very simple magma intrusion is likely to produce a heterogenous record in the crystals, with those at the closest to the newly intruding magma capturing more of the process (e.g. less delay times, Fig. 2). The proportions of the zoned crystals and the distributions of the calculated inverse times are also significantly different depending on the location of a crystal with regards to the intrusion (Fig. 10). The calculated inverse times for the different parts of the reservoir are still between 25% (for crystals from the bottom part of the crystal pile) and > 50% (crystals from the top of the crystal pile) shorter than the residence time, but the systematic behavior of the crystal data as a function of location in the reservoir has the potential to be used to reveal the dynamics and configuration of the plumbing system. Although generic by nature, our findings might be especially applicable to short -lived and small volume mafic eruptions such as those from monogenetic volcanoes or dike fed events (e.g. Valentine and Gregg., 2008). Documenting of the changes in bulk-rock and crystal chemistry with eruption sequence (e.g. Brenna et al. 2010 andMcGee et al. 2019) may reveal whether the eruption sequentially tapped the reservoir from top to bottom, whether magmas from a range of depths of the plumbing system where Red curves are crystals from the 25% upper part, crystals from the two middle areas are not shown to avoid cluttering of the figures, black curves are crystals from the 25% lowermost part, and grey curves are for all the crystals of the mixing bowl. The % next to each layer correspond to the abundance of zoned crystals. Note that there is a systematic relationship between the crystal position and how inverse time change with residence time. Crystals from the bottom of the mixing bowl (black curve) give inverse times closer to residence time than those from the top (red curve). Such relationships could be used when studying natural systems that may have sampled the complete reservoir in an eruption sequence, and for which it can be shown that the crystal zoning is due to diffusion. Although the absolute values of the time difference from crystals of the different part of the mixing bowl depend on various parameters, the relative differences are likely still applicable sampled at the same time, or whether multiple magmas from different sources intruded during the eruption (e.g. Brenna et al. 2011;Walowski et al. 2019). Depending on the duration of the magma intrusion and on the eruption itself one would expect that the crystal cargo would show progressively higher proportions of zoned crystals and of longer calculated inverse times.

Implications for comparison of crystal time scales with time series of monitoring and observational dataset
One of the objectives of this study was to illustrate and quantify, in the context of a simple example, the heterogeneity in the zoning record introduced during progressive open-system events. This delay decreases with residence time and recognition of this could be useful when trying to relate the intrusion times derived from olivine zoning with the changes in volcano monitoring parameters, such as deformation, seismicity, and gas composition (e.g. Albert et al. 2016;Kahl et al. 2011;Pankhurst et al. 2018;Ruth et al. 2018). The details of our result are however not intended to directly apply to a particular natural system, because the magnitude of the delay will depend on the geometry of the system, on whether crystals are from the most active region (like a mixing bowl), on the proportions of the resident and intruding magma, and on the duration of the intrusion itself. However, detailed comparison of the relationships between the onset and duration of seismic and deformation unrest with those of the inverse times ( Fig. 11) could be used to better understand the relationship between the processes in the plumbing system and those measured at the surface, especially if a sequence of eruption deposits can be sampled.

Conclusions
We combined the numerical results of discrete element-fluid dynamic simulations of magma mixing with the olivine zoning patterns that will be produced by a diffusion-only model. We find that the calculated inverse times from the modeling of crystals at a constant temperature and boundary composition are typically shorter (> about 25%) than the residence times, although the value of the time differences depends on many variables, including the location of the crystal with respect to the intrusion, the total residence residence time, and the zoning style. Depending on the geometry of the system and the intrusion duration, there is a systematic relation between the crystal zoning proportions, lengths, and calculated inverse times, and the location of the crystals on the crystal pile (e.g. bottom vs top). Our approach is still a large simplification of the processes that may occur during magma crystal interactions, but they incorporate part of the complexity of changing temperature and melt composition at the crystal boundary, and explore the effects of down sampling of natural systems. Notwithstanding, our statistical analysis of the processes can be used to better understand the dynamics and plumbing system configuration of, for example, monogenetic volcanoes if crystal data exist for an eruption sequence. Our results can also be used to improve the conceptual models that relate the initiation and duration of magmatic intrusions from monitoring network data (seismic, deformation, gas) and those recorded by the crystal cargo.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Fig. 11
Illustration of the relationship that can be expected between the times scales and signals of volcanic unrest such as deformation and seismicity with those obtained from model of crystals. Deformation and seismicity are adapted from Tilling (2008) and are from changes in tilt and numbers of short period earthquakes of PuuOo volcano in Hawaii one month between January-June 1986. Note that the mean of the crystal times is delayed with regards to the beginning of the magma intrusion as recorded in changes in deformation and seismicity