## Abstract

The observation that basic igneous rocks most commonly are holocrystalline under a wide spectrum of cooling regimes implies that cooling and crystallization can be uncoupled and considered separately. This is tantamount to realizing that the Avrami number is large in most igneous systems. Crystallization automatically adjusts through nucleation and growth to the cooling regime, and all aspects of the ensuing crystal population reflect the relative roles of nucleation and growth, which reflect the cooling regime. The characteristic scales of crystal size, crystal number and crystallization time are intimately tied to the characteristic rates of nucleation and growth, but it is the crystal size distributions (CSDs) that provide fundamental insight on the time variations of nucleation and growth and also on the dynamics of magmatic systems. Crystal size distributions for batch systems are calculated by employing the Johnson–Mehl–Avrami equation for crystallinity related to exponential variations in time of both nucleation and growth. The scope of the CSD is set by the difference a – b, where a and b are exponential constants describing, respectively, nucleation and growth. The batch CSD has constant slope and systematically migrates to larger crystal size (L) with increasing crystallinity. The diminution in nucleation with loss of melt is reflected in the CSD at late times by a strong decrease in population density at small crystal sizes, which is rarely seen in igneous rocks themselves. Observed CSDs suggest that a – b ∼6–10 and that b ∼0). That is, growth rate is approximately constant and nucleation rate apparently increases exponentially with time. Correlations among CSD slope, intercept, and maximum crystal size for both batch and open systems suggest that certain diagnostic relations may be useful in interpreting the CSD of comagmatic sequences. These systematics are explored heuristically and through the detailed examination of comagmatic CSDs in a number of igneous and industrial systems including, amongst other, Makaopuhi lava lake, Atka volcanic center, Peneplain sill, Dome Mountain lavas, Shonkin Sag lacolith, and Kilauea Iki lava lake. None of these systems shows CSDs typical of purely batch or purely open systems, even when the system itself is known on independent grounds to be a batch system. Instead, the CSDs of each system reflect a combination of kinetic and dynamic influences on crystallization. Heterogeneous nucleation and annexation of small crystals by larger ones, entrainment of earlier grown and ripened crystals, rate of solidification front advance, and protracted transit of a well-established mush column are some of the effects revealed in the observed CSDs. There may be an overall CSD evolution, reflecting the maturity of the magmatic system, from simple straight nonkinked CSDs in monogenetic systems to multiply kinked, piecewise continuous CSDs on well-established systems such as Hawaii and Mount Etna. This is not unlike the evolution of CSDs in some industrial systems. Finally, the fact that comagmatic CSDs are not often captured evolving systematically through large changes in nucleation rates, even in low crystallinity systems, may suggest that magma is always laced with high population densities of nuclei, supernuclei, and crystallites or clusters that together set the initial CSD at high characteristic population densities. Further evolution of the CSD occurs through sustained heterogeneous nucleation and rapid annealing at all crystallinities beginning at the liquidus itself and operating under more or less steady (not exponentially increasing) rates of nucleation.

## Introduction

A long-standing obstacle to understanding igneous rocks rests in relating experiment, theory, and process to the textures and crystal sizes in comagmatic sequences. That crystals can be readily grown, even in the simplest of experiments, mimicking in size and texture those of rocks, has been known since the experiments of James Hall (1798; Dawson, 1992). Detailed kinetic models are available that go a long way towards producing reasonable crystal populations and crystal sizes, albeit often under somewhat restrictive conditions (e.g. Dowty, 1980; Kirkpatrick, 1983; Brandeis et al., 1984; Hort & Spohn, 1991a). The method of crystal size distributions (CSDs) is also available as a means of quantitatively describing crystal content as a function of size and of directly modeling crystallization in magmatic processes (Marsh, 1988). Real time observations of crystallization in systems analogous to silicates are also available (Means & Park, 1994). There is, however, a fundamental need and opportunity to intricately link these separate areas of study using the rocks themselves.

The present work analytically explores the links between simple kinetic models of crystallization and the observed and predicted CSDs in batch and open systems. These results are then used to interpret the crystallization history of some typical magmatic sequences. The overall intent is to expose the intimate interplay between rates of nucleation and growth in establishing the common crystal sizes of basaltic rocks and to evaluate the role of magmatic processes on rock crystallinity.

This work builds on a single, simple axiom, hereafter called the Crystallization Axiom. Given a time scale imposed by cooling, nucleation and growth automatically adjust locally to attain full crystallinity. That is, the great preponderance of igneous rocks crystallize essentially completely regardless of the cooling regime. This is especially so for intrusive rocks. Although crystal numbers and sizes may vary enormously across sills, dikes, and plutons, by and large, the rocks themselves are holocrystalline. This is also mainly true for lavas. The clear exceptions are the glassy rocks (e.g. Carmichael, 1979), and an understanding of this class of rocks is also enlarged once the holocrystalline rocks are understood. A corollary of the Crystallization Axiom is that cooling and crystallization can, to a good first approximation, be decoupled.

This decoupling is equivalent to assuming, in the nondimensional form of the conservation of energy equation, that the Avrami number, which is a measure of the relative importance of the thermal (e.g. conductive cooling) time scale to the kinetic time scale (each to the fourth power), is large (Hort & Spohn, 1991a). For large values of the Avrami number (>∼1010) crystallization is controlled by the rate of removal of heat from the system. That is, the time scale (and hence heat production) associated with crystal nucleation and growth is locally short relative to that for cooling of the body as a whole. For bodies larger than ∼10 m in characteristic size (e.g. half-thickness) and commonly observed growth and nucleation rates, the Avrami number is always large (»1010). That crystallization can generally keep up with cooling allows crystallization to be analyzed independent of choosing a specific model for heat transfer, which is the course taken here. For the other extreme of Avrami number (i.e. «1010), kinetic effects dominate cooling and the magma vitrifies. That is, when the time scale for cooling is short relative to the kinetic time scale, crystallization is exceedingly ineffective. This overall behavior is also well known from experimental studies of crystallization as a function of cooling rate (e.g. Gibb, 1974; Lofgren, 1980).

I begin below by introducing the basic scales associated with crystal number, crystal size, and crystallization time. This introduces three characteristic constants that can only be found by employing specific kinetic models of crystallization. The Avrami approach, by way of the Johnson–Mehl–Avrami (JMA) equation, is employed to find these constants under a number of different assumptions of the variation of nucleation and growth as a function of time. Next the CSD equation for a batch system is introduced and discussed in relation to Avrami crystallization. Insights from the CSDs of common basaltic rocks are used to constrain kinetic models of crystallization. There are several major conclusions of this analysis: (1) the characteristic time of crystallization is only weakly dependent on the exact kinetics of crystallization, but typical crystal number and size are more sensitive to kinetics; (2) nucleation and growth may each be exponentially related to time, but the exponent of the time variation of nucleation is almost always some 6–8 times larger than that for growth; (3) crystal size is mainly the result of heterogeneous nucleation and rapid, continual grain annexation and boundary migration, as in Ostwald ripening, during all stages of crystallization beginning at the liquidus itself. The CSDs in a number of igneous systems are discussed relative to the insights gained from the analytical results.

## The Characteristic Scales of Crystallization: Number, Size, and Time

Although it is undeniable that crystallization results from an intricate interplay of kinetics, time, and temperature, the last of which couples phase equilibria to magma size, shape, and dynamics, it is crucial to focus on crystallization as a function of time and not temperature. Models explicitly involving temperature require, perforce, deciding at the outset on an appropriate cooling regime, the details of which often dominate and heavily color the results. In effect, crystallization must be viewed through the cooling regime. Instead, I focus here explicitly on a simple model of batch crystallization as a function of time so that the principles and results will be internally consistent, general, and not be clouded by the considerably intricate and often debatable functional relations between time, temperature, and phase equilibria. The advantage of this approach will shortly become clear. I emphasize that this is a model against which to compare rocks. In essence, the following furnishes a standard state which is, by definition, unattainable.

Motivated by the Crystallization Axiom, I assume that crystallization of an entire unit volume of magma can be described by characteristic rates of nucleation (Jo) and growth (Go) associated with a characteristic time (tc) of crystallization. This does not mean that these rates are taken to be constant during crystallization, but simply that these values are associated with more general functions describing the rates of nucleation and growth. The functions themselves will be considered shortly. Moreover, I ignore the fact that each solid phase may have specific kinetic characteristics, and treat the crystallization process as a bulk process. The differences in crystal number and size between, say, ilmenite and plagioclase in tholeiitic basalt are taken to be second-order effects. Hence, these results pertain to a characteristic population of crystals, but are also applicable to any single phase.

### Typical crystal numbers (No)

Purely on dimensional grounds the total number of crystals in a rock is given by

(1)
$No=CN(J0G0)3/4$
where CN is a constant, which is of order one in numerical size and depends weakly on the exact crystal shape and model of crystallization (i.e. exponential increase in nucleation with time, constant growth rate, etc.). This relation makes sense physically. A high rate of nucleation requires only limited crystal growth, whereas the greater the rate of crystal growth the fewer crystals needed to solidify the melt.

### Typical crystal size (Lo)

Similarly, the typical crystal size (e.g. diameter) for crystallization of this unit volume of magma is

(2)
$L0=CL(G0J0)1/4$
where CL is a constant similar to, but distinct in size and behavior from CN. Typical crystal sizes are directly proportional to growth rate and inversely proportional to nucleation rate. Although the size of any single crystal may depend only on the product of growth rate and total growth time (see below), nucleation affects the total growth time by using up the available melt. It is important to realize that crystal size and population depend solely, phase equilibria willing, on the single parameter (Jo/Go), which functions essentially as a single variable. Neither Jo nor Go can be found explicitly from measures of Lo and No alone. It is also clear that (1) and (2) are related, because the total number of crystals that can occupy any given volume depends on how large they are. Both of these relations have been known for some time (Winkler, 1949; Shaw, 1965; Brandeis & Jaupart, 1987a, 1987b; among others).

### Characteristic crystallization time (tc)

The time to fully crystallize a given volume of melt in which crystals are nucleating and growing at characteristic rates of, respectively, Jo and Go, is given by

(3)
$tc=Ct(G03J0)−1/4$
where Ct is another constant similar to those introduced already. Crystallization time is shortened with increases in either growth or nucleation rate and vice versa. But if both Jo and Go are similar functions of time, tc will depend more strongly on Go than on Jo, because of the exponent. In rocks, however, Jo has a much larger dynamic range (Cashman, 1993), which makes tc more or less equally sensitive to both Jo and Go. As each of these relations involves only Jo and Go, it will be seen later that the previous two constants (CN and CL) must also be determined by Ct.

This time scale is not explicitly dependent on a cooling model, which usually dominates such treatments, but arises purely from the kinetic processes of nucleation and growth. As mentioned at the outset, because many rocks are, to a very good first approximation, holocrystalline, nucleation and growth evidently automatically adjust to the thermal regime to ensure complete solidification. This kinetic time scale may thus be set independently by an appropriate thermal model containing the effect of latent heat. And because of the large dynamic range of nucleation, the coupling of kinetics and cooling is largely through adjustment of the nucleation rate.

## Crystallinity with Time

The most popular and straightforward model for predicting crystallinity as a function of time from given rates of nucleation and growth is the so-called Avrami method (Avrami, 1939, 1940). Although the original derivation pays special attention to the issue of nucleation in regions already occupied by crystals, called phantom nuclei and crystals (e.g. Kirkpatrick, 1983), the same equation can be found simply from a strict consideration of conservation of mass. The variation in crystal fraction (ϕ) with time is given by the integration over time of the product of crystal volume fraction and nucleation rate. Because crystals can appear at any time t′ during the crystallization history (tc), the integrations are nested, giving

(4)
$φ(t,J,G)=1−exp{−4π3∫otJ(t′)[∫t′tG(t)dt3]dt′}$
where now the rates of nucleation and growth are given by general functions of time as, respectively, J(t) and G(t), and the shape of the crystals is assumed to be spherical (hence the factor 4π/3). An alternative crystal shape can be chosen without adding any complexity unless shape is a function of time, which is clearly of secondary importance for the present consideration.

The cubed term in (4) represents the volume of a single crystal as a function of time, which stems from the change in crystal size with time

(5)
$L(t)=CA∫t′tG(t)dt$
where CA is a shape factor. For a spherical crystal, for example, CA = 1 and the volume of this crystal is
(6)
$V′s=43π[∫t′tG(t)dt3]⋅$
The total number of these crystals appearing at each time t′ is dictated by the rate of nucleation.

It is especially useful to consider some specific solutions to (4), and for future convenience it is good to consider some fairly general representations for the rates of nucleation and growth. Considering the exponential functions common to kinetics, two useful functions are

(7)
$J(t)=Joexp(at/tc)$

(8)
$G(t)=Goexp(bt/tc)$
where a and b are constants, which can be of either sign, and tc is as introduced already by (3) and further defined below.

There are two features of these equations that may at first seem unsatisfactory. First, they do not allow for maxima, when from kinetic theory maxima might be expected. Second, there is no explicit allowance for a time lag between nucleation and growth as might be expected from considerations of the effect of undercooling. As fundamental as these effects might be there is no real sign that they are important in the crystallization of magma. I will show later that almost all igneous CSDs indicate an apparent exponential increase in nucleation with time without encountering any maximum. And although this is harder to prove for growth rate, the evidence from Makaopuhi lava lake does not suggest a maximum or minimum in growth with time or degree of crystallinity (Kirkpatrick, 1978; Cashman & Marsh, 1988). The role of undercooling in delicately regulating the relative timing of the onset of nucleation and growth is probably of minor importance in real magmas. This may be because magmas are seldom, if ever, superheated, always possess abundant nuclei, and are often multiply saturated, which leads to heterogeneous nucleation induced by the rate of cooling and not simply undercooling (Cashman, 1993). This is not to say that undercooling does not exist, but only that it is commonly asymptotically small. Moreover, we will see below that the results are not sensitive to the exact shape of these exponential functions.

Equation (4) can be presented in a more convenient and compact form by substituting (7) and (8) and defining a relative or nondimensional time x = t/tc, which gives

(9)
$φ(x,a,b)=1−exp {−43πJoGo3tc4[∫0xexp (ax′)(∫x′xexp (bx)dx)3dx′]}.$
Substituting for tc from (3) and representing the two integrals as
(10)
$f(x,a,b)=[∫0xexp (ax′)(∫x′xexp(bx)dx)3dx′]$
allows (9) to be written as
(11)
$φ(x,a,b)=1−exp[−43πCt4f(x,a,b)]⋅$
The constant Ct can be found by defining the crystal fraction φc at which crystallization is considered complete; for example, at φc = 0.95 nondimensional time is: x ∼ 1 [i.e. when f(x,a,b) = f(1,a,b)], such that from (11)
(12)
$Ct=[−ln(1−φc)4/3πf(1,a,b)]1/4$
and
(13)
$φ(x,a,b)=1−exp{[ln(1−φc)f(1,a,b)]f(x,a,b)}⋅$
Once f(x,a,b) is known both φ(x) and Ct are also known.

The constants in (1) and (2) can now be found by employing the definition of nucleation relating the change of the total number of crystals (N) with time to the continuing effective nucleation rate J′(t), which is the nucleation rate defined relative to the volume of the whole system (i.e. crystals plus liquid). That is,

(14)
$dNdt=J′(t)$
In terms of the nucleation rate [J(t)] based only on the fraction of liquid present at any time, J′(t) = [1 – φ(t)]J(t), and (14) becomes
(15)
$dN(t)=[1−φ(t)]J(t)dt.$
In laboratory studies nucleation is commonly measured relative to the liquid content (i.e. 1 − φ), whereas in CSD studies of actual rocks it is often more convenient to report nucleation rate relative to the total (local or specific) volume of the sample. Using (11), (15) can be integrated
(16)
$N=Ct(J0G0)3/4∫01[1−φ(x)]exp(ax)dx$
where the characteristic time tc has been replaced by (3).

The constant CN in (1) is found to be

(17)
$CN=Ct∫01[1−φ(x)]exp(ax)dx.$
The constant CL in (2) can be found by considering a unit volume containing N crystals, each of a volume (4π/3)Lo3, where each crystal is taken to be a sphere. Then, using (1)
(18)
$CL=(43πCN)−1/3$
where CN is as given by (17). [Given the full distribution of crystal sizes and numbers (see below), there are clearly other measures of CL, but this is perhaps the most direct.]

All the constants and characteristic scales are now known as long as the integral in (10) can be found, which is always possible numerically. Thus each constant (i.e. CN and CL) is some multiple of Ct. Using the functions (7) and (8) various useful analytical results are possible.

### Constant rates of nucleation and growth (a = 0, b = 0)

These are the conditions leading to the well-known Johnson–Mehl–Avrami (JMA) equation (e.g. Kirkpatrick, 1983). The exponential constants a and b in equations (7) and (8) are both zero and (10) can be integrated to give

(19a)
$f(x,0,0)=x44.$
At complete crystallinity x = 1 and f(1,0,0) = 1/4, from which all the characteristic constants (CL, CN, and Ct) can be found (see Table 1; all symbols used are listed in Table 2). The increase in crystal fraction φ(x) with nondimensional time x is shown by Fig. 1. The results for a number of other kinetic models are also given by Table 1 and shown by Fig. 1; dimensional time t is recovered from x through t = xtc, where tc is given by equation (3). Ultimately, of course, the crystallization time is set by the rates of nucleation and growth. The general form of the increase in crystallinity is sigmoidal, which is broadly similar to the variation in crystallinity found for actual magmas (e.g. Marsh, 1981).

Fig. 1.

(Upper) the variation in crystallinity (reported as crystal fraction φ) with dimensionless time (x = t/tc). The inset table gives the values of the exponential constants a and b for, respectively, nucleation and growth. (Lower) crystallinity as a function of the characteristic time constant (Ct) for complete solidification. The numbers against each curve relate to the model numbers from Table 1 and the inset of the upper figure. It should be noted that the solidification time shortens and the solidification interval narrows as nucleation rate and/or growth rate increases more strongly with time.

Fig. 1.

