| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Journal of Petrology Volume 42 Number 5 Pages 999-1018 2001
© Oxford University Press 2001
Energy-Constrained Open-System Magmatic Processes I: General Model and Energy-Constrained Assimilation and Fractional Crystallization (EC-AFC) Formulation
INSTITUTE FOR CRUSTAL STUDIES AND DEPARTMENT OF GEOLOGICAL SCIENCES, UNIVERSITY OF CALIFORNIA, SANTA BARBARA, CA 93106, USA
Received January 10, 2000; Revised typescript accepted August 17, 2000
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
Geochemical data for igneous rocks provide definitive evidence for the occurrence of open-system processes in magma bodies, including Replenishment by intrusion of primitive magma, Assimilation of anatectic wallrock melt and cumulate formation by Fractional Crystallization. A general class of models (Energy Conserved-RAFC or EC-RAFC) can be constructed, which allow simulation of geochemical paths for trace elements and isotopic ratios for magma undergoing simultaneous replenishment, assimilation and fractional crystallization during the approach to thermal equilibration. The general problem of EC-RAFC is formulated as a set of 3 + t + i + s coupled nonlinear ordinary differential equations, where the number of trace elements, radiogenic and stable isotope ratios simultaneously modeled are t, i and s, respectively. Partial melting of wallrock, modeled as fractional melting, is incorporated, as are sensible and latent heat effects. Temperature-dependent partition coefficients are used to describe trace element distributions. Solution of the set of differential equations, with magma temperature ( Tm ) as the independent variable, provides values for the average temperature of wallrock ( Ta ), fraction of melt within the magma body ( Mm ), mass of cumulates ( Mc ) formed, mass of wallrock involved in the thermal interaction ( Ma°), mass of anatectic melt assimilated ( Ma*), the concentration of t trace elements (Cm ) and i + s isotopic ratios (
m) in magma (melt plus cumulates) and in anatectic melt delivered to the evolving magma body. Input parameters include a user-defined equilibration temperature ( Teq ), the initial temperatures and compositions of magma, recharge magma, and wallrock, distribution coefficients and their temperature dependences, heats of fusion and crystallization, and isobaric specific heat capacities of the various constituents. For a priori defined magma recharge mass and Teq, the mass of wallrock heated to Teq is computed and the geochemical trajectory of melt is determined as magma temperature approaches Teq from its starting value, Tm°. The effects of imperfect extraction of anatectic wallrock melt may be accounted for by introduction of an extraction efficiency factor. Mathematical details of a simpler model, EC-AFC (no replenishment) are provided. Energy-constrained models have the advantage of linking thermal and chemical properties of magma chambers. Compared with ‘classical’ models that conserve mass and species only, they represent more complete assessments of the complex physicochemical dynamics governing the geochemical evolution of open-system magma bodies. Results of EC-AFC simulations demonstrate that geochemical trends can differ significantly from predictions based on ‘classical’ AFC even when recharge plays no role. Incorporation of energy conservation and partial melting into geochemical models allows important coupled effects to play their natural role. KEY WORDS: assimilation–fractional crystallization; geochemical model; isotope; magma chamber; trace element
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
Magmatism represents an important mechanism for transport of energy and matter within the mantle and between the mantle and crust on the terrestrial planets. On Earth, roughly 30 km3 of magma is generated and emplaced into or onto the crust each year (Shaw, 1973
; Crisp, 1984). By this reckoning, a volume of mafic magma equal to 15% of the mantle has been extracted since the beginning of the Archean 4·0 Gyr ago. When allowance is made for magma generation during the Hadean, an inescapable conclusion is that a significant fraction, perhaps exceeding half, of the mantle has undergone at least one cycle of partial melting. The transport of material and heat by magmatic processes has played an important if not defining role in the geological (tectonic and geochemical) evolution of Earth.
Magma acts as an efficient heat-transfer agent. Advection of heat is geologically rapid (e.g. Spera, 1984) and proportional to the amount and rate of magma transport as well as to its calorimetric properties. On a longer timescale, the redistribution of incompatible radiogenic elements (K, U and Th) during partial melting, a process controlled by atomic scale crystal–melt cation exchange reactions, also contributes to heat transfer and the geotherm.
Melts generated in the mantle are the dominant contributors to the formation of continental crust by accretion of subduction-derived arc terranes (e.g. Ernst, 1999) or mafic oceanic plateaux onto pre-existing cratons. Intrusion of mafic magmas beneath, into or onto pre-existing crust (mafic plating) may be an equally important mechanism for crustal growth. The processes of anatexis (partial melting of pre-existing crystalline rock) and fractional crystallization typically leave geochemical fingerprints in an evolving magma body, betraying both chemical and thermal interaction with surroundings. This is especially apparent when magma and host are compositionally and thermodynamically distinct. An aim of modern geochemistry is to interpret these fingerprints in light of petrologic processes and thereby constrain the geologic history of a magma body. This history can then be fitted into the larger tectonic context and afford better comprehension of the thermal and chemical evolution of Earth.
Although many facets of the compositional diversity problem have been studied—indeed, this issue lies at the very heart of igneous petrology—there is no general formulation for open magmatic systems that self-consistently treats energy, momentum and material conservation during country rock partial melting, magma replenishment and intra-chamber cumulate formation. There are several classes of extant models. Thermodynamically based transport models (e.g. Spera et al., 1982, 1995; Sparks, 1986; Ghiorso & Kelemen, 1987; Huppert & Sparks, 1988; J. S. Marsh, 1989; Kelemen, 1990; Barboza & Bergantz, 1997, 1998; Raia & Spera, 1997; B. D. Marsh, 1998; Weinberg & Leitch, 1998) account for energy balance, phase equilibria and rheological effects but generally ignore isotopic and trace element constraints. In transport models where geochemical details are explicitly considered, the physical model is often very specific in terms of process or configuration. Generality is sometimes compromised for the sake of tractability. It is often difficult to recognize when the model applies to real systems. On the other hand, ‘classical’ replenishment, assimilation and fractional crystallization geochemical models (e.g. DePaolo, 1981; O’Hara, 1995, 1998), although characterizing the geochemical evolution of a magma in detail, neglect important constraints imposed by energy conservation, the thermodynamics of partial melting, and the dynamics of melt extraction and magma mixing.
In the model developed here and applied to a number of systems in a companion paper (Bohrson & Spera, 2001), we attempt to bridge the gap between these approaches by combining constraints from energy conservation in partially molten systems with material balance relations for trace elements and isotopic ratios. The general model also incorporates constraints from momentum conservation in an approximate way. The model is zero dimensional in the sense that spatial fields are not explicitly considered. The parameterization involves energy, mass, and isotopic and trace element conservation applied to arbitrary-sized open magmatic systems undergoing concurrent magma replenishment, assimilation of anatectic melts and formation of cumulates by fractional crystallization. The main purpose is to develop and make available to the petrologic community an energy-constrained model for the evolution of trace elements and isotopic ratios in a magma body. Ease of application to natural systems is an important objective. There are thousands of magmatic systems that have been compositionally characterized. A purpose of this study is to provide a general tool, easily applied to such systems, to gain further geological insight.
|
ENERGY-CONSTRAINED OPEN-SYSTEM MAGMATIC PROCESSES: SYNOPSIS |
|---|
|
TOPABSTRACTINTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
Natural systems can be modeled as energy-conserving open composite systems sealed off adiabatically from the rest of their surroundings. The assumption of adiabaticity implies a regional length scale; obviously, it is only at such scales that the adiabatic approximation is valid. Adiabaticity facilitates the analysis of energy-conserved recharge, assimilation, and fractional crystallization (EC-RAFC) processes and allows for application to a wide variety of natural systems. The sub-systems making up the composite system include country rock, magma body, and a reservoir of recharge magma that are separated from one another by diathermal and open (permeable) boundaries that allow complete or partial heat and species exchange. The sub-systems themselves can be viewed as composite systems so that compositionally zoned magma bodies can be treated. Within this hierarchical framework, the energy-constrained formulation presented here for a (homogeneous) magma body consists of a set of 3 + t + i + s coupled ordinary differential equations expressing conservation of energy (enthalpy), total mass, t trace element species and i + s isotopic ratios. The independent variable is conveniently taken as the temperature of melt (Tm) within the magma sub-system, which is composed of melt plus cumulates. Tracking the geochemical path during RAFC demands an accurate rendering of phase relations that govern trace element mineral–melt and fluid–melt distribution coefficients including variations with composition and temperature as well as the thermal properties and melt productivity functions for magma, crust and recharge magma. Results from both experimental phase equilibria and thermodynamic modeling (e.g. MELTS, Ghiorso, 1997
) may be used to constrain the thermodynamics of melting, and these constraints can then be included explicitly into the general EC formulation thereby accounting for major element constraints. The energy-constrained formulation incorporates geological and geochemical information along with thermophysical and thermochemical parameters relevant to a specific magmatic system, thereby providing a holistic view of an evolving magma body. A schematic representation of EC-RAFC processes is pictured in Fig. 1.
|
Unfortunately, accurate accounting of momentum, chemical species and heat exchange during open-system processes demands precise knowledge of the size and configuration of the system in addition to a host of other details unique to a particular system. Development of a general three-dimensional thermodynamic and fluid dynamic model to account for the variability of natural systems is bewilderingly complex and not attempted here. We adopt the more limited goal of establishing a class of models that self-consistently accounts for mass, species and energy conservation within a ‘simple’ composite system framework. Predictions of this model may then be compared with voluminous data from natural systems.
In this paper, the conceptual framework, assumptions and approximations of the general EC-RAFC model are first discussed. Following this, the more specific application, EC-AFC, is presented in quantitative terms. Although recharge is undoubtedly an important process in many natural systems, it is useful to first consider the implications of energy conservation and partial melting to systems dominated by AFC processes. The full differential equations for EC-RAFC will be presented elsewhere. A computer program for performing EC-AFC computations is available at http://magma.geol.ucsb.edu/research/recharge.html. The program may also be obtained from the Journal of Petrology Web site at http://www.petrology.oupjournals.org. In the companion paper (Bohrson & Spera, 2001), the geochemical implications of EC-AFC are discussed and the model is applied to a number of natural systems.
|
ENERGY-CONSTRAINED OPEN-SYSTEM MAGMATIC PROCESSES: CONCEPTUAL FRAMEWORK AND BACKGROUND |
|---|
|
TOPABSTRACTINTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
Historical perspective on RAFC models
Crustal assimilation and reaction of assimilated material with magma together with fractional crystallization have been recognized as important petrogenetic processes for over a century. Studies on the Tertiary volcanic centers of western Scotland (Harker, 1904
; Thomas, 1922; Bailey et al., 1924; Tyrrell, 1928), on Hakone volcano, Japan (Kuno, 1950) and on Paricutin, Mexico (Wilcox, 1954) emphasized the necessity of invoking assimilation of anatectic melts to explain petrographic observations and chemical trends. Early treatments relied mostly on phase equilibrium arguments. The classic treatment by Bowen (1922a, 1922b, 1928) enumerates possible consequences of assimilation of igneous and sedimentary rocks based on phase equilibria and remains vital to this day. Kelemen (1990) reviewed the assimilation–reaction concept with a special focus on isenthalpic constraints and accurate rendering of phase relations by minimization of appropriate thermodynamic potentials. He has argued that assimilation–reaction of peridotite by ascending olivine tholeiitic magmas combined with crystal fractionation can give rise to calc-alkaline and basaltic andesitic liquids so common in island-arc environments.
Attention to the effects of assimilation on trace element and isotopic signatures helped usher in the modern era of study of assimilation. Investigations of Pb isotopes (e.g. Doe, 1967; Doe et al., 1968; Church, 1976) and Sr–Nd isotopes (Carter et al., 1978) underscored the importance of crustal contamination as a mechanism for imparting ‘enriched’ isotope signatures in continental volcanic and plutonic suites. Noting the existence of magmatic systems characterized by radiogenic Sr and 18O enrichments, Taylor and colleagues (Taylor et al., 1979; Taylor, 1980), James (1981) and others (e.g. Hoefs et al., 1980) documented cases where assimilation was clearly at work (e.g. Roccamonfina). The pairing of Sr and O isotopes is especially telling because processes that lead to enrichments in each can be very different. Radiogenic Sr comes from decay of high Rb/Sr sources, typically of continental origin, whereas elevated 18O is commonly the result of low-temperature interaction between rock and fluid, a common upper-crustal process. Thus, identification of high 87Sr/86Sr–18O/16O rocks led those workers to recognize the contribution that crust and crustal fluids have in the isotopic and trace element characteristics of many igneous suites.
The process of crustal assimilation was initially viewed as a two-component problem involving mixing between magma and wallrock melt. Taylor and colleagues (Taylor et al., 1979; Taylor, 1980) underscored the point made by Bowen, that assimilation of crustal material was inextricably linked to the process of fractional crystallization. The lack of superheated magmas implies that an important source of heat for melting wallrock is the latent heat of crystallization (cumulate formation). Taylor et al. (1979) illustrated variations in the oxygen isotope composition of crystallizing magma by calculating an approximate mass of wallrock melted per unit mass of crystallized cumulate, in the absence of recharge. His calculation emphasizes recognition of the energy balance inherent in the process of assimilation–fractional crystallization by identifying the ratio of mass of assimilated anatectic melt to the mass of cumulates formed during an AFC ‘event’. Heat associated with magma cooling as well as heat advected to the magma body by influx of recharge magma can be added to the latent heat generated by magma crystallization.
DePaolo (1981) described a quantitative approach to modeling trace elements and isotope ratios in AFC by solving differential equations expressing conservation of mass and species with fraction of melt (F ) as the independent variable. Energy balance was not explicitly considered in this early form of AFC. To obtain closed-form analytical solutions to the classical AFC expressions, the ratio of the mass assimilated to the mass crystallized (defined by the symbol r) was set equal to a constant. The tractable nature of these solutions was a critical factor in the emergence of numerous, quantitative treatments of AFC by many petrologists.
In recognition of the complex behavior that characterizes magma chambers, the effects of recharge were added to the total mass and species conservation expressions describing AFC, thereby generating total mass and species conservation expressions for AFC augmented by recharge (RAFC; DePaolo, 1985). Other variations of ‘classical’ AFC models have been proposed; among the most useful of these is the ‘multiple element’ approach developed by Aitcheson & Forrest (1994). Noting the weakness of somewhat arbitrarily setting the value of r, those workers solved mass and species conservation equations (including recharge) so that r and 
(defined as the mass of material assimilated/original mass of magma) may be determined. Solutions using a variety of trace element concentrations (or isotope ratios) were then plotted, with the intersection of the r– curves apparently representing the ‘true values’ of these variables. The advantage of this approach is that r is not arbitrarily defined but instead is independently inferred from an over-determined system of mass balance equations. The disadvantages are that the value of r determined by this procedure depends on the particular trace elements used to constrain its value, and no energy conservation constraint is imposed.
Although effort has mostly been directed at investigating solutions to mass and species conservation equations, some attention has been paid to energy balance. For example, Grove & Kinzler (1986) recognized that the amount of assimilant required to explain the composition of Medicine Lake lava flows exceeded the amount theoretically available, based simply on heat available from latent effects. The discrepancy was reconciled in the context of geologic evidence that suggested the proportion of contaminated material is small relative to the total volume of erupted lava. Thus, a large volume of basalt was available as a source of heat to melt wallrock but only a small volume of basalt was chemically contaminated. DePaolo (1985), DePaolo et al. (1992) and Perry et al. (1993) derived equations that consider the importance of recharge, assimilation, fractional crystallization and loss of heat from the magma body by a rough accounting of conductive heat transfer. This was accomplished by defining an effective ambient crustal temperature and assuming the ratio of the rate of heat loss to the rate of cumulate formation was a constant in the range 0·5–1. More recently, Reiners et al. (1995), building on the earlier work of Ghiorso & Kelemen (1987) and Kelemen (1986, 1990), examined the consequences of isenthalpic AFC using the MELTS algorithm. Their model balances the enthalpy released by cooling and solidification against that required for digestion of assimilant. Reiners et al. found that for basaltic magma invading and assimilating felsic crust, assimilation can lead to suppression of crystallization of plagioclase and clinopyroxene, which in turn can potentially lead to large values of r, in some cases >1. After crystallization of plagioclase and clinopyroxene commences, r values can drop significantly, commonly to <1. Their calculation clearly demonstrates that r is not, in general, constant during an AFC ‘event’. By utilizing the DePaolo equations in an iterative way, Reiners et al. further concluded that relatively mafic magmas that do not appear contaminated on the basis of major elements (e.g. <51 wt % SiO2, >6 wt % MgO) can, in fact, be contaminated by continental crust. The trace element and isotopic ratios of such a mafic melt therefore may not be representative of the mantle source from which they were derived. This possibility has important implications for characterizing chemical heterogeneity of mantle source regions based on the composition of basaltic lava, as the usual assumption is that primitive MgO-rich, low-crystallinity melts are not substantially contaminated. One concludes from these considerations that application of an energy conservation principle can have a direct and major effect on the inferences drawn regarding petrologic processes based on geochemical data, which, in retrospect, is not surprising.
Edwards & Russell (1998) traced thermal and chemical changes in a magma body undergoing AFC by determining dissolution rates of wallrock phases and balancing enthalpy. By experimentally determining mineral dissolution rates, those workers determined the rate at which assimilant is added to a magma chamber. As a function of time, they tracked the amount of crystallization, the amount of assimilant incorporated (and therefore r), the temperature change, and the major element variations (liquid lines of descent) in the magma body. They found that, depending on the mineral dissolution rates, high values of r could be achieved, particularly early in the crystallization history of the body. In addition, their analysis predicts the timescale of AFC to be weeks to years.
Magma chamber–country rock: a composite thermodynamic system
A simplified thermodynamic description of a magmatic system and its surroundings serves to outline the conceptual framework of the energy-constrained approach (Fig. 1). A system undergoing EC-RAFC is taken as a composite of three sub-systems: (1) a magma body of initial mass Mo; (2) a finite mass (defined self-consistently based on energy conservation) of surrounding country rock (Ma°); (3) a reservoir of recharge magma that is added to the magma body sub-system during the course of an RAFC ‘event’. The cumulative amount of magma injected during the RAFC ‘event’ is Mr°. The three-part composite system is separated from its surroundings (i.e. the rest of the lithosphere) by a boundary that is closed and adiabatic. The country rock–magma body boundary is diathermal and open, allowing exchange of enthalpy and anatectic melt. The replenishment sub-system reservoir allows recharge magma, of specific composition and enthalpy, to flow into the magma body. The chamber is filled with standing magma (a mixture of cumulates and melt in thermal and chemical equilibrium) that is stirred by the release of chemical and thermal buoyancy, by the settling of dislodged pieces of country rock, and by momentum transfer between recharge and standing magma. The magma body may have a complex compositional structure: regions may be homogeneous, exhibit continuous gradients (e.g. a compositionally stratified magma cupola), preserve compositional discontinuities or horizons, or be mixed at any of a variety of spatial scales [see, e.g. Oldenburg et al. (1989) for quantification of mixing]. Mixing of magmas is favored in cases where steep gradients in temperature, composition and velocity are established. In systems with relatively high recharge rates, magma may become well mixed except within relatively thin marginal chemical and rheological mixed-phase boundary layers, typically several meters to tens of meters thick (Trial & Spera, 1990). The momentum transfer induced by possible magma fountaining during recharge will also aid in mixing. When the mean recharge rate is small, the tendency for mixing is lower, and stratified (thermally and compositionally) magma horizons may develop and persist.
Enthalpy transport from magma to country rock leads to compositional modification of resident melt by fractional crystallization. Enthalpy loss as a result of magma cooling and crystallization will heat country rock and induce anatexis when country rock temperatures exceed the local solidus, which is a function of composition and pressure. The extent of anatectic melt formation depends on the country rock solidus, the amount of heat provided by magma cooling and crystallization, and the amount of heat available from recharge magma. Assimilation of anatectic melt will alter the elemental and isotopic composition of melt within the magma body. By writing statements for conservation of enthalpy, mass, trace species and isotopic ratios, one can compute changes in the composition of melt as each sub-system approaches thermal and chemical equilibrium. Although in nature reaction may proceed between country rock and magma without melting, as illustrated in the work of Kelemen and co-workers (Ghiorso & Kelemen, 1987; Kelemen, 1990; Kelemen et al., 1990), it is assumed in the EC-RAFC model that melting of country rock and the mixing of anatectic liquids with ‘pristine’ magma is a necessary prelude to magma contamination.
Thermal interaction between magma and country rock: thermal heuristics
Perhaps the most fundamental interaction between magma and its surroundings is a thermal one. The main irreversibility involves the production of entropy during heat exchange between all sub-systems as thermal equilibrium is approached. Developing a general set of rules for predicting quantitative aspects of magma body thermal evolution is complicated by many factors. In this section, we discuss some of the more problematic issues as a backdrop for how these matters are dealt with in the energy-constrained models.
As a point of departure, we consider a model detailed by Bejan & Anderson (1983) that serves as an illustration of a coupled country rock–magma body heat transfer system. Similar models are described elsewhere (Spera, 1979; Carrigan, 1988; Bergantz, 1989, 1992; Bowers et al., 1990). In this conjugate heat transfer problem, the ambient temperature of the country rock ‘far away’ from the contact (Ta°) and the magma temperature in the interior of the magma body (Tm°) are the characteristic temperatures of the problem. The average country rock–magma boundary temperature (Tb) is determined by the relative vigor of hydrothermal convection in the country rock (which transports heat away from the contact) compared with magma convection within the magma body itself (which transports heat to the contact). The ratio B = kaRd1/2/kmRa1/4 (all parameters defined in Table 1) is the dimensionless parameter that governs the boundary temperature, Tb. When B is large, hydrothermal convection efficiently carries away heat brought to the boundary and, as a consequence, the boundary temperature remains relatively low. That is, in the limit of large B, the country rock–magma body boundary acts as a good thermal conductor and Tb 
Ta°. An environment where this condition might prevail is in the upper crust in a region undergoing extensional shear failure as a result of tectonic forces where a vigorous hydrothermal system is established. In contrast, for small B, magmatic heat cannot be transported rapidly away from the country rock–magma contact and hence Tb Tm°. In this case, the boundary is more like an adiabatic one than a diathermal one. In low-permeability catazonal rocks (20–40 km depth), hydrothermal flow may be less vigorous (e.g. Schoofs et al., 2000). In this case, the boundary margin temperature is controlled by the balance between heat convected to the magma chamber and heat conduction into country rock. In the general case, the boundary temperature Tb is given by
 |
