# Molecular modeling of the effects of 40Ar recoil in illite particles on their K-Ar isotope dating

M. Szczerba, A. Derkowski, Andrey G. Kalinichev, J. Środoń

► To cite this version:

M. Szczerba, A. Derkowski, Andrey G. Kalinichev, J. Środoń. Molecular modeling of the effects of 40Ar recoil in illite particles on their K-Ar isotope dating. Geochimica et Cosmochimica Acta, 2015, 159, pp.162-176. 10.1016/j.gca.2015.03.005. in2p3-01166255

HAL Id: in2p3-01166255

https://hal.in2p3.fr/in2p3-01166255

Submitted on 28 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Molecular modeling of the effects of <sup>40</sup>Ar recoil in illite particles on their K–Ar isotope dating

Marek Szczerba<sup>a,*</sup>, Arkadiusz Derkowski<sup>a</sup>, Andrey G. Kalinichev<sup>b</sup>, Jan Środoń<sup>a</sup>

The radioactive decay of <sup>40</sup>K to <sup>40</sup>Ar is the basis of isotope age determination of micaceous clay minerals formed during diagenesis. The difference in K–Ar ages between fine and coarse grained illite particles has been interpreted using detrital-au- thigenic components system, its crystallization history or post-crystallization diffusion. Yet another mechanism should also be considered: natural <sup>40</sup>Ar recoil. Whether this recoil mechanism can result in a significant enough loss of <sup>40</sup>Ar to provide observable decrease of K–Ar age of the finest illite crystallites at diagenetic temperatures – is the primary objective of this study which is based on molecular dynamics (MD) computer simulations.

All the simulations were performed for the same kinetic energy (initial velocity) of the <sup>40</sup>Ar atom, but for varying recoil angles that cover the entire range of their possible values. The results show that <sup>40</sup>Ar recoil can lead to various deformations of the illite structure, often accompanied by the displacement of OH groups or breaking of the Si–O bonds. Depending on the recoil angle, there are four possible final positions of the <sup>40</sup>Ar atom with respect to the 2:1 layer at the end of the simulation: it can remain in the interlayer space or end up in the closest tetrahedral, octahedral or the opposite tetrahedral sheet. No sim- ulation angles were found for which the <sup>40</sup>Ar atom after recoil passes completely through the 2:1 layer. The energy barrier for <sup>40</sup>Ar passing through the hexagonal cavity from the tetrahedral sheet into the interlayer was calculated to be 17 kcal/mol. This reaction is strongly exothermic, therefore there is almost no possibility for <sup>40</sup>Ar to remain in the tetrahedral sheet of the 2:1 layer over geological time periods. It will either leave the crystal, if close enough to the edge, or return to the interlayer space. On the other hand, if <sup>40</sup>Ar ends up in the octahedral sheet after recoil, a substantially higher energy barrier of 55 kcal/mol prevents it from leaving the TOT layer over geological time.

Based on the results of MD simulations, the estimates of the potential effect of <sup>40</sup>Ar recoil on the K–Ar dating of illite show that some of <sup>40</sup>Ar is lost and the loss is substantially dependent on the crystallite dimensions. The <sup>40</sup>Ar loss can vary from 10% for the finest crystallites (two 2:1 layers thickness and <0.02 µm in diameter) to close to zero for the thickest and largest (in the ab plane) ones. Because the decrease of the K–Ar estimated age is approximately proportional to the <sup>40</sup>Ar loss, the finer crystallites show lower apparent age than the coarser ones, although the age of crystallization is assumed equal for all the crystallites. From the model it is also clear that the lack of K removal from illite fringes (potentially Ar-free) strongly increases the apparent age differences among crystallites of different size.

## 1. INTRODUCTION

The radioactive decay of <sup>40</sup>K to <sup>40</sup>Ar is the basis of the isotope age determination of micaceous clay minerals: illite, glauconite, aluminoceladonite, sericite (fine variety of

<footnote><sup>a</sup> Institute of Geological Sciences, Polish Academy of Sciences, Kraków, Poland</footnote>

<footnote><sup>b</sup> Laboratoire SUBATECH, Ecole des Mines de Nantes, Nantes, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: ndszczer@cyf-kr.edu.pl (M. Szczerba).</footnote>
muscovite), and their mixed-layered interstratifications with
smectite. These dating methods are based on the assump-
tion that the mineral of interest remained a closed system
since crystallization (e.g., Clauer et al., 1997). Illite is the
most common K-bearing clay mineral. Its structure is typi-
cal for a 2:1 layer silicate: an octahedral Al-dominated sheet
(O) sandwiched between two tetrahedral sheets (T), where
SiO4⁻ tetrahedra are partially substituted by AlO3⁻ and
the excess of negative charge is compensated by the pres-
ence of potassium cations in the interlayer. Radiogenic
⁴⁰Ar is assumed to occupy the position of parent potassium,
trapped mechanically between the TOT layers (Dong et al.,
1995; Dahl, 1996; Sletten and Onstott, 1998). Therefore, the
radiogenic ⁴⁰Ar retention in the structure is attributed
exclusively to the electrostatic cohesion forces between the
adjacent negatively charged TOT layers mediated by the
positively charged interlayer K⁺ cations, OH group
orientation, and the interlayer ionic porosity (interlayer
unit-cell volume not occupied by cations; Dahl, 1996).

Several decades of extensive research have proven the methods of K–Ar dating and its ⁴⁰Ar/³⁹Ar analogue to be successful and accurate tools for dating various alterations using mixed-layered illite–smectite clay. There is a general consensus that the difference in the measured K–Ar ages between fine- and coarse-grain illite particles can be explained using detrital-authigenic components system (e.g., Środoń, 1999; van der Pluijm et al., 2001; Elliott and Haynes, 2002; Clauer et al., 2012). In the case of K–Ar age variability among finest subfractions that should represent purely authigenic clay (i.e., <0.2 µm in ben tonites), their crystallization history and illite–smectite crys tallite growth may play a major role (Clauer et al., 1997; Środoń et al., 2002; Osborn et al., 2014). On the other hand, Lerman and Clauer (2005) and Lerman et al. (2007) have suggested post-crystallization diffusion in the diagenetic conditions as a mechanism responsible for lower K–Ar iso tope age estimates of finer particles, as compared to coarser ones. However, the diffusion model is not consistent with the experimental evidence: Hassanipak and Wampler (1996) did not find any increased Ar loss in the finer grain size fractions in preheated illite and glauconite samples. The Ar release pattern found by Hassanipak and Wampler (1996) was clearly related to the dehydroxylation reaction of clays. Nevertheless, even if the diffusion model of Lerman et al. (2007) is not correct, their concern is still valid: could there be other reasons for the ⁴⁰Ar loss in the diagenetic temperature zone over the geological time span?

While the K–Ar methodology measures K and radio-
genic ⁴⁰Ar concentrations, ⁴⁰Ar–³⁹Ar method involves prior
neutron irradiation in order to transform ³⁹K into ³⁹Ar for
⁴⁰Ar/³⁹Ar ratio measurement (Dalrymple and Lanphere,
1971). During the irradiation-induced ³⁹K to ³⁹Ar trans-
formation, high recoil energy of ³⁹Ar, ~177 keV (Onstott
et al., 1995), allows ³⁹Ar penetration through the clay struc-
ture and its escape from the fine crystallites (Smith et al.,
1993; Dong et al., 1995). Due to the high ³⁹Ar recoil energy
and the resulting substantial ³⁹Ar loss, the application of
the ⁴⁰Ar/³⁹Ar method to clay minerals requires using a
microampoule technique of encapsulation (Dong et al.,
1995; Onstott et al., 1997), and thus the use of step-heating

and plateau-age methods seem problematic (Clauer et al.,
2012). The fraction of ³⁹Ar lost from the illite structure dur-
ing irradiation is a function of its crystallite thickness and
grain size equivalent (Dong et al., 1995; see Clauer et al.,
2012 for extensive discussion), suggesting that ³⁹Ar recoil
not only stimulates the ³⁹Ar nuclei diffusion along the inter-
layer, but also its penetration through the TOT layers.

