A NEW ESTIMATE OF THE CHONDRULE COOLING RATE DEDUCED FROM AN ANALYSIS OF COMPOSITIONAL ZONING OF RELICT OLIVINE

We simulate the olivine growth in a manner consistent with a thermodynamical model of olivine. We consider a simple binary system for olivine composed of two crystalline end members, forsterite (Mg₂SiO₄) and fayalite (Fe₂SiO₄). Although natural ferromagnesian olivines contain minor amounts of additional divalent cations including Ca, Mn, Ni, and Co (Hirschmann 1991), we ignore these minor elements as a first step. The fayalite contents in solid and liquid in equilibrium are given by an ideal solution model in good approximation (Bowen & Schairer 1935). The concentration c is defined by an atomic ratio Fe/(Mg + Fe). The equilibrium concentrations in liquid and solid, $c_{\rm L}^{\rm e}(T)$ and $c_{\rm S}^{\rm e}(T)$, at temperature T are given by

$$\begin{equation} c_{\rm L}^{\rm e} (T) = \frac{e_{\rm A}(T) - 1} {e_{\rm A}(T)-e_{\rm B}(T)}, \quad c_{\rm S}^{\rm e} (T) = e_{\rm B}(T) \, c_{\rm L}^{\rm e} (T), \end{equation} \tag{ 1 }$$

where

$$\begin{equation} e_{\rm A,B} (T) = \exp \left[ \frac{L_{\rm A,B}}{R} \left(\frac{1}{T} - \frac{1}{T_{\rm A,B}} \right) \right]. \end{equation} \tag{ 2 }$$

Here, subscripts A and B stand for pure forsterite (Mg₂SiO₄) and pure fayalite (Fe₂SiO₄), respectively; L_{A, B}, their molar latent heats of crystallization; T_{A, B}, their absolute melting points; and R, the gas constant. The parameters are set as L_A = L_B = 58.6 kJ mol⁻¹ (Bowen & Schairer 1935), T_A = 2163 K, and T_B = 1490 K. The phase diagram is shown in Figure 2. Values of the equilibrium concentrations $c_{\rm L}^{\rm e}$ and $c_{\rm S}^{\rm e}$ are listed in Table 1. Note that the equilibrium partitioning coefficient, $k_0 \equiv c_{\rm S}^{\rm e} / c_{\rm L}^{\rm e} = e_{\rm B}$, is less than unity for all T.

Zoom In Zoom Out Reset image size

Table 1. Parameter Values Deduced from the Phase Diagram for Various Temperatures T

T	$c_{\rm L}^{\rm e}$	$c_{\rm S}^{\rm e}$	k0	ΔTc
(K)				(K)
1500	0.990	0.960	0.969	7.6
1600	0.886	0.640	0.722	80
1700	0.764	0.426	0.557	144
1800	0.625	0.277	0.443	194
1900	0.471	0.170	0.360	224
2000	0.303	0.091	0.299	220
2100	0.121	0.031	0.253	142

Notes. Here, $c_{\rm L}^{\rm e}$ and $c_{\rm S}^{\rm e}$ are equilibrium concentrations in liquid and solid, respectively; $k_0 \equiv c_{\rm S}^{\rm e} / c_{\rm L}^{\rm e}$, equilibrium partitioning coefficient; ΔT_c, temperature difference between liquidus and solidus for given concentrations $c = c_{\rm L}^{\rm e}$.

Download table as:  ASCIITypeset image

Crystal growth associated with solute diffusion in a liquid has been formulated as a moving boundary problem to describe the compositional zoning (Tiller et al. 1953; Smith et al. 1955). Let us consider a crystal growing toward x direction. We also take another coordinate x' comoving with the solid–liquid interface; the origin x' = 0 is fixed at the interface. We neglect the solute diffusion in solid (x' < 0) because the diffusivity in the solid is much smaller than that in the liquid (Dohmen et al. 2007). The diffusion equation for the concentration c_L(x'; t) in the liquid is given by

