Igneous microstructures from kinetic models of crystallization

Journal of Volcanology and Geothermal Research 154 (2006) 34 – 47 www.elsevier.com/locate/jvolgeores

Igneous microstructures from kinetic models of crystallization Taber G. Hersum ⁎, Bruce D. Marsh Johns Hopkins University, Morton K. Blaustein Department of Earth and Planetary Sciences, 34th and North Charles Streets, Baltimore, MD 21218, USA Received 29 October 2004; accepted 30 September 2005 Available online 2 March 2006

Abstract During crystallization, a variety of competing kinetic processes determine the evolution of the igneous microstructure, yet, the relative contribution of each process remains elusive. To this end, a stochastic algorithm is developed to yield a detailed spatial representation of the igneous microstructure during progressive crystallization. This algorithm is used to test a variety of kinetic models for nucleation and crystal growth that yield realistic igneous microstructures. The algorithm itself relies on a theoretical model of crystallization (the Avrami method) and can therefore be internally validated in addition to being verified against natural microstructures using crystal size distributions (CSDs). The most realistic simulated igneous microstructure, in which the CSDs remain approximately log-linear throughout the crystallization interval, is produced using a kinetic model involving exponential nucleation rate and constant growth rate. An extension of the algorithm is used to simulate multiply saturated mineral growth, emulating basalt crystallization through simultaneous nucleation and growth of plagioclase and clinopyroxene. A fundamentally important feature of this procedure is the numerical production of detailed 3D representations of the microstructure that can be interrogated to study many physical and chemical processes. For example, the critical crystallinity necessary for the onset of a finite yield strength and the interstitial melt flow within the igneous microstructure depend on the percolation of the solid and melt phases, respectively. Bounds on the percolation thresholds, the critical crystallinity at which the phase of interest is connected or contiguous such that it extends across the full microstructure from one face to the opposite face, are determined for the simulated igneous microstructure and are comparable with previously published results for the percolation threshold of both the melt phase and solid network. The simulated igneous microstructure can also be used as input into other physical models to ascertain the physical properties of partially molten magmas that are otherwise difficult to estimate by experiment. In a real sense, these computed microstructures are equivalent to 3D microtomographic images of partially molten basalt within a solidification front. © 2006 Elsevier B.V. All rights reserved. Keywords: crystallization kinetics; stochastic simulation; basaltic microstructure; percolation

1. Introduction The extreme variability of natural igneous microstructures in terms of crystal size, shape, number, and position is due to the interplay of competing kinetic ⁎ Corresponding author. Fax: +1 410 516 7933. E-mail address: [email protected] (T.G. Hersum). 0377-0273/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jvolgeores.2005.09.018

processes operating during crystallization (e.g. nucleation, growth, aggregation, and grain boundary annealing). It has proven challenging to decipher the contribution of each kinetic process to the final microstructure, including the role of the intensive parameters (e.g. temperature, time) governing crystallization. Present knowledge of the crystallization kinetics of solidifying magmas is mainly limited to relatively

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

short duration melting experiments that are most applicable to rapidly cooled magmas (i.e., lavas) and not coarse-grained, slowly cooled igneous rocks that record a fuller and more dynamic history of magmatic activity. To augment empirical results, some of which have been available from the late 1700s (e.g. Hall, 1798), a theory of crystallization kinetics was adapted from classical kinetic theory (e.g. Kirkpatrick, 1981), but this has proven too broad based to be of immediate practical use in predicting the styles of igneous microstructures. That is, this atomic scale theory is inherently incomplete for silicate melts and, in the light of sparse kinetic information from actual magmatic systems, is exceedingly awkward for direct application as a forward model of microstructural development. An alternative approach is to view igneous microstructural development as a grain-scale process with fundamental rules describing the transient appearance and morphological evolution of crystals, which is consistent with but does not resolve atomic scale processes of diffusion and bonding during crystallization. This rule-based kinetic model of crystallization relies on long established general observations from crystallization experiments, from post-mortem studies of solidified magmas, and the results of both crystal size distribution (CSD) theory and CSD measurements on rocks. The thrust of the present work is to demonstrate that a stochastic numerical algorithm for forward modeling of igneous microstructural development can be used to delineate the relative contributions of individual kinetic processes in forming realistic igneous microstructures using various kinetic models of crystallization. In brief, a 3D space of melt containing randomly orientated and randomly positioned crystals is represented by an ensemble of discrete volumetric elements; changes in time, and increasing crystallinity, are achieved by discrete temporal steps. In assuming that cooling and crystallization can be decoupled, as discussed below, the design of the algorithm is greatly simplified and crystallinity can be expressed solely as a function of time. Functions describing the change of crystallinity with time are derived using the Avrami method by employing any number of exponential variations (including constant) in functions describing the rates of nucleation and crystal growth (Marsh, 1998). An investigation of the sensitivity of the algorithm to variations in the input parameters precedes the evaluation of the microstructures generated by the various kinetic models. Three basic kinetic models are investigated: constant rates of nucleation and growth, exponential nucleation rate with constant growth rate, and exponential nucleation rate with crystal-dependent

35

