Tracking Crystals Evolution in Episodes of Magma Mixing

1. Introduction

Eulerian formulations are widely employed in magma fluid-dynamics because they naturally solve the Navier–Stokes equations on fixed/moving computational meshes, provide robust and scalable simulations over large domains such as magma reservoirs, and efficiently handle multiphase and multicomponent mixtures. (Batchelor 1967,Garg Papale 2022,Ishii Hibiki 2011,Prosperetti Tryggvason 2007). While the evolution of velocity, pressure, temperature and composition are successfully obtained over the whole domain, track is lost of the pathways of magma parcels. Lagrangian trajectories of dispersed passive tracers allow to determine the locations of magma “particles” that move along the Eulerian flow. As a consequence, for each instant of motion and for each coordinate of the main flow, it can be known the original position of magma. Further, along each Lagrangian path, there can be traced the chemical and thermodynamic local conditions to which each magma particle is subject.

During magma mixing, suspended crystals can act just as passive Lagrangian particles, if their volume fraction is such that they do not disturb the overall flow by changing magma density and viscosity. While they transverse different environments, they record the pressure, temperature and composition by accretion of successive crystallizing layers.

Petrology and geochemistry reconstruct the underground magmatic conditions from such crystals’ accretionary layers (John S. Armstrong-Altrin 2022,Leonid L. Perchuk 1991). They can discern the original compositions that undergo mixing, and the depth at which it occurred. However, they register the final product of crystallization, but they cannot define precisely the locations of these crystals inside the reservoir and follow continuous ambient changes during time. This work is aimed at filling the gap between the backward reconstruction of magmatic conditions from eruptive products and the processes occurring at depth.

Previous numerical studies accounted for crystals redistribution and grow during episodes of gas-driven overturn in a single magma chamber (Ruprecht et al. 2008), and mixing in cases of injection of fresh melt into a basal crystal mush (Bergantz et al. 2015,Bergantz et al. 2017,Carrara et al. 2020,Schleicher Bergantz 2017,Schleicher et al. 2016). Ruprecht et al. (2008) solve at first magma fluid-dynamics, and adopt the Bassinet-Boussinesq-Oseen equation to describe the force-balance for crystals’ spherical particles to determine their motion. Afterwards, they study the overall motion of the seeded crystals, arguing on their grow and dissolution on the basis of the Damköhler number. Carrara et al. (2020) and references therein solve the Navier-Stokes equation for the fluid part, coupled with the second Newton’s law for crystals. They study the reorganization of the crystal network during the intrusion.

The innovative contribution of this study relies in the application of Lagrangian trajectories to magma mixing. And to do this in a complex underground system, made of deep and shallow reservoirs connected through a vertical dike, with their respective deep and shallow resident magmas. Indeed, mixing of magmas typically occurs between a less evolved, high-pressure, volatile-rich end-member, and a more evolved, low-pressure, volatile-poor end-member. These magmas reside respectively in a deep and a shallow reservoir (Araya et al. 2019,Foroozan et al. 2011,Gudmundsson 2020,Iguchi et al. 2013,Murakami 2001,Suzuki et al. 2013,Tomiya et al. 2010), which are connected through one or more conduits or dikes. In such systems, mixing occurs inside the shallow chamber and the upper part of the dike. In this work, the initial locations of magma involved in this mixing are determined, and their precise recirculation trajectories are followed.

Petrological studies modelled crystals kinetic as well as melt and intracrystalline diffusion with thermodynamic theory (Brandeis Jaupart 1987,Kirkpatrick et al. 1976), high-pressure experiments (Bonechi 2020), and from observations (Agostini et al. 2013,Cashman 2020,Cashman Marsh 1988). But environment conditions typical of magmatic reservoirs are statically reconstructed in these studies.

Here, for the first time, the continuous change in solidifying crystals’ phases is tracked by the use of the MELTS code (Ghiorso Gualda 2015) along the Lagrangian trajectories of motion inside a physically consistent underground magmatic system. The feasibility of growing layers and dissolution is discussed on the bases of the rate of change of chemical and thermodynamic conditions along the trajectories.

2. Materials and Methods

Deep-shallow mixing in the pre-eruptive evolution at Campi Flegrei is chosen as representative case (Longo et al. 2023). However, results for this particular volcanic complex go beyond the specific case and apply to general mixing episodes. Petrological studies highlight that a deep, volatile-rich shoshonite mixed with a shallow volatile-poor phonolite (Arienzo et al. 2009,Arienzo et al. 2010,Fourmentraux et al. 2012,Mangiacapra et al. 2008). The previous study by Longo et al. (2023) on magma mixing will be enriched in the present work by calculating Lagrangian trajectories of crystals as well as their evolving size and composition.

The geometry of the two-chambers system is designed on the basis of geophysical and petrological studies (Arienzo et al. 2010,de Siena et al. 2010,Di Renzo et al. 2011,Tizzani et al. 2024,Zollo et al. 2008). The shallow and deep chambers are connected through a vertical dike. The shallow chamber hosts the volatile-poor phonolite, while the dike and the deep chamber host the volatile-rich shoshonite (Figure 1). Such setting is representative of an ascending dike that suddenly intersect the shallow chamber, causing mixing of the two end-members both in the shallow reservoir and in the upper part of the dike.

Magma mixing dynamics is computed in an Eulerian reference system , by solving the conservation equations of mass of the two end-members, momentum and energy of the compressible mixture (Longo et al. 2023). Such equations are closed by means of the Lange & Carmichel law for density (Lange Carmichael 1990), a non-ideal solubility model for water and carbon dioxide (Papale et al. 2006), and viscosity that accounts for magma composition, dissolved and exsolved volatiles (Giordano et al. 2008,Ishii Zuber 1979). The two end-members mixture is dealt with as an homogeneous mixture and its properties are computed as weighted means. This implies that, at the aim of Eulerian simulations, mingling occurs.