$$\begin{equation} \frac{\partial c_{\rm L}}{\partial t} = \frac{\partial }{\partial x^{\prime }} \left(D_{\rm L} \frac{\partial c_{\rm L}}{\partial x^{\prime }} \right) + V \frac{\partial c_{\rm L}}{\partial x^{\prime }} \quad (x^{\prime } > 0), \end{equation} \tag{ 3 }$$

where V(t) is the time-dependent growth velocity and D_L is the diffusivity in the liquid. We neglect convective solute transportation in the liquid in Equation (3). A coordinate fixed with the solid, x, is then

$$\begin{equation} x = x^{\prime } + \int _0^t V(t^{\prime }) dt^{\prime } = x^{\prime } + d(t), \end{equation} \tag{ 4 }$$

where d(t) is the position of the interface on the coordinate x.

Equation (3) must be solved subject to the initial condition given by

$$\begin{equation} c_{\rm L}(x^{\prime };0) = c_{\rm L0}, \end{equation} \tag{ 5 }$$

the boundary conditions given by

$$\begin{equation} c_{\rm L}(\infty; t) = c_{\rm L0} \end{equation} \tag{ 6 }$$

at x = ∞, and

$$\begin{equation} \frac{\partial c_{\rm L}}{\partial x} + \frac{V}{D_{\rm L}} (c_{\rm Li} - c_{\rm Si}) = 0 \end{equation} \tag{ 7 }$$

at x' = 0 for any t (Lasaga 1982), where c_L0 is the initial concentration in the liquid and c_Li and c_Si are the concentrations at the interface on the liquid and solid sides, respectively. Tiller et al. (1953) and Smith et al. (1955) set c_Si = k₀c_Li in Equation (7), where k₀ is a constant partitioning coefficient. However, k₀ depends on the temperature T in general and varies with T in the course of cooling. Furthermore, they set V to be constant. But V actually varies with T and this determines the zoning profile we are interested in. Thus, we must determine the time variation of V as well as the concentration c_L. In the narrow regions on both sides of the solid–liquid interface, we can expect equilibrium partitioning to be realized. Hence, we assume local equilibrium partitioning, namely,

$$\begin{equation} c_{\rm Li} = c_{\rm L}^{\rm e}(T) \quad {\rm and} \; c_{\rm Si} = c_{\rm S}^{\rm e}(T). \end{equation} \tag{ 8 }$$

With the use of Equation (8), the boundary condition (7) is written as

$$\begin{equation} \frac{\partial c_{\rm L}}{\partial x} + \frac{V}{D_{\rm L}} \left[ c_{\rm L}^{\rm e} (T) - c_{\rm S}^{\rm e} (T) \right] = 0 \quad {\rm at} \; x^{\prime }=0 \; {\rm for \; all}\; t. \end{equation} \tag{ 9 }$$

We consider the spatial domain $0 < x^{\prime } < x^{\prime }_{\infty }$ for the numerical simulation, where $x^{\prime }_{\infty }$ is set to be so large that $c_{\rm L}(x^{\prime }_\infty;t)$ is not affected by the diffusive boundary layer. The domain is divided into N grid cells of uniform width, $\Delta x^{\prime } = x^{\prime }_{\infty }/N$. The time direction is also discretized by a constant interval Δt. $c_i^{(n)}$ denotes the concentration at the grid cell i at the time step n. From Equation (3), we obtain the difference equation:

$$\begin{eqnarray} && c_i^{(n+1)} = \nonumber\\ && \left\lbrace \begin{array}{@{} ll} c_i^{(n)} + \eta \left(c_{i+1}^{(n)} - c_i^{(n)} \right) + v \left(c_{\rm L}^{\rm e}-c_{\rm S}^{\rm e} \right) + v \left(c_{i+1}^{(n)} - c_i^{(n)} \right) &\! {\rm for} \; i=1, \\ c_i^{(n)} + \eta \left(c_{i+1}^{(n)} - 2 c_i^{(n)} + c_{i-1}^{(n)} \right) + v \left(c_{i+1}^{(n)} - c_i^{(n)} \right) & \! {\rm for \; others,} \end{array} \right. \nonumber\\ \end{eqnarray} \tag{ 10 }$$

