Документ взят из кэша поисковой машины. Адрес оригинального документа : http://geo.web.ru/~ariskin/docs/papers/p06_071.pdf Дата изменения: Tue May 16 18:12:05 2006 Дата индексирования: Sat Dec 22 04:31:32 2007 Кодировка: ISO8859-5

ISSN 0016-7029, Geochemistry International, 2006, Vol. 44, No. 1, pp. 72-93. ? Pleiades Publishing, Inc., 2006. Original Russian Text ? A.A. Ariskin, A.A. Yaroshevsky, 2006, published in Geokhimiya, 2006, No. 1, pp. 80-102.

Crystallization Differentiation of Intrusive Magmatic Melt: Development of a Convection-Accumulation Model
A. A. Ariskin* and A. A. Yaroshevsky**
* Vernadsky Institute of Geochemistry and Analytical Chemistry, Russian Academy of Sciences, ul. Kosygina 19, Moscow, 119991 Russia ** Division of Geology, Moscow State University, Vorob'evy gory, Moscow, 119992 Russia
Received July 7, 2005

Abstract--The paper summarizes the principal results obtained over the past three decades at the Vernadsky Institute and the Department of Geochemistry of the Moscow State University by the computer simulation of basaltic magma differentiation in magma chambers. The processes of diffusion-controlled mass transfer in a chamber are demonstrated to be principally limited by the heat resources of the cooling magma and cannot play any significant role during the large-scale partitioning of melt components. The leading mass-transfer mechanism is the settling of crystals from convecting magma in the form of suspension flows that are enriched and depleted in the solid phase. The physical prerequisite for the onset of this concentration convection is the existence of boundary layers, which are characterized by volume crystallization and a gradient distribution of the suspended phases. Considered in detail are the principles used in the development of algorithms with regard for the occurrence of a boundary layer and the "instantaneous" stirring of the crystallizing magma that does not hamper the settling of mineral grains forming the cumulus. The plausibility of the convection-accumulation model is illustrated by the example of the reconstructed inner structure of differentiated Siberian traps. In application to these bodies, it is demonstrated that the solutions of the forward and inverse simulation problems with the use of geochemical thermometry techniques are identical. This is a convincing argument for the predominance of convection-accumulation processes during the formation of thin tabular magmatic bodies. The further development of the computer model for the differentiation dynamics should involve the processes of compositional convection related to the migration and reactivity of the intercumulus melt. DOI: 10.1134/S0016702906010083

INTRODUCTION The problem of basaltic magma differentiation in magmatic chambers and the origin of layered mafic- ultramafic intrusions is pivotal in magmatic geochemistry and petrology because of the following two reasons. The "compositional" aspect of this problem is related to the scale of basaltic magmatism and the parental role of mafic melts as a source of magmatic associations. The genetic significance of layered intrusions is determined by their internal structures, first of all, the systematic spatial variations in the compositions of the rocks and minerals, which provide information on the magmatic evolution, mechanisms responsible for the differentiation of material in the chamber, and the physicochemical parameters under which the massifs crystallized. The multifactor character of intrusive magmatism manifests itself as the superposition of and interaction between various processes of mass transfer in the magma and country rocks, with these processes occurring within broad ranges of thermodynamic parameters and in different geodynamic environments. This makes it extremely difficult to systematically analyze differentiation in magmatic chambers and partly explains why these processes are still discussed with reference to different, often mutually exclusive, schemes and concepts.
72

A constructive way out of this situation, which stems from the multifactor character of the problems, is evident: this is the application of a mathematically formalized theory of transport processes to magmatic systems. This approach enables the transformation of hypothetic "scenarios" of differentiation into quite rigorous physicochemical models, which can be correlated with observational results in terms of a great variety of mineralogical and geochemical parameters [1]. In formulating the mathematical basis of such models, it is necessary to bear in mind that the chemical differentiation of magmas can take place only if the following two conditions are satisfied. First, the material affected by the transformations should be heterogeneous (consist of several phases), with the magmatic melt serving as a transporting medium. The second condition concerns the differences between the chemical compositions of the coexisting crystals and melt, which is a thermodynamic prerequisite for crystallization differentiation. A combination of these factors inevitably manifests itself during the cooling of magmatic masses, and this is why mechanical motions in the form of the Stokes settling of crystals, the infiltration and expelling of an intercumulus melt, or the movements of suspension flows can generate a spatial


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

73

chemical heterogeneity in an originally homogeneous magmatic system. This conclusion, which is not questioned any more by most researchers dealing with the physics of magmatism, did not seem as obvious in the early 1970s, when the first computer simulations of the crystallization differentiation of basaltic magmas were conducted at the Vernadsky Institute. The first results of this research were published in a series of papers by Frenkel and Yaroshevsky [2-5] (Fig. 1). These papers dealt with the theoretical aspects of the cooling and crystallization of magmatic systems, presented the principles formulated to develop models involving the settling of solid phases, and considered the thermal and chemical dynamics of the origin of layered stratiform intrusions. In the wake of this research, Koptev-Dvornikov et al. [6] analyzed data obtained by imitating the inner structures of thin (200 m) intrusions according to the mechanism of the Stokes sedimentation of crystals and arrived at the conclusion that most of the characteristics of the modeled bodies and the principal mineralogical features of intrusive traps are consistent. The results thus obtained were of pioneering character and made the foundation for the modeling of the more realistic convection-accumulation differentiation mechanism [6, 7]. The construction of a numerical model for the accumulation of crystals from a convecting magma and the possibility of comparing the petrological-geochemical consequences of this process for various starting and boundary conditions heralded the beginning of a qualitatively new stage in studying crystallization differentiation. The analysis of concepts concerning the genesis of layered intrusions was underlain by a quantitative basis [8]. Starting in the mid-1900s, the researchers arrived at the consensus that the phase and modal layering of mafic-ultramafic complexes was produced by the origin and evolution of the cumulus, which filled the magmatic chamber from below and whose composition systematically evolved from high- to low-temperature mineral assemblages (see the reviews [8-10]). The variations in the rock and mineral compositions were thought to result from fractional crystallization that occurred simultaneously with the efficient stirring of magma in the chamber. At the same time, the concepts concerning the mechanisms responsible for material differentiation, the effect of the thermal regime, and the origin of cumulus minerals remained principally different. The discussions were mostly centered on two hypotheses: crystal settling and in-situ crystallization. According to the former concept, crystalline material nucleates mostly near the intrusion roof, the crystals sink within the main magma volume, and are then redistributed into the lower portion of the chamber to form layers of cumulus minerals [11-13]. The concept of in-situ crystallization denied that extensive gravitational separation can play any significant role and focused on processes near the bottom of the chamber or in the vicinity of its walls [14-16].
GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

Fig. 1. M.Ya. Frenkel (left) and A.A. Yaroshevsky developing a model for the crystallization differentiation of intrusive magmas. Vernadsky Institute of Geochemistry and Analytical Chemistry, Russian Academy of Sciences, 1978.

It is pertinent to recall that alternative hypotheses of differentiation were formulated when magmatic chambers were interpreted as "pools" of a practically homogeneous melt, which contained no suspended crystals and was originally overheated with respect to its liquidus. This brought about the opinion that the settling of minerals can occur only during homogeneous (volume) nucleation of the crystalline phase, which is, in fact, impossible because of the difference between the adiabatic gradient (controlling the onset of temperature convection [14, 17]) and the actual baric dependence of the magma temperature.1 This view on the rheology and temperature profiles of intrusive magmas underlie the concepts of melt fractionation at the predominant role of crystallization near the lower solidification front or at the walls of the chamber [15, 19]. This scenario of differentiation implied that the bulk of the crystals should be formed in thin boundary layers that separate the front of complete solidification from the main volume of the magmatic liquid. It was postulated that the crystallizing solid phase remains mechanically connected with the crystallization zone (in-situ growth), and the chemical composition of the magma is modified by diffusion exchange with the residual melt adjacent to the surface of the crystallizing material. RESULTS OF THEORETICAL ANALYSIS AND THE FIRST NUMERICAL MODELS FOR MAGMA DIFFERENTIATION: 1976-1979 The concept of in-situ crystallization explained some features of the compositional evolution of minerals in intrusive bodies, whose analogues were found in technological and metallurgical processes [20]. Because of this, the first steps in the development of a
1

This physical effect was used in developing the model of the zone refining of the Earth [18].

2006


74 (a) E C

ARISKIN, YAROSHEVSKY (b)
0

(c)
s FA

B TE B+L

A Solidus T Boundary layer I x II III
s FA = 0

Solid phases (A + B) A+L Melt (L) TB

TL TA

Melt

T

x

Fig. 2. Schematic structure of a boundary layer near the complete solidification front [3]. (a) Modeled diagram for a binary A-B system with eutectic point E (TL is the liquidus temperature of the initial melt of the composition C0); (b) representation of the phase heterogeneity of the boundary layer; (c) distribution of solid phase A in the vertical section (results of theoretical analysis [3]).

realistic model for differentiation in magma chambers involved the theoretical analysis of the diffusion mechanisms responsible for material differentiation in magmatic chambers and the significance of the diffusion mechanisms of material differentiation in magmatic processes [2]. This problem was solved in application to a plane-parallel chamber, which approximated intrusive traps. The successive analysis of heat transfer and concentration diffusion at the boundaries of a tabular intrusive body made it possible to evaluate the proportionality coefficient ki that relates the diffusion flow of melt component j within the boundary layer ( I M , Fick's law) with the heat flux through the crystallization front (JQ, Fourier's law) I
j M j

affected the further research even more significantly than the conclusion about the inefficiency of diffusionrelated mass-transfer mechanisms. The possibility of volumetric (homogeneous) nucleation near the solidification front invalidated the main arguments of opponents of the settling of crystals and made it possible to start the examination of other physical schemes that allowed for the relative motions of the solid phases and melt. These concepts were developed in a later paper by Frenkel and Yaroshevsky [3] concerning the relations between the structure of the upper boundary layer and the conditions under which convection flows can be formed and serve as a factor of mass transfer for the crystalline phase. The first things to be analyzed were the consequences of the delocalization of the nucleation temperature near the chamber roof during the early crystallization stages. This problem was solved analytically, in application to a binary system of the eutectic type (A-B) for a parental melt C0 that can crystallize within a finite temperature range, from a liquidus temperature TL to an eutectic temperature TE (Fig. 2a). A simplified (linearized) form of the phase diagram made it possible to derive a semiempirical equation to relate the mass of crystals (of phase A) newly formed in the boundary layer with the temperature-compositional parameters of the original melt and the eutectic. The analysis of the dependencies thus obtained led to the conclusion that the initial crystallization of the "magma" between the liquidus and solidus proceeds unevenly and is, in fact, restricted to the vicinities of two points: TL and TE. Hence, the whole crystallization zone can be subdivided into three parts: (i) a relatively thin layer with the "freezing" of the low-temperature eutectoid material; (ii) an intermediate unit with insignificant variations in the melt crystallinity; and (iii) a significantly heterogeneous layer x in which the temperature reaches the liquidus values and the amount of the suspended solid phase decreases to zero (Fig. 2c). The second and
Vol. 44 No. 1 2006