If the high-energy irradiation-induced ³⁹Ar recoil leads to the ³⁹Ar loss from fine illite particles, perhaps the natural ⁴⁰Ar recoil can also result in some ⁴⁰Ar loss, lower than in the case of ³⁹Ar, but still significant enough to provide observable decrease of K–Ar age in the finest illite crystallites? A random walk model of Hall et al. (2000) assumes no ⁴⁰Ar loss from fine illite crystallites, which serves as reference for the irradiation-induced ³⁹Ar loss. Despite numerous successful applications of the model in ⁴⁰Ar/³⁹Ar dating of illite, the assumption about ⁴⁰Ar has not been independently validated and the ⁴⁰Ar loss, if occurring, may be covered by a strong effect of ³⁹Ar loss, thus not detectable by a model calculations.

The illite–smectite crystallites commonly have thick-
nesses of their fundamental illite particles from two to sev-
eral 2:1 layers (Nadeau et al., 1984; Środoń et al., 1992;
Dudek et al., 2002). If the ⁴⁰Ar recoil energy can make some
⁴⁰Ar atoms penetrate through a single 2:1 layer, the broadly
accepted assumption that all radiogenic ⁴⁰Ar remains
mechanically trapped in the interlayer may be incorrect.
It is also not clear whether ⁴⁰Ar after recoil can be incorpo-
rated somewhere else into the clay structure or its recoil
energy is not sufficient to break the chemical bonds of the
clay layer, so ⁴⁰Ar can remain only in the crystallographic
positions of potassium. Similar questions were hypothe-
sized earlier by Goncharov et al. (2007). Hetherington
and Villa (2007) argued that a minor portion of ⁴⁰Ar in
coarse muscovite can occur in a position different than
the interlayer: within the T–O–T framework. Due to the
recent developments of molecular dynamics methodology
it is now possible to address these questions quantitatively
and to draw specific conclusions relevant for K–Ar geologi-
cal age estimates.

## 2. METHODS

**2.1. Molecular dynamics simulations**

In order to calculate the effect of argon recoil, the kinetic energy and velocity of ⁴⁰Ar after recoil were calculated on the basis of the energy of emitted gamma quantum after the ⁴⁰K–⁴⁰Ar electron capture transition: 1460.822 keV (Bé et al., 2010). This transition produces a photon with a mass of 1.56825 × 10⁻³ amu. Employing the principle of momentum conservation we can calculate the ⁴⁰Ar atom velocity as 0.11765 Å/fs, which corresponds to the kinetic energy of 28.66 eV (the same value as given by Dong et al., 1995). All the simulations were performed for the same initial velocity but with varying recoil angles (Fig. 1) that randomly cover all possible recoil angles within a hemisphere. The angle θ is defined as an angle between the recoil vector and the direction perpendicular to the clay layer (“vertical”). The angle φ is defined as an angle between
<figure><img src="image_2.png" /><figcaption>Fig. 1. Definition of the recoil angles $\theta$ and $\varphi$. Green – $^{40}$Ar, red – oxygen, yellow – $SiO_4$ tetrahedrons, violet – $AlO_6$ octahedrons. Black arrow indicates the recoil vector, while the crystallographic direction *a* is indicated by another black arrow with letter *a* above it. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the crystallographic direction *a* and projection of the recoil vector on the *ab* plane of the clay crystal (“horizontal”; Fig. 1).

Due to limitations of the available force fields for molecular simulation of clay materials, two different approaches were used. An explicit atomistic model of illite fundamental particles was built out of two TOT layers and simulated using the CLAYFF force field (Cygan et al., 2004). Because CLAYFF was developed as a non-reactive force field for an undistorted clay structure (Cygan et al., 2004), this approach was used only to study the $^{40}$Ar recoil effect for angles close to the *ab* plane ($\theta \approx 90^{\circ}$). In this case the interaction of Ar with the TOT layer is negligible and interlayer $K^+$ ions play the crucial role. Reliability of such an approach is based on the energy of interaction between interlayer $K^+$ and the TOT layer. The validity of the model parameters is supported by the calculations of Heinz et al. (2005) that show that cleavage energies for mica, and, therefore, the attraction forces between interlayer potassium and the TOT layer, calculated with CLAYFF are relatively close to the experimental values.

In the second approach, the interactions between the TOT layer and the Ar atom were tested for the range of angles $\theta < 80^{\circ}$ while the interactions with the interlayer potassium were assumed to be negligible and were tentatively ignored. Bonds breaking, structure deformations, and Ar penetration into the TOT layer, were assumed as dominant reactions at these low $\theta$ angles, and were investigated using an implementation of the reactive force field REAXFF (van Duin et al., 2001) in the LAMMPS molecular simulation program (Plimpton, 1995). Contrary to the classical non-reactive force fields, such as CLAYFF, REAXFF allows modeling chemical reactions in a semi-empirical way, and does not require explicit definition or preservation of all interatomic bonds. In our REAXFF simulations one layer structure of trans-vacant pyrophyllite was used as the best available approximation of the mica/illite TOT layer structure. This approximation was necessary because the REAXFF parameterization for potassium is not yet available. Moreover, potassium can be ignored at all in TOT–Ar interactions at low $\theta$ angles, thus justifying the use of pyrophyllite layer as a fully valid model. In our REAXFF approach, the force field parameters optimized for clay-zeolite structures (Pitman and van Duin, 2012)

were used along with the parameters for noble gases (van Duin, personal communication). These simulations employed a single TOT layer of pyrophyllite structure $Al_2Si_4O_{10}(OH)_2$, consisting of 8 × 4 unit cells in the *a* and *b* crystallographic directions, respectively. All crystallographic octahedral positions were filled only with Al atoms, and all the tetrahedral positions – with Si atoms. Except for the electrostatically neutral $^{40}$Ar atom in the position of interlayer potassium as in illite/mica structure, no other interlayer species were allowed. Prior to the simulation, the layer structure was optimized allowing to relax the crystallographic dimensions *a* and *b*. The recoil angles were set to cover the entire range of statistical possibilities of recoil toward the TOT layer. In the REAXFF simulations the vertical angle $\theta$ was set in the range from $0^{\circ}$ to $80^{\circ}$ with a step of $10^{\circ}$, while the horizontal angle $\varphi$ was set in the range from $0^{\circ}$ to $360^{\circ}$ with a step of $20^{\circ}$ (Fig. 1). All the possible angle combinations resulted in the total of 145 individual MD simulation runs.

The time step during REAXFF simulations was set to 0.25 fs, and the MD runs were performed for 12.5 ps each (after sufficient pre-equilibration for structure optimization) using the NVE-ensemble with the initial velocity of the selected $^{40}$Ar set at 0.11765 Å/fs. It is important to note, that all simulations were intrinsically non-equilibrium, because the additional kinetic energy was suddenly introduced into the system and eventually redistributed into the clay structure by increasing a little the potential energy of all atoms in the system due to the lattice deformation. In spite of the simplification (one layer and no interlayer potassium), such a model is expected to give relatively clear picture of the structural rearrangements with general conclusions applicable to illite fundamental particles. It is because all of the studied recoil vectors were directed toward the TOT layer neighboring the interlayer $^{40}$Ar atom and interactions with interlayer potassium are of much less importance here. For the angles $\theta$ higher than $60^{\circ}$, the REAXFF model might give results more strongly dependent on the neighboring interlayer potassium cations, which is in turn simulated using the CLAYFF models. However, as will be shown below, in such a case the $^{40}$Ar usually does not penetrate the TOT layer and the resulting structure remains mostly undistorted.
In order to study the recoil distance of the <sup>40</sup>Ar atom at θ angles higher than 60°, i.e. along the illite interlayer, with the CLAYFF force field, a model consisting of two TOT layers with interlayer potassium was built. The chemical composition was set to mimic a hypothetical structure of an illite fundamental particle build of two 2:1 layers and having 50% of smectite as outer surfaces (Środoń et al., 2009). For simplicity, the number of positions filled with interlayer potassium was larger than in natural illites (mica composition). Isomorphic Al/Si and Mg/Al substitutions (in the tetrahedral and octahedral sheets, respectively) were distributed randomly, avoiding direct neighboring of two Mg and two Al<sup>IV</sup> substituted sites. The recoil effect in different directions (angles φ from 0° to 60°, varying by a step of 10°) was studied assuming that the interlayer positions surrounding the tested <sup>40</sup>Ar atom were either all occupied by potassium cations or one position was left vacant (Fig. 2). However, H<sub>2</sub>O molecules may occur in the K-vacant positions (Vidal et al., 2010). In order to further test the fate of the <sup>40</sup>Ar atom after recoil at the edges of crystallites, calculations assuming <sup>40</sup>Ar-neighboring vacancies (dashed black lines in Fig. 2) were also performed to mimic the crystallite fringe edge with angles φ varying from -60° to 120°, with a step of 10° and θ from 60° to 90°, with a step of 10°. Although the CLAYFF force field has not been designed to take into account crystal edges (Cygan et al., 2004) such a trick allows us to model the Ar recoil at crystallite edges and fringes. All CLAYFF simulation runs were performed similarly to the REAXFF runs described above. Prior to each MD run, the system structure was carefully optimized by energy minimization procedure with energy tolerance, i.e., the energy change between successive iterations divided by the total energy magnitude, of 10<sup>-6</sup>.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Definition of the recoil angle φ in plane and positions of an assumed K-free vacancy. Numbers indicate the location of potassium cations removed in different simulations. In other series of calculations linear set of neighboring potassium ions marked with dashed black lines were removed to mimic a crystallite edge environment.</figcaption></figure>