where η = D_LΔt/Δx'² and v = VΔt/Δx'. The term $v(c_{\rm L}^{\rm e} - c_{\rm S}^{\rm e})$ denotes the boundary condition given by Equation (9). For i = N, we use $c_{N+1}^{(n)} = c_{\rm L0}$. We use the first-order upwind differential scheme for the advection term under the assumption of v > 0. Note that v is not given a priori in our model but is determined to satisfy a local equilibrium condition at x' = 0:

$$\begin{equation} c_1^{(n+1)} = c_{\rm L}^{\rm e} (T^{(n+1)}). \end{equation} \tag{ 11 }$$

We need an iterative scheme to obtain v as described below. $\tilde{v}^{(k)}$ denotes a trial value of v at kth iterative step. We obtain the trial solution of concentration $\tilde{c}_i^{(k)}$ by

$$\begin{eqnarray} {{ \tilde{c}_i^{(k)} = \left\lbrace \begin{array}{@{} ll} c_i^{(n)} + \eta \left(c_{i+1}^{(n)} - c_i^{(n)} \right) + \tilde{v}^{(k)} \left(c_{\rm L}^{\rm e}-c_{\rm S}^{\rm e} \right) + \tilde{v}^{(k)} (c_{i+1}^{(n)} - c_i^{(n)}) & {\rm for}\; i=1, \\ c_i^{(n)} + \eta \left(c_{i+1}^{(n)} - 2 c_i^{(n)} + c_{i-1}^{(n)} \right) + \tilde{v}^{(k)} \left(c_{i+1}^{(n)} - c_i^{(n)} \right) & {\rm for \; others.} \end{array}\right. }} \nonumber\\ \end{eqnarray} \tag{ 12 }$$

If $\tilde{c}_1^{(k)}$ satisfies the condition given by Equation (11), we can proceed to the next time step with $c_i^{(n+1)} = \tilde{c}_i^{(k)}$. If $\tilde{c}_1^{(k)} > c_{\rm Li}$, which suggests that large $\tilde{v}^{(k)}$ causes over-partitioning, we should choose $\tilde{v}^{(k+1)}$ smaller than $\tilde{v}^{(k)}$. In contrast, if $\tilde{c}_1^{(k)} < c_{\rm Li}$, we should choose $\tilde{v}^{(k+1)}$ larger than $\tilde{v}^{(k)}$. We use

$$\begin{equation} \tilde{v}^{(k+1)} = \frac{c_{\rm L}^{\rm e}}{\tilde{c}_1^{(k)}} \tilde{v}^{(k)} \end{equation} \tag{ 13 }$$

for the next trial. The condition for convergence is taken to be

$$\begin{equation} \left| \tilde{c}_1^{(k)} - c_{\rm L}^{\rm e} \right| < 10^{-6}. \end{equation} \tag{ 14 }$$

The change in the threshold value does not affect the solution significantly. In this paper, we use $x^{\prime }_{\infty } = 100 \ {\rm \mu m}$, N = 10⁴, Δx' = 10 nm, and Δt = 10⁻⁶ s.

We consider a solid–liquid equilibrium at T = T₀ at the beginning. The initial liquid concentration is given by $c_{\rm L0} = c_{\rm L}^{\rm e}(T_0)$. We decrease T at a constant rate R_c and continue the numerical procedure until the zoning profile in the solid forms obviously. We carry out the numerical simulation for different sets of parameters T₀ and R_c to investigate how the zoning profile depends on these parameters. For simplicity, we use a constant diffusivity of D_L = 10⁻⁵ cm² s⁻¹, which is a typical value for silicate melt (Brady 1995).