= kj JQ,

(1)

where kj is close to ~10-7 mol/cal. This led to the conclusion that the transfer of components by the mechanism of concentration diffusion is controlled by thermal fluxes, and the possibility of efficient fractionation in compliance with the in-situ crystallization scheme is principally limited by the amount of heat "stored" in a magma with a low kj [2]. This is the fundamental difference between the spontaneous solidification of magmatic melts and technological processes based on in-situ crystallization, in which the main liquid volume receives heat from an external heater [21]. Another important result of this analysis pertains to the relations between the amount of crystals remaining at the crystallization front (heterogeneous nucleation) and formed within this layer (homogeneous nucleation). It was determined for an actively stirred melt that only an insignificant portion of this melt (a few percent of the overall amount of the crystals) can crystallize in situ at interfaces with the host rocks, and the rest of the solid phase is produced within the boundary layer and is distributed over the volume of the residual liquid in the form of suspended crystals [2]. This result

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

75

third units of the boundary layer are characterized by the homogeneous nucleation and growth of the crystalline phase. The most heterogeneous interval x, in which the temperature varies from TL - T to TL, is of principal importance (Figs. 2a, 2b) because, first, it is responsible for the origin of the bulk of the crystals whose composition is the closest to equilibrium with the liquidus melt, and, second, the boundary layer occurs in immediate contact with the main volume of the magmatic reservoir and controls the development of concentration convection flows [3]. The point is that the appearance of crystals suspended in melt (even in relatively thin layers) causes a gradient distribution of density near the roof and walls of a magmatic chamber. This, in turn, gives rise to the hydrodynamic instability of the magmatic system and results in the extension of the density heterogeneity due to suspension flows (flows of suspended crystals) and/or the sinking (or ascent) of the crystals relative to the liquid [3, 8]. The superposition and interaction of these two principal types of convection flows in the chamber controls the general structure of the magmatic body and the variability of manifestations of magmatic layering [22]. M.Ya. Frenkel emphasized that, regardless of the predominant mechanism of the transport of the crystalline phase to the bottom of the intrusion, the thickness of the high-temperature part of the boundary layer affects the convection structure and the whole process of magma differentiation in the chamber. The possible thickness of the transition layer was first evaluated in [3] for an approximation of the problem of the heating of a semi-infinite medium [23] x
T

t/ T L ,

(2)

where t is the time, and is the thermal conductivity of the host rocks. For the situation T = 10Аы (which corresponds to the difference TL - TE of about 100Аы), = 5 з 10-3 cm2/s and TL 103Аы, it was determined that x becomes close to 10 cm within one year after the emplacement and increases to ~1 m over the following 100 years. It is known from the hydrodynamic theory that the horizontal size of convective cells (the diameter of the ascending and descending flows) is commonly close to the thickness of the layer in which the density inversion is generated. Because of this, the obtained estimates led to the conclusion about the small-scale character of concentration convection during the initial solidification stages of intrusive bodies. During this stage, the convective structure of the crystallizing magma should look like as an array of thin heterogeneous flows that merge near the roof of the intrusion [3, 8]. It is reasonable to expect that magma cooling should be associated with the widening of the convective cells and the transition to the regime of large suspension flows.
GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

Of course, this analysis and the evaluations were of provisional character. They were not extended to multicomponent systems and involved a limited list of thermal parameters that were available in the early 1970s. However, it should be mentioned that these studies attracted attention to details in the structure of the boundary layer and the importance of this structure for the differentiation mechanisms of magmatic melts. Over the past two decades, the analysis of the systematic distribution of the crystalline phase in transitional layers was one of the key problems of the dynamics of magmatism [24-28]. The plausibility of the predictions in [3] follows from the data of the direct observations of lava lakes at Hawaiian volcanoes. The results obtained by the drilling of the chill crusts and the direct sampling of layers (which are up to 5-8 m thick) in immediate contact with magma confirmed the continuous character of the variations in the temperature (~920-1130Аы) and the compositions of the mesostasis and glass (76- 50% SiO2), as well as the significant variations in the amounts of phenocrysts in the melt (from a few to a few dozen percent [29-30]). Hence, the conclusion about the possibility of homogeneous nucleation near a crystallization front has been confirmed experimentally. The answer to the question as to whether crystals that appear at the upper contact can sink in the magma depends on the ratio of the thickness of the contact crystallization zone x (Fig. 2c) and the overall crystallization rate within the boundary layer. To solve this problem, it is not sufficient to elucidate the physical nature of the boundary layers or to determine their general inner structure in dependence on the cooling regime. The systematic analysis of the dynamics of magma differentiation required the transition to spatiotemporal models with regard for the simultaneous displacement of the crystallization fronts and the arrival of the crystalline material to the volume of the convecting magma. The complex composition of naturally occurring melts and the necessity of taking into account phase transitions left practically no possibility of using analytical relations describing the processes of mass and heat transfer (useful reviews of these problems are given, based on the data of foreign researchers, in [10, 31]). The difficulties in the solution of even some of the required differential equations (for example, during the consideration of the thermal problems [9]) were not compensated by the realism of the simulations in terms of their comparison with the structure of natural magmatic bodies. This situation called for the construction and examination of a numerical model for the thermal history and differentiation dynamics of magma in a chamber. SEDIMENTATION MODEL The difficulties were related to the fact that magmatic petrology thirty years ago had no experience in solving such complicated problems. It was necessary to correlate the physical concepts concerning the driving
2006


76 Pl

ARISKIN, YAROSHEVSKY

1400АC
1200
1200

TE = 1110АC II I

1300

1400АC

Ol

Px

Fig. 3. Crystallization diagram for the modeled olivine- pyroxene-plagioclase system that was used to develop the early models for basaltic magma differentiation in a chamber [4-6]. Crystallization proportions (mol %): Ol-Pl (35/65); Ol-Px (8/92); Px-Pl (40/60); Ol-Pl-Px (3/60/37). Circles correspond to the compositions of the "magmatic melts" used to test the sedimentation and convection-accumulation models. Solid circles correspond to the calculation variants displayed in Figs. 5a (composition I) and 5b (composition II).

forces of magma differentiation with reconstructions of the phase diagram for basalt crystallization, the motion laws of solid phases with the conditions of their accumulation, and the calculated temperatures with the evaluated proportions of minerals and melt in different parts of a magmatic chamber. Because of this, the initial stage of the development of the model involved certain principal assumptions that made it possible to formalize the relations between the dynamic and thermodynamic parameters without distorting the general relations between the physical processes and their influence on the structure of the intrusive bodies [4, 5]. The first magmatic body selected to be modeled was a relatively thin (h0 = 200 m) tabular intrusion of "basaltic" composition, which was hosted by sedimentary rocks and occurred at a depth of 1-2 km. This choice implied subsequent comparison with differentiated traps in eastern Siberia, whose simple structure and clear geologic setting clarified some problems related to the initial conditions of their origin and the dependence of the contamination of the host rocks on the shape of the intrusion. The other simplification concerned the phase diagram of basalt crystallization. The calculations were conducted with linearized formalism of the Di-An-Fo system, whose principal points (the melting points of the end members and the compositions of the cotectics and eutectic) were corrected with regard for the anchieutectic character of trap magmas (Fig. 3). The diffusion mechanisms of mass transfer

were taken into account when the basic algorithm was tailored but were then "turned off" and rejected from further considerations because of the low efficiency and low petrogenetic significance for the actual simulations [2]. The principal physical assumptions were related to the settling mechanism of crystals and the procedure for evaluating the effect of the boundary layer. During the first stage of the development of the algorithm, the accent was placed on the examination of the Stokes mechanism of crystal settling [4, 5]. This model (which was named the sedimentation model) allowed for the movements of the solid phases in the liquid but not for the convective constituent related to the possibility of the development of suspension flows. Although this approach reduced the significance of thermal and concentration convection, it was an inevitable step in studying the principal effects of the sedimentation of solid phases as a constituent of convection flows. The paper [4] presents a detailed description of the system of partial differential equations that characterize the thermal history of a tabular magmatic body with regard for the multicomponent composition of the melt (Fig. 3) and the possibility of crystal settling. To solve these equations using a computer, a finite difference scheme was proposed that involved the subdivision of a 200-m-thick magma layer into 100 intervals (Fig. 4), with the balance of the heat and material in- and outflows monitored in each of them. The corresponding time scale was subdivided into a number of conjugate segments, and this made it possible to trace the evolution of the temperature of the magma and its modal and chemical composition over the vertical section of the chamber in the course of the cooling of the body. A principal feature of this scheme was the recognition of a transitional layer, which was adjacent to the upper crystallization front (Fig. 4) and modeled the occurrence of a boundary layer in which homogeneous nucleation takes place. The marginal conditions of this model were related to the more definite, thermal part of the problem of the cooling of a magmatic tabular body (see above) and were commonly not changed when a series of simulations was conducted. It was principally interesting to assay whether it is possible to vary the starting parameters, including the starting composition of the parental magma and its crystallinity and the settling rates of crystals of different minerals. Fifteen simulations were conducted with BESM-4 and BESM-6 computers to model the solidification of "trap magmas" of different compositions (Fig. 3) at varying settling rates of Ol, Pl, and Px [5, 6, 8]. A typical example of the structure and temperature field of a crystallizing magma at a given moment in time is presented in Fig. 4. As can be seen, the system consists of a number of layers with mobile phase boundaries, which is manifested in the fact that the mineral assemblages and their proportions vary from point to point with time. These characteristics and the amount of melt were calculated depending on the bulk composition of the subsystem and the temperature
Vol. 44 No. 1 2006

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

77