dispersive, but individually constant, growth rate. An extension of the algorithm is also demonstrated by modeling multiply saturated mineral growth, which simulates basalt crystallization involving simultaneous plagioclase and clinopyroxene nucleation and growth. Aside from understanding the kinetics of crystallization, the computed microstructures are invaluable as 3D samples of test materials necessary for input into models developed to obtain macro-scale physical properties under conditions unobtainable by experiment. Failure behavior, elastic moduli, permeability, and grain-scale melt flow patterns including flow channelization are critical for continuum scale models of magmatic processes and can, thus, be determined throughout the melting interval over a wide range of pressure. 2. Background Previous studies of igneous crystallization predicting crystal size and number were hampered by having to solve simultaneously for the thermal regime in conjunction with either scaling laws (e.g. Brandeis and Jaupart, 1987) or classical kinetic theories of nucleation and crystal growth (Spohn et al., 1988; Tomaru, 2001). These models are useful to compare with laboratory measurements, but because they are unable to predict the spatial representation of the evolving microstructure they cannot test kinetic models of crystallization for either accuracy or completeness of the final rock texture. Some studies have also modeled the evolution of igneous microstructure during crystallization using stochastic approaches and simple kinetic laws (Jerram et al., 1996; Amenta, 2001, 2004; Cheadle et al., 2004). The present model differs from previous models, however, in being based on a theoretical model of crystallization (the Avrami method) and is, thus, open to internal validation by allowing the algorithm to be used as a tool to test different analytical kinetic models of crystallization. That is, the numerical simulations can be tested against known analytical solutions. 3. Decoupling cooling and crystallization Most igneous rocks are holocrystalline (not glassy) regardless of the cooling regime in which they have formed. This suggests that for any given time scale imposed by cooling, nucleation and growth automatically adjust locally to attain full crystallinity (Marsh, 1998). Although this is not always strictly true, it is a very good approximation, especially for basaltic compositions. This directly implies that crystallization can be decoupled from cooling, which can be demonstrated

36

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

through scaling. That is, the time scale for crystallization of a volume of melt by nucleation and crystal growth is (e.g. Marsh, 1998):

 

tc ¼ Ct ðJo G 3o Þ−1=4

ð1Þ

where Jo˙ and Go˙ are, respectively, characteristic rates of nucleation and growth, and Ct is a constant of order one. The time scale for an igneous body to cool is: tco ¼ cðH 2 j−1 Þ

ð2Þ

where H is the half-thickness of the body, κ is thermal diffusivity, and c is a constant of order one. The ratio of the two time scales to the fourth power is the Avrami number, Av = (tco / tc)4 (Spohn et al., 1988), and for reasonable kinetic and cooling parameters Av is ∼1010, which means that the characteristic time for local crystallization is much, much smaller than the time for heat transfer. That is, crystallization can always immediately adjust, through changes in nucleation and growth, to any thermal regime. Our approach is, thus, greatly simplified by neglecting specific cooling regimes in favor of modeling crystallinity as an explicit function of time and not temperature. Several functional forms describing the increase in crystallinity with time due to phase equilibria are presented. 4. Kinetic model During crystallization, a variety of kinetic processes can influence the final igneous microstructure. The present approach is to consider two of the most basic processes, heterogeneous nucleation and crystal growth, and to test simple models of these processes in reproducing reasonable microstructures. The essential feature of this approach is the imposition of a function describing the variation of bulk crystallinity with time. The utility of employing a crystallinity function, in addition to the inherent simplicity introduced by decoupling cooling and crystallization, is that only nucleation or growth needs to be explicitly defined initially and the remaining rate is an outcome of the numerical simulation itself. In effect, the one unspecified parameter becomes a closure condition for the crystallization process. The obvious choice of the closure parameter is the rate of nucleation; it is enormously flexible in this regard. Over a wide range of crystallization conditions (e.g. total crystallization time) and even during the crystallization process itself, nucleation rate can vary over many (perhaps 15 or more) orders of magnitude; whereas it is unlikely that growth rate varies by more than a few orders of magnitude during any single

crystallization history (Cashman, 1993; Marsh, 1998). The present procedure uses the Avrami equation (Avrami, 1940) to obtain crystallinity over time by prescribing functions describing the rates of nucleation and growth (Marsh, 1998). These same functions are used as input into the algorithm to generate an igneous microstructure through the combined effect of nucleation and growth. The change in number of crystals is tracked during the simulation. The algorithm is validated by comparing the number of crystals found in the simulation against that found using the analytical function for nucleation rate. The evolution of the solid fraction over time is found by integrating over time the product of crystal volume fraction and nucleation rate. Although any particular crystal shape can be included in the Avrami formulation, rectangular prisms with fixed aspect ratios are used to approximate crystal shape. The derived crystallinity function, however, is much less sensitive to changes in crystal shape than to changes in nucleation and growth rate. Crystals can appear at any time t′ during the crystallization interval and therefore the length, D, of the long axis of a crystal with time is: Z t  DðtÞ ¼ 2 G ðtÞdt ð3Þ tV

Given the aspect ratios A1 and A2 that represent the long axis length of the crystal over the short and intermediate axis lengths of the crystal, respectively, the crystal volume is: Z t 3 8  V ðtÞ ¼ G ðtÞdt ð4Þ A1 A2 t V The nucleation rate determines the total number of crystals appearing at time t′ and the Avrami equation becomes: ( Z t 3 ) Z t 8     /ðt; J ; G Þ ¼ 1−exp − J ðt VÞ G ðtÞdt dt V A1 A2 0 tV ð5Þ Many functions for nucleation and growth rate can be accommodated into the Avrami equation, but particularly flexible ones are exponential functions that are common to kinetics:   at   J ðtÞ ¼ Jo exp ð6Þ tc





G ðtÞ ¼ Go exp

  bt tc

ð7Þ

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