Figure 3 shows a numerical result for T₀ = 1900 K and R_c = 10³ K s⁻¹. Panel (a) shows the concentrations in solid and liquid phases. The dashed, dotted, or dash–dotted curve shows a snapshot of concentration c_L(x; t) in the liquid phase at every 50 K decrease in temperature. The horizontal axis indicates x, not x', so the solid–liquid interface is not fixed at the origin. As the temperature decreases, the interface moves rightward. Figure 3(a) shows that c_L rises ahead of the moving interface as a result of the partitioning. The diffusive boundary layer appears when the solute diffusion is insufficient to homogenize the concentration in the liquid. The appearance of the diffusive layer plays an important role in determining growth velocity and a zoning profile in the solid.

Zoom In Zoom Out Reset image size

The solid curve shows a series of c_Si formed as a result of the motion of the interface. This is the zoning profile recorded in the solid immediately after its formation; in fact, a solid phase has non-zero diffusivity, so the zoning profile will be modified thermally (Miyamoto et al. 1986). In this study, we neglect the solute diffusion in the solid. Note that c_Si increases monotonically with x and exceeds the initial liquid concentration c_L0 = 0.471 at x = 16.5 μm. The region in which the composition changes steeply is called an initial transient in alloy solidification (Smith et al. 1955). Here, we define the region in which c_S < c_L0 as the initial transient; its width is d_f = 16.5 μm in this case.

Figure 3(b) plots c_Si and c_Li on the phase diagram during solidification. Note that both c_Si and c_Li vary respectively along the solidus and liquidus curves with the decrease in temperature. This confirms that the local equilibrium condition is satisfied in the calculation.

Figure 4 shows the interface position d(t) for T₀ = 1900 K and R_c = 10³ K s⁻¹ as a function of time, indicating that the crystal growth is accelerated with time. Let us assume that the growth velocity $V = \dot{d}$ is expressed by

$$\begin{equation} V(t) = V_0 \, e^{t/t_0}, \end{equation} \tag{ 15 }$$

where V₀ is the initial growth velocity and t₀ is a time constant. The physical grounds of Equation (15) will be discussed in Section 4.1. Integration of Equation (15) leads to

$$\begin{equation} d(t) = V_0 t_0 \, (e^{t/t_0} - 1), \end{equation} \tag{ 16 }$$

where we put d(0) = 0. We fit the numerical result using Equation (16) and obtain best-fit values of V₀ = 26.2 μm s⁻¹ and t₀ = 0.124 s for T₀ = 1900 K and R_c = 10³ K s⁻¹. Note that Equation (16) reproduces the numerical result excellently in almost the entire region of the initial transient.

Zoom In Zoom Out Reset image size

We carry out calculations for a wide range of T₀ and R_c. Numerical results are summarized in Table 2. Figure 5 shows dependence of the initial transient parameters of d_f (panel a), V₀ (panel b), and t₀ (panel c) on R_c. Results for different T₀ values are shown by using different symbols (see the legend in panel a). We find that all these parameters show a power-law dependence on R_c, namely, $d_{\rm f} \propto R_{\rm c}^{-1/2}$, $V_0 \propto R_{\rm c}^{1/2}$, and $t_0 \propto R_{\rm c}^{-1}$. The power-law dependence is observed for a wide range of T₀.

Zoom In Zoom Out Reset image size

The numerical simulation yielded a remarkable result that the interface velocity increases exponentially with time for a constant rate of temperature decrease. We investigate its physics briefly in this subsection. Figure 6(a) shows the motion of the solid–liquid interface followed by the change in the concentration profile during a short duration Δt. The black and red curves are the profiles before and after the motion of the interface, respectively. The concentration excess (see Appendix A.1) in the liquid is evaluated to be

$$\begin{equation} A_{\rm L}(t) = \int _{x_{\rm i}}^\infty [ c_{\rm L}(x;t) - c_{\rm L0}]\, dx^{\prime } \sim \left[ c_{\rm L}^{\rm e}(T(t)) - c_{\rm L0} \right] \delta, \end{equation} \tag{ 17 }$$

