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

The radioactive decay of 40 K to 40 Ar is the basis of isotope age determination of micaceous clay minerals formed during diagenesis. The diﬀerence in K–Ar ages between ﬁne and coarse grained illite particles has been interpreted using detrital-au-thigenic components system, its crystallization history or post-crystallization diﬀusion. Yet another mechanism should also be considered: natural 40 Ar recoil. Whether this recoil mechanism can result in a signiﬁcant enough loss of 40 Ar to provide observable decrease of K–Ar age of the ﬁnest 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 40 Ar atom, but for varying recoil angles that cover the entire range of their possible values. The results show that 40 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 ﬁnal positions of the 40 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 simulation angles were found for which the 40 Ar atom after recoil passes completely through the 2:1 layer. The energy barrier for 40 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 40 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 40 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 eﬀect of 40 Ar recoil on the K–Ar dating of illite show that some of 40 Ar is lost and the loss is substantially dependent on the crystallite dimensions. The 40 Ar loss can vary from 10% for the ﬁnest crystallites (two 2:1 layers thickness and <0.02 l 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 40 Ar loss, the ﬁner 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 diﬀerences among crystallites of diﬀerent size.


INTRODUCTION
The radioactive decay of 40 Kt o 40 Ar is the basis of the isotope age determination of micaceous clay minerals: illite, glauconite, aluminoceladonite, sericite (fine variety of muscovite), and their mixed-layered interstratifications with smectite.These dating methods are based on the assumption 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 typical for a 2:1 layer silicate: an octahedral Al-dominated sheet (O) sandwiched between two tetrahedral sheets (T), where SiO 4 4À tetrahedra are partially substituted by AlO 4 5À and the excess of negative charge is compensated by the presence of potassium cations in the interlayer.Radiogenic 40 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 40 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 40 Ar/ 39 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., S ´rodon ´, 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 lm in bentonites), their crystallization history and illite-smectite crystallite growth may play a major role (Clauer et al., 1997;S ´rodon ´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 isotope 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 40 Ar loss in the diagenetic temperature zone over the geological time span?
While the K-Ar methodology measures K and radiogenic 40 Ar concentrations, 40 Ar-39 Ar method involves prior neutron irradiation in order to transform 39 K into 39 Ar for 40 Ar/ 39 Ar ratio measurement (Dalrymple and Lanphere, 1971).During the irradiation-induced 39 Kt o 39 Ar transformation, high recoil energy of 39 Ar, $177 keV (Onstott et al., 1995), allows 39 Ar penetration through the clay structure and its escape from the fine crystallites (Smith et al., 1993;Dong et al., 1995).Due to the high 39 Ar recoil energy and the resulting substantial 39 Ar loss, the application of the 40 Ar/ 39 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 39 Ar lost from the illite structure during 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 39 Ar recoil not only stimulates the 39 Ar nuclei diffusion along the interlayer, but also its penetration through the TOT layers.
If the high-energy irradiation-induced 39 Ar recoil leads to the 39 Ar loss from fine illite particles, perhaps the natural 40 Ar recoil can also result in some 40 Ar loss, lower than in the case of 39 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 40 Ar loss from fine illite crystallites, which serves as reference for the irradiation-induced 39 Ar loss.Despite numerous successful applications of the model in 40 Ar/ 39 Ar dating of illite, the assumption about 40 Ar has not been independently validated and the 40 Ar loss, if occurring, may be covered by a strong effect of 39 Ar loss, thus not detectable by a model calculations.
The illite-smectite crystallites commonly have thicknesses of their fundamental illite particles from two to several 2:1 layers (Nadeau et al., 1984;S ´rodon ´et al., 1992;Dudek et al., 2002).If the 40 Ar recoil energy can make some 40 Ar atoms penetrate through a single 2:1 layer, the broadly accepted assumption that all radiogenic 40 Ar remains mechanically trapped in the interlayer may be incorrect.It is also not clear whether 40 Ar after recoil can be incorporated somewhere else into the clay structure or its recoil energy is not sufficient to break the chemical bonds of the clay layer, so 40 Ar can remain only in the crystallographic positions of potassium.Similar questions were hypothesized earlier by Goncharov et al. (2007).Hetherington and Villa (2007) argued that a minor portion of 40 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 geological age estimates.

Molecular dynamics simulations
In order to calculate the effect of argon recoil, the kinetic energy and velocity of 40 Ar after recoil were calculated on the basis of the energy of emitted gamma quantum after the 40 K-40 Ar electron capture transition: 1460.822keV (Be ´et al., 2010).This transition produces a photon with a mass of 1.56825 Â 10 À3 amu.Employing the principle of momentum conservation we can calculate the 40 Ar atom velocity as 0.11765 A ˚/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 h is defined as an angle between the recoil vector and the direction perpendicular to the clay layer ("vertical").The angle u is defined as an angle between 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 (h % 90°).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 h <80°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 h 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 semiempirical 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 h 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 2 Si 4 O 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 h was set in the range from 0°t o 80°with a step of 10°, while the horizontal angle u was set in the range from 0°to 360°with a step of 20°( 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°A ˚/fs.It is important to note, that all simulations were intrinsically nonequilibrium, 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 h higher than 60°, 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 40 Ar atom at h 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 (S ´rodon ´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 IV substituted sites.The recoil effect in different directions (angles u from 0°to 60°, varying by a step of 10°) was studied assuming that the interlayer positions surrounding the tested 40 Ar atom were either all occupied by potassium cations or one position was left vacant (Fig. 2).However, H 2 O molecules may occur in the K-vacant positions (Vidal et al., 2010).In order to further test the fate of the 40 Ar atom after recoil at the edges of crystallites, calculations assuming 40 Ar-neighboring vacancies (dashed black lines in Fig. 2) were also performed to mimic the crystallite fringe edge with angles u varying from À60°to 120°, with a step of 10°and h 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 À6 .

Possible recoil reactions
The recoil energy of 40 Ar is 28.66 eV, which is equal to 660.9 kcal for 1 mol of 40 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 40 Ar in the clay structure or the final residual deformation of the clay structure.

Recoil effects from the perspective of final position of 40 Ar in the structure (REAXFF)
Four possible final positions of 40 Ar can be identified at the end of the REAXFF simulation when the initial structural perturbation is dissipated (Fig. 3).For h >40°, the most common effect observed is that the 40 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 h <40°, if the direction of the recoil vector is pointing to a Si atom of one of the SiO 4 tetrahedra.No bond breaking was then observed (Fig. 3a, b; Supplementary video 1).
The most common effect observed for h <40°is the penetration of 40 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 40 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, 40 Ar Electronic Annex: Video 1 can occupy the empty cavity after the OH group dissociated from the tetrahedral sheet (Fig. 3c, d; Supplementary video 2).
More precisely, 40 Ar finds its new position close to the position of the displaced OH group.Depending on the extent of the OH group displacement, 40 Ar is finally located closer to the surface or deeper within the layer and closer to the octahedral position.
Only for three combinations of h and u angles we observe the penetration of the octahedral sheet by the 40 Ar atom upon recoil.In one case, 40 Ar is recoiled directly into the vacancy in the octahedral sheet, avoiding any energy transfer to the OH group.In two other cases, 40 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).
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 40 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).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 40 Ar atom after recoil to pass completely through the TOT layer to the interlayer on the opposite side.

Layer structure deformation due to 40 Ar recoil (REAXFF)
Although the major layer structure deformation must result from a b-decay ( 40 K ! 40Ca) that occurs $89.3% of the time, with energy similar to the electron capture (1.33 MeV), understanding the behavior of 40 Ar after recoil and its recoil pathway requires determining the 2:1 layer structure deformations related exclusively to 40 Ar recoil.The possible effects of 40 Ar recoil on the deformation of the clay structure vary depending on the recoil angle.The most common reaction at low vertical angles (h <40°)i s associated with kinetic energy transfer to the closest OH group and a subsequent recoil of this group.Depending on the 40 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.After the water molecule is formed, 40 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 40 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 4 (OH) 2 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).
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 (h P 50°), the recoiled 40 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 40 Ar bouncing back to the interlayer: accompanied by simultaneous bond breaking and reformation in the octahedral sheet (Fig. 5; Supplementary video 8).
Bouncing of the incoming 40 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 h >60°.If the recoil vector is pointing to a Si atom of a SiO 4 tetrahedron, the recoil energy of the 40 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 h % 30°where 40 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 h irreversible deformations were observed within the simulation time of our studies.

