Theoretical Models of Decompression-Induced Plagioclase Nucleation and Growth in Hydrated Silica-Rich Melts

Magma ascending from the storage region toward the surface may crystallize small and rapidly-grown crystals, designated as microlites. Upon decompression, the rapid changes of the microlite textures, such as number density and crystal size, directly impact the rheology of the magma in the volcanic conduit, and eventually control the effusive vs. explosive style of the subaerial eruption. This is of prime importance for volcanic risk assessment involving highly viscous silica-rich magmas that are prone to fragmentation. To this aim, we present a theoretical modeling of plagioclase nucleation and growth in rhyolitic-like melts, based on previous decompression experiments. Thus, the modeling is valid for plagioclase crystallization in rhyolitic melts decompressed at 875°C from 200 MPa to final pressures of 50, 75, or 100 MPa, which represents effective undercooling (ΔTeff) of 110, 80 and 55°C, respectively. Our results of nucleation-rate calculation using the Classical Nucleation Theory (CNT; case of homogeneous nucleation) and values of crystal-melt interfacial energy (σ) either from literature or empirically calculated, strongly disagree with previously-determined experimental nucleation rates. By inverting the CNT calculation using the experimentally-determined nucleation rates, we propose plagioclase-liquid interfacial energies from 0.041 to 0.059 J.m−2 with ΔTeff increasing from 55 to 110°C, which is about 2–3 times lower than σ determined empirically and macroscopically. We modeled plagioclase growth rates by atom diffusion in melt (following Fick's second law adapted to multicomponent systems), considering a crystal-melt interface advancing over time (as crystal grows) and a chemically-closed finite reservoir (no component supply). The component limiting plagioclase growth has been determined to be CaO, for which we calculated diffusion coefficients from 10−14 to 10−15 m2/s (conditions of a silicic melt, 875°C, and H2O saturation pressures from 50 to 100 MPa). The overall good agreement between the model and the experimental growth laws validates diffusion as the main process controlling isothermal decompression-induced plagioclase growth under moderate ΔTeff < 110°C. Combining both calculations, crystal number density and size, may be relevant to predict plagioclase microlite textures in ascending rhyolitic magmas, which may now be incorporated in conduit flow models to assess the magma rheological behavior that controls eruption dynamics.


INTRODUCTION
The complex eruptive behavior of island-arc silicic volcanoes, i.e., effusive vs. explosive eruptions, is commonly attributed to the magma ascent conditions in the volcanic conduit (e.g., Jaupart and Allègre, 1991;Woods and Koyaguchi, 1994;Barclay et al., 1996). In particular, magma ascent rate controls the time available for syn-eruptive degassing and crystallization, so that a magma that has time and capability to outgas commonly erupts effusively whereas magma that still contains overpressurized gases may give rise to explosive eruptions (e.g., Eichelberger et al., 1986;Melnik et al., 2005;Mueller et al., 2005;Martel and Iacono-Marziano, 2015). One of the major consequences of decreasing melt volatile content through degassing is the increase of the liquidus temperature and subsequent melt crystallization (Tuttle and Bowen, 1958). Crystal nucleation and growth modify the flow properties of magma and gases, and may control the eruptive dynamics by posing the limits of a ductile vs. brittle behavior of the magma in response to the imposed stress. Therefore, investigating decompression-induced crystallization may eventually provide links to eruption dynamics.
Numerous studies have been conducted on crystallization kinetics in silicate melts via cooling (Fenn, 1977;Kirkpatrick, 1981;Lasaga, 1987, 1988;Lattard and Partzsch, 2001;Zieg and Lofgren, 2006;Iezzi et al., 2008;Pupier et al., 2008) and for compositionally simple systems (e.g., Davis et al., 1997;Potapov et al., 1999). These studies have shown that crystallization kinetics exert the primary control on crystal number density, proportion, size, and shape. However, the main driving force for crystallization in ascending hydrous magmas in volcanic context is more decompression (and volatile loss) than cooling, since temperature is roughly constant upon ascent (Wilson et al., 1980). Such decompression-triggered crystals have common sizes below ∼100 µm and are designated as microlites, i.e., small crystals that grow rapidly under thermodynamic disequilibrium.
Several studies have dealt with decompression-induced crystallization in silica-rich melts (dacite to rhyolite), for which the microlite phase is mainly plagioclase, with the aim of understanding the magma ascent conditions in the volcanic conduit and the link with the eruptive behavior. Besides studies of microlite characteristics in natural volcanic pyroclasts (Cashman, 1992;Klug and Cashman, 1994;Wolf and Eichelberger, 1997;Hammer et al., 1999Hammer et al., , 2000Nakada and Motomura, 1999;Clarke et al., 2007;Martel and Poussineau, 2007;Noguchi et al., 2008), experimental simulations of magma decompression have been carried out (Geschwind and Rutherford, 1995;Hammer and Rutherford, 2002;Couch et al., 2003;Martel and Schmidt, 2003;Larsen, 2005;Brugger and Hammer, 2010;Cichy et al., 2011;Martel, 2012;Mollard et al., 2012;Riker et al., 2015;Waters et al., 2015;Befus and Andrews, 2018). On the contrary, theoretical modeling of decompression-induced crystallization is rare. Plagioclase nucleation in a rhyolitic melt has been numerically studied by Hammer (2004). Syn-eruptive crystallization during magma ascent has been considered in numerical models of conduit flow and lava dome extrusion (Melnik and Sparks, 1999;Melnik et al., 2005;De'Michieli Vitturi et al., 2010). These models consider crystallization as a crystal volume fraction at thermodynamic equilibrium, which evolves with time and pressure. One model investigated the effect of crystal growth kinetics on the conduit flow properties via crystal size distributions (Melnik et al., 2011). Theoretical modeling of syneruptive crystal nucleation and growth provides a mean to infer the time evolution and extent of groundmass crystallization for a given volcanic system. Predicting groundmass crystallization is a pre-requisite step to further investigate the physicochemical and flow properties of the magma that control the eruptive behavior.
In this paper, we present two theoretical models. The first one concerns the rates of plagioclase homogeneous nucleation, for which the use of the CNT equation is not new [i.e., similar approach to that used in Hammer (2004)], but allows to refine the plagioclase-liquid surface tensions for a new rhyolitic-like composition. The second model concerns plagioclase crystal growth using the equations of diffusion, which first permits the time-evolution assessment of crystal size (and not only bulk volume) in decompressing magmas; note that the simplified melt composition used here (haplotonalite) is an asset for the diffusion modeling purpose. The proposed crystal growth model is thus suitable for decompression-triggered crystallization from a hydrous rhyolitic melt at constant pressure and temperature.
Rapidly-decompressed magmas that crystallize at constant pressure may not be the first scenario invoked to explain explosive eruptions (which are commonly interpreted as being linked to magma ascent rate). However, two main reasons can be highlighted in favor of the interest of such a scenario. First, such a scenario is possible in nature. Indeed, explosive eruptions generated from lava domes, such as blasts (e.g., the eruptions of Montagne Pelée in 1902 and of Mount Saint Helens in 1980), are interpreted to be triggered by gas overpressures in a shallow-level residing magma. Martel (2012) demonstrated experimentally that the blast magmas ascend from the reservoir fast enough (rapid decompression) to prevent crystallization during ascent. Microlite crystallization actually occurs at shallow/dome level, which drastically increases viscosity and has a major control on the degassing process (i.e., capable of generating a second bubble nucleation event). The second reason for considering such a scenario is that it provides the only way to link textural parameters observable in natural erupted products to thermodynamic parameters controlling crystal growth. Crystal growth is driven by undercooling that is the difference between the crystal liquidus temperature (T L ) and the real temperature of crystallization. In ascending H 2 O-saturated silicic magmas, the driving force for crystallization is not cooling, but melt dehydration as pressure decreases. Anhydrouscrystal T L increases with decreasing melt H 2 O content, so that isothermal decompression of the hydrated magma leads to an increasing difference between T L and magma temperature, referred as to the effective undercooling, T eff (Hammer and Rutherford, 2002). For a given magma temperature, T eff is thus dependent on melt H 2 O content, here controlled by pressure. The scenario of decompression-induced crystallization at constant pressure is thus ideal to properly assign one T eff to one set of crystal textures (e.g., number density, fraction, size, and morphology). This clearly constitutes a prerequisite at thermodynamic equilibrium to further crystallization modeling in which kinetic factors are involved.
After a brief summary of the experimental results of Mollard et al. (2012), from which the modeling aims at reproducing the empirical laws of plagioclase nucleation and growth, we present below (i) new constraints on the plagioclase-melt interfacial energy that controls the nucleation kinetics and (ii) a model of crystal growth-rate limited by element diffusion that may be used to predict final sizes of plagioclase microlites crystallizing from rhyolitic melts.

PREVIOUS RESULTS ON EMPIRICAL LAWS OF PLAGIOCLASE CRYSTALLIZATION
This section briefly presents the experimentally-determined laws of plagioclase crystallization of Mollard et al. (2012), which the modeling aims at reproducing, in order to validate the processes controlling decompression-induced crystallization in rhyolitic melts.
The experimental data set of the decompression-induced plagioclase crystallization and the experimental and analytical methods are reported in Mollard et al. (2012). The starting glass composition is haplotonalitic (HTND; SiO 2 : 78.7; Al 2 O 3 : 14.1; CaO: 1.8; Na 2 O: 5.4, in wt %), which represents a simplified analog composition for natural rhyolitic glasses. Such a simplified composition has the advantage to ease modeling. The HTND glass was saturated with H 2 O at 875 • C and 200 MPa (initial pressure, P i ; 6.7 ±0.3 wt % H 2 O) and decompressed from these conditions using cold-sealed pressure vessels. The final pressures (P f ) ranged from 160 to 25 MPa and decompression rates were of 1,200, 150, and 30 MPa/h (decompression durations of 0.25, 1.5, and 5 h, respectively). For the purpose of this study, only the P f of 100, 75, and 50 MPa will be considered. No crystals were observed before and during decompression, thus restricting crystallization, of plagioclase exclusively, to the dwell stage at P f during time periods up to 17 days. However, degassing was fully completed before the dwell stage at P f , based on direct measurements of H 2 O content in the experimental glasses. The driving force for crystallization, T eff , was calculated to 55, 80, and 110 • C for P f of 100, 75, and 50 MPa, respectively. As pointed out in Mollard et al. (2012), the three decompression rates did not lead to systematic differences in phase composition, proportions, and crystallization rate for a given P f , except for the runs decompressed to 100 MPa at 1,200 MPa/h, which show nucleation difficulties and have thus been ruled out from the calculation. The experimental conditions are reported in Table 1.
The area number density (N A ), proportion (Φ) and size of the plagioclases were determined by image analysis (using SPO2003 software; Launeau and Robin, 1996) as a function of experimental duration ( Table 2). As defined in Mollard et al. (2012), the crystal size was determined on SEM pictures by the ten longest semilong axes of the crystals (L/2), in order to best approach the true stereological size of the crystals by minimizing cross-cutting effects (Higgins, 1994).  Mollard et al. (2012). MPa in Mollard et al. (2012) and Mt Pelée rhyolitic liquidus in Mollard et al. (2012). d Effective undercooling ( T eff = T L -T). e Melt H 2 O content calculated at P f -T after Newman and Lowenstern (2002) and given with an error of 0.2 wt%. f η is melt viscosity calculated after Hess and Dingwell (1996).  (2000), which considers a monodisperse population of tabletshaped objects (raw data inputs are the longest axes of ellipses with areas equivalent to the 2D sections and "shape" is set to 1.0:1.0:a/b for the short, intermediate and long axis, with a/b being the aspect ratio (long to short axis) of the ellipses; the CSD calculations were corrected for volumetric phase abundance). c Maximum rates of plagioclase nucleation (Imax) determined from time differentiating the Nv functions.

Textural and Analytical Data
For the modeling purpose, we had to perform complementary analyses to the experimental study of Mollard et al. (2012).
Regarding plagioclase textures, we had to convert N A into volume number densities (N v ). Because plagioclase microlites show strongly elongated tablet shapes, we determined N v from the method of crystal size distribution (CSD) using the spreadsheet of Higgins (2000), which considers a monodisperse population of tablet-shaped objects ( Table 2). Plagioclase crystals grow anisotropically, but we lack growth rate data for each of the three crystalline axes, so that the crystal size to model was the average of the ten longest semi-lengths of the crystal long axes. The time evolution of plagioclase N v , Φ, and L/2 are presented in Figure 1. We had also to determine the mean distance between two crystal centers from the SEM images used in Mollard et al. (2012), via the "center-center" function of the SPO2003 software (Launeau and Robin, 1996) (Table 3).
Regarding glass compositions, we performed analytical profiles in SiO 2 , Al 2 O 3 , CaO, and Na 2 O contents from selected decompression run products. The analyses have been performed FIGURE 1 | Experimental data of Mollard et al. (2012) showing the time-evolution of the plagioclase crystal (A) volume number density calculated after Higgins (2000) as in Table 2; (B) plagioclase proportion, and (C) average of the ten longest semi-axes of the crystals as defined in Mollard et al. (2012). Temperature was 875 • C, initial pressure was 200 MPa, and final pressures were 100, 75, and 50 MPa, corresponding to T eff of 55, 80, and 110 • C, respectively. Note that the curves start at different times, corresponding to an increasing delay in plagioclase nucleation with decreasing T eff . However, for growth modeling, 0 time corresponds to the onset of nucleation.
using a Cameca SX50 electron microprobe (EMP; BRGM-ISTO, Orléans), with the conditions of 15 keV accelerating voltage, 6 nA beam current, 10 s counting time on the element peaks, and a defocused beam (∼10 µm in diameter) to minimize Na migration during analysis (Na analyzed first; Pichavant, 1987;Devine et al., 1995). Analytical errors on these analyses were calculated at ∼1% (relative) for SiO 2 , Al 2 O 3 and CaO, and 5% for Na 2 O. a Characteristic length of the system corresponding to the semi-distance between two crystal centers, determined from image analysis on the experimental samples. b Diffusion model following Fick's equations for a finite reservoir and an advancing boundary layer, with i) fixed a (from the experiments) and ii) adjusted D CaO and a calc in order to best fit the experimental data.

Classical Nucleation Theory (CNT)
Nucleation may be homogeneous or heterogeneous. Heterogeneous nucleation occurs if any interfaces, such as pre-existing crystals or bubbles, are capable of decreasing the energy barrier required to trigger crystal nucleation. The starting products used in the experiments of Mollard et al. (2012) were aphyric and the SEM images did not show pre-existing gas bubbles that could have acted as preferential sites for plagioclase nucleation (i.e., no systematic crystal attachment on bubbles). Therefore, we assume that nucleation in the experiments was homogeneous. The CNT has been developed by Volmer and Weber (1926), Farkas (1927), Kaischew and Stranski (1934), Becker and Döring (1935), improved by Turnbull and Fisher (1949), Frenkel (1946), Kashchiev (1969), and compiled in Kirkpatrick (1981) and Fokin et al. (2006). The CNT, for the case of homogeneous nucleation, considers a submicroscopic process in which random fluctuations in a homogeneous one-component liquid held below its liquidus form crystal-like atom clusters, in order to reach thermodynamic stability (low free energy). This results in the creation of a cluster-liquid interface, the maintaining of which causes clusters below a critical size to shrink whereas the largest grow to macroscopic sizes (the critical size decreasing with increasing undercooling). The temperature-dependence of the nucleus size distribution is given by the Boltzmann function: where n is the number of spherical nuclei, n o is the total number of atoms, G is the excess free energy, and k B the Boltzmann's constant. According to the CNT, the steady-state rate of homogeneous nucleation has the following form James (1985): where N c = A/V m , with N c being the volumetric concentration of reactant atoms, A the Avogadro's constant and V m the molar volume of the crystallizing phase, h is the Planck's constant, G * is the free energy required to form a spherical critical nucleus having the properties of the macroscopic crystal, and G D is the activation energy required for atom attachment. G * can be expressed as: where σ is the free energy associated with the crystal-melt interface and G v the free energy change per volume of the transformation that can be defined as G v σ G/V m .
To compare experimental data with the CNT, several simplifications are necessary. The first one sets that the activation energy of atomic jumps across the nucleus-melt interface equals that of the shear relaxation time of the melt (Stokes-Einstein approximation), allowing G D to be expressed in terms of melt viscosity (η). The second assumption sets that heat capacities of the crystal and melt compare, reducing G to Turnbull's equation: where H is the enthalpy of formation and T is the undercooling which is replaced by T eff in the present case of decompression-induced isothermal crystallization. Combining (2) and (3), the steady-state nucleation rate becomes: where a o is the atomic jump distance. This shows that for small T eff , the exponential term controls the shape of the curve and I increases with increasing T eff . In contrast, for high T eff , η increases faster than the exponential term and I passes per a maximum before decreasing. Thus, I should show a bell-shape and an asymmetrical curve around the maximum depending on T eff . To resolve (5), thermochemical and viscosity data are available and I is empirically determined from the crystal number density evolution with time in the experiments, whereas σ remains the most unknown parameter.

Diffusion-Driven Crystal-Growth Theory
A nucleus that reached a critical size starts growing until the thermodynamic equilibrium is reached or the supply in crystal components stops. The theory of crystal growth has been extensively developed by Christian (1965), Dowty (1980), and Kirkpatrick (1981). In principle, crystal growth involves two main mechanisms: (1) atom incorporation at the crystal/melt interface which locally leads to melt depletion in species consumed by the crystal, and thus leads to a concentration gradient that drives (2) transportation of the plagioclasecomponent atoms from melt to crystal interface. The kinetics of either mechanism control crystal growth. In the following, however, interface processes are not considered as the growth limiting factor, practically meaning that the atoms are assumed to be instantaneously inserted into the crystalline structure.
Although this hypothesis should be demonstrated, our first approach of plagioclase growth modeling is to test atom diffusion in melt as the limiting factor on the entire time period of crystal growth.

Adaptation of Fick's Law to Multicomponent Systems
Fick's second law classically describes diffusion in binary aqueous solutions. Therefore, its application to multicomponent silicate melts requires the assumption that the silicate melt can behave as an effective binary system (Cooper, 1968;Zhang, 2010). Given the high silica content of the melt (∼80 wt%), we defined SiO 2 as the solvent solution, in which solutes of Al, Ca, and Na, diffuse. As the crystal/melt interface processes are not considered, the progression of the crystallization interface only depends on the flux of the various components arriving at the interface. It is reasonable to assume that growth is limited by only one such flux, i.e., the smallest one. Note that in the present case of a finite reservoir (chemically closed system), the flux at the interface depends on the diffusion coefficient of the considered component but also on the concentration gradient at that point. Thus, comparing fluxes does not necessarily mean comparing diffusion coefficients (in contrary to the case of an infinite reservoir with continuous component supply). The species with the lowest diffusion flux may as well be the one with the highest diffusion coefficient, because it starts depleting the reservoir faster than the other species.

Dimensional System
Diffusion in binary aqueous solutions is classically described by Fick's second law (Bird et al., 2007) which gives the solute concentration c(x,t) at distance x from the interface and time t as follows: where D is the diffusion coefficient of the considered species. At t = 0, the initial concentration is C o and the concentration at the crystal-liquid interface is the equilibrium concentration C eq (below which the component cannot precipitate anymore into the crystal). The position of the crystal/melt interface at time t is noted x B (t) and the growth rate v B (t) is its time derivative: The flux J at the interface is as follows: The fact that growth rate (i.e., interface velocity) is controlled by diffusion is expressed by the Stefan condition (Ghez, 1980): with C s being the species concentration in the crystal. The righthand term of the equation is the quantity of matter (per unit square per second) that is consumed by the crystal as the interface extends with velocity v B . The left-hand term is the quantity of matter brought at the interface by diffusion. The two are forced to be equal following the principle of mass conservation. Therefore, estimating plagioclase growth velocity is equivalent to computing a diffusive flux at the interface, which requires solving Fick's law (6) with the appropriate boundary conditions. Considering a finite reservoir requires the determination of the characteristic length of the system (a) over which diffusion operates, which we retrieved from the average semi-distance between neighboring crystal centers (image analysis on the SEM pictures of the experimental samples using the SPO2003 software; Table 3). The position of the interface x verifies 0 < x < a. The boundary conditions are: (10) postulates that at t = 0 (when nucleation starts), the concentration in the system is equal to the initial concentration c o . (13) gives the boundary condition at the interface (at position x B ) where the concentration is fixed at the value c eq . (12) is the boundary "closed system" condition at the other end of system, which states that the flux through this extremity is zero. This hypothesis reflects the fact that there is another crystal growing, opposite to this one, on the other side of this boundary. However, because of the diffusion coefficient issue addressed below, we had to circumvent this dimensional system of equation into a non-dimensional one.

Non-dimensional System
To solve (6), one requires knowing D of the growth-limiting component (CaO in our case; see below). Yet, because of the absence of CaO diffusion data relevant to our conditions (see below), we deduce D from the experimental data using a non-dimensional form of (6). The solution is also in a non-dimensional form (i.e., independent of D) and must be transformed back into a dimensional form using an adequate D.
The non-dimensional form of this system is obtained via the change of the following variables: with their respective differentials (dx'=a dx; dt'=a²/D dt; dc'=(c o -c eq ) dc), the D-independent form of (6) writes: with boundary conditions as follows: This system is solved numerically using a fixed finite difference scheme and a front tracking method to follow the movement of the interface (i.e., giving the position of the interface). The details are available in Crank (1984). The method works as follows: 1. solve the non-dimensional form of Fick's law (16), which returns a D-independent concentration field, a flux at the interface, and an interface velocity [v B (t,D)] Crank (1984); 2. a D-independent interface position is obtained by time integration [x B (t,D)], which is exactly the half-size of the crystal; 3. find D that best fits this computed crystal size with the experimental measurements [crystal size vs. time].
It has to be noted that for comparison, the same procedure was performed using the Fick's equations for a concentric growth [finite reservoir and advancing boundary layer; equations in Crank (1984)]. The results turned out to be very similar to those of the model described above (i.e., planar interface), so that only this latter model is presented here.

Nucleation Rate and Interfacial Energy
The energy required to create interfaces between the solid nucleus and melt (σ) plays a critical role in the crystallization process by drastically impacting the nucleation rate (I), with tiny σ variations generating several orders of magnitude differences in I (Figure 2; Hammer, 2004). Indeed, strong σ forces the magma to become highly supersaturated before crystals can nucleate. Should σ become great enough, crystals would not nucleate, even if the driving force for nucleation, T eff , is high. Solid-liquid σ has been determined independently from the CNT, mostly in pure metals (Howe, 1997) and from indirect estimations in magmatic system (Zanotto, 1987). For instance, macroscopic σ values have been obtained by measurements of crystal-melt wetting angles, giving values of 0.3 J.m −2 for quartz-granite (Laporte and Provost, 2000) and 0.5 J.m −2 for olivine-basalt (Cooper and Kohlstedt, 1982). Thus, we focus below on defining σ for the plagioclase-haplotonalite (rhyolite) system. The following empirical equation gives macroscopic σ estimations for different compositions (Turnbull, 1950;Kirkpatrick, 1981;James, 1985): where β is an experimentally-determined constant that ranges from ∼0.3 for non-metals to ∼0.5 for metals (Matusita and Tashiro, 1973;Neilson and Weiberg, 1979). In the case of our experimental data set, P i is 200 MPa and the plagioclase composition varies from An 27 to An 37 depending on P f (Mollard et al., 2012), so that H varies from 74300 to 79180 J.mol −1 (Burnham and Nekvasil, 1986;Muncill and Lasaga, 1988;Hammer, 2004) and V m , calculated as a linear albite-anorthite solution and assumed to be pressure independent, is ∼1.003 × 10 −4 m 3 .mol −1 (Klein, 1994;Hammer, 2004). From (20) and with β = 0.3, σ emp varies from 0.122 J.m −2 for An 27 to 0.130 J.m −2 for An 37 (Table 4). However, these values may still be overestimated with respect to a microscopic critical nucleus (Fokin et al., 2000).
In the absence of additional methods to determine plagioclaserhyolite σ independently from the nucleation rate data, we calculated I with σ emp (20) using the CNT rate Equation (5), giving log Imax below 20 m −3 .s −1 .

Glass Compositional Gradients
The glass compositions at initial and final conditions are reported in Table 5 and their compositional evolutions as a function of distance from the plagioclase interface for durations of 7 and 17 days at P f (100, 75, and 50 MPa) are shown in Figure 3.
For T eff = 55 • C (P f of 100 MPa) and a crystallization duration of 7 days, the SiO 2 content of the glass next to the plagioclase is similar to the one expected at P f (i.e., 80.8 wt %) and gradually deviates toward the initial concentration (78.7 wt %). This trend is confirmed by the Al 2 O 3 , CaO, and Na 2 O profiles, which all start at the plagioclase interface close to the value expected at P f and join the P i concentration 80-100 µm away from the crystal. The slopes of the Al 2 O 3 , CaO, and Na 2 O profiles are inverted with respect to SiO 2 , since plagioclase is SiO 2 -poorer than melt, but richer in Al 2 O 3 , CaO, and Na 2 O. The presence of compositional gradients in the glasses after 7 days suggests that none of the oxides chemically re-equilibrated at 100 MPa over this period of time.
For T eff = 80 • C (P f of 75 MPa) and crystallization duration of 7 days, the SiO 2 , Al 2 O 3 , and Na 2 O profiles in glasses show similar trends than for T eff = 55 • C. After 17 days, all profiles are flat and close to the expected concentrations at P f , suggesting complete chemical re-equilibration between 7 and 17 days.
For T eff = 110 • C and crystallization duration of 7 days, the profiles are close to the expected values at equilibrium, but they show gradients, so that the glass did not fully reequilibrate. The Na 2 O profiles at 7 days are richer than the starting glass concentration, an inconsistency that might be attributed to the large analytical errors associated with Na microprobe analysis (although precautions have been taken for the analysis, as mentioned above). It has also to be noted that the small difference between the initial and equilibrium compositions do not favor accurate estimations of the gradients. For crystallization duration of 17 days, the profiles are more or less flat (at least they start and end with similar values), but they show differences in bulk oxide contents. For instance, profile #7 shows flat profiles close to equilibrium compositions, whereas profile #8 shows 2-3 wt% SiO 2 less and 1-2 wt% Al 2 O 3 more, that is, closer to the initial concentrations. This suggests that melt areas may be chemically re-equilibrated, whereas others are still in progress after 17 days. Such chemical heterogeneities may result from (i) spatial heterogeneities in plagioclase crystal distribution as described in Mollard et al. (2012), (ii) a nucleation event that is not strictly instantaneous, possibly due to the slow kinetics prevailing at low melt H 2 O contents; this would result in transient plagioclase crystallization from chemically evolving melts, or (iii) the vicinity of neighboring crystals that contaminate some glass analyses (mixed analyses between glass and crystal).

Crystal-Growth Limiting Component
The component with the smallest flux that limits crystal growth can be identified from comparing the compositional profiles in the experimental melts as a function of distance from plagioclase interface (Figure 3). The source of a given component is assumed to be exhausted when the component profile shows a flat slope at the equilibrium concentration expected at P f . Al 2 O 3 did not reach equilibrium at P f after 7 or 17 days (positive slopes and profile #8 far from equilibrium concentration). The Na 2 O profile looks flat at 17 days, but the large analytical scatter makes it inherently untrustworthy. In contrast, CaO is the only component showing equilibrium concentrations after 7 days (and also after 17 days) at P f = 75 MPa. CaO is also the component with the lowest concentrations (∼1-2 wt%) and the smallest slopes of the profile lines (which give an indication of the magnitude of the flux), explaining why it gets depleted first. Therefore, we chose CaO as the plagioclase growth limiting component.
Thus, solving (6) requires knowing D CaO . Ca diffusivity values are particularly variable in dry rhyolitic melts, with logD CaO ranging from about −11 m 2 /s at 1,400 • C to −13 m 2 /s at 900 • C [data compiled in Zhang (2010)]. In H 2 Orich (3.6-6.3 wt%) melts, Baker and Bossanyi (1994) measured logD CaO of the order of −11 to −10 m 2 /s for temperatures from 1,100 to 1,400 • C, highlighting the strong effect of H 2 O on the diffusion coefficient of CaO. In a Ca-doped (<0.1 wt%) haplogranite, Mungall et al. (1999) measured logD CaO from −14 m 2 /s at 1,137 • C to −11 m 2 /s at 1,600 • C in dry conditions, reaching −11 m 2 /s at 1,300 • C and −10 m 2 /s at 1,600 • C with 3.7 wt% added H 2 O. In conclusion, none of the above D CaO have been determined for pressure, temperature and H 2 O content conditions relevant to the experiments of Mollard et al. (2012), Therefore, we deduce D from the experimental  (Mollard et al., 2012). b σemp is empirically determined at macroscopic scale using (6) with β = 0.3, Vm =1.003x10 −4 m 3 .mol −1 , and H calculated for An (see text). c σ calc is calculated at microscopic scale by rearranging the CNT rate Equation (5), using the experimental Imax given in Table 2. σ rug(α;ds) is a σ calculation assuming non-spherical and rough-faceted nuclei (with a shape parameter α > 4π and a rugosity parameter ds > 2), which requires the use of Equation (7) in (5)   data using a non-dimensional form of (6). The solution is also in a non-dimensional form (i.e., independent of D) and must be transformed back into a dimensional form using an adequate D.

Results of the Plagioclase Growth Simulations
Following the methodology given in section § Non-Dimensional System, we performed simulations at P f = 100, 75, and 50 MPa, the results of which are given in Table 3 and  (Figure 4C).
Fixing a to the measured values does not allow a good fitting of the experimental data, but the calculated D CaO are of the same order of magnitude. In the worst case, i.e., for T eff of 110 • C, the calculation returns D CaO of 0.7 × 10 −15 m²/s instead of 1.1 × 10 −15 m²/s when calculated for both parameters being adjusted (i.e., about twice as slow; Figure 4).

DISCUSSION
In the section below, we discuss the results of the crystallization modeling in terms of validation and implications for volcanic systems.

New Constraints on the Plagioclase-Rhyolite Interfacial Energy (σ ) σ Values Deduced From the CNT
The calculation of I with σ emp (20) using the CNT rate Equation (5) gives extremely low log Imax (<20 m −3 .s −1 ), which represents inconsistent values with respect to the experimentallydetermined log Imax of 8 to 12 m −3 .s −1 with T eff increasing from 55 to 110 • C (P f from 100 to 50 MPa, respectively) obtained by time differentiation of the N v functions ( Table 2). This suggests that either the CNT or the calculation of σ is not valid in this case.
Turning things around, we assumed the basic form of the CNT Equation (5) and back-solved for σ using our  Mollard et al. (2012). Each column corresponds to a given run product at T eff of 55 • C (blue), 80 • C (green), or 110 • C (red), and a crystallization duration of 7 or 17 days. The symbols represent different analytical profiles performed in a given run product. The black dotted lines represent the initial concentrations (C o at 875 • C and 200 MPa) and the colored dashed lines give the expected equilibrium melt compositions at crystallization pressure (C eq at 875 • C and P f ) as given in Table 5. The vertical bars on the right give the analytical error.
experimentally-determined nucleation density data (I max ). The calculations suggest σ calc from 0.041 to 0.059 J.m −2 with T eff increasing from 55 to 110 • C (Figure 2), which is about 2-3 times lower than σ emp (Table 4). Following a similar approach (i.e., assuming the basic form of the CNT and back-solving for σ using the experimental nucleation data of Hammer and Rutherford, 2002), Hammer (2004) suggests σ variations from 0.045 to 0.068 J.m −2 with H 2 O content decreasing from 3.8 to 2.5 wt% (i.e., P f from 100 to 50 MPa and T eff from 55 to 110 • C). Both results are in very good agreement, giving confidence in the use of such an approach to estimate σ values.
However, these calculated σ may deviate from the real σ because it incorporates all the assumptions inherent in CNT (Hammer, 2004;Shea et al., 2009). In addition, another source of variance is the dwell time during which textural ripening occurs, which may lead to a decrease of the number density during the period between nucleation (first hours to days after perturbation) and the growth of crystals to observable size (7-17 days); such an effect of dwell time may influence the comparison between experimental studies of otherwise similar systems (feldspar in highly evolved silica-rich liquids).

Compositional Dependence of σ
The T eff dependence of σ (Figure 2) is mainly ascribed to the melt H 2 O content (Hammer, 2004). The decrease of melt H 2 O content with decompression (increasing T eff ) has three main effects. The first one is to increase melt viscosity and thus the kinetic barrier to the atomic attachment (Heslin and Shelby, 1993). The second is the influence of devolatilization on the free energy term in the nucleation rate equation (Hammer, 2004;Shea et al., 2009). The third effect concerns the progressive divergence between the liquid and plagioclase compositions in the course of crystallization that leads to an increase of the crystal-melt σ (Takei and Shimizu, 2003). This effect results from intra-crystalline chemical bonds less energetic than at crystal-liquid interface, so that the crystal-liquid σ will increase as the concentration of the crystal component dissolved in the liquid decreases. In other words, the nonideality of chemical bonding causes a strong compositional  Table 3. The red dots represent the experimental data, as given in Figure 1C. The plain blue curves show the results for the model adjusting both, the characteristic length of the system (a calc ) and D CaO , in order to best fit the experimental data. The dashed green lines show the results for the model adjusting D CaO with a fixed to the experimental value.
FIGURE 5 | Dependence of the interfacial energy (σ) on the nucleus shape (α) for a roughness d s of 2.3 (An 30 plagioclase), as a function of effective undercooling ( T eff ). σ calc(4π;2) and σ emp (as determined in Table 4) define a range of σ bracketed between microscopic and macroscopic values, respectively.
dependence of the crystal-liquid σ (Takei and Shimizu, 2003;Shimizu and Takei, 2005). All three effects may act together to increase σ as the H 2 O-saturated rhyolitic melt depressurizes and crystallizes plagioclase.

Refinement of σ Through Shape-Dependent Nucleus Parameters
Because the CNT assumes smooth-faceted and spherical nuclei, σ calc may be wrong in case of rough and elongated nuclei. Indeed, a rough surface increases the surface area, which modifies the crystal-melt σ and thus nucleation rate (Sen and Mukerji, 1999). Although the nanometric size of the critical nucleus does not permit determination of its true shape characteristics, an alternative consideration to the CNT is to view the critical nucleus as the unit cell shape of the macroscopic crystal, i.e., triclinic and rough-faceted for plagioclase (Holdren and Speyer, 1987;Brantley and Mellot, 2000). Thus, assuming that the general morphology of the nucleus compares to the one of the macrocrystals, the surface fractal dimension (d s ) may be extracted from the surface area (S A ), using S A = r ds−3 , where r is the radius of the nucleus (Mandelbrot, 1982;Van Damme et al., 1987). Measurements on albite and anorthite crystals suggest d s of ∼2.2 and ∼2.7, respectively (Holdren and Speyer, 1987;Brantley and Mellot, 2000). Assuming an ideal solid solution between albite and anorthite gives d s value of ∼2.3 for An 30 (average plagioclase composition analyzed in the experiments).
The nucleus shape, expressed as the shape-dependent parameter α, has the value of 4π for a perfect sphere, but there are no α data available for a plagioclase shape. The most relevant d s -α data to our magmatic system concern congruent nucleation in Li 2 Si 2 O 5 liquids (that have often been taken as melt analogs), which show a positive correlation between d s and α, with α values of ∼400π where d s = 2.3 (Sen and Mukerji, 1999). The influence of d s and α on σ is expressed through the following equation: (21) Assuming d s = 2.3 for An 30 and using our experimental I max , we calculated σ values for rough and non-spherical nuclei, σ rug(α;ds) , using (5), in which (21) replaces (3), and for α varying up to 2,000 π. Bracketing σ rug between σ emp and σ calc (calculated for a smooth sphere) suggests α between 400 π and 1,500 π ( Figure 5). σ emp would be compatible with σ rug having d s = 2.3 and α of 483, 635, and 792 π at P f = 100, 75, and 50 MPa, respectively ( Table 4). Although there are no available data on a correlation between α and a given crystal aspect ratio, the elongation and the increasing surface area of the nucleus as T eff increases would agree with the observed evolution of the plagioclase macrocrystal characteristics.

Validation of the Process Controlling Crystal Growth
The theoretical results and the experimental decompressioninduced plagioclase growth laws agree on three points (Figure 4). First, the experimental data agree with the theoretical curve profiles defined for a finite model with an advancing crystalmelt interface limited by CaO diffusion. Second, the returned a calc , i.e., the distance between two crystals that defines number density, are not far from those measured in the experimental samples, i.e., differing by ∼60 µm at P f = 100 MPa, ∼80 µm at P f = 75 MPa, and ∼10 µm at P f = 50 MPa (Table 5). Third, the returned D CaO from 1.1 × 10 −15 to 8.4 × 10 −14 m 2 /s stand in the lower limit of the range of those from the literature (full range of 10 −14 to 10 −10 m 2 /s), although none of them were determined for conditions (melt composition, H 2 O content, temperature, pressure) comparable to those in Mollard et al. (2012).
Therefore, this overall good agreement between the theoretical and empirical models validates CaO diffusion as the main process controlling isothermal decompression-induced plagioclase growth under moderate T eff (<110 • C).

Limitations of the Model
Constraining D CaO by an independent method would of course help to refine the model, which, in its present form, possibly compensates different kinds of discrepancies by adjusting D CaO . The discrepancies may come from limitations on Fick's laws and/or crystal-melt interface barriers, as explained below.
First, the aqueous solutions for which the Fick's laws were determined are often chemically dissociated, whereas melts contain polymeric transient units, possibly leading to coupled diffusions of some species. Also, aqueous solutions are dilute, thus permitting to approximate diffusion of particles against a constant background of water, whereas the notion of particles diffusing against a constant background of a silica melt has not been discussed so far. The viscosities in aqueous solutions and melts differ by many orders of magnitude, which may introduce some bias in the formulation. Although these issues need to be addressed in the future, the application of Fick's law to multicomponent systems is a reasonable starting point of work, as already used by many authors willing at determining diffusion coefficients in complex systems such as polymers, glasses, and melts.
Second, plagioclase growth may not be limited by kinetics (diffusion) over the entire crystallization period of time. In particular, Roskosz et al. (2005a,b) demonstrated that at large undercooling, for which melt viscosity is high and the crystalmelt Al/Si ratio is constant, crystal growth is thermodynamically controlled by crystal-melt interface processes. Yet, the calculation of the Al/Si ratios in the experiments used in the present study does not consolidate this hypothesis (Mollard et al., 2012). In the opposite side, interface-controlled growth may be the regime prevailing for small degree of undercooling in case of incongruent crystallization, since plagioclase-melt interfacial energy is high and the component diffusion is often fast compared to the rate of atom attachment. If true for the early crystallization, the starting conditions used for Fick's diffusion modeling may not be correct. However, the equations governing plagioclase growth controlled by interface processes are not time-dependent and consequently do not allow to assess when diffusion actually takes over as the controlling mechanism. Interface-controlled mechanisms may likely affect the early plagioclase crystallization at T eff = 55 • C (P f = 100 MPa), but the overall shape of the experimental plagioclase growth is satisfactorily simulated by diffusion-controlled processes.
Thus, potential biases, such as the application of Fick's law to melts and the possible interface-controlled growth, need further evaluation to fully validate the present model and the CaO diffusion coefficients.

Implication for Volcanic Systems
The present study provides a mean to retrieve crucial parameters for the understanding of syn-eruptive volcanic processes. Indeed, the simulations may be relevant to phenocryst-poor rhyolitic magmas that rapidly ascend from a storage region to a given shallower depth where they crystallize microlites for a given amount of time. Although rapidly-decompressed magmas that crystallize at constant pressure may not be the most common scenario triggering explosive eruptions (more commonly thought to follow polybaric crystallization), dome-generated blasts have been recognized to follow such a decompression-induced isobaric crystallization path (e.g., Castro and Dingwell, 2009;Martel, 2012). The driving force for crystallization, T eff , may be deduced from the plagioclase liquidus for the bulk composition of interest and the crystallization pressure/depth that is either directly determined (e.g., from geophysical measurements) or retrieved post-eruption from the pyroclasts [i.e., from number density, proportion, and shape of the plagioclase microlites as defined in Mollard et al. (2012); their Figure 11]. This study has two major implications for the understanding of the volcanic systems.
Firstly, predicting crystal nucleation rates in volcanology is crucial for two reasons. First, high nucleation rates may lead to an extensive crystallization of microlites at low pressures (<50 MPa), which drastically enhances the risk of dome explosivity (Sparks, 1997;Martel, 2012). Second, nucleation rate directly controls crystal growth kinetics by fixing crystal number density that defines the characteristic distance on which diffusion proceeds between two growing crystals. The shorter the characteristic distance, the quicker the reservoir depletion in the fastest species, and thus the slower the growth rate. Therefore, magmas with high crystal number densities must show slower crystal growth rates than poorlycrystallized magmas, all other conditions being comparable. Our study suggests that homogeneous nucleation rates in rapidly decompressed natural rhyolites may be deduced from the CNT equation, using a plagioclase-melt surface tension from 0.04 to 0.06 J.m −2 with T eff increasing from 55 to 110 • C. Yet, natural rhyolitic magmas may contain impurities that could trigger heterogeneous nucleation by lowering crystalliquid surface tensions and the energetic barrier, for instance on pre-existing crystals or gas bubbles. So far, however, neither experimental or natural products have brought evidence of preferential sites favoring plagioclase nucleation during magma decompression.
Secondly, the size evolution of the plagioclase microlites may be predicted as a function of time at various crystallization pressure (or T eff ) by the present model of growth limited by the diffusion of CaO. Our model predicts that a H 2 Osaturated rhyolitic melt at 875 • C that ascends from about 6-8 km deep (∼200 MPa) to 1.5-2.0 km (∼50 MPa) will nucleate about 10 13 m −3 plagioclase microlites ( Table 2) with final half-lengths of ∼60 µm ( Figure 4C). Additional data at lower pressures would be helpful to assess the crystallization conditions at dome level. However, coupling data on crystal number density to crystal size evolution is a key parameter to assess magma flow properties. Indeed, many small needle-shaped crystals do not behave the same as few large crystals under shear. Integrated numerical modeling of magma crystallization during ascent in volcanic conduits (e.g., Melnik and Sparks, 1999;De'Michieli Vitturi et al., 2010) consider crystallization as a volume fraction at thermodynamic equilibrium. Such models would simulate magma rheology more realistically by incorporating time evolutions of textural parameters, such as crystal number density and size.
For a fully confident extrapolation to natural systems, the calculations of the crystallization kinetics still need further investigations. In particular, the chosen composition represents the SiO 2 -richest end-member of the typical interstitial melts from natural andesites to rhyolites (70-79 wt % of SiO 2 ; Hammer, 2008). A validation of the calculations for SiO 2 -poorer melts would be necessary, especially to verify whether Fick's diffusion law is still acceptable with a SiO 2 -poorer background solution. Additionally, the calculations are performed for a four-components simplified analog of a rhyolitic melt. Although such analog compositions have been successfully used to determine various physicochemical properties (including diffusion) of natural magmas, they have to be confronted to more realistic melts, the compositional complexity of which may induce deviations of the crystallization processes and kinetics. Finally, crystallization at fixed T eff needs to be extended to polybaric crystallization, for which crystal nucleation and growth experience time-evolving T eff .

CONCLUSIONS
Modeling the experimental results of decompression-induced crystallization of plagioclase from a rhyolite-like melt permits to refine the plagioclase-rhyolite surface tensions that control crystal nucleation rates and validates plagioclase growth limited by diffusion. In particular: -Plagioclase-rhyolite interfacial energies (σ) have been determined from 0.041 to 0.059 J.m −2 with T eff increasing from 55 to 110 • C. -Considering the shape characteristics of the plagioclase nucleus as similar to those of the macrocrystal, σ can be refined by adjusting parameters of nucleus roughness (d s ) of ∼2.3 and shape (α) between 500 and 1,350 π depending on T eff . -Plagioclase growth with time can be numerically simulated by a diffusion-limited model (Fick's second law adapted to multicomponents), considering a (i) finite reservoir, (ii) an advancing crystal-melt interface, and (iii) CaO as the component limiting plagioclase growth. -The CaO diffusion coefficients, D CaO , that best fit the experimental data range from 8.4 × 10 −14 to 1.1 × 10 −15 m 2 /s at 875 • C and H 2 O-saturated pressures from 50 to 100 MPa. -Predicting crystal nucleation (number density) and growth (time-evolving size) is of prime importance for assessing the rheological properties of magmatic flows that decipher eruptive dynamics, with a direct application to plagioclase microlite crystallization in ascending silica-rich magmas.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/supplementary material.

AUTHOR CONTRIBUTIONS
EM participated to the implementation of the model based on her own experimental data. CM participated to the model implementation and overall writing of the manuscript. EL wrote the program using the diffusion equations. GR performed tests using the program. All authors contributed to the article and approved the submitted version.

FUNDING
This work was part of EM Ph.D. (French MNRT grant) and was funded by the Agence Nationale de la Recherche (contract No. ANR-05-CATT-0004 to CM).