Thermodynamic Model for Energy-Constrained Open-System Evolution of Crustal Magma Bodies Undergoing Simultaneous Recharge, Assimilation and Crystallization: the Magma Chamber Simulator

The Magma Chamber Simulator quantifies the impact of simultaneous recharge, assimilation and crystallization through mass and enthalpy balance in a multicomponent–multiphase (melt + solids ± fluid) composite system. As a rigorous thermodynamic model, the Magma Chamber Simulator computes phase equilibria and geochemical evolution self-consistently in resident magma, recharge magma and wallrock, all of which are connected by specified thermodynamic boundaries, to model an evolving open-system magma body. In a simulation, magma cools from its liquidus temperature, and crystals ± fluid are incrementally fractionated to a separate cumulate reservoir. Enthalpy from cooling, crystallization, and possible magma recharge heats wallrock from its initial subsolidus temperature. Assimilation begins when a critical wallrock melt volume fraction (0·04–0·12) in a range consistent with the rheology of partially molten rock systems is achieved. The mass of melt above this limit is removed from the wallrock and homogenized with the magma body melt. New equilibrium states for magma and wallrock are calculated that reflect conservation of total mass, mass of each element and enthalpy. Magma cooling and crystallization, addition of recharge magma and anatectic melt to the magma body (where appropriate), and heating and partial melting of wallrock continue until magma and wallrock reach thermal equilibrium. For each simulation step, mass and energy balance and thermodynamic assessment of phase relations provide major and trace element concentrations, isotopic characteristics, masses, and thermal constraints for all phases (melt + solids ± fluid) in the composite system. Model input includes initial compositional, thermal and mass information relevant to each subsystem, as well as solid–melt and solid–fluid partition coefficients for all phases. Magma Chamber Simulator results of an assimilation–fractional crystallization (AFC) scenario in which dioritic wallrock at 0·1 GPa contaminates high-alumina basalt are compared with results in which no assimilation occurs [fractional crystallization only (FC-only)]. Key comparisons underscore the need for multicomponent–multiphase energy-constrained thermodynamic modeling of open systems, as follows. (1) Partial melting of dioritic wallrock yields cooler silicic melt that contaminates hotter magma. Magma responds by cooling, but a pulse of crystallization, possibly expected based on thermal arguments, does not occur because assimilation suppresses crystallization by modifying the topology of multicomponent phase saturation surfaces. As a consequence, contaminated magma composition and crystallizing solids are distinct compared with the FC-only case. (2) At similar stages of evolution, contaminated melt is more voluminous (∼3·5×) than melt formed by FC-only. (3) In AFC, some trace element concentrations are lower than their FC-only counterparts at the same stage of evolution. Elements that typically behave incompatibly in mafic and intermediate magmas (e.g. La, Nd, Ba) may not be ‘enriched’ by crustal contamination, and the most ‘crustal’ isotope signatures may not correlate with the highest concentrations of such elements. (4) The proportion of an element contributed by anatectic melt to resident magma is typically different for each element, and thus the extent of mass exchange between crust and magma should be quantified using total mass rather than the mass of a single element. Based on these sometimes unexpected results, it can be argued that progress in quantifying the origin and evolution of open magmatic systems and documenting how mantle-derived magmas and the crust interact rely not only on improvements in instrumentation and generation of larger datasets, but also on continued development of computational tools that couple thermodynamic assessment of phase equilibria in multicomponent systems with energy and mass conservation.

INTRODUCTION

Volcanic eruptions are one manifestation of the thermal and mass flux from Earth’s interior. They pose risks to population and property, and thus an important pursuit in Earth Science is improved prediction of the timing, magnitude, intensity and style of eruptions. It is equally important to understand the origin and evolution of magmas that ultimately solidify to generate new crust. Such intrusive bodies form oceanic and continental crust, are the reservoirs of many natural material and energy resources, transport heat from the mantle to the crust, and preserve a geochronological record of crustal and mantle evolution over Earth history (e.g. Condie et al., 2011; Sawyer et al., 2011).

Whether stored in the crust or erupted onto the surface, terrestrial magmas undergo variable processing en route to yield the physical and chemical diversity evident on Earth. Progress in understanding this diversity has relied upon fundamental developments in igneous petrology research that include marked improvements in high-precision, small spatial resolution geochemical analysis of rock components (e.g. Feldstein et al., 1994; Thirlwall & Walder, 1995; Davidson & Tepley, 1997; Reid et al., 1997; Bindeman et al., 2001; Ramos et al., 2005; Schmitt, 2011; and numerous others), laboratory studies and modeling of the thermodynamic and transport properties of melts and magmas (e.g. Spera, 2000; Stixrude & Lithgow-Bertelloni, 2010; Zhang & Cherniak, 2010), application of thermodynamic and fluid mechanical models to magmatic systems (e.g. Sigurdsson, 2000; Annen et al., 2002; Jellinek & DePaolo, 2003; Dufek & Bergantz, 2005; Fowler & Spera, 2010) and experimental phase equilibria [see Ghiorso & Sack (1995) and the Library of Experimental Phase Relations: lepr.ofm-research.org].

Among the most important discoveries made by the integration of this array of research tools is the prevalence of open-system magmatic behavior at a variety of scales (e.g. Sparks, 1986). Magma bodies in the crust and upper mantle are dynamic open systems that exchange heat, chemical components and momentum with their surroundings via a variety of processes [see Spera & Bohrson (2001) for a short summary]. A plethora of case studies of volcanic and plutonic rocks provides indisputable evidence that magma mixing (e.g. O’Hara & Mathews, 1981; Clynne, 1999; Murphy et al., 2000; Izbekov et al., 2004; Ginibre & Worner, 2007; Waight et al., 2007; Salisbury et al., 2008; Eichelberger et al., 2013) and crustal assimilation (e.g. Grove et al., 1988; Barnes et al., 2004; Wanless et al., 2010; Cebriá et al., 2011), along with finite-increment fractional crystallization, are the key processes contributing to magma compositional diversity.

Despite an abundance of data documenting recharge, assimilation and fractional crystallization (RAFC) processes, a continuing challenge involves utilizing the abundance of geochemical data to quantify and systematize these phenomena in multicomponent–multiphase open-system magma bodies. A step toward this goal is the energy-constrained (EC) modeling approach (Bohrson & Spera, 2001, 2003, 2007; Spera & Bohrson, 2001, 2002, 2004) that balances energy, trace elements and isotopes within open-system RAFC magma bodies. Limitations in the energy-constrained approach include the inability to track major element compositions and phase abundances. In contrast, MELTS and rhyolite-MELTS (Ghiorso & Sack, 1995; Asimow & Ghiorso, 1998; Ghiorso et al., 2002; Gualda et al., 2012) models yield phase equilibria and major element results relevant to crustal pressures but do not routinely incorporate trace elements and isotopes (although see capabilities of alpha-MELTS; Smith & Asimow, 2005), and are not currently configured to handle an arbitrary number of recharge events or direct coupling of magma sensible and latent heat for wallrock partial melting. In this study, the capabilities of EC-RAFC have been combined with those of rhyolite-MELTS to produce the Magma Chamber Simulator (MCS), a computational tool that rigorously tracks enthalpies, compositions and temperatures of melt, fluid and solids in a (resident) magma body–magma recharge–wallrock system undergoing simultaneous RAFC. MCS output includes masses and compositions (major and trace elements and isotopes) of melt, fluid and solids in wallrock and contaminated magma, as well as the thermal evolution of both, with allowance made for simultaneous heat and magma input by recharge. MCS results, which are thermodynamically self-consistent so that trace element and isotopic solutions are quantitatively coupled to relevant phase results, can be compared with whole-rock data as well as in situ analyses of crystals and melt inclusions for specific magmatic systems by specifying geologically relevant boundary and initial conditions.

As a thermodynamic tool, MCS can be utilized in two ways: the first is to generate forward models to illustrate the influence input parameters (e.g. mass of wallrock involved in AFC interaction, critical melt fraction that remains in wallrock, initial compositions and temperatures, etc.) have on evolving open-system magma bodies. The second use involves direct comparison of MCS models with high-resolution geochemical data and engagement of an iterative process by which a ‘best’ model of a specific natural system is derived. Both approaches provide an equilibrium reference state for open-system magma body behavior that can be used to explore the importance of model assumptions and limitations and to help elucidate the parameters that most influence the compositional, thermal and mass evolution of open systems. The MCS is not a dynamical model per se, although information from ancillary models, such as a heat transport model or a model of anatectic melt delivery by percolative flow, can be integrated with MCS results to obtain kinetic, dynamic, and timescale insight. MCS version 1.0 (v1.0), introduced in this work, represents the first stage of model development, and we anticipate that future code enhancements will respond to modeling needs that emerge as data from natural systems are compared with simulated results, as forward modeling defines the host of parameters that are critical to understanding open-system evolution and as constraints from other types of inquiry such as timescales of magmatic processes are explored.

OVERVIEW OF THE MAGMA CHAMBER SIMULATOR

The geochemical and petrological evolution of the composite system, which is composed of subsystems (1) magma body (resident magma melt ± crystals ± fluid), (2) cumulate reservoir, (3) wallrock, and (4) a set of recharge reservoirs, is linked to the temperature evolution of hotter magma and cooler wallrock as the two approach and eventually achieve thermal equilibrium. As the magma body cools, crystallizes and (possibly) receives additions of recharge magma (each of distinct bulk composition, temperature and phase state), it delivers its enthalpy to wallrock, which heats up, potentially to temperatures that exceed its solidus temperature. Wallrock partial melts may contaminate the magma body; both wallrock and magma compositionally evolve as a result of mass and enthalpy exchange and the state of both along the path to thermal equilibrium is determined by attainment of local (internal) chemical equilibrium subject to appropriate constraints (described below).

MCS retains all of the rigorous constraints imposed by rhyolite-MELTS phase equilibria and the EC-RAFC composite system formulation, allowing the computation of internal states of equilibrium in each subsystem. Subsystems can gain or lose energy or mass but total composite system energy and mass are strictly conserved. In MCS v1.0, once a major element, multiphase solution is found for the composite system, trace element and isotopic implications for all subsystems are explored via separate calculations based on the trace element and isotope mass balance approach used in EC-RAFC. A given MCS solution provides a continuous record of output as a function of magma temperature that includes the composition (major and trace elements, isotopic ratios), masses and temperatures of all relevant phases (melt + solids ± fluid) within the magma body, cumulate reservoir, recharge magma and wallrock. The utility of the model is emphasized by the range of quantitative results it produces. For example, magma crystal compositions that account simultaneously for RAFC are calculated, as is a record of the composition of melt potentially trapped within cumulates as melt (glass) inclusions. In addition, the compositional sequence and masses of anatectic melts generated by wallrock partial melting, as well as the compositions and masses of wallrock restite phases, are determined as a function of wallrock temperature. With such a broad array of computed quantities, MCS results can be compared directly with data from well-studied natural systems, providing evolutionary context for whole-rock, mineral and melt inclusion data.

The goal of this contribution is to introduce the Magma Chamber Simulator (v1.0). The conceptual model is presented together with details about initial and boundary conditions, thermodynamic potentials, output and model assumptions. Major and trace element, Sr and Nd isotope, mass and thermal results of a single AFC case in which high-alumina basalt intrudes upper crustal diorite illustrate the diversity of predicted quantities and highlight the utility of the tool as a framework for understanding the evolution of a particular open magmatic system. We discuss model limitations and briefly outline future enhancements to v1.0. This is followed by a discussion of heat transfer and partial melt mobility to provide a dynamic context for the MCS. We close by providing some logistical details about the code and information about how to access it.

MAGMA CHAMBER SIMULATOR: MODEL DESCRIPTION, INITIAL CONDITIONS AND COMPUTED RESULTS

Definitions of subsystems and initial conditions

The MCS defines a composite thermodynamic system that includes four subsystems: (1) a resident magma body that includes melt ± solids ± fluid; (2) cumulates ± separate fluid phase that form in equilibrium with magma melt and are fractionated incrementally (at each temperature step of imposed cooling) to a cumulate reservoir; (3) wallrock; (4) a set of recharge reservoirs. These subsystems are separated by boundaries described by classical thermodynamics as elucidated below. The composite system is modeled as isobaric [i.e. magma, cumulates, wallrock and recharge magma(s) interact at identical pressure], and the progress variable for the calculations is magma temperature, which does not necessarily monotonically decrease for a given simulation. For MCS v1.0, all melts in any given subsystem are assumed to be homogeneous, and enthalpy is assumed to be instantaneously delivered with its associated mass. Initial conditions and subsystem parameters required for an MCS simulation are listed in Table 1, and the path to thermodynamic equilibrium between magma body and wallrock is schematically illustrated in Fig. 1.

Fig. 1.

Schematic illustration of the four subsystems and thermal and mass exchanges among them in an RAFC event computed by the MCS. Steps that involve the use of rhyolite-MELTS to achieve a new equilibrium state are highlighted. Simulation stops when magma body and wallrock reach thermal equilibrium.

Table 1:

Initial conditions and subsystem parameters required for an MCS simulation

Pressure, P, for composite system
Ratio of initial WR mass to initial M mass, Λ = /
Critical melt fraction in WR,
Temperature decrement to subsystem M during approach towards thermal equilibrium of WR and M, ΔT
M subsystem melt temperature for jth recharge event, ⁠, etc.
Ratio of mass of jth recharge event to initial mass of M, Φ = /

Magma body, wallrock and recharge magma subsystem input for MCS simulation:
Subsystem	Initial bulk major oxide composition (for i oxide components)	Temperature	constraint	Mass	Initial trace element concentration, initial isotopic ratio1 (for x elements/isotope ratios)	Solid–melt, solid–fluid partition coefficients for each trace element
Magma body (M)		Initial temperature of subsystem	Buffer (e.g. QFM or Fe2+/Fe3+ initial ratio	Initial mass of subsystem (100% melt),	⁠,	⁠,
Wallrock (WR)		Initial temperature of subsystem	Buffer or Fe2+/Fe3+ initial ratio	Initial mass of subsystem,	⁠,	⁠,
Recharge, j events (Rj)			Buffer or Fe2+/Fe3+ initial ratio	Mass of jth recharge increment,	⁠,

¹It should be noted that WR requires a bulk isotopic ratio for the equilibrium isotope case. In the case of isotopic disequilibrium, each phase can be defined by a distinct initial isotope ratio. (See text for discussion.)

The resident magma body subsystem, denoted by M in what follows, is of specified starting bulk composition (⁠⁠) where i refers to the ith component (e.g. SiO₂, Al₂O₃, etc.), initial mass (⁠⁠), pressure (P), initial oxygen fugacity, and initial temperature (⁠⁠). M is initialized at its liquidus temperature, and thus is completely molten.

The second subsystem is solids ± fluid (collectively called cumulates, denoted C) that form in thermal and compositional equilibrium with M melt and are then transferred adiabatically to C in incremental batches, the size (mass) of which is governed by the temperature decrement imposed on M. (In any given simulation, the temperature decrement is set at the outset and kept constant. Extensive calculations show that temperature decrements in the range 0·2 to ∼20 K produce nearly identical results, except in certain situations, easily identified, as explained in Supplementary Data; Supplementary Data is available for downloading at http://www.petrology.oxfordjournals.org.) No initial conditions are required for C. Once solids ± fluid are put into C, no further chemical exchange is allowed with the melt of M; thermal exchange may or may not occur, as described below.

The third subsystem is the wallrock (WR) surrounding the magma body that thermally interacts with and transfers anatectic melt to M, provided certain conditions are met (see below). WR is initialized by specification of its initial mass (⁠⁠), bulk composition (⁠⁠), initial oxygen fugacity, and initial temperature (⁠⁠), which is always below its solidus. The ratio of initial WR mass to initial M melt mass, defined as ⁠, is constant during a single simulation and is set as an initial condition. Thus, Λ specifies the mass ratio of initially colder WR, ⁠, that interacts with ⁠; Λ will vary depending on the geological system of interest. For example, if one were interested in studying the effects of stoping, then Λ would be set to a small value to emulate complete ‘digestion’ or reaction of stoped blocks with M melt. In contrast, to study model systems where regional-scale metamorphism and assimilation of anatectic melt are important, Λ would be set to a larger value. Sensitivity studies elucidate how small vs large Λ affect the compositional, mass and thermal evolution of an AFC system, and in the case where MCS is used to model natural systems, estimates of Λ will derive from geological knowledge of the magmatic system in question. In fact, one can perform a variety of simulations with all initial conditions fixed except Λ to understand what the effective Λ might actually have been in a particular system. Because this ratio is widely variable, we have chosen to make Λ an initial condition rather than setting it a priori according to some specific heat transfer model. This maintains the applicability of MCS to a variety of systems.

The fourth subsystem is a set of recharge reservoirs (R_j). Addition of magma from R to M occurs via j distinct mass increments, ⁠, each of specified bulk composition, ⁠, for i oxide components, oxygen fugacity, and temperature, ⁠. Either melt or magma (i.e. melt + crystals + fluid) can be added. In either case, the constraint is that R is in internal (i.e. local) thermodynamic equilibrium at before addition to M. It should be noted that for each recharge event, the temperature of M (⁠⁠) must also be defined. The dimensionless parameter Φ ( = ⁠) is the mass ratio of the jth recharge mass increment to the initial mass of M and is the metric used to quantify the mass of a given recharge event. Φ is set as an initial condition of the simulation and, as for Λ, estimates will be based on knowledge of the system of interest.

For all subsystems, the oxygen chemical potential or constraint can be set to a specific buffer [e.g. quartz–fayalite–magnetite (QFM), nickel–nickel oxide (NNO), etc.], thus providing a fixed as a function of temperature at the pressure of the MCS calculation. Alternatively, the concentrations of FeO and Fe₂O₃ are specified in the subsystem, and is computed as part of the thermodynamic potential minimization or maximization, depending on the process of interest.

For the remainder of this work, M is used synonymously with (resident) magma body (or magma), C with cumulates, WR with wallrock, and R with the set of recharge reservoirs. The terms magma body melt or magma melt refer exclusively to the melt phase of the magma body. Cumulates represent crystals ± fluid that form in equilibrium with magma melt and are incrementally fractionated to C. Crystals (and fluid phase, if present) that remain as part of the WR assemblage are called restite, and the term anatectic melt is used solely in reference to the partial melt formed in wallrock, some of which may contaminate the magma body. A ‘step’ in the simulation reflects the new state of the composite system after some change has been imposed on one or more subsystems (e.g. temperature decrement imposed on magma body, crystals removed from magma body, anatectic melt removed from wallrock, enthalpy added to wallrock). Each simulation involves c. 50–500 temperature decrements, and the simulation ends when the magma body (M) and wallrock (WR) are in thermal equilibrium, which is dependent upon Λ.

Boundary conditions and thermodynamic potentials

M and WR are coupled by a diathermal and ‘osmotic’ (semi-permeable) boundary. The osmotic condition allows anatectic melt to enter the magma body when the fraction of melt in wallrock exceeds a critical threshold (described below). Once crystallization begins in M, crystals (± fluid) form in local equilibrium with melt and are removed to C. In v1.0, once in C, no further chemical interaction is permitted, and thus the boundary is closed. In terms of the energetics, two possibilities exist for cumulates: adiabatic or isothermal. In the adiabatic condition, cumulates do not stay in thermal equilibrium with M melt. That is, once formed, crystals (± fluid) are added to C, which is separated from M by a closed and adiabatic boundary. The alterative case is for crystals (± fluid) to exchange heat (but not mass), remaining in thermal equilibrium with M. Thus, C is closed but shares a diathermal boundary with M. In v1.0, C is closed and adiabatic. This assumption minimizes the amount of anatectic melt generated in WR because heat stored in C is not available for wallrock heating and partial melting. Future work will explore the impact of imposing a diathermal boundary and allowing cumulate crystals to remain in chemical equilibrium with melt.

The R subsystem is isolated from both M and WR except during a recharge event, when the boundary between M and R is open with respect to mass (including phase transfer) and energy (diathermal); thus, enthalpy and mass associated with the addition of recharge magma, itself in internal equilibrium before the mixing event, are added to M. After a recharge event, R is again isolated from M. WR interacts with R only indirectly; recharge enthalpy is added to M, and thus is available to WR through subsequent cooling and crystallization of M.

At each step in the MCS computation, when the state of equilibrium is disturbed in the M, WR or R subsystem during AFC or RAFC evolution, the new state of equilibrium is computed via minimization or maximization of an appropriate thermodynamic potential (Tisza, 1978). Two thermodynamic potentials are utilized in MCS.

1. For computational steps in which the temperature is known (pressure is always constant for a single simulation), the new equilibrium state is found via minimization of the Gibbs energy. For example, during cooling (e.g. M temperature decrement from 1200 to 1198°C) and crystallization within M, the Gibbs energy is minimized subject to constraints of fixed composition, temperature and pressure provided the oxygen fugacity is determined internally by ferrous/ferric redox. If, instead, oxygen fugacity is defined by a buffer, then the Khorzhinskii potential, defined as ⁠, is minimized. In either case, the potential minimization returns the mass and composition of each phase in M at the new temperature (1198°C).

2. For computational steps in which enthalpy has been exchanged, and the temperature of the subsystem is not known, entropy is maximized (Ghiorso & Kelemen, 1987). For example, as WR temperature increases owing to enthalpy transfer from cooling and crystallizing M, WR phase abundances and compositions change. To compute its new (local) equilibrium state, WR is evaluated by maximizing its entropy at fixed enthalpy and pressure.

The stoichiometric phases, liquid solution, fluid phase and solid solutions in the rhyolite-MELTS thermodynamic data/model have been documented by Ghiorso & Sack (1995; see additional references cited therein) and Gualda et al. (2012), and the solid phases treated by rhyolite-MELTS are listed in the Supplementary Data.

THERMODYNAMIC PATH TO THERMAL EQUILIBRIUM

Assimilation–fractional crystallization (AFC)

In MCS, cooling and crystallization in the magma body are accomplished by imposing a sequence of temperature decrements of constant value (e.g. ΔT = 2 or 5°C, etc.; Fig. 1, Step 1), beginning at the liquidus temperature of M and ending when magma and wallrock reach thermal equilibrium. At this stage, the simulation ends because the temperatures of wallrock and magma are equal, and thus there is no thermodynamic driving force to produce additional anatectic melt. (Initialization of M at a superliquidus temperature is possible and would enhance wallrock heating and the potential for partial melting. However, as a general condition, there is little evidence of superheated magmas in the crust and upper mantle.) Crystals ± fluid form in equilibrium with melt over the specified temperature decrement (e.g. from 1200 to 1198°C) (Fig. 1, Step 1). After melt + crystals ± fluid achieve a new equilibrium state, crystals ± fluid are fractionated to C (Fig. 1, Step 2). Because this process is not infinitesimal (i.e. not ‘perfect’ Rayleigh fractional crystallization), we define it here as incremental batch crystallization. The size of the incremental batch is governed by the chosen M temperature decrement; when it is small, crystallization very closely approaches perfect fractional crystallization, whereas if a larger temperature decrement is imposed, then crystallization approaches equilibrium crystallization. (For detailed discussion, see Supplementary Data.)

For each step in which M decreases in temperature and crystals form, the appropriate amount of enthalpy (i.e. 100% of sensible and latent heat associated with the mass of magma present at this temperature step and the mass of crystals that form in each magma temperature decrement step) is transferred through the diathermal M–WR boundary, and the multiphase–multicomponent phase equilibria of both M and WR and the local temperature of the wallrock are computed (Fig. 1, Step 1). Heat transfer continues as magma temperature decreases (and crystals ± fluid form; Fig. 1, Step 2), until at some point, anatectic melt is generated in WR (Fig. 1, Step 3). The attainment of this condition depends on the WR bulk composition, initial temperature, pressure, and the magnitude of sensible and latent heat generated in M owing to the imposed temperature decrement. We utilize the rhyolite-MELTS calibration to perform these WR melting calculations because we have confidence that phase equilibria computed on melting are modeled with reasonable accuracy; in contrast, rhyolite-MELTS should be used with caution (see Gualda et al., 2012) to examine phase relations coming down in temperature to near solidus conditions. The distinction between the two cases, which may at first appear contradictory, has to do in part with the accuracy with which the bulk composition of the system is known. If the bulk composition of the primary liquid is a bit off, as crystallization advances to low melt fractions, the liquid just does not partition correctly to the solid phases. Additionally, rhyolite-MELTS does not model the phase proportions well enough to do the solid–liquid partitioning correctly when the given bulk composition cannot be exactly expressed in terms of the more limited composition space of the solid phases in the model. On the other hand, if starting with a known rhyolite-MELTS phase assemblage below the solidus, then partially melting such an assemblage is a more straightforward problem that we have found generally yields useful and intuitive results.

When the wallrock solidus temperature is exceeded, the associated fraction of melt is computed by

(1)

where ⁠, ⁠, and represent the mass of melt, restitic solid and fluid in WR. Based on the rheological dynamics of melt segregation, the following simple criterion is adopted to specify at what temperature and how much anatectic melt flows through the osmotic boundary to contaminate the magma body melt: when the anatectic melt fraction exceeds a defined critical fraction in the wallrock (⁠⁠), then the mass of anatectic melt in excess of the amount represented by is removed from WR and added to M (Fig. 1, Step 4). By this action, the WR melt fraction after removal is precisely equal to the critical value, ⁠. Conceptually, the critical melt fraction represents a percolation threshold required for anatectic liquid mobility and hence transfer to M (see ‘Thermomechanical constraints’ discussion below). As an example, let us consider a case where is = 0·08. We consider that the total mass of WR is 1 kg and that 0·05 kg of this is anatectic melt. Obviously, in this state, the wallrock melt fraction is below the critical value of 0·08, and no anatectic melt is allowed to cross the WR–M boundary. After the next increment of heat is passed from M to WR, we assume the mass of anatectic melt increases to 0·1 kg. Because exceeds the critical fraction, anatectic melt removal is triggered. The mass of the anatectic melt, Ω, removed from WR and delivered to M is

(2)

where is the mass of anatectic melt in the wallrock subsystem at the current wallrock temperature (i.e. before any melt removal to M). In this particular case with = 0·08, a mass of anatectic melt (⁠⁠) of 0·1 kg, and a mass of WR restite (⁠⁠) of 0·9 kg, the mass of anatectic liquid added to the M subsystem is ∼0·0217 kg. The mass of melt that remains in the wallrock after the mobile portion is removed is 0·0783 kg, which is consistent with the ‘new’ wallrock melt fraction of 0·08. The assumption is made that anatectic melt is completely homogenized when it is brought into M; thus, MCS v1.0 does not allow for compositional gradients to form in M. It should be noted that if a fluid phase is present in WR, v1.0 assumes it remains in the wallrock restite. With the removal of anatectic melt, an updated initial condition for wallrock is defined by a new bulk composition, mass, and enthalpy.

Because addition of anatectic melt changes the state of the magma body, a new equilibrium condition is computed. This yields a new temperature and composition for M. Crystals ± fluid may form in response to assimilation; these form in equilibrium with contaminated melt and are then removed to C (Fig. 1, Step 5). Hence, the cumulate compositional record will convey information regarding the assimilation process. Once M has achieved its new equilibrium state and crystals ± fluid have been fractionated, the M temperature decrement is (again) imposed (Fig. 1, Step 1), generating another increment of heat that passes through the diathermal boundary into WR. In this fashion, the calculation is continued until M and WR reach the equilibration temperature.

Recharge–assimilation–fractional crystallization (RAFC)

In RAFC, allowance is made for the addition of recharge magma of specified mass (⁠⁠), bulk composition (⁠⁠) and temperature (⁠⁠) at any point in the evolution of M (Fig. 1, Step 3). To illustrate how recharge is handled, let us imagine a scenario in which M is undergoing cooling and incremental batch crystallization and WR, although heating up, remains below its solidus. M melt of specific composition is at some temperature ⁠. A recharge event now takes place. That is, a mass of recharge magma of bulk composition (⁠⁠) and temperature (⁠⁠) is added to M; pressure is that of M and WR. Recharge mass can be melt or melt ± crystals ± fluid. Because the recharge event has added mass of particular bulk composition and its enthalpy to M, the previous state of internal equilibrium in M is upset. A new state of equilibrium is achieved for this ‘recharged’ magma body. Because and at the instant of recharge are not generally identical, computation of the new equilibrium state defines a new magma body temperature, which may be above or below ⁠. At the end of these calculations, M possesses a new temperature, bulk composition, and mass of magma (melt ± solids ± fluid). Crystals ± fluid induced by recharge form in equilibrium with the ‘recharged’ magma body and are then removed to C. In this way, a record of recharge event(s) will be recorded in cumulate rocks. At the completion of the recharge event, another temperature decrement is imposed on M and the evolution continues with transfer of heat between (newly recharged) M and WR. There is no limit (except practical) on the number of recharge events or of the characteristics of recharge magma(s).

In summary, the following sequence of calculations models simultaneous RAFC. The equilibrium state of the magma in the recharge reservoir is evaluated at the local recharge temperature prior to the recharge event, and this state illuminates what is transferred to the magma body (i.e. melt ± crystals ± fluid). A recharge event occurs, M achieves a new ‘recharged’ equilibrium state and crystals ± fluid are fractionated to C. An amount of enthalpy appropriate to the new recharged state of the magma (including latent heat of crystallization that results from crystal formation in response to recharge) is transferred to WR, and the new equilibrium state of WR is assessed. If anatectic melt is transferred, M achieves a new ‘contaminated’ state of equilibrium and any crystals ± fluid that form in equilibrium with magma melt are removed to C by incremental batch crystallization. The wallrock equilibrium state is re-evaluated upon removal of anatectic melt. The relevant set of calculations continues for RAFC until M and WR are in thermal equilibrium.

MCS TRACE ELEMENT AND ISOTOPE MODELING

Overview

Once a major element/phase equilibria simulation has been computed (wherein M and WR achieve thermal equilibrium), isotopic ratios and trace element concentrations are calculated for wallrock and magma at each step. New input required to complete these calculations includes solid–melt and solid–fluid partition coefficients for each relevant phase and element, and initial whole-rock trace element concentrations and isotopic compositions for magma and wallrock (Table 1). For each distinct recharge event, only initial whole-rock trace element concentrations and isotopic values are required.

Solid–melt partition coefficients (K_sm = C_s/C_m where C_s is the concentration of element in solid and C_m is concentration of element in melt) are functions of composition and temperature at the fixed pressure of the MCS evolution. Either constant or variable partition coefficients can be employed in the MCS trace element calculation, depending on the range of magma and wallrock compositions and temperatures and on the availability of appropriate partition coefficient information. In many cases, the effects of temperature and phase (crystal, melt or fluid) composition on partition coefficients are inadequately known so a constant partition coefficient is assumed. For some treatments, it might be useful to calculate trace element partition coefficients based on crystal chemical models involving elastic strain in host crystals although this is not part of the MCS model per se (e.g. Wood & Blundy, 1997; van Westrenen et al., 1999; Law et al., 2000).

Some elements commonly utilized in trace element modeling exercises, including Cs, Rb, Ba, Pb, Sr, Mo, Th, Tb, Yb, U, Ce, Be and B, exhibit relatively low solid–fluid partition coefficients (K_sf = C_s/C_f where C_f is concentration of element in fluid) for common crustal and mantle phases; these elements are relatively soluble in hydrothermal supercritical fluids [see Spera et al. (2007) for a compilation of K_sf values]. Consequently, element partitioning among fluid, melt and solids is incorporated into the trace element mass balances (described below). Solid–fluid partition coefficients for each element of interest are therefore required as input.

For all trace elements in MCS, ⁠, ⁠, ⁠, and ⁠, where the overbar represents the bulk partition coefficient for each element of interest, are calculated for each step where crystals form in M and are then incrementally removed to C and for each step where anatectic melt forms, is removed from WR and is added to M. The bulk partition coefficient calculations require the solid–melt and solid–fluid partition coefficients as noted above and the proportion of solid and fluid phases in magma and wallrock, as output by MCS phase equilibria results. Thus, calculation of bulk partition coefficients tracks not only the changing solid ± fluid phase assemblage, but also changes in the solid–melt and solid–fluid partitioning behavior as the magma body cools and the wallrock heats up.

In the section below, we first treat trace elements and isotopes in the AFC case and follow that with the modifications needed to handle RAFC. We discuss two possibilities for the isotopic state of wallrock: equilibrium and disequilibrium. Isotopic equilibrium assumes that all phases in WR are in isotopic equilibrium, and thus at all stages of partial melting, anatectic melt isotopic ratios are constant. In contrast, (radiogenic) isotopic disequilibrium should be considered when wallrock exhibits a significant non-zero geological age. Because the MCS gives a self-consistent thermodynamic picture of partial melting in WR, the contribution that each phase with its distinct isotopic composition makes to the isotopic signature of partial melt can be determined.

Assimilation–fractional crystallization

For AFC, in M, each trace element has an initial concentration, ⁠. Because the simulation begins at the liquidus, the initial composition is also the bulk composition of M (i.e. M is 100% melt). Once the simulation begins but before anatectic melt crosses the boundary between WR and M, trace element concentrations in all of the phases of M are determined by the equilibrium crystallization relations [see AIV-2, 3, 4 of Spera et al. (2007) for the expressions used in MCS]. It should be noted that in rhyolite-MELTS, certain accessory phases such as zircon and allanite are not treated because of the lack of thermodynamic data. Although these phases are not directly incorporated into a trace element assessment using the phase data from MCS, consideration of how these phases affect trace element concentrations can be made using, for example, zircon saturation estimates from Watson & Harrison (1983) and incorporation of these into the relevant trace element calculation.

Once the trace element mass balance among crystals, melt and fluid in M is known, the mass of trace element removed via incremental batch crystallization (including fluid) is calculated algebraically. For each subsequent M temperature decrement, these calculations are repeated.

WR has an initial bulk trace element concentration, ⁠. For each WR temperature, the concentration of trace elements in coexisting phases is found from the equilibrium melting equations [see AIV-6, 7, 8 of Spera et al. (2007)], which are the same as the equilibrium crystallization equations noted above, with the exception that the initial conditions allow for the presence of a separate fluid phase at the solidus. We recall that if such a fluid phase does exist in the wallrock, it remains as part of the WR restite when anatectic melt is extracted and added to M. (Of course, any H₂O dissolved in the anatectic melt travels with it.) Once is exceeded, assimilation of a portion of WR anatectic melt occurs. The transfer of anatectic liquid may or may not lead to crystallization and/or formation of a separate fluid phase in M. In the case where it does not, the concentration of each trace element in contaminated magma melt is a function of the mass of trace element added by anatectic melt, the mass of trace element in the magma body melt just before the contamination step, and the (total) mass of contaminated melt in M (which reflects addition of a specified mass of anatectic melt). In the case where addition of anatectic liquid leads to an increment of crystallization or formation of a separate fluid phase in M, the calculation is slightly more involved. First, the new trace element concentration in M melt is found as just described. Because crystals or a fluid phase form in response, the distribution of the trace element among these newly precipitated solids, melt and fluid (if present) in M is calculated utilizing the equilibrium crystallization equations. The relevant mass of the trace element in crystals ± fluid is then subtracted from M because these phases are fractionated to C. At this point, the melt phase in M has a trace element composition reflecting contamination by anatectic melt and removal of associated crystals ± fluid. This procedure is followed for every M temperature decrement until WR and M achieve thermal equilibrium.

Isotopic ratios, radiogenic as well as stable, in WR and contaminated M (i.e. after coupling) may be calculated. Required input includes initial isotope ratios for M (⁠⁠) and WR (⁠⁠) (Table 1). For most radiogenic isotope systems, radiogenic ingrowth is negligible during relatively short-lived (roughly ⋚1 Myr) MCS evolution of M and is neglected in the calculations; the most extreme situations of parent element enrichment (e.g. extremely Rb-rich phases) or cases in which the half-life is short enough to yield ingrowth are not considered in v1.0. For WR, two general cases for radiogenic isotopes are possible. In the case of isotopic equilibrium, all phases are in internal isotopic equilibrium (e.g. identical ⁸⁷Sr/⁸⁶Sr). The isotopic composition of anatectic melt added to the M subsystem is therefore identical to ⁠. The isotopic composition of contaminated M melt is computed by isotopic mass balance using the concentration and isotopic ratio of M melt before the current assimilation increment and those of the relevant anatectic melt. Obviously, all crystals that form from M melt in response to assimilation are in isotopic equilibrium with that melt. In the case of isotopic disequilibrium, each WR phase possesses a distinct isotopic composition as an initial condition. The contribution each phase (including residual melt in WR) makes to the anatectic melt can be calculated from MCS phase equilibria results. With concentration and the isotope ratio information for each phase, the isotope mixing equation can be employed to determine the bulk anatectic melt element concentration and isotopic ratio. Mixing between this anatectic melt (not in isotopic equilibrium with its associated restite) with M melt is treated as in the isotopic equilibrium case described above.

Recharge–assimilation–fractional crystallization

To calculate trace element changes in the magma body owing to each episode of recharge, the bulk trace element concentrations for each recharge magma, ⁠, are required, as are the initial (bulk) isotope ratios, ⁠. For each distinct recharge event, melt and associated crystals are mixed with M melt, and the new trace element concentration for the ‘recharged’ or hybridized magma is calculated using the mass of trace element in melt of M (prior to recharge), mass of trace element added by R, and the (new) mass of M (which has changed as a result of addition of recharge melt ± crystals ± fluid). A new equilibrium state for the hybridized bulk composition and temperature of M is determined by MCS, and the equilibrium crystallization concentration mass balance expressions provide the distribution of trace element among melt ± solids ± fluid. The mass of each trace element is then debited from hybridized M according to the appropriate mass of solid(s) ± fluid removed. Isotopic mass balance is subject to the same limitations as addition of anatectic melt (i.e. no ingrowth is allowed during an MCS event). No radiogenic isotope disequilibrium among melt ± solid(s) ± fluid is permitted in R, so for each recharge event, each isotopic ratio is constant (e.g. a single ⁸⁷Sr/⁸⁶Sr or ¹⁴³Nd/¹⁴⁴Nd is defined for each recharge event).

RECAP OF MODEL ASSUMPTIONS AND FEATURES OF MCS v1.0

The Magma Chamber Simulator represents the merging of two extant models for treating compositional diversity in magma bodies: rhyolite-MELTS and EC-RAFC. A number of model features and assumptions have been described that pertain to MCS v1.0. We recognize that some of these assumptions are simplifications of how an open-system magma body behaves, but the current features of the code lend themselves to verification and exploration that will allow us to conduct forward modeling exercises that examine the importance of particular initial conditions and subsystem variables (e.g. ⁠, Λ, Φ), thoughtfully augment MCS in response to these modeling results, and model natural systems. Below we recapitulate the critical assumptions of MCS v1.0.

MCS is a zero-dimensional model; it therefore does not track time nor does it consider the complexity associated with thermal gradients in the wallrock. Total enthalpy and mass for the composite system are conserved, and thus an adiabatic boundary surrounds the composite system. Enthalpy exchange between subsystems is allowed as prescribed above (e.g. diathermal, adiabatic, etc.), and enthalpy is assumed to be instantaneously delivered with its associated mass. In MCS v1.0, complete homogenization of different magmas or melts (e.g. anatectic and magma melt, recharge magma and magma melt) is assumed, precluding development of chemical gradients. When crystals ± fluid phase form in equilibrium with melt in the magma body, they are incrementally removed to the cumulate reservoir (C); once removed, they do not interact chemically or thermally with the magma body. Thus, in v1.0, zoned crystals do not develop, but instead, the cumulate pile represents the range of open-system magma crystal compositions that form in equilibrium with M melt at distinct temperatures. Before anatectic melt transfer occurs, a critical melt fraction must develop in the wallrock, and transfer is based on the concept of a percolation threshold. In the case where anatectic melt forms in the wallrock and there is a separate fluid phase, MCS v1.0 transfers only the melt phase to the magma body; the fluid phase remains as part of the restite. Prior to addition, recharge magma is in internal thermodynamic equilibrium at its temperature and may contain melt ± crystals ± fluid. Upon addition, recharge crystals can be resorbed. Trace element and isotopic analysis requires solid–melt and solid–fluid partition coefficients to be defined for each element and phase, and different K_sm and K_sf can be employed in response to changes in M and WR. Radiogenic isotopic equilibrium between anatectic melt and wallrock restite is assumed in v1.0; thus, regardless of any isotopic heterogeneity in wallrock mineral phases, anatectic melt isotope ratios remain constant during a simulation. The simulation progress variable is magma temperature, and the simulation ceases once wallrock and magma are in thermal equilibrium. Although temperature is used as a progress variable, it does not necessarily monotonically decrease. For example, if hot recharge magma is added to cooler resident melt in M, the temperature after mixing could be higher in M than just before mixing.

MAGMA CHAMBER SIMULATOR AFC RESULTS: HIGH-ALUMINA BASALT INTRUDES UPPER CRUSTAL DIORITIC WALLROCK

In this section, we present results for a single AFC case to illustrate the range of predictions generated by MCS and to specifically highlight the behavior of a typical AFC scenario. We focus on intrusion of high-alumina basaltic melt into dioritic wallrock in the upper crust. Our intent is to illustrate the kinds of information available from an MCS solution, and we emphasize that this case in not intended to test specific hypotheses or model a particular magmatic system. Detailed analysis of the effects of particular parameters, such as the critical wallrock melt fraction, pressure, or initial wallrock and magma body temperatures and compositions, cases with significant episodes of recharge and application to specific magmatic systems (e.g. layered intrusions and volcanic successions) will be presented elsewhere.

AFC case overview

We present an AFC case in which high-alumina basaltic magma intrudes dioritic wallrock at 0·1 GPa (∼3 km depth). These compositions and the pressure were chosen to represent a typical magma–wallrock interaction in a continental arc. Initial major and trace element and radiogenic isotopic data for magma (Brophy & Marsh, 1986; Wilson, 1989) and wallrock (Wilson, 1989; Rudnick & Gao, 2004) and other input parameters are given in Table 2; selected results are included in Supplementary Data. Trace elements included in this analysis, Sr, Rb, Ba, La, Nd, Yb, and Ni, reflect a range of elemental behavior (large ion lithophile, rare earth, transition metal). and are presented in Supplementary Data and were selected based on a survey of values available at http://earthref.org/KDD/. Rarely, there were no data in the database (e.g. quartz); these were estimated to be 0·01. and were set to 50 for all elements. This implies that a negligible amount of solute is dissolved in the coexisting fluid phase when melt becomes volatile (H₂O) saturated. Although this value does not accurately reflect the correct solid–fluid partitioning for some of the modeled elements, its use simplifies analysis of the impact of AFC processes as determined by the MCS algorithm. Because rhyolite-MELTS does not currently include any volatile components other than H₂O, mixed volatiles in the system H–O–C–S–Cl are not modeled in the MCS v1.0. The impact of solid–fluid partitioning based on published partition coefficients and cases that involve mixed volatiles (e.g. H₂O–CO₂) will be analyzed elsewhere.

Table 2:

Major and trace element and Sr and Nd isotopic characteristics of magma and wallrock and selected input parameters for AFC and FC-only cases

	Magma (high-alumina basalt)	Wallrock (diorite)
SiO2	51·44	58·74
TiO2	0·60	0·77
Al2O3	16·31	17·80
Fe2O3	1·24	0·07
FeO	6·78	5·29
MgO	10·74	2·64
CaO	9·63	5·89
Na2O	2·23	4·02
K2O	0·42	2·38
P2O5	0·11	0·42
H2O	0·50	1·98
Rb	17	49
Sr	518	320
Ba	364	456
Ni	97	59
La	23	20
Nd	18	20
Yb	8	1·9
87Sr/86Sr	0·7035	0·7200
143Nd/144Nd	0·5131	0·5120

	AFC: Magma	Wallrock	FC: Magma
Pressure (GPa)	0·1	0·1	0·1
Magma temperature
decrement, ΔT (°C)	2		2
Initial wallrock T (°C)		500
		0·04

Trace element concentrations taken from Brophy & Marsh (1986) for magma and Rudnick & Gao (2004) for wallrock. Isotope ratios estimated based on knowledge of typical mantle and crustal values (e.g. Wilson, 1989). All major oxides are reported in wt %, and all trace elements are reported in ppm.

The simulation begins at the magma liquidus of ∼1279°C; the WR initial and solidus temperatures are 500°C and ∼765°C, respectively. The simulation ends at ∼945°C (Fig. 2a), when the temperature of WR and M are equal. The value Λ is unity, and the magma temperature decrement is 2°C. Discussion of results is divided into four stages that reflect changing conditions in the composite system, as follows.

1. Stage 1 (M temperature between ∼1279 and 1187°C): the WR subsystem is thermally but not chemically interacting with the magma body, and olivine is the sole precipitating phase in M.

2. Stage 2 (M temperature between ∼1187 and 1145°C): the WR subsystem is still heating up, reaches, and then surpasses its solidus temperature. However, the wallrock melt fraction remains below the critical value = 0·04; thus, no anatectic melt is transferred to M. In the magma, olivine, plagioclase, orthopyroxene (opx), high-Ca clinopyroxene (high-Ca cpx) and low Ca-pyroxene (low-Ca cpx) crystallize over different temperature intervals.

3. Stage 3 (M temperature between ∼1145 and 1066°C; WR temperature between ∼770 and 772·3°C): mass transfer between wallrock and magma occurs because WR contains a melt fraction above = 0·04. Wallrock partially melts over (only) an ∼2·3°C interval, although the cumulative per cent of anatectic melt added to the magma body is 18% (per cent is relative to initial magma body mass as described below). The very small temperature interval for stage 3 partial melting is due to pseudoeutectic melting that is dominated by alkali feldspar, plagioclase, quartz, H₂O, and apatite in dioritic wallrock.

4. Stage 4 (M temperature between ∼1066 and 945°C; WR temperature between ∼772·3 and 945°C): assimilation continues owing to mass and thermal exchange between M and WR. Over ∼167°C temperature drop in M, an additional 25% anatectic melt is added to the magma body.

Fig. 2.

Comparison of thermal and mass results for magma (blue) and anatectic melt (red) for AFC case and magma (gray) for comparative FC-only case. Stages of AFC case are labeled and identified on the AFC magma melt curve by tick marks. (For description of stage characteristics, see text.) (a) Per cent melt in magma body vs magma temperature (°C). (Per cent is relative to initial magma body mass.) Tick marks labeled with initial and final magma temperatures for each stage. Wallrock and magma are in thermal equilibrium at ∼945°C (end of stage 4). (b) Wallrock temperature vs magma temperature for stages 3 and 4. Final magma temperatures for stages shown. (c) Cumulative per cent cumulate minerals vs magma temperature. Each trend represents a mineral assemblage, starting with the liquidus phase, olivine, followed by olivine + plagioclase, etc. Thick lines are AFC case; thin lines are FC-only. It should be noted that the AFC and FC-only trajectories for olivine are identical. (For order of crystallization and distinctions between AFC and FC-only cumulate assemblages, see text.) Ol, olivine; Plag, plagioclase; Opx, orthopyroxene; Pyro, high- and low-Ca clinopyroxene; Sp, spinel. (d) Per cent cumulative anatectic melt added to magma body for stages 3 and 4. (e) Per cent plagioclase restite (in wallrock) vs decreasing magma temperature. (f) Per cent orthopyroxene (Opx), alkali feldspar, quartz, high-Ca clinopyroxene (Cpx), H₂O (phase), rhombohedral oxides (Rhm-oxides), apatite, and olivine that are in wallrock restite vs magma temperature. It should be noted that quartz and alkali feldspar are completely melted out of wallrock at magma temperatures of ∼1066 and 1022°C, respectively. Olivine joins the wallrock residual assemblage at a magma temperature of ∼1022°C.

Below, we present selected details of the phase equilibria, major and trace element, isotopic, and thermal evolution of these stages for M, WR and C (Tables 2–4, Figs 2–4). Because this case involves significant assimilation of wallrock partial melt into the magma body, a comparison case in which no assimilation occurs (i.e. fractional crystallization only, labeled as FC-only) isolates the distinctive effects of assimilation. For these two cases, all relevant parameters are the same. Conditions that lead to FC-only include the wallrock being of large mass (e.g. Λ = 10), having a low initial temperature and/or setting the wallrock critical melt fraction to a large number (e.g. 0·8).

Fig. 3.

Major oxide (wt %) vs magma temperature (°C) trends for magma body melt and anatectic melt in AFC case and magma melt in FC-only case. Stage labels, ticks and trend colors are the same as in Fig. 2. (a) SiO₂; (b) MgO with cumulate (blue) and wallrock (red) olivine Fo content and high-Ca clinopyroxene (cpx) and orthopyroxene (opx) Mg# annotated (note that the Mg# values for opx and cpx are labeled in stage 2 at the magma temperatures at which they first precipitate); (c) Al₂O₃; (d) CaO with cumulate plagioclase An content for AFC (blue) case annotated. The plagioclase-absent temperature interval for the AFC case should be noted. Wallrock plagioclase Ab content and alkali feldspar Or content (red) annotated. It should be noted that alkali feldspar is no longer a part of the wallrock residual assemblage at magma temperature of ∼1022°C (AF out). (See text for discussion.) (e) FeO; (f) TiO₂; (g) K₂O; (h) Na₂O; (i) P₂O₅; (j) H₂O (component).

Fig. 4.

Trace element (ppm) and Sr and Nd isotope ratios vs magma temperature (°C) trends for magma body melt and anatectic melt in AFC case and magma melt in FC-only case. Stage labels, ticks and trend colors are the same as in Fig. 3. (a) Ni; (b) Sr; (c) ⁸⁷Sr/⁸⁶Sr; (d) Rb; (e) Ba; (f) La; (g) Yb; (h) Nd; (i) ¹⁴³Nd/¹⁴⁴Nd.

Table 3:

Selected output for AFC and FC-only cases

	AFC: Magma	AFC: Wallrock	FC: Magma
Magma liquidus (initial) temperature (T)1	1279		1279
Coupling T	∼1145	∼770
Equilibration T	∼945	∼945
Olivine precipitation T2	1277		1277
Plagioclase precipitation T	1187		1187
Orthopyroxene precipitation T	1169		1169
High-Ca clinopyroxene precipitation T	1165		1165
Low-Ca clinopyroxene precipitation T	1157		1157
Spinel precipitation T	1048		1081
Cumulative per cent cumulates (magma), anatectic melt transferred to M at equilibration T3	72	43	81
Total per cent melt in magma at equilibration T3	72		19
Per cent residual wallrock solids (restite) + melt remaining at equilibration T3		57

¹All temperatures reported in °C. Liquidus temperature is determined by rhyolite-MELTS.

²Precipitation T refers to the temperature at which mineral first appears in magma body.

³Mass per cent calculated using initial magma body mass as denominator.

Table 4:

Selected output for stages of AFC case

	Stage 1	Stage 2	Stage 3	Stage 4	At equilibration T
Magma temperature (°C) range for stage1	1279–1187	1187–1145	1145–1066	1066–945	945
Magma melt mass (%) range1	100–91	91–60	60–69	69–71	71
Magma melt SiO2 (wt %) range1	51·5–52·5	52·5–54·2	54·2–59·4	59·4–64·4	64·4
Magma melt 87Sr/86Sr range1	0·7035	0·7035	0·7035–0·7053	0·7053–0·7084	0·7084
Magma melt 143Nd/144Nd range1	0·5130	0·5130	0·5130	0·5130–0·5125	0·5125
Cumulative per cent of crystals fractionated to C range1	0–9	9–40	40–49	49–72	72
Magma body cumulate mineral assemblage throughout stage2	olivine	ol + pl + opx + high-Ca cpx + low-Ca cpx	pl3 + high-Ca cpx + low-Ca cpx	pl + opx4 + high-Ca cpx + low-Ca cpx + spinel	pl + opx + spinel
Magma olivine Fo range1	88–84	84–83	–	–	–
Magma plagioclase An range1	–	80–75	75–723	72–53	53
Magma Opx Mg#	–	74–724		52 at 1013°C4	42
Magma high-Ca Cpx En–Fs–Wo5	–	53–10–37	51–11–38	45–15–40
Magma high-Ca Cpx Mg# range1	–	70–66	66–58	58–47
Wallrock temperature (°C) range1	500		770–772·3	772·3–945	945
Cumulative mass of anatectic melt added to magma body (%) range1	0	0	0–18	18–43	43
Anatectic melt SiO2 (wt %) range1	–	–	72	72–61	61
Residual wallrock solid + fluid assemblage6	–	–	pl opx ≈ af > q7 > high-Ca cpx > H2O > rhb-ox > apa	pl opx > af8 > high-Ca cpx > rhb-ox ≈ H2O > apa, (ol)9	pl > opx > ol > rhb-ox > high-Ca cpx ≈ H2O > apa
WR plagioclase An, Ab range1	–	–	38·5–37, 56–58	37–50, 58–50	50, 50
WR alkali feldspar Or, Ab	–	–	74–72, 23–26	72, 26 at start of stage; 68, 30 at magma T 1022°C8
WR Cpx Mg# range1	–	–	48	48–52	52
WR Cpx En–Fs–Wo	–	–	34–20–46	34–20–46	38–19–43
WR Opx Mg# range1	–	–	34	34–45	45
WR olivine Fo range1	–	–		29–469	46
WR H2O wt % (i.e. fluid phase in WR restite) range1	–	–	1·8–1·3	1·3–1·1	1·1

¹Reported range represents start of stage and end of stage or at equilibration temperature, unless otherwise noted.

²It should be noted that not all phases are present over the same magma temperature ranges (see text for discussion). (For magma temperatures at which phase initially forms, see Table 3.)

³Plagioclase out at magma temperature of ∼1142°C and begins precipitating again at ∼1100°C.

⁴Opx out at magma temperature of ∼1159°C and begins precipitating again at 1013°C.

⁵For clarity, high-Ca cpx En–Fs–Wo provided for start of stage only.

⁶It should be recalled that wallrock always has residual melt to satisfy = 0·04.

⁷Quartz out at magma temperature of ∼1073°C, wallrock temperature of ∼772°C (end of stage 3).

⁸Last occurrence of alkali feldspar at magma temperature of 1022°C, wallrock temperature of ∼852°C.

⁹Olivine joins stable wallrock assemblage at magma temperature of 1022°C, wallrock temperature of ∼852°C.

ol, olivine; pl, plagioclase; opx, orthopyroxene; high-Ca cpx, high-Ca clinopyroxene; low-Ca cpx, low-Ca clinopyroxene; af, alkali feldspar (i.e. sanidine); q, quartz; rhb-ox, rhombohedral oxides; apa, apatite.

Several conventions are applied below. Comparisons between the AFC and FC-only cases are always at the same magma temperature. Masses are reported as per cent relative to the initial mass of the M subsystem. For example, for the case described here, if the initial magma subsystem mass is 1 kg, then the total cumulative mass of anatectic melt added to the M subsystem when thermal equilibrium is achieved is 0·43 kg; therefore, the relative mass of assimilated anatectic melt is reported as ∼43% (Fig. 2d).

Stage 1

Evolution during stage 1 (M temperature = 1279–1187°C) involves cooling and crystallization in M. The enthalpy generated is transported to WR, but no delivery of wallrock partial melt has (yet) occurred because wallrock remains below its solidus temperature. Thus, phase equilibria and geochemical characteristics of the AFC and FC-only cases are identical. Major element trends plotted versus M temperature behave as anticipated for a basaltic system in which olivine is the only fractionating phase (Figs 2c and 3b). The Fo content of olivine for this stage is between 88 and 84 (Fig. 3b). The mass of melt in the magma body at the end of stage 1 is ∼91%, and the mass of the dunitic (actually pure olivine) cumulate in C is 9% (Fig. 2c, Table 4).

Precipitation of olivine leads to the marked depletion of only Ni in M melt during stage 1 (Fig. 4a); Sr, Rb, Ba, La, Nd, and Yb all behave incompatibly in olivine and therefore increase in concentration. Isotope signatures remain constant.

Stage 2

Stage 2 (M temperature = 1187–1145°C) is discussed in two parts. The first (stage 2a), between 1187 and 1169°C, involves a fractionating assemblage of olivine + plagioclase, thus producing a troctolitic cumulate; olivine ceases to be a part of the assemblage at 1171°C. The second (stage 2b), beginning at 1169°C, involves an assemblage of pyroxenes + plagioclase (gabbroic cumulate). The mass of M melt at the end of stage 2 is ∼60% and thus the cumulate reservoir, which constitutes ∼40% of the starting melt mass, is composed of dunite, troctolite and gabbro. Wallrock continues to heat up, but no anatectic melt has been produced. Additional details of stage 2 are provided in Tables 3 and 4.

Stage 2a initiates at 1187°C when plagioclase (An₈₀) joins olivine (Fo₈₄) as a precipitating phase (Fig. 3b and d). Wallrock continues to increase in temperature but no anatectic melt has yet been transferred to M. Plagioclase composes ∼75% of the solid assemblage (by mass), and olivine ceases to crystallize at 1171°C.

Major oxide trends for stage 2a reflect the dominance of plagioclase crystallization (Fig. 3) and the continued crystallization of magnesian olivine. Compared with stage 1, FeO vs M temperature slope changes, becoming positive because the proportion of fractionating plagioclase is relatively high compared with olivine. The slopes of SiO₂, K₂O, Na₂O, P₂O₅, TiO₂ and H₂O vs M temperature are steeper than those in stage 1 because for each 2°C decrement in magma temperature the total mass of solids removed is 5–8 times greater than for each 2°C decrement in stage 1.

Slope changes in plots of concentration vs M temperature are also observed for most of the trace elements (Fig. 4) because the mass of crystals removed per temperature decrement is larger as just noted and/or because the phase assemblage has changed. Sr concentration begins to decrease in response to its compatible behavior in plagioclase. for Ni decrease in response to the changing magma body mineral assemblage; Ni remains compatible and therefore its concentration continues to decrease. The trends for Ba, Rb, La, Nd, and Yb vs magma body temperature become steeper for the same reason as noted above for SiO₂, etc.

Inflections in magma body temperature vs oxide trends are again observed in stage 2b when pyroxenes become fractionating phases starting at ∼1169°C (opx) and 1165°C (high-Ca cpx). Opx is stable for a narrow temperature interval, and ceases to be part of the assemblage when low-Ca cpx begins to crystallize at 1157°C. Mg# values for opx and high-Ca cpx are shown in Fig. 3b, and En–Fs–Wo components are listed in Table 4.

Inflections in concentration vs M temperature slope are seen among the oxides and trace elements at ∼1169°C because of the changing fractionating assemblage. With the exception of Ni and Sr, the trace elements act incompatibly. Inflections in slope are visible when pyroxenes begin to fractionate and/or the total relative mass of fractionated crystals changes appreciably. These inflections illustrate the sensitivity melts have to changes in the identity and relative masses of the fractionating assemblage. As in stage 1, for both stages 2a and 2b, Sr and Nd isotope ratios remain constant.

Stage 3

Stage 3 begins when the magma temperature reaches ∼1145°C and WR is ∼770°C. This is the first step at which the fraction of anatectic melt exceeds the critical melt fraction set a priori at = 0·04, and, consequently, addition of anatectic melt to M begins. During stage 3, WR temperature changes from ∼770 to 772·3°C (Fig. 2b), and ∼18% anatectic melt is added to M (by the start of stage 4) (Table 4). The WR equilibrium phase assemblage (in addition to melt) includes (in order of abundance) plagioclase opx ≈ potassium feldspar > quartz > high-Ca cpx > H₂O (phase) > rhombohedral oxides > apatite (Fig. 2e and f, Table 4). During this stage, partial melting of wallrock is pseudoeutectic (Fowler et al., 2007) as one might anticipate from the ternary minimum model system SiO₂–NaAlSi₃O₈–KAlSi₃O₈ (e.g. Carmichael et al., 1974). As wallrock temperature increases, relative proportions of plagioclase, high-Ca cpx, opx, and rhombohedral oxides that remain in the WR restitic assemblage increase. Of these, high-Ca cpx actually increases in mass (because of the changing phase equilibria as melt is removed from WR) and plagioclase, opx and rhombohedral oxides decrease (Fig. 2e and f). In contrast, relative proportions and masses of alkali feldspar, quartz, H₂O (phase), and apatite decrease, indicating that they are contributing proportionally more to anatectic melt. (Plagioclase Ab, alkali feldspar Or, and Mg# of opx and high-Ca cpx are shown in Fig. 3b and d.) At this stage, the anatectic melt composition is felsic, with ∼72 wt % SiO₂, and its major element composition does not change appreciably until the end of stage 3 (Fig. 3). Quartz is no longer part of the restite at M temperature of ∼1073°C (and WR temperature of 772·3°C), which is very close to the end of stage 3 and cessation of pseudoeutectic melting of wallrock. (The quartz-out compositional evolution of anatectic melt dominates stage 4 and is therefore described in the next section.) By the end of stage 3, WR is composed of 78·4% solid and ∼3·3% (residual) melt (yielding = 0·04).

In distinct contrast to major elements, trace element concentrations in the anatectic melt change, in some cases dramatically (Fig. 4). La, Nd, Yb, and Rb behave incompatibly (⁠ < 1) in WR during melting. Such behavior mimics complete melting of accessory phases by the temperature at which the melt fraction in wallrock exceeds or lack of stability of accessory phases. Although we recognize that important trace element-bearing accessory phases may be present in silicic rocks and magmas, the trace element analysis presented here provides a reference state highlighting the role that major phases can play in trace element mass balance. Incorporation of accessory phases [e.g. zircon using Watson & Harrison (1983)] into MCS is indeed important and will be treated elsewhere. Because of their concentrations in average crust, upon the first step of transfer of anatectic liquid to the magma body, Rb, La, and Nd in anatectic melt have concentrations that are higher than those in M melt at 1145°C, whereas Yb has a concentration very similar to that in M melt at 1145°C. For all of these elements, as WR melting progresses (and M temperature decreases), concentrations decrease (Fig. 4) because these elements are preferentially stripped from the wallrock restitic phases. Ba is interesting because it behaves differently. Because alkali feldspar is a relatively abundant phase in the WR solid assemblage (∼13%) at the beginning of stage 3, is slightly greater than unity (∼1·2). The initial Ba concentration in anatectic melt is lower than that of M melt because of this compatible behavior and the initial concentration in average crust (Rudnick & Gao, 2004; see Table 2). As alkali feldspar melts, its abundance in WR restite decreases to ∼4%, and the bulk solid–melt partition coefficient decreases to ∼0·7. As a consequence, Ba concentration in anatectic melt increases. Because of the dominance of plagioclase in WR restite, for Sr changes little during stage 3, and thus the concentration in anatectic melt increases only very slightly. The Ni concentration of anatectic melt does not change dramatically; stays close to two because the chosen partition coefficients for opx, high-Ca cpx, and/or rhombohedral oxides are greater than unity and the relative abundances of these minerals do not change significantly during this stage.

Figures 2–4 illustrate the impact of assimilation on the temperature, phase equilibria and geochemical evolution of contaminated magma melt, compared with the FC-only case. During stage 3, the magma body temperature drop of 79°C (1145–1066°C) is dominated by addition of a relatively large mass of lower enthalpy anatectic melt; of the 79°C drop, ∼85% of that decrease is due to addition of colder anatectic melt. Differences in the mass of crystals that form and the identity of the crystallizing phases are significant between the AFC and FC-only cases. Crystallization of plagioclase ceases between ∼1142 and 1100°C, whereas it continues in FC-only (Fig. 2c). By the end of stage 3, the mass of low-Ca cpx is lower in the AFC case (∼7 vs 10·5% for AFC vs FC-only), whereas high-Ca cpx abundances are very similar (∼8·5 vs 9% respectively). Thus, in stage 3, the per cent of crystals that precipitate by 1066°C is less than that in the FC-only case (∼49 vs 66%, respectively). This difference leads to differences in the per cent of crystals formed per °C. During stages 1 and 2, M temperature dropped ∼134°C and M crystallized about 40% whereas, during stage 3, M temperature decreases about 79°C but M undergoes only 9% crystallization. These yield per cent crystallization per °C rates of 0·30 for stages 1 and 2 versus 0.11 for stage 3. Thus, although adding lower temperature anatectic melt to a hotter magma body generates significant cooling of the body, enhanced crystallization, which might be expected based on thermal arguments, does not occur. Phase relations respond even more dramatically to the change in magma bulk composition because the magma body subsystem moves away from rather than towards the phase saturation surface. This phenomenon was identified by Reiners et al. (1995) in an analysis that involved use of MELTS to explore isenthalpic AFC. The difference in the mass of fractionated crystals coupled with addition of anatectic melt explains why in the MCS AFC case, the mass of melt in the magma body is much higher compared with FC-only: 69% vs 33%, respectively (Fig. 2a).

Because anatectic melt is being added and the precipitating phase assemblage is different, the magma melt composition is also different between the AFC and FC-only cases (Fig. 3). At the same M temperature, SiO₂ is slightly higher in the AFC magma melt in response to relatively high SiO₂ in anatectic melt. Because during AFC, plagioclase crystallization in M ceased between 1142 and 1100°C, Al₂O₃ concentration in M melt increases slightly during this temperature interval, yielding a change from a negative (stage 2) to positive (stage 3) Al₂O₃ vs M temperature slope. K₂O and P₂O₅ concentrations in M melt increase because they are not abundant oxides in any of the crystallizing phases. In addition, because their concentrations in anatectic melt are high relative to M melt, they are higher in AFC than in FC-only. Concentrations of FeO, TiO₂ and Na₂O in AFC melt are lower than their counterparts in FC-only, and the slopes are opposite: the AFC oxide versus magma temperature slopes are negative because these oxides are diluted by assimilation of anatectic melt, whereas slopes in the FC-only cases are positive. MgO and CaO concentrations, which are virtually identical in the AFC and FC-only magma melts, highlight the combined effects of addition of anatectic melt, which has relatively low concentrations of these oxides (thus diluting them), and removal of different masses of crystallizing phases. Viewed in isolation, one might anticipate that a slightly smaller mass of pyroxene fractionating in the AFC case would yield higher concentrations of CaO and MgO in the magma body melt. However, these ‘enrichments’ are offset by assimilation of anatectic melt that has relatively low concentrations of CaO and MgO. The opposite is the case for H₂O (component) in the magma melt. Although it shows similar behavior for the AFC vs FC-only case, because it acts incompatibly (as there are no H₂O-bearing phases forming in the magma), less crystallization in the AFC case should yield lower concentrations, but the H₂O-‘enriched’ anatectic melt offsets this ‘depletion’. These examples highlight an important conclusion of this work: simple ‘petrological rules of thumb’ may lead to erroneous conclusions about the behavior of major elements in an open system.

Trace element concentrations in contaminated melt vary (Fig. 4) depending on the mass of the element in anatectic melt delivered to the magma body, ⁠, the mass of element in M melt, and the (total) mass of magma melt after contamination. For example, the concentration of Sr in contaminated melt is typically slightly higher than it is for the same M temperature in the FC-only case (Fig. 4b). This observation suggests that the dilution afforded by assimilation of relatively low Sr concentration anatectic melt is compensated by the reduction in the mass of fractionating plagioclase; for Sr in the AFC case is significantly lower than that in the FC-only case (∼0·3 vs 1·5, respectively). The proportion of Sr from anatectic melt added to M melt, calculated as a percentage of the initial mass of Sr in the magma body, is relatively small (up to 6% by the end of stage 3), but still imparts a distinctly higher Sr isotope signature of ∼0·7053 to magma melt by the end of stage 3 for the equilibrium case presented.

Ba in stage 3 is less abundant in magma melt of AFC than FC-only. Although for the two cases are similar, the concentration of Ba in anatectic melt is lower than that of M melt throughout the stage. Thus, addition of anatectic melt dilutes Ba in magma melt. Yb is also diluted in the AFC magma melt compared with the FC-only melt.

Concentrations of La and Nd in contaminated magma compared with the FC-only case illustrate what appears initially as a paradoxical result. These elements behave incompatibly in the magma and wallrock subsystems, and one might anticipate that such elements should be highly enriched in contaminated melt compared with the FC-only case. Indeed, this might be the expectation particularly given that the concentrations of both in anatectic melt through most of stage 3 are higher than that of magma melt. Instead, Fig. 4 shows that these elements have concentrations that are lower than those in the FC-only case. The explanation for this lies in the mass balance relationships between mass of melt and mass of trace element in the contaminated magma compared with the same masses in the FC-only case. It is straightforward to recognize that the FC-only case has a monotonically decreasing mass of melt, whereas the mass of melt in the contaminated case increases during stage 3. Figure 2a shows that by the end of stage 3, the AFC melt mass per cent is ∼69%, compared with ∼33% for the FC-only case. Because concentration is defined as the mass of element/mass of melt, the very different melt masses have a profound impact on trace element concentrations. For La and Nd in stage 3, the ratio of mass of M melt in AFC to mass of M melt in FC-only (herein called the melt mass ratio) is greater than the ratio of mass of trace element in M melt for AFC to mass of trace element in M melt in FC-only (herein called the trace element mass ratio). These differences yield concentrations of La and Nd in magma melt for the AFC case that are lower than for the FC-only case. The per cent Nd contributed by anatectic melt is 44% (by the end of stage 3) of the initial mass of Nd in the magma body, and the ¹⁴³Nd/¹⁴⁴Nd attains a value of ∼0·5127 for the equilibrium case. It should be noted that the per cent of ‘crustal’ Nd is very different from that of Sr (up to 6%), and the total mass transfer by the end of stage 3 is ∼18% (using the total initial mass of melt in the magma body as denominator). These value differences highlight the ambiguity inherent in defining crustal ‘indices’ that are based on assessment of single elemental contributions and underscore another significant conclusion: crustal and mantle contributions (i.e. indices) should be cast as functions of total mass rather than masses of single elements.

Perhaps equally non-intuitive is the concentration of Ni in M melt of the AFC case, which is higher than that of the FC-only case (Fig. 4a). is >1, in part because Ni is compatible in pyroxenes and spinel. Its concentration in the AFC case is higher than that in the FC-only case because the melt mass ratio is less than the trace element mass ratio. A primary influence on the mass relationships is the initial concentration for average crust (59 ppm, Table 2); at the magma temperature at which assimilation begins, anatectic melt has similar or higher concentrations of Ni compared with the magma, and thus it is ‘enriched’ by contamination. In summary for the case examined here, which involves modeling of typical high-alumina basalt intruding average continental crust, the classification of some elements as ‘enriched’ or ‘depleted’ is not a simple function of the element behavior (incompatible vs compatible) and concentration in the contaminant. Element concentrations are inextricably tied to the total mass of the magma subsystem via the energetics and self-consistent phase equilibria, and thus changes in magma melt mass that occur because of anatectic melt addition must be strictly tracked to correctly model open-system trace element and isotope behavior. Simpler models (e.g. DePaolo, 1981) that do not incorporate major element mass conservation and consistent phase equilibria may yield incorrect predictions or at best are incomplete and should be used with caution.

One of the fundamental conclusions of this work is that the trace element mass balance impacts of assimilation on the magma subsystem are non-linear functions of the trace element mass added by assimilation, the trace element mass in the magma body, the mass removed by a potentially distinct fractionating assemblage (compared with FC-only) and the mass of magma melt. Magma melt mass is influenced by the mass of anatectic melt delivered to M for each step of the assimilation process, which is, in turn, a function of ⁠. Such complex AFC mass balance relationships require rigorous phase equilibria and associated self-consistent mass and enthalpy (sensible and latent heat) balance calculations.

Stage 4

Stage 4 begins when M is ∼1066°C and WR is ∼772·3°C; both subsystems evolve to ∼945°C, the equilibration temperature. Enthalpy delivered to wallrock between magma temperatures 1066 and 1064°C (first step of stage 4) raises the WR temperature from 772·3 to ∼777°C, a larger temperature increase than all stage 3 steps combined. This highlights the distinction between a system’s enthalpy and its temperature. In stage 3, enthalpy delivered to WR affected a relatively small change in temperature but a relatively large mass was melted; thus, a large proportion of the enthalpy (latent heat) delivered to WR was applied to cause melting at essentially isothermal conditions. Stage 4 is different in that proportionally more enthalpy is used to heat wallrock up, and for each step, proportionally less anatectic melt is generated. Also in contrast to stage 3, during stage 4, anatectic melt changes composition during progressive partial melting.

At the beginning of this stage, the WR solid assemblage is plagioclase opx > alkali feldspar > high-Ca cpx > rhombohedral oxides ≈ H₂O (phase) > apatite (Table 4). Quartz disappeared from the assemblage at the end of stage 3. Plagioclase (Ab₅₈ at start of stage; see also Fig. 3d) and rhombohedral oxides continue to increase proportionally in wallrock restite throughout stage 4. Opx (Mg# 34 at start of stage) continues to increase until ∼1024°C, when its relative proportion begins to decrease. The relative abundance of high-Ca cpx (Mg# 48 at start of the stage) is approximately constant from 1066 to 1024°C, but below 1024°C, its relative abundance decreases. Alkali feldspar (Or₇₂ at the start of the stage) continues to decrease in proportional abundance; magma temperature of ∼1022°C and wallrock temperature of ∼852°C mark the last appearance of alkali feldspar in WR restite. At approximately this same M temperature, olivine (Fo₂₉) joins the stable WR assemblage and increases in proportional abundance until the end of the simulation, at which point it has become slightly more magnesian (Fo₄₆); it should be noted that this compositional change reflects the increasingly mafic nature of the wallrock restite, consistent with removal of anatectic melt that is more silicic than the initial bulk composition of the wallrock. H₂O (phase) and apatite continue to decrease in proportional abundance, but remain as restitic phases at the equilibration temperature.

The change in anatectic melt composition observed at the end of stage 3 characterizes much of stage 4 and is therefore discussed in this section (Figs 2–4). Because quartz is no longer a residual solid, the SiO₂ concentration of anatectic melt begins to decrease. MgO, Al₂O₃, CaO, FeO, TiO₂, K₂O, Na₂O and H₂O (component) begin to increase in concentration in the anatectic melt in part because SiO₂ is proportionately less concentrated. P₂O₅ concentration decreases because the proportional role played by apatite in the melting process is decreasing, and, as with SiO₂, P₂O₅ is being diluted in the anatectic melt. At M temperature of ∼1022°C and wallrock temperature of ∼852°C, the last step at which alkali feldspar is a part of the WR restite, K₂O and Al₂O₃ are at their peak concentrations. Above this WR temperature, K₂O concentration begins to decrease by dilution. The contribution of Al₂O₃ from melting of both plagioclase and alkali feldspar is greater than plagioclase alone, and therefore, its concentration begins to decrease as well. The flat to positive slopes of MgO, CaO, FeO, TiO₂, Na₂O and H₂O (component) vs magma temperature reflect the variable contributions that high-Ca cpx, opx, plagioclase, H₂O (phase), and rhombohedral oxides make to anatectic melt.

Although Rb, La, Nd, and Yb all have < 1 (i.e. incompatible), concentrations in anatectic melt consistently decrease through stage 4 (Fig. 4) because the mass of trace element in WR restite is decreasing more quickly than the total mass of restite. Between the start of stage 4 and ∼1022°C, Ba concentration in anatectic melt increases because alkali feldspar is melting. The proportion of this mineral is decreasing in the restite, so the bulk solid–melt partition coefficient gradually decreases. The maximum Ba concentration in anatectic melt is approximately coincident with the disappearance of alkali feldspar, at which point the concentration begins to decrease. In contrast to these elements, Ni initially shows rather flat concentration vs M temperature trajectories. At ∼1022°C, its concentration in anatectic melt increases slightly to moderately because the relative proportion of high-Ca cpx in restite begins to decrease. As a consequence, for Ni gradually becomes slightly less compatible, imparting increased concentration of trace element to anatectic melt. Sr concentration in anatectic melt increases slightly throughout stage 4 in response to the slight decrease in as the relative proportions of phases change in WR restite.

At the start of stage 4, the magma crystallizing assemblage includes plagioclase, high-Ca cpx, and low-Ca cpx. Although plagioclase returns to the stable assemblage at 1098°C, close to the end of stage 3, the total per cent crystallized was only 0·4 between that temperature and 1066°C. Thus, the impact its fractionation has on the magma in stage 3 is minor and discussion was delayed to this section. The cumulative mass per cent of plagioclase that fractionates is lower in the AFC than in the FC-only case at the same magma temperature; thus, plagioclase crystallization is clearly suppressed by assimilation. Low- and high-Ca cpx continue to crystallize during stage 4; high-Ca cpx crystallizes to within a few degrees of the equilibration temperature, and low-Ca cpx ceases to crystallize at 1022°C. In the FC-only case, both minerals crystallize to 945°C. Assimilation reduces the temperature range over which low-Ca cpx precipitates, and the total mass formed is smaller. By the end of the stage, the compositions are also slightly different: Mg# 47 and En–Fs–Wo 37–20–43 for AFC vs Mg# 35 and En–Fs–Wo 29–29–42 for FC-only. Orthopyroxene begins to crystallize again at ∼1013°C and precipitates until the equilibration temperature. In contrast, it does not crystallize again in the FC-only case, indicating that addition of anatectic melt enhances its stability. Spinel joins the stable assemblage at 1048°C and continues to crystallize until 945°C; compared with the FC-only case, spinel begins crystallizing at a slightly lower magma temperature (1048 in AFC vs 1081°C in FC-only). In summary, at the equilibration temperature, it is apparent that assimilation has suppressed the total mass of plagioclase, low-Ca cpx, and spinel crystallization and enhanced slightly crystallization of high-Ca cpx and opx. Olivine is not affected by AFC.

Magma major oxide compositions obviously respond to the complex interplay of changing anatectic melt composition and mass, and mass and identity of fractionating phases (Figs 2 and 3). The impact of AFC on SiO₂ is subtle. For the part of the trajectory in which anatectic melt has higher SiO₂ than the magma melt, SiO₂ in the AFC case is slightly higher than in the FC case. Eventually, the SiO₂ concentration in anatectic melt is equal to or lower than that in M melt, and SiO₂ in the AFC case is slightly lower than that in the FC-only case. Al₂O₃ decreases throughout the stage, but it is always higher than in the FC-only case at the same magma temperature. The slope of Al₂O₃ vs M temperature changes from slightly positive to negative at ∼1066°C, when the mass of plagioclase crystallizing per M temperature decrement approximately doubles. At this same temperature, Al₂O₃ in anatectic melt starts to increase in concentration, but despite this, its concentration is still lower than that of M melt, so it reinforces the declining concentration in magma melt. At the same magma temperature, the mass of plagioclase fractionating in the FC-only case is always greater, and this contributes to lower Al₂O₃ in the FC-only case. The concentrations of K₂O and P₂O₅ in M melt increase throughout this stage because there is no significant crystallization of a K- or P-bearing phase. Both have M melt concentrations that are higher than their FC-only equivalents. Although there is less cumulative crystallization, suggesting that K₂O and P₂O₅ should be lower in the AFC case, the significantly higher concentrations in anatectic melt for most of the stage render the contaminated melt concentrations greater than in the FC-only case. Na₂O changes very little through most of stage 4, rising slightly at the stage’s end when the concentration in anatectic melt is much higher than that of magma melt. FeO concentration decreases throughout the entire stage, in response to crystallization of spinel and pyroxenes. At the beginning of stage 4, TiO₂ increases slightly compared with its concentration at the end of stage 3. The most likely explanation for this is the reduction in the proportional mass of high-Ca cpx that fractionates per M temperature decrement. The decrease in TiO₂ at ∼1048°C is in response to spinel as a (new) fractionating phase. AFC TiO₂ is lower than that in FC-only because anatectic melt has a significantly lower concentration of the oxide, and thus dilutes the magma melt. (The major change in trajectory for the FC-only case occurs at the temperature at which spinel stabilizes.) MgO, CaO, and H₂O (component) in the AFC versus FC-only magma melts are similar for reasons explained in stage 3.

Like the major elements, trace elements also respond to the changing magma fractionating assemblage, and composition and mass of anatectic melt. Ba, La, and Nd generally have positive concentration vs M temperature slopes and are less abundant than their FC-only counterparts. During stage 4, the slopes of each vary slightly in response to changing bulk solid–melt partition coefficients and the mass and composition of anatectic melt. Nd isotopes smoothly decrease to ∼0·5125 by the equilibration temperature. The contribution of crustal Nd is significant by the end of stage 4; it is ∼90% (relative to initial magma body Nd). Toward the end of stage 4, La evolves to a negative concentration vs M temperature slope, probably through dilution by anatectic melt. As in stage 3, Rb mostly increases in concentration with decreasing magma temperature and is always greater than that in the FC-only case. Sr concentration decreases for all of stage 4 because plagioclase is the dominant magma crystallizing phase, and thus >1. The final ⁸⁷Sr/⁸⁶Sr is ∼0·7084, and the total per cent Sr from anatectic melt relative to initial magma body Sr is ∼16%. The slope of Ni concentration vs magma temperature is initially (slightly) negative because is well above unity as pyroxenes and spinel are crystallizing. At the same approximate M temperature at which anatectic melt begins to increase its Ni concentration, the slope of the Ni–M temperature curve becomes positive, and thus AFC Ni in magma melt is greater than that of the FC-only case.

The above analysis highlights the interplay evident among magma element or oxide concentration, magma fractionating assemblage and composition and mass of anatectic melt. For example, Al₂O₃ might be expected to decrease in concentration, particularly at the onset of assimilation because anatectic melt has a concentration that is significantly lower than magma melt, but cessation of plagioclase crystallization for ∼40°C leads to higher concentration than the FC-only case. The concentration of H₂O (component) might be expected to be lower in the AFC case because at all stages after assimilation begins, less total crystallization (of an anhydrous assemblage) occurs. Instead, addition of relatively H₂O-rich anatectic melt yields AFC and FC-only cases that have similar H₂O concentrations.

Recap: the utility of the Magma Chamber Simulator

The example described above illustrates the rich array of predicted quantities generated by the MCS, and underscores the critical importance of correctly evaluating energetics, multicomponent–multiphase phase equilibria and mass exchange. Insights derived from this single analysis include phenomena such as suppression of crystallization despite addition of cooler anatectic melt, concentrations of incompatible trace elements in magma contaminated by continental crust that are lower than their FC-only counterparts, and different crustal contamination ‘indices’ for each element evaluated. These insights and others require a rethinking of the geochemical and petrological impacts of assimilation on a magma body and the conventions by which assimilation is identified in rock suites. Although changing isotopic characteristics remain an indisputable clue for open-system evolution, trace element behavior will vary because of not only the addition of specific trace elements but also the addition of total mass to the contaminated body. Hence, the concentration of a given trace element in contaminated magma melt can be lower than in the FC-only case despite incompatible behavior in both magma and wallrock. Anatectic melt compositions that change as assimilation progresses illustrate that formulations using constant composition assimilant are flawed; although pseudoeutectic melting produces major oxide concentrations that are generally constant, trace element variations can in fact be significant. Some elements that act incompatibly in WR decrease in concentration in anatectic melt as wallrock partial melting progresses, whereas some compatible elements actually increase in concentration. Thus, the impact of addition of such melts to the magma body can be the opposite of that assumed for addition of continental crust. The conclusions drawn from this single AFC simulation illustrate the illuminating insights that can be derived by systematic confrontation of MCS predictions and underscore the potential MCS has in modeling natural systems.

THERMOMECHANICAL CONSTRAINTS APPLIED TO THE MCS

Overview

Because the MCS is a thermodynamic, not transport, model, time and spatial scales relevant to mass, momentum and heat transfer are not explicitly considered, nor are magma transport properties part of the input specifications. Consequently, the applicability of the MCS to natural systems hinges critically on the validity of implicit MCS assumptions regarding heat and mass transfer. For example, in MCS v1.0, once a threshold melt fraction is present in wallrock, partial melt above this limit is instantaneously transferred to the resident magma M. Is there a rheological basis for this assumption? Similarly, what is the timescale for percolative flow of partial melt from WR to M relative to the solidification timescale of M? If percolative flow is slow relative to M solidification, contamination cannot take place on the appropriate magma solidification timescale, and MCS results become moot by a metastability argument. Similarly, the extent of partial melting in WR is ultimately controlled by the temperature field along and near the wallrock–magma body boundary. How thick is this region, and how does it evolve in time? How does that timescale compare with timescales associated with magma recharge? For open systems with significant recharge, what limits can be placed on rates of recharge so that the system can be considered dominated by recharge relative to cooling and solidification? Based on these considerations, a fundamental question can be posed: Are the implicit assumptions intrinsic to the MCS consistent with the characteristics of natural magmatic systems? The detailed answer to this question provides explicit information on the limitations of the MCS in its present form (v1.0) and serves to illustrate how future modifications can enhance its applicability.

Because rates of heat and mass transfer between M, WR and R subsystems in natural settings exhibit considerable variation depending on magma material properties, boundary and initial conditions, wallrock thermophysical properties and the replenishment rate and properties of recharge magma, the goal of this section is to investigate by scale analysis the timescales for heat and mass exchange implicitly intrinsic to the MCS model. In particular, heat transport between resident magma and wallrock is analyzed for both closed and open systems including calculation of thermal timescales for systems dominated by conduction, convection and conjugate heat transfer. The timescale of Darcy percolative flow is also addressed and compared with magma solidification times to discover if the assumption of rapid percolative flow is generally valid in typical natural systems. The rheological basis for selection of the threshold melt fraction is also presented. The overall conclusion is that implicit timescales for both heat and matter transport in the MCS are consistent with those in natural systems in general. In addition, a guide to how the MCS parameters ⁠, Λ and Φ relate to system dynamics and may be chosen to model particular petrogenetic scenarios is provided. Finally, it should be noted that the utility of the MCS does not lie solely in the generation of ‘reference’ RAFC models applied to better understand particular magmatic systems. The MCS can also be applied to more general problems by sequential application to emulate the long-term (10⁴–10⁷ years) consequences of protracted episodes of crustal heating, partial melting, injection of mantle melts and consequent crustal growth and by application to shallow volatile-saturated hydrothermal–magmatic systems in island arc and continental settings with implications for eruptive mechanisms, volcanic hazards and economic mineralization. These problems will be presented in detail elsewhere.

Heat transfer

For a magma body to crystallize, heat must be dissipated. It is important to distinguish for timescale analysis between open and closed systems and for systems governed strictly by conductive cooling versus mixed conductive–convective cooling. Convection may be important both within the magma body and within wallrock where hydrothermal supercritical and subcritical advection of heat may dominate.

For closed systems dominated by conduction, the relevant timescale is the time required for sufficient heat to be extracted so that a fixed mass of impulsively or quasi-impulsively emplaced melt, initially at its liquidus, undergoes phase transformation to a crystalline solid or rock at its solidus. Allowing for both sensible and latent effects, ∼1 MJ kg⁻¹ must be dissipated to accomplish the required cooling and phase transition for a mafic bulk composition and ∼0·7 MJ kg⁻¹ for a silicic bulk composition. There is a long history of thermal calculations that model the cooling and crystallization of magma bodies of various geometries including latent heat effects and variable transport properties. Classic calculations of the type pioneered by Jaeger (1957, 1968; see also Carslaw & Jaeger, 1959) mainly handle heat conduction and apply various boundary and initial conditions (fixed temperature or fixed thermal gradient). A characteristic of most conduction models is that the wallrock–magma body temperature quickly attains a value equal to the average and slowly decays towards the far-field wallrock temperature as time progresses. In more sophisticated but still conduction-dominated models, conjugate solutions have been obtained in which conditions within the magma and wallrock far afield are specified as boundary conditions and the solution provides an estimate of the boundary temperature as a function of time. In virtually all these types of solutions, the solidification cooling timescale (τ_S) is related to a characteristic length of the magma body (⁠ where V is the magma volume) and the thermal diffusivity of the wallrock, κ = k/ρ according to τ_S ≈ /κ, where the characteristic length scale of the problem is based on the resident magma volume. Characteristic solidification times for magma volumes of 0·1, 1, 10 and 100 km³ are 67 kyr, 310 kyr, 1·4 Myr and 6·7 Myr, respectively, with nominal values of k = 1 W m⁻¹ K⁻¹, ρ = 2800 kg m⁻³ and an effective isobaric specific heat capacity that accounts for latent and sensible heat of = 3500 J kg⁻¹ K⁻¹. Accounting for variable properties and different magma body geometries, timescales will vary. Under most circumstances, however, τ_S remains within the range 10⁴–10⁷ years for closed conductive systems. Of the important factors, geometric effects (size and shape) are the most critical. For example, in conductive models, solidification times decrease as the wallrock–magma body contact area increases. For example, the surface area ratio for a cubical body to that a sheet-like body of volume ab² and aspect ratio A = b/a goes as A^2/3. Because the solidification time is proportional to wallrock/magma body surface area, for the same volume and wallrock–magma properties, cooling times for a sheet of aspect ratio A = 10 are ∼5 times shorter. For extremely thin sheets (say A ∼ 100), cooling times are faster by a factor of ∼20 compared with an equal volume of magma in an equant body: a 100 km³A = 100 sheet-like magma body solidifies in 135 kyr.

Remaining within the framework of closed systems, models that account for magma convection and/or hydrothermal convection in surrounding wallrock have also been extensively developed in the past 30 years (Cathles, 1977; Norton & Taylor, 1979; Spera, 1980; Huppert & Sparks, 1988; Marsh, 1989). Here the essential issue is the ratio of thermal resistances between wallrock and magma body. In purely conductive models, these are essentially equal except for small differences in the thermal conductivity of wallrock and magma. In models allowing for magma convection and wallrock porous medium convection, thermal resistances can be substantially different. The prototypical magma convection–wallrock hydrothermal advection model of Bejan & Anderson (1983) serves as an illustration of a conjugate country rock–magma body heat transfer system. Similar models have been described elsewhere (Spera et al., 1982; Carrigan, 1988; Bergantz, 1989, 1992; Bowers et al., 1990). In this conjugate heat transfer problem, an impermeable vertical surface separates a porous wallrock of temperature at great distance from the contact from a magma reservoir of interior temperature ⁠. The boundary temperature, the one relevant to the MCS model, is the result of the interaction of two buoyancy-driven flows, one in the wallrock treated as a porous medium and the other in the magma body. The wallrock–magma boundary temperature (T_b) is determined by the relative vigor of hydrothermal convection in the country rock compared with magma convection within the magma. The thickness scale of the porous medium boundary layer is δ^WR = Rd^–1/2 whereas the magma body thermal boundary layer scale is δ^M = Ra^–1/4. Here Rd is the Darcy Rayleigh number,

for the wallrock treated as a porous medium with permeability K and

is the Rayleigh number for the magma body treated as a viscous fluid. Other parameters are defined in table 1 of Spera & Bohrson (2001) and are not repeated here. The continuity of the heat flux across the boundary gives rise to a dimensionless parameter B = (k^WRRd^1/2)/(k^MRa^1/4) that controls the average wallrock–magma body boundary temperature as well as the magnitude of the heat flux across the boundary. Physically, B is the effective conductivity ratio of wallrock to magma when conjugated boundary layers, an upwelling one in wallrock and a downwelling one in the magma, are present. The value of B increases with increasing temperature difference between wallrock and magma as B ∝ ΔT^1/4 but decreases with the magma body characteristic length according to B ∝ ⁠. When B is large, hydrothermal convection efficiently carries away heat at the boundary, and the boundary temperature assumes the temperature of the porous wallrock (T_b → ⁠). The solidification time of the magma is significantly shorter (by a factor of 10–100) than the standard conduction solidification time owing to rapid advection of heat away from the contact. The heat flux between magma and wallrock scales as q ≈ k^MΔTRa^1/4/⁠, which leads to vigorous hydrothermal convection in wallrock but relatively little wallrock partial melting owing to the efficient advection of heat away from the wallrock–magma body boundary. In the MCS formulation, this limiting case is accomplished by setting the wallrock/magma mass ratio Λ 1, which ensures low wallrock temperatures. This condition might prevail in the upper crust in regions undergoing extensional shear failure as a result of tectonic forces where a vigorous and deep hydrothermal system is established perhaps transiently or episodically (e.g. see Ingebritsen & Manning, 2010). In the opposite limit where B → 0, magmatic heat cannot be transported rapidly away from the country rock–magma contact, and the boundary temperature approaches the magma temperature (T_b → ⁠). Significant partial melting of wallrock occurs in this case, and the heat flux from magma to wallrock scales as q ≈ k^WRΔTBRa^1/4/⁠. This condition might be met in low-permeability catazonal rocks (20–40 km depth) where hydrothermal flow is negligible. In this case, very steep thermal gradients exist in the wallrock and significant partial melting occurs in a zone that laterally propagates at velocity v ≈ a_c(κ/t)^1/2, where a_c is a constant of order unity (e.g. equal to 1·6 for a planar vertical wall). In the MCS formulation, employing Λ in the range 0·5–2 can attain this limit. An extreme end-member that can be modeled in MCS is the assimilation of stoped blocks. In this case Λ should be set to a small value.

In the general case, the boundary temperature is where the function f(B) approaches ½ as B → 0 and −1/2 for B → ∞. Values of f(B) for intermediate values of B have been given by Bejan & Anderson (1983). When both sides of the boundary are lined by boundary layers, the heat flux from magma to wallrock measured perpendicular to the wallrock–magma boundary is

It is important to note that the extent of contamination by country rock is related not only to the contact temperature but also to wallrock bulk composition, pressure and thermodynamic factors such as heat capacity, and fusion enthalpies and phase equilibria, all factors that are taken into account in the MCS.

For open systems in which magma recharge is important, the relevant timescale is determined from the balance between heat flow to wallrock and heat delivered to resident magma via recharge of magma of higher specific enthalpy. Where a system resides along the open-system spectrum can be characterized by defining a dimensionless ratio Π of the solidification time of resident magma of volume V_M to the time required to replenish melt lost by crystallization given a fixed rate of recharge ⁠; that is,

For Π > 1, the system is recharge dominated; solidification cannot keep pace with recharge, and the magma body mass grows. In contrast, if Π < 1, recharge cannot compete with solidification and, in a strict RFC process, the magma body shrinks. As an illustration, let us consider a magma body of 1000 km³ to which recharge is being added at rate of 3 × 10⁻⁵ km³ a⁻¹. With a nominal thermal diffusivity of 10⁻⁷ m² s⁻¹, Π = 1, which indicates approximately a volumetrically steady-state magma body. The Π number scales with Φ, an MCS input parameter defined as the ratio of recharge increment mass to initial resident magma mass according to

where V_R is the increment of recharge magma added to resident magma during an episode of recharge. Because the MCS is a thermodynamic and not transport model per se, to convert from Π to Φ, one must define the total recharge magma added during a given replenishment event.

The conditions for magma body growth can be addressed in more detail by noting that growth occurs when the magma body–wallrock boundary temperature remains at or above the wallrock solidus during the time interval between discrete replenishment events. The simple Laplace-type conduction model of Hardee (1982) gives the wallrock–magma body boundary temperature at the onset of each successive recharge event for different values of the volumetric (time-averaged) recharge flux to the magma body. When the wallrock boundary temperature reaches its solidus, the condition for magma body formation, growth or steady behavior is roughly satisfied. In the MCS, the relative importance of recharge is measured by the ratio of recharge mass to initial resident magma mass, Φ. An illustration of this model that also demonstrates the role of magma body shape can be addressed. For a diapiric pipe intrusion [V = (π/4)ba²] or rectangular parallelepiped (V = a²b), a detailed calculation shows that if the rate of basaltic recharge is >2 × 10⁻³ km³ a⁻¹ (∼6 Tg a⁻¹), typical wallrock solidi temperatures (a function of bulk composition and pressure) of ∼700–900°C are reached after about a dozen intrusive pulses. The calculated volumetric rate is approximately the same as observed eruption rates integrated over decades at typical active volcanic centers. For example, in the period 1961–2001, the average volumetric flux of magma erupted at the active volcano Pacaya, Guatemala, where geochemical evidence indicates that RAFC is important, was ∼5·3 × 10⁻³ km³ a⁻¹ (Durst, 2008). For volume rates <5 × 10⁻⁴ km³ a⁻¹ (∼1·5 Tg a⁻¹), the wallrock solidus is not reached even after several hundred intrusive pulses because conductive heat loss rates exceed advective, latent and sensible heat gains, given typical magma and wallrock properties. In contrast, for a dike or sheet-like intrusion of volume ab² and aspect ratio b/a = 10, higher intrusion rates of order ∼0·1 km³ a⁻¹ are required for wallrock partial melting. The transition from an early, more brittle style of magma emplacement (e.g. swarms of sheet-like propagating dikes) to a later more diapiric (cylindrical) structure might be anticipated in response to heat input.

In summary, the most critical factor affecting the solidification time of a magma body is whether the system is closed or open. Closed-system timescales vary as with a proportionality coefficient that varies with the mode of heat transport (e.g. conductive, convective, conjugate, etc.) in magma and in wallrock. A continuum of behaviors exists with a strong dependence on depth of magma–wallrock interaction, as depth controls the ambient crustal temperature and the thermophysical properties of wallrock such as permeability and the effective thermal conductivity. The MCS can accommodate a variety of system characteristics by appropriate choice of Λ. In closed systems or systems with limited recharge emplaced in the middle or deep crust, boundary temperatures are high and the effective Λ is of order 0·2–0·6 in strict AFC scenarios. The extreme limit of stoped blocks corresponds to small Λ of order 0·1. In contrast, at shallow levels where magmatic heat can be rapidly advected away from the magma body rock, boundary temperatures remain relatively low and assimilation of wallrock partial melts is limited. In this case, resident magma solidification is relatively short and assimilation of wallrock partial melts is limited. In the MCS formalism, Λ lies in the range 0·6–3 in this case. In RAFC scenarios, where a particular system lies on the open–closed continuum depends most critically on the rate of recharge supply and the thermal conductance of the wallrock surrounding the resident magma body. In general, Π scales as Φ^1/3. Typical values of Φ lie in the range 0·2–1 in the MCS model. Finally, it should be noted that there is an implicit coupling between recharge and assimilation. That is, in scenarios when the recharge magma has a higher specific enthalpy than the resident magma, additional enthalpy is available for heating and partial fusion of wallrock. Because phase relations change as the bulk composition of the resident M changes, the coupling between phase equilibria, assimilation and heat transfer is nonlinear.

Stability field model for magma chamber growth

Within the context of the MCS, it is useful to consider the stability and longevity of magma bodies. The main variables governing the growth and geochemical evolution of a magma body include its initial composition, volume (mass), surface area (A) to volume (V) ratio (geometric form), composition, enthalpy and flux of magma recharge delivered to the body, the ambient state of wallrock (e.g. its composition, temperature, state of stress, thermophysical properties) and the mean depth (pressure) of wallrock–magma interaction. Two regional-scale parameters, with connections to wallrock rheology and net magma buoyancy, include the initial (undisturbed) geotherm and the crustal density structure prior to the onset of magmatism. To obtain a bird’s-eye view of how these factors influence magma body development, Karlstrom et al. (2009, 2010) have mapped out a three-parameter magma chamber ‘stability field’ model pertinent to emplacement of both ‘wet’ and ‘dry’ magmas in arc and continental settings. The regime diagram was constructed by locating the boundaries between four distinct temporal magma system behaviors including (1) in situ solidification, (2) runaway magma body growth, (3) magma body growth leading to eruption, or (4) stable magma body with approximately constant melt volume. This stability field approach can be elucidated and quantified by application of the MCS using the Karlstrom and coworkers parameters of (1) initial magma volume, (2) recharge magma flux, and (3) mean pressure (or depth) of magma storage and evolution. The factors upon which the regime diagram is built are linked by dynamic feedback loops. For example, the orientation, magnitude and spatial variation of the principal stresses exert a first-order constraint on the geometric form magma bodies assume (e.g. pipes, plutons, dikes, sills, sheets) and the depth to which magma may intrude. In turn, the magma body A/V ratio (or average A/V for a collection of bodies such as a swarm of magma-filled cracks) exerts a strong influence on the rate of heat transfer to wallrock and hence its thermo-rheological properties, the timescale for magma crystallization, and the relative role of assimilation of wallrock partial melts into magma. The recharge flux is itself related to the orientation of the principal stresses, the local geotherm, the crustal density structure and possible subsolidus convective upwelling in the upper mantle. Simply put, the act of intrusion affects the thermo-rheological properties of the wallrock, which, in turn, determine the rate of magma cooling (or heating) and crystallization or growth of the magma body. The ‘stability field’ approach developed by Karlstrom and coworkers rests upon the thermomechanical coupling between the local stress field and magma pressure, the thermal balance between recharge heat advection, magma cooling and crystallization versus heat loss to wallrock via conduction or hydrothermal fluid advection, and the viscoelastic relaxations that occur in wallrock as it heats up and possibly undergoes partial melting. MCS simulations use these parameters either implicitly or explicitly, and thus results can be related to the ‘stability field’ representation of Karlstrom et al. (2010). In this way, multiple sequential MCS solutions can be utilized to focus upon long-term crustal evolution issues. An example application is modeling the incipient stage of island arc volcanism in regions where subduction has begun geologically recently (e.g. Leat et al., 2003). For MCS calculations relevant to systems with disturbed geotherms a sequence of MCS calculations of increasing ⁠, Λ and Φ and decreasing pressure can be conducted to approximate the long-term effects of magma underplating into the crust. Because the extent of thermal interaction between wallrock and magma can be highly variable, a series of such MCS sequential calculations can be used to emulate many different petrological scenarios.

Wallrock critical melt fraction and application to MCS

In the MCS formulation, it is assumed that when the fraction of melt in wallrock exceeds a critical fraction (⁠⁠), partial melt is ‘instantaneously’ removed and added to the resident magma M such that the melt fraction remaining in wallrock remains precisely at the critical value. This procedure is based on the concept of a percolation threshold; that is, partial melt is immobile until a threshold melt fraction is reached, after which the melt is mobile relative to the solidification timescale. In this section, we examine the validity of this assumption. In particular, we compare the rate at which melt can percolate from wallrock to resident magma with the solidification timescale. Although wallrock may contain partial melt, resident magma contamination requires material transport across the wallrock–magma boundary and therefore the percolation timescale must be short relative to the solidification timescale.

Many factors that govern the properties of multiphase, partially molten wallrock mediate the transport of wallrock partial melt into an adjacent magma body (e.g. Tait & Jaupart, 1992; Wark & Watson, 1998; Cheadle et al., 2004). For a fixed bulk composition, the wallrock melt fraction is a function of temperature and depth (pressure). Mineralogical heterogeneity in the protolith wallrock leads to a patchy distribution of melt on scales defined by the mineral abundance heterogeneity length scale, typically 0·01–1 m. At temperatures below the brittle–ductile transition, melt is collected in swarms of cracks or fractures approximately orthogonal to the minimum principal stress. Fracture formation is triggered by anisotropic regional tectonic stress and the release of thermal stress near wallrock–magma contacts where temperature gradients are steep. The magnitude of lateral pressure gradients in brittle wallrock can be estimated assuming contact heating is approximately isochoric in the brittle regime. In this case, the lateral pressure gradient is

where α and β are the isobaric expansivity and isothermal compressibility respectively of wallrock and dT/dy is the lateral temperature gradient along the wallrock–magma boundary. The y-component velocity of wallrock partial melt into resident magma owing to Darcy flow is given by

where η is the dynamic viscosity of partial melt in wallrock and the y-coordinate is orthogonal to the wallrock–magma body contact. When thermal stresses are the primary cause of fracture, the migration of anatectic melt can be related to the wallrock–magma thermal gradient and the timescale for percolative (Darcy) flow can be estimated according to

where δ is the thickness of the wallrock–magma boundary thermal layer and ΔT is the temperature drop. The condition that percolative flow is rapid on the timescale τ_solid is that

That is, the percolation timescale is short in comparison with the heat loss or solidification timescale. Adopting the typical parameters η_m = 10⁴ Pa s, β = 10⁻¹¹ Pa⁻¹, δ = 10 m, κ = 10⁻⁷ m² s⁻¹, α = 10⁻⁵ K⁻¹, ΔT = 100 K and V_M = 10 km³ gives Ξ ≈ 10⁻⁴ for a permeability of K = 10⁻¹² m² (Norton & Knight, 1977). Hence the MCS assumption that percolation can deliver anatectic melt to resident magma rapidly is easily met. Because Ξ ∝ ⁠, the required condition is more easily achieved the larger the resident magma body. At higher temperatures, where wallrock is ductile, stresses can be relieved by thermally activated viscoelastic relaxation that depends on grain size, temperature, pressure, and mineral type and abundance. In this case, a viscous relaxation time τ_η = η^WR/σ can be defined and compared with the Darcy timescale τ_Darcy. Again, the condition for validity of the MCS assumption is that

as fractures can develop in viscoelastic materials. With typical values for the deviatoric stress and wallrock viscosity (σ = 1 MPa and η^WR = 10¹⁸ Pa s), the condition that Ψ is <1 is met. Hence Darcy percolative flow can outpace both magma solidification and crustal relaxation, and the validity of the MCS assumptions is sustained.

The final issue involves the validity of the percolation threshold assumption. As noted above, the efficiency of melt extraction from partially molten wallrock is controlled by the permeability (K) of the wallrock, which in turn depends on the volume, connectivity and topology of partial melt, its viscosity (η_m), and the pressure gradient driving percolative flow. The porosity topology is determined by the melting behavior, lithological texture and grain size of the crystals making up the wallrock restite. The existence of a percolation threshold, defined as the volume fraction of melt at which wallrock becomes permeable, is of special interest in MCS calculations as achieving this threshold is a necessary condition for magma contamination. Once the percolation threshold is reached during partial melting, porous media flow in response to buoyancy and pressure forces will commence and allow magma contamination by percolative flow at rates sufficient to realistically apply the MCS formulation. The MCS is configured such that when the wallrock melt fraction exceeds this threshold, anatectic liquid in excess of the critical amount is removed and mixed with melt in the adjacent magma body. Detailed analysis of the dynamics of melt segregation in multiphase rocks undergoing partial melting is beyond the scope of this study. However, recent developments in the dynamics of solid–liquid systems subject to pressure gradients are directly applicable to the MCS formulation (see Burg & Vigneresse, 2002; Vigneresse & Burg, 2004; Cheadle et al., 2004). Consistent with percolation theory and experimental studies, two critical melt fractions have been identified in solid + melt mixtures as melt fraction varies from zero to unity. Starting at subsolidus conditions, as melt fraction increases the first critical fraction is reached at a melt fraction of f_crit ≈ 0·08 ± 0·04. In silicate systems, because melt density is not very different from crystal density, volume fractions and mass fractions are essentially identical. For example, a volume fraction of 0·08 corresponds to a mass fraction of melt of 0·07, given typical densities. The physical picture is that the threshold melt fraction corresponds to the presence of a melt film along crystal grain boundaries in the solid–melt mixture owing to grain boundary wetting. This film forms a three-dimensional tubular network that accommodates strain partitioning, allowing for the segregation and flow of melt in response to lateral stress gradients and the local distribution of buoyancy and viscous forces. Precise values of this critical melt fraction depend on the distribution of phases in the wallrock at the centimeter scale and cannot be categorized a priori without petrographic information. Presumably, the critical value can itself vary spatially. In the MCS, the critical melt fraction for partial melt mobility is typically set in the range 0·04–0·12, although other values may be justified if additional information is available.

CODE LOGISTICS AND ACCESS

MCS v1.0 utilizes rhyolite-MELTS for all thermodynamic computations. A simulation requires appropriate compositional, mass and thermal input for magma body, wallrock and recharge reservoir(s), which are entered through an interface written in Visual Basic (Microsoft Excel for Mac 2011). Output related to the rhyolite-MELTS results (i.e. major element and phase equilibria, mass and thermal data) is collated in an Excel workbook. Appropriate simulation results are transferred into the trace element and isotope calculator, which is a separate Excel workbook. The user is required to enter solid–melt and solid–fluid partition coefficients, and because full phase information is provided by the rhyolite-MELTS results for each temperature decrement in the magma body, the user has the option of entering different solid–melt and solid–fluid partition coefficients for each temperature step. Bulk partition coefficients are calculated using the well-described bulk partition coefficient equations for solid–melt and solid–fluid. Information about access to the MCS v1.0 is available at http://melts.ofm-research.org/mcs.html.

CONCLUSIONS

1. The Magma Chamber Simulator is a computational tool that allows forward modeling of the geochemical and petrological evolution of a composite system composed of four subsystems: a magma body consisting of melt ± solids ± fluid, a cumulate reservoir, wallrock, and a set of recharge reservoirs. Magma and wallrock exchange material and heat during cooling and crystallization of the magma body. As it heats up from its initial subsolidus temperature, surrounding wallrock may generate anatectic melts that are mixed with magma body melt. An arbitrary number of distinct recharge events is accommodated, allowing realistic modeling of simultaneous RAFC processes self-consistently. Thermodynamic solutions are provided by rhyolite-MELTS, and EC-RAFC informs the trace element and isotope conservation approach. Input parameters include compositions, masses and temperatures for all subsystems and require the specification of initial conditions such as pressure and the ratio of initial magma body melt mass to initial wallrock mass (Λ).

2. An MCS solution provides a continuous record of the composition (major and trace elements, isotopic ratios), masses and temperatures of all relevant phases (melt ± solids ± fluid) within the magma body, cumulate reservoir, recharge magma and wallrock as a function of magma temperature. Thus, the geochemical and petrological trajectory through composition and temperature space during concurrent RAFC can be rigorously explored.

3. Detailed results presented for a case of intrusion of high-alumina basalt into dioritic upper crust illustrate that formulations that do not correctly assess the mass impacts of assimilation may lead to incorrect conclusions. For example, the expectation that typically incompatible elements are more enriched as crustal contamination progresses (compared with an FC-only case) is not necessarily correct. Likewise, three- or four-component phase diagram analysis may not successfully predict the suppression or enhancement of phase stability during AFC. Realistic open-system behavior cannot be correctly predicted by formulations, such as binary mixing, because the chemical and phase evolution of magma and wallrock is governed by the complex interplay of mass exchange among magma melt, cumulate reservoir and wallrock. Rigorous thermodynamic solutions constrained by mass and energy balance are required to accurately map the impact of AFC on melt, crystals, and fluid in wallrock and magma body.

4. Petrologists and geochemists can document igneous suites in exquisite detail using ever-improving technology. Progress in understanding how these remarkable sets of data can be used to make predictions that will continue to shape our understanding of the origin and evolution of magmas, and address societal problems such as volcanic risk mitigation and acquisition of vital natural resources, requires concomitant improvements in computational modeling. The Magma Chamber Simulator is one step toward the goal of an integrated tool that realistically models the complex thermal, chemical, energetic and dynamical evolution of open-system crustal magmatic systems.

ACKNOWLEDGEMENTS

We are indebted to Gene Wilder for inspiration. We are grateful to Paul Asimov, Paula Smith and an anonymous reviewer for helpful reviews. We also thank Executive Editor Marjorie Wilson and Editorial Manager Alastair Lumsden for editorial handling and comments that improved this work.

FUNDING

This work was supported by the National Science Foundation (EAR 0440010, 0810086, 0838182) and Central Washington University.

SUPPLEMENTARY DATA

Supplementary Data for this paper are available at Journal of Petrology online.