Role of partial molar enthalpy of oxides on Soret effect in high-temperature CaO–SiO2 melts

The Soret effect or thermodiffusion is the temperature-gradient driven diffusion in a multicomponent system. Two important conclusions have been obtained for the Soret effect in multicomponent silicate melts: first, the SiO2 component concentrates in the hot region; and second, heavier isotopes concentrate in the cold region more than lighter isotopes. For the second point, the isotope fractionation can be explained by the classical mechanical collisions between pairs of particles. However, as for the first point, no physical model has been reported to answer why the SiO2 component concentrates in the hot region. We try to address this issue by simulating the composition dependence of the Soret effect in CaO–SiO2 melts with nonequilibrium molecular dynamics and determining through a comparison of the results with those calculated from the Kempers model that partial molar enthalpy is one of the dominant factors in this phenomenon.

The Soret effect was discovered by C. Ludwig 1 and tested by C. Soret 2 . Although it has been over 150 years, and the Soret effect can be quantitatively explained by the Chapman's theory 3 for the case of molecular gases, the mechanism of the Soret effect in liquid remains controversial [4][5][6] . The Soret coefficient indicates whether the components diffuse toward hot or cold region and provides the separation degree of the components between hot and cold regions. Neglecting convection, the flux under a temperature gradient in a binary system can be written as follows 7 : where x is the position, n 1 is the mole fraction of the species 1, ρ is the mass density, T is the temperature, and D M and D T are the mutual and thermal diffusion coefficients, respectively. In the steady state, the flux J = 0 and the Soret coefficient can be written as follows: should be valid in the binary system. The Soret effect in silicate melts is important in the field of glass and earth sciences, where silicate is a representative component, because it causes the spatial inhomogeneity of the composition in a glass melting container and in the Earth's interior. Many reports on the Soret effect of silicate melts were released [8][9][10][11][12] , but the dominant factor to determine the Soret coefficients remained controversial. In 2010, Huang et al. suggested that the Soret coefficient in silicate melts is expressed by an additive function of the mass-and chemical-effect terms. The mass-effect term is expressed by the variable of mass, charge, and radius of ion. However, the explicit function for the chemical-effect term is unclear 10 . In 2011, Dominiques et al. suggested that the electric energy barrier and the vibrational zero-point energy in the activation process of ion diffusion are important factors in determining the Soret coefficient 11 . When tackling this issue, determining the diffusion species in silicate melts is difficult because silicate melts generally contain many structure types (e.g., bridging, nonbridging, and free oxygen, Q n unit (n is the number of bridging oxygen per SiO 4 unit), and network ring size 13 ). Furthermore, the chemical reaction of the Si-O network 14 and the electric neutrality constraint among ions 15 make the diffusion process complicated. The diffusion and viscous flow units in silicate melts are also unclear. Therefore, it is difficult to employ the kinetic approach, as similar was done for the Artola 4 and Eslamian 5 models, where the Soret coefficient was discussed using the activation energy of the diffusion or viscous flow process of molecular liquid. There has been an approach from a different angle; Kempers model 16,17 is a thermodynamic model mainly used for the molecular gas and liquid system. This model will be effective in dealing with the problem of silicate melt complexity because it does not need to determine the diffusion species. Kempers said that factors contributing the Soret effect are a function of thermodynamic parameters, such as the partial molar enthalpy, partial molar volume, and chemical potential.
The present study discusses the dominant factor contributing to the Soret effect in a CaO-SiO 2 system. We use the nonequilibrium molecular dynamics (NEMD) simulation since we can neglect the convection effect caused by gravity and surface tension. Inspired by Kempers 16,17 , we also employ a thermodynamic approach to the Soret effect in silicate melts. We then compare the simulation result with that of a Kempers model, where we take the segregation limit as the standard state of the thermodynamic parameters to adjust the original Kempers model 16,17 to the silicate melts.