### 3. RESULTS

#### 3.1. Possible recoil reactions

The recoil energy of <sup>40</sup>Ar is 28.66 eV, which is equal to 660.9 kcal for 1 mol of <sup>40</sup>Ar. This energy is much higher than the typical dissociation energies of covalent bonds in clay minerals: Si-O – 148 kcal/mol, Mg-O – 94 kcal/mol, Fe-O – 98 kcal/mol and Al-O – 122 kcal/mol (Cottrell, 1958; Darwent, 1970). A bond can be broken if the energy transfer leads to the bond elongation in a particular direction beyond a certain length threshold. Therefore, based on the observed behavior of the recoiled Ar atoms for different recoil angles, different effects on the silicate layer structure may occur. They can be described from a perspective of either the final position of <sup>40</sup>Ar in the clay structure or the final residual deformation of the clay structure.

#### 3.2. Recoil effects from the perspective of final position of <sup>40</sup>Ar in the structure (REAXFF)

Four possible final positions of <sup>40</sup>Ar can be identified at the end of the REAXFF simulation when the initial structural perturbation is dissipated (Fig. 3). For θ > 40°, the most common effect observed is that the <sup>40</sup>Ar atom bounces back from the basal oxygen surface and returns to its initial position in the interlayer. Depending on the recoil angle, some deformations of the structure can accompany this process with the breaking and formation of new chemical bonds. Ar bouncing back also takes place even for θ < 40°, if the direction of the recoil vector is pointing to a Si atom of one of the SiO<sub>4</sub> tetrahedra. No bond breaking was then observed (Fig. 3a, b; Supplementary video 1).

<figure><img src="image_4.png" /><figcaption>Electronic Annex: Video 1</figcaption></figure>

The most common effect observed for θ < 40° is the penetration of <sup>40</sup>Ar into the tetrahedral sheet, which occurs along with the deformation (displacement) of the neighboring OH groups. Most often, the OH group that was initially the closest to the <sup>40</sup>Ar, changes its position forming new bonds with other atoms of the TOT layer (usually on the opposite side to the initial Ar location). Therefore, <sup>40</sup>Ar
can occupy the empty cavity after the OH group dissociated from the tetrahedral sheet (Fig. 3c, d; Supplementary video 2).

<figure><img src="image_5.png" /><figcaption>Electronic Annex: Video 2</figcaption></figure>

More precisely, <sup>40</sup>Ar finds its new position close to the position of the displaced OH group. Depending on the extent of the OH group displacement, <sup>40</sup>Ar is finally located closer to the surface or deeper within the layer and closer to the octahedral position.

Only for three combinations of $\theta$ and $\phi$ angles we observe the penetration of the octahedral sheet by the <sup>40</sup>Ar atom upon recoil. In one case, <sup>40</sup>Ar is recoiled directly into the vacancy in the octahedral sheet, avoiding any energy transfer to the OH group. In two other cases, <sup>40</sup>Ar bounces from the tetrahedral basal oxygen and the octahedral aluminum to eventually end up in the octahedral vacancy (Fig. 3e, f; Supplementary video 3).

<figure><img src="image_6.png" /><figcaption>Electronic Annex: Video 3</figcaption></figure>

In these three cases the recoil is also accompanied by some deformation of the structure that includes only small distortion of the OH groups.

Under certain recoil angles pointing approximately toward the vacancy in the octahedral sheet, the <sup>40</sup>Ar atom can pass through the octahedral vacancy and embed into the tetrahedral sheet on the opposite side of the TOT layer (Fig. 3g, h; Supplementary video 4).

<figure><img src="image_7.png" /><figcaption>Electronic Annex: Video 4</figcaption></figure>

This effect is accompanied by a deformation of the neighboring OH groups, with bond breaking and the formation of new bonds between the OH groups and the octahedral or tetrahedral species. However, none of the recoil angle combinations allowed the <sup>40</sup>Ar atom after recoil to pass completely through the TOT layer to the interlayer on the opposite side.

### 3.3. Layer structure deformation due to <sup>40</sup>Ar recoil (REAXFF)

Although the major layer structure deformation must result from a $\beta$-decay (⁴⁰K → ⁴⁰Ca) that occurs ~89.3% of the time, with energy similar to the electron capture (1.33 MeV), understanding the behavior of ⁴⁰Ar after recoil and its recoil pathway requires determining the 2:1 layer structure deformations related exclusively to ⁴⁰Ar recoil. The possible effects of ⁴⁰Ar recoil on the deformation of the clay structure vary depending on the recoil angle. The most common reaction at low vertical angles ($\theta$ < 40°) is associated with kinetic energy transfer to the closest OH group and a subsequent recoil of this group. Depending on the ⁴⁰Ar recoil angle, the following reactions were observed:

(1) OH group was recoiled but still remained connected to an aluminum atom of the octahedral sheet (Fig. 3c, d). The newly formed structure had some octahedral Al atoms in defect penta-coordination and/or hepta-coordination.

(2) The recoiled OH group created new bonding with another OH group, resulting in the formation
of a water molecule and a residual oxygen (Fig. 4a, b; Supplementary video 5).

<figure><img src="image_8.png" /><figcaption>Electronic Annex: Video 5</figcaption></figure>

The H₂O molecule can then behave in three different ways:

- leave the TOT structure (as in Fig. 4a, b);

- remain as a weakly-bonded H₂O molecule within the TOT layer;

- form a new bonding with two Al atoms of the octahedral sheet (oxygen of this H₂O molecule is at the position of the former oxygen of the OH group, Fig. 4c, d; Supplementary video 6).

<figure><img src="image_9.png" /><figcaption>Electronic Annex: Video 6</figcaption></figure>

After the water molecule is formed, <sup>40</sup>Ar can either remain in the TOT layer structure or bounce back to the interlayer space. Depending on the initial relative positions of the O and H atoms to which the <sup>40</sup>Ar recoil energy was transferred to, different final structures can be created. A configuration resembling the locally dehydroxylated structure (Drits et al., 1995) was formed after Ar struck the

two OH groups belonging to the same AlO₄(OH)₂ octahedron (Fig. 4a, b). Alternatively, if the OH groups shared the same vacancy, the hydrogen of one of the OH groups was attached to another OH group and the remaining oxygen was bound to silica tetrahedrons after the recoil (Fig. 4c, d).

(3) The recoiled OH group created new bonds with Si atoms of the tetrahedral sheet. In this newly formed structure the recoiled OH group could form bonds with one or two Si atoms, leading to the formation of penta-coordinated Si sites with some bonds visibly elongated compared to the bonding in tetra-coordinated silicon (Fig. 4e, f; Supplementary video 7).

<figure><img src="image_10.png" /><figcaption>Electronic Annex: Video 7</figcaption></figure>

Depending on the recoil angle, the OH group could be reattached either within the TOT layer or on its outer surface.