0 and -1/2 as B 50. Values of f(B) at intermediate values of B may be found in table 1 of Bejan & Anderson (1983).
|
It is important to note that the extent of magmatic contamination by country rock is not a simple function of the boundary temperature but depends also on the initial country rock temperature (Ta°), melt productivity functions fa(T) and fm(T), and thermodynamic parameters such as the heat capacity and phase change enthalpies of relevant materials.
Magmatic–hydrothermal systems: transport modeling
Knowledge of the dynamics of magmatic–hydrothermal systems has expanded enormously in the past quarter-century. The role of chemical as well as thermal convection in magma dynamics has been explored by a variety of analytical, laboratory and numerical studies. The effects of crystals and bubbles on magma rheology and convection, as well as inferences of magma residence time based on particle (crystal and bubble) size distributions and radiogenic isotope systematics have also been addressed. Multiphase effects have been considered by incorporation of phase relations in model two- or three-component systems with the equations of motion, energy and material balances to account for marginal mush layers along magma body–country rock contacts, the generation of chemically induced buoyancy within such layers, and the growth of porous cumulates (Campbell & Turner, 1987; Huppert & Sparks, 1988; Oldenburg & Spera, 1991, 1992; Bergantz & Dawes, 1994; Brown et al., 1995; Spera et al., 1995; Barboza & Bergantz, 1997, 1998; Raia & Spera, 1997). Similarly, models of hydrothermal convection have also progressed in the last quarter-century. Complexities associated with medium anisotropy, heterogeneous permeability fields, salinity and thermally driven buoyancy, the effects of topography, variations in the thermodynamic and transport properties of H2O–CO2–NaCl fluids, and the role played by metamorphic fluids from prograde reactions have been addressed. Based on extensive geochemical investigations, Taylor and colleagues (e.g. Criss & Taylor, 1986, and references cited within) have pointed out important differences in the behavior of meteoric–hydrothermal systems depending upon magma composition and temperature.
It is clear from collective study of these works that any attempt to track the chemical evolution of magma bodies that neglects thermal budgets and the energetics of partial melting is inchoate. The extent of chemical interaction and chemical evolution of magma is inextricably bound to the thermal interaction amongst the composite sub-systems. The irreversible production of entropy implied by heat exchange between sub-systems in RAFC interactions allows one to define a temperature, the equilibration temperature (Teq), as a natural ‘progress’ variable to describe the approach to equilibrium during the irreversible process of heat exchange (Prigogine, 1962). Because the eruptive temperature of magma can be directly measured or estimated from geothermometry, it is sensible to utilize magma temperature as the independent variable in EC-RAFC modeling.
Extraction of anatectic melt
How efficiently is anatectic melt extracted from its source host and added to the standing magma body? This question is central to petrogenetic models that attempt to predict the geochemical fingerprint of anatexis. Heat and material exchange can be decoupled if anatectic melts are generated but not mixed into the magma body sub-system. Hence it is essential to discuss the efficacy of melt extraction before detailing EC-RAFC and EC-AFC models.
The processes whereby melt is extracted from partially molten host rock are complex and remain imperfectly known despite considerable research in the past quarter-century (Arzi, 1978; van der Molen and Paterson, 1979; Shaw, 1980; McKenzie, 1984, 1985; Richter & McKenzie, 1984; Fowler, 1985; Ribe, 1985, 1987; Richter, 1986; Sleep, 1987; Scott, 1988; Bergantz, 1989; Scott & Stevenson, 1989; Stevenson, 1989; Riley et al., 1990; Riley & Kohlstedt, 1991; Takahashi, 1992; Hart, 1993; Bergantz & Dawes, 1994; Daines & Kohlstedt, 1994, 1997; Sawyer, 1994; Brown et al., 1995; Petford, 1995; Rushmer, 1995; Rutter & Neumann, 1995; Kohlstedt & Zimmerman, 1996; Connolly et al., 1997; Barboza & Bergantz, 1998; Petford & Koenders, 1998; Zimmerman et al., 1999). There are two general models for the description of melt segregation: percolative and melt-vein. These are ‘end-member’ models and reality is likely to incorporate elements of both. In percolation models, the matrix or residual solid is considered a granular, deformable material of isotropic and uniform permeability. Melt is segregated from the deformable matrix at rates that depend on the balance amongst capillary, gravity and viscous forces. The compaction length 
c (McKenzie, 1984) gives the length scale over which the deformation must occur in order that Darcy pressure balances compaction forces:
 |
