Modeling the Mechanisms of Clay Damage by Molecular Dynamic Simulation

We carried out MD simulations to review the mechanisms of clay damage under hydration, the manners clays interact with water and cations, and the structures they form in between aluminosilicate sheets and offered an alternative explanation as to how clay structures become instable after certain water concentration. Water molecules form one, two, and sometimes three layers parallel to the clay surfaces. The cations already present in montmorillonites (the most representative of the smectites) or coming with low salinity water will either be adsorbed by the clay surfaces or become fully hydrated, detached from the surfaces, and located in the midst of the interlayer space between clay surfaces, while water molecules that do not hydrate the cations invade the clay structures, making them unstable. The driving force of this phenomenon happens to be the hydration energy of the cations in question.


Introduction
Clay minerals are encountered in all types of petroleum bearing formations.Under interaction with water, clays undergo two phenomena that have a major and direct impact in permeability impairment of formations: (1) swelling effects and (2) instability of the structure.The former has been widely studied and rather well characterized theoretically [1], experimentally [2][3][4], and using Monte Carlo (MC) simulations [5][6][7][8][9], molecular dynamics (MD) simulations [10,11], and ab initio molecular dynamics simulations [12].Clay instability is still to be explained properly.
Among the group of clay minerals, smectites comprise the subset of clays that readily swell in the presence of water.Smectites are a natural product of the weathering and decomposition of igneous rocks [13] or other rocks, including shales, consisting of negatively charged two-dimensional aluminosilicate layers.Each layer is composed of tetrahedral silicate sheets surrounding octahedral aluminum sheets.Isomorphic substitutions can occur in either the tetrahedral (Si 4+ with Al 3+ ) or octahedral (Al 3+ with Mg 2+ or Fe 2+ ) sheets, which yields a net negative charge that is balanced by cations located in the interlayer region.
Barshad [2] studied the factors affecting the interlayer expansion of vermiculite and montmorillonite with organic substances and determined that the extent of interlayer expansion was affected primarily by the size, charge, and total amount of the cations and by the magnitude of the dipole moment and the dielectric constant of the immersion liquid.Norrish [1] carried out experiments focused on the swelling of montmorillonite and confirmed that this effect exhibits two regimes: crystalline and osmotic swelling.Crystalline swelling can occur in all types of clay minerals [14], being the process in which adsorbed water increases to approximately 0.5 g H 2 O/g clay while the interlayer spacing increases from approximately 9.5 Å (for dry material) to roughly 20 Å.Several other studies have shown that the swelling process occurs by increasing the water content through the formation of one, two, and perhaps three layer hydrates; Sun et al. [15] summarized several works in which the majority of experimental studies reported the formation of one-layer (1W) and two-layer (2W) hydrates, whereas the formation of the three-layer (3W) hydrate was observed only in a few studies (see, e.g., [5,6,10,16]).Osmotic swelling occurs when montmorillonite is placed in contact with water and takes up 10 g H 2 O/g clay, increasing its volume by about twenty times.
In this paper we present the results from MD simulations for the basal spacing of Na-, K-, and Ca-montmorillonites and the density profiles and snapshots of cations and oxygenwater atoms in the interlayer region.The results for the basal spacing indicate that the K-montmorillonite expands to a higher extent than both Na-and Ca-montmorillonite when water content increases; the density profiles provide insights into how the hydration of cations takes place while the interlayer space increases and how the structure itself becomes unstable.Although Na + and K + are both monovalent cations, and Na + and Ca 2+ have similar ionic radii, their behavior under hydration is different [17].

