Mechanical Properties of Chemically Modified Clay

Serpentine clay minerals are found in many geological settings. The rich diversity, both in chemical composition and crystal structure, alters the elastic behavior of clay rocks significantly, thus modifying seismic and sonic responses to shaley sequences. Computation of the elastic properties is a useful tool to characterize this diversity. In this paper we use first principles methods to compare the mechanical properties of lizardite Mg3(Si2O5)(OH)4, a polymorph of serpentine family, with the new compounds derived by substituting Mg ions with isovalent elements from different chemical groups. New compounds are first selected according to chemical and geometrical stability criteria, then full elastic tensors, bulk and shear modulii, and acoustic velocities are obtained. Overall, the new compounds have a lower anisotropy and are less resistant to mechanical deformation compared to the prototype, thus providing valuable information regarding chemical composition and mechanical properties in these systems.

Serpentine clay minerals are Mg-rich hydrous phyllosilicates found predominantly in the Earth's mantle, and their physical properties have major implications in diverse area of geophysics. Valuable information on the tectonic history can be gained by analyzing the structure and elastic constants of these minerals 1 . Knowledge of velocities of sedimentary rocks is mainly used in seismic imaging, hydraulic fracturing, and drilling. Since the Earth's interior is not homogeneous, the propagation velocity of internal seismic waves depend on the physical properties of the rocks that they cross. Measuring the elastic constants directly in fine grained materials is very difficult due to the impossibility of isolating a single grain of clay; for this reason, theoretical simulations combined with experimental measures or extrapolations are useful tools to investigate mechanical properties of these systems.
The structures of the serpentine minerals are based on uncharged tetrahedral Si 2 (O b ) 3 (O a ) 2 and octahedral Mg 3 (OH) 4 (O a ) 2 units arranged in sheets in a 1:1 ratio, where O a and O b are the apical and the basal oxygen atoms, respectively (Fig. 1). The serpentine group is composed of three clay minerals with the same chemical formula, but different crystal structures. The arrangement of these sheets is responsible for the different species: lizardite and antigorite (planar shape) and chrysotile (tubular form). Lizardite is Mg-rich 1:1 trioctahedral layer mineral, ideally Mg 3 (Si 2 O 5 )(OH) 4 , with space group P31m and trigonal crystal structure. The first crystal structure refinement was reported by Mellini 2 in 1982, then other researches have performed structural studies on two polytypes of lizardite 3,4 . See Fig. 1 for the crystal structure.
A few experimental measurements on the mechanical properties have been reported by Mellini and Zanazzi 3 , Tyburczy et al. 5 and Hilairet et al. 6 and some theoretical studies have been done using density functional theory (DFT) calculations [7][8][9][10] . Lizardite was also readily synthesized in laboratory 11 . In addition to Mg, numerous samples of serpentine rocks contain traces of other elements like Co, Mn, Fe, Al and Zn 12 . Understanding the effect of chemical variations on the physical properties of prototype lizardite is useful in synthesis and design of new functional clay minerals and this is the focus of the present work. Here we perform a systematic first principles study of the chemical, structural and mechanical properties of isostructural new compounds derived by the chemical substitution in lizardite, where Mg in the octahedral sheet is replaced with one of the twelve elements from different chemical groups (alkaline earth metals, Be and Ca, transition metals, Ni, Mn, Fe, and Zn, post-transition metals, Al and Sn, metalloids, Ge and Te, and the nonmetal, S and Se) with the same oxidation state as Mg 2+ . We compute the full elastic tensor and the derived criteria for stability, including bulk modulus, shear modulus and acoustic velocities for the new structures and a detailed comparison is performed with the existing prototype.
The paper is organized as follows. In Sec. 1, we describe the computational method, stability criteria and key quantities used to analyze mechanical properties. In Sec. 2, we first validate the method on the prototype material, then we discuss our results on the formation energies and elastic properties of the derived systems. Sec. 3 closes the paper with a summary of the results.