As regards crystals, it is assumed that their volume fraction is such that it does not affect neither density, neither viscosity. Crystals are therefore treated as passive dispersed particles advected by the main flow. These tracers are uniformly distributed in the fluid domain and their Lagrangian trajectories are computed. Changes in pressure, temperature and composition along these trajectories are also tracked. The equilibrium crystals path for descending pressure is computed with MELTS (Ghiorso Gualda 2015) for steps in intermediate compositions between the two end-members, and the actual crystal content for each point along the trajectories is obtained by interpolation in pressure and in the amount of the end-members. This procedure implies that perfect mixing occurs. At the sub-grid scale, local stirring and diffusion can indeed hybridize the two end-members leading to equilibrium crystallization from intermediate compositions (Insolia et al. 2024,Petrelli et al. 2011). Furthermore, the interpolation over intermediate compositions will always be more accurate than the interpolation between the two extreme end-members.

Eulerian and Lagrangian computations are made with the in-house GALES finite element code (Garg et al. 2018,Garg et al. 2021,Longo et al. 2012) available at https://gitlab.com/dgmaths9/gales, while the GALES version for Lagrangian tracers is available at https://github.com/antonellalongo/gales_lagrangian/releases. A comprehensive description of the code is given at https://gales.pi.ingv.it.

The growth of a crystal may be interface-controlled (attachment and detachment of atoms at the interface), melt diffusion-controlled (diffusion of components in the melt), or solid-state diffusion-controlled (intracristalline diffusion). The timescales or length scales of the three process has to be compared to correctly compute the crystals’ growth.

The kinetic law for growth/dissolution at the interface can be written as:

$$v={\displaystyle \frac{dR}{dt}}=kS$$

(1)

where v is the growth/dissolution velocity, R is the radius of the growing/dissolving crystal, k is the interface-controlled growth rate, and S is the supersaturation or other normalized growth driving force. This driving force should be the activity of the crystal with respect to the background melt phase. This quantity, even though exact, is arduous to be obtained, and a more simple option is taken:

$$S={\displaystyle \frac{m_{eq}-m}{m_{eq}}}$$

(2)

where $m_{eq}$ is the equilibrium mass of crystals in the melt, and m is the actual existing mass before reaching the equilibrium.

Growth is obtained for $S>0$ when the equilibrium mass is greater than the actual mass, and dissolution for $S<0$ when it is lower. In case of growth the equilibrium composition forms a new rim.

According to our simulations, the maximum absolute value of S is ∼1. For typical melts, k varies between 10⁻¹² and 10⁻⁶ m/s (Figure 2b in Moschini et al. 2023, and references therein). Taking as reference a variation in radius of the crystal of $L={10}^{-4}$ m, the timescale of interface-controlled growth is:

$${\tau }_{int}={\displaystyle \frac{L}{kS}}\approx {10}^2-{10}^8\mathrm{s}$$

(3)

The timescale of diffusion in magma can be derived from the flow state, whether it is slow (diffusion dominated) or fast (replacement of impoverished diffusive layer by the main flow). An estimate of the diffusion length $\delta$, the thickness of the compositional boundary layer around a growing crystal is

$$\delta ={\displaystyle \frac{D}{k_m}}$$

(4)

where D is the diffusion coefficient in the melt, and $k_m$ is the mass transfer coefficient set by the flow. For a sphere (Bird et al. 2002,Fogler 2017,Green Southard 2008)

$$k_m={\displaystyle \frac{ShD}{2R}}$$

(5)

where Sh is the Sherwood number. A widely used correlation for the flow past a sphere is (Fogler 2017,Green Southard 2008,Ranz Marshall 1952) is

$$Sh=2+0.6R{e}^{1/2}S{c}^{1/3}$$

(6)

where Sc is the Schmidt number $Sc=\nu /D$ and the Reynolds number $Re=U\left(2R\right)/\nu$. Therefore,

$$\delta \approx {\displaystyle \frac{D}{k_m}}={\displaystyle \frac{2R}{Sh}}$$

(7)

With the typical flow values stemming from Eulerian simulations, $U=0.1$ m/s, $\rho =2400$ kg/m³, $\mu =3100$ kg/ms, with characteristic diffusion coefficient in melt of $D={10}^{-11}$ m²/s, and with the reference crystal size of $R={10}^{-4}$ m, it is obtained $Re=1.5·{10}^{-4}$, $Sc=1.29·{10}^{11}$, and $Sh=14$. Therefore, $\delta \approx 1.5·{10}^{-5}$ m. The timescales of diffusion and of interface-controlled growth are thus

$${\tau }_{diff}={\displaystyle \frac{{\delta }^2}{D}}=20\mathrm{s}$$

(8)

$${\tau }_{int}={\displaystyle \frac{R}{kS}}={10}^2-{10}^8\mathrm{s}$$

(9)

From which it is evident that the process of melt diffusion is negligible.

The role of solid-state diffusion has to be carefully considered through the comparison of the thickness of the growing rim and of the diffusive layer that smooths it. The maximum values of S resulting from simulations is ∼1. The coefficients of diffusion in the feldspar solid-state are in the range 10⁻²¹ – 10⁻¹⁵ m²/s (Brady Cherniak 2010). Considering a reference interval of time $\Delta t=1$ s and $k={10}^{-12}$–${10}^{-6}$ m/s (Moschini et al. 2023), the size of the accretionary rim is

$$\Delta R_{growth}=kS\Delta t\approx {10}^{-12}-{10}^{-6}\mathrm{m}$$

(10)

while the range of possible diffusive layers is

$$l_{diff}=\sqrt{D\Delta t}={10}^{-10.5}-{10}^{-7.5}\mathrm{m}$$

(11)

The diffusive layer may be of the same order of magnitude of the width of the accretionary rim and cannot be neglected.

The size and compositional evolution of crystals is obtained from the equilibrium assemblages computed with MELTS, by the kinetic law (1) for crystals growth or dissolution and from diffusion of components in spherical crystals.

Let us consider a unitary magma parcel that moves along the Lagrangian trajectory from state 1 at $(t-dt)$ to state 2 at t, it carries a mass $m(t-dt)$ of a solid solution (e.g., feldspars). In state 2 the equilibrium mass should be $m_{eq}\left(t\right)$. As stated before, a proxy for the instability of the solid solution in state 2 can be the ratio

$$S={\displaystyle \frac{m_{eq}\left(t\right)-m(t-dt)}{m_{eq}\left(t\right)}}$$

(12)

Therefore, the kinetic law of growth/dissolution can be written as:

$${\displaystyle \frac{dR\left(t\right)}{dt}}=k{\displaystyle \frac{m_{eq}\left(t\right)-m(t-dt)}{m_{eq}\left(t\right)}}=k{\displaystyle \frac{x_{eq}\left(t\right)-x(t-dt)}{x_{eq}\left(t\right)}}$$

(13)

where x are weight fractions of crystals.

Considering the high value of interface-controlled timescale ${\tau }_{int}$ (3), and considering that the growth/dissolution velocity $dR/dt$ in the kinetic law (13) decreases linearly as the equilibrium is approached, there is the actual possibility that the mass of growing/dissolving crystals do not reach $m_{eq}$.

Therefore, while the equilibrium mass $m_{eq}$ is given by the MELTS equilibrium computations, the mass of crystals m carried by the unitary magma parcel has to be calculated from the history of the whole crystals population radius. Assuming dilute spherical crystals, the change in mass in the time interval $(t-dt,t)$ while approaching equilibrium is:

$$m\left(t\right)-m(t-dt)=N{\displaystyle \frac{4}{3}}\pi \left[{\left(R(t-dt)+L\left(t\right)\right)}^3-R{(t-dt)}^3\right]$$

(14)

where $m\left(t\right)$ is the total mass of the new crystals population after growth/dissolution, $m(t-dt)$ was their previous mass, N is the number of crystals, R is their radius, and L is the size of the growth/dissolution rim while reaching the equilibrium. Transforming the expression into crystals weight fractions, it is obtained:

$$x\left(t\right)-x(t-dt)=n{\displaystyle \frac{{\rho }_c}{\rho }}{\displaystyle \frac{4}{3}}\pi \left[{\left(R(t-dt)+L\left(t\right)\right)}^3-R{(t-dt)}^3\right]$$

(15)

where n is the volumetric number density, ${\rho }_c$ the crystals’ density, and $\rho$ the density of magma.

The algorithm that allows to track the crystals’ radius and weight fractions while growing or dissolving towards equilibrium is:

starting at $t=0$ from:

$$\left\{\begin{array}{cc}\hfill L(t=0)& =0\hfill \\ \hfill x(t=0)& =x_{eq}(t=0)\hfill \\ \hfill R(t=0)& =R_{eq}=\left({\displaystyle \frac{x_{eq}}{n\frac{{\rho }_c}{\rho }\frac{4}{3}\pi }}\right)\hfill \end{array}\right.$$

(16)

for $t>0$,

$$\left\{\begin{array}{cc}\hfill L\left(t\right)& =k{\displaystyle \frac{x_{eq}\left(t\right)-x(t-dt)}{x_{eq}\left(t\right)}}dt\hfill \\ \hfill x\left(t\right)& =x(t-dt)+n{\displaystyle \frac{{\rho }_c}{\rho }}{\displaystyle \frac{4}{3}}\pi \left[{\left(R(t-dt)+L\left(t\right)\right)}^3-R{(t-dt)}^3\right]\hfill \\ \hfill R\left(t\right)& =R(t-dt)+L\left(t\right)\hfill \end{array}\right.$$

(17)

where as initial condition it is assumed that the crystals are at equilibrium, as the most simple assumption.

As new rims growth with the equilibrium composition dictated by MELTS, intracristalline diffusion tends to homogenize them. Solid state diffusion follows the standard diffusive equation

$${\displaystyle \frac{\partial x}{\partial t}}=D{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial x}{\partial r}}\right)$$

(18)

that can be solved by separating the flux of composition J (19) and the conservation of mass (20)

$$J=-D{\displaystyle \frac{dx}{dr}}$$

(19)

$${\displaystyle \frac{dxV}{dt}}=-{\left(AJ\right)}_{out}+{\left(AJ\right)}_{in}$$

(20)

where x is composition, V is the volume of a shell, A its internal and external surface area.

Equations (16), (17), (19) and (20) are solved via a Finite Volume scheme with a MATLAB script.

3. Results

This study allows to discern about:

• overall motion of magma parcels;
• Lagrangian trajectories;
• time-evolution of the size of the crystals;
• time-evolution of the composition of the crystals.

Referring to the work by Longo et al. (2023), simulations where active stirring occurs are considered, i.e., P0, P3.5, P5 and P6.5. These involve pure buoyancy of the deep magma with respect to the shallow one (P0), or buoyancy and deep overpressure, where the deep magma is overpressurized ad pushes magma upwards (P3.5, P5, and P6.5). Table 1 reports the strength of buoyancy by the Archimede’s number and the strength of overpressure by the Hagen number. In overpressurized cases, buoyancy, which drives convection, decreases from P5 to P3.5 to P6.5, while pressure, which inhibits convection, increases.

In each of the chosen cases, about 2500 magma parcels are seeded in the shallow chamber, and 500 parcels in the upper conduit where convection is highly active. No particles are set in the deep chamber as no relevant motion is present there (Longo et al. 2023). The motion of each particle is followed, along with the pressure and the magma composition it transverses.

3.1. Overall Magma Parcels Location

Thanks to the Lagrangian computation, the overall motion of unitary magma parcels can be found. In particular, it is interesting to note the initial position of the particles in the shallow chamber that undergo recirculation inside the dike (Figure 2). It is evident that the majority of particles involved in this recirculation are located next to the lower and side walls of the shallow chamber. Instead, the particles marked as down-going, that are initially in the center of the shallow reservoir, are transported by internal convection towards the right wall and from there they descend into the dike (Videos S1-S4 in SI). Several parcels of magma after descending in the dike, rise upwards, and are subject to convection in the shallow chamber again.