At larger incoming angles (θ ≥ 50°), the recoiled <sup>40</sup>Ar atom interacts not with the OH group located in the center of the hexagonal cavity, but with the tetrahedral sheet. Three effects were then observed, each associated with the <sup>40</sup>Ar bouncing back to the interlayer:

(1) If the recoil vector is pointing to a Si atom of the tetrahedron, the <sup>40</sup>Ar atom bounces back from the tetrahedron and the TOT layer structure stays intact (Fig. 3a, b);

(2) If the recoil vector is pointing to a basal O atom, the energy transfer from <sup>40</sup>Ar to that O atom at the induced maximal deformation makes the H atom of the neighboring octahedral OH group become attached to that O atom and the residual O remains in the octahedral sheet at the position of the former OH group.

(3) The <sup>40</sup>Ar recoil leads to a substantial deformation of the tetrahedral sheet, followed by the formation of new OH groups attached to Si tetrahedra, and
accompanied by simultaneous bond breaking and
reformation in the octahedral sheet (Fig. 5;
Supplementary video 8).

<figure><img src="image_11.png" /><figcaption>Electronic Annex: Video 8</figcaption></figure>

Bouncing of the incoming <sup>40</sup>Ar atom back from the TOT layer structure without any bond breaking appears uncommon and occurs only at certain angles. This effect is generally observed for angles θ > 60°. If the recoil vector is pointing to a Si atom of a SiO<sub>4</sub> tetrahedron, the recoil energy of the <sup>40</sup>Ar atom is distributed among all four Si–O bonds and other neighboring parts of the TOT structure, and no bond breaking is observed (Fig. 3a, b). There were also four cases for θ ≈ 30° where <sup>40</sup>Ar penetrated deeply into the structure but the recoil energy was distributed without any bond breaking. This was, however, exceptional, because in most cases for low angles θ irreversible deformations were observed within the simulation time of our studies.

**3.4. Spatial distribution of the effects of <sup>40</sup>Ar recoil on the**
**TOT layer structure and final Ar position**

Quantitative estimation of all possible consequences of
$^{40}$Ar recoil, based on all REAXFF simulations, as a func-
tion of angles θ and φ (in a hemispherical projection) is
necessary in order to assess the $^{40}$Ar behavior over geologi-
cal time span and possible effects of this process on K–Ar
dating of clay minerals. A plot depicting the recoil effect
for different positions of the intersection of the $^{40}$Ar trajec-
tory with the TOT layer surface (a plane formed by the cen-
ters of basal oxygen atoms), is shown in Fig. 6a, b.

The final position of $^{40}$Ar after recoil (Fig. 6a) varies
depending mostly on the angle $\theta$, but some dependence
on the angle $\varphi$ is also observed. For $\theta \le 40^\circ$, $^{40}$Ar is most
often positioned in the tetrahedral sheet. Combining this
information with the data of Fig. 6b, which depicts the dis-
tortions of the clay structure, it becomes clear that placing
$^{40}$Ar in the tetrahedral sheet is always coupled with the
deformation of OH groups. The displaced OH group is usu-
ally bonded simultaneously to the tetrahedral and

octahedral sheet, and very often more than one OH group
is distorted. This substantial qualitative difference both in
$^{40}$Ar location and in the deformation of the TOT layer
structure between the recoil effect simulated for $\theta > 40^\circ$
and $\theta \le 40^\circ$ angles is due to the fact that in the latter case
$^{40}$Ar first collides with the OH group and the recoil results
in its displacement, while in the former case (i.e., for higher
$\theta$) $^{40}$Ar initially hits the basal oxygen surface of the tetrahe-
dral sheet and the recoil usually leads to the deformation of
the tetrahedral sheet. For $\theta \approx 70^\circ$, a reattachment of the
knocked off hydrogen to a basal oxygen was observed,
but the lack of distortion of the TOT layer structure is
the dominant situation.

The calculation of the total probability for $^{40}$Ar to be
incorporated into the TOT layer structure, gives the value
of 20.5%, with a 16.5% chance for its incorporation into
the tetrahedral sheet, 1.5% – into the octahedral sheet,
and 2.5% – into the tetrahedral sheet on the opposite side
of the TOT layer. There is also a 5% probability to form
a water molecule during this process. All these values were
calculated from the expression of a spherical dome as frac-
tions of the sphere surface corresponding to each combina-
tion of the two recoil angles, $\theta$ and $\varphi$.

**3.5. Maximum recoil distance of <sup>40</sup>Ar along the interlayer**
**(CLAYFF; θ close to 90°)**

As shown in the previous sections, if the angle θ
approaches 90°, the effect of $^{40}$Ar recoil on the TOT struc-
ture is negligible and $^{40}$Ar never penetrates the TOT layer.
The CLAYFF simulations of the $^{40}$Ar recoil distance along
the interlayer for a model consisting of two TOT layers
(Fig. 2) show that if there is no vacancy next to the initial
$^{40}$Ar position (i.e., in the positions marked as 1, 2, and 3
in Fig. 2), the recoiled $^{40}$Ar atom always bounces back from
the neighboring K⁺ and returns to its initial interlayer posi-
tion, independent of the recoil angle φ.

If the interlayer vacancy is an immediate neighbor of the
$^{40}$Ar atom (position 1 on Fig. 2), the recoiled $^{40}$Ar ends up
occupying the vacancy for most of the studied angles
($0 \le \varphi \le 60^\circ$). Only for $\varphi \approx 0^\circ$ (or, by analogy, for
$\varphi \approx 120^\circ$) $^{40}$Ar bounces back to the initial position. For
any vertical angle $\theta < 60^\circ$ at $\varphi = 30^\circ$, $^{40}$Ar bounces back
and returns to its initial position.

For a vacancy located in the position of a second neigh-
bor (position 2 on Fig. 2), the range of angles φ leading to
an eventual jump of the recoiled $^{40}$Ar to the vacant position
is approximately between 20° and 50°. However, if $\theta < 80^\circ$,
$^{40}$Ar bounces back and returns to its initial position regard-
less the φ angle.

If a vacancy is located in the farther neighboring posi-
tion 3 on Fig. 2, the only effect observed was pushing the
K⁺ cation from the position 2 to the position 3, while the
$^{40}$Ar ends up embedded in the position 2, initially occupied
by that potassium ion. This effect was observed for the
range of angles $20^\circ < \varphi < 40^\circ$. For a vacancy located in the
position 4 of Fig. 2, the $^{40}$Ar always bounces back to its ini-
tial position independently of the initial recoil angle φ.

The results of simulations performed for linear series of
neighboring vacancies as a model mimicking the crystallite
<figure><img src="image_12.png" /><figcaption>Fig. 3. Recoil effects resulting in different final position of the <sup>40</sup>Ar atom in the structure (final simulation frames): (a, b) the largest deformation amplitude when <sup>40</sup>Ar atom bounces from the TOT layer without deforming the structure. (c, d) <sup>40</sup>Ar atom is embedded in the tetrahedral sheet. (e, f) <sup>40</sup>Ar atom is embedded in the octahedral sheet. (g, h) <sup>40</sup>Ar atom is embedded in the tetrahedral sheet on the opposite side of the TOT layer.</figcaption></figure>

fringe edge show substantial dependence on angles $\theta$ and $\phi$. For $\theta < 80^{\circ}$, <sup>40</sup>Ar located behind the row of potassium cations bounces back and returns to its initial position regardless of the value of $\phi$. For $\theta$ equal to 80° or 90° the fate of <sup>40</sup>Ar is very similar and with important dependence on the recoil angle $\phi$ (Fig. 7). In some cases, if <sup>40</sup>Ar recoil