where t is time, tc is the total crystallization time as defined in Eq. (1), and a and b are constants of either sign. The CSDs of many rocks resulting from bulk crystallization show a log-linear trend of population density with increasing crystal size. From CSD theory, a kinetic model assuming a constant crystal growth rate and an exponential nucleation rate can reconcile the observed CSD trends. This is not to suggest that this is the only combination to explain the observed CSD trends, but that a kinetic model of this nature is simple, the numerical results are, therefore, straightforward to interpret, and the result is consistent with observed CSDs. A model with constant crystal growth rate and exponential nucleation rate is also easily implemented into the algorithm, although more complex models can also be accommodated, including time-dependent growth rate (Burkhard, 2002), boundary layer thickening, and nuclei clustering. Eq. (5) can be recast using the Eqs. (6) and (7) and introducing dimensionless time, x = t / tc:    8J o G 3o tc4 /ðx; a; bÞ ¼ 1−exp − A1 A2 Z x 3  Z x  expðax VÞ expðbxVÞdx dx V xV

0

ð8Þ By substituting Eq. (1) into Eq. (8) and representing the two integrals as: (Z Z x 3 ) x f ðx; a; bÞ ¼ exp expðaxVÞ expðbxVÞdx dxV 0

xV

37

Fig. 1. Crystallinity versus dimensionless time given by Eq. (12) for constant growth rate (b = 0) and various degrees of exponential nucleation rate.

ation rate with various a values, therefore Eq. (9) becomes: f ðx; a; 0Þ ¼

1 f6 expðaxÞ−½ðaxÞ3 þ 3ðaxÞ2 þ 6ax þ 6g a4 ð13Þ

Eq. (12) is plotted with various a values in Fig. 1. The sigmoidal shape of the crystallinity function is consistent with experimental studies and with observations of magma crystallizing in lava lakes (Wright and Okamura, 1977; Marsh, 1981). 5. The stochastic algorithm

ð9Þ Eq. (8) becomes:   8Ct4 /ðx; a; bÞ ¼ 1−exp − f ðx; a; bÞ A1 A2

ð10Þ

By assigning the crystal fraction at which crystallization is complete (e.g. ϕc = 0.99), Ct as defined in Eq. (1) can be found as a function of a and b (Marsh, 1998):   −A1 A2 lnð1−/c Þ 1=4 Ct ða; bÞ ¼ 8f ð1; a; bÞ and Eq. (10) can be simplified to:    lnð1−/c Þ /ðx; a; bÞ ¼ 1−exp f ðx; a; bÞ f ð1; a; bÞ

ð11Þ

ð12Þ

In this study, we only consider constant growth rate for individual crystals (b = 0) and exponential nucle-

The crystallizing domain is a cubic volume of length L on a side that is discretized into a number of smaller cubic volume elements (voxels) each of width W, containing either solid or liquid (see Fig. 2). A crystal is composed of multiple adjacent solid voxels and the number of voxels per side of the domain is m. The crystallization time is discretized into a number of time steps nt of equal time duration Δt. At the beginning of a crystallization simulation no nuclei or crystals are present in the domain and all voxels are liquid. At each time step and given the current crystallization time, crystallinity is calculated from the crystallinity function and converted into a number of solid voxels out of the total number of voxels in the domain (see Fig. 3). For the first time step, a voxel is chosen at random as the nucleation site for a new crystal. The orientation of the new crystal in 3D space is also determined by three independent random numbers. A crystal length for the long axis is calculated by multiplying twice the

38

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

Fig. 2. 3D crystallizing domain discretized into a number of voxels that are either solid or liquid. A 2D slice of the 3D domain demonstrates how solid voxels are mapped onto a crystal boundary box.

instantaneous growth rate G˙ by Δt. Given the aspect ratios that characterize the crystal shape, a boundary box is created centered in the middle of the chosen nucleation site that represents the crystal boundary. The crystal volume is then discretized into a number of solid voxels by scanning through the nearby domain and assigning voxels whose centers fall within the boundary box of the new crystal (see Fig. 2). The creation of other new crystals continues in the same manner until the requisite crystallinity is achieved for the current crystallization time. After the first crystal is formed,

nucleation sites are only eligible in liquid voxels. During crystal growth, no previous solid voxels of one crystal type can be converted into another crystal type, which simulates the impingement of crystals and, thus, allows the entrapment of one crystal within another. Liquid voxels that fall within a crystal's boundary box but do not share any common vertices with the solid voxels of the crystal are not assigned to the growing crystal, which prevents a crystal from growing through another crystal. For subsequent time steps, crystallinity must first attempt to be satisfied by growth of all pre-existing

Fig. 3. Diagram of stochastic algorithm procedure. At t = 0 there is no crystallinity and after the first time step, t = tc / 4, a crystallinity of 11% is achieved through a nucleation event. During the second time step, all pre-existing nuclei are grown before new nuclei are formed to attain a crystallinity of 50% g at t = tc / 2.

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

crystals that are not entrapped. This is then followed by the creation (nucleation) of new crystals. If the requisite crystallinity for the current crystallization time is achieved before all of the pre-existing crystals are grown, no more crystals are grown, no new crystals are created, and the algorithm continues to the next time step. Pre-existing crystals cannot be grown simultaneously during a time step because the computational method requires crystal-by-crystal growth. However, simultaneous growth is approached by minimizing Δt (and, hence, the growth distance of a crystal) and by randomly choosing the order in which pre-existing crystals are grown at each time step. The new length of a growing crystal is determined by appending the previous crystal length with a calculated growth distance for the most recent time step. Crystallization does not occur across the boundaries of the overall initial melt domain. Changes in system volume associated with the phase change from a liquid to solid state are ignored and the size of the domain remains constant throughout the crystallization experiment. Crystallization occurs in a zero-shear environment and crystal movement during crystallization through mutual impingement or buoyancy forces is not incorporated into the algorithm. A realization refers to a single complete computation or experiment stemming from a given set of input parameters with a unique initial random seed; a simulation refers to a collection of realizations with the same input parameters. Because of the inherent randomness built into the algorithm, no two realizations will yield the same microstructure. Any simulation will therefore be characterized by a set of statistical measures.