On the contrary, all but a few magma parcels initially located in the dike rise upwards in the shallow reservoir and there are subject to local convection (Figure 3). Some of these particles descend again back into the dike and occupy the upper 25–50 m of the dike (P0 and P3.5). In P5 case, where dynamic is more active, such recirculation reaches depths of -3900 m (500 m down into the dike). Finally, simulation P6.5 does not present any of this motions. In simulation P0 the removal of particles from the dike and the start of recirculation downwards starts at about 4400 s. In P5 at 2950 s the downward motion initiates while still particles are being removed from the dike. The same for simulation P3.5, but at 2800 s. The complete drag of particles form the dike in simulation P6.5 is reached at 5760 s.

3.2. Lagrangian Trajectories

Lagrangian trajectories of magma parcels seeded in the shallow chamber are depicted in Figure 4. It is evident the strong recirculation within the shallow chamber and between it and the dike. Convective lines of motion are strongly entangled one with the other. As far as the buoyancy contribution decreases, the shallow chamber convection is more and more regular. Trajectories of the magma parcels rising from the dike are depicted in Figure 5. Their trajectories inside the shallow reservoir are similar to those of the parcels already residing in the shallow chamber, indicating good level of mixing.

The Finite-Time Lyapunov Exponent measures how fast nearby particles separate over a finite time interval. It would be a useful indicator, but it assumes exponential growth of the distance between particles. This can happen only in unbounded domains. In a restricted domain like the magmatic system it cannot be used. Therefore, a useful simple quantity to be investigated is the evolution in time of the distance between particles. To exemplify the divergence or proximity between magma parcels, five neighborhood particles are seeded in different initial locations inside the shallow chamber, and their respective distance is calculated. Depicted are the initial positions in the center of the reservoir (Figure 6), in the center of the left recirculating patterns (Figure 7), and along the bottom wall (Figure 8). The maximum distance reached between particles, along with the flow characteristics are reported in Table 2.

Going through the columns of Table 2 it is evident the dispersive behavior of simulation P5, while P0, P3.5 and P6.5 show nearly always coherent flow. Starting from several initial locations, there is detachment of one or two particles from the rest of the neighbor cluster. Simulation P6.5 distinguishes itself for its low buoyancy contribution and high overpressure, the upward flow of magma forms highly coherent cells on the left and right portions of the shallow chamber. Parcels that descend into the dike are always subject to chaotic motion.

Lagrangian trajectories and reciprocal distances between magma parcels initially located in the center of the dike 200 m deep (Figure 9) show that some parcels fall down into the dike. In these latter cases the reciprocal distances reach ∼1 km.

3.3. Parametric Study on Interface-Controlled Growth Rate and Crystals Volumetric Number Density

Algorithm (16-17) depends on the interface-controlled growth rate k, the crystal volumetric number density n, and the densities of crystals and magma, ${\rho }_c$ and $\rho$. The choice of the interface-controlled growth rate k and of the volumetric number density n is critical due to their uncertainty and wide range of variability. Therefore, a parametric study is performed on these two variables to state their effects on the crystal growth behavior.

The solid solution of feldspar is taken as reference, so that ${\rho }_c=2600$ kg/m³ (MELTS), while from simulations, $\rho =2400$ kg/m³.

According to Moschini et al. (2023), k varies between ${10}^{-12}$ and ${10}^{-6}$ m/s.

Furthermore, Schiavon et al. (2023) report a pre-eruptive chamber-resident crystal cargo of 7 – 35 mm⁻³, values around 10¹⁰ m⁻³ are also confimed by Salisbury et al. (2008); Suhendro et al. (2021).

The parametric study is thus conducted varying n in the range ${10}^8-{10}^{12}$ m⁻³ taking k = ${10}^{-8}$ m/s constant (Table 3 and Table 5), and on $k={10}^{-12}-{10}^{-6}$ m/s, fixing $n={10}^{10}$ m⁻³ (Table 4 and Table 6). Magma parcels initially residing next to the bottom walls of the chamber are considered at first (Table 3 and Table 4), lately, magma parcels initially residing in the dike are investigated (Table 5 and Table 6).

Table 3. Results of the parametric study on the volumetric number density n and $k={10}^{-8}$ m/s for magma parcels initially residing next to the bottom walls of the shallow chamber. Values of initial radius $R_{eq}$, minima and maxima of equilibrium proxy S and crystal weight fraction x.

Sim.	n	${\mathit{R}}_{\mathit{e}\mathit{q}}$ (m)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^8$	5.9419e-04	-1.8685	0.1640	0.0707	0.0952
	${10}^{10}$	1.2802e-04	-1.8682	0.4690	0.0448	0.0952
	${10}^{12}$	2.7580e-05	-1.8672	0.5729	0.0337	0.0952
P5	${10}^8$	5.9378e-04	-1.8641	0.2587	0.0609	0.0950
	${10}^{10}$	1.2793e-04	-1.8115	0.5318	0.0369	0.0950
	${10}^{12}$	2.7561e-05	-1.7175	0.5648	0.0330	0.0950
P3.5	${10}^8$	6.0159e-04	-1.9205	0.1369	0.0732	0.0988
	${10}^{10}$	1.2961e-04	-1.8309	0.4205	0.0451	0.0988
	${10}^{12}$	2.7923e-05	-1.7094	0.5640	0.0346	0.0988
P6.5	${10}^8$	6.7361e-04	-2.5993	0.2246	0.0950	0.1387
	${10}^{10}$	1.4512e-04	-2.0513	0.5404	0.0517	0.1387
	${10}^{12}$	3.1266e-05	-1.5446	0.5970	0.0378	0.1387

Table 3. Results of the parametric study on the volumetric number density n and $k={10}^{-8}$ m/s for magma parcels initially residing next to the bottom walls of the shallow chamber. Values of initial radius $R_{eq}$, minima and maxima of equilibrium proxy S and crystal weight fraction x.