vector is pointed to the center of a potassium ion (1a, 1b and 2b in Fig. 7) Ar bounces back to its initial position. There are, however, some values of $\phi$ for which <sup>40</sup>Ar ends up outside of the “crystallite” by passing between either two (positions 2b and 2c in Fig. 7) or three (positions 1a, 2a and 2b in Fig. 7) potassium atoms. There is also a
<figure><img src="image_13.png" /><figcaption>Fig. 4. Recoil effects depending on the final deformation of the TOT layer structure due to the recoil of the initially closest OH group to the incoming $ ^{40}\text{Ar} $. A final simulation frames represent: (a, b) a new $H_2O$ molecule is formed and ejected from the TOT sheet, while $ ^{40}\text{Ar} $ is bounced back from the TOT surface; (c, d) a new $H_2O$ molecule is formed in the position of one of the initial OH groups, $H_2O$ oxygen forms two penta-coordinated Si sites, while $ ^{40}\text{Ar} $ atom is embedded in the tetrahedral sheet; (e, f) the OH group forms two penta-coordinated Si sites, while$^{40}\text{Ar}$ atom is bounced back from the oxygen surface.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 5. Recoil effect under high angle $\theta$ leading to the deformation of the tetrahedral sheet and, subsequently, of the octahedral sheet. Ar (not shown) has bounced back to the interlayer.</figcaption></figure>

substantial fraction of simulations in which $ ^{40}\text{Ar} $ ended up at the direct edge of the “crystallite” in the position formerly occupied by potassium (which, in turn, is pushed

out of the structure) – at position 2b, 2c and 2d in Fig. 7. In one case $ ^{40}\text{Ar} $ switched its position with a neighboring potassium ion (position 1b in Fig. 7). The “crystallite edge”
<figure><img src="image_15.png" /></figure>

Fig. 6. Projection of the hemisphere of recoil angle combinations on the crystallographic *ab* plane presenting a map of the recoil effects. Left: final position of the <sup>40</sup>Ar atom in the structure: orange – <sup>40</sup>Ar in the tetrahedral sheet on the same side of the TOT layer, yellow – <sup>40</sup>Ar in the octahedral sheet, green – <sup>40</sup>Ar bounced back from the TOT surface, brown – <sup>40</sup>Ar in the tetrahedral sheet on the opposite side of the TOT layer. Right: various types of structural distortions due to recoil: red – misplaced OH reconnected to the tetrahedral and/or octahedral sheet, light blue – newly formed H<sub>2</sub>O molecule stays in the TOT layer, blue – newly formed H<sub>2</sub>O molecule leaves the TOT structure, orange – distortion of the tetrahedral sheet, yellow – hydrogen attached to a basal oxygen; green – the structure remained undistorted. The distance from the center of the circle characterizes the *θ* angle. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

simulation results extrapolated to illite/smectite crystallites indicate that if <sup>40</sup>Ar is located in the position of a second neighbor to the edge, it will end up out of the crystallite in 3.8% of all recoils (normalized to the hemispherical projection) and will end up at the edge of a “crystallite” in 4.8 % of all cases.

**3.6. Potential energy of <sup>40</sup>Ar leaving the tetrahedral and octahedral sheet**

Fig. 8 shows the variation of the short-term moving average value of the potential energy in two cases:

- <sup>40</sup>Ar passing through the hexagonal (ditrigonal) cavity, in the direction from the tetrahedral sheet to the interlayer (Fig. 8a).

- <sup>40</sup>Ar passing the substantial part of TOT layer from octahedral position within TOT layer to the interlayer (Fig. 8b).

Average values of the energies were calculated for periods of 400 fs sliding along the entire generated MD trajectory. The $f_x$, $f_y$ and $f_z$ components of the force acting on the <sup>40</sup>Ar atom were set to zero and its velocity was set to a constant value – hundred times smaller than the velocity of the recoil. Additionally, to keep the clay layer in position, the $z$ coordinate for one of the oxygen atoms in the structure was fixed. In the first case the initial structure corresponded to a particular case of recoil with $\theta = 0$, i.e., the recoil vector pointing perpendicularly to the TOT layer. In the second case the structure corresponded to <sup>40</sup>Ar located in the octahedral sheet with the initial recoil angles $\theta$ and $\varphi$ set to 20°.

The average potential energy of the structure with <sup>40</sup>Ar located in the tetrahedral sheet is higher by about 26 kcal/-mol compared to the case when <sup>40</sup>Ar is eventually located out of the TOT structure (Fig. 8a). The energy barrier for

<figure><img src="image_16.png" /></figure>

Fig. 7. Map of the recoil effects along the interlayer ($\theta = 90°$). Solid arrows – <sup>40</sup>Ar transferred to the pointed positions, short lines – <sup>40</sup>Ar returned to the initial position, dashed arrows – marked $K^+$ ions displaced to the pointed positions.

<sup>40</sup>Ar passing through the hexagonal cavity is ≈17 kcal/mol. For comparison, the recoil energy of <sup>40</sup>Ar is equal to 660.9 kcal/mol, which indicates that <sup>40</sup>Ar can pass easily through the hexagonal cavity at the angles close to $\theta = 0$.

In the case of the structure with <sup>40</sup>Ar located in the octahedral sheet, the average potential energy is higher of about ≈30 kcal/mol compared to a position outside of the TOT structure (Fig. 8b). The energy barrier of <sup>40</sup>Ar passing through the structure is much higher than for the situation of passing only through the hexagonal cavity, i.e., ≈55 kcal/mol. This is related to the fact that on leaving the octahedral sheet structure <sup>40</sup>Ar needs to deform it substantially and to pass between an OH group and one of the basal oxygens, and then also through hexagonal (ditrigonal) cavity.
<figure><img src="image_17.png" /><figcaption>Fig. 8. Moving averages of the potential energy of (a) $ ^{40}\text{Ar} $ passing through the hexagonal cavity from an initial location in the tetrahedral sheet after the recoil perpendicular to the surface; (b) $ ^{40}\text{Ar} $ passing the TOT layer from initial octahedral position.</figcaption></figure>

# 4. DISCUSSION

## 4.1. Deformation of the clay structure associated with $ ^{40}\text{Ar} $ recoil in natural conditions

The extent of structural deformation clearly depends on the recoil angle. The deformation for low $\theta$ angles mostly involves the recoil of the OH groups. The formation of an H$_2$O molecule and its passage through the hexagonal cavity found in this study is analogous to the OH removal during dehydroxylation, as evidenced by extensive modeling studies (Sainz-Díaz et al., 2004; Molina-Montes et al., 2008). The OH removal reaction should be considered as reversible, as it is shown by the dehydroxylation-rehydroxylation experiments (Derkowski et al., 2012) and by molecular modeling (Molina-Montes et al., 2010). A dehydroxylated diocahedral 2:1 layer resorbs water molecules from the environment and uses them to reconstruct the OH groups, thus leading to structural rearrangements. In the case of the studied system, because of the presumed isolation from ambient sources of water, either the OH/H$_2$O removed by recoil or the interlayer H$_2$O that may occur in the positions unoccupied by K$^+$ in K-deficient mica (Vidal et al., 2010) are the potential sources of water used to rehydroxylate the TOT layer structure. Therefore, it can be also assumed that under real geological conditions and time span, the structural deformation caused locally by $ ^{40}\text{Ar} $ recoil leading to the deformation of OH groups will be eventually healed back to the initial structure. The reconstruction occurs due to the fact that all the deformed structures are thermodynamically less stable than the initial state.

## 4.2. Probability of retaining the $ ^{40}\text{Ar} $ atom in the layer structure over geological time periods

Taking into account the energy barrier of $ ^{40}\text{Ar} $ passing through the hexagonal cavity of the 2:1 clay structure, which is ca. 17 kcal/mol, and $ \Delta E $ of this reaction (-26 kcal/mol) there is hardly any possibility of $ ^{40}\text{Ar} $ staying within the tetrahedral sheet of the TOT layer structure during a time span of geological scale. $ \Delta E $ of the reaction was calculated for the $ ^{40}\text{Ar} $ transfer from TOT layer to vacuum; it will be lower if $ ^{40}\text{Ar} $ moves to the interlayer space due to a partial overlap of electron clouds of Ar