(heat content) within each of the intervals. The crystallization of the last melt portions in the cooling magma and cumulus layers corresponded to the complete solidification of the intrusive body, a stage when the solving of the thermal problem was terminated. The final result of these simulations was the structure of a complex magmatic body with a certain set of features of modal layering (Fig. 5). The significance of these simulations can hardly be overestimated. First of all, it was demonstrated for the first time that crystals homogeneously nucleating in a boundary layer and sinking in the liquid should dissolve at the front of the underlying higher temperature melt (Fig. 4). However, because of the general cooling of the system, the dissolution of minerals in the melt cannot preclude the systematic displacement of the "cloud of the crystalline phase" toward the chamber bottom and the origin of cumulates [5]. The most significant consequence of the accumulation of crystals is the progressive enrichment of the rocks in olivine towards the inside of the intrusive body from its lower contact and the depletion of this mineral during the development of the main cumulate sequence ("layered series"), as well as the return to the original content near the upper contact (Fig. 5a). This distribution shows features of the so-called S-shaped profile of the variations in the modal composition of rocks, which is typical of most layered intrusions and differentiated sills [25, 27, 31]. The obvious outcome of the research was the establishment of the fact that the changes in the assemblages of cumulus phases in the vertical sections of the modeled bodies were controlled by the crystallization sequence of the parental magma. For example, the crystallization sequence of composition I (Fig. 3) was Ol Ol + Px Ol + Px + Pl and predetermined the appearance of a "peridotite" unit (Ol + Px) and the origin of "anorthosite" with approximately 80 vol % Pl in the upper part of the magmatic rock succession (Fig. 5a). This corresponded not only to the crystallization sequence but also to the specified low velocity of plagioclase sedimentation (0.4 m/yr) relative to that of pyroxene (4.0 m/yr). For composition II (Ol Ol + Pl Ol + Pl + Px, Fig. 3), the appearance of the Ol-Pl assemblage was followed by the appearance of cumulus pyroxene (modeled gabbroids), in spite of the fact that plagioclase sank more slowly than pyroxene. Another principal point was related to the fact that the sedimentation of crystals was a more powerful mechanism of heat exchange with the environment than conductive heat transfer [5, 8]. A detailed analysis of the petrological consequences of the sedimentation model was presented by KoptevDvornikov et al. [6]. The comparison of the results of the computer simulations with data on the structure of differentiated sills in the Siberian Platform indicates that the models were able to reproduce the principal features of the intrusive traps. These analogies included olivine accumulation in the lower units, the presence of an upper and a lower contact zone, and the enrichment
GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

JQ

t Host rocks Boundary layer S L + Ol + Px + Pl L + Ol + Px

t + t

T

L + Ol

L + Ol + Px S Host rocks JQ x

Fig. 4. Typical distribution of temperature and modal composition in the vertical section of a cooling intrusion in which crystal sedimentation is not coupled with convective stirring [4]. The left-hand column demonstrates the principle of the fragmentation of the initial melt layer, as well as the solidified and host rocks, in intervals for which the mass and heat balance problem is to be solved. The right-hand column demonstrates the variations in the phase composition of the magma owing to a heterogeneous temperature distribution. S are the zones of complete solidification.

of low-melting components in the upper part of the vertical section. These features pointed to the determining role of crystal settling as a factor of magma differentiation in a chamber. The problem was the much more extensive layering of the modeled bodies, a feature that was expressed in the contrasting mineral compositions of the rocks at the transitions from one cumulus unit to another (Fig. 5a). This feature persisted regardless of the crystal settling velocities assumed in the models: s even at low j (from a few tenths to a few meters a year), the boundaries between the modeled cumulates were sharp and step-shaped. The simulation results implied that the contrasting variations in the modal composition were related to the heterogeneity of the temperature profile and the localization of phase boundaries in the vertical section of the intrusive magma (Fig. 4). This distribution was predetermined by the differences in the sinking velocities of crystals relative to the liquid and by the absence of global stirring in the magmatic system. The way out of this situation was thought to be provided by accounting for the
2006


78 Height, m 200 Ol 150

ARISKIN, YAROSHEVSKY

Px

(a)

Pl

UCZ

Ol + Pl + Px + L 100 Ol + Pl + L 50 Ol + L LCZ 200 Ol 150 Ol + Pl + Px + L 100 Px (b) Pl UCZ

50 Ol + Pl + L Ol + L LCZ

0

20

40

0

20

60

80 vol %

100

Fig. 5. Distribution of mineral phases in the vertical section of a modeled intrusion that was "produced" by (a) the sedimentary and (b) convection-accumulation mechanisms [6, 8]. The starting simulation conditions were as follows: the thickness of the body h0 = 200 m, and the body occurred at a depth of 2 km. Starting compositions (mol %): (a) 50% Pl, 43% Px, and 7% Ol; (b) 58% Pl, 30% Px, and 12% Ol (see solid circles in Fig. 3). None of the magmas contained intratelluric phenocrysts. Liquidus temperatures: (a) 1179АC, (b) 1149АC. Sedimentation velocity of mineral crystals (m/yr): Pl--0.4, Px--4.0, Ol--7.6. Complete solidification durations: (a) 302 years, (b) 240 years. Successions of modeled "cumulates": (a) Ol Ol + Px Ol + Pl + Px; (b) Ol Ol + Pl Ol + Pl + Px. Contact zones: UCZ--upper, LCZ--lower.

convective (hydrodynamic) constituent of the solid phase transport [5, 6]. CONVECTION-ACCUMULATION MODEL OF MAGMA DIFFERENTIATION It was mentioned above that the physical prerequisite for the onset of concentration convection in magma is an inversion of the density of the medium, which arises as a consequence of melt crystallization in boundary layers, first of all, near the roof of an intrusion. The density heterogeneity leads to the development of differently oriented convection flows (of material relatively enriched and depleted in the solid phase) in the magma. This situation was first qualitatively described by Hess [12], who analyzed the conditions under which convective flows were induced in the Stillwater magma. Hess noted that

the cooling of the boundary layer and its enrichment in suspended crystals can give rise to "blebs" of heavier liquid, which should then sink and entrain the magma layer from the vicinity of the chamber roof. It was not until forty years had elapsed that these ideas were verified by hydrodynamic calculations. Trubitsyn and Kharybin conducted the numerical simulation of convection in a two-phase system and determined that a cloud of an evenly sinking denser phase should be quickly modified into plumes of "crystalline mush," whose integral velocity becomes higher than the velocity of the Stokes settling of individual crystals [32].2
2

The authors of [32] named this sinking regime of denser material not concentration convection, but sedimentation convection, which can lead to misunderstanding. In [4-6, 8] and in this paper, the term sedimentation model is applied exclusively to the Stokes mechanism of crystal settling relative to the liquid. Vol. 44 No. 1 2006

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

79

The inevitable development of vertical heterogeneous flows in magma was also evident for the authors of [3-6], but the algorithm for the simulation of convective stirring was then far from evident. The systematic solution of this problem should have involved at least a two-dimensional computation scheme that would have taken into account ascending and descending material flows. The structure and sizes of the circulation cells should be controlled by the gradients of the mechanical stress, temperature, and chemical potentials, which are related to the reaction and modal transformations of components at all points of the medium and at any moment of time [8]. It is known that the simulation of hydrodynamic models even for systems not complicated with chemical transport and phase transitions requires the use of the most advanced computers. The calculation capability of computers available in the 1970s could hardly allow the systematic simulation of convection in multicomponent systems with phase transitions, and these simulations were practically inefficient. The long times of the calculations needed to construct such algorithms would have not allowed the systematic inspection of a model whose results are comparable with natural data. This situation called for the modification of the sedimentation model in a way that made it possible to take into account (within the scope of the one-dimensional calculation scheme) the principal mineralogical and geochemical consequences of stirring in a plane-parallel chamber. The one-dimensional model "with convection" was implemented under the following two assumptions. One of them is underlain by the presumption that the sinking velocities of crystals relative to the liquid are low compared with the average displacement velocity of the heterogeneous melt (see above). The other assumption followed from the theoretical analysis of the physical nature of the concentration convection [3, 8]. A principal feature of cooling magmas is the circulation character of their convective flows. This is explained by the absence of additional heat sources and by the cooling of the magma throughout the whole perimeter of the magmatic chamber. The origin of concentration flows under these conditions cannot be of stationary character, in spite of the continuous generation of a solid phase near the boundaries of the magmatic reservoir. Descending heterogeneous flows can appear only when a certain critical density of the material is reached, along with a critical amount of crystals in the boundary layer (Fig. 2). This delay gives rise to a nonstationary character of convection, which can be visualized as a succession of fading cycles of the origin and sinking of suspension flows during the continuous cooling of the magma (Fig. 15 in [8]). For a numerical model of differentiation, a principal issue is the duration of these circulations. Assuming that the number of convection cycles corresponding to a cellular subdivision in terms of time is much greater than the number of intervals specified for the model, one can use the approximation of the infinitely high velocity of convecGEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

tion flows. This is equivalent to the "instantaneous" homogenization of the magma and elucidates the meaning of the second assumption utilized to construct the model "with convection." The high efficiency of convection processes implies the homogeneity of the temperature and phase composition of the melt at any moment of time. In this situation, the problem of the assessment and simulation of general convection can be reduced to the summing of the bulk compositions and heat contents of the distinguished melt layers (Fig. 4) and the subsequent calculation of the weighted mean parameters (including the temperature) for a single common convection layer [6]. It should be thereby taken into account that convective stirring does not hamper either the passage of the sinking crystals from the melt into precipitate or the capture of the magmatic mush by the advancing crystallization front. This conjugacy of concentration convection and crystal sedimentation is the pivotal point of the convection-accumulation model for differentiation, which was justified in [6-8]. The genesis of rocks during this combined regime was examined with the use of the same basic algorithm that was utilized in the sedimentation model [4, 5]: the changes involved only an additional procedure of the averaging of the physicochemical characteristics of the magma after each cycle of nucleation and the sedimentation of the solid phase. This approach to the simulation of the settling of crystals oversimplified the physical scheme of mass transfer within a convecting melt layer and proved to be efficient as a tool for examining the petrological consequences of the model. The principal goals of such simulations are not the reconstruction of the hydrodynamic characteristics of the medium but the study of the spatial structure of magmatic bodies that could be formed with the participation of different magma differentiation mechanisms [1]. The plausibility of the model "with convection" was confirmed by the first simulation results for thin tabular magmatic bodies [6]. When these results were compared with the structure of Siberian traps, it was determined that natural bodies show more similarities to bodies modeled in compliance with the convection- accumulation mechanism (Fig. 5b) than to layered bodies simulated as produced by the simple sedimentation of crystals (Fig. 5a). For example, if the factor of stirring is taken into consideration, the transitions between different cumulates become smooth and gradual. The final distribution of minerals in the vertical section (the succession of cumulus phases) is also controlled by the crystallization sequence in the transitional layer, and the proportion of cumulus phases depends on the ratios of the settling velocities of the crystals of different minerals [5, 8]. Another corollary of general convection is an elevated heat flux from the intrusion. This is associated with a relatively quick leveling of the temperature profile in the inner parts of the chamber and the simultaneous shift of magma characteristics toward the eutectic, in our case, the stability field of the Ol-Pl-Px
2006


80

ARISKIN, YAROSHEVSKY