Sim.	n	${\mathit{R}}_{\mathit{e}\mathit{q}}$ (m)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^8$	5.9419e-04	-1.8685	0.1640	0.0707	0.0952
	${10}^{10}$	1.2802e-04	-1.8682	0.4690	0.0448	0.0952
	${10}^{12}$	2.7580e-05	-1.8672	0.5729	0.0337	0.0952
P5	${10}^8$	5.9378e-04	-1.8641	0.2587	0.0609	0.0950
	${10}^{10}$	1.2793e-04	-1.8115	0.5318	0.0369	0.0950
	${10}^{12}$	2.7561e-05	-1.7175	0.5648	0.0330	0.0950
P3.5	${10}^8$	6.0159e-04	-1.9205	0.1369	0.0732	0.0988
	${10}^{10}$	1.2961e-04	-1.8309	0.4205	0.0451	0.0988
	${10}^{12}$	2.7923e-05	-1.7094	0.5640	0.0346	0.0988
P6.5	${10}^8$	6.7361e-04	-2.5993	0.2246	0.0950	0.1387
	${10}^{10}$	1.4512e-04	-2.0513	0.5404	0.0517	0.1387
	${10}^{12}$	3.1266e-05	-1.5446	0.5970	0.0378	0.1387

Table 4. Results of the parametric study on interface-controlled growth rate k and $n={10}^{10}$ m⁻³ for magma parcels initially residing next to the bottom walls of the shallow chamber. Values of minima and maxima of equilibrium proxy S and crystal weight fraction x.

Sim.	k (m/s)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^{-12}$	-1.8761	0	0.0952	0.0952
	${10}^{-8}$	-1.8682	0.4690	0.0448	0.0952
	${10}^{-6}$	-1.8434	0.6091	0.0331	0.0952
P5	${10}^{-12}$	-1.8875	0	0.0950	0.0950
	${10}^{-8}$	-1.8115	0.5318	0.0369	0.0950
	${10}^{-6}$	-1.3048	0.4256	0.0329	0.0950
P3.5	${10}^{-12}$	-1.9642	0	0.0988	0.0988
	${10}^{-8}$	-1.8309	0.4205	0.0451	0.0988
	${10}^{-6}$	-0.9057	0.3608	0.0333	0.0988
P6.5	${10}^{-12}$	-2.8770	0	0.1387	0.1387
	${10}^{-8}$	-2.0513	0.5404	0.0517	0.1387
	${10}^{-6}$	-0.8598	0.3358	0.0358	0.1387

Table 4. Results of the parametric study on interface-controlled growth rate k and $n={10}^{10}$ m⁻³ for magma parcels initially residing next to the bottom walls of the shallow chamber. Values of minima and maxima of equilibrium proxy S and crystal weight fraction x.

Sim.	k (m/s)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^{-12}$	-1.8761	0	0.0952	0.0952
	${10}^{-8}$	-1.8682	0.4690	0.0448	0.0952
	${10}^{-6}$	-1.8434	0.6091	0.0331	0.0952
P5	${10}^{-12}$	-1.8875	0	0.0950	0.0950
	${10}^{-8}$	-1.8115	0.5318	0.0369	0.0950
	${10}^{-6}$	-1.3048	0.4256	0.0329	0.0950
P3.5	${10}^{-12}$	-1.9642	0	0.0988	0.0988
	${10}^{-8}$	-1.8309	0.4205	0.0451	0.0988
	${10}^{-6}$	-0.9057	0.3608	0.0333	0.0988
P6.5	${10}^{-12}$	-2.8770	0	0.1387	0.1387
	${10}^{-8}$	-2.0513	0.5404	0.0517	0.1387
	${10}^{-6}$	-0.8598	0.3358	0.0358	0.1387

Table 5. Results of the parametric study on the volumetric number density n and $k=10-8$ m/s for magma parcels initially residing in the dike. Values of initial radius $R_{eq}$, minima and maxima of equilibrium proxy S and crystals weight fraction x.

Sim.	n	${\mathit{R}}_{\mathit{e}\mathit{q}}$ (m)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^8$	5.9419e-04	-0.0845	0.6497	0.0331	0.0403
	${10}^{10}$	1.2802e-04	-0.4048	0.6491	0.0332	0.0632
	${10}^{12}$	2.7580e-05	-1.3234	0.6462	0.0331	0.0922
P5	${10}^8$	5.9378e-04	-0.0900	0.6452	0.0331	0.0437
	${10}^{10}$	1.2793e-04	-0.3662	0.6390	0.0330	0.0736
	${10}^{12}$	2.7561e-05	-0.6348	0.6220	0.0329	0.0921
P3.5	${10}^8$	6.0159e-04	-0.0452	0.6600	0.0331	0.0413
	${10}^{10}$	1.2961e-04	-0.1879	0.6509	0.0331	0.0677
	${10}^{12}$	2.7923e-05	-0.5056	0.6230	0.0331	0.0952
P6.5	${10}^8$	6.7361e-04	-1.0022e-09	0.7518	0.0331	0.0437
	${10}^{10}$	1.4512e-04	-0.1322	0.7355	0.0331	0.0829
	${10}^{12}$	3.1266e-05	-0.4686	0.7024	0.0331	0.1346

Table 5. Results of the parametric study on the volumetric number density n and $k=10-8$ m/s for magma parcels initially residing in the dike. Values of initial radius $R_{eq}$, minima and maxima of equilibrium proxy S and crystals weight fraction x.

Sim.	n	${\mathit{R}}_{\mathit{e}\mathit{q}}$ (m)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^8$	5.9419e-04	-0.0845	0.6497	0.0331	0.0403
	${10}^{10}$	1.2802e-04	-0.4048	0.6491	0.0332	0.0632
	${10}^{12}$	2.7580e-05	-1.3234	0.6462	0.0331	0.0922
P5	${10}^8$	5.9378e-04	-0.0900	0.6452	0.0331	0.0437
	${10}^{10}$	1.2793e-04	-0.3662	0.6390	0.0330	0.0736
	${10}^{12}$	2.7561e-05	-0.6348	0.6220	0.0329	0.0921