and TOT atoms. However, it will still be an exothermal reaction because of the substantial difference in the steric constrains affecting the $ ^{40}\text{Ar} $ atom, when comparing the interlayer space environment and the interior of the TOT layer. The low energy barrier can be overcome in low-temperature diagenetic reactions in geological environment. For example, a similar range of activation energies has been found for minerals weathering and dissolution at temperatures corresponding to a land surface or shallow burial before the onset of illitization ($t < 60$ °C; Plummer et al., 1978; Palandri and Kharaka, 2004). The activation energy for diffusion of carbon atoms in ferrite is a little higher, 20.1 kcal/mol, while the reaction occurs at temperatures close to ambient (Cermák and Král, 2010). The calculated energy barrier thus allows predicting that all $ ^{40}\text{Ar} $ atoms placed after recoil in the tetrahedral sheet will be eventually moved back to the interlayer. Rehydroxylation, mentioned in Section 4.1, may also play a role in the removal of $ ^{40}\text{Ar} $ from the tetrahedral sheet back to the interlayer.

One important conclusion based on the simulations is that there is 2.5% chance of $ ^{40}\text{Ar} $ to end up recoiled to the tetrahedral sheet on the other side of the TOT layer. If the other side of the TOT layer is a smectitic, external surface of illite crystallite, with no strongly fixed cations, within geological time the $ ^{40}\text{Ar} $ embedded in such a position can escape from the TOT layer structure, leading to the decrease of the apparent K–Ar age. The effect of Ar loss will depend on the thickness of illite crystallites ($c^*$ dimension). It will be negligible for thick mica crystals, but for the illite-smectite crystallites consisting of just a few TOT layers the effect can be significant.

On the other hand, there is non-negligible probability of ≈1.5% that $ ^{40}\text{Ar} $ ends up within the octahedral sheet. Much higher energy barrier of 55 kcal/mol indicates that after entering the octahedral sheet, $ ^{40}\text{Ar} $ will not be easily removed due to thermal vibrations and therefore will be kept within TOT layer over geological time.

## 4.3. Crystal size and crystallite thickness vs. the effect of $ ^{40}\text{Ar} $ loss after recoil

In order to apply the obtained MD simulation results to the model of illite particles and thus to estimate potential
effects of $^{40}$Ar recoil on the K-Ar dating of these minerals,
the following assumptions were made:

- illite crystallites have hexagonal shape and consist of a relatively small number of TOT layers (Środoń et al., 1992);

- $^{40}$K atoms that decay to $^{40}$Ar are distributed randomly;

- $^{40}$Ar can be lost due to recoil if its position before recoil is close enough to the crystallite side edges (external position in Fig. 9); $^{40}$Ar can escape from the structure also due to thermal motions if its original position is right at the edge of a crystallite (external position in Fig. 9). In the geochronologic calculations, the sum of these two effects was taken into account because they are undistinguishable over a geologic time span.

No direct loss of $^{40}$Ar during recoil was assumed along
the interlayer for Ar atom initial locations in the internal
part of the crystallites (blue color in Fig. 9). If $^{40}$Ar was
located at the most external part of the crystallite (red color

<figure><img src="image_18.png" /></figure>

Fig. 9. Different positions of $^{40}$Ar probed to study the $^{40}$Ar loss after recoil. Blue – internal, yellow – next to external, red – external part of the crystallite. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

in Fig. 9) substantial amount of $^{40}$Ar can be lost due to the
recoil outside of the layer. The exact value of this loss from
recoil has not been determined, because the most external
parts of a crystallite are substantially different in their
structure and composition from the interior due to the for-
mation of OH groups on the crystallite edges and the weak-
ening of attraction between the neighboring layers (e.g.,
Churakov, 2007; Liu et al., 2012; Tazi et al., 2012). In the
next-to-external position in the layer (yellow color in
Fig. 8), the $^{40}$Ar recoil model is analogous to the situation
for the series of K-vacancies (Fig. 7). In the illite structure
not all interlayer positions are occupied by potassium.
These vacant positions are, however, usually considered
as occupied by H$_2$O molecules (Vidal et al., 2010). In a very
approximate calculation, water molecules were considered
to block the escape of $^{40}$Ar similarly to potassium ions.
Taking into account the results from Section 3.5 above,
the loss of $^{40}$Ar from these positions due to recoil will hap-
pen for ≈3.8% of all cases.

$^{40}$Ar escape from the TOT structure due to thermal motions can happen both from the tetrahedral sheet on the external side of the crystallite, where the recoiled $^{40}$Ar was embedded on the other side of the TOT layer, and from crystallite side edges. In the first case, this loss would account for ≈2.5% of $^{40}$Ar produced, if the illite crystallite consists of only two TOT layers. For particles with more layers, this value would correspond to the ratio of the external smectitic surfaces to the internal illitic surfaces (e.g., Środoń et al., 1992). In the second case, we assume that over geological time, 100% of $^{40}$Ar is lost from the external part of the edge (red color in Fig. 9) due to the thermal motions and substitution by other atoms/molecules (ions or water), because there is no bonding that can keep Ar in interlayer. From the next-to-external part of the crystallite (yellow color in Fig. 9), the total Ar loss is a sum of the direct effect of Ar recoil, passing between most external potassium ions (3.8%) and the Ar embedding in most external positions after bouncing of the external $K^+$ ions out of the crystallite (4.8%). It was assumed again that all $^{40}$Ar is then lost from the final external positions due to thermal motions.

<figure><img src="image_19.png" /></figure>

Fig. 10. Percentage of $^{40}$Ar loss from an illite particle as a function of the number of illite layers (2–10) for different particle fractions. Left plot represents the total calculated $^{40}$Ar loss while the right one represents the theoretical $^{40}$Ar loss after removing the potassium cations from low-Ar retention sites at crystallite edges. Note that the vertical scales are different between the plots.
<figure><img src="image_20.png" /><figcaption>Fig. 11. Apparent K–Ar age differences for crystallites of different grain size fraction (ab plane) and thickness (number on a bar), computed for the crystallization events of 20, 250, and 500 Ma. Gray bar – apparent K–Ar age before K removal from low-Ar retention sites at crystallite fringes; Black (total) bar – the K–Ar age corrected for the K present at illite outermost fringes.</figcaption></figure>

The calculations in which $ ^{40}Ar $ can escape only through recoil are geologically irrelevant, because they describe only a short period of time after each consecutive recoil. The transformation of $ ^{40}K $ to $ ^{40}Ar $ is a continuous process, and dating of a mineral occurs at a certain moment of time, therefore the effect of recoil is always overlapped by the effect of thermal motions. The calculations are shown thus only for such case, i.e., the percentage of $ ^{40}Ar $ lost due to both the recoil effect and the thermal motions (Fig. 10). The net $ ^{40}Ar $ loss shows a substantial dependency on the size and number of layers: it can vary between 10% for the finest crystallites to close to zero for the thickest and largest (in the *ab* plane) ones. When virtually removing the potassium that occupies the most external positions at illite crystallite edges (red zone in Fig. 9), i.e., the potassium atoms parent for Ar lost due to thermal motions, the total Ar loss becomes significantly corrected and drops by a factor of 2–4, or even more (Fig. 10).

## 5. CONCLUSIONS AND IMPLICATIONS FOR AR-BASED GEOCHRONOLOGY

Present results strongly support and provide a firm atomic scale background to the illite K–Ar and $ ^{40}Ar/^{39}Ar $ dating methodology. They clearly prove that a predominant fraction of $ ^{40}Ar $ remains in its initial interlayer position and that the recoil-induced $ ^{40}Ar $ loss has little contribution. This is experimentally supported by numerous K–Ar dating studies on the finest illite–smectite fractions showing no considerable Ar loss (Clauer et al., 1997; Środon et al., 2002; Osborn et al., 2014) and by $ ^{40}Ar/^{39}Ar $ illite dating model of the $ ^{39}Ar $ recoil random walk (Hall et al., 2000), which has been remarkably successful in explaining stereotypical age spectral shapes.

The “low-Ar retention sites” identified by Aronson and Douthitt (1986) in illite–smectite structure are, in fact, the exchange sites on smectitic surfaces of an illite–smectite crystallite (interlayers and crystallite outer surfaces), i.e., sites where cations are not locked, thus there is no mechanism to keep Ar in place. If such sites are occupied by potassium, this can lead to a significant underestimation of the age (Thompson and Hower, 1973; Środon et al., 2002; Derkowski et al., 2009). Therefore, appropriate cation exchange procedure is always required prior to