where δ is a width of the diffusive boundary layer and is approximately constant. A_L increases monotonically with decreasing T because of the ejection of a portion of the solute from the solidified layer. On the other hand, the amount of solute ejected from the solidified layer during Δt is given (see Appendix A.1) by

$$\begin{equation} \Delta A_{\rm S}(t) \simeq \left[ c_{\rm L0} - c_{\rm S}^{\rm e}(T(t)) \right] V \Delta t. \end{equation} \tag{ 18 }$$

The solute conservation requires ΔA_S(t) = A_L(t + Δt) − A_L(t), hence, we obtain

$$\begin{equation} V \sim \frac{ R_{\rm c} \delta }{ k_0 } \cdot \frac{ \left| dc_{\rm S}^{\rm e}/dT \right| }{ c_{\rm L0} - c_{\rm S}^{\rm e}(T(t)) }, \end{equation} \tag{ 19 }$$

where R_c = −dT/dt and $c_{\rm S}^{\rm e} = k_0 c_{\rm L}^{\rm e}$. Here, k₀ can be regarded as a constant around T ∼ T₀. Since solidus can be approximated by a straight line around T ∼ T₀ (see Figure 6(b)), we obtain

$$\begin{equation} {-} \frac{dc_{\rm S}^{\rm e}}{dT} \simeq \frac{(1-k_0) c_{\rm L0}}{\Delta T_c}, \end{equation} \tag{ 20 }$$

where ΔT_c is the temperature difference between liquidus and solidus at a given concentration. Using Equation (20), the time variation of V is evaluated to be

$$\begin{equation} \frac{\dot{V}}{V} \sim \frac{ R_{\rm c} }{\Delta T_c} \cdot \frac{ (1-k_0) c_{\rm L0}}{c_{\rm L0} - c_{\rm S}^{\rm e}(T)} \sim \frac{R_{\rm c}}{\Delta T_c}, \end{equation} \tag{ 21 }$$

where $c_{\rm S}^{\rm e}(T) \sim k_0 c_{\rm L0}$ around T ∼ T₀. Equation (21) indicates that V increases exponentially with t at least in the early stage of the growth, which is in harmony with Equation (15). It is also indicated that the time constant t₀ is on the order of ΔT_c/R_c.

Zoom In Zoom Out Reset image size

We showed t₀ ∼ ΔT_c/R_c ≡ t_sol, where t_sol is called a solidification time, that is, a period required for the temperature to decrease from the liquidus to the solidus for a given c_L0. In Figure 7(a), we compare t_sol with t₀ calculated for all input parameters. Each symbol corresponds to a result for a set of T₀ and R_c. It is found that all symbols are aligned on a straight line with a slope of unity, suggesting that t₀ is proportional to t_sol for wide range of T₀ and R_c. The symbols align along the t₀ = t_sol/2 line. Hence, including the numerical factor, t₀ is expressed by

$$\begin{equation} t_0 \approx \frac{1}{2} \frac{\Delta T_c}{R_{\rm c}}. \end{equation} \tag{ 22 }$$

Note that t₀ is inversely proportional to R_c.

Zoom In Zoom Out Reset image size

Taking the exponential acceleration of the growth velocity into account, the solute distribution in the crystal in the initial transient is expressed (see the Appendix) by

$$\begin{equation} c_{\rm S}(x) \approx k_0 c_{\rm L0} \left[ 1 + (1-k_0) \sqrt{ \frac{2 \pi R_{\rm c}}{D_{\rm L} \Delta T_c} } \cdot x \right], \end{equation} \tag{ 23 }$$

where we used Equation (22). Equation (23) shows that the solute concentration increases linearly with distance x in contrast to the Rayleigh fractionation (Rayleigh 1896) and Smith's model (Smith et al. 1955). From c_S(d_f) = c_L0, we immediately obtain the width of the initial transient expressed by