Spatial distribution of the effects of 40 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 function of angles h and u (in a hemispherical projection) is necessary in order to assess the 40 Ar behavior over geological 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 trajectory with the TOT layer surface (a plane formed by the centers 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 h, but some dependence on the angle u is also observed.For h 6 40°, 40 Ar is most often positioned in the tetrahedral sheet.Combining this information with the data of Fig. 6b, which depicts the distortions 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 usually 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 h >40°a nd h 6 40°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 h) 40 Ar initially hits the basal oxygen surface of the tetrahedral sheet and the recoil usually leads to the deformation of the tetrahedral sheet.For h % 70°, 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 fractions of the sphere surface corresponding to each combination of the two recoil angles, h and u.

Maximum recoil distance of 40 Ar along the interlayer (CLAYFF; h close to 90°)
As shown in the previous sections, if the angle h approaches 90°, the effect of 40 Ar recoil on the TOT structure 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 position, independent of the recoil angle u.
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 6 u 6 60°).Only for u % 0°(or, by analogy, for u % 120°) 40 Ar bounces back to the initial position.For any vertical angle h <60°at u =30°, 40 Ar bounces back and returns to its initial position.
For a vacancy located in the position of a second neighbor (position 2 on Fig. 2), the range of angles u leading to an eventual jump of the recoiled 40 Ar to the vacant position is approximately between 20°and 50°.However, if h <80°, 40 Ar bounces back and returns to its initial position regardless the u angle.
If a vacancy is located in the farther neighboring position 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 < u <40°.For a vacancy located in the position 4 of Fig. 2, the 40 Ar always bounces back to its initial position independently of the initial recoil angle u.
The results of simulations performed for linear series of neighboring vacancies as a model mimicking the crystallite Electronic Annex: Video 8 fringe edge show substantial dependence on angles h and u.For h <80°, 40 Ar located behind the row of potassium cations bounces back and returns to its initial position regardless of the value of u.For h equal to 80°or 90°the fate of 40 Ar is very similar and with important dependence on the recoil angle u (Fig. 7).In some cases, if 40 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 u for which 40 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 substantial fraction of simulations in which 40 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 Ar switched its position with a neighboring potassium ion (position 1b in Fig. 7).The "crystallite edge"  simulation results extrapolated to illite/smectite crystallites indicate that if 40 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.