K–Ar or Ar–Ar dating. Fringes of illite crystallite, however, also belong to the “low-Ar retention sites”. They can be formed by more than one outermost layer of unit cells (Fig. 9) and may even propagate deep toward the center of a crystallite (McKinley et al., 2004). If the fringe is open wide enough that local $ ^{40}Ar $ is removed due to thermal motions, potassium can also be removed from these sites by either natural or laboratory cation exchange, that is supported by e.g., $ Cs^+ $ fixation in illite (Wampler et al., 2012; Fuller et al., 2014). Removal of potassium from the illite fringes is expected to be more difficult by common cation exchange procedures than from the surface of illite crystallites due to strong electrostatic forces retaining potassium in the site. The cation exchange procedure with low-hydration (high attraction to the $ K^+ $ position) cations that can easily substitute interlayer $ K^+ $, i.e., $ Cs^+ $ or $ Rb^+ $ (Wampler et al., 2012; Fuller et al., 2014) may thus be suggested for the most complete potassium removal from the low-Ar retention sites. Such experiments could verify the results of this study.

The apparent K–Ar age decrease is proportional to the fraction of $ ^{40}Ar $ loss. Therefore, the difference in $ ^{40}Ar $ loss among crystallites of different *ab* and $ c^* $ dimensions can manifest itself as an apparent age difference, simulating the real K–Ar age measurement. The finer crystallites should show lower apparent age than the coarser ones, if the age of crystallization is assumed equal for all crystallites (Fig. 11).

From the model presented in Fig. 10 it is clear that the lack of $ K^+ $ removal from illite fringes (potentially Ar-free) strongly increases the apparent age differences among crystallites. The gray and black bars, showing the K–Ar apparent ages with the illite crystallite fringes occupied by $ K^+ $ and free from $ K^+ $ not linked to the $ ^{40}Ar^* $, respectively, represent the end-members of the real cases (c.f. Derkowski et al., 2009). Some potassium is always removed from the illite crystallite fringes, but there is always some K left, strongly attached, despite the procedure involved (Aronson and Douthitt, 1986; Wampler et al., 2012).

Fine equivalent grain size fractions of <0.02 µm and 0.02–0.05 µm or <0.05 µm, as presented in Figs. 10 and 11, can be separated as concentrates of fundamental illite particles of different sizes from shales and bentonites, in order to determine the age and pathway of diagenesis from the perspective of illite crystal growth mechanism (Clauer et al.,
1997; Šrodoń et al., 2002). However, in the case of K–Ar age increase from the finest to coarser illite subfractions, which is a pattern commonly found in shales and even in bentonites, the detrital contamination has to be considered as equally possible as a certain pattern of illite crystal growth (Šrodoń, 1999; Šrodoń et al., 2002). As demonstrated in this study, the pattern of decreasing K–Ar age in the finer fractions can represent the Ar loss due to natural Ar recoil. While the absolute difference among three cation-exchanged, fringe-K-removed crystallites is acceptably low for the crystallization model of 20 Ma (0.4 Ma), it reaches up to 9 Ma for an event 0.5 Ga old. In such a case, the coarse crystallite fraction represents the lowest Ar loss thus the age closest to the true age. Additional criteria will have to be applied to distinguish the Ar recoil loss effects from the detrital contamination effects.

The ⁴⁰Ar loss found in this study is an order of magnitude lower than that suggested by Lerman et al. (2007) for fine illite particles. Because the activation energy of Ar release from the illite interlayer is high, from 47 to 66 kcal/mol (Hassanipak and Wampler, 1996), similar to muscovite 63 ± 7 kcal/mol (Harrison et al., 2009), there is very little probability that tens of % of ⁴⁰Ar could have been lost from fine 20 Å particles due to diffusion. Argon recoil, however, is a high-energy reaction that can easily drive the ⁴⁰Ar atoms into the TOT layer and, eventually, out of the illite structure. Besides the obvious necessity of removing exchangeable K from smectitic surfaces, a correct technique to remove potassium from low-Ar retention sites at illite crystallite fringes seems to be the next crucial step toward correct measurement and interpretation of the K–Ar age of illite–smectite.

## ACKNOWLEDGEMENTS

This work was inspired by the work of Goncharov et al. (2007) as well as by the discussions with K. Balogh to whom the authors are very grateful.

This work was supported by the collaborative project ATLAB (FP7-REGPOT-2011-1 No. 285989) funded by the European Union. A.G.K. acknowledges the support of the industrial chair “Storage and Disposal of Radioactive Waste” at the Ecole des Mines de Nantes, funded by ANDRA, Areva, and EDF.

