| ログオフ・・・ゲスト様 |
コラム
CSD theory
以下は私が３０年以上前に書いたCSD理論に対するコメントの全てである．論文化されなかったため，著作権はすべて私にあるのでここに公開しておく．
Some comments on applying the CSD (Crystal Size Distribution) theory to magmatic rocks.
Abstract The CSD theory, which was originally developed in the field of chemical engineering, provides a useful way of relating the crystal size of rocks to crystallization kinetics. The first geologic application of the theory was made by Marsh (1988b) and Cashman and Marsh (1988); however, it contained some errors and problems. One issue is the solving method of a differential equation. Another concern is the treatment of system volume and the definition of crystal number density in relation to the volume problem. In this paper, the CSD theory was reconsidered according to the original method of Randolph and Larson (1971), and a correct way of applying the theory to geologic problems is discussed.
Introduction
Among several parameters that describe rock textures, crystal size is the most popular and best studied (Cashman, 1990). The crystal size is the time-integrated product of nucleation and growth of crystals; thus, it should bear important information about rock formation and crystallization processes. A full description of crystal size must be made in terms of size distribution or size frequency diagrams. In 1988, Marsh first introduced a theory that deals with such distribution quantitatively - the crystal size distribution (CSD) theory (Randolph and Larson, 1971), which was originally developed to deal with particulate processes in chemical engineering. The theory is based on considerations of population dynamics, i.e., the population balance of crystals of a particular size range during crystallization processes. It enables the observed CSD to be interpreted to yield quantitative information on dynamic parameters such as crystal growth rate, nucleation rate, and nucleation density, without knowing or making assumptions about the thermodynamic data of minerals or melts. The theory was very attractive because it appeared to provide a powerful tool to invert observations (CSD) into physical processes in geological problems. Indeed, Marsh (1988b) and Marsh and Cashman (1988) deduced many physical parameters from natural CSD by applying this method to igneous and metamorphic processes and demonstrated the usefulness of the method. Unfortunately, however, their application of the theory contained some serious errors, and thus some of the results they obtained had no geological significance and thus were misleading. Since Marsh (1988b) and Cashman and Marsh (1988), CSD was discussed by other authors with somewhat different approaches and different emphases (Maaloe et al., 1989; Armienti et al., 1994). but no major progress has been made since then in this field and the errors originally contained in Marsh (1988b) has not yet been corrected.
The purpose of this paper is to point out the errors of Marsh (1988b) and to propose a correct way of applying the CSD theory to geological problems. Because the CSD theory was formulated for engineering purposes, cares must be taken when applying it to geological problems, in which setting of chemical system is often arbitrary and thus must be clearly defined in order to avoid ambiguities and confusion. Marsh\textquotesingle s error partly stemmed from his erroneous definition of the system volume in his "batch crystallization ". He also started with a master equation, different from the original one of Randolph and Larson (1971), which became another source of error.
I first briefly review the basics of the theory of CSD, mostly following the way of Randolph and Larson (1971), with slight modification and expansion for geological application purposes. Then I analyze the method of Marsh (1988b) and point out the difficulty and errors in it. It is hoped that present review and correction help readers to sort out confusions in understating the CSD theory, that may have been present among readers, and promote sound studies in this field of petrology.
CSD theory
Size of minerals in rocks can be measured, if the shape of minerals is approximated to geometrical shape; circle, square, and so on. By measuring large number of crystals, we obtain accumulative crystal number $ N\left(L \right) $ - the total number of crystals of size $L$ and smaller in per unit volume. Function $N \left( L \right) $ is a monotony-increase function. Crystal population density $N \left( L \right) $ is defined, as
\[ n(L) = \frac{dN \left( L \right)}{dL} \tag{1} \]Total number of crystals between size $L_1$ and $L_2$ ($L_2 < L_1$) is therefore given as
\[ N \left( L_2 \right) - N \left( L_1 \right) = \int_{L_1}^{L_2}{n \left( L \right) dL}. \tag{2} \]Measured crystal population density for rocks represents that recorded at the time of solidification of magmas. In actual crystallization processes, each crystal migrates, being carried by flowing magma; therefore, crystal population density changes with time at fixed positions. In general, population density is a function of five variables: crystal size $L$, its position $x$, $y$, $z$, and time $t$. We consider an arbitrary small size range from $L$ to $L+\Delta L$ in a small cube, which occupies a region from $x$ to $x+\Delta x$, $y$ to $y+\Delta y$, and $z$ to $z+\Delta z$, for the $x$, $y$, and $z$ axes, respectively. Volume $\Delta V \left(=\Delta x \Delta y \Delta z \right)$ is a volume of the small cube. The population density $n \left( L \right)$ and $n \left( L+\Delta L \right)$ indicate those at size $L$ and $L+\Delta L$, respectively. Let the growth rate of crystals of size $L$ be $G \left( L \right)$ and that of size $ L+ \Delta L$ be $G \left( L+\Delta L \right)$. For an increment of time $\Delta t$, the number of crystals entering this size range due to the growth is $n \left( L \right) G \left( L \right) \Delta V \Delta t$. Similarly, the number of crystals leaving from this size range is $n \left( L+ \Delta L \right) G \left( L+\Delta L \right) \Delta V \Delta t$. The change of crystal number by growth is, therefore, $\left[ n \left( L \right) G \left( L \right) n \left( L+\Delta L \right) G \left( L+\Delta L \right) \right] \Delta V \Delta t$. The crystal flow rate per unit area in the $x$ direction at $x$ and $x+\Delta x$ are $v_x \left( x \right)$ and $v_x \left( x+\Delta x \right)$, respectively. The net flow rate for small area with dimensions $\Delta y$ and $\Delta z$ is $ \left[ v_x \left(x \right) n \left( x \right) - v_x \left( x + \Delta x \right) n \left( x + \Delta x \right) \right] \Delta y \Delta z \Delta L \Delta t$. The net flow rates along the $y$ and $z$ axes are similarly defined. The total changes of crystals number by growth and crystal migrations are balanced with the crystal number change from time $t$ to $t+\Delta t$, that is, $ \left[ n \left( t+ \Delta t \right) - n \left( t \right) \right] \Delta V \Delta L$.
\[ \begin{array}{lcl} \left\lbrack n(t + \Delta t) - n(t) \right\rbrack \Delta V \Delta L = & & \\ & & \left\lbrack v_x n \left( x \right) - v_x n \left( x + \Delta x \right) \right\rbrack \Delta y \Delta z \Delta L \Delta t \\ &+& \left\lbrack v_y n \left( y \right) - v_y n \left( y + \Delta y \right) \right\rbrack \Delta x \Delta z \Delta L \Delta t \\ &+& \left\lbrack v_z n \left( z \right) - v_z n \left( z + \Delta z \right) \right\rbrack \Delta x \Delta y \Delta L \Delta t \\ & +& \left\lbrack Gn \left(L \right) - Gn \left( L + \Delta L \right) \right\rbrack\Delta V\Delta t \\ \end{array} \tag{3} \]Dividing Equation (3) by $\Delta V$, $\Delta L$ and $\Delta t$, and setting these increments infinitesimally small, we obtain,
\[ \frac{\partial n}{\partial t} + \frac{\partial v_x n}{\partial x} + \frac{\partial v_y n}{\partial y} + \frac{\partial v_z n}{\partial z} + \frac{\partial Gn}{\partial L} = 0 \tag{4} \]where $v_x$, $v_y$, $v_z$, $G$ are defined as
\[ v_{x} = \frac{\partial x}{\partial t}, v_{y} = \frac{\partial y}{\partial t}, v_{z} = \frac{\partial z}{\partial t}, G = \frac{\partial L}{\partial t} \tag{5} \]$G$ is referred to the crystal growth rate. Equation (4) may be expressed using a vector notation
\[ \frac{\partial n}{\partial t} + \nabla \cdot \mathbf{v}n + \frac{\partial Gn}{\partial L} = 0 \tag{6} \]where, $\mathbf{v}$ is velocity vector for crystals. In order to consider crystal number change in a certain volume \emph{V}, we integrate Equation (6) with respect to volume,
\[ \int_{V}^{} \left( \frac{\partial n}{\partial t} + \nabla \cdot \mathbf{v} n + \frac{\partial Gn}{\partial L} \right) dV^{'} = 0 \tag{7} \]We assume that the crystal population density does not depend on the position within the system. In this case, care must be taken with the second term in equation (7). The second term can be transformed into a surface integration as
\[ \int_{V}^{}{\nabla \cdot \mathbf{v} n}dV^{'} = \int_{S}^{}{\mathbf{v}_{N}n}dS^{'} \tag{8} \] where $S$ is the surface of the system, $\textbf{v}_N$ is the normal velocity of crystals on the system surface. Equation (8) indicates the change in the number of crystals caused by crystals passing the system surface. Equation (8) is divided into two terms: volume change due to a moving surface and due to a flowing mixture. The volume change by the moving surface that moves to the crystal-free area is \[ n\int_{S_e}^{} \mathbf{v}_{e} dS^{'} = n\frac{dV}{dt} \tag{9} \] where $\mathbf{v}_e$ is the velocity of the moving surface, and $S_e$ is the free surface of the liquid. Crystals enter the system being carried by the flowing mixture. If we assume the population density of the mixture is $n_i$, and the volume flux of the flowing mixture into the system is $Q_i$, then the rate of increase of the number of crystals in the system is $n_i Q_i$. Similarly, we consider the number of crystals that leave the system with the flowing mixture. The volume flux of the flowing mixture leaving the system is assumed to be $Q_o$. Because the crystal population density in the mixture leaving the system is equal to the population density in the system, $n$, the number of crystals leaving the system is also represented as $n Q_o$. The change in the number of crystals due to the flowing mixture and the moving surface becomes, \[ \int_{S}^{}{\mathbf{v}_N n} dS^{'} = - n_{i}Q_{i} + n Q_{o} + n\frac{dV}{dt} \tag{10} \] substituting this into the equation (7) \[ \left( \frac{\partial n}{\partial t} + \frac{\partial Gn}{\partial t} \right) V + n\frac{dV}{dt} = n_{i}Q_{i} - nQ_{o} \tag{$11.a$} \] or \[ \begin{array} \frac{\partial n}{\partial t} + \frac{\partial Gn}{\partial L} &= - \frac{n}{V} \cdot \frac{dV}{dt} + \frac{\left( n_{i}Q_{i} - nQ_{o} \right)}{V} \\ & = - n\frac{d\left( \text{ln} V \right)}{dt} + \frac{\left( n_{i}Q_{i} - nQ_{o} \right)}{V} \end{array} \tag{11.b} \] where $Q_i$, $Q_o$, and $V$ should be related to each other as follows, \[ \frac{dV}{dt} = Q_{i} - Q_{o} \tag{12} \] In the following, we use Eqs. (11) and (12). Note that $V$ represents the volume of the mixture, which consists of liquid and solid phases. The region in which nucleation and crystal growth are actively taking place may be treated as a "liquid phase." Conversely, the "solid phase" may be defined as a region in which no nucleation or crystal growth occurs. Then volume $V$ is \[ V = V_{l} + V_{s}\ \tag{13} \] where $V_l$ and $V_s$ are the volumes of the liquid and solid, respectively. Even if the system volume $V$ in Eq. (11) is exchanged for the liquid or solid volume, Eq. (11) still holds. For example, let us express Eq. (11) in terms of the liquid volume $V_l$. The population density $n_l$ is defined to be calculated based on $V_l$. \[ V_n = V_{l}n_{l} \tag{14} \] This equation indicates that the population density for liquid phase becomes $V / V_l$ times larger than that for volume $V$. $Q_i^l$ is also defined to represent the volume of liquid in the inflowing mixture. $n_i^l$ is also defined to be calculated based on $Q_i^l$. \[ Q_{i}n_{i} = Q_{i}^{l}n_{i}^{l} \tag{15} \] Similarly, $Q_o^l$ gets \[ nQ_o = n_{l}Q_{l}^{o} = \frac{V}{V_{l}}nQ_{l}^{o} \tag{16.a} \] or \[ Q_o = \frac{V}{V_{l}}Q_{l}^{o} \tag{16.b} \] Substituting Eqs. (14), (15) and (16) into Eq. (11), we have \[ \left( \frac{\partial n_{l}}{\partial t} + \frac{\partial Gn_{l}}{\partial t} \right)V_{l} + n_{l}\frac{dV_{l}}{dt} = n_{i}^{l}Q_{i}^{l} - n_{l}Q_{o}^{l} \tag{17} \] This equation is also correct. Population density may be defined on the basis of either of the following three kinds of volumes: it is possible to use $V$, $V_l$ or $V_s$. If we use $V_l$ for a crystallization occurring system, we should take care of the nucleation density, because nucleation density exponentially increases with decreasing liquid volume. Combining Eq. (11) and (12), we get \[ \frac{\partial n}{\partial t} + \frac{\partial Gn}{\partial L} = \frac{\left( n_{i} - n \right)Q_{i}}{V} \tag{18} \] We consider a system in which $n$ does not depend on its position but only on crystal size and time. Eq. (18) can be solved under certain initial and boundary conditions. As an initial condition, the population density at time 0 is expressed as $f \left( L \right)$. There are two kinds of boundary conditions to be considered: one is concerned with the crystal flows across the system boundaries, $Q_i \left( t \right)$, $Q_o \left( t \right)$, and $n_i \left( L, t \right)$. The other is the crystal population density at size 0, $n^0$, which may be related to the crystal growth rate $G \left( t \right)$ and the nucleation rate $J \left( t \right)$ as \[ J \left(t \right) = \left. \frac{\partial n}{\partial t} \right|_{L = 0} = \left. \frac{\partial n}{\partial L} \cdot \frac{\partial L}{\partial t} \right|_{L = 0} = n \left(0,t \right) G \left( 0,t \right) \tag{19} \] or \[ n^{0} = \frac{J \left(t \right)}{G \left( 0,t \right)} \tag{20} \] where \( n^0 \) is the nucleation density, which is the population density at size 0. Eq. (20) indicates that the nucleation density is obtained from the nucleation rate divided by the growth rate. In the following, I will show some solutions to Eq. (18).1. The CSD at steady state
We consider a magma chamber the volume of which does not change with time ($dV / dt = 0$). Such a situation may be realized by the presence of a steady inflow and outflow of magma in the system. From Eq. (12), we have \[ Q_i = Q_o \tag{21} \] We assume that the crystal growth rate $G$ does not depend on crystal size (McCabe's $\Delta L$ law). We assume that the environment in the magma chamber does not change, so the crystal growth rate and nucleation rate do not change with time, and that the inflowing magma does not contain crystals ($n_i = 0$). With these conditions being considered, Eq. (18) becomes \[ \frac{\partial n}{\partial t} + G \frac{\partial n}{\partial L} = - n\frac{Q_o}{V} = - \frac{n}{\tau} \tag{22} \] where $\tau = Q_o / V$ is called the residence time. Eq. (22) can be solved using Laprace transformation as, \[ n = n^{0} \cdot \text{exp}\left( - \frac{L}{G\tau} \right) \left( L < Gt \right) \tag{23.a} \] \[ n = f\left( L - Gt \right) \cdot \text{exp}\left( - \frac{t}{\tau} \right) \left( L > Gt \right) \tag{23.b} \] where $f \left( L \right)$ is the initial crystal population density function. The CSD approaches to a steady distribution that is expressed by Eq. (23.a), in which the nucleation density $n^0$ is $J / G$, and the slope of exponential distribution is $ -1 / G \tau $.2. The CSD for a closed system
In a closed system, there is no inflow or outflow of magma (i.e., $Q_i = Q_o = 0$). If the system volume does not change during crystallization, and if the crystal growth rate $G$ does not depend upon grain size (McCabe's $\Delta L$ law), Eq. (18) becomes \[ \frac{\partial n}{\partial t} + G \left( t \right) \frac{\partial n}{\partial L} = 0 \tag{24} \] This equation is applied to the system in which magma is well-mixed and the population density is independent of its position. Such conditions may be satisfied in a lava lake in which crystal settling does not occur. Eq. (24) is similar to that which expresses the propagation of waves, suggesting that the CSD pattern moves with a velocity $G$ in the positive direction of the size axis. The general solution of Eq. (24) can be obtained, in which the growth rate $G$ and the nucleation rate $J $ are time - dependent as follows. We assume a CSD function $f \left ( L \right)$ ( $ L > 0$) at time 0, and we next consider its time development. If the crystal size is $ L_0 $ ($> 0$) at time 0, the size $L$ at time $t$ is, \[ L = L_{0} + \int_{0}^{t} {G \left(t \right) \text{dt}} = L_0 + L_G \tag{25} \] Here, $L_G$ is the amount of crystal growth from time 0 to $t$. Therefore, the CSD at time $t$ is written in terms of the initial CSD function $f \left( L \right)$ as \[ \begin{array}{lr} n \left( L \right) = f \left( L - L_G \right) & \left( L > L_G \right) \tag{26} \end{array} \] If the crystal size is smaller than $L_G$ at time $t$, this crystal must have nucleated and grown after the eruption. If this crystal nucleated at time $t_n$ ( $0 < t_n < t$), the size of this crystal is 0 at time $t_n$. Then, the crystal population density $n \left( L \right)$ becomes, using Eq. (20), \[ \begin{array}{lr} n \left( L \right) = \frac{J \left( t_n \right)}{G\left( t_n \right)} & \left( L < L_G \right) \tag{27} \end{array} \] The size $L$ and $t_n$ are related as \[ L = \int_{t_n}^{t} {G \left( t \right) dt} \tag{28} \] Eq. (27) indicates that the crystal population density is governed by the crystal growth rate and nucleation rate at time $t_n$ when the crystal is nucleated. Eqs. (26) and (27) are a general solution for Eq. (24) for an initial CSD, $f \left( L \right)$ (at time 0) and for a boundary condition of the nucleation rate $J \left( t \right)$ and the growth rate $G \left( t \right)$. In the simplest case in which the growth rate $G$ and nucleation rate $J$ are constant, being independent of time, Eqs. (26) and (27) become \[ n \left( L \right) = f \left( L - Gt \right) \left( L > Gt \right) \tag{29} \] \[ n \left( L \right) = \frac{J}{G} = \text{const.} \left( L < Gt \right) \tag{30} \] \textbf{An interpretation of the CSD theory by Marsh} The master equation that Marsh (1988b) used for the CSD was \[ \frac{\partial\left( Vn \right)}{\partial t} + \frac{\partial\left( GVn \right)}{\partial L} = Q_{i}n_{i} - Q_{o} \tag{31} \] Marsh (1988b) and Cashman and Marsh (1988) used a concept of batch crystallization and applied Eq. (31) to analyze the CSD of a Hawaiian lava lake. This model assumes no inflow or outflow of magma, therefore, $Q_i = Q_o = 0$, and that the volume of "actively crystallizing magma" decreases. If the crystal growth rate does not depend on the crystal size, Eq. (31) becomes \[ \frac{\partial\left( Vn \right)}{\partial t} + G\frac{\partial\left( Vn \right)}{\partial L} = 0 \tag{32} \] Marsh (1988b) solved the equation for the following initial condition. The initial CSD function $F \left( L \right)$, which is equal to $Vn$, at the time of eruption ($t = 0$) is defined as follows: \[ F \left(L \right) = V_0 n^0 \cdot \text{exp} \left( - \frac{L}{\tau G} \right) \tag{33} \] where $V_0$ is the volume of the actively crystallizing magma at time 0. The exponential distribution is the steady CSD that was derived in the previous section. If this steady CSD is already established in the magma chamber before magma eruption, it is appropriate to use this CSD as an initial condition for the consideration of the erupted magma. Knowing that the growth rate $G$ is constant, Marsh obtained the solution for Eq. (32) as \[ \begin{array}{lcl} V \left( t \right) n & = & F \left(L - Gt \right) \\ & = & V_0 n^0 \cdot \text{exp} \left( - \frac{L}{\tau G} + \frac{t}{\tau} \right) \end{array} \tag{34} \] The CSD therefore becomes \[ n = \left\lbrack \frac{V_0}{V \left(t \right)} \right\rbrack n^0 \cdot \text{exp} \left( - \frac{L}{\tau G} + \frac{t}{\tau} \right) \tag{35} \] We now discuss problems of Eq. (35).The boundary conditions for size axis
The CSD that is smaller than size $Gt$ at time t must have grown from nuclei that were formed at some time after the eruption. Therefore, Eq. (35) should be valid only for $L$ greater than $Gt$. The size $L$ must always be positive, but the initial condition of Marsh, i.e., Eq. (33), does not include this condition. An erroneous application of Eq. (35) to those crystals of size $Gt$ and smaller would imply that negative values were supposed for $L$ at $t = 0$ in Eq. (33), a function of the initial condition. Therefore, Eq. (35) expresses that the exponential distribution moves in the positive direction with a velocity $G$. By substituting Eq. (35) at size $L = 0$ into Eq. (20), it was concluded that the nucleation rate $J$ also increases exponentially with time. Marsh considered the growth rate to be constant but did not mention the nucleation rate as a boundary condition. Nevertheless, the nucleation rate was obtained uniquely. This erroneous result comes from their misapplication of Eq. (35) for size $L$ smaller than $Gt$.Discussion
We now discuss other problems included in other geologic applications of Marsh's method. Cashman and Marsh (1988) reported that the plagioclase CSD of the Kilauea lavas has an exponential distribution. They applied Marsh's theory, particularly Eq. (35) of steady state for lava lakes, and obtained the crystal growth rate and nucleation rate. They concluded that the growth rate of plagioclase is constant and that the nucleation density increases exponentially with time. They, however, solved Eq. (35) at a constant growth rate, and they did not consider the boundary condition at size 0, which determines both the growth rate and nucleation rate. Eq. (35) also indicates that the exponential CSD at the initial condition moves at a constant growth rate, and the nucleation rate increases exponentially. Their result for the lava lake is equal to the assumption made in Eq. (31). Moreover, Cashman and Marsh (1988) used Eq. (35), which was derived from the steady state CSD theory, and estimated the growth rate of plagioclase in the lava lake from the slope of the CSD. The assumption they made was that for a magma chamber at a steady state, which cannot be justified. Exponential distribution indicates that the magma chamber was at a steady state before the eruption and does not indicate the condition for the lava lake. The followings are remaining problems of CSD that have to be considered in further work.Summary
CSD theory was reconsidered according to Randolph and Larson's method. I also considered the rate of volume change and derived Eq. (18) as a master equation. Eq. (18) was solved for two cases: one for a steady state in an open system and another for a closed system. For the former case, CSD converges to exponential distributions. In this case, the slope of the logarithmic CSD is related to the growth rate and the residence time. For the latter case, CSD is expressed by the equation of wave propagation. The CSD, therefore, is determined by the nucleation and growth rate after magma eruption. These two cases may be applicable to a magma chamber and a lava lake, respectively. The problems with Marsh (1988b) and Cashman and Marsh (1988) may be summarized as follows: (1) Volume $V$ in Marsh (1988b) only includes the "actively crystallizing magma volume." It was interpreted that the increase in population density in Marsh's equation was mostly caused by a decrease in the liquid volume. Population density should be calculated for system volume including liquid and solid phases, so that a change in population density is not caused by a change in volume. (2) Marsh (1988b) did not consider conditions at size 0. In his method, the number of crystals whose sizes are smaller than $Gt$ is determined by reflecting the initial CSD. The CSD in this size range should be obtained by using the crystal growth rate and nucleation rate at size 0 during time $t$. (3) Cashman and Marsh (1988) concluded from their study of plagioclase in Hawaiian lava lakes that the growth rate was constant and the nucleation rate increased exponentially with time. This result only supports the assumption of Marsh (1988b). Cashman and Marsh (1988) also calculated the growth rate of plagioclase from the slope of the CSD. However, this method is only possible for an open magma chamber in a steady state, and not for a closed-system lava lake.