Conventional Mechanisms of Clay Damage
It has been widely accepted that Na-montmorillonite swells more than Ca-montmorillonite because Ca 2+ cations are strongly adsorbed (by the clay surfaces) compared to Na + cations [18].Accordingly, under contact with water, Camontmorillonite platelets remain practically intact and close to each other, whereas the Na-montmorillonite aggregates readily swell and the platelets separate widely.As a result, water can easily invade the gaps between the platelets and form thicker water envelopes around the Na-montmorillonite platelets than the Ca-montmorillonite platelets.However, recent studies using MD simulations show that Na-montmorillonite and Ca-montmorillonite basal spacing are fairly similar except in that window in which the 2W layer is formed in Ca-montmorillonite but not in Na-montmorillonite, where Ca-montmorillonite exhibits stronger swelling than Na-montmorillonite (see [15]).Moreover, Ca 2+ cations exhibit larger hydration energies relative to Na + cations, which leads to higher water coordination numbers and more pronounced association of water molecules with Ca 2+ cations.These new results using MD simulations demonstrate that the underlying mechanisms of clay swelling are still open to discussion.
It has also been accepted that clay damage can be prevented by maintaining high concentrations of K + cations in aqueous solutions.The proposed explanation is that, due to the small size, the K + cations can readily penetrate the interlayers of the clay and hold the clay platelets together [18,19].To further support this statement, Reed [20] conducted laboratory core tests by flowing deionized water, 3% NaCl brine, and 3% CaCl 2 brine through cores extracted from micaceous sand formations to determine permeability reduction, hypothesizing that mica alteration is a result of the exchange of K + cations with cations of larger sizes, such as Na + , Li + , Ca 2+ , and Mg 2+ , as depicted in Figure 1.Mica alteration generates fines that later deposit in porous rocks.When clays are exposed to brines containing either no or small amounts of K + cations or larger cations, K + cations diffuse out of the clay platelets according to Fick's law due to the concentration gradient between clay and brine, while larger cations diffuse into clays.Since larger cations cannot fit into the interlayer region, the edges of the friable mica flakes break off into small pieces.It is now known that the ionic radii (for a coordination of VI) of the most common ).Reed assumed that the ionic radius of Na + was larger than that of K + .cations in clays are (in pm) Mg 2+ (72) < Li + (76) < Ca 2+ (100) < Na + (102) < K + (138) < Rb + (152) < Cs + (167) [21], so that there must be an alternative mechanism that explains the generation of fines.

Models and Methods
We considered a model of Wyoming-type montmorillonite with octahedral and tetrahedral substitutions, as schematically shown in Figure 2. The substitutions throughout each layer are subject to Loewenstein's substitution rule that the substitution sites are connected by neither one hydroxyl group nor an oxygen atom [22].Al 3+ atoms in the octahedral layers are replaced with Mg 2+ , and Si 4+ atoms in the tetrahedral layers are replaced with Al 3+ .The position of the tetrahedral substitution is not linked to any of the octahedral substitutions by an oxygen atom.Due to the isomorphic substitutions, no structural data with explicit location of substitutions is available in X-ray crystallographic databases.It is customary to build up a model cell for montmorillonite starting with the unit cell of pyrophyllite, which has identical aluminosilicate layers to montmorillonite but exhibits no substitutions (see, e.g., Mignon et al. [12]).A different approach is to start with the atomic positions determined experimentally by Tsipursky and Drits [23] for a smectite sample (see, e.g., Minisini and Tsobnang [24]).We followed the first approach.
The simulation supercells used have 64 unit cells (8 × 4 × 2), with a composition of M 3/ + (Si 31 Al) (Al 14 Mg 2 ) O 80 (OH) 16 ⋅nH 2 O, where M represents the cation, x is the charge of the cation, and  varies from 0 to 10.5 water molecules/unit cell.Figure 2 shows a schematic representation of the montmorillonite layers and interlayer species (H 2 O and Na + ).The unit cell of pyrophyllite is triclinic (belongs to space group C1) with parameters  = 5.160 Å,  = 8.966 Å,  = 9.347 Å,  = 91.18∘ ,  = 100.46∘ , and  = 89.64∘ , as reported by Lee and Guggenheim [25].Cations were placed randomly in the region between the clay layers.Then, 0.5 water molecules per unit cell were introduced into the simulation cell for each simulation, using Packmol [26], up to a maximum of 10.5 water molecules per unit cell.This paper uses a supercell of 64 unit cells to minimize any potential side effects.
MD simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [27].Clayff force field was used to describe the interactions between atoms [10,28].It is based on the single point charge (SPC) water model of Berendsen et al. [29] to represent water, hydroxyl, and oxygen-oxygen interactions.The SPC model has partial charges centered directly on each of three atoms, and the short-range interactions are represented by a Lennard-Jones (LJ) 12-6 term.Bond stretch and bond angle terms are introduced into the SPC model using the expressions determined by Teleman et al. [30] to ensure full flexibility for the water and hydroxide components.It is worth noting that the clayff force field was not fitted to describe the hydrogen bonding network explicitly but include the contributions in the force field parameters.The Lorentz-Berthelot mixing rule was used to obtain the Lennard-Jones parameters for interactions between unlike atoms [31,32].The simulations were performed under periodic boundary conditions, with the long-range electrostatic term treated by the standard Ewald method [31].The cutoff distance for the nonbonded van der Waals interactions was 15 Å and the electrostatic interactions were calculated by the Ewald summation method.All simulations were carried out using an isobaric-isothermal (NPT) ensemble at  = 300 K and  = 1 atm.Pressure was controlled by the Parrinello-Rahman barostat [33] while temperature was controlled by the Nosé-Hoover thermostat [34][35][36].A 10 ns MD run was performed to relax the systems studied, namely, the no water layer (0W), one water layer (1W), and two water layers (2W).Another 1 ns production MD run was then carried out for data analysis.