*, K and a represent the effective viscosity of the matrix, the permeability (a function of grain size and porosity) and the viscosity of anatectic melt. If the region of partial fusion is significantly larger than c, melt segregation is governed by the balance between buoyancy and Darcy friction. The permeability is notoriously difficult to estimate. Experiments by Riley et al. (1990) suggest an expression of the form K = do2fan/b where do is a typical grain size, fa the melt fraction, n 
1 and b 
5 x 103. For values appropriate to crustal anatexis ( fa 0·2, do = 3 mm, g = 9·8 m/s2, 
= 300 kg/m3, * 1018 Pa s and a 105 Pa s), K 4 x 10-11 m2 and c 20 m. The relative vertical velocity between melt and matrix (u) is then given by
 |
represents the density difference between solid and melt. The velocity of melt relative to restite is of the order of 150 m/Ma. According to the percolative model, a thin (tens of meters on average) selvage of country rock could be drained on a 105 year timescale. In the mantle where the melt viscosity is rather low (a 1 Pa s), the percolative model gives extraction rates 105 times higher, essentially geologically instantaneous. Short-lived radioisotope studies are consistent with such high melt extraction efficiencies for partial melting of peridotitic sources to generate basaltic magma.
The melt-vein model for melt segregation, perhaps more applicable to melting of non-peridotitic sources, gains support from field, laboratory and fluid dynamical studies. This scenario may be more geologically realistic in light of the heterogeneous and anisotropic nature of most source regions. The essential insight is that melt accumulates in discrete veins or pods, of the order of centimeters in length, oriented parallel to the axis of minimum effective stress during shear deformation associated with intrusion of partially molten material. There is clear field and laboratory evidence that extensional shear failure is associated with melt mobilization (e.g. Shaw, 1980; Spera, 1980). Once veins grow large enough, an interconnected critical drainage network may develop so that melt can be transported at rates of 1 mm/s (30 km/yr) or greater in response to small differential stresses. These rates are much greater than those characterizing Darcy percolative flow. Laboratory experiments on crystal–melt suspensions often exhibit small-scale melt redistribution in response to surface forces and small deviatoric stresses (van der Molen & Patterson, 1979; Cooper & Kohlstedt, 1986). The total volume of melt that can exist in a drainage network depends upon the volume of the partially molten region and the details of the distribution of crack apertures, lengths and widths, as well as the orientation of the local stress field including but not restricted to thermal stresses. Presumably each magma–country rock interaction is unique. Once a critical level of connectivity exists amongst melt-filled fractures, relatively rapid drainage of the fracture plexus can occur as the effective (equivalent) permeability dramatically increases. Petford & Koenders (1998) proposed a phenomenological model for self-organization and fracture connectivity created by thermal stresses generated by intrusion of hot magma into cooler ambient country rock. They estimated transient permeability in the range 10-10–10-5 m2, many orders of magnitude greater than in the percolative model. The permeability (K) was related to the fracture connectivity (
) by Gueguen & Dienes (1989) in a fracture model according to
 |