(Upper) the variation in crystallinity (reported as crystal fraction φ) with dimensionless time (x = t/tc). The inset table gives the values of the exponential constants a and b for, respectively, nucleation and growth. (Lower) crystallinity as a function of the characteristic time constant (Ct) for complete solidification. The numbers against each curve relate to the model numbers from Table 1 and the inset of the upper figure. It should be noted that the solidification time shortens and the solidification interval narrows as nucleation rate and/or growth rate increases more strongly with time.

Table 1:

Kinetic models and characteristic constants

Model a b f (1,a,bCN CL Ct
1/4 0.895 0.644 1.448
0.31 1.241 0.577 1.372
−1 0.055 1.544 0.537 2.114
2.179 0.895 0.644 0.843
0.396 1.82 0.508 1.291
0.725 4.73 0.370 1.110
1.585 18.938 0.233 0.913
4.182 60.791 0.158 0.716
Model a b f (1,a,bCN CL Ct
1/4 0.895 0.644 1.448
0.31 1.241 0.577 1.372
−1 0.055 1.544 0.537 2.114
2.179 0.895 0.644 0.843
0.396 1.82 0.508 1.291
0.725 4.73 0.370 1.110
1.585 18.938 0.233 0.913
4.182 60.791 0.158 0.716
Model Function [f(x,a,b)]
f(x,0,0) = x4/4
f(x,1,0) = 6exp(x) − (x3+3x2+6x+6)
$f(x,1,1)=32(1+2x)exp (x2)+12[2+exp (x3)−6exp (x2)]exp (x3)$
f(s,1,1) = [exp(x)−1]4/4
f(x,2,0) = f(x,1,0)/24 as for model 2 above but with x = 2x
f(x,4,0) = f(x,1,0)/44 as for model 2 above but with x = 4x
f(x,6,0) = f(x,1,0)/64 as for model 2 above but with x = 6x
f(x,8,0) = f(x,1,0)/84 as for model 2 above but with x = 8x
Model Function [f(x,a,b)]
f(x,0,0) = x4/4
f(x,1,0) = 6exp(x) − (x3+3x2+6x+6)
$f(x,1,1)=32(1+2x)exp (x2)+12[2+exp (x3)−6exp (x2)]exp (x3)$
f(s,1,1) = [exp(x)−1]4/4
f(x,2,0) = f(x,1,0)/24 as for model 2 above but with x = 2x
f(x,4,0) = f(x,1,0)/44 as for model 2 above but with x = 4x
f(x,6,0) = f(x,1,0)/64 as for model 2 above but with x = 6x
f(x,8,0) = f(x,1,0)/84 as for model 2 above but with x = 8x
Table 2:

Symbols used

Symbols
C characteristic constants (see subscripts)
D fractal dimension
f(x,a,bfunction describing the Avrami integral [see equation (10)
G crystal growth rate (length/time)
J nucleation rate (number/unit volume/time)
K rate constant
N number of crystals per unit volume
n population density of crystals (number/unit volume/size)
S numerical value of CSD slope
t time; when with a subscript, a characteristic crystallization time
V volume (cubic length)
x nondimensional time (= t/ tc
φ crystal fraction
τ residence time or flushing time
Subscripts
shape factor
characteristic value of crystallization time or crystallinity
index
constant for crystal size
maximum crystal size
constant for crystal concentration
characteristic values of various parameters
solid
constant for crystallization time
Superscripts
a exponential constant for nucleation rate function
b exponential constant for growth rate function
o nucleation density (e.g. no')
′ nucleation rate based on system volume; also nondimensional crystal size
Symbols
C characteristic constants (see subscripts)
D fractal dimension
f(x,a,bfunction describing the Avrami integral [see equation (10)
G crystal growth rate (length/time)
J nucleation rate (number/unit volume/time)
K rate constant
N number of crystals per unit volume
n population density of crystals (number/unit volume/size)
S numerical value of CSD slope
t time; when with a subscript, a characteristic crystallization time
V volume (cubic length)
x nondimensional time (= t/ tc
φ crystal fraction
τ residence time or flushing time
Subscripts
shape factor
characteristic value of crystallization time or crystallinity
index
constant for crystal size
maximum crystal size
constant for crystal concentration
characteristic values of various parameters
solid
constant for crystallization time
Superscripts
a exponential constant for nucleation rate function
b exponential constant for growth rate function
o nucleation density (e.g. no')
′ nucleation rate based on system volume; also nondimensional crystal size

### Exponential nucleation and constant growth (a = 1, 2, 4, 6, 8, b = 0)

For these conditions only, the exponential constant b in (8) and (10) is zero, and a is nonzero. Integration of (10) gives

(19b)
$f(x,a,0)=(1/a)4{6exp(ax)−[(ax3+3(ax)2+6ax+6)]}.$
This is a flexible result, in terms of nucleation rate, and results for a series of models where a = 1, 2, 4, 6, and 8 are given by Fig. 1 and Table 1, which also gives the values of f(1,a,0) necessary for equation (13).

### Exponential nucleation and growth (a = 1, b = 1)

Here nucleation and growth both increase exponentially with time.

(20)
$f(x,1,1)=[exp(x)−1]44$
and f(1,1,1) = 2.179. The variation in crystallinity (Fig. 1) is almost identical to model 7 (Table 1) with a = 6, b = 0, but, as will be shown later, the number of crystals and crystal size are almost identical to the simple JMA result (i.e. model 1 with a = 0, b = 0).

### Exponential nucleation and exponentially decreasing growth (a = 1, b = 1)

Here the rate of nucleation increases and the rate of crystal growth decreases with time; even though there are increasing numbers of nuclei they grow slower and slower with time.

(21)
$f(x,1,−1)=3(1+2x)2exp(x2)+12[2+exp(x3)−6exp(x2)]exp(x3)$
where f(1,1,−1) = 0.055. The pattern of increasing crystallinity is distinct from all other models (Fig. 1).

### Overview of crystallization models

Although the variations in crystallinity with nondimensional time of Fig. 1 (upper) are broadly as expected given the underlying assumptions of each model, they are somewhat misleading when plotted in this fashion. It would appear that the models with high nucleation rates crystallize more slowly than those with slower growth and nucleation rates. This misperception is purely due to the nondimensional time scale, which has a different characteristic time (tc) for each model. A more meaningful and intuitively appealing time scale is that defined by the characteristic constant for crystallization time (Ct; see also Table 1), which for any given characteristic nucleation and growth rates gives the relative time to complete crystallization [Fig. 1 (lower)]. The stronger the increases in nucleation and growth the sooner crystallization is complete. Model 8 (a = 8, b = 0) fully crystallizes over a relative time of Ct = 0.716, whereas model 3 (a = 1, b = −1) takes approximately three times longer (Ct = 2.114). Across all the models a reasonable choice for crystallization time might be Ct ∼1.

The larger the nucleation rate the more crystals are produced, but as the amount of melt is finite the final crystal size will be smaller. These effects are reflected in the constants CL and CN, which are shown along with Ct in Fig. 2. As the characteristic number of crystals increases from about CN = 0.6 to 60, at the highest nucleation rates, crystal size decreases from CL = 0.75 to ∼0.15, as nucleation increases by ∼103.5.

Fig. 2.

The characteristic constants for crystal size (CL, upper), solidification time (Ct, lower), and number of crystals (CN, x-axis) for a variety of models of nucleation and growth as discussed in the text.

Fig. 2.

The characteristic constants for crystal size (CL, upper), solidification time (Ct, lower), and number of crystals (CN, x-axis) for a variety of models of nucleation and growth as discussed in the text.

Overall, strong changes in nucleation rate with time have a relatively mild effect on the ensuing crystallization time and the size of the crystals; it does have an appreciable effect on the size spectrum of crystals, which will show up in the CSD. The effect of growth rate is also not particularly strong except in determining the crystallization time. These analytical results involving no cooling model can be compared with the numerical results of Brandeis & Jaupart, (1987b), which employed fully coupled conductive cooling and crystallization. They found the shape of the function representing growth rate to have no influence on the results, whereas nucleation, when compared with crystal size in rocks, is best represented by a constant rate unless near a contact. Although this may be so in the mean, we will find below that when the actual size spectrum of crystals is considered, as measured by a CSD, nucleation must vary exponentially with time.

## Crystal Size Distributions (CSDs)

By considering the crystal size distributions associated with the kinetic models presented already, which are each a batch system, another measure is available by which to gauge the value of each model. Beginningagain with equation (15) and recalling that CSDs aremeasured by a crystal population density (n) definedby n = dN/dL, where N is crystal number and L is crystal size (Marsh, 1988; Randolph & Larson,1988)

(22)
$n(t)=dNdL=[1−φ(t)]J(t)G(t)$
where it is recalled that n measures the crystal population density of the whole system and J is the nucleation rate of the liquid. Substituting equations (7) and (8) into (22), and rearranging and taking the natural log,
(23)
$ln[n(x)n0]=ln[1−φ(x)]+(a−b)x$
where the initial (i.e. at time t or x = 0) nucleation density no defined using (22) above as
(24)
$n0=[1−φ(x=0)]J0G0=J0G0$
has been used. This convention of defining no at x = 0 is different from that used in most CSD studies where no (then denoted as No) is defined as the final (i.e. x = 1) nucleation density, which pertains to the nucleation density at the smallest (i.e. L = 0) crystal size. The present convention (or standard state) allows no to be unambiguously related to Jo and Go. Thus no relates to the initial nucleation density, in contrast to the final nucleation density no.

Given any real CSD it is a simple matter to convert from no to No, and given a nucleation rate as described by (7), ln(n/No) + a = ln(n/no). The disadvantage in using no is that it relates to the far end of the CSD, where L = Lm, and measurement uncertainty can be large. Substituting (11) into (23) gives

(25)
$ln [n(x)no]=−43πCt4f(x,a,b)+(a−b)x$
and by substituting for Ct from (12)
(26)
$ln [n(x)no]=[ln (1−φc)f(1,a,b)]f(x,a,b)+(a−b)x.$
This is the batch CSD as a function of the log of the liquid fraction (the first term on the right) and the product of the difference in the exponential constants (ab), describing the rates of nucleation and growth, and nondimensional time (x = t/tc; second term on the right). (It actually describes the rate of production of nuclei, but because these must grow we can take it also to give the batch CSD.) When plotted as a function of time both terms contribute to the CSD, but it is the last term that gives the overall slope to the CSD. If a = b, for example, this term vanishes and the CSD is defined only by the variation in liquid fraction over time; the CSD is essentially horizontal with no slope. The larger the difference in ab, the steeper the CSD slope. It is of utmost importance to realize here that the CSD slope is set by the rates of nucleation and growth (i.e. ab). For constant growth rate, for example, b = 0 and the slope is determined solely by the nucleation rate. If the nucleation rate is constant (i.e. a = 0) the batch CSD will be flat (zero slope) as long as G is also constant. And if nucleation decreases with time (i.e. a <0), the CSD slope will be positive if plotted against crystal size (see below). That individual rock CSD slopes are never flat or positive strongly implies that, for this strictly batch model, the effective nucleation rate always increases locally on the scale of the sample. Much later we shall see that the change in nucleation rate suggested by successive samples (e.g. successive lavas) sometimes increases, is sometimes contant, and sometimes decreases. This apparent inter-sample change should not be confused with the nucleation experienced by local samples, for these apparent nucleation changes more reflect differences in cooling regime. Thus if apparent nucleation changes were used to describe the local nucleation rate (e.g. a <0), the resulting CSD would be exceedingly unusual. We will return to these points again.

CSDs for rocks are always shown as a function of crystal size and not time as in equation (26). Later we will see how this arises from the CSD continuity equation considered jointly with the present calculation taken as merely a nucleation condition or side condition. But a common CSD result can also be extracted directly from (26), and this gives insight into later, more involved results. Thus, for kinetic models where growth rate is constant (i.e. b = 0) crystal size relates directly to time, but the largest crystals are those nucleated at t = 0 (or x = 0) and the smallest are those just nucleated as x reaches unity. The explicit relation between growth time and crystal size is found by considering (5) with CA = 1, which for constant G (= Go) becomes after integrating

(27a)
$L=Go(tc−t′)$
where t′ is the time within the solidification interval 0 ≤ t′ ≤ tc when the crystal nucleated and began growing. Rewriting this to put it in a nondimensional form
(27b)
$LLm≡L′=1−x$
where Lm(= Gotc) is the maximum possible crystal size. Wherever x appears in (26) it can be replaced by 1 − L′ and the resulting CSD equation can be plotted against nondimensional crystal size L′; alternatively, equation (26) can simply be plotted against 1 − x to give the CSD as a function of nondimensional crystal size L′. It is important to emphasize that (27b), strictly speaking, holds only for G = Go. When G is a function of time or size the relationship between L′ and x will be more involved, but for values of b near zero (27b) is a reasonable approximation (see below).

In general, for b ≠ 0, the maximum crystal size (Lm) is found from equations (5) and (8), which in dimensional form gives

(28)
$Lm=Gotc∫01exp(bx)dx=Gotcb[exp(b)−1]$
or substituting for tc from (3)
(29)
$Lm=[exp(b)−1]bCt(GoJo)1/4$
where Ct is as given by (12). In the special case limit of constant growth rate, b = 0 and, using L'Hospital's rule,
(30)
$Lm=Gotc=Ct(GoJo)1/4$
which in comparison with (2) shows that for the largest crystal size CL = Ct, under this restriction of constant growth rate, but this is certainly not the typical crystal size. That is, in comparison with (18), where all crystals are of the same size, for the maximum size crystals considered here CN = Ct − 3/(3/4π), which gives constants for the number of crystals different from that found from (17).

### CSDs

Batch CSDs calculated from equation (26) for the previous kinetic models with ab = 0, 2, 4, 6, 8 and 10 (where b = 0 in all cases) are shown by Fig. 3. Each CSD represents the final population density for the batch system after 95% crystallization (i.e. φc = 0.95). The diminishing liquid fraction during the last stages of crystallization causes the pronounced decrease in the CSD at the smallest crystal sizes. That is, the late-stage effective nucleation rate [see equation (15)] becomes vanishingly small because of the diminishing volume of melt even though the nucleation rate itself is increasing exponentially. The slope of each CSD reflects the numerical difference in the exponents a and b describing the functions for, respectively, nucleation and growth given by equations (7) and (8).

Fig. 3.

The final (i.e. φc = 0.95) crystal size distribution in a batch system as a function of the difference ab describing the exponential increases in nucleation (a) and growth (b) rates. Only for larger values of ab does the CSD resemble those of igneous rocks. The decrease in population density at smaller crystal sizes reflects the loss of melt with time.

Fig. 3.

The final (i.e. φc = 0.95) crystal size distribution in a batch system as a function of the difference ab describing the exponential increases in nucleation (a) and growth (b) rates. Only for larger values of ab does the CSD resemble those of igneous rocks. The decrease in population density at smaller crystal sizes reflects the loss of melt with time.

The CSDs typically observed in igneous rocks show population densities [i.e. ln(n/no)] that vary by 6–8 across the observed size range [see, for example, the summary Fig. 24 of Cashman, (1990)]. That flatter CSDs have not so far been found suggests that nucleation varies much more strongly than growth rate during crystallization, which favors the larger values of ab. This is not at all surprising judging from the strong reduction in crystal numbers commonly observed inward from the contacts of high-level intrusions. With this in mind, it is interesting to consider the maturation of the CSD with progressive crystallization.

The progression of the batch CSD during crystallization is found by varying the characteristic crystallization time through the variation of φc in equation (26) in coordination with the commensurate variation in x or t/tc. Because φc varies in a nonlinear fashion with x, for purposes of calculation it is simplest to make the approximation that x ∼φc in the last term of (26). The CSDs after 5, 25, 50, 75, and 95% crystallization and ab = 8, calculated in this manner, are shown by Fig. 4. It is emphasized that these are, strictly speaking, approximate results. The CSDs are, for the most part, parallel to one another. The maximum crystal size increases systematically with increasing nucleation rate. Only in the later stages of crystallization, once φc is greater than ∼75%, is the CSD markedly different from that at the smallest degrees of crystallization.

Fig. 4.

The development of the batch crystal size distribution with increasing crystallinity or time for a system where ab = 8. It should be noted that the slope is generally constant and that the effect of decreasing melt on nucleation becomes prominent beyond ∼50% crystallization.

Fig. 4.

The development of the batch crystal size distribution with increasing crystallinity or time for a system where ab = 8. It should be noted that the slope is generally constant and that the effect of decreasing melt on nucleation becomes prominent beyond ∼50% crystallization.

Strongly linear CSDs are characteristic of volcanic rocks where φc is always less than ∼50%. Thus in general for volcanic rocks there is little practical difference in how nucleation rate is defined; it can be defined relative to either the whole system volume or only the liquid volume with little loss of generality [see equation (15)]. It is important to realize here that these CSDs are representative for any exponential factors a and b as long as their difference is eight. That is, even though in calculating Fig. 4 the growth rate has been taken to be constant (i.e. b = 0), for any other value of b the CSD will not be too much different; the principal effect will be to change φc(t), which will change the curvature to the CSD as a result of a nonlinear growth rate.

It will be shown later that the log-linear form of the CSD is basically due to an exponential change in nucleation with time which is propagated across the CSD diagram at a rate determined by the growth rate. The overall slope of a straight [i.e. ln(n)-linear] CSD is thus determined by the overall amplitude of the effective nucleation event (ab) and the maximum crystal size Lm, which is also obvious by combining equations (7) and (38). The meaning of Lm is that it is a measure (through the growth rate Go) of the overall time of the nucleation event. The slope of the CSD is thus determined by the nucleation exponent a and the growth rate Go. With all else being equal, CSDs having similar rises in nucleation rates reflect different growth rates by showing different slopes. And successive CSDs of the same batch system that do not have the same slope cannot have a constant growth rate. There are comagmatic CSDs that do not have the same slopes, but instead fan and suggest variations in growth rate with time or crystal size. I will return to this important issue later.

### Actual crystal size

There is probably no more fundamental length scale in geology than crystal size. Any model must produce reasonable crystal sizes in reasonable times. Crystal size is the product of growth rate and crystallization time [see (29) and (30)], but crystallization time depends on nucleation rate.

Crystal size, using (29) is given by

(31)
$L(x)=Ct(GoJo)1/4[exp(bx)−1b]$
and substituting for Ct from (12)
(32)
$L(x)=[−ln(1−φc)4/3πf(1,a,b)]1/4(GoJo)1/4[exp(bx)−1b]⋅$
When growth rate is constant (b = 0), as before using L'Hospital's rule on the last term on the right,
(33)
$L(x)=[−ln(1−φc)4/3πf(1,a,b)]1/4(GoJo)1/4x⋅$
As (5) really describes the effective radius of the crystal, the characteristic crystal size (diameter or length) will be twice that given by these equations. Crystal sizes calculated over a range of values of Go and Jo using (33) for the models of Table 1 are shown by Fig. 5. For a constant growth rate of 10−9cm/s, for example, to grow a 1 mm diameter crystal calls for a nucleation rate of between 10−3 and 10−5cm−3 s−1.

Fig. 5.

The variation in crystal size (radius) with nucleation rate when growth rate is constant (i.e. b = 0) and nucleation increases exponentially with time (i.e. a ≠ 0).

Fig. 5.

The variation in crystal size (radius) with nucleation rate when growth rate is constant (i.e. b = 0) and nucleation increases exponentially with time (i.e. a ≠ 0).

Nucleation rates of this magnitude are common for volcanic rocks. The plagioclase of Makaopuhi lava lake shows typical nucleation rates of 10−3 (number/cm3 s) and the basaltic lavas of Dome Mountain, Nevada, show rates of 10−3–10−5 (number/cm3 s) (Cashman & Marsh, 1988; Resmini & Marsh, 1995). The crystals of Makaopuhi are much smaller than 1 mm because they are not yet at their final size.

The main difficulty in determining either Jo or Go from rock CSDs is that they appear in all formulae together, making it difficult to separate and determine either parameter. Even in the simplest formula for maximum crystal size [equation (30)], either Jo or tc must be independently known to estimate Go from crystal size (i.e. Go = Lm/tc). Growth rate also cannot be estimated with any reliability from kinetic theory, which involves guessing the degree of undercooling (e.g. Brandeis & Jaupart, 1987a; Solomatov, 1995). As true undercooling probably has little meaning in magmatic systems, growth rates have been correlated with cooling rate by Cashman, (1993). In the spirit of the Crystallization Axiom (see the Introduction), the cooling rate sets the local crystallization time (tc). These clear log-linear correlations do not have the slopes expected from the exponents of equations (1), (2), and (3). Suffice it to say here that other scalings are possible (e.g. nucleation not a time-dependent process) that can give more agreeable exponents. But the CSDs of the rocks themselves seem to demand a time-dependent nucleation process, which would suggest that Gotc correlations must explicitly include nucleation rate and not simply a measure of tc from heat transfer models.

Although the results of Fig. 5 are for constant growth rates (i.e. b = 0), the effect of b ≠ 0 can be readily assessed using equation (28), and some results are shown by Fig. 6. For b = ±1, the effect is not large, but for large values of |b| the effect can be large. Although there is little firm evidence of the dependence of G on time, what information there is from CSDs suggests that G is certainly not a strong function of time or crystal size (Cashman & Marsh, 1988; Marsh, 1988; Resmini & Marsh, 1995).

Fig. 6.

The variation in the scale factor for maximum crystal size as a function of the magnitude of increase in growth rate as function of time as measured by the exponential constant b.

Fig. 6.

The variation in the scale factor for maximum crystal size as a function of the magnitude of increase in growth rate as function of time as measured by the exponential constant b.

The characteristic time of crystallization, particularly if the growth rate is constant, sets both the overall time of crystal growth and of course the maximum size of the ensuing crystal. The crystallization time is determined by the product Go3Jo as given by equation (3). Even if the style of nucleation is radically different (e.g. nucleation is not a rate at all, but a fixed set of nuclei), this product is the governing parameter, although the fractional exponent of the whole group will change from −1/4 to −1/3.

Representative times calculated using (3) and (12) for constant growth rate (i.e. b = 0) are given by Fig. 7; along the upper axis nucleation rates are given for the specific growth rate of 10−9 cm/s. The smaller the nucleation rate the longer the time of crystallization. For nucleation rates of J = 10−310−6 cm−3 s−1, tc ≈1–10 years, and the effect of increasing the exponential increase in nucleation rate from a = 1 to a = 8 decreases tc by a factor of two, which is also clear from Ct in Fig. 2.

Fig. 7.

The time for solidification as a function of the factor Go3Jo for constant growth rate and nucleation functions described by, respectively, a = 1 and a = 8.

Fig. 7.

The time for solidification as a function of the factor Go3Jo for constant growth rate and nucleation functions described by, respectively, a = 1 and a = 8.

## Summary of Crystallization Models

Beginning with a set of equations for N, L, and tc based purely on scaling or dimensional arguments, it is possible to develop an internally consistent and petrologically reasonable set of results by using an Avrami approach to calculating crystallinity as a function of time. A wide range of models employing exponential functions describing the rates of nucleation and growth show fairly similar results. All of these are analytical models and represent only a small subset of a vast array of models that can give broadly similar results.

Perhaps the most useful results are those relating to the CSDs associated with these models. The CSDs of actual rocks provide an added dimension by which to discriminate between models. This is the critical bridge between Avrami-style kinetic calculations and the real world of magma crystallization. Against this background of formal models, I now turn to interpreting the general kinetic features of crystallization of igneous systems using CSDs from various magmatic systems. I begin below with a brief introduction of the basic formalism of CSDs in batch and open systems.

## CSD Systematics in Batch and Open Systems

The unusual character of many igneous CSDs is best appreciated by gauging them against theoretical CSDs where the conditions of nucleation and growth are known exactly. It is also important to see that the formal population balance for CSDs in batch systems is entirely equivalent to that arrived at already using Avrami crystallization. In a batch system of volume V the governing population balance is described by a basic continuity equation (Randolph & Larson, 1971, 1988; see p. 59):

(34)
$∂n∂t+∂(Gn)∂L+nd(lnV)dt=0$
where all symbols are as before. It should be noted that the last term on the left is the total derivative for the change in volume of the system, which when written out becomes
(35)
$nd(lnV)dt=n(∂lnV∂t)+n(∂lnV∂t)G.$
Upon substitution, (34) becomes
(36)
$∂(nV)∂+∂(GnV)∂L=0.$
Because n and V appear as a product it is convenient to define the population density relative to the volume of the whole system, which, neglecting the volume change of solidification, can be taken as constant. This is consistent with the definition of n given by (22).

Under this condition, (36) becomes

(37)
$∂n∂+∂(Gn)∂L=0$
which describes the nonsteady CSD in a batch system.

The general form of (37) when G is constant (G = Go), or only a function of time and not crystal size, is that of a wave equation, which is solvable using the method of characteristics (Randolph & Larson, 1988, p. 55). Given any initial CSD, this equation predicts the forward propagation in time of this CSD at a rate Go. On a plot of ln(n) vs L, for example, each crystal moves horizontally across the diagram according to (37) and the crystal size equations (28) and (30). Because (37) merely describes the conservation of the crystal population density, it really contains, as it stands, little useful kinetic information unless it is considered in conjunction with equations (28) and (22). That is, the set of equations describing the initiation and evolution of the CSD is

(38a)
$∂n∂±Go∂n∂L=0$

(38b)
$L(t)=CA∫t′tG(t)dt$

(38c)
$n(L,t)=[1−φ(t)]J(t)G(t)$

along with a relation describing φ(t). These are an essential difference between assuming G to be independent of time or independent of crystal size L. If G is a function of time, the growth rate of all crystals regardless of size will vary with the age of the system. If, on the other hand, G is dependent on L, the growth of each crystal depends on the age of that crystal. Both situations are possible to some degree. A changing thermal regime may ‘globally’ modify G, and hence G = G(t). Size-dependent growth G = G(L), on the other hand, could very well be due to crystal crowding, boundary-layer thickening, or any number of other effects (see later). When G is stated to be constant in the following, it is assumed that G is not a function of either t or L.

The meaning of the last equation of this set is also important to emphasize. Recalling equation (15), (38c) comes from the definitions of nucleation rate and population density. That is, for L=0.

(39a)
$J′(t)=dNdt=dNdLdLdt$
or
(39b)
$J′(t)=n(t)G(t)$
where J′(t) is the rate of nucleation per unit volume of the system. Because nucleation occurs in the melt, regardless of whether or not nucleation is heterogeneous, it is useful in batch systems to base nucleation on the amount of melt present [i.e. J(t)]. This requires a model for the loss of melt as a function of nucleation and growth rates, which is available, for example, from the preceding Avrami model. Then,
(40)
$J′(t)=[1−φ(t)]J(t)$
and (39b) becomes
(41)
$n(t)G(t)=[1−φ(t)]J(t)$
or at a time t = 0
(42)
$noGo=Jo.$
The important aspect of this equation in batch systems is that this holds for any crystal size and not just for the smallest, newly born crystals. That is, in a batch system where each crystal, once nucleated, is preserved until complete solidification, the population balance at all times and all crystal sizes obeys (42). Only if the crystals grow dispersively where the growth of each cohort of crystals does not depend only on time or if crystal aggregation occurs would this relation be violated. It is also useful to notice that for constant G [≠ f(t) or f(L)], as mentioned earlier n/no = J/Jo. In essence, at constant G, n and J are interchangeable.

The importance of (38c) and its relation to (38a) is further emphasized by rewriting the population balance (38a) in a form involving nucleation. Recalling that J = dN/dt, (38a) can be written as

(43)
$∂2N∂t∂L+∂(Gn)∂L=0$
now changing the order of differentiation in the first term
(44)
$∂J′∂L+∂(Gn)∂L=0$
rewriting gives
(45)
$∂∂L(J′+Gn)=0$
or
(46)
$J′+Gn=constant$
and if J′(L = 0) = 0, then
(47)
$J′+Gn=0.$
The relation of this to (38c) is clear, but here the sign is different. This is because (47) relates to (38a), which describes the propagation of a given CSD. That is, within a size range the population density can only change in time as a result of changes in Gn with L; only at L=0 does the sign change.

### Initial crystal-free systems

Returning to the set of equations (38), let us consider a melt free of crystals and nuclei that begins at a time t = 0 to cool and nucleate crystals at a rate J(t), which go on to grow at a constant rate Go. By rewriting the batch population balance equation [i.e. equation (38a)], it can be seen that the slope of the CSD gives a record of the nucleation history. Multiplying (38a) by Go/J′,

(48)
$1J′∂nGo∂t+Go2∂nJ′∂L=0$
and using Gon = J′ in the first and second terms, gives
(49)
$1Go∂ln(J′)∂t=−∂ln(n)∂L.$
The term on the right is the slope of the CSD on a conventional semi-log CSD ln(n) plot. Any slope on such a plot translates into an equivalent variation of effective nucleation rate with time; conversely, any variation in effective nucleation rate gives rise to a change in slope of the batch CSD. That is, any time variation in the effective nucleation rate changes the population density along the L = 0 axis, which then propagates along the L axis because of the growth rate. Some examples will clarify this point.

Let us consider, for example, a CSD exhibiting a constant slope S on an ln(n) vs L plot. Beginning with (49), let S = ∂ln(n)/∂L, and then

(50)
$dln(J′)=−GoSdt$
which can be integrated to yield
(51)
$J′=Joexp(−SGot)$
where the condition J′(t = 0) = Jo has been used. Because the actual CSD slope is negative (i.e. S <0), nucleation increases with time. As seen earlier, a linear CSD relates to an exponential increase in nucleation. Some schematic variations in nucleation with time and the effects on the associated CSD are shown by Fig. 8, which are also evident in the earlier results of Fig. 4. When the nucleation rate is exponential the resultant CSD is linear. If nucleation rate suddenly changes to a new exponential rate, the CSD, although still linear, will change. The resulting composite CSD will thus be kinked. Kinked CSDs record sudden changes in nucleation rate. Later we will also see that growth processes can also kink CSDs.

Fig. 8.

The relation of the CSD slope to the time variation in nucleation rate. The upper set describes a single exponential nucleation event and the lower, kinked CSD, reflects a sequence of two nucleation events.

Fig. 8.

The relation of the CSD slope to the time variation in nucleation rate. The upper set describes a single exponential nucleation event and the lower, kinked CSD, reflects a sequence of two nucleation events.

On the other hand, choosing the exponential variation in nucleation given earlier by equation (7) [and using relations such as (44)] and substituting this into the first term in (38a) yields

(52)
$dln(n)dL=−aG0tc$
or, because from (30)Lm = Gotc,
(53a)
$dln(n)dL=−aLm$
which shows that the slope in a batch system is a constant and is set by both the nucleation increase, as set by (38c), and the rate of growth of crystals as set by (38b). This is also true if growth rate varies with time according to (8), as has been demonstrated by the CSD of equation (23), whence the slope is given by (ab)/Lm. It is important to realize for these models that the slope is fixed for the entire batch crystallization history, although there may be curvature if b ≠ 0. The only possible serious modification of the slope comes from the eventual loss of all melt at advanced stages of crystallization where the factor [1 minus; φ(t)] becomes important (see Fig. 4), which has been ignored in these examples. As long as the nucleation rate increases smoothly with time (in an exponential sense) the slope of the associated CSD will be constant. This is illustrated by Fig. 9, which shows how the CSD develops under a constant growth rate. From this it is also clear how the continuity and nucleation rate equations arise.

Fig. 9.

The development of the CSD and CSD continuity equation as a result of the combined effects of an exponential increase in nucleation rate and crystals growing at a constant rate.

Fig. 9.

The development of the CSD and CSD continuity equation as a result of the combined effects of an exponential increase in nucleation rate and crystals growing at a constant rate.

The batch CSD for constant growth rate is thus a clear record of the history of nucleation, and, as mentioned earlier, the population density can be readily converted to nucleation rate using (39), with constant G, that is, n/no = J/Jo. Relative variations in n are equivalent to relative variations in J. The true range of variation of J is found by knowing at least a single value of J, which can be found from J = Gon, where specific values of n and Go are employed. This conversion produces a plot of ln(J/Jo) vs L, which can be converted to relative time using (27c). The age of the batch system is given by the largest crystal size, which can be used in (30) to give t = Lm/Go, and this age should correspond to the overall degree of crystallinity [i.e. φ(t)] of the system. Any CSD can be converted into an equivalent variation in the effective nucleation rate over the period of solidification.

To discriminate a batch system of this type from other systems, plots of slope vs Lm and slope vs intercept, i.e. ln[n(L = 0)] = ln(No), can be used (see Fig. 10). Because in a batch system the slope is always constant, each diagram shows a horizontal line. It will be shown presently that most other styles of systems show distinct variations and igneous rocks, even within the same system, often show a significant range of relations between slope, Lm, and intercept.

Fig. 10.

Possible diagnostic relations for batch (upper) and open systems (lower) among CSD slope, intercept, and maximum crystal size (Lmax).

Fig. 10.

Possible diagnostic relations for batch (upper) and open systems (lower) among CSD slope, intercept, and maximum crystal size (Lmax).

Before moving on, it is pertinent to inquire if the log-linear form of the CSD can be preserved when the nucleation rate is constant. With nucleation rate at steady state, the first term in (38a) drops out and we are left with ∂(Gn)/∂L = 0, or when expanded and rewritten

(53b)
$G=Goexp (−SL)$
If G is constant the usual CSD slope can only be zero (flat) as already noted. But if the CSD slope is chosen to be S, where S <0, then we find that G must vary as
(53c)
$G=Goexp (−SL)$
where Go = G(L = 0). As S < 0, the growth rate must increase exponentially with crystal size. This makes sense. That is, with a steady production of nuclei and a growth rate increasing strongly with size, the population density is increasingly thinned with L as crystals rapidly grow through the larger sizes. Although it is unlikely that single isolated crystals could ever follow such a growth rate function, it may be possible for crystals growing by annexation (see later) of smaller crystals or for crystals undergoing strong shape changes because of growth. Also, depending on the details of annexation or shape changes, the CSD slope could vary about some mean value; overall the slope would again probably be more or less constant.

In sum, in batch systems the downward sloping log-linear form of observed CSDs can be formed from either a nucleation rate that increases exponentially with time, with a constant growth rate, or a constant nucleation rate with an effective growth rate that increases exponentially with crystal size. The maximum crystal size, Lm = [ln(SGotc = 1)]/S, is still a measure of the age of the system, but the relation between crystal size and time is not linear. Because this batch system is somewhat unusual (but not necessarily unreasonable), when further discussing batch crystallization below the system in mind, unless otherwise mentioned, is that with a nonconstant nucleation rate.

### Systems with an initial distribution of crystals

If there is a hiatus in crystallization and further crystallization (i.e. nucleation) begins again after some interlude, or if the system simply inherits a population of stable ‘tramp’ crystals that remain suspended in the magma, the ensuing CSD will be a composite of the old and new CSDs. There are any number of possibilities, but a few reasonable ones are depicted schematically by Fig. 11.

Fig. 11.

The effect of various nucleation events on the resulting CSD. (a) A single nucleation event followed by annealing during no nucleation. (b) Double nucleation events (deep and shallow) separated by a hiatus in nucleation. (c) CSD curvature because of superexponential increases in nucleation rate or size-dependent growth rate. (d) The net result of the processes of (a), (b), and (c) on the resulting CSD.

Fig. 11.

The effect of various nucleation events on the resulting CSD. (a) A single nucleation event followed by annealing during no nucleation. (b) Double nucleation events (deep and shallow) separated by a hiatus in nucleation. (c) CSD curvature because of superexponential increases in nucleation rate or size-dependent growth rate. (d) The net result of the processes of (a), (b), and (c) on the resulting CSD.

If nucleation ceases because of phase equilibria, but growth is still possible, the CSD may also be affected by annealing or Ostwald ripening whereby the smaller crystals dissolve to the advantage of the larger crystals. The size range is reduced and if nucleation recommences, the succeeding CSD may show a distinct hump. This pattern might also appear if the magma were to entrain a population of ‘tramp’ crystals from a well-sorted or well-ripened cumulate bed.

These composite CSDs may still represent a batch system, but the histories of growth and nucleation may not be simple. Their first-order interpretation, however, is still possible through a segment-by-segment consideration of the CSD.

To place the batch system in further distinction, it is useful to consider briefly the CSDs in open systems.

### Nonbatch open systems

The population balance in open systems, which are the so-called MSMPR (i.e. mixed suspension, mixed product removal) systems of Randolph & Larson, (1988), considers a chamber of volume V filling and emptying at a rate Q, and the incoming magma contains no crystals. This population balance is given by (38a) with an additional term:

(54)
$∂n∂t+Go∂n∂L+nτ=0.$

The additional term (n/τ) represents the outflow containing a population of crystals that has a mean residence time τ in the chamber. Various scenarios for such systems have been previously discussed (Marsh, 1988), and the most revealing is the simple steady-state system (i.e. ∂n/∂t = 0), whereby (54) becomes

(55)
$ddL=−nGoτ.$
Upon rearranging and integrating
(56)
$ln (nno)=−LGoτ$
where the condition that n(L = 0) = no has been used.

Under the assumption of steady state, there is no time dependence and with G constant, the nucleation rate is [from (39b)] also constant. Yet the CSD is log-linear because of the recharge or flushing time (τ) of the system. That is, the mean residence time of a crystal in the system before it is evicted is τ, and its size is L = Goτ; the maximum crystal size is ∼10 times larger than this size, and this crystal is ∼10 times older. The parameter Goτ is given by the slope of the CSD, which depends on the recharge rate and the growth rate. It should be clear, however, that without discharge of crystals, with this constant nucleation rate, the CSD would be horizontal and have a slope of zero, which reflects an infinite residence time in a batch system. If the nucleation rate were to change, over a period of flushing times, depending on the characteristics of the new CSD, a new steady CSD would be established (see Marsh, 1988). Other processes such as crystal settling and crystal accumulation can also influence the CSD, but most such processes affect only a relatively limited range of crystal sizes and thus mainly distort or produce curvature in the CSD.

In terms of discriminating between a batch system and an open system, in the present steady, open system, both the slope and intercept are invariant with time, as is also the maximum crystal size. On the preceding discrimination plots, this system is a single point. Only if the system is unsteady does it form a trajectory leading to the steady-state end point (see Fig. 10).

### Summary of CSD in batch and open systems

In a batch system with exponentially increasing nucleation rate, the CSD migrates systematically up the ln(n) axis and across the L or size axis. The slope of the CSD remains constant regardless of t and L, until the melt starts becoming scarce for φ >∼0.5. Changes in slope reflect changes in nucleation rate relative to an ln-linear (in time) nucleation rate. The maximum crystal size (Lm) should increase systematically with increasing overall crystallinity.

In a batch system with a constant nucleation rate, a log-linear CSD is also possible if the effective growth rate increases exponentially with crystal size.

In a steady open system the CSD slope is determined by the product of growth rate and residence or flushing time. As long as the system is steady (i.e. ∂n/∂t = 0), nucleation rate is constant and has no influence on the slope as long as J ≠ 0. Crystal size now measures residence or growing time in the system. The age of the system is not measured by crystal size, even though the age of any particular crystal is measured by its size. If Goτ increases, the slope will decrease and vice versa. If the system is not at steady state and nucleation rate either increases or decreases with time, the slope will be affected by both Goτ and J(t).

The effects of these variations on the slope, Lmax, and intercept are shown systematically by Fig. 10. In terms of slope vs Lmax, it may be difficult to discriminate between batch and open systems, but on plots of slope vs intercept the two systems may be more differentiable (see Fig. 10). Additional information concerning the solidification history of most igneous systems is available, which may allow an easier interpretation of the CSD systematics. Before introducing further kinetic scenarios, it is essential to consider the general forms of CSDs in some natural igneous systems.

## CSD Systematics in Natural Igneous Systems

Although the majority of CSDs observed for igneous rocks are often remarkably linear on ln(n)–L plots (see Fig. 12), very few of these systems exhibit CSD systematics similar to those for either strictly batch or open systems. This is not to say that natural systems do not show systematic variations, they clearly do, and they are broadly similar for many systems. But these systematic variations are distinct from what has so far been considered, and are also distinct relative to CSDs in many industrial systems.

Fig. 12.

CSDs of plagioclase of the Peneplain sill, Antarctica. The clear linearity in these CSDs from at or near the contacts of this sill should be noted. The error bars relate to the internal measurement statistics.

Fig. 12.

CSDs of plagioclase of the Peneplain sill, Antarctica. The clear linearity in these CSDs from at or near the contacts of this sill should be noted. The error bars relate to the internal measurement statistics.

There are two features of igneous systems that make them difficult to work with as a CSD system: (1) Active systems are generally inaccessible to direct observation during crystal nucleation and dispersal. Are most systems more aptly described as local batch systems, or are fresh crystals continually sorted and deposited so that any local, say, 10 m region functions more like an open system? (2) Because lavas are never erupted carrying more than ∼55% (vol.) phenocrysts, just at the point in a strictly batch system where loss of melt begins to dampen nucleation, sampling through the natural process of eruption becomes improbable (e.g. Marsh, 1981). The magma surely goes on to complete solidification, but slow cooling with approach to the solidus allows extensive annealing and modification of the earlier CSDs. The most reasonable approach is first to gauge the applicability of the previous scaling relations using CSD results from systems where the style of the igneous process is independently known. The most notable of these systems are the Hawaiian lava lakes and diabase sills emplaced free of phenocrysts. Thus I begin below by discussing the basic features of single CSDs and follow with a discussion of sequences for successive comagmatic CSDs from a single igneous system. The overall intent is to search for connections between igneous CSDs and the CSDs of idealized batch and open systems in the hope of discerning something of the crystallization dynamics of magmas.

### Growth rate

The CSDs of lava lakes and sills are often linear (see Figs 12 and 13, and Figs 18 and 22, below) (e.g. Cashman & Marsh, 1988; Heyn et al., 1997). Because variations in growth rate with time or crystal size introduce curvature in most model CSDs, this suggests that growth rate does not vary much with crystal size or time. Much less likely is the possibility that nucleation rate changes subtly with time to offset exactly any changes in growth rate with time. Also, a review of the evidence for the magnitude of the growth rate itself suggests variations by perhaps no more than a factor of 2–5 for plagioclase in basaltic magma cooling slowly enough such that the local cooling time is long (>∼3 years) (Cashman, 1990). Moreover, other major silicate phases, such as olivine and pyroxene, seem to have overall growth rates similar to plagioclase. This is especially clear in the common equigranular texture of lava lakes and sills initially free of phenocrysts. This may reflect diffusion-controlled growth where the coupled diffusivities of the major components may be more or less similar at any temperature, as is observed experimentally (Kress & Ghiorso, 1993, 1995; Liang et al., 1996). In extrusives and thin dikes, however, Cashman, (1993) showed large variations in growth rate spanning as much as five orders of magnitude.

Fig. 13.

(Left) the CSDs of plagioclase of Makaopuhi lava lake with increasing crystallinity, depth, and time in the upper solidification front (after Cashman, 1986). (Right) the CSD calculated for a batch system like Makaopuhi, where ab = 12, as a function of increasing crystallinity.

Fig. 13.

(Left) the CSDs of plagioclase of Makaopuhi lava lake with increasing crystallinity, depth, and time in the upper solidification front (after Cashman, 1986). (Right) the CSD calculated for a batch system like Makaopuhi, where ab = 12, as a function of increasing crystallinity.

In terms of determining the actual growth rate using CSDs, the type of system (batch or open) must be known independently along with a measure of the age of the system. In a batch system with variable nucleation and constant growth rate, the largest crystals directly measure the age of the system; growth rate comes from G = Lm/t, where t is the age of the system. In an open system, the slope 1/Gτ of the CSD gives a measure of the growth rates as long as residence time τ is known. As mentioned earlier, these two determinations will give different estimates of G because the characteristic size involved in each determination (i.e. Gτ vs Lmax) differs by perhaps a factor of ∼10. Nevertheless, it is easy to convert from one to the other if the system type is known, but more important is the realization that growth rate does not seem to be highly variable in slow-cooling igneous systems. Even still, a variation by a factor of 2–5 can obviously have a significant effect on crystal size. That the size of single generation, or freshly grown, crystals does not vary much across igneous rocks of similar cooling histories also suggests relatively small variations in growth rate. It is another matter altogether for nucleation rate.

### Nucleation

The natural log of the population density [i.e. ln(n/No)] in lava lakes and sills commonly varies by ∼6–12 over the size range [i.e. ln(n/No) = 6–12], and as G does not vary strongly with crystal size, J is well determined by the CSD. Although the absolute value of J, because it depends on the magnitude of G, will be uncertain to the same degree as G, the overall pattern of variation is clearly recorded by the population density. It is thus clear that the nucleation rate varies strongly in igneous systems, and the magnitude of nucleation rate itself varies with the style of emplacement. Lavas and lava lakes have high nucleation rates, and sills and plutons have significantly lower nucleation rates (Cashman, 1990; Marsh, 1996).

The often log-linear form of igneous CSDs shows that not only is nucleation indeed a rate process [i.e. nucleation = f(t)], but also that it commonly increases effectively exponentially with time. In fact, there is little evidence in any single igneous CSD that nucleation rate is ever constant or ever decreases with time. [A possible exception to this may be in the Mount Etna lavas that show for plagioclase, olivine, and clinopyroxene unusually flat CSD spectra at large sizes (Armienti et al., 1994).] If nucleation were not a rate process, all the scalings for tc, N, and Lo derived earlier [i.e. equations (1), (2), and (3)] would be changed, as would the basic result for crystallinity [equations (9) and (11)] (Zanotto & Galhardi, 1988; Cashman, 1993). The exact functional form of J(t) does not affect the scalings, but it does influence the form of φ(t) simply through suitable choice of J(t) and G(t) with, in real systems, due concern for the appearance of more than a single solid phase.

In all of the CSDs presented henceforth, the uncertainties associated with the measurements themselves are ignored. This is not to say that the uncertainties are unimportant, for they are, but at this level of inquiry it is mainly the shapes of the CSDs and not the numbers themselves that are of interest. A concise discussion of reproducibility and minimum necessary measurement area (generally ∼1.5 cm2) has been given by Armienti et al., (1994). Methods of extracting high-quality CSD information from measurements in two dimensions have been treated by Higgins, (1994) and Peterson (1996).

## CSDs of Comagmatic Sequences

### Kinked CSDs in batch systems from growth by annexation

Crystallization within the upper solidification front of a lava lake would seem to be the ideal batch system, especially at higher crystallinities. That is, once the crystallinity in the downward advancing front reaches ∼25%, crystals are unlikely to settle from the front, or even move much relative to one another, and the system crystallizes as a collection of local batches (e.g. Mangan & Marsh, 1992; Marsh, 1996). We would expect the CSD to have a constant slope and simply migrate systematically, as in Fig. 4, to larger L and larger ln(n). CSDs for plagioclase in the upper solidification front of Makaopuhi lava lake, spanning a range of crystallinity from about 10 to 80% (vol.) are shown by Fig. 13 (Cashman & Marsh, 1988). Also shown is the CSD of an ideal batch system with constant growth rate and exponential nucleation over time, where the overall increase in nucleation (i.e. ab = 12–0) is similar to that observed.

Instead of each being parallel to one another, the observed CSDs fan out relative to the CSD of lowest crystallinity, which is from the deepest (i.e. hottest) portion of the front. Because crystallinity is a direct measure of time, it is clear that with time the largest crystals become larger faster than the smaller crystals. Crystals larger than ∼0.1 mm seem to grow at the expense of crystals smaller than ∼0.1 mm. There is a decided kink in the CSD near L ∼0.1 mm. I hasten to mention also that the CSD at the smallest sizes (i.e. 0 ≤ L ≤ 0.1 mm) is not particularly well determined as only two or three points determine the slope in this region. Nevertheless, the kink persists in all these CSDs as well as in those of broadly similar systems. The fanning out of these kinked CSDs and the fact that the earliest CSD is not kinked indicates that this kinking is not simply due to separate exponential nucleation events.

A clue to the cause of this kink in the CSD comes from studying the crystallization process itself as captured in thin sections. (I used a series of photographic transparencies of Makaopuhi thin sections kindly supplied by Dr T. L. Wright some time ago.) These same sections have been shown by McBirney, (1993). Beginning at ∼10% (vol.) crystals and continuing to the solidus at ∼95% crystals, these eight sections show a distinctive process of multiple phase heterogeneous nucleation with growth caused by both diffusion and grain boundary migration associated with annexation of smaller grains by larger ones of the same phase. More specifically, clinopyroxene, olivine, and plagioclase coprecipitate and steadily stimulate nucleation of one another through local enrichment of chemical components rejected during growth of foreign phases. Crystallization thus begins as poorly sorted clusters of all three minerals all of which coarsen with progressive crystallization. During this process larger crystals annex adjacent small crystals of the same phase through grain boundary migration. The net effect is crystal growth through both diffusion of components to the crystal surface and by annexation of smaller crystals by larger ones. The overall texture continually coarsens or ripens into an assemblage of spatially separate, relatively large optically continuous crystals set in a groundmass of smaller, later-appearing crystals. Representative photomicrographs and schematic interpretive sketches are shown as Fig. 14.

Fig. 14.

Schematic illustrations (right) of the crystallization process in tholeiitic basalt of Makaopuhi lava lake as suggested in the accompanying photomicrographs (left) of samples collected by T. L. Wright (personal communication) in 1965. The dimension of each photographed area is about 0.9 mm × 1.2 mm; the upper section contains 10% and the lower one 92% crystals. Crystallization is by local multiphase nucleation and grain boundary migration. The final result of this rapid coarsening process is a collection of optically continuous crystals made from many small crystals.

Fig. 14.

Schematic illustrations (right) of the crystallization process in tholeiitic basalt of Makaopuhi lava lake as suggested in the accompanying photomicrographs (left) of samples collected by T. L. Wright (personal communication) in 1965. The dimension of each photographed area is about 0.9 mm × 1.2 mm; the upper section contains 10% and the lower one 92% crystals. Crystallization is by local multiphase nucleation and grain boundary migration. The final result of this rapid coarsening process is a collection of optically continuous crystals made from many small crystals.

In essence, crystallization is controlled by local melt stoichiometry, with each nucleation cluster crystallizing relatively independent of neighboring clusters. The actual scale of each cluster increases with the local increase in crystal size, which continually fosters nucleation. To a large degree, this process reflects the magnitude of the difficulty of chemical diffusion relative to thermal diffusion, which is measured by the nondimensional Lewis number (Le). Le is typically ∼10−4 for basalt, which implies that chemical diffusion, relative to thermal diffusion, is a much more local process. From a standard diabase or dolerite texture, slower cooling will coarsen this texture into that of a conventional gabbro. The overall process reflects a coarsening similar to that common in dendrite crystallization experiments (e.g. Means & Park, 1994) through grain boundary migration (Hunter, 1986) and surface energy minimization (Boudreau, 1994). In effect, as the texture coarsens the kinetic time scale for further coarsening increases such that the Avrami number becomes smaller, which allows quenching of the texture as a chilled diabase or gabbro.

In an approximate fashion, this growth process by coarsening or aggregation is the reverse of fragmentation accompanied by growth of the antifragmenting crystals. Fragmentation itself is a fractal. A fractal, by definition, is a power-law relation between particle number and size. The number of objects (N) of a characteristic size greater than L is (e.g. Turcotte, 1986)

(57)
$N∼L−D$
where the exponent D is the fractal dimension. The fractal dimension for fragmentation is often 1.5 < D < 3.5. In terms of a cumulate number of crystals as a function of L, this proportionality suggests the equation
(58)
$N=1−(1+L)−D$
from which the population density is found by taking the derivative with respect to L and writing the result as a natural log:
(59)
$ln (n)=ln (D)−(D+1)ln (1+L).$
Population densities for D = 1.5, 2.5, and 3.5 are shown by Fig. 15. These CSDs show distinct upward curvature which is broadly similar to that observed, but the observed CSDs are especially linear for L > 0.1 mm. Although this process will produce larger crystals at the expense of smaller ones, there is no actual conservation of these number of particles across the CSD, nor is there any obvious way to incorporate crystal growth.

Fig. 15.

A possible CSD caused by the fractal process of antifragmentation. The constant D is the fractal dimension found for fragmentation processes. The distinct upward curvature of the CSD should be noted.

Fig. 15.

A possible CSD caused by the fractal process of antifragmentation. The constant D is the fractal dimension found for fragmentation processes. The distinct upward curvature of the CSD should be noted.

A more straightforward way to consider this process is to use the original batch CSD equation [i.e. equation (37)], which when written for size-dependent growth becomes

(60)
$∂n∂t+G∂n∂L+n∂G∂L=0.$
If small crystals are being lost from one segment of the CSD and added to another segment, terms must also be added to this equation to consider, in essence, the output and input of crystals across the CSD. Let these flux terms be given by, respectively, Kon and Kin. Then (60) becomes
(61)
$∂n∂t+G∂n∂L+n∂G∂L(Kon−Kin)=0$
where Ko and Ki are essentially rate constants.

This equation can be put in a more convenient form by dividing throughout by nG and rewriting the derivatives in logarithmic form, which upon rearranging becomes

(62)
$∂ln(n)∂L=−1G∂ln(n)∂t−n∂ln(G)∂L−1G(Ko−Ki)$
where Ko and Ki are each treated as intrinsically positive entities. This form is convenient because it shows the contributions to the CSD slope, which is the term on the left, in terms of changes in nucleation rate (second term), in growth rate as a function of crystal size (third term) and the loss (Ko) or gain (Ki) of crystals from one segment of the CSD. It should be noted that each of these terms, except for Ki, gives rise to a negative CSD slope, and also that, if G were to increase with L, the CSD slope would lessen or become concave up as in the case of antifragmentation. This equation can now be applied to understanding these kinked CSDs.

#### CSD for small sizes

Let us consider the Makaopuhi CSDs where, to a good first approximation, G is constant (i.e. no observed curvature in the CSDs), and for the small size range (0 < L < 0.1 mm) Ki = 0. The above equation then becomes

(63)
$∂ln(n1)∂L=−1Go∂ln(n1)∂t−KoGo$
where the notation n1 is remindful of the smaller size range.

Because G is constant the second term can be written in terms of nucleation rate [see equation (48), etc.], yielding

(64)
$∂ln(n1)∂L=−1Go∂ln(J1)∂t−KoGo.$
For an exponential rise in nucleation with time, as described by equation (7) [see also equation (52), etc.], this becomes
(65)
$∂ln(n1)∂=−aLm−KoGo.$
The first term on the right gives the usual batch CSD slope and the loss of crystals at a rate Ko further steepens the slope (see Fig. 16). The actual CSD is found by integrating (65) with respect to L:
(66)
$ln (n1)=−(aLm+KoGo)L+C1.$
C1 is a constant found from the condition that n(t, L = 0) = no(t), and at any specific time or crystallinity, the CSD is given by
(67)
$ln (n1/n1o)=−(aLm+KoG1)L$
which is reflected in the slopes of the CSDs of Fig. 16.

Fig. 16.

(Left) kinking in the CSD as a result of aggregation of small crystals to form larger crystals. (Right) the possible migration of a kinked CSD in a batch system.

Fig. 16.

(Left) kinking in the CSD as a result of aggregation of small crystals to form larger crystals. (Right) the possible migration of a kinked CSD in a batch system.

#### CSD for large sizes

In the large size range (i.e. L > 0.1 mm), the slope is similarly modified, but in the opposite sense of that for the small size range. Thus if Ko represents the outflux from the small sizes, then the input flux to the large sizes must be equal and opposite. The difference in signs of Ki and Ko in (62) already accounts for this opposite sense, so in magnitude Ki = Ko. The CSD for the large sizes is just (67) with a sign change, namely,

(68)
$ln (n2/n2o)=−(aLm+KoG2)L.$
Here the slope of the CSD is less steep than that of both the small size range and the ideal batch system (see Fig. 16).

It should also be emphasized that the growth rate Go will be different for each size range. In the small sizes, growth rate should be that associated with the normal sense of diffusive growth (i.e. Gi = Go). But in the large size range it will be enhanced by the aggregation effect such that G2 = Go + G′, where G′ is the growth caused by aggregation.

The effective growth rate can often be found from a linear CSD by assuming, in analogy with open systems (and see below), that the slope measures Goτo, where τo is an effective residence time. If τo is known, Go can be found from the observed slope. For the present batch system the CSD slope is given by a/Lm or, because Lm = Gotc, the slope is a/Gotc, where now tc/a is an effective residence time. In this fashion for the kinked CSD the effective residence times are, for the small size range

(69)
$τ1=tc/aKo[(1/Ko)−(tc/a)]=τo/Ko(1/Ko)−τo$
and for the large size range
(70)
$τ2=tcG2/GoaKi[(G2/GoKi)−(tc/a)]=τoG2/GoKi[(G2/GoKi)+τo]$
where τo is the effective residence time of the ideal batch system. Although τ1 and τ2 may be independently known, it is clearly challenging to extract more than simply effective measures of G and K from kinked CSDs because of aggregation.

The net effect of aggregation annealing on the batch CSD is also shown by Fig. 16, where the kinked CSD migrates across the diagram. That the calculated CSDs do not fan out at large sizes, as observed, suggests that the flux rate constants (K) are themselves changing with time. The size range supplying the small crystals for aggregation may also be increasing slightly with time. Each of these effects will cause fanning in the CSDs at larger sizes. It is also not at all certain that the entire range of small crystals (i.e. 0 < L < 0.1 mm) is simultaneously supplying crystals to the larger size range. In fact, from the earlier description of the grain annexation within the clusters, it would seem more reasonable that, once nucleated, growth may proceed normally until a certain critical size whereupon the crystal meets a larger crystal and is incorporated. If so, then the formulation leading to (62) would need to be amended to allow output only over a restricted size range. But the same analysis can be used in a piecewise sense across the CSD. Some indication of such a process may be in the detailed bulk CSD (i.e. all solid phases in one CSD) of the 1984 lava from Mauna Loa (Crisp et al., 1994). Here there is an almost discontinuous drop in the CSD at L ∼ 20 µm, with both the earlier (i.e. L > 20 µm) and later CSD slopes being much gentler. This sudden change could also, however, be simply associated with a strong nucleation event in response to the decompression accompanying eruption, which is favored by Crisp et al.

### Batch vs open system aggregation

Allowing for the flux of crystals from one part of the CSD to another [i.e. Ko and Ki in (62)] makes the batch system somewhat analogous to an open system where the crystal flux is due to a global removal of crystals from the system itself. Beginning with an exponential increase in nucleation, a typical log-linear CSD can be produced in each instance from a steady state of nucleation where ∂n/∂t = 0. That is, from (64),

(71)
$∂ln()∂L=−KoGo$
will hold for the small size range. And in the true open system (i.e. the MSMPR system) at steady state, from (54)
(72)
$∂ln(n)∂L=−1Goτ.$
The rate constant Ko is analogous to the reciprocal of the residence time in the open system. The two CSDs will be similar, but each will reflect a fundamentally different process. The distinct difference in process will be reflected in the full CSD, which will migrate in the batch system and be steady in the open system. Aggregation also produces a kinked CSD in both systems.

One such open system has been produced experimentally by Burkhart et al. [1980; see also Randolph & Larson, (1988), p. 294]. In this experiment, crystallites of ammonium polyuranate aggregate at very small, unmeasurable, sizes to form crystal clusters, which themselves further aggregate to form agglomerates. The method of analysis is as used above for the batch system except that the system is at steady state. The observed CSD is shown as Fig. 17, which resembles the Makaopuhi CSDs, but this is an open system, which Makaopuhi is not. The diagnostic difference is that this CSD is invariant with time, whereas the Makaopuhi CSDs systematically change with time for increasing crystallinity. The general rule in most volcanic systems is that the CSDs do change, albeit not necessarily systematically, with time. A remarkable exception to this general rule is the CSDs of plagioclase in the lavas of Mount Etna presented by Armienti et al., (1994). These CSDs have not changed in any noticeable detail for a period of at least 70 years.

Fig. 17.

CSD of an experiment involving aggregation of clusters of ammonium polyuranate (after Burkhart et al., 1980).

Fig. 17.

CSD of an experiment involving aggregation of clusters of ammonium polyuranate (after Burkhart et al., 1980).

### Kinked CSDs in other systems

#### Atka comagmatic sequence

Phenocryst-rich lavas often show kinked CSDs. A clear example of this is the CSDs measured by Resmini, (1993) of high-alumina basalts from Atka volcanic center, Aleutian Islands, Alaska. The CSDs of plagioclase are of 12 flows in a sequence of 22 flows from the northeast region of northern Atka. The modal mineralogy of these flows was reported by Marsh, (1981) and the general geochemistry was reported by Myers et al., (1986). Anorthitic plagioclase is by far the dominant phenocryst, typically making up ∼95% of the overall phenocryst content, which itself ranges between about 35 and 48%. Some representative CSDs for these lavas are shown by Fig. 18.

Fig. 18.

CSDs of plagioclase of high-alumina basalt from Atka, Aleutian Islands, Alaska (from Resmini, 1993).

Fig. 18.

CSDs of plagioclase of high-alumina basalt from Atka, Aleutian Islands, Alaska (from Resmini, 1993).

These CSDs are characterized by a distinct break in slope or kink, marking a small size range (L < ∼0.55 mm) and a larger size range (0.5 mm < L < ∼2.0 mm). Each size range is well represented by a straight-line fit, which defines the whole CSD by two slopes and two intercepts. In addition to these slopes and intercepts, the overall plagioclase content (vol. %) and the maximum crystal size are known for each of these basalts. The maximum crystal size is estimated from the average of the four largest L values (i.e. Lmax) of each CSD; in cases where there is a large gap between these points, the most extreme values of L were neglected.

As the maximum crystal size is a direct measure of system age or residence time, it is of interest to consider, as mentioned earlier, correlations between slope and intercept and Lmax. Plots of slope against Lmax, intercept against Lmax, and slope against intercept are shown in Fig. 19a. As might be expected, the slopes (S2) and intercepts (n2o) for the larger size range of each CSD correlate better with Lmax than do the values of S1 or n1o. And the larger sizes also correlate better on the slope-intercept plot. Referring back to the diagnostic plots of Fig. 10, these Aleutian lavas may have characteristics of a batch system with a time-decreasing nucleation exponential constant [i.e. a(t) in (7) decreases with time]. This does not necessarily rule out the possibility of an open system, for a decreasing nucleation rate in an open system would also give similar diagnostics. Overall, however, it is clear that the effective nucleation rate (i.e. found by extrapolation from the large size range to the L = 0 limit) decreases as crystal size increases. The associated CSDs thus fan upward about a pivot point near L ∼ 0.75 mm. With this fanning upward and lessening of slope not only does Lmax increase but so does the overall crystallinity.

Fig. 19.

(a) CSD diagnostic relations among slope, intercept, and maximum crystal size for the comagmatic sequence of Atka lavas. (b) (Upper) correlations between maximum crystal size and modal plagioclase content in Atka lavas. (Lower) detailed variation between slope − 2 and intercept. The larger number against each symbol is the sample number, which indicates sequence of eruption, and the smaller number is the maximum crystal size (mm). There is a clear correlation between decreasing slope, decreasing intercept, and increasing crystal size.

Fig. 19.

(a) CSD diagnostic relations among slope, intercept, and maximum crystal size for the comagmatic sequence of Atka lavas. (b) (Upper) correlations between maximum crystal size and modal plagioclase content in Atka lavas. (Lower) detailed variation between slope − 2 and intercept. The larger number against each symbol is the sample number, which indicates sequence of eruption, and the smaller number is the maximum crystal size (mm). There is a clear correlation between decreasing slope, decreasing intercept, and increasing crystal size.

The variation in maximum crystal size with plagioclase content is shown by Fig. 19b. Within this variation, there may also be a more subtle subdivision into a larger Lmax series and a smaller Lmax series. That is, for any crystallinity the maximum crystal size follows one of two variations, which are similar. These two series, for the most part, also distinguish themselves on the accompanying, more expanded slope-intercept plot for the larger size range CSDs. Both slope and intercept lessen with increasing crystal size and overall crystallinity, but the eruptive sequence (smaller to larger sample numbers) shows no correlation with slope or intercept for the S2n2o CSDs, and for the S1n1o CSDs only a weak negative correlation with slope and a weak positive correlation with intercept is evident. This may be reasonable as originally deeper-held magma is increasingly decompressed during the eruption. There is, however, no correlation between eruptive sequence and phenocryst content.

Comagmatic lava sequences so far examined for Makaopuhi lava lake, Mount St Helens dacites, Dome Mountain, Nevada, Atka, and Kilauea show a decrease in CSD slope and intercept with increasing crystal size. This may be a general characteristic of local batch systems where both chemical (i.e. annealing through aggregation) and physical (i.e. sorting by buoyancy) processes tend to produce populations of large crystals, which indicates a low effective nucleation rate. What is also especially clear in these Atka lavas is that the CSD varies systematically spatially in the pre-eruptive staging area. Exceptions to this pattern are the Mount Etna and Kameni flows reported by, respectively, Armienti et al., (1994) and Higgins, (1996). These systems seem to show a much more steady nucleation rate and a more time invariant CSD pattern in general. These systems may be more akin to a true steady-state, open system.

If the Atka lavas are thus interpreted as a batch system, considering the results of the earlier modeling leading to Figs 5, 7, and 8, each linear or quasi-linear range of the CSDs can be represented by a value of the nucleation exponent a and a value of Lmax for that segment. For the larger size range, a varies between about 4 and 7, with the larger values (i.e. steeper slopes) characteristic of the smaller Lmax series of Fig. 19b. The values of Lmax range between about 1.5 and 2.0 mm, which relates to an effective crystal radius of 0.75–1.0 mm. From the results shown by Fig. 5, for a constant growth rate (b = 0), 4 ≤ a ≤ 7, and Lmax ∼ 1 mm, suggests that JoGo × 104. Assuming Go ∼ 10−9cm/s, gives a nucleation rate of ∼10−5s−1, which means that this size range could mark an event of perhaps 25 years or less in duration.

For the smaller size range, 3 ≤ a ≤ 5 and Lmax ∼ 0.5 mm, which translates to a marked increase in nucleation rate. That is, for an equivalent change in crystal size, this latter nucleation rate is larger by a factor of ∼100 over that of the earlier nucleation event. This suggests a change in crystallization environment, which is most probably related to migration to the near subsurface pending eruption with increased nucleation as a result of high-level volatile loss. These crystals mark an event of perhaps a year or so in characteristic time.

In sum, the characteristic kinking of the Atka CSDs suggests a quasi batch system marked by two distinct nucleation event environments. Within each environment the systematic fanning of the CSDs, especially in the larger size range, suggests a systematic spatial variation in crystal population density within the magmatic body. The systematic correlations between plagioclase size and modal content and CSD characteristics further suggests that these CSDs, as a whole, are intimately linked. How they might be linked can only be judged in relation to the dynamics of the actual system.

In this context, one of the obvious difficulties in interpreting the CSDs from volcanic sequences concerns the basic physical CSD nature of the magmatic system itself. In this regard, it is insightful to consider CSDs from dynamically well-characterized systems, which is only possible in industrial crystallization studies.

#### Industrial ‘comagmatic’ systems

Mixed suspension, mixed product removal (MSMPR) magma crystallizers are exceedingly common in industry (e.g. Jancic & de Jong, 1984). In the idealized system, a continually well-mixed tank of volume V of liquid and growing crystals is fed by a flow of crystal-free liquid and emptied at the same rate (Q) by an outflow containing a true representation of the overall contents (liquid and crystals) of the tank. The CSD of the outflow at steady state is log-linear and time invariant, as given earlier by (56) and as discussed in detail by Marsh, (1988) and Randolph & Larson, (1988). The time to become steady is observed to be (12–15)τ, where τ is the residence time measured by V/Q. As mentioned already, although linear CSDs are not uncommon in igneous systems, there is little sign that they are invariant with time and that they represent a true MSMPR system. Because deviations from the log-linear classic MSMPR CSD are also common in industrial applications, it is of interest to consider in some detail the relationship between observed non-ideal MSMPR CSDs and crystal sorting within the crystallizer itself. After considering the industrial results, I will briefly turn to translating the tank crystallizer to a possible equivalent magmatic system.

A standard industrial draft-tube baffle (DTB) crystallizer and the associated CSD for potassium chloride as measured by extractions from three dynamic zones of crystal sorting are shown by Fig. 20 (Canning, 1970). A propeller at the base of the draft tube in the center of the tank keeps the fluid (i.e. crystals plus liquid) circulating. Although the incoming fluid is crystal free, nucleation occurs freely in response to supersaturation caused by evaporation, cooling, or earlier crystal growth. (Here growth is due to evaporation.) The entire central volume or active zone (Zone 2) is ideally kept well mixed and the product extracted from this region shows a linear CSD. The active zone is surrounded by a baffle settling region (Zone 1) that accepts overflow, which consists of fine particles, and allows settling. Fluid drawn from the top of this zone contains the smallest size range of crystals. This fluid is sometimes reheated or otherwise treated to dissolve these ‘fines’ and is fed back into the tank. Crystals too large to be kept in circulation settle to the base of the tank (Zone 3) where a crystal-rich slurry is removed in the product flow. These largest crystals have a lower residence time relative to smaller sizes and the CSD is kinked downward in the active zone. Overall, however, the large size of these crystals betrays their age.

Fig. 20.

The operation (left) and resulting CSD (right) of a draft-tube baffle, mixed suspension, mixed particle removal (MSMPR) industrial crystallizer (after Canning, 1970).

Fig. 20.

The operation (left) and resulting CSD (right) of a draft-tube baffle, mixed suspension, mixed particle removal (MSMPR) industrial crystallizer (after Canning, 1970).

The intent of the MSMPR system is to keep the active, central zone well mixed such that the outflow is an accurate representation of the internal dynamics, but each separate zone has a CSD reflecting local crystal fractionation (Fig. 20). In Zone 1 the smallest crystals are drawn off relative to all others, this fluid has the smallest residence time in the system, and the slope of the CSD is steep. This reflects both strong fractionation and a small residence time. The assemblage of Zone 2, the active region, has lost the smallest crystals, which gives it a low effective nucleation density, but a long residence time because of stirring. If the forced convection is not strong enough, however, to keep the largest crystals in suspension before extraction in the slurry line, they will settle to the base and be removed in the product slurry. The CSD in this region will show some fractionation over that of the active zone; the slope will be slightly steeper because of premature crystal loss. The net result is a double kinked CSD (Fig. 20) reflecting significant crystal fractionation. It must be kept clear, however, that the observed CSD comes from sampling of specific parts of the system. The bulk of the system is represented by Zone 2, which alone would give a log-linear CSD. Given the composite CSD alone, different conclusions might be drawn on the dynamic regime within the crystallizer. But it is none the less clear that the CSD of each zone coupled with the maximum crystal size of each zone is crucial in discerning the dynamics of each zone. It is this interpretive inversion that lies at the heart of understanding igneous systems.

#### Composite CSDs and mixing

The two fundamental differences between igneous and industrial systems are that (1) the product crystals generally remain somewhere in the system and (2) sampling of the system is by eruption, which is generally not deliberate. Moreover, residence time is not independently known and the system cannot be freely manipulated to yield kinetic information. Thus in the equivalent magmatic problem, crystals come from two principal sources, both of which are mainly due to cooling. One source is crystals nucleated and grown at the present location of the magma. These crystals are almost always associated with solidification fronts and are available for sorting in the magma only up to crystallinities of ∼25% (e.g. Marsh, 1996). These crystals will always be small. The second source of crystals is cumulate beds of previously grown crystals incorporated by the magma during transport and emplacement. These are the so-called ‘tramp’ crystals [e.g. forsteritic olivine in Hawaiian picrites (Wright & Fiske, 1971; Marsh, 1996)], which may have resided for some time in the feeding conduit system and annealed or ripened (see below) to a strongly non-log-linear CSD. And where the system is multiply saturated, aggregation-dominated nucleation and growth will produce a kinked CSD. Depending on the nature and vigor of the eruptive or displacement event, the CSD can take on any number of styles.

First and foremost, because the eruption event involves significant decompression and possibly degassing, a burst of nucleation will probably be the final event recorded by the CSD. Earlier transport from deep levels, interspersed with significant holding times, followed by wholesale eruption may produce a CSD not too dissimilar to that of Fig. 21a. Although locally the system may function as an open or quasi-MSMPR system, overall the system is a traveling batch system. Perhaps only within the solidification fronts themselves does the system truly approach a batch system. These concepts are summarized by Fig. 21a.

Fig. 21.

(a) CSD evolution associated with transmission of magma through a well-developed and long-standing magmatic ‘mush’ column. The CSD of the eruptive product may reflect a long series of local processes. (b) Mount Etna CSDs for (A) plagioclase from three samples of the 1992 eruption (January, March and May) and for (B) plagioclase, olivine and titanomagnetite of lavas of the 1991–1993 eruption. All are from Armienti et al., (1994). The overall shape of the CSDs suggests perhaps three crystallization regimes, possibly similar to those depicted schematically by Fig. 20. The 1991–1993 CSDs have each been fitted (by eye) with three straight lines of exactly the same slope for plagioclase, olivine, and titanomagnetite (two lines only). That the lines of the same slope apparently fit the CSDs for all three minerals suggests crystallization in similar regimes.

Fig. 21.

(a) CSD evolution associated with transmission of magma through a well-developed and long-standing magmatic ‘mush’ column. The CSD of the eruptive product may reflect a long series of local processes. (b) Mount Etna CSDs for (A) plagioclase from three samples of the 1992 eruption (January, March and May) and for (B) plagioclase, olivine and titanomagnetite of lavas of the 1991–1993 eruption. All are from Armienti et al., (1994). The overall shape of the CSDs suggests perhaps three crystallization regimes, possibly similar to those depicted schematically by Fig. 20. The 1991–1993 CSDs have each been fitted (by eye) with three straight lines of exactly the same slope for plagioclase, olivine, and titanomagnetite (two lines only). That the lines of the same slope apparently fit the CSDs for all three minerals suggests crystallization in similar regimes.

The simple addition of some types of cumulate crystals to a less crystallized magma of the same initial bulk composition may be directly detectable through geochemistry. The addition of plagioclase, for example, may produce a positive europium anomaly when the original bulk fluid had none. This independent information may be useful in recognizing a CSD produced from two distinct nucleation events over that from cumulate entrainment. In this regard, the CSD resulting from mixing of two magmas with distinct initial CSDs has been modeled by Higgins, (1996) for dacites from Thera. The log-linear CSD dictates that each initial CSD dominates the slope of the new CSD in the region where it has the highest population density regardless of the mixing ratios. As long as the two initial CSDs have significantly different slopes, with the larger size range having a more gentle slope, the final CSD will be kinked. An example from Higgins, (1996) is shown as Fig. 22. Moreover, if the added phenocrysts or megacrysts are for the most part unzoned plagioclase and the host magmas are of similar composition, they reflect either ripening in a holding region or a long growth period under steady, low nucleation density, and single saturation. The fact that both olivine and plagioclase (and perhaps also titanomagnetite) show exceedingly similar kinked or flattened CSD spectra in the Mount Etna lavas (see Fig. 21b) strongly suggests deep-seated storage as a common feature of this magmatic column (Armienti et al., 1994).

Fig. 22.

(a) Modification of the CSD as a result of mixing of magmas with distinct CSDs, and (b) a possible example from Thera (both after Higgins, 1996). (c) A possible example of a CSD arising from magma mixing, Usu volcano, Japan (after Tomiya & Takahashi, 1995).

Fig. 22.

(a) Modification of the CSD as a result of mixing of magmas with distinct CSDs, and (b) a possible example from Thera (both after Higgins, 1996). (c) A possible example of a CSD arising from magma mixing, Usu volcano, Japan (after Tomiya & Takahashi, 1995).

Another example of eviction of magmas carrying related CSD spectra, but where the parental magma types vary considerably in bulk composition, may be in the post-1663 products of Usu volcano, Japan (Tomiya & Takahashi, 1995). Three phenocryst types (A, B, C) of, respectively, low, high, and intermediate An (plagioclase) and mg-number (opx) composition are found in addition to microclots of mafic magma within rhyolitic magma. The A (∼An45) and B (∼An90) plagioclase types are large (up to 2.5 mm) and remarkably unzoned, and the plagioclases of the microclots are smaller and profuse. The presence of A and B phenocrysts in the 1663 lavas and mainly only C (∼An70) phenocrysts in the 1769 lavas, suggests early commingling with eruption followed by re-equilibration in the remaining magma with eventual eruption. These CSDs are also shown by Fig. 22. In comparison with the mixing diagram of Higgins, (1996), to match the CSDs in the small size range through mixing, the mafic magma may have contributed >50% (vol.) to the later (e.g. 1769) magma. The overgrowths on the rare A plagioclase in the 1769 dacite are generally 50–70 µm, which, if developed over a 100 year period, amounts to a growth rate of about 2 × 10−12 cm/s; this is similar to the rate found for plagioclase in the 1980–1986 Mount St Helens dacite (Cashman, 1988). Yet another example is the olivine-bearing tholeiites of Kauai, Hawaii, discussed by Maaloe et al. (1989).

In summary, industrial MSMPR systems show a CSD spectrum broadly similar to what can be expected in a generalized magmatic system. The constituent crystals are sorted by dynamic processes that in the industrial system can be individually tapped to reveal the influence of the local dynamics on the CSD. Magmatic systems are unlikely to be systematically tapped to reveal in the flow sequence a consistently recognizable CSD ensemble, but that does not mean that judicious study will not reveal the intimate kinetic details of the system.

Magmatic systems are, strictly speaking, mainly closed systems that emit a wide range of possible CSD spectra, even within a comagmatic suite. Individual suites, however, generally show a consistent CSD style that is often reflected in the detailed petrology (i.e. overall crystallinity, crystal zoning, and geochemistry) and is understandable within the framework of a generalized magmatic system consisting of serial holding regions connected by transport zones all dominated by solidification fronts and new and old phenocrysts.

### Linear CSDs in comagmatic systems

A series of 20 basalts and andesites emitted in rapid succession in the moat of the Timber Mountain caldera show consistently linear CSDs. The overall nature of the CSDs and the general petrology of the lavas have been described by Resmini & Marsh, (1995). The sequence begins with several basalts followed by a long series of andesites all of which are of fairly low crystallinity (<15%), which is almost solely plagioclase. The diagnostic indicators of these CSDs are shown by Fig. 23.

Fig. 23.

CSD diagnostic relations among slope, intercept, and maximum crystal size in lavas from Dome Mountain, Nevada. The numbers near the data points are flow numbers.

Fig. 23.

CSD diagnostic relations among slope, intercept, and maximum crystal size in lavas from Dome Mountain, Nevada. The numbers near the data points are flow numbers.

A plot of slope against intercept (Fig. 23) does show overall a typical negative slope, but it also shows a sequence (perhaps two) of flows that has a more or less constant slope and systematically increasing intercepts. In the wholesale eruption of a batch system, a sequence of CSDs, amongst others, might be expected to show constant slope and an intercept increasing with crystallinity. In half of these 20 flows there seems to be such a relation among the andesites (see Fig. 23). And within this overall constant slope sequence there is perhaps a hint of two related sequences of increasing slope with increasing flow number (numbers increase with order of eruption). The maximum crystal size is also expected to increase with increasing intercept in a batch system. Although Lmax does increase overall for this group of flows, within each sub-sequence (as delineated on the intercept-crystallinity plot) crystal size remains approximately constant at about, respectively, 1.0 and 1.5 mm (Fig. 23). On the plot of slope against maximum crystal size there is an increase in Lmax as slope remains approximately constant (Fig. 23), although the suggestion of a correlation with flow number is less certain here. The earlier flows (1–5) show more scatter in slope and crystal size is generally large. Those flows of similar slope generally group together (numbers 6–10 and 16–20).

Overall, these Dome Mountain lavas show characteristics of strong unsteadiness in the earliest lavas followed by two sequences of CSDs possibly indicating batch system evolution, where the slope remains fairly constant as crystal size and intercept both increase. There is some scatter in the slopes but the longest series of ten lavas has a slope of about −12 ± 1, which is fairly tight considering that the entire sequence shows slopes from about −5 to −17.

#### Peneplain sill, Antarctica

This 330 m thick diabase occurs within a sequence of thick sills in the Dry Valley region of Antarctica. The general nature of these sills has been described by Gunn, (1966). The Peneplain sill is particularly noteworthy in its occurrence at Solitary Rocks in being exceedingly uniform in texture and composition throughout. Only near the top is the uniformity disrupted by a series of silicic segregations caused by internal tearing of the roofward solidification front (Marsh & Wheelock, 1994; Marsh, 1996). Because the sill here carried few, if any phenocrysts, it is particularly valuable from which to gauge the effects of in situ or batch crystallization. The rocks are all holocrystalline and the plagioclase content of the diabase is more or less constant from top to bottom of the sill.

The CSDs of plagioclase in the Peneplain sill are in general linear (see Fig. 12; some show some relatively minor curvature) throughout the body (Heyn et al., 1997). There is little sign of any dramatic turn-down in the smallest size range of the CSD. On a slope-intercept plot (not shown) most of the slopes fall between −4 and −6, which is a tight batch-like range relative to that seen in the Dome Mountain lavas. Almost all these samples are from the sill interior, whereas the samples near the upper and lower contacts have steeper slopes and larger intercepts. The reason for this grouping is apparent if the CSD data are viewed against stratigraphic height within the sill itself (see Fig. 24). By and large, the slopes, intercepts, and maximum crystal sizes (Lmax) show distinctive variations inward from each contact. Slope and crystal size each increase inward whereas intercept decreases until near the sill center where each trend abruptly reverses. Each parameter returns to near contact values at the sill center. A possible reason for this variation has been suggested by Heyn et al., (1997).

Fig. 24.

The variation of the CSD parameters slope, intercept, and maximum crystal size with height through the Peneplain sill of diabase, Dry Valleys, Antarctica. The diagram on the lower left shows the variation in nondimensional isotherm velocity (VL/K) as a function of position within a sheet of thickness L cooling by conduction; K is thermal diffusivity.

Fig. 24.

The variation of the CSD parameters slope, intercept, and maximum crystal size with height through the Peneplain sill of diabase, Dry Valleys, Antarctica. The diagram on the lower left shows the variation in nondimensional isotherm velocity (VL/K) as a function of position within a sheet of thickness L cooling by conduction; K is thermal diffusivity.

It was apparently first noted by Tomkeieff, (1940) [see also Jaeger, (1968)] that the velocity of isotherms within a cooling sheet is large both at the contacts and also at the center. The actual variation in isotherm velocity is also shown by Fig. 24, where the numbers on the curves give the nondimensional temperature (Heyn et al., 1997). It is obvious why isotherms move fast near the contacts where cooling is rapid, but much less obvious why this occurs also at the center. In essence, as the upper and lower cooling fronts approach each other heat is taken from an increasingly smaller volume at the sill center. Isotherm velocity thus increases without limit as the exact center is reached. This effect has been suggested by Tomkeieff, (1940) and Jaeger, (1968) to be instrumental in the variation of joint patterns in lava flows.

In the Peneplain sill rapid isotherm advancement in the center produces CSDs that resemble those at the contacts. The slowest advance of isotherms is at intermediate horizons between the contacts and the center, for it is here that the crystals are largest, the slope is smallest, and the intercept lowest. Inward from the contacts, nucleation rate decreases and then increases near the center. Nucleation rate is directly related to isotherm velocity or cooling rate, which is not surprising (e.g. Cashman, 1993), but this nucleation rate variation alone, with more or less constant growth rate, can give rise to all these CSD variations. On the other hand, if this correlation should prove to be robust and general, CSDs from volcanic suites such as that of Dome Mountain may give direct insight into subsurface cooling regimes through the nucleation rate (Resmini & Marsh, 1995).

Similar evidence for rapid final cooling at the center of a sheet comes from prehistoric Makaopuhi lava lake (Evans & Moore, 1967). Here the glass content increases more or less symmetrically toward the center, from a few percent near the top and bottom to >15% near the center. There is also a hint of this effect in plagioclase size and composition, which was originally suggested to be due to foundered crust. Additional hints of this effect, although slight, may be apparent in the CSD of Antarctic lavas and sills measured by Wilhelm & Worner, (1996; see their fig. 8).

In closing this section, it should also be pointed out that the correlations with height shown by Fig. 24 are not in all respects perfect. Two samples at 470 feet and 950 feet do not have particularly correlative parameters. Whether this is meaningful or not remains to be decided. Overall, however, the style of crystallization in this sill seems to be clearly of the batch type and strongly related to the rate of cooling.

### Annealing vs attenuation of nucleation in batch systems

There is little sign in most igneous CSDs of the expected strong diminution in effective nucleation rate at the most advanced stages of crystallization. Few igneous CSDs turn down strongly at the smallest crystal sizes. This probably means that most igneous systems experience arrested crystallization because of relatively rapid cooling rates as the solidus is passed through. Metamorphic rocks, on the other hand, commonly show CSDs with maxima at small, to intermediate crystal sizes (e.g. Cashman & Ferry, 1988). Whether this is due to annealing or attenuated nucleation is difficult to assess in metamorphic systems as long periods at high, but subsolidus temperatures are common. Some evidence capturing the transition from straight to ‘humped’ CSDs may be available, however, in igneous systems.

#### Chromite in Stillwater

Relatively thin monomineralic layers in large igneous bodies may experience conditions similar to a metamorphic environment; high, subsolidus temperatures may prevail for long (i.e. many times individual crystal growth times) periods of time. CSDs of chromite in chromite cumulates of the Stillwater intrusion show a distinct maximum at relatively large (∼0.1 mm) crystal sizes (Waters & Boudreau, 1996). One of these CSDs is shown by Fig. 25, where it has been fitted with a batch system CSD model with constant growth rate (i.e. b = 0) and a nucleation rate exponential constant of a = 8 (i.e. ab = 8). The fit is particularly good for 99.5% crystallization, at which point further crystallization would have to be arrested to preserve the original CSD from further annealing under subsolidus conditions. Because the CSD cannot be observed at previous times or stages of crystallization, it is difficult to tell whether the final CSD is due to nucleation attenuation or to subsolidus annealing. If the CSD evolves upward with time (i.e. increasing nucleation density) then the system is truly a batch. But if the CSD evolves downward with time, with strong loss of the smallest crystals, this is annealing. Waters & Boudreau, (1996) found, however, clear petrographic evidence of grain coalescence or annexation and also a spatial variation of the CSD. This evidence points to this CSD being a result of annealing and not simply cessation of nucleation. Evidence on this issue also comes from sills.

Fig. 25.

CSDs of chromite of the Stillwater layered intrusion [left, after Waters & Boudreau, (1996)] and of clinopyroxene of Box Elder laccolith, Montana (Congdon et al., 1993). Although the chromite CSD could be interpreted as a result of melt loss, the more reasonable explanation is annealing, as grain annexation can be seen in thin section.

Fig. 25.

CSDs of chromite of the Stillwater layered intrusion [left, after Waters & Boudreau, (1996)] and of clinopyroxene of Box Elder laccolith, Montana (Congdon et al., 1993). Although the chromite CSD could be interpreted as a result of melt loss, the more reasonable explanation is annealing, as grain annexation can be seen in thin section.

#### Cpx in Box Elder laccolith

Sills formed of crystal-laden magma may furnish good evidence on the prior evolution of humped CSDs. That is, a comparison of the CSD of the chilled margin with that of the interior will indicate whether the CSD is evolving upward or downward with time. A CSD of clinopyroxene in Box Elder laccolith of north-central Montana is also shown by Fig. 25 (Congdon et al., 1993). This 140 m thick body overall has the composition of shonkinite (Pecora, 1941) and has a structure broadly similar to Shonkin Sag laccolith, which is some 100 km south (e.g. Marsh et al., 1991). The CSD from the chilled margin is steeper, much more linear at the small sizes, and has a larger nucleation density than that of the interior. The interior CSD exhibits a maximum at L/Lm ∼0.13 (L ∼0.2 mm), which is a sufficiently large size not to be confused with poor measuring resolution. At the same time, the CSD at larger sizes is higher than that of the chilled margin. It would appear that small crystals have been lost to the growth of larger crystals. Based on the size of the body, and the fact that it contained ∼25% phenocrysts upon emplacement, the interior would reach the solidus ∼80 years after emplacement. The interior sample, however, is ∼43 m above the lower contact and would have quenched within ∼20 years of emplacement. All in all, it appears that the maximum in this CSD is also due to post-emplacement annealing, which has occurred relatively rapidly.

Be this as it may, some caution is in order, because of the possibility that flow differentiation of phenocryst-laden magma may have also influenced the differences between these two CSDs. Flow differentiation concentrates larger crystals in the interior, away from the contacts. However, there is little reason to expect that the smaller crystals will be strongly depleted in the interior. On the contrary, they should be more evenly distributed. Thus although flow differentiation may have partly influenced these CSDs, the maximum is most probably due to annealing.

#### Resorption vs annealing in cpx of Shonkin Sag laccolith

This 3 km × 70 m potassic basalt laccolith in north–central Montana has long been studied as a small body showing a classic differentiation sequence: a thick pile of cumulates on the floor, a thinner zone of basaltic rock near the roof, and a syenite sandwich horizon in between. In fact, many of the basic principles thought to govern magmatic evolution in closed systems germinated through study of this body early in this century (Marsh, 1996). The geochemical nature of Shonkin Sag, however, is almost entirely due to the fact that the magma contained upon emplacement ∼35% phenocrysts of mainly clinopyroxene (Hurlbut, 1939). Settling of these crystals dominated the evolution of this body, forming a classic S-shaped profile of modal cpx. The general petrology and CSDs for cpx at seven horizons in the body were presented by Marsh et al., (1991) and these are shown here as Fig. 26a along with a stratigraphic profile. These CSDs are particularly revealing in terms of crystal settling and cpx resorption in the syenite horizon.

Fig. 26.

(a) CSDs and stratigraphy with bulk MgO (wt %) variation of Shonkin Sag laccolith, Montana (after Marsh et al., 1991). (b) CSD diagnostic relations among slope, intercept, and maximum crystal size for Shonkin Sag laccolith. The pronounced effects of resorption should be noted.

Fig. 26.

(a) CSDs and stratigraphy with bulk MgO (wt %) variation of Shonkin Sag laccolith, Montana (after Marsh et al., 1991). (b) CSD diagnostic relations among slope, intercept, and maximum crystal size for Shonkin Sag laccolith. The pronounced effects of resorption should be noted.

The upper and lower chilled margins show a phenocryst-laden magma with an unusual concentration of cpx crystals of ∼1.2 mm in size (i.e. L ∼1.2 mm). This is seen in the characteristic local maximum in the CSDs for the chill margins (Fig. 26a; see also Fig. 11d). This bump in the CSD may represent an earlier interval of annealing in a holding region or cumulate pile followed by ascent and further nucleation. Aside from this feature (and perhaps a similar less distinct one at L ∼ 2.5 mm), the initial CSDs are fairly linear. (The downturn at the smallest sizes is almost certainly due to the fact that these were measured by tracing from projected images, which underestimates the smallest crystals.) The similarity of the upper and lower CSDs suggests that the initial magma was fairly homogeneous. Once emplaced, some phenocrysts were captured during settling by the rapidly advancing upper solidification front, but eventually settling outstripped capture and a thick cumulate pile formed on the floor. [Other details of this process have been discussed by Marsh et al., (1991).]

This process of crystal capture, settling, and accumulation is clearly reflected in the CSDs. The CSDs of the cpx in the syenite become increasingly concave downward with depth, reflecting loss of crystals, and the CSDs in the cumulate pile are concave upward, reflecting crystal accumulation. In addition, resorption or annealing is reflected in the systematic decrease in nucleation density (i.e. no or intercept) within the syenite and uppermost lower shonkinite (cumulate pile). The CSD systematics are shown by Fig. 26b. The effects of resorption or annealing have been estimated by extrapolating the linear, central portion of each CSD to the L = 0 limit.

The slope vs intercept plot shows two trends beginning at the chilled margins. One trend shows increasing slope with relatively little change in intercept, which is due to crystal loss and accumulation. The other trend shows a significant decrease in nucleation density (intercept), which reflects resorption or annealing. (The slope values for both trends are the same because the slope of the CSD in the midsize range is used for both.) These trends are unlike those seen for the comagmatic sequences of Atka, Alaska, and Dome Mountain, Nevada.

When plotted against height or stratigraphic position in the laccolith, the strong steepening of the slope in the syenite is clearly evident. This is reminiscent of the slope variation in the central region of the Peneplain sill, but the cause here is different. In the Peneplain sill, the slope increase is suspected to be due to increased cooling rate, whereas in Shonkin Sag it is due to extensive resorption during protracted cooling. The residual cpx crystals in the syenite become increasingly ragged and resorbed in this region (Congdon, 1990). Not only are small grains lost but large grains are also resorbed. Based on the remaining crystal forms, Congdon, (1990) was able to restore the full crystals and show that the restored CSD resembles that of the chilled margins. The least affected crystals, however, seem to be those of intermediate size, as would be expected in annealing. The reduction in Lmax by ∼0.7 mm over the final, post-emplacement solidification time (∼1 year), suggests a resorption rate of 10−9 cm/s.

In sum, Shonkin Sag CSDs reflect an extensive history of pre-emplacement crystallization followed by extensive crystal settling, resorption, and accumulation.

#### Annealing theory

In terms of batch crystallization, annealing represents the effects of an effective reverse in nucleation rate and resorption of crystals smaller than some critical size and continued growth of crystals larger than this critical size, which itself, because of its inverse relation to the degree of 7lsquo;supersaturation’, migrates slowly (as ∼t1/3) to larger sizes (Lifshitz & Slyozov, 1961; see also Cashman & Ferry, 1988). The relative change in crystal stability is due to surface free energy. Small grains, because of their intrinsically high surface free energy, are less stable than large grains. The approach of Lifshitz & Slyozov is inherently the same as used here for a batch system. The central difference is that they based the entire result on a specific prescribed growth model dependent on surface free energy effects relative to a critical size of crystal. They solved for the final CSD and showed that the only stable distribution is a universal function (in normalized coordinates). The form of their general solutions [e.g. the fourth equation above their equation (15)] is the same as that used here [e.g. equation (39b)]. But because they prescribed the functional form of the growth rate, through the constraints of continuity [i.e. equation (37) here] and conservation of mass they obtained, in effect, a result describing the reverse nucleation rate over time and space. The universal function describing the distribution of grain size has grains no larger than 1.5 times the critical size. The CSD takes on this form regardless of the initial form of the CSD. At large times the number of crystals decreases asymptotically with time as t−1 and the critical grain size increases as t1/3, but the shape of the distribution of crystal sizes, when normalized to the most abundant and largest crystal classes remains, in this limit, invariant with time. Metamorphic rocks show distributions of grain size that approximate this ideal distribution, which is not of the log-linear form [see discussion by Cashman & Ferry, (1988)]. The log-linear part of the CSD of the Box Elder cpx and the Stillwater chromite can be analyzed for the mass lost in annealing using the straightforward CSD-moment functions as demonstrated by Cashman & Ferry, (1988). The mass of small crystals lost in the Box Elder CSD of the chilled margin does indeed approximate that gained at large crystal sizes in the interior CSD.

In sum, there are clear signs that annealing or ripening is an efficient process in igneous CSDs at all stages of crystallization. Some of the more dramatic effects occur in high-level intrusions formed of phenocryst-laden magma. It is important not to confuse these effects with the attenuation of nucleation as crystallization wanes, which is generally not seen in igneous rocks.

### CSD curvature vs composite linearity

Upward curvature in CSDs is a clear reflection of the physical addition of crystals to an evolving CSD. An example of this curvature is found in the picrite phase of the 1959 Kilauea eruption leading to the establishment of Kilauea Iki lava lake (see Fig. 27; data from Cashman, 1986). The pronounced upward concavity in the CSD for olivine suggests a multistage CSD history. Crystals may have been added to ascending magma from other local crystallization regimes as depicted by Fig. 21a. These added crystals may also accurately record the CSD systematics of these local regimes. On the other hand, this CSD may reflect a long and varied, but continuous and batch-like, evolution for the ascending magma where changing thermal regimes induce different nucleation regimes. In the former instance, the CSD may be piecewise continuous, being a collection of linear CSDs, and in the latter, the CSD may show continuous curvature. Two possible interpretations are given with the data of Fig. 27.

Fig. 27.

Two possible interpretations of the CSD of olivine in a picritic basalt of the Kilauea Iki eruption (data of Cashman, 1986), Hawaii, with inset (left) of broader CSD from Mangan, (1990).

Fig. 27.

Two possible interpretations of the CSD of olivine in a picritic basalt of the Kilauea Iki eruption (data of Cashman, 1986), Hawaii, with inset (left) of broader CSD from Mangan, (1990).

First, three straight lines have been put through the data (by eye). The first two lines fit the data well and broadly resemble some of the kinked CSDs shown already. The third straight line, at the largest sizes, goes through apparently rather noisy data. Be this as it may, this linear CSD (larger sizes) is nearly identical to that measured by Mangan, (1990) using the air quenched eruption pumices. Her measurements cover a much wider size range, up to ∼2.5 mm, and show much less unevenness (see inset to Fig. 27). These three linear CSDs could thus reflect a deep cumulate phase (largest L), nucleation during continued ascent, and finally strong nucleation attending eruption (smallest L).

On the other hand, a simple model relating crystal size to the continuous change of nucleation (and possibly growth) can be gauged from equations (2) and (24). That is,

(2)
$L=CL(GoJo)1/4=CL(no)1/4$
or
(73)
$ln (no)∼4ln (CL)−4ln (L).$
For a true batch system where Jo and Go change with time, no is an accurate measure of n. Curves calculated using (73) for a series of constants (0.25 ≤ CL ≤ 1.0) are also shown by Fig. 27 along with the Kilauea Iki data. For CL ∼ 0.5 the calculated curve resembles the data. The match could be made closer if CL also varied with size. A slightly smaller CL at small L and a larger CL at large L is needed. A comparison of these values of CL with the earlier kinetic model calculations summarized by Fig. 2 shows that these correspond to ab values of ∼1–4. The CSDs themselves, however, show a nucleation amplitude of ∼6 or more for all size ranges (considering also Mangan's data). This might suggest that a model of this type is inadequate in explaining these data. Instead, a CSD associated with the magmatic column of Fig. 21a a would seem to be more reasonable.

## Discussion

CSDs in comagmatic sequences are valuable records of magmatic behavior. As a rule, they do not closely resemble CSDs expected in either batch or open systems. Instead, they resemble the dynamic style of the system. That is, in systems where batch conditions are known on independent grounds to prevail (e.g. Peneplain sill, Makaopuhi lava lake, etc.) the CSDs reflect either the thermal regime of the solidification fronts (Peneplain) or local kinetic controls on crystallization (Makaopuhi). The contrast in CSDs of the Peneplain sill and Makaopuhi lava lake is particularly noteworthy.

In the Peneplain sill, nucleation rate apparently closely followed a conductive-style thermal regime with inward propagating solidification fronts. The holocrystalline final product clearly shows this CSD evidence because the system was unable to undergo prolonged, near-solidus annealing. Although the lava lakes probably also show this broad overall CSD variation, when examined in situ as the solidification front passes through a relatively restricted spatial dimension, the detailed kinetic effects of grain-induced heterogeneous nucleation and annealing through crystal annexation are clear. In the final product, however, these relatively subtle effects may not be at all obvious (e.g. Peneplain sill).

In systems where the dynamic style is not known, as is so for volcanic systems, there may be a dynamic distinction based on overall degree of crystallinity. In low-crystallinity systems, such as Dome Mountain, Nevada, the CSDs do resemble, at least in part, a batch system where the intercept increases as slope is more or less constant. And although crystallinity is low, it does increase with increasing intercept and modestly increasing maximum crystal size. These diagnostics, however, are by no means clear for the full lava sequence.

At high crystallinities in volcanic systems the CSDs inevitably tend to be either singly or multiply kinked. The high-alumina basalts of Atka, Alaska, show two distinct stages of plagioclase nucleation and growth. The CSD signal is coherent enough to suggest magmatic parcels dominated by separate deep-seated andnear-surface processes. In essence, these CSDs seem transitional between the Dome Mountain style of single-regime CSDs and the multi-regime systems of Hawaii (Kilauea Iki picrite) and Mount Etna, as shown so well by Armienti et al., (1994). The Mount Etna system, in particular, shows one property that essentially defines open systems. The CSD, at least for this century, seems to be invariant with time. Moreover, Etna may resemble the crystallization regimes seen in industrial crystallizers except where the output contains cumulate phases (see Fig. 21b). This style of system is thus open in a continuously intermittent sense where the underlying mush column is clearly reflected in the continuous entrainment of older, cumulate crystals. Thus the spatial effects of thermal regimes, which are so apparent in the CSDs of systems sampled in situ, are not nearly as evident in volcanic systems, although some systems must surely show them. Moreover, if sill-forming magma is crystal laden, subsequent slow cooling may lead to significant resorption or annealing.

There are three additional features of igneous CSDs that typify both igneous rocks and magmatic systems in general. The first is the apparently meaningful correlations between Lmax and CSD slope and intercept. Some reflection will convince one that Lmax is an independent measure of the coherence of the CSD. That is, given a well-developed CSD, crystals can easily be lost or gained (e.g. xenocrysts) such that the resulting CSD does not show systematic correlations with Lmax. In both batch and open systems the age of the CSD is reflected in the slope of the CSD and also in the size of the largest crystals. In purely dimensional terms the decrease of the CSD slope should always be proportional to Lmax. There are, nevertheless, a wide range of possible correlations between Lmax and inverse slope (i.e. |S1|), and it is important to consider these briefly (see Fig. 28). The two possible extremes are (1) systems where Lmax is constant and slope changes, and (2) where slope is constant and Lmax changes. In the first instance, crystals beyond a certain size could be lost from the system as it evolves as either a batch or open system (Fig. 28a and 28b). Another example is the annealing of a given population of phenocrysts where the largest crystals are mostly unaffected. In the alternative situation of constant slope and changing Lmax, the most obvious example is that of a strict batch system (Fig. 28a). A second, much less likely, example might be the approach to steady state of an open system under constant nucleation and growth rates (Fig. 28b).

Fig. 28.

Overview of possible diagnostic correlations between maximum crystal size and CSD inverse slope (absolute value) for both batch and open systems.

Fig. 28.

Overview of possible diagnostic correlations between maximum crystal size and CSD inverse slope (absolute value) for both batch and open systems.

The most reasonable case, however, is a progressive evolution of both Lmax and slope (Fig. 28d and 28e), which is exhibited by most igneous suites (see Fig. 29). These latter plots show fairly good correlations between Lmax and inverse slope with indications of two extremes. The one extreme of increasing inverse slope with constant Lmax is shown by Shonkin Sag. This reflects the post-emplacement annealing or resorption of the initial population of cpx phenocrysts. On the other hand, the Dome Mountain lavas show the steepest increase of Lmax with inverse slope, which again suggests batch behavior for this system. The other examples (Peneplain and Atka) show progressively changing Lmax and inverse slope, which may well represent a fundamental nucleation feature of igneous systems (see below).

Fig. 29.

Correlations between maximum crystal size and inverse CSD slope for four comagmatic systems.

Fig. 29.

Correlations between maximum crystal size and inverse CSD slope for four comagmatic systems.

The second distinctive feature of igneous CSDs is the apparent variation of nucleation rate. It is clear that the ln-linear (or segmented ln-linear) form of igneous CSDs comes at face value from an effective exponential rise in nucleation, regardless of cause, and simultaneous crystal growth. What is not seen, at least so far, are signs of either a constant rate of nucleation or a sustained hiatus in nucleation, both of which, on purely thermal grounds, might be expected. Recalling the examples of Fig. 11, either of these processes should be readily observable in igneous CSDs. What is seen instead, as exhibited, for example, by the Makaopuhi CSDs (Fig. 13), is a progressive fanning in the CSDs over the most easily measured size range of the CSD. As Lmax increases slope lessens. The net result is a progressive decrease in the effective nucleation rate with increasing crystal size. And because crystal size increases with time, this implies a progressive decrease in effective nucleation rate with time. This is, of course, aside from the possible overriding effect of a strongly increasing local cooling rate that may increase nucleation. But even for Makaopuhi, where increased cooling rate is occurring, increasing Lmax with decreasing effective nucleation rate dominates the system.

The fact that these are effective nucleation rates, uncorrected for melt fraction, must also be considered. If the Makaopuhi effective nucleation rates are corrected for decreasing melt fraction, nucleation rate decreases and reaches an approximate constant level only for advanced crystallinities (75%; Cashman & Marsh, 1988). That the effects of decreasing melt content are not large in terms of discussing the effective nucleation rates is shown by Fig. 4, where the CSD is essentially unaffected until crystallinity is well beyond 50%. This is the reason why these effects are not seen routinely in CSDs of volcanic systems. These effects, are, however, evident in the Peneplain sill CSDs where the rocks are holocrystalline. The CSDs at the smallest crystal sizes apparently reach a maximum and turn downward. The overall drop in population density in the Peneplain CSDs suggests that ab ∼12, whence, from the ideal batch CSDs of Fig. 4, the magnitude of the decrease in nucleation seems appropriate. This is especially so in light of the fact that plagioclase is but one of several minerals present.

The third distinctive feature of these CSDs is the apparent absence of a strong record of progressive increase in nucleation rate with time from one comagmatic CSD to another. For all intents and purposes, nucleation seems to begin almost full blown at high nucleation rates. This is particularly evident in CSDs in systems with pre-eruption low crystallinity and small crystals, such as the Mount St Helens dacites studied by Cashman, (1988, 1992). Once established, nucleation seems to sustain itself through perhaps heterogeneous nucleation effects; local saturations as a result of growth of foreign phases sustain nucleation and promote growth by annexation. If so, this style of nucleation may simply mean that igneous systems, by their very nature, are always laced with nuclei. Igneous systems have little chance of becoming superheated and destroying their nuclei. Igneous systems are intrinsically dirty in this respect by always being in contact with crystalline materials, indigenous or exotic. A spectrum of sizes of nuclei and clusters may always exist in magmatic systems such that, under suitable phase equilibria and cooling, crystal growth and further nucleation sets in. Because of the great sensitivity of nucleation to cooling rate, this possible characteristic may be difficult to discern in silicate glasses.

It is therefore possible that the typical ln-linear CSD form for igneous rocks does not typically come from simply exponential rises in nucleation rate, but instead from secondary processes, such as crystal annexation or growth rate dispersion, operating under more or less steady rates of nucleation. This issue is fundamental to understanding the crystallization of magmas. For we have seen [e.g. equation (53b)] that under constant nucleation, the negative log-linear CSD slope can be preserved if growth rate increases greatly with crystal size. Although it is clearly unreasonable to expect any isolated crystal to grow in this fashion, it may not be unreasonable to find that growth by annexation may approximate such a relation (i.e. G ∼ eL). Thus individual neighboring crystals may grow by a diffusive process, but coalescence of the crystals follows a dramatically different growth law. This is a fertile area for further study.

In closing, it is interesting to compare this wide family of observed igneous CSDs with those calculated by Hort & Spohn [1991a, 1991b; see also Spohn et al., (1988)] for coupled cooling and crystallization in thin dikes and sills. These studies assumed crystallization in a binary eutectic governed by Avrami-style kinetics with both an initial and continuous undercooling and growth and nucleation functions with maxima. The earlier study assumed homogeneous nucleation throughout, whereas the more recent work employed heterocatalytic nucleation (i.e. nucleation on the surfaces of earlier crystals), which is only effective at subliquidus temperatures. The critical difference between these works and the present is in the choice of Avrami number (see the Introduction), which Hort and coworkers set at ∼105. This choice assures that crystallization and heat transfer are equally important in crystallization, which makes analytical calculations of CSDs involved. Nevertheless, these numerical studies are valuable for elucidating the detailed, nonlinear interaction between kinetics and cooling. In this respect, the CSDs are particularly revealing. Although there is a good deal of variety depending on nucleation style, relative to natural CSDs, none of the CSDs show strong linearity. Instead, they more resemble annealed CSDs, such as shown by the Stillwater chromites (Fig. 25); there is commonly a pronounced downward trend in the CSD at small sizes. When nucleation includes hetercatalytic effects there are two distinct CSD populations: one produced at the liquidus, where nucleation is more homogeneous, and another at lower temperatures. Although these CSDs do not generally closely resemble natural CSDs they are important gauges for what to expect when the Avrami number is in this range.

## Conclusions

Log-linear or near log-linear CSDs (crystal size distributions) are commonplace in igneous rocks and would appear to represent an exponentially increasing effective nucleation rate with time accompanied by steady crystal growth. Nucleation rates are evidently sensitive to cooling rates, whereas growth rates are far less sensitive. The consequence of this is that crystals in igneous rocks do not vary extremely in size (roughly ∼0.1–10 mm) whereas nucleation rate, reflected in numbers of crystals, varies enormously. Because both nucleation and growth produce solid material, crystallization proceeds efficiently by either means. Conversely, all scales (i.e. time, length, and number) of crystallization depend on both nucleation and growth rates, and cannot be uniquely separated from one another except with additional information. Moreover, a wide variety of combinations of models of growth and nucleation rates show broadly similar sigmoidal increases of crystallinity with time. The characteristic crystallization times (tc) are, for the most part, similar regardless of the rates of nucleation and growth; that is, tc ≈ (Go3Jo)−1/4. Thus a cooling rate prescribed by external factors (i.e. body size, emplacement depth, etc.) is, within reason, easily adhered to by crystallization.

Direct insight into the relative roles of nucleation and growth comes from the CSDs, which also often allow a critical evaluation of a multitude of chemical and physical processes. From CSD studies made so far where some independent measure of cooling time is known, it appears that the effective growth rate for individual crystals on a local scale may be approximately constant in slowly cooling systems. The expected curvature in the CSD from nonlinear growth does not appear to be there. The inherent multiply-saturated nature of many igneous systems, however, strongly suggests, in concert with the CSDs, that growth is greatly aided by annexation of small grains by larger ones with subsequent grain boundary migration to produce optically continuous final minerals. When coupled with the effects of grain sorting, annealing, and heterogeneous nucleation, caused by multiple saturation, cooling and decompression, distinctive signatures of the CSDs suggest a continuum of magmatic evolution from single-stage fissure eruptions to well-established mush columns. Overall, the basic ln-linear form of typical igneous CSDs may owe itself less to exponential changes in nucleation rate than to secondary processes, such as grain annexation, operating under a fairly steady nucleation rate set by the local cooling regime.

## Acknowledgements

This work has benefited from intermittent discussions with many people interested in CSD analysis. Among these people are O. M. Phillips, K. V. Cashman, M. Hort, R. G. Resmini, M. T. Mangan, A. T. Lasaga, K. A. McCormick, and G. Bergantz. Helpful reviews byD. Jerram, K. V. Cashman, J. Crisp, A. Boudreau, M. Higgins, F. Gibb, M. Hort, R. Resmini, K. Schwindinger and K. McCormick are much appreciated. Special thanks are due to C. Weeks and M. Zieg for help with the manuscript and figures. This work is supported by grants from the National Science Foundation (DPP-9418513) and the Nuclear Regulatory Commission (NRC-04-94-96).

## References

Armienti
P.
Pareschi
M. T.
Innocenti
F.
Pompilio
M.
Effect of magma storage and ascent on the kinetics of crystal growth
Contributions to Mineralogy and Petrology
,
1994
, vol.
115
(pg.
402
-
414
)
Avrami
M.
Kinetics of phase change I
Journal of Chemical Physics
,
1939
, vol.
7
(pg.
1103
-
1112
)
Avrami
M.
Kinetics of phase change II
Journal of Chemical Physics
,
1940
, vol.
8
(pg.
212
-
224
)
Boudreau
A. E.
Mineral segregation during crystal aging in two-crystal, two-component systems
South African Journal of Geology
,
1994
, vol.
4
(pg.
473
-
485
)
Brandeis
G.
Jaupart
C.
Mysen
B. O.
Crystal sizes in intrusions of different dimensions: constraints on the cooling regime and the crystallization kinetics
Magmatic Processes: Physicochemical Principles
,
1987
University Park, PA
Geochemical Society
(pg.
307
-
318
)
Brandeis
G.
Jaupart
C.
The kinetics of nucleation and crystal growth and scaling laws for magmatic crystallization
Contributions to Mineralogy and Petrology
,
1987
, vol.
96
(pg.
24
-
345
)
Brandeis
G.
Jaupart
C.
Allègre
C. J.
Nucleation, crystal growth and the thermal regime of cooling magmas
Journal of Geophysical Research
,
1984
, vol.
89
(pg.
10161
-
10177
)
Burkhart
L. E.
Hoyt
R. C.
Oolman
T.
Kuczynski
G. C.
Control of particle size distribution and agglomeration in continuous precipitations
Sintering Processes
,
1980
New York
Plenum
(pg.
23
-
38
)
Canning
T. F.
Interpreting population density data from crystallizers
Chemical Engineering Progress
,
1970
, vol.
67
(pg.
74
-
80
)
Carmichael
I. S. E.
Yoder
H. S.
Crystal sizes in intrusions of different dimensions: constraints on the cooling regime and the crystallization kinetics
The Evolution of the Igneous Rocks
,
1979
Princeton, NJ
Princeton University Press
(pg.
233
-
244
)
Cashman
K. V.
Crystal size distributions in igneous and metamorphic rocks
The Evolution of the Igneous Rocks
,
1986
Princeton, NJ
Johns Hopkins University
pg.
345

Ph.D. Dissertation
Cashman
K. V.
Crystallization of Mount St. Helens 1980–1986 dacite: a quantitative textural approach
Bulletin of Volcanology
,
1988
, vol.
50
(pg.
194
-
209
)
Cashman
K. V.
Nicholls
J.
Russell
J. K.
Textural constraints of the kinetics of crystallization of igneous rocks
Modern Methods of Igneous Petrology: Understanding Magmatic Processes. Mineralogical Society of America, Reviews in Mineralogy
,
1990
, vol.
24
(pg.
259
-
314
)
Cashman
K. V.
Groundmass crystallization of Mount St. Helens dacite 1980–1986: a tool for interpreting shallow magmatic processes
Contributions to Mineralogy and Petrology
,
1992
, vol.
109
(pg.
431
-
449
)
Cashman
K. V.
Relationship between plagioclase crystallization and cooling rate in basaltic melts
Contributions to Mineralogy and Petrology
,
1993
, vol.
113
(pg.
126
-
142
)
Cashman
K. V.
Ferry
J. M.
Crystal size distribution (CSD) in rocks and the kinetics and dynamics of crystallization III. Metamorphic crystallization
Contributions to Mineralogy and Petrology
,
1988
, vol.
99
(pg.
401
-
415
)
Cashman
K. V.
Marsh
B. D.
Crystal size distribution (CSD) in rocks and the kinetics and dynamics of crystallization II: Makaopuhi lava lake
Contributions to Mineralogy and Petrology
,
1988
, vol.
99
(pg.
292
-
305
)
Congdon
R. D.
The solidification of the Shonkin Sag laccolith: mineralogy, petrology, and experimental phase equilibria
The Evolution of the Igneous Rocks
,
1990
Princeton, NJ
Johns Hopkins University
pg.
362

Ph.D. Dissertation
Congdon
R. D.
Resmini
R. G.
Marsh
B. D.
Differentiation style in the Box Elder and Shonkin Sag Laccoliths; dependence on initial conditions
EOS Transactions, American Geophysical Union
,
1993
, vol.
73
pg.
336

Crisp
J.
Cashman
K. V.
Bonini
J. A.
Hougen
S. B.
Pieri
D. C.
Crystallization history of the 1984 Mauna Loa lava flow
Journal of Geophysical Research
,
1994
, vol.
99
(pg.
7177
-
7198
)
Dawson
J. B.
First thin sections of experimentally melted igneous rocks: Sorby's observations on magma crystallization
Journal of Geology
,
1992
, vol.
100
(pg.
251
-
257
)
Dowty
E.
Hargraves
R. B.
Crystal growth and nucleation theory and the numerical simulation of igneous crystallization
The Physics of Magmatic Processes
,
1980
Princeton, NJ
Princeton University Press
(pg.
419
-
485
)
Evans
B. W.
Moore
J. G.
Olivine in the prehistoric Makaopuhi tholeiitic lava lake, Hawaii
Contributions to Mineralogy and Petrology
,
1967
, vol.
15
(pg.
202
-
223
)
Gibb
F. G. F.
Supercooling in the crystallization of plagioclase from a basaltic magma
Mineralogical Magazine
,
1974
, vol.
39
(pg.
641
-
653
)
Gunn
B. M.
Modal and element variation in Antarctic tholeiites
Geochimica et Cosmochimica Acta
,
1966
, vol.
30
(pg.
881
-
920
)
Hall
J.
Experiments on whinstone and lava
Transactions of the Royal Society of Edinburgh
,
1798
, vol.
5
(pg.
43
-
76
)
Heyn
J.
Marsh
B. D.
Wheelock
M. M.
Crystal size and cooling time in the Peneplain sill, Dry Valley Region, Antarctica
Antarctic Journal of the United States
,
1997
, vol.
30
(pg.
50
-
51
)
Higgins
M. D.
Numerical modeling of crystal shapes in thin sections: estimation of crystal habit and true size
American Mineralogist
,
1994
, vol.
79
(pg.
113
-
119
)
Higgins
M. D.
Magma dynamics beneath Kameni volcano, Thera, Greece, as revealed by crystal size and shape measurements
Journal of Volcanology and Geothermal Research
,
1996
, vol.
70
(pg.
37
-
48
)
Hort
M.
Spohn
T.
Crystallization calculations for a binary melt cooling at constant rates of heat removal: implications for the crystallization of magma bodies
Earth and Planetary Science Letters
,
1991
, vol.
107
(pg.
463
-
474
)
Hort
M.
Spohn
T.
Numerical simulation of the crystallization of multicomponent melts in thin dikes or sills 2. Effects of heterocatalytic nucleation and composition
Earth and Planetary Science Letters
,
1991
, vol.
96
(pg.
485
-
499
)
Hunter
R.
Parsons
I.
Textural equilibrium in layered igneous rocks
Origins of Igneous Layering
,
1986
Dordrecht
D. Reidel
(pg.
473
-
503
)
Hurlbut
C. S. J.
Igneous rocks of the Highwood Mountains, Montana. Part I. The laccoliths
Geological Society of America Bulletin
,
1939
, vol.
50
(pg.
1043
-
1112
)
Jaeger
J. C.
Hess
H. H.
Poldervaart
A.
Cooling and solidification of igneous rocks
Basalts: The Poldervaart Treatise on Rocks of Basaltic Composition
,
1968
New York
Interscience
(pg.
503
-
536
)
Jancic
S. J.
de Jong
E. J.
Industrial Crystallization 84
,
1984
Amsterdam
Elsevier
(pg.
498
-
536
)
Kirkpatrick
R. J.
Lasaga
A. C.
Kirkpatrick
R. J.
Kinetics of crystallization of igneous rocks
Kinetics of Geochemical Processes
,
1983
Washington, DC
Mineralogical Society of America
(pg.
321
-
398
)
Kress
V. C.
Ghiorso
M. S.
Multicomponent diffusion in MgO–Al2O3–SiO2 and CaO–MgO–Al2O3–SiO2 melts
Geochimica et Cosmochimica Acta
,
1993
, vol.
57
(pg.
4453
-
4466
)
Kress
V. C.
Ghiorso
M. S.
Multicomponent diffusion in basaltic melts
Geochimica et Cosmochimica Acta
,
1995
, vol.
59
(pg.
313
-
324
)
Liang
Y.
Richter
F. M.
Davis
A. M.
Watson
E. B.
Diffusion in silicate melts: I. Self diffusion in CaO–Al2O3–SiO2 at 1500°C and 1 GPa
Geochimica et Cosmochimica Acta
,
1996
, vol.
60
(pg.
4353
-
4367
)
Lifshitz
I. M.
Slyozov
V. V.
The kinetics of precipitation from supersaturated solid solutions
Journal of Physics and Chemistry of Solids
,
1961
, vol.
19
(pg.
35
-
50
)
Lofgren
G.
Hargraves
R. B.
Experimental studies on the dynamic crystallization of silicate melts
Physics of Magmatic Processes
,
1980
Princeton, NJ
Princeton University Press
(pg.
487
-
551
)
Maaloe
S.
Tumyr
O.
James
D.
Population density and zoning of olivine phenocrysts in tholeiites from Kauai, Hawaii
Contributions to Mineralogy and Petrology
, vol.
101
(pg.
176
-
186
)
Mangan
M. T.
Crystal size distribution systematics and the determination of magma storage times: the 1959 eruption of Kilauea volcano, Hawaii
Journal of Volcanology and Geothermal Research
,
1990
, vol.
44
(pg.
295
-
302
)
Mangan
M. T.
Marsh
B. D.
Solidification front fractionation in phenocryst-free sheet-like magma bodies
Journal of Geology
,
1992
, vol.
100
(pg.
605
-
620
)
Marsh
B. D.
On the crystallinity, probability of occurrence, and rheology of lava and magma
Contributions to Mineralogy and Petrology
,
1981
, vol.
78
(pg.
85
-
98
)
Marsh
B. D.
Crystal size distribution (CSD) in rocks and the kinetics and dynamics of crystallization I. Theory
Contributions to Mineralogy and Petrology
,
1988
, vol.
99
(pg.
277
-
291
)
Marsh
B. D.
Solidification fronts and magmatic evolution
Mineralogical Magazine
,
1996
, vol.
60
(pg.
5
-
40
)
Marsh
B. D.
Wheelock
M. M.
The vertical variation of composition in the Peneplain sill and Basement sills of the Dry Valleys: the null hypothesis
Antarctic Journal of the United States
,
1994
, vol.
29
(pg.
25
-
26
)
Marsh
B. D.
B.
Congdon
R.
Carmody
R.
Hawaiian basalt and Icelandic rhyolite: indicators of differentiation and partial melting
Geologische Rundschau
,
1991
, vol.
80
(pg.
481
-
510
)
McBirney
A. R.
Hawaiian basalt and Icelandic rhyolite: indicators of differentiation and partial melting
Igneous Petrology
,
1993
2nd
Boston, MA
Jones and Bartlett
pg.
508

Means
W. D.
Park
Y.
New experimental approach to understanding igneous texture
Geology
,
1994
, vol.
23
(pg.
323
-
326
)
Myers
J. D.
Marsh
B. D.
Sinha
A. K.
Geochemical and strontium isotopic characteristic of parental Aleutian Arc magmas: evidence from the basaltic lavas of Atka
Contributions to Mineralogy and Petrology
,
1986
, vol.
94
(pg.
1
-
11
)
Pecora
W. T.
Structure and petrology of the Boxelder laccolith, Bearpaw Mountains, Montana
Geological Society of America Bulletin
,
1941
, vol.
52
(pg.
817
-
854
)
Peterson
T. D.
A refined technique for measuring crystal size distributions in thin section
Contributions to Mineralogy and Petrology
,
1996
, vol.
124
(pg.
395
-
405
)
Randolph
A. D.
Larson
M. A.
Theory of Particulate Processes
,
1971
New York
pg.
251

Randolph
A. D.
Larson
M. A.
Theory of Particulate Processes
,
1988
2nd
New York
pg.
369

Resmini
R. G.
Dynamics of magma within the crust: a study using crystal size distribution
1993
New York
Johns Hopkins University
pg.
329

Ph.D. Dissertation
Resmini
R. G.
Marsh
B. D.
Steady-state volcanism, paleo-eVusion rates, and magma system volume inferred from plagioclase crystal size distributions in mafic lavas: Dome Mountain, Nevada
Journal of Volcanology and Geothermal Research
,
1995
, vol.
68
(pg.
273
-
296
)
Shaw
H. R.
Comments on viscosity, crystal settling, and convection in granitic magmas
American Journal of Science
,
1965
, vol.
263
(pg.
120
-
152
)
Solomatov
V. S.
Batch crystallization under continuous cooling: analytical solution for diffusion limited crystal growth
Journal of Crystal Growth
,
1995
, vol.
148
(pg.
421
-
431
)
Spohn
T.
Hort
M.
Fischer
H.
Numerical simulationof the crystallization of multicomponent melts in thin dikes or sills1. The liquidus phase
Journal of Geophysical Research
,
1988
, vol.
93
(pg.
4880
-
4894
)
Tomiya
A.
Takahashi
E.
Reconstruction of an evolving magma chamber beneath Usu Volcano since the 1663 eruption
Journal of Petrology
,
1995
, vol.
36
(pg.
617
-
636
)
Tomkeieff
S. I.
The basalt lavas of the Giant's Causeway district of Northern Ireland
Bulletin of Volcanology
,
1940
, vol.
6–7
(pg.
89
-
146
)
Turcotte
D. L.
Fractals and fragmentation
Journal of Geophysical Research
,
1986
, vol.
91
(pg.
1921
-
1926
)
Waters
C.
Boudreau
A. E.
A reevaluation of crystal size distributions in chromite cumulates
American Mineralogist
,
1996
, vol.
81
(pg.
1452
-
1459
)
Wilhelm
S.
Worner
G.
Crystal size distribution in Jurassic Ferrar flows and sills (Victoria Land, Antarctica): evidence for processes of cooling, nucleation, and crystallisation
Contributions to Mineralogy and Petrology
,
1996
, vol.
125
(pg.
1
-
15
)
Winkler
H. G. F.
Crystallization of basaltic magma as recorded by variation of crystal sizes in dikes
Mineralogical Magazine
,
1949
, vol.
28
(pg.
557
-
574
)
Wright
T. L.
Fiske
R. S.
Origin of the differentiated and hybrid lavas of Kilauea volcano, Hawaii
Journal of Petrology
,
1971
, vol.
12
(pg.
1
-
65
)
Zanotto
E. D.
Galhardi
A.
Experimental test of the general theory of transformation kinetics: homogeneous nucleation in a NaO2.2CaO2.3SiO2 glass
Journal of Non-Crystalline Solids
,
1988
, vol.
104
(pg.
73
-
80
)