Results
Soret coefficient calculated from the NEMD simulation. We simulated the Soret effect of mCaO-(1 − m)SiO 2 (m = 0.5, 0.6, 0.7, 0.8, and 0.9) melts with the NEMD simulation using approximately 12000 particles at pressure of about 100 MPa(see Method section). The hot region and cold region in the simulation box was kept to be 2200 K and 1800 K, respectively, and linear temperature gradient was obtained. We conducted 22 simulations with different compositions and initial ion positions as shown in Table 1. First, we focus on Run No. 1~21. We waited for the time period of θ under temperature gradient before sampling the mole-fraction distribution, where θ is characteristic time described in Methods section. The sampling period was 2θ. Figure 1 shows mole-fraction distribution and fitted line under the temperature gradient. The variation of the mole-fraction distribution even in the same composition of the system is mainly due to the small number of particles in the system (apporoximately 12000 particles). Every fitted line for 0.9CaO-0.1SiO 2 has a negative gradient, while every fitted line for 0.7CaO-0.3SiO 2 , 0.6CaO-0.4SiO 2 , and 0.7CaO-0.3SiO 2 has a positive gradient. This means that the SiO 2 -concentrated region is changed with SiO 2 content of the system; the SiO 2 concentrates in hot side in a SiO 2 -rich melt, while SiO 2 concentrates in cold side in a SiO 2 -poor melt. As shown in Fig. 1, the deviation of mole fraction from linear relationship against temperature increases with decreasing SiO 2 content, which should be due to the small number of Si ions in the simulation box in the SiO 2 -poor melt. To improve statistics for 0.9CaO-0.1SiO 2 , we conducted a simulation with longer sampling times as shown in Fig. 2. The mole-fraction distribution becomes linear with increasing the sampling time. The Si ions will move around the simulation box during the long-time sampling, which will be the reason why we obtained linear relationship. The small difference of plotted data between the 33.25 ns and 42.75 ns in Fig. 2(b) indicates that the concentration distribution is almost converged. This supports the negative gradient of mole fraction of SiO 2 in the 0.9CaO-0.1SiO 2 system.
In this study, we take the SiO 2 and CaO as component 1 and 2, respectively. Since the gradient of the fitted line corresponds to the ∂n/∂T, by using the Eq. (2), we calculated the Soret coefficients of SiO 2 component and summarized them in Table 1. The positive and negative values of the Soret coefficient mean the SiO 2 concentrates in cold region and hot region, respectively. We have reported the Soret effect in 50CaO-50SiO 2 by a laser irradiation experiment 18 , where the SiO 2 component concentrated in cold side under temperature gradient. This is qualitatively consistent with the result of NEMD simulation shown in Fig. 1 Kempers proposed a model to calculate the thermodynamic effect of the Soret effect 16,17 . Kempers assumed the two-bulb apparatus, which have equal and constant volumes, as shown in Fig. 3. The system has two components and the bulb A and bulb B are kept to be homogeneous temperature of T A = T + ΔT/2 and T B = T − ΔT/2, respectively. T is the average temperature, and ΔT is positive value and indicates the temperature difference between bulb A and bulb B.

Derivation of Kempers model and adjustment to the binary silicate melt.
Kempers assumed that the canonical partition function of the whole system(Z total ) will be maximum at steady state under a temperature gradient: where z A and z B are the canonical partition function of the partial system of bulb A and bulb B, respectively. By using the thermodynamic relationship, where Z is canonical partition function, F is the Helmholtz free energy, R is gas constant, T is temperature, we can modify condition (3) to Constraint conditions are as follows: where N is number of the component, v is the partial molar volume, i is the number of component. The (6) and (7) express law of conservation of mass and equal volume of each bulb, respectively. By using the method of Lagrange multiplier for (5) under the condition of (6) and (7) and approximation v 1 where v 1 and v 2 are the average partial molar volume of bulb A and bulb B, we can obtain the following condition for steady state under a temperature gradient: where μ is chemical potential. By using Taylor expansion, we obtain the following relationship:  The Gibbs-Duhem equation in this system is expressed as: By the Eqs (9-12), we can obtain the relation between Soret coefficient and thermodynamic parameters: In this study, we assume that the mixing thermodynamic parameters of the two components of binary glass melts should be important for the Soret effect. Under this assumption, we use the departure of thermodynamic parameters from the segregation limit of liquid mixture. We modify the Eq. (13) to the following one:    Table 2. ΔH Mix was calculated by H − H SL . ΔS Mix was calculated using a model proposed by P.L. Lin 19 . The H SL can be calculated by simple addition of enthalpy of the pure liquid state: mH CaO + (1 − m)H SiO2 , where m, H CaO , and H SiO2 are mole fraction of CaO, enthalpy of pure CaO liquid and pure SiO 2 liquid, respectively. By using V, ΔH Mix , and ΔG Mix , we calculated v(partial molar volume), h − h°(departure of partial molar enthalpy from pure liquid state), and μ − μ° (departure of chemical potential from pure liquid state) of the mixture, respectively. Finally, we obtained Soret coefficient of Kempers model(σ SiO2 Kempers ) using Eq. (14). Figure 4 shows the mole-fraction dependence of thermodynamic factors which appear in Eq. (14). The h − h° of the SiO 2 increases monotonically with increase of SiO 2 content, whereas that of the CaO decreases. The h − h° has cross over point of CaO and SiO 2 around 0.39SiO 2 . The v and x SiO2 {∂(μ SiO2 − μ SiO2°) /∂x SiO2 } changes gradually with the SiO 2 mole fraction.