39

least a distance of one voxel width. By assigning values for these dimensionless numbers, discretization of space and time are automatically coupled such that the ratio of domain length in voxels (m = L / W) over the number of time steps (nt = tc / Δt) equals 2β / α. The effect of spatial and temporal discretization on model accuracy is tested with monominerallic simulations involving fixed values of α and β but variable values for m and nt and comparing the numerically and analytically determined number of crystals. The analytical number of crystals during crystallization is found by integrating Eq. (6) and using Eq. (1) to replace Jo˙ with known quantities: 8 Ct4 L3 > > > x; for a ¼ 0 >  < ð Go tc Þ3 ð14Þ no V¼ > C 4 L3 > t > > expðaxÞ; for aN0 :  að Go tc Þ3 where no′ is the number of crystals in the domain and recall that x is dimensionless time. Fig. 4 shows for increasing m and nt that the final ensemble-averaged number of crystals determined numerically approaches the analytical result for simulations with α = 8, α = 0.3125, β = 1, and r = 3, where r is the number of realizations; although the relationship is more ambiguous for intermediate time steps. Unfortunately, memory restrictions of the CPU place an upper limit on m and therefore nt. The results presented here are based on discretization parameters that are a trade-off between the memory constraints of the CPU, which limit the grid fineness, and a decrease in step size to increase model accuracy.

6. Sensitivity analysis The simulated microstructure is, clearly, sensitive to the degree of spatial and temporal discretization used in the stochastic algorithm. Judicious choice of the discretization parameters is required for the simulation to approximate the analytical solutions for numbers and lengths of crystals (validation) and also to minimize error when the microstructure is input into other physical models. The number of crystals in the domain and the spatial resolution of crystals are controlled by two dimensionless numbers. The first, α = 2G˙tc / L, is the length of a crystal growing unhampered throughout the entire crystallization interval relative to the length of the domain itself. ˙ / W, is the growth distance of a The second, β = GΔt crystal during a single time step relative to the width of a voxel. Growth during each time step must be at

Fig. 4. Analytical (solid line) and numerical (ensemble-averaged) number of crystals for monominerallic simulations with m equal to 32, 64, and 96 but fixed α and β.

40

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

Model validation, using analytical solutions for numbers and lengths of crystals, and comparison of computed versus analytical CSDs, is sensitive to the number of crystals in the domain or domain size in two respects, which are investigated independently. The cataloguing of crystal number and lengths occurs for crystals that do not intersect the domain boundary. Therefore, the error in the total number of crystals is expected to be greater for a small domain than a large domain with crystals of equivalent size (i.e. differences in surface area to volume). A larger domain is also necessary to capture the randomness of crystal length, shape, location, and orientation. The former effect is investigated by dividing a single realization of a biminerallic simulation with α = 8, m = 96, nt = 15, α = 0.46875, and β = 1.5, discussed in detail later, into a number of cubic subdomains ranging in length from 8 to 48 voxels. The mean crystal length, based on converting crystal volume into an equivalent euhedral shape and determining the long axis length, is determined for each subdomain using only crystals that do not intersect the subdomain boundary. Ensemble-averages of the mean crystal length are determined for groups of subdomains of equal domain size. With decreasing domain size there is a significant decrease in the ensemble-average of mean crystal length relative to the mean crystal length of the original domain of 96 voxels in size (see Fig. 5). By extrapolating the monotonic curve to larger domain sizes and estimating an asymptotic value, it is likely that the mean crystal ¯ / L equal length for a domain of 96 voxels in length, or D

Fig. 5. Ensemble-averaged mean crystal length of groups of equally sized cubic subdomains taken from a single realization of a biminerallic simulation all normalized by the mean crystal length of the original domain. Error bars represent a one-sigma standard deviation in mean crystal length between groups of equal sized subdomains.

to 0.069, is underestimated by less than 10% due to the effect of ignoring crystals that intersect domain boundaries during cataloguing. Error associated with insufficient domain size due to not capturing the inherent randomness of the igneous microstructure is considered negligible if the coefficient of variation of mean crystal length, the ensemble standard deviation normalized by the ensemble-average, for a group of realizations approaches zero or is smaller than some arbitrary small number (such as 20%, Dagan, 1989; Zhang et al., 2000). The coefficient of variation of mean crystal length calculated for three realizations based on the biminerallic simulation mentioned above is 9.3%. The domain is therefore considered an REV (representative elementary volume) at least for the property of mean crystal length, provided ¯〉 / L is less than 0.073, where the angled brackets 〈D denote ensemble-average. As mentioned previously, memory limitations of the CPU limit the number of voxels in the domain and, therefore, limit the domain size if a reasonable crystal resolution is maintained. To relax error associated with insufficient domain size on the numbers and lengths of crystals, multiple realizations of the simulated microstructure are created with the same crystallization parameters and ensemble-averages are used. 7. Monominerallic crystallization As a simple first example of the stochastic algorithm, we simulate monominerallic crystallization of rectangular prisms of aspect ratio 1 : 1.5 : 3.5 in three realizations (r = 3) per simulation. The domain and crystallization times are discretized with the following parameters: m = 96, nt = 15, α = 0.46875, and β = 1.5. Three kinetic models are tested: constant nucleation and growth rate (a = 0), exponential nucleation and constant growth rate (a = 4 and a = 8), and exponential nucleation and dispersive growth rate (a = 8). Growth dispersion has been observed during crystallization experiments (e.g. Randolph and Larson, 1988) and is modeled by constant growth rate per individual crystal but with a varied growth rate over the bulk population of crystals. The growth rate is determined by, G˙ = G˙o10yσ, where σ is the standard deviation in the growth rate exponent, and y is a random deviate obeying a Gaussian distribution. All numerical results are averaged by the number of realizations (ensemble-averages) and only crystals that do not intersect the domain boundary at any time during the crystallization interval are considered. Crystal lengths are determined by converting the crystal volume in voxels to an equivalent euhedral shape and

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