P3.5	${10}^8$	6.0159e-04	-0.0452	0.6600	0.0331	0.0413
	${10}^{10}$	1.2961e-04	-0.1879	0.6509	0.0331	0.0677
	${10}^{12}$	2.7923e-05	-0.5056	0.6230	0.0331	0.0952
P6.5	${10}^8$	6.7361e-04	-1.0022e-09	0.7518	0.0331	0.0437
	${10}^{10}$	1.4512e-04	-0.1322	0.7355	0.0331	0.0829
	${10}^{12}$	3.1266e-05	-0.4686	0.7024	0.0331	0.1346

Table 6. Results of the parametric study on the interface-controlled growth rate k and $n={10}^{10}$ m⁻³ for magma parcels initially residing in the dike. Values of minima and maxima of equilibrium proxy S and crystals weight fraction x.

Sim.	k (m/s)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^{-12}$	-0.0060	0.6501	0.0333	0.0333
	${10}^{-8}$	-0.4048	0.6491	0.0332	0.0632
	${10}^{-6}$	-1.0993	0.6355	0.0331	0.0951
P5	${10}^{-12}$	-0.0061	0.6473	0.0331	0.0331
	${10}^{-8}$	-0.3662	0.6390	0.0330	0.0736
	${10}^{-6}$	-0.4962	0.4928	0.0329	0.0937
P3.5	${10}^{-12}$	0	0.6634	0.0331	0.0331
	${10}^{-8}$	1.2961e-04	0.6509	0.0331	0.0677
	${10}^{-6}$	-0.5110	0.4797	0.0331	0.0982
P6.5	${10}^{-12}$	-4.3730e-13	0.7585	0.0331	0.0331
	${10}^{-8}$	-0.1322	0.7355	0.0331	0.0829
	${10}^{-6}$	-0.3988	0.5275	0.0331	0.1368

Table 6. Results of the parametric study on the interface-controlled growth rate k and $n={10}^{10}$ m⁻³ for magma parcels initially residing in the dike. Values of minima and maxima of equilibrium proxy S and crystals weight fraction x.

Sim.	k (m/s)	min S	max S	min x (w.f.)	max x (w.f.)
P0	${10}^{-12}$	-0.0060	0.6501	0.0333	0.0333
	${10}^{-8}$	-0.4048	0.6491	0.0332	0.0632
	${10}^{-6}$	-1.0993	0.6355	0.0331	0.0951
P5	${10}^{-12}$	-0.0061	0.6473	0.0331	0.0331
	${10}^{-8}$	-0.3662	0.6390	0.0330	0.0736
	${10}^{-6}$	-0.4962	0.4928	0.0329	0.0937
P3.5	${10}^{-12}$	0	0.6634	0.0331	0.0331
	${10}^{-8}$	1.2961e-04	0.6509	0.0331	0.0677
	${10}^{-6}$	-0.5110	0.4797	0.0331	0.0982
P6.5	${10}^{-12}$	-4.3730e-13	0.7585	0.0331	0.0331
	${10}^{-8}$	-0.1322	0.7355	0.0331	0.0829
	${10}^{-6}$	-0.3988	0.5275	0.0331	0.1368

The initial radius $R_{eq}$, which stems from MELTS equilibrium mass of feldspars, falls inside the range of radii of microphenocrysts, $5·{10}^{-5}$–$2.5·{10}^{-4}$ m (Hansen Grönvold 2000,Salisbury et al. 2008,Wang et al. 2003,Zellmer 2021). This is in accordance with the hypothesis on the simulated system, which involves pre-eruptive resident crystals and not the much smaller sin-eruptive microlites.

The range of weight fraction of crystals is between ∼3% and ∼13%, in accordance with the hypothesis of dilute crystals.

The minimum crystals weight fraction decreases with increasing n and increasing k (Table 3 and Table 4). Dissolution is more intense for smaller particles (increasing n) and for greater interface-controlled growth rate k. On the contrary, maximum crystals weight fraction increases with increasing n and k, as growth is enhanced by smaller particles and greater k (Table 5 and Table 6).

The minimum and maximum value of S (Table 3 and Table 4) show that there is an accentuated tendency towards dissolution in magma parcels initially residing next to the bottom walls of the shallow chamber. Crystals evolve towards dissolution as soon as they are advected into the dike. The growth interval corresponds to their successive re-entrainment inside the shallow chamber.

The maximum crystal weight fractions for magma parcels initially residing in the shallow chamber remain equal to the initial equilibrium values, as the particles growth, rising in the shallow chamber after dissolution into the dike, up to a size lower then the initial radius.

The magma parcels that lie initially in the dike have an opposite behavior: crystals growth with no dissolution as soon as they rise inside the shallow chamber. The overall tendency is towards growth as maxima of S is greater than the absolute value of minima of S (Table 5 and Table 6).

Minimum weight fraction stays equal to the minimum initial x, while the maximum increases above the initial maximum of x. This is due to the net increase of radius of particles initially residing in the dike.

3.4. Kinetic Proxy

The crystal particles that undergo the maximum values of S manifest the maximum tendency to growth successive rims. On the contrary, those that are subject to the minimum values of S show stronger dissolution. The distribution in the shallow chamber of crystal particles and their respective maximum and minimum values of S over time are depicted in Figure 10 for the reference case of $k={10}^{-8}$ m/s and $n={10}^{10}$ m⁻³. Such particles represent a small subset of the particles inside ths shallow chamber, They are located next to the bottom walls and are the ones that approach the dike or descend into it (minimum S) and return back into the shallow chamber (maximum S). The same quantities, but for the dike, are shown in Figure 11. The particles having the maximum S fill the entire dike, as all or nearly all magma parcels in the dike rise upwards inside the shallow chamber and there increase in size. No dissolution develops in the dike.

3.5. Kinetic Proxy and Crystals Radius