three-mineral assemblage at 1100АC (Fig. 3). The possibility of the convective transport of the solid phase is also reflected in the fact that the modeled bodies bear no "anorthosite" units, in spite of the notably lower sedimentation velocity of plagioclase than those of mafic minerals. Conversely, much plagioclase is distributed in the mass of olivine-plagioclase "cumulates" and "gabbroids" (Fig. 5b), which is consistent with the cotectic (Ol + Pl) nature of the tholeiitic magma of trap provinces. MODELING OF THE INNER STRUCTURE OF DIFFERENTIATED SILLS IN THE SIBERIAN PLATFORM: 1985-1989 The conclusion that crystal settling from a convecting magma plays a leading role in relatively thin tabular bodies made it possible to proceed from examining the petrological consequences of the convection-accumulation mechanism to the numerical simulation and reconstruction of the structure of natural intrusions such as differentiated sills. The determining role in the further progress of this research belonged to two issues. One of them was fieldwork at intrusive traps in eastern Siberia. A group of researchers from the Vernadsky Institute and the Department of Geochemistry of Moscow State University participated in three expeditions to the Podkamennaya Tunguska River in 1969-1976 and to the upper reaches of the Vilyui River in 1980. These expeditions were aimed at the systematic searches for and detailed sampling of as full as possible a succession of dolerite sills, whose thicknesses varied from a few dozen to 160 m. The expeditions provided unique materials on the structure of these tabular intrusions, some of which were originally described as undifferentiated. A few sills were sampled over their whole vertical sections, from the lower to upper contacts, with sampling sites spaced 3-5 m apart. These data were then utilized in simulations that verified the convection-accumulation model [7, 8]. The other aspect was related to the beginning of the systematic computer simulations of basalt crystallization [33-35], which was possible thanks to the accumulation of a representative massif of experimental data on the distribution of components between the main rock-forming minerals and a mafic silicate liquid. Experience in the calibration of more and more precise mineral-melt geothermometers [36-38] stimulated the application of these temperature dependences to the construction of multicomponent silicate systems with phases of variable composition (see the review of models [39]). The thermodynamic principles of the development of computer programs for the calculation of the equilibrium and fractional crystallization of basaltic magmas were described in [33, 39]. The development of this approach resulted in a realistic fractionation model, which was later named COMAGMAT [40, 41]. By the mid-1980s, prerequisites were created for the development of a physicochemical model for magma

differentiation in a chamber that could be applied to the simulation and reconstruction of the conditions under which trap magmas of different compositions solidify and evolve. The principal difficulties concerned the techniques for a combination of the modeling procedures for convectional mass transfer and the algorithm for the calculation of liquid lines of descent for natural systems. The capabilities of the finite-difference scheme [4] were limited by the complexity of the calculation of the phase characteristics of basaltic magmas as functions of their enthalpy. It should be mentioned that the first achievements in the development of the thermal models involving the sedimentation of crystals were underlain by the simple form of the crystallization diagram (Fig. 3), which enabled the derivation of an algebraic equation to relate the temperature and modal composition of a system to the stored heat for each interval of the cell subdivision [4, 5]. This made it possible to maximally efficiently utilize the capabilities of computers during multiple (tens and hundreds of thousands) repetitions of the procedure calculating the proportions of crystals and melt in the mixture. The earliest and later versions of the COMAGMAT program package were based on other principles. The problem of thermodynamic equilibrium was solved here for a known proportion of crystals and melt but not for a given heat content of the system [39, 40]. This modified the strategy for the further development of the convection-accumulation model. The accent was now placed on the assessment of the chemical and mineralogical effects of differentiation, which was made possible with the use of the COMAGMAT model. The evolution of the physical characteristics of processes in the chamber were taken into account by means of analytical approximations. EQUATIONS OF DYNAMICS OF INTRA-CHAMBER DIFFERENTIATION In 1982 Frenkel [42] proposed a system of equations for the process of magma cooling connected with the transport and deposition of the crystallizing phases in compliance with the convection-accumulation mechanism. These equations specified the relations of the heat fluxes across the crystallization fronts with the thermodynamic and dynamic parameters, including the accumulation rate of the crystalline material, the modal composition of the cumulus, and the thickness of the "rock" units produced during discrete time intervals. The derivation of these equations for the dynamics of the convection-accumulation model was based on the results of simulations by means of the successive solution of the thermal problem [6, 8]. The most important consequence of these studies was the reproducibility of the physical characterization of the situation near the upper crystallization front. Figure 6 demonstrates the evolution of the positions of the complete solidification and crystal accumulation boundaries for an intrusive body whose cooling is conVol. 44 No. 1 2006

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT t1 t2 t3 t4 Time t5 t6

81

1.0

Complete solidification

Magma 0.5

Cumulus

Complete solidification Relative height, m

0 1.0

t1

t2

t3

t4

t5

t6 VI

VII V

L + Ol 0.5

L + Ol + Pl

L + Ol + Pl + Px

Ol + Pl + Px + L

Ol + Pl + Px + L

IV

Ol + Pl +L Ol + L 0

Ol + Pl +L III II I

Fig. 6. Movement dynamics of the solidification and crystal accumulation fronts during magma differentiation according to the convection-accumulation mechanism [4-6]. The starting moment of time t0 = 0 corresponds to the emplacement of a magma layer with a homogeneous distribution of temperature and suspended (intratelluric) phases. The time intervals subdivide the cooling and differentiation history into the following stages: t0-t1--development of chilled zones; t1-t2 --opposite movements of the upper and lower crystallization fronts; t2-t3-t4--the upper front stops moving, and cumulus layers are filled (the chamber is filled with a crystalline precipitate); t4-t5-t6--closing cooling stage of the crystal-melt mush. The plots at the bottom illustrate the evolution of the chamber structure (vertical sections) with time. By the time of complete solidification t6, the vertical section can be subdivided into seven zones: I, VII--chilled zones; II--olivine "porphyrites" (LCZ); III--Ol-Pl "cumulates"; IV--Ol-Pl-Px "cumulates; V--sandwich unit; VI--UCZ.