measuring the long axis. Maximum length, taken as the maximum distance between any two vertices of an unshared voxel common to the crystal, increases nonmonotonically from 19% to 32% greater than the long axis length over the crystallization interval, however because the CSDs are qualitatively similar, only one crystal length metric, the long axis length, is considered further. The ensemble-averaged number of crystals for each simulation is shown compared to the analytical number of predicted crystals in Fig. 6 for the kinetic models involving a equal to 0, 4, and 8. The kinetic models with exponential nucleation recover more closely the analytical number of crystals than does the kinetic model with constant nucleation rate. The root-mean-square relative error between the numerical and analytical numbers of crystals is 1.53, 0.89, and 0.82 for kinetic models with, respectively, a equal to 0, 4, and 8. The reason for decreasing relative error with increasing degree of exponential nucleation is caused by the crystallization rule implemented in the algorithm. That is, at each time step all existing crystals must be grown before nuclei are formed. For some time steps, particularly in the simulations with small a, no nuclei are formed because all of the requisite crystallinity is taken up by growth of crystals. At subsequent time steps, however, the number of nuclei overshoots the predicted number. The effect of a discontinuous change in the number of crystals is minimized, as in the case for a = 8, if a sufficiently large change in crystallinity occurs between time steps. The constant growth rate model is not responsible for this effect, even though a change in unit length increases the unit volume by a factor of about eight, because this is taken into account in the derivation of the Avrami

Fig. 6. Ensemble-averaged number of crystals for monominerallic simulations with m = 96 and a equal to 0, 4, and 8.

41

Fig. 7. Ensemble-averaged crystal resolution of monominerallic simulations with a equal to 0, 4, and 8 as measured by ζ, the number of voxels in a crystal over the length of the crystal in voxels.

equation. Instead, this effect might derive from the Avrami formulation itself, which makes assumptions regarding removal of the effect of over-lapping crystal volumes when determining the crystallinity versus time curve. Therefore, there is reason to suspect that the results of the numerical model that reflects the spatial representation of crystallization may not necessarily identically match the theoretical model. An increase in ensemble-averaged crystal resolution as a function of crystallization time is shown in Fig. 7 for kinetic models with a equal to 0, 4, and 8. Crystal resolution is defined as the number of voxels in a crystal over the length of the crystal in voxels. The decrease in ensemble-averaged crystal resolution with increasing degree of exponential nucleation reflects the increase in number of crystals within the domain. The dimensionless mean crystal length of groups of crystals nucleated at the same time step is shown in Fig. 8 for a single simulation realization with a = 8. Simulations with different kinetic parameters yield a similar pattern and are not shown. The slope of each curve can be regarded as the instantaneous dimensionless mean growth rate for a group of crystals. All slopes should never exceed the input dimensionless growth rate (slope of the dotted line), slightly larger slopes, however, reflect the coarse mapping of crystals into the discretized domain and not an actual increase in the dimensionless growth rate. The mutual impingement of crystals terminates growth during the crystallization interval and all slopes eventually decrease to zero. The temporal evolution of the ensemble-averaged dimensionless CSDs for the simulation with a = 8 is plotted for certain time steps in Fig. 9. The schematic

42

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

Fig. 8. Dimensionless mean crystal length for groups of crystals nucleated at the same time step for a single monominerallic simulation realization with a = 8. The dotted line represents the hypothetical length of an unimpinged crystal.

diagram on the upper right shows how the CSD would evolve with time for a kinetic model involving exponential nucleation rate and constant growth rate. The numerical CSD is qualitatively similar to the theoretical CSD in that the slope remains approximately the same with time and the y-intercept increases with time. Local fluctuations in the slope of the numerical CSDs reflect the discontinuous production of nuclei during the crystallization interval. The ensemble-averaged dimensionless CSDs for simulations with a equal to 0, 4, and 8 at complete crystallinity are shown in Fig. 10. With increasing degree of exponential nucleation the CSDs are more linear and less concave down, have steeper negative slopes, and, if extrapolated, have increasing y-

Fig. 9. Temporal evolution of ensemble-averaged dimensionless CSD for monominerallic simulation with a = 8. Crystal length is made dimensionless using the maximum crystal length in the simulation.

Fig. 10. Ensemble-averaged dimensionless CSDs of monominerallic simulations at complete crystallinity with a equal to 0, 4, and 8. Error bars represent a one-sigma standard deviation in population density based on three realizations.

intercepts. They also show a fanning pattern which is both predicted and observed in igneous systems (Zieg and Marsh, 2002). The ensemble-averaged dimensionless CSDs at complete crystallinity for simulations involving dispersive crystal growth with σ equal to 0, 0.2, and 0.3 are shown in Fig. 11. The slopes of all CSDs are linear and approximately the same over most of the length dimension. At the largest lengths, however, the