1/
3) is the fracture number density (where is the fracture spacing), and c is the fracture length. Theory allows estimation of , which changes dramatically as the ratio of connected to unconnected nodes in the percolation network varies from zero to unity. As a scale permeability, let us consider A 0·1, 0·5 m, c 0·1 m and 0·2. This gives K 10-8 m2, considerably larger than the classical Darcy estimate (K = do2fan/b) of 10-11 m2.
For the purposes of this study, we parameterize the extraction of anatectic melt from its source (country rock) by defining a melt extraction efficiency (
) defined as the fraction of melt delivered to the magma body relative to mass generated in the country rock (Ma*). Dimensional analysis shows that depends on the melt viscosity, magnitude and orientation of deviatoric stresses, dimensions of typical cracks (apertures, crack lengths and widths) and the crack number density. Wholesale digestion of stoped country rock (Walker & Kiefer, 1985) is another way anatectic melt and hydrothermal fluids can enter a magma body. Further discussion of the efficiency of anatectic melt extraction may be found in the studies by Ashworth & Brown (1990) and Rubie & Brearley (1990). Wickham (1990) has made some important arguments regarding evidence from Sr and O isotopes and their redistribution during anatexis. Vielzeuf & Holloway (1988) discussed melt production in pelitic sources. Hart (1993) concluded that melt generated by anatexis of upwelling peridotitic mantle experiences only limited diffusive re-equilibration, on the basis of a fractal network–melt channel model similar to the one discussed here.
Because of the difficulty in assessing the extent of melt extraction in a general way, is set equal to unity in the EC-AFC expressions below and in the examples presented by Bohrson & Spera (2001). By setting = 1, the geochemical effects of contamination are maximized, as all melt produced during anatexis acts as a chemical contaminant and not just an energy sink.
|
EC-RAFC AND EC-AFC MODEL FORMULATIONS |
|---|
|
TOPABSTRACTINTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
EC-RAFC model formulation
In the EC-RAFC model, a mass of magma Mo at initial temperature Tm° is injected into crust at temperature Ta° (Fig. 1). Both magma and crust have unique initial trace element (Cm°, Ca°) and isotopic (
m°, a°) compositions. Recharge magma of mass Mr with distinct initial trace element (Cr°) and isotopic (r°) composition and initial temperature (Tr°) is added to original magma according to some a priori defined relation Mr = f(Tm) such that total mass injected during the RAFC ‘event’ equals Mr°. The constraint imposed by enthalpy balance is that heat associated with cooling and crystallization of original magma and recharge magma is matched with the enthalpy requirements for heating and partial melting of country rock. The amount of country rock involved in the interaction (Ma°) uniquely follows upon specification of the equilibration temperature, Teq, which is set a priori. Country rock partial melts are partially or completely extracted, added to and thermally homogenized with the molten portion of the standing magma body. Conservation of enthalpy enables one to track the mean temperature of country rock restite (Ta) as a function of the melt temperature (Tm) during the approach to the common equilibrium temperature (Teq). At each step in the approach to Teq, local thermal equilibrium is maintained among cumulates, magma in the chamber, recharge magma, and the newly added mass of assimilant. This means that assimilated country rock melt and recharge magma are brought to the local and common temperature Tm. The calculation is continued until Tm and Ta equal Teq.
Latent heat effects are treated by adopting a single effective fusion (hfus) or crystallization (hcry) enthalpy computed as a weighted average for the mixture of phases involved in the phase transition. Values of fusion enthalpies for a number of phases are given in Table 2. It is assumed that enthalpy is released uniformly throughout the melting or crystallization interval, which is defined by the temperature difference between the liquidus (Tl,m or Tl,a) and solidus temperature (Ts). This is an approximation; representative calculations accomplished using the MELTS algorithm (Ghiorso & Sack, 1995; Ghiorso, 1997) indicate that crystallization intervals do not vary widely during AFC. Sensible contributions to the enthalpy balance are computed using average specific isobaric heat capacities calculated from available thermodynamic data (Table 3). Although the temperature dependence of the isobaric heat capacity of silicate materials is well known, this level of detail is likely to be of second-order importance.
|
|
The solution of the balance equations representing conservation of energy, mass and species provides self-consistent values for the mass of heated country rock (restite) plus partially molten country rock involved in the RAFC process (Ma°), the amount of assimilant partial melt (Ma*) (some portion of Ma°), the mass of melt in the chamber (Mm), the total amount of recharge magma (Mr) and the mass of cumulates (Mc) as a function of the magma temperature Tm along the path Tm° Teq. Simulation results also include the trace element (Cm) and isotopic (m) composition of contaminated magma as well as the trace element composition (Cc) of cumulates. Because of the hierarchical nature of the EC-RAFC formulation, compositionally distinct melt reservoirs within the melt portion of the magma body can be accounted for. However, in most cases so little is known regarding the structure of melt compositional heterogeneity that this level of detail seems unwarranted.
The equilibration temperature, Teq, is used to parameterize a particular RAFC evolution in geochemical space. Teq characterizes the overall extent of interaction between magma and country rock and is related to Tb discussed earlier. For Tb approaching Tm°, the Teq of the RAFC ‘event’ is relatively high, and the mass of country rock (Ma°) involved in the RAFC process is relatively small (although the fraction Ma*/Ma° may be large). In this case, the thermal effects of injected magma do not extend great distances from the contact into country rock because of poor effective heat conductivity. In contrast, if Tb is low (e.g. rapid hydrothermal heat advection in an epizonal environment), Ma° is larger but a smaller fraction Ma*/Ma° is formed. In the case Teq 
Ts, the country rock temperature (Ta) does not exceed the solidus (on average), and assimilation of anatectic melts is not important. An appropriate Teq can often be estimated by geological considerations, thus linking the model with geologically relevant data. For example, in reconstruction of petrogenetic processes in a volcanic system characterized by a spectrum of compositions, one would naturally choose Teq less than or equal to the eruptive temperature of the most evolved composition. In contrast, if the magmatic system under study is characterized by limited compositional variation, Teq may be close to but somewhat less than Tm°. Teq is used as the independent parameter in the EC-RAFC model because it provides an intuitive measure of the extent of thermal interaction of country rock and magma. Additionally, the irreversible production of entropy, generated mainly by dissipation of magmatic heat, is directly related to Teq.
EC-AFC model formulation
In the remainder of this paper, we consider only the process of assimilation–fractional crystallization; recharge has been set identically to zero to highlight the constraints imposed purely by energy conservation and partial melting on assimilation of anatectic melts. Future work will delineate the mathematics and applications of EC-RAFC. In the following sections, the differential equations expressing conservation of energy, mass and species as well as initial conditions and parameters of the model are presented. A derivation of these equations along with a description of the algorithm used to model an EC-AFC ‘event’ is presented in the Appendix. All symbols are defined in Table 1.
EC-AFC mathematical model
The description of the EC-AFC model comprises two parts. The first is an integral calculation defining a relationship between Teq and Ma° for a particular set of thermal parameters (heat capacity, latent heats, liquidii of magma and country rock, melting functions) and initial conditions (Ta°, Tm°). The second describes the energy-constrained thermal and geochemical path followed by a batch of magma during approach to Teq. For a particular set of thermal parameters and initial conditions, each Teq is associated with a unique Ma°, which is required to solve the path-dependent differential equations describing the chemical and thermal evolution of the magma body.
Integral balance. The integral enthalpy balance provides an essential constraint on the geochemical path followed as initially pristine melt cools, crystallizes and assimilates anatectic liquid. The balance can be represented as an algebraic relationship between the mass of country rock involved in the thermal interaction (Ma°) and the equilibration temperature (Teq), the independent variable in the EC-AFC formulation. It should be noted that not all of Ma° may actually melt, although all of it does eventually come to thermal equilibrium (Teq) with magma. The extent of country rock melting depends on the temperature–melt productivity relation, fa(Ta) = Ma*/Ma°, where Ma* is the mass of country rock that partially melts and is incrementally added to the batch of pristine magma, of initial mass Mo, undergoing simultaneous fractional crystallization. Energy conservation incorporates heating and partial melting of country rock, magma cooling and formation of cumulates self-consistently. The parameters include the melt production functions for country rock and magma ( fa and fm, respectively), enthalpy of crystallization for magma (hcry), the fusion enthalpy for country rock (hfus) and the average isobaric specific heat capacity of magma (Cp,m) and assimilant (Cp,a). The relationship may be expressed in terms of Ma° as a function of Teq:
 |
) simulations and are chosen to account for differences in bulk composition independent of the EC-AFC model per se.
Geochemical and thermal path. To compute the thermal, trace element and isotopic path of melt within the magma (melt plus cumulate) sub-system, a system of coupled ordinary differential equations expressing conservation of enthalpy, mass, species and isotopic composition is solved. The temperature of the country rock restite (Ta), the mass of assimilant that partially melts (Ma*), the fraction of melt (Mm) within the magma body, the mass of cumulates (Mc) crystallized and the concentration of trace elements in melt (Cm) and crystalline cumulates (Cc) are computed self-consistently. The isotopic compositions of radiogenic elements (m) and of oxygen (m) in the melt are determined as well. The independent variable is the melt temperature (Tm), and the calculation ends when the equilibration temperature (Teq), set a priori by the investigator, is attained. We emphasize that Teq must be specified, to compute the path in temperature–composition space, because Ma°, the mass of country rock involved in the EC-RAFC ‘event’, is a function of Teq. Although the reason for setting the calculation this way may appear arcane, a windfall arises from the flexibility of the scheme when applied to real systems, which, after all, are highly variable in conformation.
The model consists of 2 + t + i + s differential equations, where t is the number of trace elements, i the number of radiogenic isotopic ratios and s the number of stable isotopes considered in the calculation. There are no limitations on t, s or i. Pressure is not accommodated explicitly in the current generation of EC-AFC models although adjustment of liquidi and the solidus temperatures can be made to simulate variations in pressure. A full treatment of major phase relations assuming local thermodynamic equilibrium would improve the predictive capabilities of the model.
The first differential equation expresses conservation of energy along the path as country rock heats up, partially melts and thermally equilibrates with melt from the magma body. Energy conservation leads to a differential expression for the country rock restite temperature (Ta) as a function of the magma temperature Tm. Anatectic melt, generated as country rock temperature (Ta) exceeds the solidus (Ts), is assumed to rapidly mix both compositionally and thermally with standing magma. It should be noted that along the equilibration path, restite (i.e. the crystalline part of the country rock) is not immediately brought into thermal equilibrium with magma whereas anatectic melt is brought to local thermal equilibrium with standing magma. It is possible to allow for incomplete extraction of anatectic melt by introduction of an extraction efficiency factor (). In this study, we adopt the approach that all anatectic melt finds its way into the magma body (i.e. = 1). For cases where Teq is less than Tl,a, the thermal equilibrium condition is Ta = Tm = Teq. When Teq > Tl,a, the equilibrium condition is that Tm = Teq, as no restite remains (i.e. all wallrock has been melted). The differential equation gives the derivative of the restite temperature Ta with respect to the magma temperature Tm along the path towards thermal equilibrium at Teq:
 |