Computer time allocations made available within the Distributed European Computing Initiative (project DEC07_NUWCLAY and DEC11 COMPCLAY by the PRACE-2IP receiving funding from the European Community's FP7/2007–2013 under grant agreement R1-283493) and within PLGRID infrastructure are also gratefully acknowledged. M.S. also acknowledges NCN (grant 2012/05/B/ST10/01948), which provided IGS PAS with high-performance computing server.

## APPENDIX A. SUPPLEMENTARY DATA

Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.gca.2015.03.005.

## REFERENCES

* Aronson J. L. and Douthitt C. B. (1986) K/Ar systematics of an acid-treated illite/smectite: implications for evaluating age and crystal structure. *Clays Clay Miner.* **34**, 473–482.

* Bé M. M., Chisté V., Dulieu C., Mougeot X., Browne E., Chechev V., Kuzmenko N., Kondev F., Luca A., Galán M., Nichols A. L., Arinc A. and Huang X. (2010) *Table of Radionuclides (Vol. 5-A = 22 to 244)*. Monographie BIPM-5. Bureau International des Poids et Mesures, Sèvres, France, p. 256.

* Čermák J. and Král L. (2010) Low-temperature diffusion of carbon in iron measured by SIMS technique. In *METAL 2010 Conference Proceedings*, 18–20.

* Churakov S. V. (2007) Structure and dynamics of the water films confined between edges of pyrophyllite: a first principle study. *Geochim. Cosmochim. Acta* **71**, 1130–1144.

* Clauer N., Šrodoń J., Franců J. and Šucha V. (1997) K–Ar dating of illite fundamental particles separated from illite–smectite. *Clay Miner.* **32**, 181–196.

* Clauer N., Zwingmann H., Liewig N. and Wendling R. (2012) Comparative ⁴⁰Ar/³⁹Ar and K–Ar dating of illite-type clay minerals: a tentative explanation for age identities and differences. *Earth Sci. Rev.* **115**, 76–96.

* Cottrell T. L. (1958) *The Strengths of Chemical Bonds*, second ed. Butterworth, London, p. 317.

* Cygan R. T., Liang J. J. and Kalinichev A. G. (2004) Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. *J. Phys. Chem. B* **108**, 1255–1266.

* Dahl P. S. (1996) The crystal-chemical basis for Ar retention in micas: inferences from interlayer partitioning and implications for geochronology. *Contrib. Miner. Petrol.* **123**, 22–39.

* Darwent B. deB. (1970) Bond dissociation energies in simple molecules. *National Standard Reference Data Series, National Bureau of Standards*, no. 31, Washington, p. 48.

* Dalrymple G. B. and Lanphere M. A. (1971) ⁴⁰Ar/³⁹Ar technique of K–Ar dating: a comparison with the conventional technique. *Earth Planet. Sci. Lett.* **12**, 300–308.

* Derkowski A., Šrodoń J., Franus W., Uhlik P., Banaś M., Zieliński G., Čaplovičová M. and Franus M. (2009) Partial dissolution of glauconitic samples: implications for the methodology of K–Ar and Rb–Sr dating. *Clays Clay Miner.* **57**, 531–554.

* Derkowski A., Drits V. A. and McCarty D. K. (2012) Nature of rehydroxylation in dioctahedral 2:1 layer clay minerals. *Am. Mineral.* **97**, 610–629.

* Dong H., Hall C. M., Peacor D. R. and Halliday A. N. (1995) Mechanisms of argon retention in clays revealed by laser ⁴⁰Ar–³⁹Ar dating. *Science* **267**, 355–359.

* Drits V. A., Besson G. and Muller F. (1995) An improved model for structural transformation of heat-treated aluminous dioctahedral 2:1 layer silicates. *Clays Clay Miner.* **43**, 718–731.

* Dudek T., Šrodoń J., Eberl D. D., Elsass F. and Uhlik P. (2002) Thickness distribution of illite crystals in shales. Part I: XRD vs HRTEM measurements. *Clays Clay Miner.* **50**, 562–577.

* Elliott W. C. and Haynes J. T. (2002) The chemical character of fluids forming illite in the Southern Appalachian Basin (in: Reynolds Commemorative Volume). *Am. Mineral.* **87**, 1519–1528.

* Fuller A. J., Shaw S., Peacock C. L., Trivedi D., Small J. S., Abrahamsen L. G. and Burke I. T. (2014) Ionic strength and pH dependent multi-site sorption of Cs onto a micaceous aquifer sediment. *Appl. Geochem.* **40**, 32–42.

* Goncharov G. N., Zaitseva T. S. and Aref'ev I. A. (2007) Crystallochemical peculiarities and the retention of K-Ar ages of Upper Riphean globular phyllosilicates (the Uk Formation, the South Urals). In *The Youth Conference “Geology is Our Future”. Geological Faculty of St. Petersburg State University, St. Petersburg, Book of Abstracts*, 18–19.

* Hall C. M., Kesler S. E., Simon G. and Fortuna J. (2000) Overlapping cretaceous and Eocene alteration, Twin Creeks Carlin-type deposit, Nevada. *Econ. Geol.* **95**, 1739–1752.
Harrison T. M., Celerier J., Aikman A. B., Herman J. and Heizler
M. T. (2009) Diffusion of <sup>40</sup>Ar in muscovite. *Geochim. Cosmochim. Acta* **73**, 1039–1051.

Hassanipak A. A. and Wampler J. M. (1996) Radiogenic argon released by stepwise heating of glauconite and illite: the influence of composition and particle size. *Clays Clay Miner.* **44**, 717–726.

Heinz H., Koerner H., Vaia R. A., Anderson K. L. and Farmer B.
L. (2005) Force field for phyllosilicates and dynamics of
octadecylammonium chains grafted to montmorillonite.
*Chem. Mater.* **17**, 5658–5669.

Hetherington C. J. and Villa I. M. (2007) Barium silicates of the Berisal Complex, Switzerland: a study in geochronology and rare-gas release systematics. *Geochim. Cosmochim. Acta* **71**, 3336–3347.

Lerman A. and Clauer N. (2005) Losses of radiogenic <sup>40</sup>Ar in sediment fine-clay size fractions of sediments. *Clays Clay Miner.* **53**, 233–248.

Lerman A., Ray B. M. and Clauer N. (2007) Radioactive production and diffusional loss of radiogenic <sup>40</sup>Ar in clays in relation to its flux to the atmosphere. *Chem. Geol.* **243**, 205–224.

Liu X., Lu X., Meijer E. J., Wang R. and Zhou H. (2012) Atomic-scale structures of interfaces between phyllosilicate edges and water. *Geochim. Cosmochim. Acta* **81**, 56–68.

McKinley J. P., Zachara J. M., Heald S. M., Dohnalkova A.,
Newville M. G. and Sutton S. R. (2004) Microscale distribution
of cesium sorbed to biotite and muscovite. *Environ. Sci.
Technol.* **38**, 1017–1023.

Molina-Montes E., Donadio D., Hernandez-Laguna A. and Sainz-Díaz C. I. (2008) DFT research of the dehydroxylation reaction of pyrophyllite 2. Characterisation of reactants, intermediates, and transition states along the reaction path. *J. Phys. Chem. B* **112**, 6373–6383.

Molina-Montes E., Donadio D., Hernandez-Laguna A. and Sainz-díaz I. (2010) Exploring the rehydroxylation reaction of pyrophyllite by ab initio molecular dynamics. *J. Phys. Chem. B* **114**, 7593–7601.

Nadeau P. H., Wilson M. J., McHardy W. J. and Tait J. M. (1984) Interstratified clay as fundamental particles. *Science* **225**, 923–925.

Onstott T. C., Miller M. L., Ewing R. C., Arnold G. W. and Walsh
D. S. (1995) Recoil refinements: implications for the
<sup>40</sup>Ar/<sup>39</sup>Ar dating technique. *Geochim. Cosmochim. Acta* **59**,
1821–1834.

Onstott T. C., Mueller C., Vrolijk P. J. and Pevear D. R. (1997)
Laser <sup>40</sup>Ar/<sup>39</sup>Ar microprobe analyses of fine-grained illite.
*Geochim. Cosmochim. Acta* **61**, 3851–3861.

Osborn S. G., Duffield L. T., Elliott W. C., Wampler J. M.,
Elmore R. D. and Engel M. H. (2014) The timing of
diagenesis and thermal maturation of the Cretaceous Marias
River Shale, Disturbed Belt, Montana. *Clays Clay Miner.*
*62*(2), 112–125.

Palandri J. L. and Kharaka Y. (2004) A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling. In *U.S. Geological Survey Open File Report 2004-1068*. U.S. Department of the Interior, p. 64.

Pitman M. C. and van Duin A. C. T. (2012) Dynamics of confined reactive water in smectite clay-zeolite composites. *J. Am. Chem. Soc.* **134**, 3042–3053.

Plimpton S. (1995) Fast parallel algorithms for short-range molecular dynamics. *J. Comput. Phys.* **117**, 1–19.

Plummer L. N., Wigley T. M. L. and Parkhurst D. L. (1978) The kinetics of calcite dissolution in CO₂-water systems at 5 to 60 °C and 0.0 to 1.0 atm CO₂. *Am. J. Sci.* **278**, 179–216.

Sainz-Díaz I., Escamilla-Roa E. and Hernandez-Laguna A. (2004) Pyrophyllite dehydroxylation process by First Principles calculation. *Am. Mineral.* **80**, 1092–1100.

Sletten V. W. and Onstott T. C. (1998) The thermal instability of muscovite during in vacuo heating: implications for Ar loss systematics in the laboratory. *Geochim. Cosmochim. Acta* **62**, 123–141.

Smith P. E., Evensen N. M. and York D. (1993) First successful <sup>40</sup>Ar-<sup>39</sup>Ar dating of glauconies: argon recoil in single grains of cryptocrystalline material. *Geology* **21**, 41–44.

Šrodoń J. (1999) Extracting K–Ar ages from shales: a theoretical test. *Clay Miner.* **33**, 375–378.

Šrodoń J., Elsass F., McHardy W. J. and Morgan D. J. (1992) Chemistry of illite/smectite inferred from TEM measurements of fundamental particles. *Clay Miner.* **27**, 137–158.

Šrodoń J., Clauer N. and Eberl D. D. (2002) Interpretation of K–Ar dates of illitic clays from sedimentary rocks aided by modelling. *Am. Mineral.* **87**, 1528–1535.

Šrodoń J., Zeelmaekers E. and Derkowski A. (2009) Charge of component layers of illite-smectite in bentonites and the nature of end-member illite. *Clays Clay Miner.* **57**, 650–672.

Tazi S., Rotenberg B., Salanne M., Sprik M. and Sulpizi M. (2012) Absolute acidity of clay edge sites from ab-initio simulations. *Geochim. Cosmochim. Acta* **94**, 1–11.

Thompson G. and Hower J. (1973) An explanation for the low radiometric ages from glauconite. *Geochim. Cosmochim. Acta* **37**, 1473–1491.

van der Pluijm B. A., Hall C. M., Vrolijk P., Pevear D. R. and Covey M. (2001) The dating of shallow faults in the Earth's crust. *Nature* **412**, 172–174.

van Duin A. C. T., Dasgupta S., Lorant F. and Goddard W. A. (2001) ReaxFF: a reactive force field for hydrocarbons. *J. Phys. Chem. A* **105**, 9396–9409.

Vidal O., Dubacq B. and Lanari P. (2010) Comment on “the role of H₃O⁺ in the crystal structure of illite” by F. Nieto, M. Melini, and I. Abad. *Clays Clay Miner.* **58**, 717–720.

Wampler J. M., Krogstad E. J., Elliott W. C., Kahn B. and Kaplan D. I. (2012) Long-term selective retention of natural Cs and Rb by highly weathered coastal plain soils. *Environ. Sci. Technol.* **46**, 3837–3843.