Fig. 11. Ensemble-averaged dimensionless CSDs of monominerallic simulations at complete crystallinity with dispersive growth rate and σ equal to 0 (constant growth rate), 0.2, and 0.3. The two images are 2D slices through the 3D simulated microstructure showing the crystal fabric at complete crystallinity for σ equal to 0 (lower left corner) and 0.3 (upper right corner). Specific grayscale values for individual crystals are assigned arbitrarily for aid in visualization and do not represent relative crystal age. Error bars represent a one-sigma standard deviation in population density based on three realizations.

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

CSDs become more concave up with increasing degree of dispersive growth. Such a slight variation in the shape of the CSD with an increase in degree of dispersive growth is manifested more clearly by visual inspection of the resulting microstructure. An arbitrary 2D slice through a single 3D microstructure realization at complete crystallinity is shown in the lower left for σ equal to 0, and lower center for σ equal to 0.3 in Fig. 11. The microstructure with a high degree of dispersive growth (σ = 0.3) appears unnaturally porphyritic to the eye and therefore unrealistic in comparison to the simulated microstructure with no dispersive growth (σ = 0). Error bars in Figs. 10 and 11 represent a onesigma standard deviation in population density based on three realizations. The cumulative number of crystals is calculated based on subdividing the crystal length population after the final time step into ten evenly spaced bins. The population density is calculated by approximating the derivative as the slope of the line connecting two adjacent bin values and made dimensionless by normalizing to the volume of the domain and maximum crystal length.

43

8. Polyminerallic crystallization and percolation The stochastic algorithm can be easily extended to include crystallization of more than one mineral phase. It is common in most magmas for several distinct solid phases to grow simultaneously. Basalt crystallization is simulated with a simple model of simultaneous plagioclase and clinopyroxene (or olivine) crystallization in which both phases equally share the total crystallinity with time. Plagioclase is represented by a lath shape of aspect ratio, 1 : 2 : 5 and clinopyroxene by a stubby prism shape of aspect ratio, 1 : 1 : 2. The domain volume and crystallization time are discretized with the same parameters as in monominerallic crystallization and an exponential nucleation model of degree, a = 8, is used. A graphical display of the simulated microstructure in two and three dimensions is shown as Fig. 12. For brevity, plots of crystal numbers, sizes, and CSDs are not shown because they are very similar to the results of monominerallic crystallization. The ensemble-average of the final dimensionless mean crystal length, ¯〉 / L, is ∼1 / 14 and the ensemble-average of the final 〈D

Fig. 12. Microstructure simulation of biminerallic crystallization involving rectangular prisms of different but fixed aspect ratio and the same constant growth and exponential nucleation rate. This sequence is a simple model of basalt crystallization where the cool colors (blues) represent plagioclase and the hot colors (reds to yellows) represent clinopyroxene. Specific colors for individual crystals are assigned arbitrarily for aid in visualization and do not represent relative crystal age. The bottom eight figures are 2D slices represented by the plane indicated in the top right 3D microstructure by the bold dotted black lines.

44

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

of both melt and solid phases is tested using the results of a ‘mathematical opening’ procedure, shown schematically in two dimensions in Fig. 13. The procedure calculates the volume of the phase of interest in the domain, VS(d), that is swept out, or opened, by a sphere of given diameter, d (Serra, 1984; Hilpert et al., 2003). In this manner, percolation of VS(d) can be tested for a range of sphere diameters and a cumulative morphological phase-size distribution of the phase of interest (either melt or solid), F(d) = VS(d) / VP, can be constructed where VP is the volume of the phase of interest. The ensemble-averaged cumulative morphological phase-size distribution for the melt and solid phases is shown in Fig. 14 for crystallinities of 3.4% to 99.0%. Percolation of the melt and solid phases is

Fig. 13. ‘Mathematical opening’ of both melt and solid phases of a schematic 2D microstructure by two circles of different diameters shown on the left. Cumulative morphological phase-size distribution, F(d), is the ratio of the area swept across by a circle of diameter, d, shown in gray in the phase of interest divided by the total area of the phase of interest.

number of crystals not intersecting the domain boundary is 4255. For comparison, the ensembleaverage of the final dimensionless mean crystal length, ¯〉 / L, based on measuring crystal length as the 〈D maximum crystal dimension and not the equivalent long axis length is ∼1 / 10. As mentioned earlier, once crystallinity is spatially known in 3D, it becomes a valuable ‘specimen’ for further testing. An immediate application of the simulated microstructures, for example, is in determining the percolation threshold for either the melt or solid phase. The percolation threshold of a given phase is the critical volume fraction in which the phase of interest is connected or contiguous such that it extends across the full domain from one face to the opposite face. The percolation threshold of the solid network establishes the onset of finite yield strength in the partially molten magma and hence marks a change in bulk rheology from that of a Newtonian to a Bingham fluid (Saar et al., 2001). Conversely, the percolation threshold of the melt phase is a necessary but not sufficient condition for the onset of interstitial porous flow of melt (Cheadle et al., 2004). Percolation

Fig. 14. Ensemble-averaged cumulative phase-size distribution, 〈F(d)〉, for melt and solid phases over a range of crystallinities and sphere sizes. Cumulative phase-size distribution is the ratio of volume swept by a sphere in the phase of interest over the total volume of the phase of interest. The probability of percolation for a given phase, crystallinity, and sphere size is given by the pattern of the filled circular data points and is based on the percolation results of three realizations.

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