Results and Discussion
4.1.Swelling Curves.The basal spacing for Na-montmorillonite after hydration obtained from our simulations is presented in Figure 3(a) along with the experimental data obtained by Fu et al. [4].The error bars represent two standard deviations.The simulation results are in very good agreement with experiment.
The results exhibit the expansion of the clay through two well-defined expanded layer structures (perhaps three layer structures if we consider the incipient plateau in the region between 18.7 Å and 19.3 Å), displaying plateaus corresponding to formation of monolayer and bilayer water in the interlayer region.Norrish [1] had noticed that the crystalline swelling proceeded only to 19 Å.It is accepted that the coexistence of different hydration states (0W, 1W, and 2W) in a smectite sample is common even under controlled conditions.It is also accepted that fractional hydration states correspond to a sample with different integer hydration structures.X-ray diffraction patterns provide evidence that hydration structures evolve gradually from one hydration state to the other through mixed-layer structures composed of discrete hydration states (Ferrage et al., 2005).It could perfectly be the explanation for experimental data outside of the regions defined by the dotted lines in Figure 3(a), which might have been originated from a mixture of montmorillonites at different hydration states.The basal spacing for K-montmorillonite after hydration obtained from our simulations is presented in Figure 3(b) along with the experimental data obtained by Calvet [3].The error bars represent two standard deviations.The results deviate from the experimental observations.However, similar simulations carried out by Suter et al. [11] reveal the same trend.
Figure 4 compares the swelling and hydration energy curves for the Na-, K-, and Ca-montmorillonite models.Comparing the swelling curves of K + and Na + , both monovalent cations, K-montmorillonite expands more than Namontmorillonite.This result is consistent with K + having a larger ionic radius than Na + , 1.02 Å compared to 1.38 Å.In turn, comparing the swelling curves of Ca 2+ and Na + , whose ionic radii are similar (1.00 Å compared to 1.02 Å), one can see that their expansion is fairly similar but there is a window in which Ca-montmorillonite swells more than Namontmorillonite.This result is consistent with Ca 2+ having a larger hydration energy than Na + .
Na-montmorillonite presents well-defined transitions between 0W-1W and 1W-2W layers, although the 2W layer is less evident than the 1W, being characterized by a slightly smaller slope than that of the expanded region at higher water content.K-montmorillonite also presents well-defined transitions between 0W-1W and 1W-2W layers.Ca-montmorillonite expands at lower water content than those with monovalent cations, which is consistent with predicted behavior on the basis of cation hydration energies.Experimentally, Ca-montmorillonite is known to form only two-layer hydrates [37].MD simulations, however, show a virtual 1W state for Ca-montmorillonite, state that has been also reported in previous works [38].

Distribution of Cations and Water
Molecules.Figure 5 presents the swelling curves and density profiles for selected water contents, illustrating the formation of 0W, 1W, and 2W.For dry montmorillonites, the distribution of Na + and Ca 2+ is similar, with cations attached to the clay surfaces, owing principally to their small size.K + , on the other hand, is located in the middle plane of the inlayer space because of its larger ionic radius, large enough to not fit within the hexagonal structure formed by the tetrahedral arrangement.It would be the explanation as to why dry Kmontmorillonite has a larger basal -spacing than dry Naand Ca-montmorillonite.
The formation of a virtual 1W layer for Ca-montmorillonite occurs at a water content of 0.048 g H 2 O/g clay and corresponds to a basal -spacing of 11.8 Å. Due to its larger hydration energy, the early onset of the 2W layer is favored, being completely distinguishable at a water content of 0.135 g H 2 O/g clay.For both Na-and K-montmorillonite, the 1W layer is completely formed at a water content of 0.085 g H 2 O/g clay.The distribution of cations and water molecules is quite similar for the three types of montmorillonites considered.The water layer is located in the middle plane of the interlayer space, while some cations are intercalated among water molecules, although the majority are still located between the lower clay surface (the one that holds tetrahedral substitutions) and the water layer, which evidences the formation of inner-sphere surface complexes; that is, no water molecules are interposed between the clay surface and the cations it binds [14].
The profiles for the 2W layer of Na-and Ca-montmorillonite are similar, with most of the cations located between the water layers, in the middle plane of the clay interlayer space.Some cations are still between the lower clay surface and the lower water layer, although some water molecules are now in contact with the lower clay surface.It evidences the formation of outer-sphere surface complexes.The water and cation peaks in the density profile are sharper than those formed in K-montmorillonite; K + cations are distributed in a wider region in the clay interlayer space, with approximately half of the total located in the middle plane, between the water layers, being fully hydrated, while the other half is still located between the lower clay surface and the lower water layer, meaning that approximately half of the cations are hydrated but not fully hydrated yet.This is consistent with the smaller hydration energy of K + .Water molecules are also spread in a wider region in the clay interlayer space, with some molecules in contact with the clay surfaces, demonstrating the formation of outer-sphere complexes.