$$\begin{equation} d_{\rm f} \approx \frac{1}{k_0} \sqrt{ \frac{ D_{\rm L} \Delta T_c }{ 2 \pi R_{\rm c} } }. \end{equation} \tag{ 24 }$$

Note that d_f is proportional to $R_{\rm c}^{-1/2}$ in agreement with the numerical results shown in Figure 5(a). Figure 7(b) shows that Equation (24) well reproduces the numerical results.

As an example of the application of the present theory, we estimate a cooling rate by analyzing a zoning profile recorded in chondrules. Wasson & Rubin (2003) observed zoning profiles of olivines in type II chondrules in Y-81020 carbonaceous chondrite. They found that some large (>40 μm) olivine phenocrysts have low-FeO relicts near the centers and surrounding FeO-rich margins. In the margin, the fayalite composition linearly increases from c ≃ 0.22–0.39 to ≃ 0.38–0.77 within the thickness of ≃ 2–12 μm. In other words, the compositional gradient is ∼10²–10³ cm⁻¹. Wasson (2004) assumed that the narrow FeO-rich margins were produced by diffusion-controlled overgrowth of the phenocrysts during the final melting event and estimated the cooling rate to be ∼10³ K s⁻¹. However, his estimation includes an uncertain parameter of "blocking" temperature at which the overgrowth stops; it must depend on the cooling rate. Furthermore, he did not explain the origin of the linear zoning profile.

Our model naturally explains the linear zoning profile (see Figure 3(a)) and can estimate the cooling rate without introducing an uncertain parameter. From Equation (23), we obtain

$$\begin{equation} R_{\rm c} \approx \frac{D_{\rm L} \Delta T_c}{2 \pi k_0^2 (1-k_0)^2 c_{\rm L0}^2} \left(\frac{dc_{\rm S}}{dx} \right)^2, \end{equation} \tag{ 25 }$$

where dc_S/dx is the compositional gradient in the linear zoning profile. From this, the cooling rate is estimated to be R_c ≈ 200–2000 K s⁻¹ for D_L = 10⁻⁶–10⁻⁵ cm² s⁻¹, ΔT_c = 200 K, k₀ = 0.4, c_L0 = 0.5, and dc_S/dx = 300 cm⁻¹. The R_c value is on the same order of magnitude as that estimated by Wasson (2004) incidentally, although our theory clearly identifies the origin of the linear zoning profile. These estimates are higher by orders of magnitude than the cooling rates of 0.01–1 K s⁻¹ inferred from furnace-faced experiments (Jones et al. 2000; Hewins et al. 2005).

Aerodynamic drag heating induced by nebula shocks is one of the plausible mechanisms to melt chondrule precursor silicate grains (Wood 1984). This mechanism accounts for various properties of chondrules such as appropriate cooling rates inferred by furnace-based experiments (Hood & Horanyi 1991; Hood & Ciesla 2001; Desch & Connolly 2002; Miura & Nakamoto 2006; Morris et al. 2009; Morris & Desch 2010; Morris et al. 2012), maximum temperature up to the melting point (Iida et al. 2001), size distribution (Susa & Nakamoto 2002; Miura & Nakamoto 2005; Kato et al. 2006; Miura & Nakamoto 2007), oblate-shaped chondrules (Miura et al. 2008a), and formation of compound chondrules (Ciesla 2006; Miura et al. 2008b; Yasuda et al. 2009).