Potential energy of 40 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: -40 Ar passing through the hexagonal (ditrigonal) cavity, in the direction from the tetrahedral sheet to the interlayer (Fig. 8a). -4 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 40 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 h = 0, i.e., the recoil vector pointing perpendicularly to the TOT layer.In the second case the structure corresponded to 40 Ar located in the octahedral sheet with the initial recoil angles h and u set to 20°.
The average potential energy of the structure with 40 Ar located in the tetrahedral sheet is higher by about 26 kcal/mol compared to the case when 40 Ar is eventually located out of the TOT structure (Fig. 8a).The energy barrier for 40 Ar passing through the hexagonal cavity is %17 kcal/mol.For comparison, the recoil energy of 40 Ar is equal to 660.9 kcal/mol, which indicates that 40 Ar can pass easily through the hexagonal cavity at the angles close to h =0.
In the case of the structure with 40 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 40 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 40 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.The extent of structural deformation clearly depends on the recoil angle.The deformation for low h 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 dioctahedral 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 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.

Probability of retaining the 40 Ar atom in the layer structure over geological time periods
Taking into account the energy barrier of 40 Ar passing through the hexagonal cavity of the 2:1 clay structure, which is ca.17 kcal/mol, and DE of this reaction (À26 kcal/mol) there is hardly any possibility of 40 Ar staying within the tetrahedral sheet of the TOT layer structure during a time span of geological scale.DE of the reaction was calculated for the 40 Ar transfer from TOT layer to vacuum; it will be lower if 40 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 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 (C ˇerma ´k and Kra ´l, 2010).The calculated energy barrier thus allows predicting that all 40 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 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 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 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 Ar ends up within the octahedral sheet.Much higher energy barrier of 55 kcal/mol indicates that after entering the octahedral sheet, 40 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 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 (S ´rodon ´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 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 are substantially different in their structure and composition from the interior due to the formation of OH groups on the crystallite edges and the weakening 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 happen 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., S ´rodon ´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.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 Kt o 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).

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 subfractions showing no considerable Ar loss (Clauer et al., 1997;S ´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;S ´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 lm and 0.02-0.05lm or <0.05 lm, 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.,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.1997; S ´rodon ´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 (S ´rodon ´, 1999; S ´rodon ´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 40 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 40 Ar could have been lost from fine 20 A ˚particles due to diffusion.Argon recoil, however, is a high-energy reaction that can easily drive the 40 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.

Fig. 1 .
Fig. 1.Definition of the recoil angles h and u. 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.)

Fig. 2 .
Fig. 2. Definition of the recoil angle u 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.
(2) The recoiled OH group created new bonding with another OH group, resulting in the formation Electronic Annex: Video 2 Electronic Annex: Video 3 Electronic Annex: Video 4 of a water molecule and a residual oxygen (Fig. 4a, b; Supplementary video 5).The H 2 O molecule can then behave in three different ways: -leave the TOT structure (as in Fig. 4a, b); -remain as a weakly-bonded H 2 O molecule within the TOT layer; -form a new bonding with two Al atoms of the octahedral sheet (oxygen of this H 2 O molecule is at the position of the former oxygen of the OH group, Fig. 4c, d; Supplementary video 6).

( 1 )
If the recoil vector is pointing to a Si atom of the tetrahedron, the 40 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 40 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 40 Ar recoil leads to a substantial deformation of the tetrahedral sheet, followed by the formation of new OH groups attached to Si tetrahedra, and Electronic Annex: Video 5 Electronic Annex: Video 6 Electronic Annex: Video 7

Fig. 3 .
Fig. 3. Recoil effects resulting in different final position of the 40 Ar atom in the structure (final simulation frames): (a, b) the largest deformation amplitude when 40 Ar atom bounces from the TOT layer without deforming the structure.(c, d) 40 Ar atom is embedded in the tetrahedral sheet.(e, f) 40 Ar atom is embedded in the octahedral sheet.(g, h) 40 Ar atom is embedded in the tetrahedral sheet on the opposite side of the TOT layer.

Fig. 4 .
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 Ar.A final simulation frames represent: (a, b) a new H 2 O molecule is formed and ejected from the TOT sheet, while 40 Ar is bounced back from the TOT surface; (c, d) a new H 2 O molecule is formed in the position of one of the initial OH groups, H 2 O oxygen forms two penta-coordinated Si sites, while 40 Ar atom is embedded in the tetrahedral sheet; (e, f) the OH group forms two penta-coordinated Si sites, while 40 Ar atom is bounced back from the oxygen surface.

Fig. 5 .
Fig. 5. Recoil effect under high angle h leading to the deformation of the tetrahedral sheet and, subsequently, of the octahedral sheet.Ar (not shown) has bounced back to the interlayer.

Fig. 6 .
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 40 Ar atom in the structure: orange -40 Ar in the tetrahedral sheet on the same side of the TOT layer, yellow -40 Ar in the octahedral sheet, green -40 Ar bounced back from the TOT surface, brown -40 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 2 O molecule stays in the TOT layer, blue -newly formed H 2 O molecule leaves the TOT structure, orangedistortion 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 h angle.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Deformation of the clay structure associated with 40 Ar recoil in natural conditions

Fig. 8 .
Fig. 8. Moving averages of the potential energy of (a) 40 Ar passing through the hexagonal cavity from an initial location in the tetrahedral sheet after the recoil perpendicular to the surface; (b) 40 Ar passing the TOT layer from initial octahedral position.

Fig. 10 .
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.

Fig. 9 .
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.)