Proposed Mechanisms of Clay Damage
Clay damage can be expressed not only as clay swelling but also as instability of the structure.Clay flakes can swell and still be stable.The swelling process occurs while cations are hydrated in the interlayer space.Water molecules interact first with the cations, forming inner-sphere complex surfaces in all cases, as shown in Figure 6.It means that water molecules locate in the middle plane of the interlayer region, while pushing cations towards the clay surfaces, preferentially those negative charged sites.Na + and Ca 2+ are small enough to be able to bury within the hexagonal structure formed by the tetrahedral arrangement, whereas K + is larger so that it cannot bury.It causes the K-montmorillonite to swell more than Na-and Ca-montmorillonite at lower water contents.Water content keeps increasing and the higher hydration energy of Ca-montmorillonite favors the onset of the 2W layer at a lower water content than that of Na-and Kmontmorillonite.This is seen in the range 0.06 to 0.15 g H 2 O/g clay.The hydration energy is not as high for Na-and K-montmorillonite; therefore the 1W state lasts longer, which means water content increases without notorious increase in basal -spacing; in other words, water density increases in the clay interlayer space.At this point, only inner-sphere complex surfaces have formed, although the higher the water content, the more cations detach from of the clay surface.The 2W state begins to form, and cations travel to the middle plane of the clay interlayer space, becoming hydrated.Depending on the hydration energy of the cations, some water molecules attach to the clay surface, indicating that outer-sphere complex surfaces are forming.The higher the hydration energy of a cation, the more affinity it has for water molecules and therefore the smaller the amount of water molecules that are left to attach to the clay surface.Nonetheless, water molecules that bind to the clay surfaces are able to interact with them, diffuse into them, and ultimately produce unstable structures that break down into little pieces called fines.

Conclusions
The MD simulation results in this work aid in the understanding of the underlying mechanisms of clay damage.The MD simulation results are in very good agreement with experiment for Na-montmorillonite but not for K-montmorillonite, which is probably due to the clayff force field.We propose that clay damage occurs through (1) swelling and (2) instability of the structure.Swelling of clays under contact with water occurs in a stepwise manner, forming one, two, and perhaps three layers.The density distribution of cations and water in the clay interlayer space depends on hydration energy and ionic radius of the cations.The high hydration energy of Ca 2+ favors the formation of the second water layer at water contents much smaller than those for Na + and K + .Water molecules interact first with cations, hydrating them, and then the excess water molecules interact with the clay walls.Water layers screen the electrostatic interaction between dissimilar charges.In addition, water molecules that attach to the clay surfaces, helping form outer-sphere complex surfaces, are able to interact with the clay surfaces, diffuse into them, and ultimately cause their damage.

Figure 1 :
Figure 1: Schematic explanation of the Reed[20] mechanism for particle generation by mica alteration during exposure to lowpotassium brine (modified from Civan[39]).Reed assumed that the ionic radius of Na + was larger than that of K + .

Figure 2 :
Figure 2: Schematic representation of the clay/water unit cell.The production calculation uses a supercell of 64 unit cells.Notice the location of the tetrahedral and octahedral substitutions.Color code: yellow: silicon; red: oxygen; magenta: aluminum; green: magnesium; white: hydrogen; blue: sodium.

Figure 3 :
Figure 3: Swelling behavior of montmorillonite clays upon hydration.Comparison of calculated and experimental results.

Figure 4 :
Figure 4: Montmorillonite swelling (basal -spacing) (a) and hydration energy (b) curves as a function of water content from MD simulations at 1 atm and 300 K.

Figure 5 :
Figure 5: Swelling curves and density profiles for Na-(a), K-(b), and Ca-montmorillonite (c) as a function of water content from MD simulations at 1 atm and 300 K. Color code: black: clay surface; red: water oxygen; blue: sodium; ochre: potassium; green: calcium.