determined using a sphere with a diameter of one voxel and is given as a probability based on three realizations. Percolation is a scale-invariant property for randomly oriented objects therefore if percolation is achieved for one realization and not another realization at a given crystallinity, this suggests that the domain is not quite large enough to be an REV. Nevertheless, given this uncertainty, bounds on the percolation threshold can be estimated for both melt and solid phases. Percolation of the melt occurs at a crystallinity as high as 92.9% for a single realization and at greater than 78.0% for all three realizations. The percolation threshold for the melt phase around a crystallinity of 92.9% is comparable with a value of crystallinity around 92% for percolation threshold of the melt phase obtained by Cheadle et al. for a model of igneous microstructure evolution during crystallization but involving rectangular prisms of aspect ratio 1 : 1 : 3 (Cheadle et al., 2004). Percolation (or connectedness) of the solid network occurs at a crystallinity as low as 12.9% for a single realization and at less than 22.9% for all three realizations. The percolation threshold for the solid network around a crystallinity of 12.9% is consistent with the numerical modeling results for percolation of rectangular prisms of Saar et al. with the same crystal aspect ratios (Saar et al., 2001). Percolation of VS(d) over a range of spheres of diameter greater than one voxel in width could be useful in estimating certain physical properties of partially molten magmas such as electrical restivity or, in the case of melt percolation, for estimating sieve sizes for crystals settling or a secondary fluid phase either settling or rising (e.g. gas bubbles) through the microstructure as a function of crystallinity. Although the percolation of VS(d) for a given sphere diameter does not necessarily imply that a sphere of diameter d will be able to pass through the phase of interest, it is likely that the true diameter of a sphere able to pass is probably no less than half of d. As can be seen from the ensemble-averaged cumulative phase-size distribution of the melt in Fig. 14, percolation is achieved for all realizations at a maximum sphere diameter of 2, 3, 7, and 8 voxels or an equivalent dimensionless diameter, ¯〉, of 0.28, 0.43, 0.99, 1.46, at crystallinities of d / 〈D 78.0%, 57.5%, 37.9%, and 22.9%, respectively. That is, kinematically, spherical crystals or gas bubbles of ¯〉, could diameter no less than half of 0.28 times 〈D settle or rise completely through the igneous microstructure in a percolation channel at a crystallinity of 78%, provided no settling crystals or rising gas bubbles clog the percolation channel.

45

9. Conclusions A stochastic algorithm is developed to predict the 3D spatial representation of the igneous microstructure during progressive crystallization involving simultaneous nucleation and crystal growth. Of the kinetic models of crystallization tested, exponential nucleation rate with an exponential factor, a, of eight and constant crystal growth rate yields a microstructure in which the CSDs remain approximately log-linear (as often found for rocks) throughout the crystallization interval. For smaller values of a, the CSDs at complete crystallinity are increasingly concave down and therefore these kinetic models are considered to be less realistic in emulating Nature. Kinetic models involving increasing degrees of dispersive growth rate have minor but increasingly concave up CSDs at large crystal lengths. The effect of dispersive growth is much more diagnostic through visual inspection of the microstructure where the texture becomes markedly and unnaturally porphyritic. The kinetic models tested here are certainly simplified and abbreviated forms of the actual kinetic processes occurring during natural crystallization, but the reproduction of log-linear CSDs and reasonable microstructures attests to a first order realistic attempt at simulating natural microstructures. The onset of a finite yield strength in partially molten magma and the necessary condition for interstitial flow of melt within partially molten magma depends intimately on the percolation threshold for both the crystalline network and melt. The percolation threshold for the crystalline network is at a crystallinity less than ∼23% and for the melt is at a crystallinity greater than ∼78% based on the results of a ‘mathematical opening’ procedure of the simulated biminerallic microstructures. This ‘mathematical opening’ procedure can also be used to estimate the percolation threshold of spheres of varying size within either the crystalline network or melt. The implication for the melt phase is the determination of sieve sizes within the partially molten magma, which is a fundamentally important property for understanding sorting of crystals of varied size during development of layered igneous rocks (Marsh, 2004). Sieve sizes in the partially molten magma are also important for understanding upward percolation and stagnation of gas bubbles as a mobile secondary fluid phase through a crystalline network and could also be relevant to thermal rejuvenation of silicic magma bodies (Bachmann and Bergantz, 2003), chemical differentiation within compacting crystal piles (Boudreau, 2003), and possibly the explosivity of erupting volcanoes. An advantage of generating a simulated representation of the igneous microstructure is that, although serial

46

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47

sectioning and volume reconstruction (Marschallinger, 1998; Mock and Jerram, 2005) or microtomography (Ketcham and Carlson, 2001) are available to image the internal igneous microstructure of natural rocks, samples of partially molten magma arrested in a natural state are rare. Samples that fully capture the in situ spatial distribution of the different phases during crystallization that are free of any deformation due to the sampling process (i.e. crystal translation, reorientation, fracturing, or abundance of quench crystals) may be nearly impossible to obtain. Simulated igneous microstructures are invaluable as input into other physical models developed to obtain macro-scale physical properties of partially molten magma under conditions unobtainable by experiment. For example, permeability of simulated microstructures has been estimated using a networksimulation method by Cheadle et al. for porous melt flow at crystallinities greater than 70% (Cheadle et al., 2004). However, the lower limit of crystallinity for which porous melt flow in partially molten basalts is possible is near the percolation threshold of the crystalline network. In a companion paper, therefore, a lattice-Boltzmann model for grain-scale, single-phase flow is used to estimate permeability within a crystallinity range of 20% to 80% using the simulated igneous microstructures developed here (Hersum et al., 2005). Acknowledgements Reviews by R. Amenta and A. Mock are greatly appreciated. This work is supported by NSF Grant OPP0440718 to Johns Hopkins University. Appendix The stochastic algorithm is written in Fortran 90 and the simulations described in this paper were run on a single processor of a dual 1 GHz PowerPC G4 Apple Macintosh computer. A single realization of the polyminerallic crystallization simulations took 8.4 h of CPU time. The random number generator used in the algorithm is RAN1 of Numerical Recipes (Press et al., 1986). Post-processing of data and graphics were done using MATLAB. References Amenta, R.V., 2001. Three-dimensional computer modeling of fabric evolution in igneous rocks. Comput. Geosci. 27, 477–483. Amenta, R.V., 2004. Computer modeling of igneous hypidiomorphic textures using crystal prisms and plates with comparisons of measured crystal size distributions. Eos Trans. AGU, 85(17), Jt. Assem. Suppl., Abstract V43C-03.

Avrami, M., 1940. Kinetics of phase change: I. J. Chem. Phys. 8, 212–224. Bachmann, O., Bergantz, G.W., 2003. Rejuvenation of the Fish Canyon magma body: a window into the evolution of largevolume silicic magma systems. Geology 31, 789–792. Boudreau, A., 2003. Vapor separation and migration in a solidifying and compacting crystal pile. Eos Trans. AGU, 84(46), Fall Meet. Suppl., Abstract V51I-0390. Brandeis, G., Jaupart, C., 1987. The kinetics of nucleation and crystal growth and scaling laws for magmatic crystallization. Contrib. Mineral. Petrol. 96, 24–34. Burkhard, D.J., 2002. Kinetics of crystallization: example of microcrystallization in basaltic lava. Contrib. Mineral. Petrol. 142, 724–737. Cashman, K.V., 1993. Relationship between plagioclase crystallization and cooling rate in basaltic melts. Contrib. Mineral. Petrol. 113, 126–142. Cheadle, M.J., Elliott, M.T., McKenzie, D., 2004. Percolation threshold and permeability of crystallizing igneous rocks: the importance of textural equilibrium. Geology 32, 757–760. Dagan, D., 1989. Flow and Transport in Porous Formations. SpringerVerlag, 465 pp. Hall, J., 1798. Experiments on whinstone and lava. Trans. R. Soc. Edinb. 5, 43–76. Hersum, T.G., Hilpert, M., Marsh, B.D., 2005. Permeability and melt flow in simulated and natural partially molten basaltic magmas. Earth Planet. Sci. Lett. 237, 798–814. Hilpert, M., Glantz, R., Miller, C.T., 2003. Calibration of a porenetwork model by a pore-morphological analysis. Transp. Porous Media 51, 267–285. Jerram, D.A., Cheadle, M.J., Hunter, R.H., Elliott, M.T., 1996. The spatial distribution of grains and crystals in rocks. Contrib. Mineral. Petrol. 125, 60–74. Ketcham, R.A., Carlson, W.D., 2001. Acquisition, optimization and intrepretation of X-ray computed tomographic imagery: applications to the geosciences. Comput. Geosci. 27, 381–400. Kirkpatrick, R.J., 1981. Kinetics of crystallization of igneous rocks. In: Lasaga, A.C, Kirkpatrick, R.J. (Eds.), Kinetics of Geochemical Processes. Min. Soc. Am., 398 pp. Marschallinger, R., 1998. A method for three-dimensional reconstruction of macroscopic features in geological materials. Comput. Geosci. 24, 875–883. Marsh, B.D., 1981. On the crystallinity, probability of occurrence, and rheology of lava and magma. Contrib. Mineral. Petrol. 78, 85–98. Marsh, B.D., 1998. On the interpretation of crystal size distributions in magmatic systems. J. Petrol. 39, 553–599. Marsh, B.D., 2004. A magmatic mush column Rosetta Stone: the McMurdo Dry Valleys of Antarctica. Eos Trans. AGU 85 (47), 497–502. Mock, A., Jerram, D.A., 2005. Crystal size distributions (CSD) in three dimensions: insights from the 3D reconstruction of a highly porphyritic rhyolite. J. Petrol. doi:10.1093/petrology/egi024. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1986. Numerical Recipes. Cambridge University Press, 818 pp. Randolph, A.D., Larson, M.A., 1988. Theory of Particulate Processes. Academic Press, 369 pp. Saar, M.O., Manga, M., Cashman, K.V., Fremouw, S., 2001. Numerical models on the onset of yield strength in crystal-melt suspensions. Earth Planet. Sci. Lett. 187, 367–379. Serra, J., 1984. Image Analysis and Mathematical Morphology. Academic Press, 610 pp.

T.G. Hersum, B.D. Marsh / Journal of Volcanology and Geothermal Research 154 (2006) 34–47 Spohn, T., Hort, M., Fischer, H., 1988. Numerical simulation of the crystallization of multicomponent melts in thin dikes or sills 1. The liquidus phase. J. Geophys. Res. 93, 4880–4894. Tomaru, A., 2001. A numerical experiment of crystallization for a binary eutectic system with application to igneous textures. J. Geophys. Res. 106, 4037–4060. Wright, T.L., Okamura, R.T., 1977. Cooling and crystallization of tholeiitic basalt, 1965 Makaopuhi Lava Lake, Hawaii. Professional Paper 1004. U.S. Geological Survey.

47

Zhang, D., Zhang, R., Chen, S., Soll, W., 2000. Pore scale study of flow in porous media: scale dependency, REV, and statistical REV. Geophys. Res. Lett. 27, 1195–1198. Zieg, M.J., Marsh, B.D., 2002. Crystal size distributions and scaling laws in the quantification of igneous textures. J. Petrol. 43, 85–101.