Summary of the NEMD simulation and the Kempers model.
We summarize the results of NEMD and Kempers theoretical model in Fig. 5. Both values monotonically decreases with the increasing mole fraction of SiO 2 and obtained the turning point of the sign of the Soret coefficient, with the difference of 0.23 for the NEMD and 0.36 for the Kempers model. Both Soret coefficients seem to change in parallel against the mole fraction of SiO 2 . As shown in Fig. 5, the value of lowest edge of the error bar and the value of long-time sampling in 0.9CaO-0.1SiO 2 are located above zero, which indicates the SiO 2 concentrates in the cold region. The standard deviation increases with decreasing SiO 2 content. The large variation of the Soret coefficient in low-SiO 2 content should be due to the small number of Si ions in the simulation box.

Discussion
Comparison of the NEMD result with the Kempers model. The difference between the NEMD simulation and the theoretical model was almost not dependent on the composition, that is, the values changed in parallel. Kempers said that there are two contributions to the Soret effect: thermodynamic contribution from attraction/repulsion and kinetic contribution from collision interaction between components in Soret effect. Based on this idea, the difference between the NEMD simulation and the theoretical model may come from the kinetic factor because we have already considered the thermodynamic contribution in the Kempers model (Eq. (14)). The kinetic factor comes from the collision behavior, and should be expressed as a function of the factors of diffusion species, such as size, mass, and bond strength. Lacks 12 discussed the contribution of ion mass to the Soret effect caused between isotopes, called isotope fractionation, in silicate melts, and suggested that this contribution to the Soret effect is quantified by scaling relation based on the Chapman-Enskog theory 3    explanation for this is: the penetration depth of the heavier ion from the hot region to the cold region is longer than that of the lighter ion because the former can easily scatter lighter ions. In this study, the diffusion units in the system should be Si−O network, Ca ion, and free O ion. The Si-O network and may act as a heavier diffusion unit than Ca and free O ions, which in turn may contribute to the positive shift to the Soret coefficient of SiO 2 . However, this cannot explain the deviation tendency of the NEMD from the Kempers model. In contrast, analogous to this discussion, the Si-O network behaved as a large-diffusion species compared to Ca and free O ions. The penetration depth of the Si-O network from the hot region to the cold region was shorter than that of the Ca and free O ions because collision frequently occurred in the case of the large-diffusion species. Therefore, the Si-O network will easily concentrate in the hot region through this contribution, and the Soret coefficient of SiO 2 obtained from the NEMD is shifted to a negative direction from that predicted by the Kempers model. We think that the size effect may be larger than the mass effect in the system. Important role of partial molar enthalpy. The sign of the Soret coefficient in Kempers model (Eq. (14)) is determined by the term of (h 2 − h 2°) /v 2 − (h 1 − h 1°) /v 1 . As shown in Fig. 3, the h − h° is drastically changed with composition and cross over point is observed around 0.39SiO 2 14), we obtain the conclusion that the sign change mainly comes from the magnitude relationship of h − h° between the two components. In other words, the main factor to determine the diffusion direction of SiO 2 component is departure of partial molar enthalpy of mixture from pure liquid state. This discussion may be applied to other binary silicate system. In conclusion, both of NEMD and Kempers model showed a monotonic decrease of the Soret coefficient of SiO 2 and a sign change of that at low SiO 2 content. The difference between the two results may be caused by a kinetic factor. According to the Eq. (14) of the Kempers model, mole-fraction dependence of partial molar enthalpy of SiO 2 and CaO should be the cause of the sign change. As a future work, we should confirm the sign change in Soret coefficient at low SiO 2 content binary silicate system, experimentally.

NEMD simulation.
The simulation method was similar to that in our previous report 20 , but the details had some differences. We employed the potential proposed by Seo 21 . The potential was developed for CaO-SiO 2 system and is an extended version of the potential for SiO 2 proposed by Tsuneyuki 22 to CaO-SiO 2 system. The potential function is expressed as: