Adsorption of Cs⁺ Ion into Di- and Tri-Octahedral Vermiculites as Demonstrated by Classical Molecular Dynamics Simulation

Abstract

The intrinsic adsorption of a Cs⁺ ion into di- and tri-octahedral vermiculites without the presence of K⁺ ions was demonstrated by a classical molecular dynamics (MD) simulation. The calculation conditions included Coulomb and Born–Mayer–Huggins potentials, assisted by Lennard–Jones potentials under a constant pressure ensemble and valences from a force field for clays (CLAYFF) mainly as well as conventional valences. A monoclinic di-octahedral vermiculite crystal with a 6 × 3 × 1 supercell was created using crystallographic data from a monoclinic tri-octahedral vermiculite, followed by conversion to a rectangular supercell with periodic boundary conditions along the x-axis. The simulated rectangular supercell of the di- and tri-octahedral vermiculite maintained its crystalline structure for 1 ps at 298 K using a constant step of 0.1 fs. Vacancies with diameters of 0.15 nm, which is nearly equal to the ionic size of Cs⁺, or larger were found at the octahedral (O)-sheet only in the di-octahedral vermiculite simulated with valences from CLAYFF. The further MD simulations were performed by placing a Cs⁺ ion at a vacancy at the O-sheet of the simulated state of the di-octahedral vermiculite, revealing that a vacant site can be a candidate of adsorbing Cs⁺ ion. The low degree of crystallinity of the di-octahedral vermiculite because of the octahedral cationic vacancy and the tilting of hydroxyl (OH) group from perpendicular to (001) provided an additional site for absorbing Cs⁺ ion in the O-sheet. The simulation result of the di-octahedral vermiculite simulated with valences from CLAYFF suggested a novel mechanism for Cs⁺ ions to firmly adsorb into vermiculite without being desorbed again.

Fig. 5 (a), (b) Polyhedral views of the di-octahedral vermiculite simulated with CLAYFF valences without drawing the O ion and containing Cs⁺ ion in the vacancy at its states of t = 0 and 5 ps. (c), (d) The trajectories of Cs⁺, drawn with red curves on the Cs⁺ as a ball view, together with those of Si and Al ions drawn with red and green curves, respectively, of the di-octahedral vermiculite at t = 0 and 5 ps.

1. Introduction

The present paper describes classical molecular dynamics (MD) simulations for vermiculites from a class of soils to investigate their adsorption of radioactive cesium, Cs-137 (${}_{\phantom{0}55}^{137}\text{Cs}$) as well as Cs-134 (${}_{\phantom{0}55}^{134}\text{Cs}$), based on crystallographic and thermodynamic data. Previous theoretical and computational studies regarding the adsorption of Cs⁺ ions in soils include several models that explained the mechanisms of Cs⁺ ion adsorption into soils. For instance, representative models include frayed-edge-sites (FES),^1–7) hydroxyl-intercalated vermiculite (HIV),^6,8–10) and partially vermiculitized biotite, also termed weathered biotite (WB). The WB model recently reported⁴⁾ data of soil contaminated with radioactive cesium collected from Fukushima. Also, a collective adsorption model¹¹⁾ was proposed based on the examination of vermiculites taken from Fukushima. It should be noted that FES is often associated with the presence of K⁺ ion in the inter-layer at the third shrinkage step, which makes it difficult to derive the intrinsic nature of the Cs⁺ ion adsorbed into vermiculite. Furthermore, few studies demonstrate the process of Cs⁺ ions adsorbing into the vermiculite using kinetic models such as MD simulations, using precise crystallographic data of the vermiculites.

Several reports dealt with the adsorption behavior of Cs⁺ ions by MD schemes for clays as represented by first principle (ab initio) MD^7,12,13) and classical MD^6,14,15) as well as experiments.^2–5,16) In a class of classical MD simulations, Lammers et al.¹⁴⁾ performed classical MD simulation for illite nanoparticles and clarified Cs adsorption nature of the basal and edge surface sites. The significance of the results¹⁴⁾ is focusing on the surface of the clay particles, which would include considerable defects due to the surface effect, leading to the particle surface apart from the ideal crystallographic structure. The partial change in the crystallographic structure of clays took place at their interface contacting with water as reported by the author’s work,¹⁵⁾ in which classical MD simulations were performed for vermiculite as well as muscovite dipped into water containing 200 pieces of Cs⁺ ions by referring to crystallographic data of clays. These early studies^14,15) are somewhat common in focusing on the effect of the defects of the clays on the Cs⁺ ion adsorption.

The present study focused on the di-octahedral vermiculites. A reason for this was due to experimental data¹⁷⁾ that the vermiculite in Fukushima is low in Mg²⁺ ions, suggesting a preference of di-octahedral over tri-octahedral vermiculate. Another reason for this was due to a report indicating¹⁸⁾ that vermiculites with a relatively fine grain size contain considerable stacking disorder and are of the di-octahedral type containing Al ions, whereas the tri-octahedral vermiculites with large grain size are rich in Al ions at the Tetrahedral (T)-sheet, including Mg ions that are mainly concentrated at the O-sheet and inter-layers, which vermiculites are called Mg-vermiculite. In addition, it is stated¹⁹⁾ that the vermiculites produced in the Fukushima area are extremely poor in MgO. Furthermore, previous reports indicated^2,20) that the micaceous minerals in the clay of most soils are predominantly di-octahedral, particularly in developed soils. These reports recommended to use the precise crystallographic data of the di-octahedral vermiculite utilizing the existing tri-octahedral Mg-vermiculite^21,22) data for the analysis.

The present study aims to clarify the intrinsic process and mechanism of Cs⁺ ion adsorption into di- and tri-octahedral vermiculites without the presence of K⁺ ions using classical MD simulations. An appropriate and reasonable amount of Cs⁺ ions is considered from experiments and crystallographic and thermodynamic knowledge incorporated from materials science. The purpose of the present study is to perform MD simulations for vermiculite and clarify its absorption behavior of Cs⁺ ion through a simplified MD method as possible.

2. Methods

The present MD method was almost the same as the scheme in the author’s previous study,¹⁵⁾ excepting for (i) the type of vermiculite, (ii) the charge of the ions, (iii) inclusion of environmental water, and the number of Cs⁺ ions in the supercells. The previous study¹⁵⁾ dealt with (i) only the tri-octahedral vermiculite as well as di-octahedral muscovite, (ii) the charge of the ions were acquired from conventional valences, (iii) the supercells dipped into water included 200 species of Cs⁺ ion. In contrast, the present study performed MD simulations for (i) both the di- and tri-octahedral types with (ii) the values of charges by referring to a force field for clays, CLAYFF,²³⁾ and for (iii) the supercells with only one Cs⁺ ion without dipping into environmental water.

Classical MD simulations, which solves Newton’s equations of motion for many-particle systems interacting through pairwise potentials simply and do not deal with the charge transfer of ions as can be done in first principle MD, were performed using the commercial software SCIGRESS ME Ver. 2.3 and SCIGRESS Ver. 2.9 (Fujitsu).²⁴⁾ The Coulomb and Born–Mayer–Huggins (B–M–H) potentials^25–27) from a class of potentials for ionic bonding were both included in the software and primarily adopted for the simulations. The details of the B–M–H potential are described in Appendix A. In addition, the conventional Lennard–Jones (LJ) potential²³⁾ was used as a supplement for the O²⁻–Cs⁺ ion pairs where the O²⁻ ions were constituent ions in vermiculite and H₂O. Moreover, an especial LJ potential was adopted for the O²⁻–Cs⁺ ion pairs of the hydroxyl group included in the T–O–T structure of the vermiculite. The details of the coefficients of the MD potentials are summarized in Table A.1.^24,28) Furthermore, the inter-molecule potentials for the O–H and H₂O (H–O–H) components were not provided in the software, they were therefore treated as rigid bodies in the present study. This treatment (rigid body) does not allow the lengths of O–H and the angle of H–O–H with O at its center to vary but does permit turns and rotation.

Di-octahedral vermiculite was created by referring to a Mg-vermiculite^21,22) from a class of tri-octahedral vermiculite. This Mg-vermiculite was selected to obtain precise crystallographic data to perform MD simulations and to clarify the intrinsic nature of the Cs⁺ ion adsorption into the vermiculite without the presence of K⁺ ions. The Mg-vermiculite crystallographic data were acquired from literature.^21,22) The data summarized in Table 1²⁹⁾ were selected from two sets of crystallographic vermiculite data^29,30) that include H₂O in the inter-layer. Specifically, di-octahedral vermiculite was created from Mg-vermiculite by replacing Mg²⁺ ions of the octahedral sheets with Al³⁺. Two of the three Mg²⁺ ions at the Mg1–3 sites, Mg1 and Mg2 in Table 1, were replaced with Al³⁺, and the Mg3 site was set to be vacant in the di-octahedral vermiculite. This replacement of two Mg²⁺ with three Al³⁺ did not change the total charge within the di-octahedral sheets. It should be noted that a part of Mg²⁺ ions at the inter-layer of the vermiculite occupies the Mg4 site with a site occupancy of 0.41, together with the octahedral-sheet (O-sheet) sites of Mg1–3 with a site occupancy of unity. The constituent ions without occupancy of unity were randomly set in the commercial software. Vermiculite has a monoclinic structure with an axial angle β ≠ 90 degrees. For simulation, the monoclinic structure was converted to a 6 × 3 × 1 rectangular supercell to efficiently adopt the pressure ensemble under periodic boundary conditions. The ionic configurations of the di- and tri-octahedral vermiculite are illustrated in Figs. 1 and 2, respectively. Comparing of Fig. 1(d) left and right shows the lack of octahedral sites in di-octahedral vermiculite. At the initial ionic configuration, the OH group had an orientation perpendicular to the z-axis where the H heads to the direction of the T-sheet as shown in Figs. 1.

Table 1 Crystallographic data of a vermiculite (Al_0.57H_1.4Mg_1.705O_7.86Si_1.43) acquired from literature²⁹⁾ where the space group number is 15 (C12/c1). The lattice constants (a, b, c) and angles (α, β, γ) in the units of nm and degrees are (a, b, c) = (0.533, 0.9255, 2.889) and (α, β, γ) = (90, 97.12, 90), respectively.

Fig. 1

Ionic arrangements of the 6 × 3 × 1 rectangular supercells of (left) the di- and (right) tri-octahedral vermiculites drawn in y = 0–0.49 and y = 0.51–1, respectively. (b), (d), (e), (f), (g) polyhedral views without drawing O²⁻ and (a), (c) ball-and-stick views. Figures 1(a), (e) are exactly equivalent as well as Figs. 1(c), (g). Figure 1(d) depicts the octahedral sheet partially derived to x-direction from Fig. 1(b), indicating a vacancy at the ratio of 1/3 in the di-octahedral vermiculite. The star-shaped symbol indicates the six-membered ring in the tetrahedral sheet, which is composed of tetrahedral SiO₄ heading a base perpendicular to the page.

Fig. 2

Ball-and-stick and polyhedral views of the ionic arrangements of the 6 × 3 × 1 rectangular supercells. (a), (b) di- and (e), (f) tri-octahedral vermiculite simulated with valences from CLAYFF for 1 ps at 298 K, whereas (c), (d) and (g), (h), respectively, depict di- and tri-octahedral vermiculites simulated with conventional valences for comparison.

While selecting the scheme of the classical MD method, attention was also paid to another widely accepted method called LAMMPS^31,32) (large-scale atomic/molecular massively parallel simulator) for its charge of ions in units of coulomb. The following values were acquired from the CLAYFF: Si = +2.1, Al = +1.575, Mg = +1.05, O in H₂O = −0.82, H in H₂O = −0.41, O in vermiculite = −1.1688, O in hydroxyl OH group = −0.95, and H in hydroxyl OH group = +0.425. There is another charge for the Mg ion (+1.36), but this value was not selected because it collapsed the supercells in preliminary simulations for tri-octahedral vermiculite. The ionic radii of the constituent elements were acquired from literature:³³⁾ r_O²⁻ = 0.132 nm, r_Mg²⁺ = 0.078 nm, r_Al³⁺ = 0.057 nm, r_Si⁴⁺ = 0.039 nm, and r_Cs⁺ = 0.165 nm. The ionic radii provided in literature¹⁸⁾ and software differ. These data are summarized in Table 2.

Table 2 Atomic number (N_A), species, and Ionic radius (r_ion) acquired from literature,^18,33–35) together with conventional valences and valences from clay force field (CLAYFF).²³⁾

A characteristic of CLAYFF²³⁾ is the nonbond and bond parameters. In the nonbond parameters, ions are subdivided by combination states, such as hydroxyl oxygen (oh), bridging oxygen (ob), bridging oxygen coordinated to tetrahedral substituted sites (obts), bridging oxygen coordinated to two octahedral substitutions, or tetrahedral-octahedral substitutions (obss). However, the present MD scheme to focus on the simplification of the method did not distinguish the ions in the subdivided states and every O ion in a nonbond state in the T–O–T was dealt with the same O ion. Similarly, octahedral and tetrahedral aluminum with symbols of ao and at, respectively, were the same Al ion as did in the author’s previous study.¹⁵⁾ On the other hand, the author’s previous study¹⁵⁾ distinguished octahedral and hydroxide magnesium of mgo and mgh, respectively, but they were the same Mg ion in the present study. Besides, the present scheme as well as the previous study¹⁵⁾ did not deal with bond parameters due to the treatment of rigid body for O–H and H₂O molecule because of the simplification.

The present MD simulations with constant pressure (p) ensemble, including a constant number of atoms/ions (N) and absolute temperature (T), were performed to examine the validity of applying mainly the B–M–H potential in this study. The simulation conditions included a simulation time up to 1 (or 5) ps with a constant step (Δt) of 0.1 fs. Before performing the MD simulation, the rectangular supercell was subjected to structural relaxation at a constant volume (V) at T = 100 K. The simulation results were analyzed for their ionic arrangements with commercial software CrystalMaker Ver. 10.3.2,³⁴⁾ and freeware VESTA Ver. 3.4.0 (64-bit Edition).³⁵⁾ The analysis of the vacancies for the simulation results was performed with CrystalMaker by using the conventional ionic radii for the constituent ions of the vermiculite in the software. Besides, the number of Cs⁺ ions included in the supercell was determined to be only one ion per supercell by considering experimental data in previous studies^36–39) as described in Appendix B.

3. Results and Discussion

3.1 Di- and tri-octahedral vermiculites simulated for 1 ps at 298 K

The ion configurations of the di- and tri-octahedral vermiculites simulated for t = 1 ps at T = 298 K using the NTp ensemble are presented in Fig. 2 where Figs. 2(a), (b), (e), (f) and Figs. 2(c), (d), (g), (i) exhibit the results simulated with valences from CLAYFF and conventional valences, respectively. Figure 2 exhibits that MD simulations were successfully performed without causing overshoot under the given MD conditions. The reproducibility of the simulation results for the tri-octahedral vermiculite shown in Figs. 2(a), (e) was analyzed in terms of the lattice constants, density, and basal d spacing of the tri-octahedral vermiculite, as examined the validity of the CLAYFF²³⁾ for boehmite γ-AlO(OH), portlandite Ca(OH)₂, kaolinite Al₂Si₂O₅(OH)₄, and pyrophyllite AlSi₂O₅(OH) in terms of lattice constants (a, b, c), lattice angles (α, β, γ), density, basal d spacing, and specific distances such as O–H, Al–OH, Si–O_basal, Si–O_apical. The lattice constants of the supercells (a, b, c) of the tri-octahedral vermiculite in units of nm varied from (3.24, 2.81, 2.90) to (3.60, 3.14, 2.70) during the simulation. Thus, the lattice constants expanded approximately 11.1% and 12.0% in x- and y-directions, whereas shrank 6.6% in z-direction. The resultant density (ρ) varied from 2.140 to 1.782 Mg/m³, exhibiting a 17% decrease during the MD simulation. A referential simulation results shown in Figs. 2(c), (g) performed with conventional valences revealed that (a, b, c) of the tri-octahedral vermiculite varied from (3.24, 2.81, 2.90) to (3.18, 2.75, 2.72) and the resultant lattice constants shrank approximately 1.9%, 2.1% and 6.2% in x-, y- and z-directions, respectively, accompanied by the variation of ρ from 2.140 to 2.302 Mg/m³ corresponding to 7.6% increase during the MD simulation. These different tendencies between the simulation results with valences from CLAFYY and conventional ones can be confirmed by focusing on the aspect ratio of the edges of the supercells in Fig. 2. The above analysis revealed that the referential simulation with conventional valences reproduced the experimental data²⁹⁾ better than that with CLAFF due to larger Coulomb interactions for the former than the latter cases. The vermiculite simulated with conventional valences was more rigid (higher ρ) than that with CLAYFF, suggesting that the former can be regarded as a bulk body whereas the latter as a surface part. Below, the present study focuses on the latter, in particular, to the peculiar possibility of the O-sheet that can accommodate Cs⁺ ion in the di-octahedral vermiculite, although the former exhibited the better reproducibility to the experimental data.

Thus, the present MD simulation for the tri-octahedral vermiculite did not reproduce the crystallographic data^21,22) from experiments accurately within errors less than a couple of percentage. However, the present results barely kept the crystalline structure of the tri-octahedral vermiculites, which are sometimes classified into 14-angstrom interlayered clays. By considering the simplest MD simulation scheme of the present study, the present MD simulation for the tri-octahedral vermiculite was barely valid in terms of satisfying the minimum necessary to exhibit the crystalline structure of the tri-octahedral vermiculite. The features of the simulation results are presented in Fig. 2 and characterized in terms of crystallinity, orientation of the OH group, and the presence of vacancy.

It appears from Fig. 2 that the supercells remained crystalline structures with slight random orientations of the polyhedra retaining the crystalline structures. In more detail, the polyhedra in the simulated vermiculites displayed in Fig. 2 were compared to the ideal crystalline state irregularly connected due to distortion, sharing vertices and edges. Comparing the irregularity between the polyhedra with the ideal state, the degree of crystallinity of the di-octahedral vermiculite seems to be lower than that of the tri-octahedral vermiculite. With respect to the OH group orientation, the angle δ between the hydroxyl groups and the (001) of the vermiculite was initially 90 degrees as shown in Figs. 1. In the di-octahedral vermiculite δ decreased while it remained almost unchanged in the tri-octahedral vermiculite as shown in Fig. 2. The different behavior of δ agreed with earlier reports,^40,41) where it was explained¹⁸⁾ that the vacancy of an octahedral site affects the orientation of the OH groups in the di-octahedral vermiculite due to the inhomogeneity of the electrical charge. It was also reported in previous studies^42,43) for talc and pyrophyllite from a class of soil that the tri-octahedral talc exhibits hydroxyl groups configured perpendicular to the clay layer (001) while the dioctahedral pyrophyllite shows the hydroxyls situated sub-parallel to the layer in response to the presence of octahedral vacancies. In addition, it was observed that the di-octahedral vermiculite simulated with CLAYFF valences for 1 ps at 298 K gained additional space inside the T–O–T layer as shown in Fig. 3(a). Specifically, the analysis of the excess space of the vermiculites simulated at a given temperature revealed that di- and tri-octahedral vermiculites contain in the simulated states 59 and 32 vacant sites, respectively, with diameters larger than 0.15 nm. These vacant sites were distributed in the di-octahedral vermiculites at both T–O–T and inter-layers, whereas those vacancies were only observed at inter-layers in tri-octahedral vermiculites simulated with CLAYFF valences. On the other hand, both di- and tri-octahedral vermiculites simulated with conventional valences contained 33 vacant sites, in which no vacant sites were observed at T–O–T layers even in the di-octahedral sites. Thus, it was anticipated that the vacant sites observed at T–O–T layer in the di-octahedral vermiculite were affected by an intrinsic characteristic of the di-octahedral vermiculite caused by its low crystallinity including large amounts of defects that could take place at a surface area.

Fig. 3

Ball-and-stick views of the ionic arrangements of the 6 × 3 × 1 supercells of (a) di- and (b) tri-octahedral vermiculites simulated with CLAYFF valences for 1 ps at 298 K, together with vacancies with radii >0.15 nm drawn in transparent black spheres denoted as Zz. The analysis for (c) di- and (d) tri-octahedral vermiculites simulated with conventional valences are also shown for comparison. Every vacancy was drawn as a transparent circle because of the projection to (001), where some of the overlapped area of the vacancies are drawn with increased darkness due to the projection. The ionic radii of the constituent ions of the vermiculite and those of the vacancy drawn in ball view are relatively proportional to each other. Most of the vacancies were observed in the inter-layer and the place near the T-sheet. The vacancies in the O-sheet were observed only in the di-octahedral vermiculite in Fig. 3(a).

The changes of Hamiltonian ($\mathcal{H}$) of the di- and tri-octahedral vermiculite supercells simulated with valences from CLAYFF are depicted in Fig. 4, indicating that the simulated di- and tri-octahedral vermiculites were in a state that would be close to the equilibrated state as a result of $\mathcal{H}$ steadily decreasing with t at 298 K from each initial state.

Fig. 4

Changes in Hamiltonian ($\mathcal{H}$) of (a), (c) tri- and (b), (d) di-octahedral vermiculite supercells at 298 K. Figures 4(a) and (b) were the results simulated with valences from CLAYFF, whereas Figs. 4(c) and (d) were the ones with conventional valences.

3.2 Behavior of Cs⁺ ion at O-sheet simulated for di-octahedral vermiculite for 1 ps at 298 K

As the thermally stabilized di-octahedral vermiculite contained a vacancy site in the T–O–T layers, hypothetical MD simulations were performed in advance to analyze the behavior of the Cs⁺ ion at the O-sheet. For this simulation, the MD potential of UFF (Universal Force Field) in a class of LJ potential was adopted for O–Cs ion pairs, as shown in Table A.2,^24,28) to confirm if the Cs⁺ ion as large as approximately 0.180 nm can truly be accommodated at the vacancy. Specifically, of the 59 vacant sites with diameters larger than 0.15 nm, four sites had a radius larger than r_Cs⁺ ∼ 0.180 nm. Thus, a Cs⁺ ion was hypothetically placed as the initial state at a vacancy in the O-sheet as shown in Figs. 5(a) and (c), where the size of the vacancy radius was taken as 0.19 nm. The MD simulation results at t = 5 ps, as shown in Fig. 5(b) and (d), indicate that the Cs⁺ ion remained in the O-sheet. The trajectory of the Cs⁺ ion during the MD simulation is presented as red curves in Figs. 5(c) and (d) and is drawn on those of the Cs⁺ ion as a ball view, including the other trajectories of Si and Al ions drawn as curves. The trajectory of the Cs⁺ ion shows that its position remained unchanged and just thermal vibration around a site took place, which was similar to the other Si and Al ions. Figure 5 suggests that a vacancy in the O-sheet in the di-octahedral vermiculite can be a candidate, although the diffusion kinetics of the Cs⁺ ion from the inter-layer were not validated in the hypothetical MD simulation shown in Fig. 5. The average distance between O–Si, O–Al, and O–Cs at t = 0 ps was 0.182 (O–Si: tetrahedra), 0.2061 (O–Al: octahedra, and 0.3084 (O–Cs) nm, which varied to be 0.182 (O–Si: tetrahedra), 0.2068 (O–Al: octahedra, and 0.2885 (O–Cs) nm at t = 5 ps. The decrease in the distance of O–Cs suggested that ionic bonding took place between O–Cs pairs according to Fig. A.1, mainly because of the Coulomb potential.

Fig. 5

(a), (b) Polyhedral views of the di-octahedral vermiculite simulated with CLAYFF valences without drawing the O ion and containing Cs⁺ ion in the vacancy at its states of t = 0 and 5 ps. (c), (d) The trajectories of Cs⁺, drawn with red curves on the Cs⁺ as a ball view, together with those of Si and Al ions drawn with red and green curves, respectively, of the di-octahedral vermiculite at t = 0 and 5 ps.

3.3 Consideration about kinetics of a Cs⁺ ion into O-sheet

In the previous Sub-section, it was demonstrated that a Cs⁺ ion placed at an O-sheet in the di-octahedral vermiculite simulated with valences from CLAYFF can be stable under given MD potentials. Now, there remains a question regarding the actual pathway of Cs⁺ ions to be absorbed into the O-sheet kinetically. To clarify this actual pass, the author tried to performing MD simulations by placing Cs⁺ at an interlayer, failing to demonstrate the Cs⁺ to intrude into inside T–O–T layers. The following can be candidates to explain the reasons for this failure: shortage of computation time and shortcoming of Cs–O potentials. As far as the former candidate is concerned, the author thinks that adding only one Cs⁺ ion into the di-octahedral vermiculite did not cause a significant decrease in $\mathcal{H}$, which can be anticipated from a saturated state of $\mathcal{H}$ in Fig. 4. This suggested that the driving force of the Cs⁺ ion to intrude into T–O–T layers was not large enough, and thus, an extremely long time may be required for a Cs⁺ to move from interlayer to an O-sheet. As for the latter candidate, further tentative simulations were performed by differencing O(OH)–Cs potentials to the O–Cs by applying the following additional condition: a larger E₀ value of the LJ potential for O(OH)–Cs by a couple of hundred times than that for O–Cs. This additional condition originated from an idea that O ions in the T-sheet at the surface of the T–O–T sheet directly contact to H₂O molecules, whereas O ions of OH group of the vermiculite were fully surrounded by constituent ions of the vermiculite due to ionic bonding. Thus, O–Cs interactions of the O ion at the T-sheet would be weakened than O(OH)–Cs interactions inside the T–O–T layers. The tentative simulation results, which are not shown, tended Cs⁺ ion to move from interlayer to an O-sheet in the di-octahedral vermiculite. However, this tendency was not validated because of the author’s intentional variation of the O(OH)–Cs potential from O–Cs potential. In general, such a change in MD potential should be accompanied by accurate energy analysis, such as ab initio calculations. Thus, it was tentatively concluded that the present study was not able to clarify the actual pathway of the Cs⁺ to intrude into O-sheet. Further researches about the kinetical pathways of the Cs⁺ ion to intrude into O-sheet in the di-octahedral vermiculite simulated with valences from CLAYFF will be presented based on the above tentative simulations in near future.

4. Conclusions

Classical MD simulations with Coulomb and Born–Mayer–Huggins potentials primarily combined by Lennard–Jones potentials under a constant pressure ensemble were performed to analyze the intrinsic adsorption behavior of a Cs⁺ ion in vermiculites without K⁺ ions to shrink the inter-layer. Comparisons of the results performed without including Cs⁺ ion and with valences from CLAYFF and conventional valences to the experimental data revealed that conventional valences led to better reproducibility than that using CLAFF. It was found that the di-octahedral vermiculite simulated with valences from CLAYFF only exhibited vacancies with diameters of 0.15 nm, which is nearly equal to the ionic size of Cs⁺, or larger at O-sheet. Further MD simulations after placing a Cs⁺ ion at O-sheet demonstrated that the Cs⁺ ion could be accommodated at an O-sheet of the di-octahedral vermiculite in the case of the simulations with valences from CLAYFF. The actual pathways of Cs⁺ ion to intrude into the O-sheet was not clarified in the present MD simulations. The characteristics of the di-octahedral vermiculite in terms of the lower crystallinity compared to the tri-octahedral vermiculite, the presence of vacancies in the octahedral sheet, and the tilting of the hydroxyl group toward the vacancy direction, provided novel adsorption sites for Cs⁺ ions. The present MD simulations suggested that adsorption of the Cs⁺ ion takes place at an O-sheet of the di-octahedral vermiculite, because of crystallographic defects promoted by the surface effect. The simulations suggest that defects inherent to di-octahedral vermiculite are key factors when considering the adsorption of Cs⁺ ions, where Cs⁺ ions are strongly adsorbed but hardly desorbed from the soil.

Acknowledgments

The author thanks Dr. Hiromu Arai and Prof. Akira Hasegawa from REER (Research Center for Remediation Engineering of Living Environments Contaminated with Radioisotopes), Graduate School of Engineering, Tohoku University, and Dr. Keizo Ishii, professor emeritus of Tohoku University, for their fruitful discussion on the experimentally observed features of vermiculite. I specifically thank Dr. Keizo Ishii, professor emeritus of Tohoku University, who advised the author to include the contents of the number of Cs⁺ ions from an experimentally observed specific activity. This work was partially supported by JSPS KAKENHI Grant Number JP17H03375.

REFERENCES

• 1)   B.L.  Sawhney: Clays Clay Miner. 20 (1972) 93–100.
• 2)   A.  Nakao,  Y.  Thiry,  S.  Funakawa and  T.  Kosaki: Soil Sci. Plant Nutr. 54 (2008) 479–489.
• 3)  “Environmental Behavior of a Minute Amount of Radioactive Cesium –Cs-Sorption/Desorption Behavior of Clay Minerals Considering the Actual Contamination Condition in Fukushima–”, https://rdreview.jaea.go.jp/review_en/2016/e2016_1_13.html, accessed July 6, 2020.
• 4)   H.  Mukai,  T.  Hatta,  H.  Kitazawa,  H.  Yamada,  T.  Yaita and  T.  Kogure: Environ. Sci. Technol. 48 (2014) 13053–13059.
• 5)   A.J.  Fuller,  S.  Shaw,  M.B.  Ward,  S.J.  Haigh,  J.F.W.  Mosselmans,  C.L.  Peacock,  S.  Stackhouse,  A.J.  Dent,  D.  Trivedi and  I.T.  Burke: Appl. Clay Sci. 108 (2015) 128–134.
• 6)   L.K.  Zaunbrecher,  R.T.  Cygan and  W.C.  Elliott: J. Phys. Chem. A 119 (2015) 5691–5700.
• 7)  “Supercomputer Solves Mystery of Strong Affinity of Soil to Radioactive Cesium –First-Principles Calculation Reveals Mechanism of Cesium Adsorption to Soil–”, https://rdreview.jaea.go.jp/review_en/2013/e2013_1_5.html, accessed July 6, 2020.
• 8)   M.  Goto,  R.  Rosson,  W.C.  Elliott,  J.M.  Wampler,  S.  Serkiz and  B.  Kahn: Clays Clay Miner. 62 (2014) 161–173.
• 9)   A.  Nakao,  S.  Funakawa and  T.  Kosaki: Eur. J. Soil Sci. 60 (2009) 127–138.
• 10)   A.M.  Elprince,  C.I.  Rich and  D.C.  Martens: Water Resour. Res. 13 (1977) 375–380.
• 11)   R.  Motokawa,  H.  Endo,  S.  Yokoyama,  S.  Nishitsuji,  T.  Kobayashi,  S.  Suzuki and  T.  Yaita: Sci. Rep. 4 (2014) 6585.
• 12)   S.V.  Churakov: J. Phys. Chem. B 110 (2006) 4135–4146.
• 13)   M.  Okumura,  H.  Nakamura and  M.  Machida: J. Phys. Soc. Jpn. 82 (2013) 033802.
• 14)   L.N.  Lammers,  I.C.  Bourg,  M.  Okumura,  K.  Kolluri,  G.  Sposito and  M.  Machida: J. Colloid Interface Sci. 490 (2017) 608–620.
• 15)  A. Takeuchi: Adsorption behavior of Cs⁺ ions to vermiculite demonstrated by molecular dynamics simulations, Proc. The 9th edition of the BioPIXE symposium (BioPIXE 2018), (World Scientic Publishing Company, 2018) pp. 1–5.
• 16)   L.  Dzene,  E.  Tertre,  F.  Hubert and  E.  Ferrage: J. Colloid Interface Sci. 455 (2015) 254–260.
• 17)   T.  Ohnuma and  K.  Ishii: Cogent Eng. 4 (2017) 1326200.
• 18)   M.  Sambonsugi: J. Jpn. Assoc. Min. Petr. Econ. Geol. 39 (1955) 76–87 (in Japanese).
• 19)  H. Shirozu: Introduction to Clay Mineralogy — Fundamentals for Clay Science —, (Asakura Publishing Co., Ltd., Tokyo, 1988), pp. 185 (in Japanese).
• 20)  D.L. Sparks and P.M. Huang: Physical Chemistry of Soil Potassium, ed. by D.L. Sparks and P.M. Huang, (ASA–CSSA–SSSA Books, Hoboken, 1985).
• 21)  “Crystallography Open Database”, http://www.crystallography.net/cod/cod/index.php, accessed July 6, 2020.
• 22)  “American Mineralogist Crystal Structure Database”, http://rruff.geo.arizona.edu/AMS/amcsd.php, accessed July 6, 2020.
• 23)   R.T.  Cygan,  J.J.  Liang and  A.G.  Kalinichev: J. Phys. Chem. B 108 (2004) 1255–1266.
• 24)  “FUJITSU Technical Computing Solution SCIGRESS”, http://www.fujitsu.com/global/solutions/business-technology/tc/sol/scigress/, accessed March 4, 2016.
• 25)  A.M. Stoneham and E. Atomic Energy Research: Handbook of Interatomic Potentials: 1: Ionic Crystals, (UKAEA Atomic Energy Research Establishment, London, 1981).
• 26)   F.G.  Fumi and  M.P.  Tosi: J. Phys. Chem. Solids 25 (1964) 31–43.
• 27)   M.P.  Tosi and  F.G.  Fumi: J. Phys. Chem. Solids 25 (1964) 45–52.
• 28)  K. Hirao and K. Kawamura: Material Design Using Personal Computer, (Shokabo, Tokyo, 1994), pp. 217 (in Japanese).
• 29)   H.  Shirozu and  S.W.  Bailey: Am. Mineral. 51 (1966) 1124–1143.
• 30)   A.M.  Mathieson and  G.F.  Walker: Am. Mineral. 39 (1954) 231–255.
• 31)   S.  Plimpton: J. Comput. Phys. 117 (1995) 1–19.
• 32)  “LAMMPS Molecular Dynamics Simulator”, http://lammps.sandia.gov/, accessed July 6, 2020.
• 33)  Metal Databook, (Maruzen Co. Ltd., Tokyo, 2004) p. 8 (in Japanese).
• 34)  “CrystalMaker SOFTWARE”, http://www.crystalmaker.com/index.html, accessed July 6, 2020.
• 35)  “VEATA”, https://jp-minerals.org/vesta/en/, accessed July 6, 2020.
• 36)   T.  Ohnuma and  K.  Ishii: Cogent Eng. 5 (2018) 1429985.
• 37)  K. Ishii, H. Arai, T. Ohnuma and A. Takeuchi: Analysis of strontium in soil and plant samples contaminated during the 2011 Fukushima Daiichi Nuclear Power Plant accident, Proc. The 9th Edition of the BioPIXE Symposium (BioPIXE 2018), (World Scientic Publishing Company, 2018) pp. 21–27.
• 38)   K.  Tagami,  S.  Uchida,  Y.  Uchihori,  N.  Ishii,  H.  Kitamura and  Y.  Shirakawa: Sci. Total Environ. 409 (2011) 4885–4888.
• 39)   M.  Komori,  A.  Shozugawa,  N.  Nogawa and  M.  Matsuo: Bunseki Kagaku 62 (2013) 475–483 (in Japanese).
• 40)   R.F.  Giese: Clays Clay Miner. 27 (1979) 213–223.
• 41)   P.A.  Schroeder: Am. Mineral. 75 (1990) 983–991.
• 42)   R.T.  Cygan,  J.A.  Greathouse,  H.  Heinz and  A.G.  Kalinichev: J. Mater. Chem. 19 (2009) 2470–2481.
• 43)   S.  Tesson,  M.  Salanne,  B.  Rotenberg,  S.  Tazi and  V.  Marry: J. Phys. Chem. C 120 (2016) 3749–3758.
• 44)   A.K.  Rappe,  C.J.  Casewit,  K.S.  Colwell,  W.A.  Goddard and  W.M.  Skiff: J. Am. Chem. Soc. 114 (1992) 10024–10035.
• 45)   Z.W.  Peng,  C.S.  Ewig,  M.J.  Hwang,  M.  Waldman and  A.T.  Hagler: J. Phys. Chem. A 101 (1997) 7243–7252.

Appendices

Appendix A

The total energy in the simulations is determined by considering Coulombic (electrostatic) interactions and short-range interactions for oxides containing ions from the Born–Mayer–Huggins (B–M–H) type.   

\begin{equation} E_{\text{total}} = E_{\text{B-M-H}} + E_{\text{Coul.}}. \end{equation} (A.1)

The following explains the details of the B–M–H potential,^25–27) which is a type of two-body molecular dynamics (MD) potential. The target materials of the potentials are oxides containing H, Li, Be, B, C, N, O, F, Na, Mg, Al, Si, P, S, Cl, K, Ca, Sc, Ti, Ga, Ge, Br, Rb, Sr, Y, Zr, Cs, or Ba, and the electrical charges should be assigned according to the valence of the periodic table. The formula of the potential is given by eq. (A.2), where E is the potential energy, r the distance between species i and j, and A_ij, B_ij, C_ij, and D_ij are the constants for the i–j interactions:   

\begin{equation} E = A_{\text{ij}}\exp (- B_{\text{ij}}r) - \frac{C_{\text{ij}}}{r^{6}} - \frac{D_{\text{ij}}}{r^{8}}, \end{equation} (A.2)

where A_ij represents the Pauling factor determining the magnitude of the repulsive force and B_ij is the parameter determined by the ionic radius and softness, and C_ij and D_ij are dipole interaction parameters. The values summarized in Table A.1 were acquired from literature.²⁴⁾

In contrast, the Coulomb energy was evaluated based on eq. (A.3), where e is the charge of the electron, ε₀ the dielectric permittivity of vacuum (8.85 × 10⁻¹² F/m), q_i and q_j the partial charges, and r_ij the distance of separation of i–j ions.   

\begin{equation} E_{\text{Coul.}} = \frac{e^{2}}{4\pi \varepsilon_{0}}\sum_{i \neq j}\frac{q_{\text{i}}q_{\text{j}}}{r_{\text{ij}}}. \end{equation} (A.3)

Figure A.1 depicts the representative energies of O–Cs, O–Mg, O–Al, and O–Si ion pairs calculated based on eqs. (A.1)–(A.3). Figure A.1 shows that E_total for O–Cs has a minimum value of approximately r = 0.08 nm, which is considerably shorter than that of O–Mg, O–Al, and O–Si. Figure A.1 depicts the B–M–H, Coulomb, and total energies calculated as a function of distance for O–Mg, O–Al, O–Si, and O–Cs ion pairs and LJ-12, universal force field (UFF), and constant force field (CFF) for O–Cs atomic pairs. It should be noted that the ionic radius of Cs⁺ in the Cs⁺–O²⁻ ion pair from the B–M–H potential is artificially shorter than that of the Mg²⁺, Al³⁺ and Si⁴⁺–O²⁻ ion pairs. A similar tendency is seen in the CFF potential, but LJ-12 and UFF potentials provide a rational size of Cs⁺, which is larger than that of Mg²⁺, Al³⁺, and Si⁴⁺. The parameters of the UFF potentials for the O–Cs ion pair are listed in Table A.2, where the total energy (E) is formulated in eq. (A.4) as a function of distance (r) from conventional Lennard–Jones potential. The value of E₀ for the O(OH)–Cs ion pair was estimated through its optimization in the present study.   

\begin{equation} E = E_{0}\left\{\left(\frac{R_{0}}{r} \right)^{12}{} - 2\left(\frac{R_{0}}{r} \right)^{6} \right\}, \end{equation} (A.4)

The actual values of E₀ and R₀ used in the present study are summarized in Table A.2. As shown in Fig. A.1, the actual UFF potential as a function of r for the O–Cs provides almost the same magnitude of E as that of LJ potential, and this difference is rather smaller than E_Coulomb.

Appendix B

The number of Cs⁺ ions involved in the present MD simulations was selected from experimental data of highly contaminated soils. Specifically, contaminated soil acquired 30 km away from the Fukushima Daiichi Nuclear Power Plant (FDNPP) had a specific activity of 25000 Bq/kg, with soil particle diameter <0.5 mm, as indicated in an earlier study.³⁶⁾ In addition, other experimental data indicate that the specific activity of the contaminated soils was approximately 817 kBq/kg (= 817 Bq/g),³⁷⁾ collected at a location approximately 3 km from the FDNPP. Furthermore, the concentration ratios of the concentrations of radioactive cesium (Cs-134 and Cs-137) involved in the disaster were approximately 1:1.^38,39) Thus, it appeared that a specific activity in the order of approximately 400 Bq/g was an appropriate value for performing the present MD simulations. These experimental data directly denied the grouping of Cs⁺ ions in the vermiculite. Specifically, the specific activity of 400 Bq/g resulted in approximately 1 Cs⁺/µm³ according to the following simple calculations using the decaying quantity (λ), which is related to the mean lifetime of the decaying quantity (τ). More precisely, two species of radioactive cesium (Cs-137 and Cs-134) exist and their τ values are 2.0652 and 30.07 years, respectively. The corresponding values of λ for Cs-137 and Cs-134 are 7.309 × 10⁻¹⁰ and 1.064 × 10⁻⁸/sec, respectively. This fundamental knowledge together with an approximate specific vermiculite gravity of ρ ∼ 2 g/cm³ helped to convert the value of B = 400 Bq/g into the number of Cs⁺ ions (Cs-137) per cubic micrometers (µm³) as approximately unity, using the formula Bρλ⁻¹·10⁻¹². Here, it should be noted that a highly contaminated soil (400 Bq/g), led to one ion per µm³ from a microscopic point of view. This result does not agree with previous models such as the FES model and others that simulated higher amounts of Cs⁺ ions in the vermiculite. Consequently, the present study considered a small Cs⁺ ion amount as reasonable in the MD simulation. In more detail, the initial supercells had an approximate volume of 31.3 nm³ (from a = 3.53 nm, b = 3.18 nm, and c = 2.79 nm) in which the di- and tri-octahedral vermiculites possess 1600 and 1620 constituent ion species, respectively, including the OH groups and H₂O molecule at the inter-layer. The application of the periodic boundary conditions compensated for the small volume (31.3 nm³ compared to approximately 1 µm³) with respect to the experimental data mentioned above. Hoverer, the resultant number of Cs⁺ ions applied for the periodic boundary was approximately 10⁸ times higher than the approximate value of 1 Cs⁺ ion/µm³. This disagreement in terms of the number of Cs⁺ ion per µm³ cannot be solved in conventional MD simulations. Thus, as an alternative solution to solve this disagreement in this study, the cutoff distance of the MD potentials was set to be a 1 nm ion, which prevented the Cs⁺ ions from interacting with the other Cs⁺ ions through the periodic boundary conditions in the supercell. At least, the interference of the Cs⁺ ions through the periodic boundary conditions was avoided, therefore, the Cs⁺ ion was regarded as being isolated.

Table A.1 Parameters for the Born–Mayer–Huggins potential used in the present study.^24,28)

Table A.2 The parameters of the UFF (Universal Force Field) in a class of LJ potentials for O–Cs²⁴⁾ and those of a LJ potential.²³⁾

Fig. A.1

Born–Mayer–Huggins (B–M–H) potential,^25–27) Coulomb potential, and total potential calculated using eqs. (A.1)–(A.3) and with (a) valences from CLAYFF as well as (b) conventional valences. The Lennard–Jones (LJ-12),²³⁾ universal force field (UFF),⁴⁴⁾ and constant force field (CFF)⁴⁵⁾ potentials are also shown for comparison. The asterisk indicates the total potential used in the present MD simulation. The UFF potential is a type of LJ potential, whereas the CFF potential simulates the diffusion of hydrocarbon in zeolite. The charge of ions in the units of coulomb were acquired from the CLAYFF.