The aerodynamic heating requires the chondrule-forming region to be optically thick for the dust thermal radiation to satisfy the cooling rate constraint inferred by the furnace-based experiments (≲ 10³ K hr⁻¹); inefficient radiative losses in the dusty cloud result in a much slower cooling rate than that expected from the radiative cooling in vacuum (Hood & Ciesla 2001; Desch & Connolly 2002). However, our new constraint given by Equation (25) suggests an extremely rapid cooling as fast as that expected from the free radiation. The extremely rapid cooling supports shock waves produced in the optically thin environment: planetesimal bow shocks (Hood 1998; Ciesla et al. 2004; Hood et al. 2005), nebula shocks in less-dusty environments (Iida et al. 2001), and shock waves induced by X-ray flares (Nakamoto et al. 2005). It has been believed that most of chondrules were melted multiple times in the early solar nebula (Jones et al. 2000) and the FeO-rich overgrowth on the relict porphyritic olivines records the final melting event in the multiple heating episodes (Wasson & Rubin 2003). The extremely rapid cooling on chondrules at the final melting event is consistent with the fact that the solar nebula gradually became optically thin due to accumulation of dust particles onto the forming planets.

Lightning discharges are one of the leading candidates satisfying the extremely rapid cooling constraint. Most models of the nebular lightning focused on the mechanism of charge separation and initiation of electrical breakdown (Gibbard & Levy 1997; Pilipp et al. 1998; Desch & Cuzzi 2000). On the other hand, Horányi et al. (1995) examined the thermal history of chondrule precursor dust particles in the optically thin expanding discharge channel quantitatively. According to Horányi's model, the precursor particle heated above its melting point solidifies within about 1 s, which corresponds to a cooling rate of a few 100 K s⁻¹ and meets the constraint given by Equation (25).

Consider a system composed of solid and liquid phases in equilibrium. Figure 8(a) shows solute concentration distributions in solid and liquid phases during solidification. Initially, the region of x < 0 was occupied by the solid phase and x > 0 by the liquid phase. The liquid was homogeneous with a constant concentration c_L0. As the solidification proceeds, the interface moves rightward and a portion of the solute is ejected from the solidified region and added to the adjacent liquid region. The solute added to the liquid is transported away from the interface afterward by diffusion. When the interface reaches x = x_i, the total amount of solute A_S ejected from the solidified region is given by

$$\begin{equation} A_{\rm S} = \int _0^{x_{\rm i}} (c_{\rm L0} - c_{\rm S}(x)) dx. \end{equation} \tag{ A1 }$$

This must be equal to that added to the liquid given by

$$\begin{equation} A_{\rm L} = \int _{x_{\rm i}}^{\infty } (c_{\rm L}(x) - c_{\rm L0})\, dx \end{equation} \tag{ A2 }$$

from the conservation of solute content. Smith et al. (1955) derived analytic solutions of c_S(x) and c_L(x) assuming that the interface moves at a constant speed.

Zoom In Zoom Out Reset image size

We derive a solution c_S(x) in the initial transient for the case that the velocity changes with time. Let us consider an infinitesimal solute element dA_S ejected from the solidified region during the movement of the interface from x = x'' to x = x'' + dx'' (see Figure 8(b)). The solute element dA_S transferred to the liquid diffuses toward x > x_i and has a distribution dc_L(x) in the liquid as a result. From the mass conservation, dA_S must be equal to

$$\begin{equation} dA_{\rm L} = \int _{x_{\rm i}}^{\infty } dc_{\rm L} \, (x)\, dx. \end{equation} \tag{ A3 }$$

The distribution dc_L(x) has a typical diffusion length x_D given by $x_D = \sqrt{D_{\rm L} t^{\prime \prime }}$, where t'' is the time elapsed after the ejection; the solute element ejected earlier has longer t''. By determining dA_S, dA_L, dc_L(x), x_D, and t'' properly, we obtain an expression for c_S(x).

We assume that dc_L is given by

$$\begin{equation} dc_{\rm L}(x^{\prime }) = B \, e^{- x^{\prime } / \sqrt{D_{\rm L} t^{\prime \prime }} }, \end{equation} \tag{ A4 }$$

where B is a constant and x' = x − x_i is the distance measured from the interface to the point concerned in the liquid. B is determined from the mass conservation dA_S = dA_L. From Equation (A4), we obtain