The second constraint is a statement of conservation of total mass. The derivative of the mass of melt in the magma body (Mm) with respect to the magma temperature is expressed in terms of melt production functions and their derivatives and uses the constraint imposed by energy conservation. The mass of melt in the magma body along the path is related to the amount of melt initially present before the AFC ‘event’ (Mo), the amount added by assimilation of anatectic melt (Ma*) and the amount removed by cumulate formation (Mc). The expression has the differential form
 |
Conservation of species provides the basis for defining the mass concentration of a trace element in the melt as a function of the melt temperature, Tm. Concentration of the trace element in the country rock is denoted by Ca. We assume that partial melting of country rock is described by fractional melting so that the concentration of trace element in anatectic melt is given by
 |
 |
For an isotopic ratio m, the differential equation becomes
 |
a° and m represent the isotopic ratio in anatectic melt (identical to country rock) and in the melt (at temperature Tm), respectively. Radiogenic in-growth has been neglected and no isotopic fractionation upon melting or crystallization is allowed.
Finally, the differential equation expressing oxygen isotope balance may be written
 |
. In cases where magma and country rock have nearly the same oxygen isotopic ratio, temperature effects may be important to consider and (11) should be modified.
Initial conditions. Equations (6), (7), (9), (10) and (11) constitute an initial value problem with Tm as the independent variable. This set of 2 + t + i + s coupled ordinary differential equations is subject to the following initial conditions: at T = Tm°, Ta = Ta°, Mm = Mo, Cm = Cm° for t trace element species, m = m° for i isotope species including m = m° for oxygen. The system of differential equations is numerically solved by a fourth-order Runge–Kutta method once cast into dimensionless form (see the Appendix).
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
The geochemical evolution of a batch of magma is intimately tied to interactions with its environment. Open-system behavior (e.g. assimilation and magma replenishment) is more the rule than the exception as magma moves between source and surface. The open-system processes of magma replenishment, assimilation and fractional crystallization can be formulated as a coupled set of 3 + t + i + s non-linear ordinary differential equations with magma temperature as the independent variable where t is the number of trace elements, i is the number of isotopic ratios and s is the number of stable isotope ratios treated in the simulation of the geochemical path. In this paper, explicit formulation of EC-AFC has been derived. Inclusion of an energy conservation principle constrains the trace and isotopic evolution of magma system more tightly than ‘classical’ AFC models because the initial thermal state of the surroundings and the thermodynamic properties of silicates are self-consistently accommodated. In the accompanying paper (Bohrson & Spera, 2001
) the EC-AFC formulation is applied to a number of natural systems. Detailed quantitative comparison between ‘classical’ AFC and EC-AFC geochemical paths highlights some remarkable and significant differences in trace element–isotope ratio coordinates for many instances. The single most important and in some sense obvious conclusion of this work is that the energetics of petrologic processes must be explicitly accounted for in study of complex RAFC phenomena. RAFC without self-consistent application of energy conservation is like an automobile without an engine. A user-friendly program that runs on a PC or a Mac and requires Excel 98 is available at http://magma.geol.ucsb.edu/research/recharge.html. |
APPENDIX: DERIVATION OF EC-AFC EQUATIONS |
|---|
|
TOPABSTRACTINTRODUCTIONENERGY-CONSTRAINED OPEN-SYSTEM...ENERGY-CONSTRAINED OPEN-SYSTEM...EC-RAFC AND EC-AFC MODEL...CONCLUSIONSAPPENDIX: DERIVATION OF EC-AFC...REFERENCES |
|---|
Integral energy balance
The model posits that the enthalpy required for heating and partial melting of country rock (the heat absorbed) is balanced by the heat liberated by cooling and crystallization (cumulate formation) of magma. It is necessary to select a final or equilibration temperature (Teq), to define an EC-AFC ‘event’. Once Teq is picked, application of energy conservation provides a relationship between the total mass of country rock (Ma°) involved in the magma–host thermal interaction and the chosen equilibration temperature, Teq. Enthalpy available for heating and partial melting comes from the sensible heat associated with cooling of magma of mass Mo from initial temperature Tm° to Teq and the latent heat (
hcry) liberated by formation of cumulates of mass Mc = [1 -fm(T)]Mo. The enthalpy of crystallization and isobaric heat capacity are taken as constants and the temperature–melt fraction relation (also known as the melt productivity), fm(T), is a function of temperature and is obtained by independent phase equilibria experiments or MELTS modeling. We have found the MELTS simulator to be especially handy for determining liquidi and solidi and bulk distribution coefficients, Da and Dm, for different starting materials. The integral expression for heat liberated is
 |
Enthalpy is absorbed by country rock of mass Ma° as it heats up and undergoes partial fusion, provided the local country rock temperature (Ta) exceeds the solidus (Ts). The initial temperature of country rock is Ta°, and the entire mass is brought to the equilibration temperature Teq. The fraction of partial melt in the country rock is defined by the melt production function fa(T), which is analogous to fm(T) defined above. The mass of anatectic melt, a quantity that grows during the EC-AFC ‘event’, is defined as Ma*. The heat of fusion (hfus) required for melting is taken as a constant. The integral expression for the heat absorbed by country rock becomes
 |