Figure 12 shows the trajectories of reference magma parcels initially seeded along the bottom wall of the shallow chamber, along with the time-evolution of composition, pressure, albite fractions (albite is taken as the reference phase for the further computation of zoning), kinetic proxy and crystal radii for $k={10}^{-8}$ m/s and $n={10}^{10}$ m⁻³. All these magma parcels enter the dike to further rise in the shallow chamber again.

The kinetic proxy S is a complex function of composition and pressure, but, on the whole, it loosely resembles the trend of composition.

When the magma parcel descends into the dike composition increases from zero (shallow chamber resident magma: phonolite) to one (dike resident magma: shoshonite). The equilibrium mass fraction of feldspar in shoshonite is lower than the actual mass fraction beard by the phonolitic magma, with a lower amount of albite (MELTS). The kinetic proxy is negative, i.e., dissolution develops and the crystal radii decreases of about 30 $\mu$m. On the contrary, when the magma parcel rises back into the shallow chamber, composition decreases towards zero, a greater mass fraction of feldspar is at equilibrium (MELTS), the kinetic proxy becomes positive, the crystal composition is a little less albitic than the initial one and the radius increases. However, this increase in radius is such that the crystal size remains always below the initial value.

Trajectories, time-evolution of composition, pressure and albite phase for the dike are reported in Figure 13. Here the initial composition is shoshonitic, for an initial lower content of equilibrium feldspars and albitic composition (MELTS). As soon as the magma parcel enters the shallow chamber, the equilibrium feldspar weight fraction is greater (MELTS), the kinetic proxy is positive, composition in more albitic and the radius increases of ∼20 $\mu$m.

In P0 case the magma parcel is involved in recirculation inside the upper dike and its composition oscillates between one and zero; the kinetic proxy oscillates accordingly.

3.6. Parametric Study: Kinetic Proxy, Crystals Weight Fraction and Radius

Figure 14, Figure 15, Figure 16 and Figure 17 report the kinetic proxy, the weight fraction of crystals and crystal radii for the reference magma parcel next to the bottom walls of the shallow chamber in simulation P5 and P3.5 (see Figure 12). It is considered $k={10}^{-8}$ m/s and different values of crystal number density ($n={10}^8,{10}^{10},{10}^{12}$ m⁻³), (see Table 3, and Figure 14 and Figure 15). Also, it is fixed $n={10}^{10}$ m⁻³ and varied interface-controlled growth rate ($k={10}^{-12},{10}^{-8},{10}^{-6}$ m/s), (see Table 4, Figure 16 and Figure 17).

The effect of increasing the number of crystals n is to make a larger decrease in the crystals weight fraction x for a given negative kinetic proxy S, as can be seen around 1500 s for P5, and around 1000 s for P3.5. At the same time, there is a lower decrease in radius, but the larger amount of particles allows a larger decrease in weight fraction as observed above.

The lower crystal weight fraction gives rise to less negative kinetic proxy, as depicted in the interval 2000–5000 s for P5 and in 2000-3500 s for P3.5, simply because $S=(x_{eq}-x)/x_{eq}$. Given that for a negative S, $x_{eq}$ is lower than x, if the crystals weight fraction x is lower, the mass of actual crystals will be more next to the mass of equilibrium crystals, and S will be less negative.

Variations in k gives the most differences in S, x and r. Crystallization for $k={10}^{-12}$ m/s is nearly absent, while for $k={10}^{-6}$ m/s even for the small kinetic proxy there is a substantial change in x and radius.

3.7. Zoning

Solution of algorithm (16-17), and of Eqs. (20-19) allows to obtain the evolution in time of the crystals zoning by growth/dissolution and intracrystalline diffusion.

Figure 18 and Figure 19 report the final state of evolution for all simulations (Table 1) of reference magma parcels initially located next to the bottom walls of the shallow chamber for the range of solid-state diffusion coefficient of ${10}^{-15}-{10}^{-21}$ m²/s (Brady Cherniak 2010), $k={10}^{-8}$ m/s and $n={10}^{10}$ m⁻³.

The radius decreases during the whole interval of dissolution spent in the dike, and increases when the magma parcel rises into the shallow chamber. After the crystal has been dissolved, the first composition that precipitates is less albitic than the interior. The crystal was initially at equilibrium within phonolite inside the shallow chamber (higher Ab composition). After dissolution in the dike, the growth of new rims starts from a less Ab rich composition next the the exit of the dike, and occurs in a more and more phonolitic (Ab rich) environment as the crystals rise upwards. This remains almost constant during the successive advection inside the shallow chamber. The final albitic composition is still below the initial value as the shallow chamber now contains a mixture of phonolite and shoshonite.

The size of the dissolved rim is between 24 and 33 $\mu$m, while the accretionary rim is ∼ 6–8 $\mu$m, reached after 1500–1900 s of mixing inside the shallow chamber.

The case with the lower value of solid-state diffusion coefficient ($D={10}^{-21}$ m²/s) has pure interface-controlled mechanism. As D increases, it is clearly seen that the innermost accretionary layers increase in albitic composition by diffusion from the interior of the crystal. The size of the diffusive layer is ∼ 2–3 $\mu$m for $D={10}^{-15}$ m²/s. In P0 and P5 cases it is not only the accretionary rim that increases in albitic composition by diffusion from the interior, but also the interior is impoverished in albitic composition due to diffusion towards the exterior of the crystal.

The crystals zoning consists in a steep decrease in composition of 0.3–0.5 wt% of albite. With a subsequent slight increase.

The final zoning of the reference particles initially hosted by the dike for all simulations (Table 1) is reported in Figure 20 for the lowest and highest solid-state diffusion coefficients (10⁻²¹ m²/s and 10⁻¹⁵ m²/s). After no resorption, outer rims grow with more albitic fraction than the initially shoshonite albite-poor crystal in the dike when the crystal rises inside the shallow chamber.

The abrupt increase in albitic composition of crystals rising into the shallow chamber amounts to 1.7–2.5 wt% (exclunding P6.5) extending for 18–21 $\mu$m for an interval of mixing in the shallow chamber of ∼1–2 hours.

For the higher solid-state diffusion coefficient a diffusive layer of ∼5 $\mu$m develops. As for the magma parcels of the magma chamber, diffusion inside the interior of the crystals makes its composition to slightly increase.

3.8. Zoning in Parametric Study on Interface-Controlled Growth Rate

Variations of k give rise to sharp changes in the crystals radius (Figure 17 and Figure 16), which are also visible in the development of the zoning.

Zoning for magma parcels initially residing in the shallow chamber, for $k={10}^{-6}$ m/s, and solid-state diffusion coefficients of ${10}^{-21}-{10}^{-15}$ m²/s is depicted in Figure 21.

The shape of the zoning is largely different from the case of $k={10}^{-8}$ m/s. The successive growing rims are of size up to 2.5 $\mu$m, while for the lower k there is a continuous change in composition. The increase in composition after dissolution is gradual, while for the lower k it is sharp (Figure 18 and Figure 19). The increase in radius is about 40 $\mu$m. The diffusive layer is ∼5 $\mu$m. Diffusion acts all over the crystal. The internal composition decreases, the innermost accretionary layers increase in Ab fraction, while the external steps are smoothed.

The gradual $\mathit{vs}$ steep increase in composition can be explained on the basis of the radius and kinetic proxy change in time (Figure 16 and Figure 17). For $k={10}^{-6}$ m/s there is an abrupt increase in radius at 5000–6000 s in P5 and at 3500–4000 s in P3.5, followed by a very slow increase, while for $k={10}^{-8}$ m/s there is only a slow radius increase. Furthermore, the kinetic proxy becomes positive for $k={10}^{-6}$ m/s when it is still negative for $k={10}^{-8}$ m/s. Figure 22 compares different instant of crystallization for the two values of k for simulation P3.5.

At 3500 s, when for greater k ($k_g$) crystallization starts, for lower k ($k_l$) there is still dissolution in the dike. A soon as the crystal enters the shallow chamber ($t=3730$ s), $k_l$ starts crystallizing, while $k_g$ has already increased in radius of more then 10 $\mu$m with intermediate compositions. At 3970 s the innermost rim for $k_l$ is of higher composition as the crystal is already well inside the shallow chamber, and $k_g$ grown nearly to its final state. The kinetic proxy of $k_g$ decreases to zero before 4000 s, while that of $k_l$ remains positive (Figure 17). Therefore, until 5340 s there is only a little radius growth for $k_g$, with an increase of ∼10 $\mu$m for $k_l$. This latter slow growth is at almost constant composition as the crystal is advected inside the shallow chamber.

Zoning obtained for the reference magma parcels initially inside the dike is reported in Figure 23.

As for the magma parcels belonging to the shallow chamber, those inside the dike see a gradual increase in composition instead of a sharp one as in the case of $k={10}^{-8}$ m²/s. Again, the accretionary rims sum up to ∼40 $\mu$m.

4. Discussion

Three distinct crystal populations are identified. The majority of magma parcels initially hosted by the shallow chamber are involved in convective motion inside the shallow chamber itself and their composition is almost unchanged. On the contrary, magma parcels located next to the side and bottom walls are dragged into the dike where they experience dissolution, and are advected back into the shallow chamber where accretionary rims growth with the new composition of the chamber after mixing with the refilling magma. Finally, magma parcels initially hosted by the dike rise into the shallow chamber where accretionary layers of mixed composition soon develops over the refilling magma composition. The trajectories of all these magma parcels are strongly entangled, so that they tightly mix together.

Crystals seeded in the dike and crystals seeded at the bottom walls of the shallow chamber have accretionary rims over the same interval of composition: from that at the inlet of the dike to that inside the shallow chamber. However, the first do not present any reverse zoning with respect to the inner of the crystals, while the seconds show a reverse step between the core and the outer layers.

Counterintuitive sharp changes in composition appear for lower interface-controlled growth, rather then higher one. Given an interval of compositions to be covered by the growing radius, for a quick increase the steps of in radius increase are less steep in composition, while for a slower increase, the steps are steeper.

Chemical supersaturation can occur when two melts mix to produce a hybrid liquid whose equilibrium feldspar/plagioclase mode is higher than the actual crystal content. Existing feldspar then grows overgrowth rims until the melt relaxes toward equilibrium. Evidences of supersaturation during magma mixing can be found, for example, at Gran Canaria (TROLL SCHMINCKE 2002), at the roof of the Skaergaard intrusion (Naslund 1984), in the Karkonosze granite (SW Poland) (Słaby et al. 2008), and in the Bjerkreim–Sokndal intrusion (Norway) (Jensen et al. 2003). No data exist on crystal growth rate by supersaturation, but only for undercooling/decompression. These values range between 10⁻¹⁰ and 10⁻⁷ cm/s (i.e., 10⁻⁶ and 10⁻³ $\mu$m/s) (Billon 2025,Brugger Hammer 2010,Cashman 1993,Heckmann Strmić Palinkaš 2023,Moshrefzadeh et al. 2023,Mourey Shea 2019).

The growth rate of crystals among the whole parametric study is of 3·10⁻³–9.5·10⁻³ $\mu$m/s, on the higher range of what found in natural samples, even though the contrast in composition is not so sharp, but only between shoshonite and phonolite. The simulations performed in this study allow to shed light into the process of dissolution, while it is hard to measure in Nature and experiments the size of the dissolved layers. It is computed that these reach 25–40 $\mu$m at a rate of 7·10⁻³–18·10⁻³ $\mu$m/s, while experiments report a size of ∼2–3.5 $\mu$m (Mangler et al. 2024,Tsuchiyama 1985).

The simulated evolution represents an event of shallow chamber refilling. Extending the findings from one episode to repeated recharges, it can be speculated that a small population of crystals should record dissolution separating the original shallow chamber composition and the mixed one, while a substantial population should bear the mixed composition surrounding the original refilling one. The latter would be freshly produced at every repeated injection, while the former would record the variations in composition during the whole history of refilling.