$$\begin{equation} dA_{\rm L} = \int _0^{\infty } dc_{\rm L}(x^{\prime })\, dx^{\prime } = B \sqrt{D_{\rm L} t^{\prime \prime }}. \end{equation} \tag{ A5 }$$

On the other hand, dA_S is given by

$$\begin{equation} dA_{\rm S} = \int _{x^{\prime \prime }}^{x^{\prime \prime }+dx^{\prime \prime }} (c_{\rm L0} - c_{\rm S}(x)) \, dx \approx (1-k_0)\, c_{\rm L0}\, dx^{\prime \prime }, \end{equation} \tag{ A6 }$$

where we neglect the spatial variation in c_S and set c_S(x) = k₀c_L0. This assumption is valid in the early phase of the initial transient because the variation in c_S is much smaller than c_L0 − c_S(x). By equating Equations (A5) and (A6), we obtain

$$\begin{equation} B = \frac{(1-k_0)c_{\rm L0}}{\sqrt{D_{\rm L} t^{\prime \prime }}} dx^{\prime \prime }. \end{equation} \tag{ A7 }$$

The position of the interface is given by $x_{\rm i}(t) = V_0 t_0 (e^{t/t_0} - 1)$ since $V(t) = V_0 \, e^{t/t_0}$. Therefore, the period for the interface to move from x = x'' to x = x_i is given by

$$\begin{equation} t^{\prime \prime } = t_0 \ln \frac{ 1 + x_{\rm i} / (V_0 t_0) }{ 1 + x^{\prime \prime }/ (V_0 t_0) }. \end{equation} \tag{ A8 }$$

Therefore, the concentration in the liquid at the interface is calculated to be

$$\begin{eqnarray} c_{\rm Li} (x_{\rm i}) &=& c_{\rm L0} + \int _0^{x_{\rm i}} dc_{\rm L}(0) dx^{\prime \prime } \nonumber \\ &=& c_{\rm L0} \left[ 1 + (1-k_0) \sqrt{\frac{V_0^2 t_0}{D_{\rm L}}} \cdot I \left(\frac{x_{\rm i}}{V_0 t_0} \right) \right] \end{eqnarray} \tag{ A9 }$$

with the use of Equations (A4), (A7), (A8), and

$$\begin{eqnarray} &&I (u) = \int _0^{u} \frac{1}{ \displaystyle { \sqrt{\ln \frac{1+u}{1+u^{\prime }}} } } \ du^{\prime } = \sqrt{\pi } \, (1+u) \, {\rm erf}\, [ \sqrt{\ln (1+u) }]. \nonumber\\ \end{eqnarray} \tag{ A10 }$$

Here, ${ {\rm erf}(x) = ({2}/{\sqrt{\pi }}) \int _0^x e^{-\xi ^2} d\xi }$ is the error function. Finally, we obtain the concentration profile c_S = k₀c_Li in the solid in the early phase of the initial transient by

$$\begin{eqnarray} c_{\rm S}(x_{\rm i}) &=& k_0 c_{\rm L0} \left \lbrace 1 + (1-k_0) \sqrt{\frac{\pi }{D_{\rm L} t_0}} V_0 t_0 (1+x_{\rm i}/V_0 t_0) \right.\nonumber \\ &&\left.\times {\rm erf}\, [ \sqrt{\ln (1+x_{\rm i}/V_0 t_0)}] \right \rbrace. \end{eqnarray} \tag{ A11 }$$

For x_i ≫ V₀t₀, c_S is approximated by

$$\begin{equation} c_{\rm S} (x_{\rm i}) \approx k_0 c_{\rm L0} \left\lbrace 1 + (1-k_0) \sqrt{\frac{\pi }{D_{\rm L} t_0}} \, x_{\rm i} \right\rbrace. \end{equation} \tag{ A12 }$$

Although Equation (A12) gives slightly smaller values of c_S than that given by Equation (A11), the difference is less than ∼25% for any x_i for the range of the parameter values listed in Table 2. Equation (A12) is accurate enough for an order of magnitude estimation.