trolled by concentration convection. This plot, which was repeated in all of the simulations, demonstrates that the structural evolution of a tabular intrusion can be subdivided into three stages: (i) the development of chilled zones (a brief time period from the emplacement time to t1), (ii) a brief period with the opposite shifts of the upper and lower crystallization fronts ((t1-t2), and (iii) a long-lasting (with a practically zero velocity) movement of the upper front and the simultaneous thickening of the crystalline-sediment layer at the bottom of the chamber (t2-t3-t4). If the possibility of conGEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

vection flows in the cumulus is ignored (we will return to this assumption below), the final cooling stage (t4-t5-t6) is not interesting from the geochemical viewpoint. By the moment of time t4, the chamber is filled with cumulus, and, thus, if the body is cooling quickly enough (which applies to bodies 100-200 m thick [5]), no large-scale differentiation is possible for it any more. Hence, the petrologically most important stage is the third one (convection accumulation proper), when the bulk of the material crystallizing near the roof occurs in the magmatic melt and gradually sinks at the lower front. This
2006


82

ARISKIN, YAROSHEVSKY

stage corresponds to the most complete fractionation of the melt and the development of the layered series, which composes the bulk of the body. It is interesting that the cessation of the upper front movement is not assumed in the model but is its "response" to the possibility of crystal settling [5, 8]. This asymmetry is one of the principal differences between the convection-accumulation process and other petrological schemes that do not allow for the transport of the crystalline phase [24, 25] or that relate the sedimentation of the solid phase with the crystallization kinetics in the boundary layer [26]. Relations between the Thermal and Dynamic Characteristics The asymmetric movements of the solidification fronts (Fig. 6) and the general convection condition for a magmatic chamber warrant the application of analytical approximations to the modeling of the cooling history of an intrusive body [42]. On the one hand, the convective stirring of the system does not require the solution of problems related to heat transfer within a melt layer and allows the researcher to constrain the considerations to contacts with the host rocks. On the other hand, the almost unchanging position of the upper crystallization front during the main differentiation stage implies the possibility of using the solution of the problem of the heating of a semiinfinite medium with a constant (anchieutectic) temperature maintained at its boundaries [23] J
UF Q

enthalpy of the magmatic melt. This is the initial melt with an intratelluric phase and a small amount of accumulated crystals early in the course of the cooling process (over interval t0-t1-t2) and the intercumulus liquid over the interval t2-t3-t4. The situation is principally different near the roof of the chamber. The overall solid-phase flux at the upper front also depends on the melt crystallization heat IUF = J
con UF UF Q

/l (g/s),

(5)
dir UF

but is subdivided into the fluxes of directional ( I and volume ( I ) crystallization [42] I
UF

)

=I

dir UF

+I

con UF

.

(6)

The physical meaning of this subdivision follows from the characteristics of the main differentiation dir stages (Fig. 6): the value of I UF controls the velocity of the directed growth of the upper contact zone (due to con crystals mechanically linked to the roof), and I UF characterizes the rate of the flux of crystalline material into the volume of the convecting magma. The approximation I
dir LF

= ILF,

I

con LF

=0

(7)

= UF / t ,



UF

= cr r T

0

r/,

(3)

where J Q (cal/s) is the heat flux through the upper crystallization front, t is time, cr is the heat capacity (cal/g grad), r is the density (g/cm3), r is the thermal conductivity of the host rocks (cm2/s), and T0 = const is the "eutectic" temperature. In application to models with a simplified phase diagram and a constant eutectic temperature, Eq. (3) provides an accurate enough evaluation of the heat flux. The value of T0 of naturally occurring basaltic systems can vary and decreases with time. We did not attempt to take this factor into account (neither did we allow for the variations in the thermal parameters) and deliberately introduced certain simplifications in the heating history of the intrusion and focused mostly on the modal and chemical consequences of crystal settling. An analogous approximation was applied to the lower crystallization front J
LF Q

UF

is valid for the lower crystallization front over any time interval. The fluxes of directed crystallization determine the velocities of the movement of the upper (UF) and lower (LF) solidification fronts inward the magmatic chamber, which are related to the crystallinity of the initial (or residual) magma (Fmag)
UF

=I

dir UF

/ l S ( 1 - F

mag

),

F

mag

=


j=1

m

f

mag j

(8)

and the amount of crystals accumulated in the lower cum zone ( F LF )
LF

= I LF / l S ( 1 - F
dir

cum LF

),

F

cum LF

=


j=1

m

f

cum j

,

(9)

where i is the density of the melt (g/cm3), f j is the volume fraction of crystals (of mineral j, 1 j m) suscum pended in the magma, f j is the fraction of minerals in the "cumulus," and S is the unit section surface.
mag

= LF / t ,



LF

= k UF ,

(4)

where k < 1 is the proportionality coefficient with regard for the cooling of the intrusion predominantly through its upper contact. The flux of material crystallizing in situ at the lower solidification front is caused by the crystallization

Relations between the Dynamic and Thermodynamic Parameters The values of f j and the average crystallinity of the magma Fmag are extensive parameters, which reflect the P-T parameters of the mixture of crystals and melt
Vol. 44 No. 1 2006
mag

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

83

in the volume of the convecting layer. They can be calculated for any moment of time with the COMAGMAT program. The fraction of crystals captured at the lower solidification front can be evaluated by relations following from the mass balance of components at the moving boundary [42] f
s cum j

cumulates, the value of for f
cum LF cum j mag j

LF

in Eq. (10) is substituted
cum l

=f

(1 - j/
s

cum LF

),

f

= 1-F

cum

,

(12)

where Fcum is the total fraction of cumulus phases, and f l is the porosity of the cumulates. During the final stage of the process, Fmag can become greater than or equal to F LF (Fmag F LF ), which means the equalization of crystallinity in the residual magma and cumulates. This situation is equivalent to the filling of the chamber with a limiting possible mixture of crystals and signals for the termination of the simulations. Hence, the sets of thermal and dynamic equations that make it possible to evaluate the modal compositions of the rocks produced at the upper and lower crystallization and accumulation fronts are now determined for all major differentiation stages in a chamber (Fig. 6). If data on the compositions of the coexisting minerals and melt are available, one can calculate the concentrations of major and trace elements in the modeled rocks. For the upper crystallization front, which is formed by means of directed growth, the concentration of component i is determined by the equation
cr cr cum

=f

mag j

( 1 - j / LF ) ,
s

(10)

where j is the average sinking velocity of a crystal of a given mineral. Inasmuch as j and
s cum j mag j s j

LF

have different

signs, the value of f >f if > 0. Hence, relations (8)-(10) make it possible to simulate the phase composition of rocks corresponding to the stage of the opposite movements of the upper and lower crystallization fronts. During chilling and the opposite movements of the solidification fronts, these rocks can be referred to as cumulates only provisionally, because the amount of the in situ crystallizing melt is usually cum greater than the amount of the cumulus phases f l F LF . As the intrusive melt cools, the velocity of the crystallization fronts decreases, and the amount of phases suspended in the magma increases. This results in an increase in the concentration of cumulus minerals cum near the lower front f j , a feature that corresponds to the origin of the lower branch of the S-shaped distribution of mineral and geochemical characteristics in intrusive bodies [8, 39]. At a certain moment of time, the value of F
cum LF cr LF , cum

C

i UF

=f

mag l

C+

l i


j=1

m

f

mag j

Ci .

j

(13)

The composition of the rocks originating at the lower solidification or accumulation fronts can be estimated by relations that involve the amounts of minerals and the porosity of the cumulus C
i LF

reaches the maximum possible critical value F which correspond to the beginning of the origin of "normal" cumulus with a high density of crystalline phases. Thereby, the accumulation surface of crystals becomes detached from the lower crystallization front, and, as a consequence, the magmatic chamber is then relatively quickly filled with crystalline precipitate (the third stage, corresponding to the interval t2-t3-t4 in Fig. 6). To quantify the velocity of the accumulation front, Eq. (10) should be summed for crystals of differcum cr ent minerals with regard for the constraint F LF = F LF
cum LF

=f

cum l

Ci +

l


j=1

m

f

cum j

Ci .

j

(14)

=


j=1

m

f

mag s j j

/( F

cr LF

-F

mag

).

(11)

With the accumulation of crystals suspended in the magmatic melt (for example, because of the less efficient settling of plagioclase), the denominator in Eq. (11) decreases and its numerator increases. As a consequence, the accumulation front is accelerated, which corresponds to the concave dotted line in Fig. 6. In calculating the modal composition of the "actual"
GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

Calculation Succession for the INTRUSION Subroutine The system of Eqs. (3)-(14) is the basis of the "intrusive unit" of the COMAGMAT program package [39, 40]. They are appended by relations that take into account the different crystallinity of the parental magmas and by a procedure for differentiating between the fluxes of in situ and homogeneous crystallization [42]. These equations for differentiation dynamics are successively solved by the INTRUSION routine, which is periodically "switched on" in the course of the construction of the main fractionation trajectory (Fig. 7). The modeling of processes in a chamber is started with the calculation of the equilibrium crystallization of a melt of specified composition corresponding to that of the parental magma. The simulations are carried out until a specified crystallinity int, which corresponds to the assumed fraction of intratelluric phenocrysts Fint at the moment of emplacement. This is the first stage
2006


84

ARISKIN, YAROSHEVSKY
START COMAGMAT PROGRAM

Initial parameters: Crystallization parameters: chemical composition of the cr = 0, cr = 0.02, max = 0.90 cr parental magma, pressure, accuracy of temperature calculations and redox conditions Parameters of the intrusive process: thickness of the body, characteristics of the host rocks, sinking velocity of mineral grains, crystallinity of the parental magma int, minimal porosity of cumulates, and the duration of chilling
SOLUTION OF THE THERMODYNAMIC EQUILIBRIUM PROBLEM AT A GIVEN MAGMA CRYSTALLINITY

cr int Yes cr = int No

No

DYNAMIC BLOCK OF THE PROGRAM (INTRUSION SUBROUTINE)

Calculation of chilled zone thicknesses Heat fluxes through the boundaries of the intrusion Overall fluxes of in situ crystallization Fluxes of homogeneous crystallization of minerals

cum FL F

cr L

Velocity of crystallization fronts, thickness and modal composition of in situ crystallization zones Yes Thickness and modal composition of cumulates

No

Modeled rock composition and correction of the residual magma composition cr = cr + cr cr = max cr
PRINT-OUT OF SIMULATION RESULTS

Information on the state of the crystallizing magma: temperature, modal proportions, compositions of minerals and melt Characteristics of the structure of the intrusive body: position of the crystallization and crystal accumulation fronts; proportions of minerals and intercumulus melt, chemical composition of modeled rocks
Fig. 7. Scheme of successive iterations during the simulation of the liquid line of descent of the parental magma (the INTRUSION subroutine is switched on periodically). This algorithm is used in the modeling of the structure of differentiated sills by the COMAGMAT computer program [39].

when information on the equilibrium (T, phase proportions, and the compositions of the crystals and melt) is input into the dynamic block in which the thermal characteristics of the host rocks and melt are specified and calculated, and the initial magma temperature T0 at the upper boundary [Eq. (3)] is assumed to be equal to the equilibrium temperature T. This enables the calculation of the heat flux from magma layers at the upper and lower contacts, along with the corresponding fluxes of the solid phase. The nucleation rate predetermines the thicknesses of the chilled zones (see [39] for details), whose origin results in an insignificant decrease in the magma volume (the thickness of the melt layer). This requires a correction of the residual system composi-

tion, after which information from the INTRUSION block is fed into the main module of the COMAGMAT program to continue the simulation of the crystallization process. A known solid phase flux at the upper contact enables the calculation of the time interval needed for the formation of a mass of crystals corresponding to crystallization increment t / I
UF

.

(15)

This makes it possible to calculate the thicknesses of the rock units that are produced at the in situ crystallization and accumulation fronts:
Vol. 44 No. 1 2006

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

85

h

UF(or LF)

= t,

(16)

in which the velocity of the front is calculated by Eqs. (8) or (9). Upon determining the modal and chemical composition of the derivatives [see Eqs. (10)-(14)], the changes in the system bulk composition and volume as a result of the origin of new rock units are calculated. With regard for this correction, information on the current state of the magma is loaded again in the block for the solution of the equilibrium problem (Fig. 7), in which the COMAGMAT program calculates the next crystallization step. These iterations with the INTRUSION subroutine are continued until the cumulus surface reaches the upper front or until the concentration of crystals in the cr magma exceeds the specified limit F LF . After this the calculations are terminated, and the modeled variations in the compositions of the rocks, cumulus phases, and melt become available for petrological-geochemical analysis. SIMULATION OF THE STRUCTURE OF DIFFERENTIATED SILLS The intrusive version of the COMAGMAT program can be used to reproduce the structures of naturally occurring tabular magmatic bodies, including the distribution of cumulus phases in their vertical sections, evaluation of the amounts of the captured or intercumulus melts, and the variations in the concentrations of ten major and twenty trace elements. This is a typical example of "forward" modeling, which is centered on the "response" of the system to the assumed initial conditions and the fitting parameters of the system. The results of these simulations are compared with natural data, the initial conditions and fitting parameters are then corrected, and the simulations are resumed until the observed and modeled characteristics become consistent (if possible). The initial conditions include the thickness of the intrusion, the composition of the parental magma, the thermal characteristics of the host rocks (the temperature and density of the melt are calculated in the process of crystallization), pressure, and oxygen fugacity f O2 . The fitting parameters determine the semiempirical character of the model, including the fraction of intratelluric phenocrysts, the coefficients of "heat losses" UF and proportionality in Eqs. (3) and (4), the duration of the chilling regime [42], the limiting cr density of the cumulus F LF , and the effective velocities of mineral settling j . If the structure of a natural intrusion is reproduced during the simulations, some of the fitting parameters can be interpreted as physical parameters, at least within the framework of the assumption that magma differentiation proceeds in compliance with the convection-accumulation model. The first experience in modeling the structure of intrusive traps was gained by the simulations of the
s

Kuz'movskii sill near the Podkamennaya Tunguska River and two intrusions in the upper reaches of the Vilyui River: V-304 and Vavukan [7]. These data were later recalculated with the use of a revised crystallization model [35] (the results are described in detail in [8, 43]). The intrusive version of the COMAGMAT model always yielded the distributions of major and trace elements highly consistent with those observed in the natural intrusions, with the distribution lines sometimes amazingly coinciding. The modal composition of the modeled rocks was in agreement with the petrological and structural characteristics of the rocks. This can be illustrated by the simulations results for the Vavukan sill [7, 8]. The sill is an almost tabular body, which crops out along the rivers of Vilyui and Vavukan, a right-hand tributary of the Vilyui. None of the natural exposures displays the whole thickness of the body, but two representative cliff outcrops spaced about 1.5 km apart allowed us to sample long enough vertical sections of the body (with a sampling step of 5 m). Petrographic data indicate that these exposed fragments contain the rocks of the same units, and this enabled us to reproduce the whole vertical section of the Vavukan intrusion ~100 m thick. Its upper and lower inner-contact zones consist of fine-grained dolerites. The "layered series" (which composes the bulk of the body) consists (from bottom to top) of units of poikilophitic, taxitophitic, and trachytoid gabbro-dolerites, which in places grade into lenses of ferrogabbro. The results of the geochemical research made it possible to reproduce the compositional layering of the intrusion, which is manifested in the systematic variations in the contents of major and trace elements in the vertical section (Fig. 8). The structure of the Vavukan intrusive was simulated for a pressure of 1 atm and redox conditions corresponding to the wuestite-magnetite (~10 rel. % Fe2O3) buffer. The parental magma was assumed to have a composition corresponding to the weighted mean composition evaluated in the reproduced integral vertical section of the intrusion and corresponding to a saturated tholeiite with ~50% SiO2 and 7.9% MgO. The liquidus temperature of this melt is approximately 1200Аы, and its simulated crystallization succession (Ol Ol + Pl Ol + Pl + Aug Ol + Pl + Aug + Mt) is typical of tholeiitic systems. The thermal characteristics of the rocks corresponded to the parameters of the models with a simplified phase diagram [4, 5]. The variable parameters included the crystallinity of the parental magma, the limiting amount of crystals in the cumulus, and the effective settling velocities of minerals. In examining the first model, we used arbitrarily chosen (but physically sensible) values, which were corrected in the course of the reconstruction of progressively more realistic distributions of major and trace elements. The successive refining of the fitting parameters (>50 variants) resulted in the determination of the optimum characteristics of the process, when the modeled
2006

GEOCHEMISTRY INTERNATIONAL

Vol. 44

No. 1


86 Vertical section 100 UBZ 80 Height, m 60 40 20 0 100 80 Height, m 60 40 20 0 0 0.5 1.0 1.5 0 0.5 Ni Ol + Pl cumulate LBZ Ol + Pl + Aug cumulate Ilm

ARISKIN, YAROSHEVSKY Normative composition Concentrations of major elements MgO Ol Px Pl FeO Al2O3 TiO2 K2O

Na2O

1 2 0 10 20 30 40 50 vol % 6 10 10 12 14 wt % 0 1 2 Concentrations of trace elements, ы/ы0 (ы0 are the original concentrations)

31

2

wt %

Cr

Mn

Sc

V

Sr

Ba

1 2

1.0

1.5

2.0 0.5

1.0 0.5

1.0 0.5

1.0

1.5 0.5

1.0

1.0

1.5

Fig. 8. Observed and modeled structural characteristics of the Vavukan trap intrusion. (1) Natural data; (2) data calculated by the INTRUSION subroutine: P = 1 atm, WM buffer, the initial composition and fitting parameters are listed in Tables 16 and 17 in [8]; the Ol-Pl and Ol-Pl-Aug cumulates correspond to the poikilophitic dolerites and to taxitophitic dolerites and gabbro-dolerites, respectively.

parameters of the body coincided in much detail with the observational results. For example, the distributions shown in Fig. 8 require the assumption of a low content of phenocrysts in the parental magma (Fint ~ 2%) and a relatively low degree of filling with cumulus crystals (close to 42-49 vol %). The optimum settling velocities of minerals varied thereby from 10 m/yr for plagioclase to 100 m/yr for clinopyroxene [8, 43]. A distinctive feature of the convection-accumulation model is the fact that each bend of the evolutionary lines and jumps in the concentrations of major and trace elements in Fig. 8 can be explained as a consequence of the evolution of the modal composition of the two- and three-mineral cumulate, including the fraction of the intercumulus liquid. This is illustrated by the data of Fig. 9, in which the distribution of trace elements is correlated with the parameters of the "simulated cumulates." As can be seen, the results of these computer simulations perfectly reproduce the S-shaped distribution profile for Ni, which is controlled by the distribution of cumulus olivine and, to a lesser extent, augite. An analogous explanation can be given for the distribu-

tions of Cr and V, whose behavior in the lower portion of the body was controlled by the amount of the intercumulus melt, and that in the middle and upper parts was determined by clinopyroxene and magnetite. Such correlations between the chemical and modal compositions of rocks can be traced for all elements. The smooth configuration of the evolutionary trajectories for much of the vertical section is predetermined by the condition of stirring, which is inherent to the convection-accumulation model. Another important outcome of this work on the simulation of differentiated traps (not only the Vavukan sill) was related to the interpreting of the structure of dolerites. The agreement between the petrochemical and geochemical characteristics of the natural and modeled bodies made it possible to reconstruct the phase composition of the initial melt-crystal mixtures that gave rise to the natural rocks and to correlate the phase boundaries with the layering observed in the sills [7, 8, 43]. Thus, the poikilophitic dolerites composing the bottoms of the sills were interpreted as cotectic rocks, which corresponded to the mutual crystallization
Vol. 44 No. 1 2006

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT 1.0 UBZ "SH" Cumulates GD Ni 1 2

87

0.8 0.6

TOD 0.4 Lower zone POD 0.2

Chill zone

0 1.0 UBZ "SH" Cumulates GD

Cumulus Ol 0 100 200 300 400 0 Cr 10 20 30

0.8 0.6 Relative height

TOD 0.4 Lower zone POD

0.2

Chill zone

0 1.0 UBZ "SH" Cumulates GD

Cumulus Aug 0 200 400 600 800 0 V 10 20 30

0.8 0.6

Intercumulus liquid

TOD 0.4 Lower zone POD

0.2

0

Chilled zone

Melt amount

Structure of 100 200 300 400 500 40 60 80 100 intrusion Concentration in rock, ppm Modal proportions, vol %

Fig. 9. Observed and modeled relations between major rock types, concentrations of trace elements, and the original proportions of crystals and melt in the rocks of the Vavukan trap intrusion [8, 43]. (1) Natural data; (2) simulations. The lower zone includes chilled rocks and poikilophitic dolerites (POD), which crystallized with the magmatic mixture captured at the solidification front. The "cumulate" zone with 50-60% porosity of the sediment consists of taxiophitic dolerites (TOD) and gabbro-dolerites that grade into "ferrogabbro" of the sandwich horizon (SH). The upper contact zone was formed by means of the directed growth from the chamber roof. The smooth concentration trends correspond to the high efficiency of the vigorous concentration convection. GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1 2006


88

ARISKIN, YAROSHEVSKY

of olivine and plagioclase (Fig. 8). The taxitophitic dolerites mark the appearance of the first large clinopyroxene oikocrysts, and the prismatic-granular gabbro-dolerites are "mature" Ol-Pl-Aug cumulate, which can show evidence of the presence of the first grains of cumulus magnetite [43, 44]. The problematic aspect in the interpretation of these data is related to the uniqueness of the solution of the modeling problem. It can be thought that the variable parameters are too numerous, and thus, similar distributions of major and trace elements could have been produced by different combinations of the initial conditions and/or fitting coefficients of the model. This situation calls for a proof of the stability of the solution depending on the initial data. Methods for the analysis of such systems involve the introduction of random fluctuations in the parameters of the model and the inspection of the reproducibility of the simulation results [45]. This approach practically cannot be applied when the INTRUSION subroutine is utilized, because this would require not only multiple changes in the ten variables but also the examination of the simultaneous "response" of the model in terms of the concentrations of many chemical elements. Nevertheless, some evidence testifies to the stability of the solution of the modeling problem of intrusive bodies. First, in sorting out and searching for optimal fitting parameters, one should bear in mind that the geochemical data with which the model distributions are compared include the concentrations of about 30 elements in a few dozen samples. Hence, the system of the analytical data is strongly "overestimated" in terms of the number of the variable model parameters. Our experience in reproducing the structures of sills with the use of the convection-accumulation model indicates that, regardless of the predicted values, realistic schemes of the distributions of elements are consistent with the same stable combination of the fitting parameters, which vary within ~10 rel. %. Moreover, there are means for an independent test of the results of direct simulations: this is the method of geochemical thermometry, which enables the evaluation of the parameters of differentiation within a chamber by reproducing crystallization conditions for each individual rock sample. The methodological aspects of this technique were described in [46] and reviews of the results were published in [8, 39, 47]. REPRODUCING THE CRYSTALLIZATION CONDITIONS OF DOLERITES BY THE METHOD OF GEOCHEMICAL THERMOMETRY: 1987-2004 An important avenue of the development of the convection-accumulation model involves the elaboration of methods for the solution of the inverse problem in relation to the evaluation of the crystallization conditions of rocks on the basis of the spatial distribution of their petro- and geochemical characteristics. This research was carried out simultaneously with the forward modeling of intrusive bodies and was focused on

the refining of the parameters of the natural processes of the separation of crystals and melt in a chamber. The first "inverse models" were underlain by the scheme of ideal fractionation, which implied that part of the residual melt remained captured in the interstitial space of cumulates. The utilization of these models made it possible to determine the porosity of cumulates and facilitated the designing of a method for the estimation of the mineral-melt distribution coefficients during the origin of subtabular bodies [48]. A principally new stage began when these studies were conducted with the use of the thermodynamic simulations of crystallizing basalts [33-35]. In 1987-1988, two papers [46, 49] were published that justified a systematic approach to the reconstruction of the directed fractionation of an intrusive magma, including the variations in the temperature and composition of the residual liquid, the initial compositions of the cumulus phases, and the relative proportions of the crystals and melt. The essence of this method (which was named geochemical thermometry [46]) follows from the concepts of the solidification dynamics and distinctive genetic features of intrusive rocks that are illustrated in Figs. 6 and 10. The main postulate of this method is that the rock succession produced during magma differentiation in a chamber is formed from below upward, in such a manner that the volume of the residual magma lies above the solidification front or the cumulus surface (Fig. 6). An important consequence of this is as follows: by averaging the composition of the upper part of the intrusion, it is possible to evaluate, for each rock, the composition of the melt-crystal mixture (residual magma) in which crystals settled and cumulus was formed at a given level. This assumption corresponds to the known approach of Wager and Brown [13] proposed to calculate the composition of the residual magmas of the Skaergaard intrusion. A difference is that we regard magma not as a homogeneous liquid but as a suspension of crystals in a melt (Fig. 10a), whose proportions are controlled by the temperature of the mixture and the compositions of which are interrelated through the thermodynamic equilibrium conditions. Another difference concerns the composition of the cumulus that is formed by the settling of crystals from this suspension. It is natural to assume that if the melt is stirred and homogenized efficiently enough, the temperature of the melt captured in the cumulus corresponds to that of the main magma volume. In this process the compositions of cumulus minerals should not differ from the average compositions of the crystalline phases suspended in the magma. Hence, the bulk compositions of the magma and cumulus during each differentiation stage in the chamber can be expressed as a combination of the same mineral phases and melt [see Eqs. (13) and (14)] that are in equilibrium at the same temperature (Fig. 10b). This situation is promising for the application of the computer simulation of crystallization with the aim of evaluating the thermodynamic characteristics of the residual magma and the corresponding rock [8, 39].
Vol. 44 No. 1 2006

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT (a) (b) Temperature (c)

89

Cumulus Magma Ol

Pl

Rocks

Ol

Cumulus

+ Pl

Residual melt

Fig. 10. Principle of the geochemical thermometry of cumulates and the corresponding residual magmas in the hypothetical vertical section of an intrusive body. (a) Scheme of the inner structure of a chamber that implies the possibility of estimating the composition of the magma for the developing cumulus (see text); (b) phase composition of the residual magma and the equilibrium crystalline precipitate; (c) expected intersection of the liquid lines of descent for the compositions of the magma and cumulus.

The technique of geochemical thermometry involves the numerical simulation of liquid lines of descent for pairs of initial compositions, one of which corresponds to the bulk rock composition (at a given level) and the other to the weighted mean composition of the overlying part of the intrusion. The intersection point of these lines in composition-temperature space determines the crystallization conditions of the rock (proportions of cumulus phases and their compositions), the phase composition of the magma, and the chemistry of the residual melt. Figure 10b schematically shows such an intersection for a single oxide: MgO. In a multicomponent compositional space, analogous intersections exist for all oxides, and, hence, it is usually easy to visually determine the initial equilibrium temperatures (methodological aspects of the geochemical thermometry are described in detail in [47]). If data on the compositional succession of the rocks and residual magma are available, the paired geochemical thermometry of these compositions enables the researcher to reconstruct the thermodynamic characteristics of magma differentiations and the phase proportions of the cumulates produced during the solidification of the intrusive body. THERMOMETRY OF DOLERITES IN DIFFERENTIATED SILLS The approach described above was tuned up using the example of dolerites from the Vel'minskii sill,
GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

which is exposed in a series of outcrops in the valley of the Podkamennaya Tunguska River in eastern Siberia [46, 49]. The simulations were conducted with 27 selected rock compositions, with the rock samples taken in the vertical section of the body at intervals of 3-5 m. These data were used to calculate 26 compositions of the residual magmas, which corresponded to the successive compositional evolution of the taxitophitic dolerites and gabbro-dolerites [8]. The calculations were conducted at a pressure of 1 atm at an oxygen fugacity corresponding to the WM buffer, with these conditions ensuring the crystallization of magnetite and the origin of ferrogabbro during the closing differentiation stages. Geochemical thermometry yielded the temperature range of ~1160-1135Аы for the crystallization of the magma [49] and systematic variations in the composition of the residual melt. The estimations of the original modal proportions of the cumulus suggested that it contained the Ol + Pl + Aug cotectic assemblage early in the course of differentiation. Our further thermometric calculations were conducted for trap intrusions (Kuz'movskii, Vavukan, V-304) that were used previously to verify the convection-accumulation model (see above). A conspicuous feature common to these three bodies, but not the Vel'minskii sill, is their elevated MgO contents (8-9 wt %). Our research indicates that this was not a result of the original magma richness in olivine or augite crystals but suggests the more primitive composition of the parental
2006

Ol
Pl

Magma

Residual melt

Melt composition (MgO)


90 Relative height 1.0 Layered series 0.8 0.6 0.4 0.2 0 10 20 30 40 50

ARISKIN, YAROSHEVSKY

Magma Ol Pl Aug Total

for phase equilibria. Thus, the analysis of the data portrayed in Fig. 11 leads to a principal conclusion: the uniqueness of the solution of the modeling problem is predetermined by the existence of a single distribution of cumulus phases in the vertical section that describes the totality of data on the distributions of major and trace elements. CONCLUSIONS The principal outcome of our research is the development of a physicochemical model that interrelates methods for the simulation of realistic liquid lines of descent for basaltic magmas (COMAGMAT computer program) and an algorithm for the modeling of the dynamics of magma differentiation in compliance with the convection-accumulation mechanism (INTRUSION subroutine). Although the very first versions of this multiparametric model were proposed more than two decades ago, it still remains perhaps the only example of the systematic solution of the equilibrium problem and heat transfer equations, with the model simulations yielding the spatial structure of a magmatic body. The efficiency of the algorithms was proved by the simulations of the structures of sills in the Siberian Platform. The possibility of reproducing the distributions of major and trace elements and the variations in the mineralogy of the rocks leave little doubt that concentration convection in the form of suspension flows is the leading sedimentation mechanism of crystallizing material and plays a decisive role in the origin of the geochemical structure of these tabular intrusions. An important element of the work on the modeling of processes in magmatic chambers is calculations by the method of geochemical thermometry. These calculations enable the researcher to obtain complete information on the conditions under which the residual magmas evolved and the crystallization conditions of the corresponding rocks. The key parameter is the temperature of the residual magma (intercumulus melt), which is related to other thermodynamic parameters, including the modal and chemical composition of the magmatic suspensions and the proportions and chemistry of cumulus minerals. The results obtained by the thermometry of dolerites from differentiated sills have confirmed the results of the forward simulations of their structures. This leaves no doubt concerning the unambiguity of the differentiation parameters, which were quantified with the convection-accumulation model. Geochemical modeling can be regarded as the most consistent method for the genetic interpretation of intrusive rocks; it combines a variety of approaches based on approximations of the equilibrium distribution of components in melt-crystal systems. The significance of our research is determined not only by its applicability to the reconstructions of the structures of naturally occurring magmatic bodies. The simulations of differentiation dynamics has a fundamental methodological aspect related to the evaluation
Vol. 44 No. 1 2006

60 0 10 20 30 Cumulus phases, vol %

Fig. 11. Variations in the contents of cumulus minerals and crystals suspended in magma for the optimal model of the Vavukan intrusion structure displayed in Figs. 8 and 9. The results of simulations by means of forward modeling by the INTRUSION subroutine are shown as solid and dashed lines. Symbols are the results of geochemical thermometry for the rock succession and corresponding residual magmas [8, 39].

magmatic liquid. This is also confirmed by the occurrence of poikilophitic dolerites, which are the highest temperature rocks and consist of the Ol + Pl cotectic assemblage. According to thermometric data, the cumulus from which the dolerites crystallized was characterized by the systematic accumulation of solid phases, a process that reached a maximum (45-55%) in the central parts of the sections. The appearance of primary cumulus clinopyroxene coincided with petrographic evidence of the transition from poikilophitic (Ol + Pl) to taxiophitic (Ol + Pl + Aug) dolerites. This is clearly illustrated by the left-hand plot in Fig. 11, in which symbols indicate the proportions of cumulus minerals evaluated for the Vavukan intrusion. The proportions of melt and crystals in the right-hand plot characterize the evolution of the phase composition of the residual magma during differentiation. The evaluations of the amounts of olivine, plagioclase, and Cpx in the magma and cumulus may be interesting for researchers dealing with flood-basalt magmatism in eastern Siberia. In the context of this publication, we would like to emphasize another feature of the plots in Fig. 11, in which the data of geochemical thermometry are displayed together with the distributions obtained by the simulation of differentiation in an "intrusive regime" of the COMAGMAT model [39]. It can be seen that the solutions of the forward and inverse problems are practically identical. This is a convincing argument in favor of the plausibility of the convection- accumulation mechanism and the consistency of the thermodynamic and dynamic parameters of processes in the chambers. Note that the methods of the forward and inverse modeling are independent, and their only common feature is the utilization of the same models

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT

91

of the petrological significance of small- and largescale processes of mass transfer in a magmatic chamber. The main conclusion is that the diffusion and kinetic aspects (which are responsible for the origin of the texture of rocks, the sizes of their grains, and the zoning of mineral crystals) play a subordinate role in the origin of compositional layering compared with the role of processes responsible for the separations of the melt and solid phases [2-5]. This separation is controlled by the homogeneous nucleation of the crystalline phase in the boundary layers and occurs as various convection regimes, including magmatic suspension flows and the Stokes sinking of crystals in the liquid [8]. Combinations of these major mechanisms of mass transfer in a magmatic chamber give rise to the wealth of layering types in intrusive bodies, including the origin of cyclically layered series and units [22]. Similar conclusions about interactions between convection processes related to phase separations were published in [10, 31]. This viewpoint on the physical nature of differentiation in magmatic cambers highlights avenues for the further development of the convection-accumulation model. One of these avenues is related to the movement dynamics of interstitial melt in cumulus. The problem is that the method proposed for the simulation of the structures of intrusions takes into account only the processes of mass transfer that occur in the residual magma and at its boundaries. It is thereby postulated that the cooling of the crystal-melt mush (Fig. 6, interval t4-t6) does not contribute significantly to the extensive differentiation and fractionation of the magma. This approach proved efficient for the modeling of relatively thin tabular bodies (Fig. 9) but can hardly be warranted for large layered intrusions, whose rocks show evidence of the compaction and recrystallization of the cumulates [50]. A problem awaiting its solution is the structure of the lower contact zones, which show an "inverse evolutionary trend" of the composition of the cumulus and rock-forming minerals, as is manifested in an upward transition from low- to high-temperature mineral assemblages [51]. Searches for the solutions of these problems within the scope of diffusion processes (for example, by appealing to the Soret effect [52]) have no solid ground because they do not allow either for the heterogeneity and the subeutectic character of magmatic melts or for relations between the cooling rate and the scale of the possible diffusion transfer of components [see Eq. (1)]. A constructive approach should allow for the effects of crystal settling [53] or compositional convection, which is related to the migration and reaction activity of pore melt [54, 55]. These phenomena result in the reequilibration (recrystallization) of cumulus minerals and, inasmuch as the melt is forced out, can significantly contribute to the compositional evolution of the bulk of the residual magma [56, 57]. The first attempts to simulate infiltration processes in cumulus [58] and
GEOCHEMISTRY INTERNATIONAL Vol. 44 No. 1

the effects of fractionation in the boundary layer [59] represent significant steps in this direction. ACKNOWLEDGMENTS The authors thank E.V. Koptev-Dvornikov and G.S. Barmina for valuable comments concerning the content of the manuscript. G.S. Nikolaev and K.A. Bychkov are thanked for help during the preparation of the graphical materials. This study was supported by the Russian Foundation for Basic Research, project nos. 05-05-64906 and 03-05-64569. REFERENCES
1. M. Ya. Frenkel, "The Thermodynamics, Dynamics, and Mathematical Modeling of Geochemical Systems," Geokhimiya, No. 10, 1401-1411 (1992). 2. M. Ya. Frenkel and A. A. Yaroshevsky, "The Crystallization Differentiation of Intrusive Magma: The Diffusion Mechanism for Heat and Mass Transfer," Geokhimiya, No. 8, 1197-1203 (1976). 3. M. Ya. Frenkel and A. A. Yaroshevsky, "The Crystallization Differentiation of Intrusive Magma: Convection and Freezing-on Conditions," Geokhimiya, No. 11, 1624- 1632 (1976). 4. M. Ya. Frenkel and A. A. Yaroshevsky, "A Set of Equations for Heat and Mass Transfer during the Emplacement of a Tabular Intrusive Body: Problem Formulation and a Computer Algorithm for Its Solution," Geokhimiya, No. 4, 547-559 (1978). 5. M. Ya. Frenkel and A. A. Yaroshevsky, "The Crystallization Differentiation of Intrusive Magma: Mathematical Modeling of the Thermal Conditions and Differentiation of an Tabular Intrusive Body with Regard for Solid Phase Sedimentation," Geokhimiya, No. 5, 643-668 (1978). 6. E. V. Koptev-Dvornikov, A. A. Yaroshevsky, and M. Ya. Frenkel, "Crystallization Differentiation of Intrusive Magma: The Advantages and Disadvantages of the Sedimentation Model," Geokhimiya, No. 4, 488-508 (1979). 7. M. Ya. Frenkel, A. A. Yaroshevsky, E. V. Koptev-Dvornikov, et al., "The Crystallization Mechanism in the Genesis of Layered Intrusions," Zap. Vsess. Mineral. O-va CXIV (3), 257-274 (1985). 8. M. Ya. Frenkel, A. A. Yaroshevsky, A. A. Ariskin, et al., Dynamics of Basic Magma Differentiation in Chambers (Nauka, Moscow, 1988) [in Russian]. 9. V. N. Sharapov and A. N. Cherepanov, The Dynamics of Magma Differentiation (Nauka, Novosibirsk, 1986) [in Russian]. 10. C. Jaupart and S. Tait, "Dynamics of Differentiation in Magma Reservoirs," J. Geophys. Res. 100, 17 615- 17 636 (1995). 11. N. L. Bowen, The Evolution of the Igneous Rocks (Princeton Univ. Press, Princeton, 1928). 12. H. H. Hess, Stillwater Igneous Complex, Montana: A Quantitative Mineralogical Study (1960). 13. L. R. Wager and G. M. Brown, Layered Igneous Rocks (Oliver, Edinburgh, 1967; Mir, Moscow, 1970).
2006


92

ARISKIN, YAROSHEVSKY rithm of Computer-based Solution," Geokhimiya, No. 5, 679-690 (1984). M. Ya. Frenkel and A. A. Ariskin, "Computer-based Modeling of Equilibrium and Fractional Crystallization under Determined Oxygen Fugacity," Geokhimiya, No. 10, 1419-1431 (1984). A. A. Ariskin, G. S. Barmina, and M. Ya. Frenkel, "Crystallization of Basaltic Melts under Determined Oxygen Fugacity: Computer-based Modeling," Geokhimiya, No. 11, 1614-1628 (1986). P. L. Roeder and E. Emslie, "Olivine-Liquid Equilibrium," Contrib. Mineral. Petrol. 29, 275-289 (1970). M. J. Drake, "Plagioclase-Melt Equilibria," Geochim. Cosmochim. Acta 40, 457-465 (1976). R. L. Nielsen and M. J. Drake, "Pyroxene-Melt Equilibria," Geochim. Cosmochim. Acta 43, 1259-1272 (1979). A. A. Ariskin and G. S. Barmina, Modeling of Phase Equilibria for the Crystallization of Basaltic Magmas (Nauka, Moscow, 2000) [in Russian]. A. A. Ariskin, G. S. Barmina, M. Ya. Frenkel, and R. L. Nielsen, "COMAGMAT: A Fortran Program to Model Magma Differentiation Processes," Comput. Geosci. 19, 1155-1170 (1993). A. A. Ariskin, "Phase Equilibria Modeling in Igneous Petrology: Use of COMAGMAT Model for Simulating Fractionation of Ferro-basaltic Magmas and the Genesis of High Alumina Basalt," J. Volcanol. Geotherm. Res. 90, 115-162 (1999). M. Ya. Frenkel, "The Geochemical Structure of a Tabular Intrusion," in Dynamic Models in Physical Geochemistry (Nauka, Novosibirsk, 1982), pp. 19-30 [in Russian]. M. Ya. Frenkel, A. A. Yaroshevsky, A. A. Arislin, et al., "Convective-Cumulative Model Simulating the Formation Process of Stratified Intrusions," in Magma-Crust Interactions and Evolution, Ed. by B. Bonin (Theophrastus Publications, Athens, 1989), pp. 3-88. E. V. Koptev-Dvornikov, E. E. Kameneva, B. S. Kireev, and A. A. Yaroshevsky, "Evidence for the Cumulus Origin of Clinopyroxene and for the Reequilibration of Olivine in Vavukan-Sill Dolerites," Geochem. Int. 33 (1), 81-102 (1996). A. I. Shapkin, Doctoral Dissertation in Chemistry (GEOKhI, Moscow, 1998). M. Ya. Frenkel, A. A. Ariskin, G. S. Barmina, et al., "Geochemical Thermometry of Igneous Rocks: Principles and Examples," Geokhimiya, No. 11, 1546-1562 (1987). A. A. Ariskin and G. S. Barmina, "COMAGMAT: Development of a Magma Crystallization Model and Its Petrologic Applications," Geochem. Int. 42 (Suppl. 1), S1- S157 (2004). G. S. Barmina, M. Ya. Frenkel, A. A. Yaroshevsky, and A. A. Ariskin, "Crystallization-related Redistribution of Trace Elements in Sheet Intrusions," in Dynamic Models in Physical Geochemistry (Nauka, Novosibirsk, 1982), pp. 45-55 [in Russian]. G. S. Barmina, A. A. Ariskin, E. V. Koptev-Dvornikov, and M. Ya. Frenkel, "Experience in Studying the Composition of Primary Cumulus Minerals in Differentiated Traps," Geokhimiya, No. 8, 1108-1119 (1988).
Vol. 44 No. 1 2006

14. E. D. Jackson, "Primary Textures and Mineral Associations in the Ultramafic Zone of the Stillwater Complex, Montana," US Geol. Surv. Prof. Pap. 358, 1-106 (1961). 15. I. H. Campbell, "Some Problems with the Cumulus Theory," Lithos 11, 311-323 (1978). 16. A. R. McBirney and R. M. Noyes, "Crystallization and Layering of the Skaergaard Intrusion," J. Petrol. 20, 487-554 (1979). 17. R. W. Bartlett, "Magma Convection, Temperature Distribution, and Differentiation," Am. J. Sci. 267, 1067-1082 (1969). 18. A. P. Vinogradov and A. A. Yaroshevsky, "Physical Conditions of Zoned Melting in the Earth," Geokhimiya, No. 7, 779-790 (1965). 19. R. G. Cawthorn and T. S. McCarthy, "Incompatible Trace Element Behavior in the Bushveld Complex," Econ. Geol. 80, 1016-1026 (1985). 20. E. V. Sharkov, The Petrology of Layered Intrusions (Nauka, Leningrad, 1980) [in Russian]. 21. N. I. Gel'perin and G. A. Nosov, The Basis of the Mechanism for Melt Crystallization (Khimiya, Moscow, 1975) [in Russian]. 22. M. Ya. Frenkel, The Heat and Chemical Dynamics of Basic Magma Crystallization (Nauka, Moscow, 1995) [in Russian]. 23. A. N. Tikhonov and A. A. Samarskii, Equations of Mathematical Physics (Nauka, Moscow, 1966) [in Russian]. 24. M. G. Worster, H. E. Huppert, and R. S. J. Sparks, "Convection and Crystallization in Magma Cooled from Above," Earth Planet. Sci. Lett. 101, 78-89 (1990). 25. M. T. Mangan and B. D. Marsh, "Solidification Front Fractionation in Phenocryst-free Sheet-like Magma Bodies," J. Geol. 100, 605-620 (1992). 26. A. Simakin, V. Trubitsyn, and H. Schmeling, "Structure of the Boundary Layer of a Solidifying Intrusion with Crystal Sedimentation," Earth Planet. Sci. Lett. 126, 333-349 (1994). 27. B. D. Marsh, "Solidification Fronts and Magmatic Evolution," Mineral. Mag. 60, 5-40 (1995). 28. V. N. Sharapov, A. N. Cherepanov, V. N. Popov, and A. G. Lobov, "Dynamics of Basic Melt Cooling during the Filling of a Funnel-shaped Intrusive Chamber," Petrologiya 5 (1), 10-22 (1997) [Petrology 5 (1), 8-19 (1997)]. 29. R. T. Helz, "Crystallization History of Kilauea Iki Lava Lake as Seen in Drill Core Recovered in 1967-1979," Bull. Volcanol. 43, 675-701 (1980). 30. R. T. Helz, "Differentiation Behavior of Kilauea Lava Lake, Kilauea Volcano, Hawaii: An Overview of Past and Current Work," in Magmatic Processes: Physicochemical Principles, Ed. by B. O. Mysen, Geochem. Soc. Spec. Publ. 1, 241-258 (1987). 31. S. Tait and C. Jaupart, "The Production of Chemically Stratified and Adcumulate Plutonic Igneous Rocks," Mineral. Mag. 60, 99-114 (1996). 32. V. P. Trubitsyn and E. V. Kharybin, "Convection in Magma Chambers due to Inversion of Vertical Distribution of Deposited Crystals," Fiz. Zemli, No. 5, 47-52 (1997). 33. M. Ya. Frenkel and A. A. Ariskin, "The Problem of Equilibrium Crystallization for the Basaltic Melt: An Algo-

34.

35.

36. 37. 38. 39. 40.

41.

42. 43.

44.

45. 46.

47.

48.

49.

GEOCHEMISTRY INTERNATIONAL


CRYSTALLIZATION DIFFERENTIATION OF INTRUSIVE MAGMATIC MELT 50. R. H. Hunter, "Texture Development in Cumulate Rocks," in Layered Intrusions (Elsevier, New York, 1996), pp. 77-101. 51. R. M. Latypov, "The Origin of in Basic-Ultrabasic Sills with S-, D-, and I-shaped Compositional Profiles by insitu Crystallization of a Single Input of Phenocryst-poor Parental Magma," J. Petrol. 44, 1619-1656 (2003). 52. R. M. Latypov, "The Origin of Marginal Compositional Reversals in Basic-Ultrabasic Sills and Layered Intrusions by Soret Fractionation," J. Petrol. 44, 1579-1618 (2003). 53. A. A. Yaroshevsky, "Physicochemical Principles in the Behavior of a Magmatic System in the Gravitational Field with a Low Fraction of Melt: Melt Segregation and Cumulus Development," Geokhimiya, No. 6, 670-675 (2003) [Geochem. Int. 41 (6), 603-608 (2003)]. 54. T. N. Irvine, "Magmatic Infiltration Metasomatism, Double-diffusive Fractional Crystallization, and Adcumulus Growth in the Muskox Intrusion and Other Lay-

93

55. 56. 57.

58.

59.

ered Intrusions," in Physics of Magmatic Processes (Princeton Univ. Press, Princeton, 1980), pp. 325-384. S. Tait and C. Jaupart, "Compositional Convection in a Reactive Crystalline Mush and Melt Differentiation," J. Geophys. Res. 97, 6735-6756 (1992). C. H. Langmuir, "Geochemical Consequences of In Situ Crystallization," Nature 340, 199-205 (1989). R. L. Nielsen and S. E. DeLong, "A Numerical Approach to Boundary Layer Fractionation: Application to Differentiation in Natural Magma Systems," Contrib. Mineral. Petrol. 110, 355-369 (1992). A. E. Boudreau and W. P. Meurer, "Chromatographic Separation of the Platinum-Group Elements, Gold, Base Metals, and Sulfur during Degassing of a Compacting and Solidifying Igneous Crystal Pile," Contrib. Mineral. Petrol. 134, 174-185 (1999). T. Kuritani, "Thermal and Compositional Evolution of a Cooling Magma Chamber by Boundary Layer Fractionation: Model and Its Application to Primary Magma Estimation," Geophys. Rev. Lett. 26, 2029-2032 (1999).

GEOCHEMISTRY INTERNATIONAL

Vol. 44

No. 1

2006