Materials and Methods
All the calculations are performed within density functional theory (DFT) using the quantum espresso (QE) code 13 integrated in the AFLOWπ (Automatic Flow π) 14,15 -a portable framework for high-throughput first principles calculations -which allows calculation of elastic properties, dielectric functions,electronic and optical properties. The crystal structures of the prototype, lizardite, was taken from the refinement reported by Mellini 2 from American Minerologist Crystal Structure Database 16 .
After several tests (see Sec. 2), we used the revised Perdew-Burke-Ernzerhof generalized gradient approximation (PBESol) 17 as exchange-correlation functional and ultrasoft pseudopotentials for our calculations. The kinetic energy cutoffs for the plane waves and for charge density are 60 and 600 Ry, respectively, and a 5 × 5 × 5 Monkhorst-Pack k-point mesh is used for all our calculations. The structures are relaxed until the convergence for forces and total energy achieve values lower than 0.1 mRy/Å, 0.01 mRy, respectively.
Lizardite has a trigonal (space group: P31m) crystal symmetry with six independent second-order elastic constants (SOEC): = c c 11 22 , c 12 , c 13 , c 33 e c 44 in Voigt notation ( → αβ c c ijkl ). We computed these values within the stress deformation approach as implemented in ElaStic code 18 where τ i and η j are the Lagrangian strain and stress, respectively. For the given space group and the Lagrangian strain values between − .
. [ 0 005, 0 005], around 17 distorted structures are generated and each of them is then relaxed. The deformation types are defined according to Yu et al. 19 . We obtain the SOECs, bulk (B) and shear (G) moduli for Voigt (V) and Reuss (R) averaging procedure, Young's modulus (E), Poisson's ratio (ν) and SOEC eigenvalues all of which are further used to analyze the mechanical stability of the new compounds. We use Drucker's criteria 20,21 for testing the mechanical stability. Additionally, the two modulii B R and G R values should also be positive for elastic stability. As a post-processing step, we obtain Pugh's index 22 , B/G, and the universal elastic anisotropy index 23 that includes all crystal symmetries, also used as mechanical stability criterion for the new compounds. The compressional wave velocity, V p , and shear wave velocity, V s are also calculated according to Birch 24 , www.nature.com/scientificreports www.nature.com/scientificreports/ where E(A ) Mg is the total energy of the new structure with element A replacing the Mg position in lizardite chemical structure and E X is the energy of the constituent X (A, Si, O 2 and H 2 ). The reference energies for calculating the chemical potential of the constituent elements are the bulk solids. For the prototype lizardite our calculated formation energy is −6154.94 kJ/mol in the same range as reported by a few experiments 11,31 . All the new compounds have a negative formation energy suggesting chemical stability. Compared to the prototype, E f decreases for Be, Ca and is higher for the transition metal ion substitutions Mn, Fe, Zn and Ni. For verifying if the compound is indeed stable, we also analyze the geometry and elastic properties to screen for elements which show both mechanical and chemical stability.  www.nature.com/scientificreports www.nature.com/scientificreports/ The lattice parameters and the interatomic distances for all the twelve A 3 (Si 2 O 5 )(OH) 4 relaxed systems are computed. Since we substitute the atom in the octahedral site (see Fig. 1), we look at the distortions in the cage. We find that only five substitutions (Ca, Mn, Ni, Fe, and Zn) preserve the octahedral structure. Furthermore, in case of Al, Sn, Ge, Te, S, Se and Be our computed values of the elastic properties show that one or more eigenvalues of elastic tensor is negative, thus contradicting Drucker's stability criterion described in the Section 1. In the rest of the paper, we focus on the elements that fulfill the criteria of chemical, geometric and elastic stability: Ca, Mn, Fe, Ni and Zn. The values of the lattice constants and elastic constants for these cases are presented in Tables 2 and 3.
The variation in the lattice parameter for the elements is consistent with the ionic radius 32 of the substituting element with the exception of Mn. The Mn 2+ ionic radius is smaller, 0.83 Å compared to 0.89 Å in Mg 2+ , which should reduce the volume of the cell but we find an increase in both a and c by 2%. The largest change is observed for Ca 3 (Si 2 O 5 )(OH) 4     www.nature.com/scientificreports www.nature.com/scientificreports/ elastic constants and Anisotropy. Elastic constants are fundamental properties of materials and can be used to judge the mechanical stability with different chemical compositions or crystal structures. In this work we compute the full elastic tensor for the five chemically modified lizardite compounds. In Table 3, we present our calculated results of the averaged bulk modulus B, shear modulus G, Young's modulus E and Poisson's ratio (ν). Figure 3 shows the Young's modulus in direction 1 and 2 as well as the shear stiffness G 12 and G 13 . The calculated universal anisotropic index A U is also given which is a measure of the anisotropic degree of the crystal as described in Sec. 1. For the prototype, our computed value of bulk modulus = .
B 72 46 GPa and the experimental values presented in the literature are single crystal X-ray diffraction studies (XRD) ( = B 57 GPa) 3 , shock wave equation of state ( = . ± . B (63 5 3 5) GPa at low pressure) 5 and synchrotron X-ray diffraction applying to P-V equations of state ( = . B 71 0(19) GPa for volume = 180.92 Å 3 ) 6 . Both B and G do not vary much with chemical modifications. The value of Poisson's ratio, ν, which is a measure of the ductility/brittleness of materials varies from 0.28 in Mg to 0.31 in Mn, indicating that the prototype and all the derived new compounds are more brittle, in accordance with Frantsevich's criteria 34 .
One important mechanical property for modeling correctly rock and fluid processes is the anisotropy of elasticity which is determined by computing the full elastic tensor, as discussed below. As shown in Table 5, the ratio of = . c c / 250 11 33 which means that in-plane stiffness is much larger than the out-of-plane c-axis, under a uniaxial stress. Also, = . shows that the basal plane is more resistant to fracture (rigidity) than c-axis when a shear stress is applied. These values are in good agreement with experimental results 6 and indicate strong anisotropic materials. Our results are different from other DFT calculations with LDA pseudopotentials done by Mookherjee and Stixrude 9 (20% for c 11 /c 33 and 27.4% for c 66 /c 44 ) and atomistic simulations presented by Auzende et al. 10 (100% for c 11 /c 33 and 11.4% for c 66 /c 44 ). The difference arises from the choice of the LDA functional which underestimates the lattice constants and consequently affects the elastic constants.
Additionally, the anisotropy factor can also be estimated from the compliance tensor, s ij , and the Young's modulus. From our calculated elastic constants, = = . s s E E / / 2 22 33 11 1 3 , the ratio between the Young's moduli is not the same. The lizardite anisotropy can also be seen from linear compressibility, β i , defined as material's response to decrease in length of a line when the crystal is subject to hydrostatic pressure in direction = and β β = .
/ 355 3 1 6 . Poisson's ratio of lizardite-chrysotile serpentinites average is 0.36 at 1 GPa 1 and our result is 0.28 GPa at zero pressure. High anisotropy of the lizardite is also shown in the Poisson's ratio, ν ij , that describes the response in the direction orthogonal j to this uniaxial stress i. We report for lizardite ν = − = . s s / 034 12 21 11 , ν = − = . s s / 013 13 31 11 , ν = − = . s s / 006 31 13 33 . Finally, we also compute the universal anisotropy index parameter A U as described in Sec. 1. For the prototype composition, the value is around 4.22 which is higher than most inorganic compounds 23 and is further proof of strong crystal anisotropy in this mineral class. Our computational approach describes all the elastic properties for the prototype lizardite in good agreement with experiments and provides references to analyze the new isostructural compounds. Compared to the prototype, A U decreases with the chemical substitutions, with Fe-lizardite being the least anisotropic among the list (see Table 3). The relations between the compressional elastic constants > c c 11 33 and the shear elastic constants > c c 66 44 , shown in Table 5, are lower compared to lizardite, suggesting higher stiffness and lower rigidity in the new compounds.
In Table 5, we tabulate the Hill-averaged shear (G) moduli, bulk (B) moduli and Pugh's indexes (B G / ). The smaller G values for the five compounds indicate low resistance to stress deformations. For Fe, Ni and Zn, B is slightly higher which means these materials have high fracture strengths. The higher Pugh's indices (B G / ) in all the new compounds, suggest that they are more ductile than the prototype lizardite.
One other parameter that can be extracted from our first principles elastic constants is the ratio of compressional (V p ) to shear wave velocity (V s ). This information is often used as a lithology indicator and is important for petrophysical evaluation, seismic imaging and modeling geomechanical properties. We calculated the compressional, V p , and shear, V s , wave velocities using Eq. (2) and the results are presented in Table 5, together with the DFT results obtained by Mookherjee and Stixrude 9 and with atomistic simultations reported by Auzende et al. 10 .
Our calculated values for the ideal protoype lizardite composition lies in the range of experimental values ranging from 5.03 to 6.33 km/s for V p and 2.62 to 3.38 km/s for V s as reported for three different samples by Watanabe et al. 35 . With chemical variations, we see a decrease in two velocities but the ratio increases. The difference between the compressional velocities of waves in lizardite and the new compounds is around 1 km/s and between shear velocities is around 0.7 km/s, therefore, they should be detectable in practice.  Table 5. Calculated density (ρ in 10 3 kg/m 3 ), Hill-averaged bulk (B) and shear (G) moduli (in GPa), Pugh's index (B/G), elastic constant ratio c ii /c jj , compressional, V p , and stress, V s , wave velocities (in m/s) of lizardite and its chemical substitutions (A 3 (Si 2 O 5 )(OH) 4 ), separated by group. The letter A is the element substituted in the lizardite structure. a Mookherjee and Stixrude 9 with DFT.