Melting, compaction and reactive flow: Controls on melt fraction and composition change in crustal mush reservoirs

Abstract Changes in melt fraction and local bulk composition in high-crystallinity, crustal mush reservoirs are essential to produce the large volumes of low-crystallinity, silicic magma that are emplaced to form plutons, or erupted to surface. Heating (and cooling) is well understood and widely invoked in driving melt fraction change, but does not cause chemical differentiation because there is no separation of melt and crystals. Fractional crystallisation at high melt fraction is widely assumed to explain differentiation, but is inconsistent with the evidence that large-scale, long-term magma storage and evolution occurs in highcrystallinity mush reservoirs. Compaction has been suggested to explain melt fraction change and differentiation at low melt fraction, but compaction (and decompaction) causes simple unmixing (and mixing) of melt and solid crystals: to produce very refractory bulk composition by compaction, melt fraction must be driven down to very low values. Yet microstructural evidence demonstrating widespread compaction in crustal mush reservoirs at low melt fraction is lacking.

It is now widely, albeit not universally, accepted that long-term magma storage in 47 the continental crust occurs in low melt fraction (high-crystallinity) 'mush reservoirs' 48 rather than the high melt fraction (low-crystallinity) 'magma chambers' that have dom-49 inated conceptual models of magma storage and differentiation for over a century (e.g. (H-C) format is especially useful when discussing heat which, unlike temperature, is a 142 conserved quantity. We use the phase diagram in conceptual experiments to demon-143 strate the three processes considered here that impact melt fraction and local bulk 144 composition in crustal mush reservoirs: heating /cooling, compaction, and reactive flow.  (Figure 1c). Melting at constant temperature demonstrates that 173 the appropriate (conserved) quantity to describe melting is enthalpy, not temperature 174 (e.g. Figure 1b), as has been pointed out previously (Ussler and Glazner, 1992). How-175 ever, temperature is more easily measured and controlled in melting experiments. The 176 equilibrium melting process described here is reversible, such that cooling and crystalli- Cooling is an inevitable cause of melt fraction reduction in all igneous rocks, which 182 were once (partially) molten. Conversely, heating has been proposed as a mechanism 183 for melt fraction increase in crustal magma reservoirs, in a process termed 'thermal 184 rejuvenation' or 'defrosting' (Burgisser and Bergantz, 2011;Cooper and Kent, 2014; 185 Szymanowski et al., 2017). In this mechanism, a magma reservoir that is at low melt 186 fraction because it is 'cool' or 'cold' is heated by intrusion of new, hot magma, causing 187 the melt fraction in the surrounding mush to increase. 188 Thermal rejuvenation is widely suggested as a mechanism to create low-crystallinity 189 magmas in mush reservoirs, often over short timescales (e.g. Burgisser and Bergantz, In the context of our simple phase diagram, compaction can be visualised in terms of 216 a simple unmixing of a melt from its equilibrium solid residue; note we use mixing here 217 to mean mixing of melt and solid phases, not mixing at the molecular level. Consider 218 an experimental melting apparatus in which melt can be squeezed out of a partially 219 molten sample by operating a mechanical press (e.g. Philpotts and Dickson, 2000). 220 The heating system keeps the temperature constant. We place again 100 g of rock (an  The decrease in temperature caused by mixing means that, unlike compaction, the  We begin with the well known mass conservation (continuity) equation where φ, u l and u s are melt fraction, liquid (melt) and solid ( The equivalent expression for component conservation is given by where C l and C s are liquid and solid composition, respectively. Γ c is the rate of compo- The same procedure applied to equation 2(b) yields Adding equation 5 and 6 to eliminate Γ c leads to Re-arranging for Γ and substitution into equation 1(a), the rate of change of melt 395 fraction can be written as Applying a similar approach to bulk compositionC yields shows that the change in bulk composition can be expressed in terms of two different 400 contributions. We discuss the physical meaning of these terms in the next section.
At these points, melt fraction change occurs by melting or freezing at fixed temperature 405 and composition, and is a function only of latent heat exchange, given by where H is the enthalpy per unit mass and L f is the latent heat, which is assumed to where C ae is the eutectic composition (e.g. Figure 1), or in a binary solid solution when 410 C l = C s = 0 or 1.    and reactive flow to melt fraction change at a given time t are therefore given by where the subscripts H, C and R refer to melting/cooling, compaction and reactive flow respectively. The contributions to changes in local bulk composition are likewise given by Although equations (8) Table 2.

518
The experimental data demonstrate that evolved (high SiO 2 ) melt can be in equilibrium 519 with refractory (low SiO 2 ) crystals. Application of this fixed pressure phase diagram is 520 restricted to models of relatively limited vertical extent.  We next consider the same rock column which is sealed at the top and base and 572 has initial bulk composition corresponding to basalt withC = 51% SiO 2 (Fig. 3).

573
The column initially contains a uniform melt fraction (φ = 0.3; Fig (Table 2). This confirms that our nu-  We assume the mush is rigid. An analytical solution for this case can be found in As the more evolved melt injected into the base of the column flows upwards through 628 the column, it reacts with the solid phase so the local bulk composition at the top 629 of the column changes to become more evolved (Fig. 7i-l). The melt fraction also 630 increases at the top of the column, such that melt fraction no longer decreases linearly 631 with height ( Fig. 7b-d). The additional melt is created by 'chemical melting' rather  before the sill solidified (see Table 2). Using different mush material properties such as against cold country rock (Fig. 8t). As expected, heating/cooling has no effect on bulk 705 composition ( Fig. 8u-y).

878
Here the evidence of clinopyroxene formation by dissolution-reprecipitation reactions 879 with minerals in the mush (olivine + plagioclase) is more limited, probably because 880 the three phases are in a cotectic rather than peritectic relationship. In Figure 10F it 881 is apparent from the Ca X-ray map that plagioclase grains have sodic rims in contact 882 with clinopyroxene, suggesting that both minerals were precipitated during the reactive 883 flow process. This texture is not evident in plane-polarised light ( Figure 10E). In The binary phase diagram that we have adopted (Fig. 3) captures much of the im- can be written as where φ i and u i and Γ i are the melt fraction, velocity and the volume exchange rate of 1065 the ith phase, respectively. The component transport equation reads  i, respectively, we define The total change in melt fraction across all phases Φ can be decomposed into where the contributions from compaction, heating and reactive flow are given by the whereC j is the bulk composition of the jth component, which is defined as A.2. Including variable density 1082 We begin with the mass conservation (continuity) equation with variable density The source term Γ on the right hand side represents the rate of mass transfer from solid where C l and C s are liquid and solid composition, respectively. Γ c is the rate of com-1087 ponent mass transfer from solid to liquid. Applying the chain rule to equation 29(a), 1088 we have The same procedure applied to equation 29(b) yields Re-arranging for Γ and substitution into equation 28(a), the rate of change of melt 1092 fraction can be written as Thus an extra term − ∂ρ l ∂t appears when density changes are included. Since this term is 1094 independent of velocity, we suggest that density changes caused by thermal expansion 1095 or chemical reaction should be included in the heating/cooling term. Thus we have SiO 2 ) than for the liquid (glass) analysis itself (of order 0.2-0.8 wt% SiO 2 ), due to 1132 propagation of errors through all analysed phases and through the mass balance.

1133
The selected dataset of 42 experiments is given in Table 3, showing the range in predicted from our simple binary phase diagram lends further support to our approach.

1170
The good match to both fractional and equilibrium crystallisation experiments means 1171 that we can confidently use this phase diagram to explore not only simple crystallisa-1172 tion scenarios of specific starting compositions, but also reactive flow whereby the bulk 1173 composition is repeatedly modified by addition and loss of melt. 1174 We create a mathematical model that fits the experimental data for both phase 1175 diagrams ( Figs. 1 and 3), using a uniform solidus temperature, and a second order 1176 polynomial to describe the liquidus temperature, given by (Solano et al., 2014) 1177 where a k , b k , c k are fitting parameters, and k = 1 whenC <= C e and k = 2 when 1178C > C e .

1179
The anorthite-diopside binary (Fig. 1) exhibits two branches of the liquidus, one 1180 on each side of the eutectic composition C e , and we fit the data by adjusting the value 1181 of C e and the fitting parameters a, b and c (see Table 1 for the values used). The 1182 bulk compositionC varies from zero to one, with zero corresponding to An0 and one 1183 corresponding to An100.

1184
The binary phase diagram representing basalt (Fig. 3) where T s is the solidus temperature.
where φ, u, ρ, µ and ζ denote melt fraction, velocity, density, shear and bulk viscosity, where k t is the thermal conductivity.

1229
Using the relationship between enthalpy and temperature (equation 38), tempera-1230 ture can be eliminated from equation 46 and the energy equation becomes Conservation of components is given by where Fickian diffusion of components is small and is neglected here (Solano et al., We assume the viscosity of the solid phase varies with melt fraction, using a simple 1248 power-law relationship given by (e.g. Connolly and Podlachikov, 1998) 1249 where µ 0 is a reference viscosity and β is a constant that typically ranges between 0 1250 and 1 (Table 2).  (Table 2).

1255
The densities of the liquid and solid phases vary as a function of composition where ρ smin , ρ lmin , ρ smax and ρ lmax are the densities of the least and most evolved solid 1257 and liquid compositions, respectively (Table 2).
Permeability is here related to melt fraction (porosity) using the Kozeny-Carmen rela-1265 tionship, which can be simplified for small porosity to (Wark et al., 2003) 1266 where d is the grain size and A and B are adjustable constants. Coefficient c at low 1267 melt fraction is then given by At high melt fraction, we assume the separation of liquid and solid phases is described In this way, we scale the coupling between the end-members of porous media flow and 1274 hindered settling (Figure 11).
With zero velocity (no advection), the governing equations simplify to and the lever rule (equation 10). As there is no advection or diffusion of components, 1285 the bulk compositionC is constant. Substituting equations 58, 59, 10, into 60 leads to which can be solved analytically. With the insulating boundary conditions used in 1287 case I, melt fraction φ, assuming a uniform initial melt fraction φ 0 corresponding to a 1288 specified initial temperature, is given by where L is the length of the domain and The temperature T can be obtained from equation 62. with uniform initial melt fraction φ 0 . In dimensionless form, the solution is given by with boundary conditions: Equations (66) -(71) are non-dimensionalised using the compaction length and com-1298 paction time D.3. Analytical solution for reactive flow (Fig. 7) 1300 To validate the numerical method outlined in Appendix C when reactive flow domi-1301 nates (i.e. temperature remains fixed and the matrix is rigid so there is no compaction), which leads to where φ 0 is the initial melt fraction at that location and t is the elapsed time.