György Hantal

Laurent Brochard

Hadrien Laubie

Davoud Ebrahimi

Pellenq

Franz-Josef Ulm

Benoit Coasne
email: coasne@mit.edu

Atomic-scale modeling of elastic and failure properties of clays 1 2

The elastic and failure properties of a typical clay, illite, are investigated by means of molecular simulation. 2

We employ a reactive (ReaxFF) as well as a non-reactive (ClayFF) force field to assess the elastic properties 3 of the clay. As far as the failure properties are concerned, ReaxFF was used throughout the study, however 4 some calculations were also performed with ClayFF. A crack parallel to the clay layers is found to have low 5 fracture resistance (equivalent fracture toughness K Ic = 0.09 MPa.m 1/2 ) when submitted to a tensile loading 6 perpendicular to the crack (mode I). The nanoscale mechanism of both yield and fracture failures is 7 decohesion in the interlayer space, and the critical energy release rate characterizes both failures. In contrast, 8 under shear loading (mode II), the nanoscale failure mechanism is a stick-slip between clay layers. No fracture 9

propagation is observed as the clay layers slide on top of each other. The low fracture resistance in mode I and 10 the stick-slip failure in mode II are both the consequence of the lack of chemical bonds between clay layers 11 where the cohesion is provided by electrostatic interactions only. We also consider the mode I loading of a 12 crack perpendicular to the clay layers and find that, in this case, the material exhibit strain softening as the 13 result of the clay layers breaking one after the other. In this orientation illite displays a much higher fracture 14 resistance (0.61 MPa.m 1/2 ) due to the breaking of chemical bonds involved when fracturing in this direction. 15

This work, which provides a description of the failure properties of clays at the microscopic scale, is a first 16 step needed to describe the failure of clays at a larger scale where the polycrystalline distribution of clay 17 grains is a key parameter that must be taken into account. 18 19 20

Introduction 1

Oil and gas shale rocks have received particular attention in the past decades as source rocks for fossil fuel 2 production. (Vandenbroucke & Largeau, 2007) Due to the very low permeability of shale formations, the 3 common recovery technique is based on stimulating the rock by hydraulic fracturing. In particular, clay-rich 4 shales seem difficult to operate. Britt and Schoeffler report that shale with clay content higher than 40% are 5 seldom producing enough hydrocarbons to ensure profitability. [START_REF] Britt | The geomechanics of a shale play : What makes a shale prospective 1[END_REF]Various hypotheses 6

were proposed to explain those difficulties, e.g., low intrinsic permeability, incapability to activate natural 7 fractures, high in-situ stress, or swelling. Clearly, understanding the fundamental elastic and fracture 8 properties of clay-rich shale rocks is of key importance to operate in a more efficient and eco-friendly way. 9 10 Among major components in shales, clays are layered aluminosilicates that are also found in sediments, soils, 11 and argillaceous rocks. [START_REF] Bergaya | Handbook of clay science[END_REF]Clays have been used as raw material in pottery and 12 construction for millennia, but they have also been gaining importance in new technological fields over the 13 past decades. Modern applications of clays often rely on their swelling properties as they can uptake water or 14 selectively organic molecules and intercalate them between clay layers. Clays have also found applications 15 due to their large specific surface area which is responsible for a number of their unique physical and 16 chemical properties. Moreover, the fact that the adsorption and cation exchange capacity of clay minerals can 17 be controlled with the pH makes clays ideal candidates for a broad range of applications. Recently, they have 18 also been used as natural or artificial barriers in radioactive waste depositories due to their very low 19 permeability and long-term healing capacity. Clay-polymer complexes have also been used in inorganic-20 reinforced nanocomposites in fields ranging from aerospace systems to commodity plastics. (Vaia, 2002) 21 22 Owing to the numerous applications above, the structural and mechanical properties of clays have received 23 significant attention. However, because of their high affinity for water and their fine-grained polycrystalline 24 nature, measuring these properties is challenging. Indeed, clays usually appear in the form of aggregates of little crystals of a size smaller than 2 µm. Further difficulties are related to the fact that clay grains are often 1 mixed with other minerals, which renders obtaining reliable structural data, other than unit cell parameters, 2 also problematic. Even though their anisotropy is well known, very often only isotropic elastic properties of 3 clays are reported due to the limitations of experimental techniques in describing micro-scale mechanical 4 properties of clay grains. (Vanorio, Prasad, & Nur, 2003), (Faust & Knittle, 1994) 5 6

Molecular modeling provides a viable alternative to experimental techniques as it allows one to probe material 7

properties on the molecular scale. As a result, there have been several computational studies dedicated to the 8 investigation of the anisotropic elastic properties of clays. Heinz reported cleavage energies for different clay 9 The focus of this paper is on the elastic and failure properties of illite, one of the most common clay minerals 23 in Earth's crust. [START_REF] Bergaya | Handbook of clay science[END_REF] In particular, we aim to determine on the molecular scale its 24 elasticity tensor as well as its resistance to fracture and failure mechanisms under mechanical loading. We first 25 determine the elastic properties of illite with two different molecular force fields, and compare their ability to 26 predict the elasticity tensor. We also determine the failure properties of illite for different crack orientations 1 based on an empirical reactive force field. To the best of our knowledge, this study is the first to determine 2 fracture properties of clays by means of a reactive molecular force field. 3 4 2. Computational details 5

Clay model: illite 6

Illite is one of the most commonly found clay minerals in soils. [START_REF] Bergaya | Handbook of clay science[END_REF] It is also a 7 typical constituent of sedimentary rocks, especially shales. Illite, or more precisely the few distinct clay 8 minerals that are referred to as illite, constitute a sub-group in the dioctahedral 2:1 phyllosilicate family. For 9 the sake of simplicity, in what follows, we will use the term illite to refer to one typical composition. The 10 structure of illite is illustrated in Figure 1a. Illite is defined as a non-swelling clay with exclusively potassium 11 ions in the interlayer space. Pyrophyllite is considered the reference structure for dioctahedral 2:1 12 phyllosilicates as it forms charge-neutral clay layers where all tetrahedral (T) and octahedral (O) positions in 13 the structure are filled with Si and Al, respectively. Illite differs from pyrophyllite in that (i) about one 14 tetrahedral silicon position per unit cell is filled with aluminum and (ii) aluminum in one octahedral position 15 of every second unit cell is substituted in an isomorphous way, with other metal atoms. Such metal atoms can 16 be either of lower valence than Al (e.g., Mg(II) or Fe(II)) or, less typically, of the same valence (e.g., Fe(III)). 17

In case Al is replaced with a metal atom having a lower valence, the resulting charge defect is compensated 18 with potassium ions. The typical charge of clay layers in illite varies between -1e and -1.8e per unit cell 19 depending on the number of isomorphic substitutions of lower valence. In this study we used the composition 20 The starting configuration for illite used in our work is the pyrophyllite unit cell optimized from high-24 accuracy quantum chemistry calculations by Refson et al. (Refson, Park, & Sposito, 2003) In the simulation 25 cell, Si(IV) → Al(III) isomorphic substitutions were made in the top or bottom (T) layers of every replicated 1 unit cell together with Al(III) → Mg(II) substitutions in the (O) layer of every second unit cell. We ensured 2 that no O atom had more than one substituted metal neighbor in order to satisfy the Lowenstein rule (Bergaya, 3 Then, & Lagaly, 2006)). We also verified that an equal number of substitutions were made in the top and 4 bottom (T) layers to avoid highly asymmetrical charge distributions and, hence, large electrostatic fields in the 5 system. Potassium ions were distributed randomly in the clay interlayer space to counterbalance the 6 electrostatic charges induced by the isomorphic substitutions in the layers. applied to the system with an elliptical crack; while a tensile loading perpendicular to the plane of the crack 14 (xy plane here) acts on the system upon mode I ("opening mode"), a shear loading is applied parallel to the 15 plane of the crack upon mode II ("sliding mode"). The crack is infinitely long in the y direction and 16

propagation is expected to occur in the xy plane. 17

Elastic properties 1

There are different ways to obtain the elastic properties of solids from molecular simulation. A widely used 2 approach is based on applying a series of small strains to the system and calculating the resulting stress tensor. 3

If the deformation is small enough and the system remains in the linear elastic regime, the stress-strain curve 4 can be fitted with a straight line whose slope provides the stiffness parameters. The advantage of this method 5 is that is can be used at any finite temperatures, however, its application is often not possible as the result may 6 exhibit strain rate dependence. Another approach to determine the elastic properties of solids is to estimate the 7 curvature of the potential energy surface about a minimum energy configuration. Indeed, the components of 8 the symmetrical stiffness tensor C ij (in Voigt notation) can be written as the second derivative of the energy 9 density with respect to the strains ε i and ε j : 10

𝐶 !" = ! ! ! ! ! !! ! !! ! , (1) 11
where V is the volume of the simulation cell and U is the potential energy. ε i and ε j (i,j = 1…6) represent the 12 components of the second order strain tensor in Voigt notation. Here we use a Cartesian coordinate frame 13 where direction 3 is perpendicular to the clay layers. For simple potential energy functions the second 14 derivatives of the energy are often available analytically. In such a case the stiffness tensor can readily be 15 determined at any point of the potential energy surface. If the functional form is too complicated, as in the 16 case of the ReaxFF potential used in the present work, the 21 independent elements of the stiffness tensor can 17 be calculated by applying linear combinations of all elementary deformations and performing an energy 18 minimization of the resulting structures. The obtained energy profiles are then fitted with a parabola which 19 provides a numerical estimate of the curvature of the potential energy surface about the energy minimum. 20 Such a curvature is then used to estimate C ij from Eq. 1.Note that strictly speaking this method provides 21 elastic parameters at 0K. If the entropic contribution to the free energy is small C ij values determined at 0 K 22 are not very different from values at finite temperatures. This is often the case for minerals with high melting 23 point, like clays. However this assumption does not hold for solids of low melting point (e.g. polymers). 24 25

Failure properties 1

Fracture failure refers to the failure of materials where a flaw, like a micro-crack, is responsible for the 2 occurrence of a crack in a localized region of a body. By contrast, yield failure refers to the failure of flawless 3 materials. In this work, we assessed both failure properties by simulating illite with and without an initial 4 micro-crack. According to fracture mechanics, at the continuum scale, the resistance of a material to fracture 5 is characterized by a few material properties. [START_REF] Anderson | Fracture Mechanics: Fundamentals and applications[END_REF] In the simplest formulation, valid for linear 6 elastic materials under the assumption of small scale yielding, a single property is sufficient, which is either 7 the stress intensity K or the energy release rate G. When the stress intensity reaches a critical value, called 8 fracture toughness K c , the crack propagates (or, equivalently, when the energy release rate reaches a critical 9 value G c ). The usual assumption is that K c and G c are material properties. Regarding yield, the failure is 10 characterized by the maximum stress conditions, the yield strength or yield surface, that the flawless material 11 can sustain. In the general case, yield strength and fracture resistance are not correlated. In this work, we focus 12 on the particular fracture and yield properties of illite for which a peculiar yield-fracture relation is possible 13 because of the layered structure. 14 15 In the present work, we use a thermodynamic route to derive the critical energy release rate, in which the 16 variation of free energy of the system is estimated from molecular simulation (for a detailed description of the 17 method the reader is advised to refer to Ref (Brochard, Hantal, Laubie, Pellenq, & Ulm)). In isothermal 18 conditions, the thermodynamic state of a deformable solid can be described by the applied mechanical loading 19 and, if there is, the crack area. Upon loading, the mechanical work provided to the system and stored in 20 recoverable energy is released when fracture propagation occurs. Our goal is to determine by means of 21 molecular simulations the released energy when the fracture propagates, i.e. upon increasing the surface A crack 22 of an existing crack by a quantity ΔA crack . The balance of Helmholtz free energy F of the system under 23 displacement loading is: 24

𝑑𝐹 = 𝑑𝑃 + 𝐺 ! 𝑑𝐴 !"#!$ where 𝑃 = 𝜎(𝜀): 𝑑𝜀 !(r) ! dr ! (2) 25 1
The change in the free energy dF includes the variation in the mechanical energy dP and the energy released 2 by propagating the crack G C dA crack . The mechanical energy can be locally evaluated by integrating the stress 3 σ(ε) over the local strain dε(r) at any point r. A second integration over the whole volume of the system V 4 yields the total mechanical energy stored in the system. In these equations the double-underlined letters 5 correspond to second-order tensorial quantities. In isothermal equilibrium conditions, the free energy is 6 minimum dF = 0, which leads to the usual definition of the critical energy release rate G c : 7

𝐺 ! = !!" !" !"#!$ !,!"#$%&' , (3) 8
where the derivative is at constant temperature and loading. Akin to real experiments, we apply in molecular 9 simulations a load to the system and measure its mechanical response. In our study, this is done by imposing 10 the displacement at the boundary of the system (represented by the global strain 𝐸) and measuring the 11 resulting stress 𝛴 at the boundaries of the sample.The integration of 𝛴 with respect to 𝐸 provides an estimate 12 of the variation of free energy of the system per unit volume. If, at the end of the simulation, the system is in 13 the same loading state as in the beginning, the mechanical energy has been totally released by crack surface 14 creation so that the integration of the stress-strain curve yields the critical energy release rate G C : 15

𝐺 ! = 𝑉 ! ! :d! ! !"#$% ! d! !"#!$ ! !"#$% ! , (4) 16
where

d𝐴 !"#!$ ! !"#$% ! = ∆𝐴 !"#!$
is the area of crack created. It is possible to relate the energy release to the 17 stress intensity factor through Irwin's formula. In a general case, for the linear combination of all three 18 loading modes the generalized Irwin formula can be used: 19

𝐺 = 𝜋 {𝐾} ! • [𝐻] • {𝐾} (5) 20
where {K} = (K I ,K II ,K III ) is a vector containing the stress intensity factors of the three loading modes (I = 21 "opening", II = "sliding", and III = "tearing" modes). [H] is a 3×3 matrix that depends on the elastic properties 22 of the material. t {} denotes the transpose of a vector. While finding the [H] matrix is usually not 23 straightforward, its determination can be simplified in case of simple loadings or material symmetries. Here 1 we will use the solution derived for plane-strain conditions, under which all strains in the material are 2 restricted to lie in a plane. If a strain acts in the xy plane, only ε xx , ε yy and ε xy are non-zero. In the general case 3 of plane strains, there always exists a coordinate frame in which the strain tensor has only 3 non-zero 4 elements, which reduces the 3D problem to a 2D one. 5 6

Evaluation of the crack surface area 7

As shown in the previous section, evaluating the fracture energy G c requires the determination of the area of 8 the fracture created upon releasing the mechanical energy (see Eq. 4). This is far from being a trivial task in 9 molecular simulations due to the fractal-like character of atomic surfaces which is a consequence of the 10 discontinuous nature of matter at this scale. All molecular surface analysis methods necessarily have an 11 adjustable parameter that characterizes its sensitivity to the size of discontinuities and, hence, the resolution of 12 the given method. There is thus no unique way to evaluate surface area in molecular simulations .Among the 13 various surface analysis methods [START_REF] Lee | The interpretation of protein structures: estimation of static accessibility[END_REF], [START_REF] Connolly | Analytical molecular surface calculation[END_REF], (Jorge, Jedlovszky, & Cordeiro, 14 2010),we chose to apply a methodology similar to that yielding the solvent-accessible surface, (Lee & 15 Richards, 1971) also known as the specific surface. In this algorithm a probe atom of a certain size rolls over 16 all surface atoms. The specific surface is defined as the set of points visited by the center of the probe during 17 rolling. In our study, this surface is estimated by a Monte Carlo procedure; (Bhattacharya, Coasne, Hung, & 18 Gubbins, 2009) segments are generated in the system with random position and direction that may intersect 19 the surface. The specific surface area (S SS ) of the void-wall interface per unit volume is then calculated as 20

𝑆 !! = 2 ! !!! ! ! ! ! !!! ! ! ! , ( 6

) 21

where p i is the number of intersections of the segments i (=1…N) of length l i with the interface. The 22 algorithm works as follows. A random location is picked in the system, which is the initial point of the probe 23 segment. If a probe sphere of radius r p can be inserted at this point without overlapping with any atom, this 24 point is considered as belonging to the "void". In contrast, if the probe sphere overlaps with any of the atoms, 25 the point is considered to belong to the "wall". The next point analyzed of the segment is at a distance Δl from 1 the previous point in a direction chosen randomly (but kept constant for the same segment). Again, this point 2 is characterized as belonging either to the "void" or to the "wall" by checking whether a probe sphere of 3 radius r p could be accommodated at this point without overlapping with the wall atoms. An intersection (i.e., a 4 surface point) is detected when two consecutive points on a given segment belong to different domains (either 5 wall or void). If the procedure is run for a large number of segments, it provides a representative sample of 6 surface points being located at distances d∈[r p , r p +Δl] from the wall atoms. The choice of the maximum 7 length of the segment as well as of Δl affects the efficiency of the procedure but not its sensitivity. The latter 8 is controlled by the radius of the probe sphere r p . There is no unique choice of this quantity, rather a 9 reasonably defined range of possible values. Our criterion for setting r p was to allow the crack opening to be 10 detected in the clay interlayer, i.e., it was such that when the system is intact (there is no fracture in the 11 system) the probe sphere is large enough not to fit into the clay interlayer. However it fits in (and hence 12 detects surface creation) once the crack starts to propagate. This led us to set r p = 3.2 Å. Note that during the 13 analysis, molecular fragments detached from illite as well as K + ions were excluded. Figure 2 illustrates the 14 way the surface analysis method above works together with a crack and an example of a set of sampled 15 surface points. 16 available parameterizations for various kinds of molecular system. We used the parameter set which was 9 developed by Pitman and van Duin (Pitman & van Duin, 2012)to study the reactivity of water confined 10 between zeolite and smectite. As this set does not contain parameters for potassium ions, the original 11 parameters for sodium ions were modified by replacing those that describe basic physical properties such as 12 the atomic mass and sigma bond covalent radius, as well as non-bonded interaction parameters described by 13 the Morse potential. These latter parameters were obtained to best fit the Lennard-Jones potential of K ions as 14 described by the ClayFF potential. (Cygan, Liang, & Kalinichev, 2004) These changes only affect the 15 intermolecular part of the interactions. Since potassium is not expected to participate in covalent bond 16 formation and only acts as a counterion, these modifications are thought to be satisfactory and relevant. As a 17 test of our results, some simulations were also performed with the non-reactive ClayFF potential, which is 18 known to accurately reproduce structural data of various clay phases. This is due to the fact that most 19 After setting up the initial configuration, the system was first energy minimized. A series of deformations was 2 then applied; each deformation was followed by an energy minimization and then by a molecular dynamics 3 simulation at constant volume and temperature for 10 ps. Equilibrated configurations were sampled and saved 4 for further analysis. Periodic boundary conditions were used in the three directions x, y, and z. The elements of 5 the stress tensor were calculated with the virial expression and averaged for every deformation step: (Allen & 6 Tildesley, 1989) 7

𝛴 !" = - ! ! ! ! ! ! ! ! !,! ! !,! ! ! ! !!! + 𝑟 !" 𝑓 !,! ! !!! , (7) 8
where a, b, and c denote the three box vectors and the subscripts x, y, and, z refer to vectorial components 9 along the three directions in the Cartesian coordinate frame. α and β represent any of these components. The 10 index i (=1,..,N) runs over all N atoms in the system. m i , p i , and, r i correspond to the mass, momentum, and 11 position of atom i, respectively. f i is the force acting on atom i. 12 13 3. Results and discussion 14

Elastic properties of illite 15

We used the numerical method (introduced in Section 2.3) to obtain the stiffness tensor of illite with the 16

ReaxFF forcefield (presented in Table 1). Figure 3 illustrates the fitting procedure of the obtained energy 17 profiles against a parabolic function, which allows estimating for each type of deformation the stiffness 18 constant C ij using Eq. ( 6). The elasticity tensor of illite was also computed with the ClayFF forcefield (shown 19 in Table 2) using the GULP 3.1 program package [START_REF] Gale | GULP 3.1, General Utility Lattice Program[END_REF] which provides the full stiffness tensor from 20 the analytical derivatives of the intermolecular interactions after a thorough energy minimization of the 21 system. All values in Tables 1 and 2 were averaged over 5 equilibrium configurations. Both matrices are 22 positive definite, which ensures that the energy minimized systems correspond to stable structures. In this specific example, the deformation ε 3 is along the direction 3. These data were obtained at 300 K using 3

ReaxFF. 4 5

Table 1. Stiffness tensor of illite in Voigt notation as obtained with ReaxFF. Note that the tensor is shown 6 here as an upper triangular matrix due to symmetry reasons (C ij = C ji ). The shaded fields in the table show the 7 elements that are expected to be different from 0 for a transversely isotropic material (see text). 8 It is useful to consider the random polycrystal properties, i.e., the macroscopic elastic moduli, which can be 5 obtained in the framework of the Voigt-Reuss-Hill (VRH) approximation classically used in mineralogy. In 6 this approach one simply takes the average of the values obtained based on the definitions given by Voigt 7 (using the stiffness tensor) and those defined by Reuss (using the compliance tensor which is the inverse of 8 the stiffness tensor, S = C -1 ). Tables 3 Values in column 3 of Tables 1 and2 are exceptions since ReaxFF seem to predict stiffer behavior along the 3 direction perpendicular to the clay layers (direction 3). It is worth noting that, considering the stiffness tensor 4 (C 11 > C 22 > C 33 , see Table 1), the trend in the Young's moduli calculated with ReaxFF (E 2 >E 1 >E 3 , see Table 3) 5 is counterintuitive. However, this is not observed in the case of ClayFF (cf. Tables 2 and3). Note also that the 6 error computed with ReaxFF is two orders of magnitude larger than that of ClayFF. 7 8

C ij /GPa 1 2 3 4 5 6 1 216±5 76±9 29±4 -4±3 -21±2 -3±2 2 209±4 29±5 -1±3 -9±2 0±2 3 93±1 -7±1 -23±1 -4±1 4
We also examined the anisotropy of the elasticity tensor. It is common practice to consider clays transversely 9 isotropic (TI) materials. For such materials the elastic properties show a specific rotational symmetry: in 10 planes perpendicular to the axis of the rotational symmetry, elastic properties are isotropic (in Figure 1, this is 11 the z axis). Due to its high symmetry, the elastic properties of a TI material obey the following relationships: 12 relations hold for the compliance tensor). A comparison of the obtained elasticity tensors of our models to that 15 of an ideal transversely isotropic material provides a test of the reliability of the computed tensors. 16

Considering the data obtained for illite with ReaxFF ( 1 Further comparison of the elasticity tensors obtained using ReaxFF and ClayFF can be achieved by measuring 2 their distance from the tensor of a TI material using an appropriate distance metric. Despite its simplicity and 3 elegance, the well-known Euclidean metric sometimes leads to ambiguities in the interpretation of results. 4

Here we chose to use the Riemannian metric as it avoids such ambiguities. (Moakher, 2006)The Riemannian 5 metric between two positive definite n×n matrices, A and B, is defined as:

6 𝑑 R (A,B) = ln (𝐵 !/! 𝐴 !! 𝐵 !/! ) = ln ! 𝜆 ! ! !!! !/! , (8) 7
where λ i are the eigenvalues of A -1 B. For each computed stiffness tensor, a reference TI matrix was generated 8 and the distance of this reference to the original tensor was calculated using the Riemannian metric. that the Riemannian distance of the tensors calculated with ReaxFF and ClayFF to their reference matrix is 13 2.66 and 0.37, respectively. In other words, ClayFF predicts elastic properties that are closer to those of a TI 14 material. To quantify the difference between the two computed matrices and the reference TI matrices, it is 15 also possible to calculate the error one would make by using the reference TI matrix instead of the original 16 one when calculating fracture toughness. For cracks oriented parallel with the plane of transverse isotropy 17 (cracks in the 001 plane, see below), this comparison yields 108.5% error for ReaxFF but only 2.6% for 18

ClayFF. 19 20

Depending on the applied molecular force field, elastic constants, especially in-plane constants, have been 21

found to differ significantly. (Zartman, Liu, Akdim, Pachter, & Heinz, 2010) Although desirable, the 22 comparison of our results with experimental data is prohibited due to the lack of experimental values. 23 Therefore, we use the above qualitative comparisons to discriminate between the two sets of data: we 24 conclude that our molecular simulations with ClayFF give a reliable description of the elastic properties of 25 illite, due to the close resemblance of the ClayFF elasticity tensor to that of a TI material. In contrast, ReaxFF 1 gives a poor description of illite as a TI material. We will thus use, in the rest of this study, the reliable 2 elasticity tensor computed with the ClayFF potential and combine with fracture energies obtained from 3

ReaxFF when reactivity needs to be taken into account. 4 5

Failure parallel to clay layers 6

To study both yield failure and fracture failure, we simulate bulk illite both with and without an existing 7 crack. As the clay layers are not bonded covalently, the basal (001) plane is an ideal plane for cleavage for 8 both fracture failure and yield failure. Thus, we focus our investigation on the possible relationship between 9 fracture and yield. Figure 1b illustrates a crack lying in the xy plane (i.e., the basal plane). In the ideal case the 10 thickness of the crack along the z direction must be negligible compared to its dimensions in the other 11 directions. In our simulations, elliptical cracks with high aspect ratios (≈ 4-5, i.e., the ratio of the major and 12 minor axes of the ellipse) were created in the interlayer region. We consider two crack sizes of lengths 32.5 13 ands 53 Å and an intact structure (i.e., without a pre-existing crack). For the intact structure, we consider three 14 different system sizes besides the 'original' size: a system 50% longer or 50% thicker. Such we could assess 15 the dependence of the results to crack size and system size. 16 17

Mode I loading 18

The calculated stress-strain curves display the same qualitative behavior upon mode I loading, which confirms 19 that it is independent of the studied crack and system size. As an illustration, we present the stress-strain curve 20 obtained in a system without an initiated crack in Figure 4. This figure shows that illite exhibits very low 21 resistance to mechanical loading characterized by low values of the critical strain (<0.05) and critical stress (< 22 2.5 GPa). It is also seen that the material is in the linear regime at low loadings, which is followed by a 23 nonlinear strain hardening while approaching the point where the material fails and the stress drops to zero. 24

Irrespective of the presence of an initial crack in the system, all the systems fail by cleavage in the interlayer 25 space. Accordingly, the nanoscale mechanism of fracture failure and yield failure are qualitatively similar. As 1 shown in Eq.( 3), one has to estimate the crack surface area in order to estimate the fracture energy G c from the 2 integral of the strain-stress curves. To evaluate the crack area, the method introduced in Section 2.4 was 3 applied. In Figure 4, besides the stress-strain curve, the crack area in the same system is also plotted as a 4 function of the strain imposed. These data show that the chosen probe size is appropriate to describe fracture 5 propagation as the change in the crack surface area (i) is negligible in the elastic regime and (ii) increases 6 abruptly as the fracture starts propagating. For both intact and initially cracked systems, the energy released is 7

proportional to the area of uncracked material. This proves that fracture failure and yield failure originate from 8 the same mechanism: a decohesion in the interlayer space which can be characterized by an energy released 9 per unit area, that is, by an energy release rate (Eq. 4). The nanoscale heterogeneity of the clay leads to a 10 peculiar relationship between yield and fracture failures. Usually, in perfectly homogenous material, yield 11 would distribute over the volume with a homogenous strain field and the energy released would be 12 proportional to the volume; while fracture concentrates the stresses and strains in a plane. The average fracture 13 energy G c obtained by averaging over the five different systems considered in this work, is G c = 0.22 ± 0.04 14 J/m 2 . This value is quite low owing to the lack of chemical bonds between the clay layers. 15 were obtained for the smallest system without a pre-existing crack (the other systems display similar 1 qualitative results). The crack area was estimated using a probe r p = 3.2 Å (see Section 2.4). 2 3 It is worth mentioning that we found that systems thinner in the z direction than our thinnest system produced 4 unphysical stress oscillations after crack propagation. This is unrealistic as the system, once already broken, 5

shows increasing stress (and then also stress of opposite sign) upon further strain increase although there is no 6 more cohesion between the two broken parts of the fracture. Moreover, for such a small system size, the clay 7 layers were found to oscillate before and after crack propagation. We believe that this is most likely due to the 8 small size of the system considered. 9 10 As pointed out above, the pre-existence of a crack in the system is not needed for a fracture to propagate in the 11 direction parallel to the clay layers. This means that the molecular simulations above can be repeated with a 12 non-reactive force field such as ClayFF in order to confirm the failure behavior of clays. Note that this is only 13 true when the system does not contain an initial crack as a reactive forcefield is otherwise needed to account 14 for possible charge deficiencies as well as surface reconstruction. The obtained stress-strain curve calculated 15 with ClayFF is shown in Figure 5. This curve shows the same brittle mechanism as that obtained with 16

ReaxFF. However, illite simulated with ClayFF reaches higher strains before failing (ε = 0.1). Moreover, the 17 stress-strain curve obtained with ClayFF yields a fracture energy G C = 0.34 J/m 2 , which is about 50% larger 18 than the value obtained with ReaxFF. This is probably due to the lack of long range corrections of the 19 electrostatic interaction in the ReaxFF potential which leads to weaker cohesion in illite. 20 21 As shown in Section 3.1, the stiffness tensor computed with ClayFF can be very well approximated with that 22 of a TI material. Owing to this symmetry of the system, it is possible to provide a simple expression for H 11 23 (see Eq.( 5)),the matrix element needed to convert energy release rate to stress intensity (all other elements are 24 irrelevant in pure mode I).The original expression by Sih et al. (Sih, Paris, & Irwin, 1964) has been rewritten 25 in terms of the stiffness constants. [START_REF] Laubie | [END_REF]This equation reads for plane-strain cases as: 26

𝐻 !! = ! !! ! !! ! !! ! !! !! !" ! ! ! !! + ! ! !" ! ! !! ! !! (9) 1 2
Using the stiffness matrix components in Table 2 one obtains fracture toughnesses of 0.12 and 0.097 MPa.m 1/2 3 for ClayFF and ReaxFF, respectively. These fracture toughnesses are very low compared to typical values of 4 metal oxides (~ 1 MPa.m 1/2 ), (Brochard, Hantal, Laubie, Pellenq, & Ulm), (Spray, 2010)which is the result of 5 the absence of chemical bonds between the clay layers where the cohesion is provided by electrostatic forces 6 only. Indeed, as will be seen in Section 3.4, a much larger fracture toughness is obtained when a crack 7 perpendicular to clay layers is investigated i.e., when chemical bond breaking is involved in the fracture 8 propagation. 9 10 Figure 5. Stress-strain curve obtained with ClayFF when illite is loaded in mode I in a direction parallel to the 11 clay layers. 12 13

Mode II loading 14

Again, all stress-strain curves obtained for different crack and system sizes are found to display the same 15 behavior. As a result, we present in Figure 6 only one stress-strain curve obtained upon mode II loading 16 parallel to the clay layers. The system in mode II loading exhibits a behavior that is radically different from 17 what was observed for mode I loading. After an initial stable increase of the stress upon increasing the load, 18 an instability occurs and the stress sharply drops. This drop is, however, followed by a new stable stress increase until the stress drops again. These successive stable and unstable regimes alternate upon increasing 1 the shear strain. This is characteristic of a failure with stick-slip. For the systems with crack, crack does not 2 propagate. The system does not fracture, but it releases mechanical energy as the clay layers slide on the top 3 of each other. At the macroscale, we expect the local stick-slips to be averaged which would led to an 4 apparent plastic behavior. Similar failure mechanism was identified by Zartman et mechanism, but the curve displays a much greater regularity than those obtained with ReaxFF. The stick-slip 1 mechanism is almost periodic. It is worth noticing that the stress value at which the slip of clay layers occur is 2 considerably higher with ClayFF than the highest value reached with ReaxFF (1.9 GPa vs 0.8 GPa, see 3

Figures 6 and7). This can be related to the missing long range correction in ReaxFF. 4

Fracture orthogonal to clay layers 10

In addition to cracks parallel to the bedding plane of clay layers, we also considered cracks perpendicular to 11 the 001 plane (i.e., to the clay layers). In this section we investigate the mode I loading of such cracks. Given 12 that fracture propagation in this direction requires breaking chemical bonds, the molecular simulations can 13 only be performed (i) with a pre-existing crack and (ii) by using a reactive forcefield (ReaxFF in the present 14 work). Figure 8 shows the computed stress-strain curve when the crack is loaded in mode I. The crack 15 propagation exhibits some small instabilities while the general behavior is a strain softening during failure. 16

The underlying mechanism is a sequence of the brittle breaking of individual clay layers, similarly to the 17 failure of layered composites with some delamination between the layers. This conclusion is supported by the 18 evolution of the crack surface area upon loading. The increase of the crack surface area between stress drops 19 configurations show that the clay layers break in a sequential manner. The crack surface area is evaluated by 8 the method described in Section 2.4. This curve reveals that the fracture propagation proceeds by first creating 9 separate secondary cracks that eventually merge. 10 11 After having evaluated the crack surface area created upon propagation of the crack in this direction, it is 12 possible to compute the corresponding energy release rate G C . Using Eq.(4) one obtains G C = 3.86 J/m 2 which 13 is one order of magnitude larger than what has been estimated for cracks oriented along the 001 plane. One 14 can also calculate from G C the corresponding fracture toughness value using an equation developed for a crack 15 in this specific crack orientation (propagating along direction 3) and loaded in mode I. For TI materials, this 16 expression reads: (Laubie, 2013) 17 18 1 For this crack orientation the fracture toughness was found to be K C = 0.61 MPa.m 1/2 . Table 4 compares the 2 fracture toughnesses obtained in this work for the different crack orientations. We also report experimental 3 value found in the literature for muscovite mica; (Spray, 2010) although somewhat different, this material is 4 often considered as a proxy for illite in terms of mechanical properties. Indeed, illite and muscovite differ 5 slightly in their number and distribution of isomorphic substitutions, (Bergaya, (We suppose that the same trend would be found when using the corresponding experimental value of illite.) 15

𝐻 !! = ! !! ! !! ! !! ! !! !! !" ! ! ! !! + ! ! !" ! ! !! ! !! (
The trend shows that the experimental fracture toughness for muscovite is about twice the value calculated for 16 the 001 crack orientation in illite, but three times smaller than the computed value when the crack is 17 perpendicular to clay layers. To understand this result, one has to bear in mind that illite (and clays in general) 18 is a polycrystalline material at the macroscopic level and upscaling the simulated microscopic results to the 19 macroscale is not straightforward. In other words the representative elementary volume for the macroscopic 20 fracture properties is larger than what has been studied here. 21 22 The direction of the fracture propagation is always determined by the lowest fracture resistance (fracture 1 toughness) which is based locally on the grain orientation and the far field stress. As the crack can be 2 subjected to various toughening phenomena when propagating in a granular material (such as crack bridging, 3 crack deflection, front roughening, crack pinning, etc), the overall macroscopic fracture toughness will always 4 be higher than the one in the preferred fracture propagation direction. Data in Table 4 suggest that fracture 5 propagation is likely to happen in a direction parallel to the clay layers (001 direction) and will always be 6 preferred over a direction perpendicular to them. Due to the possible toughening mechanisms mentioned 7 above, it is consistent to find the measured macroscopic fracture toughness (0.2 MPa•m 1/2 ) higher than the 8 value obtained in the possible preferred direction of the fracture propagation (0.1 MPa•m 1/2 ). 9 10 This work sheds light on the fracture properties of clays at the microscopic scale and, as such, can be regarded 11 as the starting point for a mesoscopic fracture simulation and eventually to upscale fracture properties to the 12 macroscopic level. This work is also the basis for investigation of more complex systems at the molecular 13 scale, such as illite-organic composite materials. (Hantal, Brochard, Pellenq, Ulm, & Coasne) 14 15

Conclusions 16

In this paper we studied the elastic and fracture properties of a typical clay mineral, illite, by means of 17 molecular simulations. We used two force fields: the reactive ReaxFF potential which allows us to simulate 18 bond breaking and formation during crack propagation and the non-reactive ClayFF potential which is known 19 to accurately reproduce the structural properties of various clay minerals. The elasticity tensor of illite was 20 determined using the two force fields and compared to that of a generic transversely isotropic material as an 21 example of ideal layered materials. 22 23 The performance of the two force fields in terms of providing a reliable description of a transversely isotropic 24 material was assessed based on three criteria: (i) the similarity of the calculated elasticity tensors to that of a 25 transversely isotropic material; (ii) the distance of the elasticity tensors from that of an ideal TI material, and 26 (iii) the magnitude of the error one would make by using the TI approximation to calculate fracture toughness instead of the general (less symmetric) formula. All these three comparisons indicate that ClayFF provides a 1 good description of the elastic properties of illite clay. In contrast, ReaxFF is found to fail at providing a 2 reliable description of the elastic properties. As a result, we decided to use the stiffness tensor determined with 3

ClayFF throughout the study to calculate fracture properties of illite. In contrast, fracture propagation 4 involving breaking and formation was studied using ReaxFF as it necessarily requires the use of reactive 5 forcefield. 6 7 Yield and fracture failures were also investigated. Failure in the 001 plane (i.e., the most likely cleavage 8 plane) was studied under tensile and shear loading, i.e. under mode I and mode II loading for cracked systems. 9

Various system and crack sizes were studied to make sure that our results are independent of the sizes of the 10 system and crack. In mode I, it was found that yield and fracture failures originate from a single mechanism: a 11 decohesion in the interlayer space characterized by an energy released per unit area of crack created, which is 12 the usual critical energy release rate in the case of fracture failure. The corresponding fracture toughness, 13 calculated with ReaxFF, was found to be quite low: K c ReaxFF = 0.09 MPa.m 1/2 . In contrast, in mode II loading, 14 the yield failure mechanism is a stick-slip between clay layers. Fracture failure follows the same mechanism 15 and no crack propagation could be observed. Again this is due to the lack of chemical bonding between the 16 clay layers so that they slide on top of each other. We repeated the simulation of the mode I loading of such a 17 crack with the non-reactive ClayFF potential. We observed the same molecular mechanism, but we estimated 18 a somewhat higher fracture toughness K c ClayFF = 0.12 MPa.m 1/2 . The mode I loading of a crack perpendicular to 19 the clay layers was also investigated. We found that, in this case, the material exhibits a strain softening as the 20 result of the clay layers breaking one after the other. Illite displayed a much higher fracture resistance (0.61 21 MPa.m 1/2 ) in this orientation. 22 23 The highly anisotropic failure behavior of illite was explained as a result of its specific layered structure. The 24 low failure resistance under tension and the stick-slip failure under shear are both the consequence of the lack 25 of chemical bonds between clay layers where the cohesion is provided by electrostatic interactions only. The 26 mode I fracture toughness values for the 001 orientation was found to be consistent with the experimental 27 value (0.2 MPa.m 1/2 ) measured for muscovite, a clay mineral analogous to illite. The measured enhanced 1 macroscopic fracture toughness value compared to our result was interpreted with toughening mechanisms 2 arising from the granular structure of clays. We wish to thank A. Saul (CNRS/MIT) for helpful discussions. 7 8

  minerals and clay/surfactant composites based on classical molecular dynamics simulations. (Heinz, 2006) 10 This author found good agreement with available experimental measurements. Soto et al. was the first to use 11 first-principles calculations to determine the elasticity tensor of a layered silicate, kaolinite. (Sato, Ono, 12 Johnston, & A., 2005) Similar first-principles calculations were performed by Zartman et al. and Militzer et 13 al. to obtain elasticity parameters of different clay minerals including pyrophyllite, montmorillonite, kaolinite, 14 and muscovite as a function of the number and the distribution of isomorphic substitutions. (Zartman, Liu, 15 Akdim, Pachter, & Heinz, 2010) (Militzer, Wenk, Stackhouse, & Stixrude, 2011) Teich-McGoldrick et al. 16 used the classical ClayFF force field to determine the pressure and temperature dependence of the elastic 17 properties of muscovite and found better agreement with experiments than earlier first-principles calculations. 18 (Teich-McGoldrick, Greathouse, & Cygan, 2012) Ebrahimi et al. also carried out classical molecular 19 simulations with ClayFF to determine the elastic properties of Na-montmorillonite as a function of water 20 uptake. (Ebrahimi, Pellenq, & Whittle, 2012) 21 22

  given by Meunier et al., K (x+y) [Al (4-y) Mg y ](Si (8-x) Al x O 20 )(OH) 4 , (Meunier & Velde, 2004) and chose x = 1 and y 21 = 0.5 to have an overall clay charge equal to -1.5e per unit cell. 22 23

7 8

 7 Figure 1. (a) Structure of illite whose chemical formula is K 1.5 [Al 3.5 Mg 0.5 ](Si 7 AlO 20 )OH 4 . Yellow, green, red, 9 white and grey spheres correspond to silicon, aluminum, oxygen, hydrogen and potassium, respectively. Illite 10 is made up of the stacking of tetrahedral (T) and octahedral (O) layers. About one (T) silicon position per unit 11 cell is filled with aluminum and aluminum in one (O) position of every second unit cell is substituted in an 12 isomorphous way with other metal atoms. (b) Schematic illustration of the different mechanical loadings 13

17 18Figure 2 3 2

 1723 Figure 2.(a) Illustration of the surface (dashed line) accessible for the blue probe sphere of radius r p in a 19 system of the red particles (b) A representative sample of the surface points (grey) obtained in a crack by a 20

22 23

 22 interatomic interactions are treated in ClayFF as non-bonded interactions which lends a large flexibility to the 20 force field. ClayFF was also successfully used to simulate transport in clays. (Botan, Marry, Rotenberg, Turq, 21 & Noetinger, 2013) (Wang, Kalinichev, & Kirkpatrick, 2006) Molecular dynamics simulations were carried out with the LAMMPS simulation program package at 300 K. 24 (Lammps: Molecular Dynamics Simulator, http://lammps.sandia.gov) The time step for the integration of the 25 molecular equations of motion was chosen equal to 0.1 fs in the case of ReaxFF and 1fs in the case of ClayFF. 1

23 1Figure 3 .

 233 Figure 3. Determination of the the stiffness parameters C ij from the curvature of the potential energy surface. 2

C 11 =

 11 C 22 , C 13 = C 23 , C 44 = C 55 and C 66 = (C 11 -C 12 )/2. As a result, a TI material has only five independent 13 stiffness constants out of 21: C 11 , C 12 , C 13 , C 33 , C 44 and all other elements are zero (note that similar symmetry 14

16 17Figure 4 .

 164 Figure 4.Stress σ (black empty circles, left axis) and crack surface area S (black filled circles, right axis) 18

6 Heinz, 2010) 7 8 Figure 6 . 11 12 19 -

 6861119 Figure 6.Stress-strain curve upon mode II loading of a crack parallel to the clay layers in illite. The observed 9

5 6 Figure 7 .

 67 Figure 7. Stress-strain curve computed with ClayFF, when a shear deformation (mode II loading) is applied in 7

Figure 8 .

 8 Figure 8.Stress σ (filled circles, left axis) and crack surface area S (empty circles, right axis) upon mode I 6

3 45

 3 AcknowledgementThis work was funded by Royal Dutch Shell and Schlumberger in the framework of X-Shale research project. 6

Table 2 .

 2 Stiffness tensor of illite in Voigt notation as obtained with ClayFF. Note that the tensor is shown 1 here as an upper triangular matrix due to symmetry reasons (C ij = C ji ). The shaded fields in the table show the 2 elements that are expected to be different from 0 for a transversely isotropic material (see text).

	3
	C ij /GPa	1	2	3	4	5	6
	1	292.5±0.5 128.3±0.4 16.67±0.08	1.0±0.3	5.5±0.6	-17.0±0.1
	2		274.2±0.4 15.97±0.04	0.1±0.7	2.0±0.3	1.7±0.1
	3			48.9±0.1	0.37±0.03 2.89±0.06 0.87±0.01
	4				8.99±0.02	-	0.8±0.2
						0.71±0.02
	5					9.82±0.04 -0.1±0.2
	6						90.1±0.2
	4

Table 3 .

 3 summarizes the VRH polycrystal properties of illite computed with 9 Elastic moduli of illite as obtained with the ReaxFF and ClayFF forcefields. 12

	10	ReaxFF and ClayFF.
	11
				ReaxFF	ClayFF
		Bulk modulus / GPa		52±6	73.66±0.08
		Shear modulus / GPa		18±2	35.17±0.03
		Speed of sound: v s / km/s		2.6±0.1	3.517±0.001
		Speed of sound: v p / km/s		5.4±0.3	6.511±0.003
		Young's moduli / GPa	1	110±51	225.3±0.7

Table 1 )

 1 , this comparison reveals that C 11 and C 22 are 17 indeed very close as they only differ by 2%. Considering the uncertainty over the results, C 13 and C 23 can also 18 be considered equal. However, C 55 is roughly twice larger than C 44 , while C 66 is 5 times smaller than (C 11 -19 C 12 )/2. In addition, C 34 , C 15 , C 25 , and C 35 are found to be significantly different from 0. The data obtained with 20 ClayFF (Table 2) displays a much better agreement with what is expected for a TI material: C 11 and C 22 are 21 not strictly equal but their difference is only about 6%. C 13 and C 23 are found to differ by about 4% while the 22 departure of C 44 from C 55 amounts to 8%. C 66 deviates from (C 11 -C 12 )/2 by 10%. The rest of the values (white 23 cells in Table 2) are all almost zero, except for the element C 16 .

  The 9 inverse of such a reference matrix S ref is determined based on the inverse of the computed stiffness matrices 10

	(i.e., the compliance matrix S): S 11	ref = (S 11 +S 22 )/2, S 33	ref = S 33 , S 44	ref = (S 44 +S 55 )/2, S 13	ref = (S 13 +S 23 )/2, and S 12

ref 11 = S 11 (S 12 /S 11 +S 12 /S 22 )/2. The rest of the elements of the reference matrix is set to 0. Such an analysis reveals 12

Table 4 .

 4 Then, & Lagaly, 2006) 6 resulting in a charge of -2e for muscovite instead of -1.5e for illite. Both clay minerals contain, however, K + 7 Comparison of the calculated fracture toughness K C for illite with experimental data for muscovite 10 Table 4 displays a consistent qualitative trend when comparing our calculated values to experimental data. 14

	8	counter ions only.
	9
	11	(Spray, 2010)
		Simulation (MPa•m 1/2 )	Experimental (MPa•m 1/2 ) a
			(for illite)	(for muscovite mica)
		Parallel to the 001 plane	Orthogonal to the 001 plane
		0.09 (ReaxFF)	0.61 (ReaxFF)	0.2
		0.12 (ClayFF)
	12	a from (Spray, 2010)
	13


Laurent Brochard
email: laurent.brochard@enpc.fr

Gyorgy Hantal
email: ghantal@mit.edu

Hadrien Laubie
email: hlaubie@mit.edu

Franz J Ulm
email: ulm@mit.edu

.-M Pellenq
email: pellenq@mit.edu

Fracture mechanisms in organic-rich shales: role of kerogen

In this work we study role of kerogen in the fracture properties of organic-rich shales, and in particular in the ductility of shales. The presence of kerogen and clays in shale is known to increase the ductility. We propose here a multiscale approach to develop of a fine understanding of shale ductility from the molecular scale. We develop and validate a methodology at the molecular scale that can capture the toughness and ductility of a material. We apply this methodology successfully to a silica polymorph and to a kerogen analog, and we confirm the significant ductility of kerogen. Interestingly the silica-kerogen interface exhibits a similar ductility, which is central for the properties of the heterogeneous shale. Finally, we consider a tentative upscaling considering the pull out phenomenon as a likely mechanism of fracture of the shale.

INTRODUCTION

Organic-rich shales are larges reserves of unconventional hydrocarbons that have become more and more attractive with the increasing price of conventional gas and oil. The potential resources of shale gas are estimated between 50 and 100 years of current world consumption [1], but the recovery of these unconventional hydrocarbons is challenging because of the low intrinsic permeability of those reservoirs. To operate such reservoirs at reasonable costs, the oil industry performs hydraulic fracturing. Hydraulic fracturing consists in creating a network of fractures in the reservoir by injecting a fluid under high pressure [2]. The most important factors that influence the crack propagation in the reservoir are the in situ far field and local stress and the fracture properties of the rocks. A fine understanding of the fracture properties of shales is critical for the reliable modeling of hydraulic fracturing. In this work, we focus on a specific aspect of the fracture properties: how the fracture properties depend on the mineralogy and organic content of the shale? Brittleness indicators have been developed to quantify the level of ductility of the shale [3]: shales containing mostly silica and carbonate tend to be very brittle, i.e., they "shatter", leading to a vast array of small-scale induced fractures providing numerous flow paths. Conversely, shales with a significant content of clay and kerogen tend to be ductile, i.e., they deform instead of shattering, leading to relatively few induced fractures and providing only limited flow paths. The brittleness of shale is important to decide the appropriate location of the horizontal well and to estimate the productivity of the well.

In this work, we aim at developing a fundamental understanding of the toughness and ductility of shales starting at the molecular scale. The characteristic size of the elementary constituents of shales (minerals and kerogen) is as small as a micrometer. At this scale, the study of fracture properties by laboratory experiments is very challenging. Whereas molecular simulation is a convenient alternative to experiment to study the fracture properties at the nano-scale. Ultimately such an approach would enable to capture finely the impact of organic / clay content, or kerogen type / maturity on the fracture properties of shales. In the first section, we present the methodology we develop to estimate fracture properties by molecular simulation. In the second section, we apply this methodology to silica, kerogen and their interface. In the third section, we propose to upscale those fracture properties to estimate the fracture properties of the heterogeneous material.

METHODOLOGY

We propose a methodology to estimate fracture properties by molecular simulation [4]. The methodology, based on the energetic theory of fracture mechanics [5,6,7], consists in a thermodynamic integration. We consider a bulk material in which we initiate an initial elliptic through crack (see 1). Then, we load the system in the direction orthogonal to the crack until the crack is fully propagated. During the loading we impose the displacement at the top and bottom boundaries, we impose a constant temperature and we compute the derivative Σ = 1 LxLz ∂F ∂Ly T,Lx,Lz of the Helmholtz free energy F with respect to the vertical direction L y , which can be interpreted as the average stress at the top and bottom boundaries. To compute Σ, we use the virial equation [8]. This procedures corresponds to a loading of the crack in mode I, but can be easily adapted to a loading in mode II by applying a shear instead of a tension. According to the energetic theory of fracture mechanics, the Helmholtz free energy F can take the form of mechanical energy P due to the strain (r) and stress σ (r) over the volume Ω of the solid, i.e., P = Ω (r) 0 σ : d dr, and can be released during crack propagation, the critical energy release rate G c being the corresponding energy released per unit area of crack propagated. Assuming that no other form of energy transformation occur during the loading of the system, the energy balance of the system is: dF = dP +G c dA, where A is the crack area.

At the beginning of the procedure, the system is free of stress and the mechanical energy P is zero. Once we start to increase the loading, the system stores mechanical energy but the crack does not propagate due to insufficient stress intensity at the crack tip. When the loading is large enough, the crack starts to propagate and the mechanical energy stored is released during the propagation. At the end of the procedure the crack is fully propagated over all the system and the mechanical energy stored is again zero since the stress in the system is free of stress. Therefore integrating the curve Σ (L y ) over the whole process provides an estimation of the critical energy release rate G c :

G c = L x L z ∆A Lymax L y0 ΣdL y (1)
where ∆A is the total area of crack created. The thermodynamic integration we propose is not limited to linear elastic materials but is valid whatever the mechanical behavior and whatever the fracture properties (brittle or ductile), under two important assumptions: 1 -the crack propagation is assumed in isothermal conditions, 2 -the irreversible processes that occur during the loading are related to the crack propagation and no other form of dissipation take place, e.g., plasticity not related with the crack propagation. The first assumption, isothermal conditions, is valid if the heat conductivity of the material is large enough. For, the second assumption, cracks with large aspect ratios ensure a sufficient stress singularity at the crack tip and we did not observe any other form of dissipation in all the systems we simulated. In addition to those two assumptions, special care has to be taken on a few aspects [4]: size of the system, the inter-atomic energy potential, and overlap at the periodic boundary conditions.

As an alternative to the thermodynamic integration, one can also consider the usual Linear Elastic Fracture Mechanics (LEFM) to derive the critical stress intensity K Ic , also called toughness, from the maximum average stress Σ cr at the onset of crack propagation 1. According to LEFM, the toughness and the critical stress are linearly related [7]:

K Ic = Σ cr √ πaC = Σ cr √ πa ef f (2)
where a ef f = aC 2 and the coefficient C is a parameter that depend on the geometry, e.g., C equals 1 for a single crack in an infinite medium. For the doubly periodic geometry we consider, C is a function of a/W and H/W , where the geometrical parameters a, W , and H are defined in Figure 1. We used the pseudo-traction method proposed by Karihaloo et al. [9,10] to compute numerically the coefficient C (a/W, H/W ) for all the systems we considered.

Using Equation 2we can derive the toughness K Ic from the critical stress Σ cr , and the toughness is related to the critical energy release rate G c according to Irwin's formula. For an isotropic material in plane strain [6], which is the case of kerogen in the next section, we have

G c = 1-ν 2 E (K Ic ) 2
, where E is the Young's modulus and ν is the Poisson's ratio. For an anisotropic material, which is the case of the silica we consider in the next section, the Irwin's formula is more complex and depends on the different components of the elasticity tensor (see [11]).

For a linear elastic brittle material, the two approaches, i.e., thermodynamic integration and LEFM, are equivalent. However, when a material exhibits a significant ductility, e.g. due to a large plastic zone around the crack tip, the results of the two approaches differ: LEFM is not valid and underestimates the toughness whereas the thermodynamic integration is valid and lead to the actual critical energy release rate.

APPLICATION TO SILICA, KEROGEN AND THEIR INTERFACE

In this section, we apply the methodology presented in the previous section to two of the most common materials in shale: silica and kerogen. In addition, we apply the methodology to molecular reconstructions of the silica-kerogen interface. Silica is a good benchmark for the approach we propose because its mechanical behavior is linear elastic and its fracture behavior is brittle. Therefore, both LEFM and the thermodynamic integration are valid for this material. In addition, many experimental results are available for silica. In contrast, kerogen is known to increase the ductility of shale and, accordingly, we expect to identify this ductile behavior by comparing the results of LEFM and thermodynamic integration.

We considered α-cristobalite [12], a silica crystal metastable at ambient temperature and pressure. We use this specific silica, because the tetrahedral lattice of αcristobalite adapts very well to the lattice of the molecular representation of kerogen we use, which is key for the simulation of the silica-kerogen interface. For kerogen, we consider a molecular representation of microporous carbon, CS1000, representative of the matrix of a mature kerogen without mesopores [13]. The reactive energy potential we use is reaxFF [14] adapted to systems made of silicon, oxygen, carbon and hydrogen atoms [15,16].

For both α-cristobalite and CS1000, we performed the loading procedure over many systems with various sizes and various crack lengths, and varying the orientation of the crack in the case of α-cristobalite which is anisotropic. We simulated a total of 33 systems for α-cristobalite and 25 systems for CS1000 [4]. We display in Figure 2 the resulting L y -Σ curve for three of those systems, one for α-cristobalite with the crack in the (001) plane, a second one for α-cristobalite with the crack in the (100) plane (equivalent to the (010) plane by symmetry) and one for CS1000. The two curves for alpha-cristobalite cannot be compared quantitatively since the initiation of the crack is dependent on the geometry which is not the same between those two systems. But in both case the behavior is linear and brittle with a sudden propagation of the crack and a sharp drop of the average stress. By contrast, the behavior of CS1000 is much more ductile with a progressive decrease of the average stress during crack propagation.

We performed the thermodynamic integration for all the systems and we computed the critical energy release rates of CS1000 and α-cristobalite for which we distinguished between the crack propagation in the (001) plane and the propagation in the (100) plane. The results for α-cristobalite were almost independent on the crack orien- tation, crack length and size of the system (Fig. 3). However, the results for CS1000 exhibited a clear correlation with the area of crack created (Fig. 3). We attributed this trend to the overlap of the process zone at the periodic boundary. Indeed, the process zone at the crack tip is important in the case CS1000 whereas it is non-significant for α-cristobalite. A tentative estimate of the size of the process zone is the size r pl of the "plastic region" around the crack tip assuming a plastic behavior of the material. We consider the Dugdale-Barenblatt estimate: r pl = (π/8) × (K I /σ pl ), where σ pl is the plastic yield stress of the material. We estimated σ pl = 35 GPa for α-cristobalite and σ pl = 19 GPa for CS1000. Therefore we have r pl = 3.38 Å for α-cristobalite and r pl = 19.9 Å for CS1000. This size of the process zone is very small for α-cristobalite compared to the typical size of the system we considered (≈ 70 Å) and, accordingly, it is reasonable to neglect any overlap of the process zone at the periodic boundary. However, for CS1000 the overlap cannot be neglected. As a first order correction for this overlap we considered an effective crack area created ∆A ef f = ∆A -L z r pl /2 instead of the actual crack area ∆A in the integration Eq. 1. In Figure 3 we display the corrected results and no trend can be identified. According to LEFM, the critical stress Σ c r is inversely proportional to √ πa ef f (see Eq. 2). This relation is well verified for α-cristobalite but less accurate for CS1000 (see Fig. 4). We fitted the value of the toughness K Ic (lines in Figure 4) and we compare those values to that obtained previously by thermodynamic integration in Table 1. In Table 1 we converted the critical energy release rates into toughness using the Irwin's Formula for CS1000 (isotropic) and its derivative for α-cristobalite (orthotropic) [11].

We simulated bulk CS1000 to estimate its Young's modulus and Poisson's ratio. For α-cristobalite, we used the elastic tensor measured experimentally [17]. In addition, we include in Table 1 the average toughness measured experimentally for silica [18].

Effective half length of the crack πa ef f , Å Critical stress Σcr, GPa For α-cristobalite, the results obtained by thermodynamic integration, derived from LEFM and measured experimentally are very consistent. Accordingly, the methodology we propose works well for this example which is a linear elastic brittle material. For CS1000, however, the results of thermodynamic integration are significantly higher than the results derived from LEFM. This difference is due to the important ductility of this microporous carbon that is readily visible in Figure 2. LEFM is not valid when the material experiences a significant ductility at the scale of the system and leads to an underestimation of the toughness. The thermodynamic integration however is valid even when the material exhibits a significant ductility. Therefore, the value obtained by thermodynamic integration can be regarded as the actual toughness of CS1000. The toughness of CS1000 is much higher than that of silica, but this conclusion is specific to CS1000 and cannot be generalized to kerogen. Indeed, CS1000 is representative of a high density mature kerogen without its mesopores, and accordingly its toughness is in the upper range of what can be expected for kerogen. The ductility of CS1000, however, is an intrinsic property of the microporous carbon which is transferable to kerogen. This is consistent with macroscopic observations: laboratory experiments available in the literature have shown that a significant organic content in shale increases the ductility of shale, but not necessarily the toughness [19,20].

The fracture properties of a heterogeneous material depends on the fracture properties of the its elementary constituents but also on the fracture properties of the interfaces between those constituents. The crack propagation often occurs at interfaces which have a weaker resistance to fracture. In a previous work [21], we developed a methodology to reconstruct molecular models of silica-kerogen interfaces and we applied this methodology to generate a series of 11 interfaces between CS1000 and α-cristobalite assuming an average situation in terms of kerogen type and maturity. We then applied the mode I loading to these interfaces (see Fig. 5) and estimated the critical energy release rate by thermodynamic integration. For a crack lying at an interface, the theoretical stress field near the crack tip is oscillatory because of the mismatch between the elastic properties of the two materials. As a consequence, the toughness is a combination of the mode I and mode II stress intensities and to convert the critical energy release rate into toughness, one has to consider an alternative formulation of the Irwin's formula [6,22]. Thus we could estimate the toughness of the interface (see Tab. 1). We also derived the toughness from LEFM as we did for pure CS1000 and pure α-cristobalite. In theory, this approach by LEFM is derived for homogeneous materials and is not strictly exact for interface cracks. But we used this approach to have a basis of comparison with the results from thermodynamic integration and therefore quantify the ductility of the interface. The strain-stress curve for the interface system (Fig. 5) exhibits clearly some ductility like CS1000 (Fig. 2). This conclusion is supported by the discrepancy between the toughness estimated by thermodynamic integration and that derived from LEFM (Tab. 1). The critical energy release rate of a perfectly brittle material is equal to two times the surface energy γ. Therefore, a possible quantification of the ductility d is the ratio between the critical energy release rate and the surface energy: d = G c /γ. d = 1 corresponds to a perfectly brittle material, and d → +∞ corresponds to a perfectly ductile material. Here, for sake of simplicity we assume that the energy release rate derive from LEFM corresponds to the perfectly brittle material:

G LEF M c = 2γ. Therefore, we have d = G thermo c /G LEF M c = K thermo Ic /K LEF M Ic 2
. With this definition of ductility, we have d = 1.40 and d = 1.61 for α-cristobalite with the crack orientations ( 001) and (100), respectively; we have d = 4.07 for CS1000; and d = 4.75 for the interface. According to this tentative quantification of ductility, the interface is as ductile as the CS1000 itself.

EFFECTIVE FRACTURE PROPERTIES OF HETEROGENEOUS SHALE

In the previous section, we have shown that kerogen exhibits a significant ductility and therefore can turn a brittle shale into a ductile material. However, how the fracture properties of the different constituents of shale combines in the heterogeneous medium needs to be investigated. In particular, the kerogen content is usually small (≈ 5%) and therefore it is not straightforward that such a small content can significantly affect the ductility of the shale. The presence of clay also increases the ductility of shale. We propose here an approach to upscale the fracture properties of shale assuming the shale as a heterogenous material made of a brittle matrix with ductile inclusions. We apply this approach to the specific case of kerogen inclusion in a silica matrix with the results of the previous section. But this upscaling could be adapted to account for clay as the ductile inclusions or for carbonate as the brittle matrix, as long as the clay content remains small (< 20%).

The crack propagation in a heterogeneous material is subject to various unusual phenomena such as crack deflections [23,24,25,26], front roughening [27], trapping [28,29], bridging [30,31,32,33], particle pull out [30,31,34], micro-crack toughening [30,35] or crack pinning [36,37]. All these phenomena combine and should be accounted for to fully capture the effective properties of the shale. However, a model combining all these phenomena would be far too complex and we decided to focus on a few of those mechanisms that we consider the most relevant in the present situation. When a crack lying in the brittle matrix (silica) approaches a ductile inclusion (kerogen), the crack front tends to deform to account for the heterogeneity of toughness and mechanical properties between the matrix and the inclusion, and then the crack either penetrates the inclusion or is deflected at the interface. According to the works of He and Hutchinson [25], in the specific case of kerogen and silica, the crack is always deflected. Indeed, the ratio between the critical energy release rates G Interf ace c and G CS1000 c at the interface and in CS1000, respectively, is 0.143. According to the work of He and Hutchinson, whatever the angle of approach of the crack at the interface, this ratio is smaller than the ratio of energy release rates between the two scenarios deflection and penetration (at least 0.25). Therefore, following the criterion of He and Hutchinson, the crack never penetrates the ductile inclusions but bypasses eventually forming bridges between the two crack faces until the ductile particle is pulled out inside the crack space. Therefore, we focused our investigation on the particle pull out mechanism which is known to significantly affect the fracture properties of heterogeneous materials. We neglected the front roughening and crack trapping since the impact of these phenomena on the fracture properties is second order compared to the effect of the particle pull-out.

We considered the model of Bower et al. [31] which assumes a periodic distribution of circular bridges in the plane of the crack (Fig. 6). Each inclusion exerts a total force P on the crack faces which depends on the crack opening δ at the location of the inclusion. The model enables to compute the toughness ratio K ef f Ic /K matrix Ic of the heterogeneous material with respect to the matrix, in function of the pull-out law P = P (δ) and of the geometry of the system (R/L, H/L, x 0 /L). The key input of this model is the pull-out law P = P (δ). Bower et al. studied the case of a friction law which makes sense for the case of fiber reinforced materials. In the present case, we considered an other pull-out law (see Fig. 6) inspired from Figure 5, i.e., representative of the ductility at the kerogen-silica interface. The values of the slope K p.o. and of the critical crack opening δ p.o can be related to the total energy release during the complete pull-out of the particle. Assuming a spherical inclusion which is pulledout on half of its surface, we have: To fully determine the pull-out law a second condition is necessary. As a starting point, we considered an elasticity of the particle pull-out equal to the elasticity of the matrix:

K ref p.o. = (2πRE) / (2 (1 -ν 2 ))
, where E and ν are the Young's modulus and Poisson's ratio of the matrix. In the following, we vary the values of K interf ace . For these calculations, we assumed H/L = 1 and x 0 /L = 0.5. The toughness ratio K ef f Ic /K matrix Ic increases with the fraction of inclusion, i.e., the inclusions are toughening the material. The toughening can be very significant, up to a factor of three. This toughening is very sensitive to the interface toughness and to the pull out elasticity. In the case of CS1000 and α-cristobalite K interf ace Ic /K matrix Ic ≈ 0.7 and K p.o. /K ref p.o. ≈ 1. But CS1000 has a higher elasticity than usually admitted for kerogen because the mesopores are not present in this molecular model. The values of elasticity reported in the literature for kerogen correspond more to the curve K p.o. /K ref p.o. ≈ 0.125. The toughening under those conditions is very significant, more than 2 for a fraction of inclusion of 10% and more than 3 for a fraction of 20%.

As we did in the previous section we can quantify the effective ductility of the heterogeneous material as the ratio between the effective critical energy release rate and the effective surface energy: d = G ef f c /γ ef f . Considering the LEFM estimate of the previous section as estimate of the interface energy, we can compute the ductility d. We obtain a ductility of 1.61 for a surface fraction of ductile inclusions of 0% (pure silica), 9.97 for a surface fraction of 10% and 22.6 for a surface fraction of 20%. Therefore the presence of ductile inclusions in shale may increase dramatically the ductility of the rock.

CONCLUSION

In this work, we present a comprehensive approach to study the ductility of organic-rich shales starting at the molecular level and upscaling to the scale of the heterogeneous medium. We develop a methodology to capture the toughness and the ductility by molecular simulation, and apply it to a series of molecular systems representing a silica polymorph, a typical microporous carbon, representative of kerogen, and their interface. We validate the approach by comparing our results with experimental results and with predictions from LEFM valid for linear elastic brittle materials like silica. The kerogen analog and its interface with silica exhibit a significant ductility. Considering an upscaling scenario based on the pull out phenomenon, we perform a tentative upscaling to understand how the presence of ductile inclusions in shale can increase the ductility of the heterogeneous material. We apply this upscaling approach to the case of silica and kerogen and show that the ductility of shale can increase tremendously.

Figure 1 .

 1 Figure 1. Schematic representation of the proposed methodology.

Figure 2 .

 2 Figure 2. Example of loading curves Σ = Σ (L y ) for α-cristobalite (left) and CS1000 (right).

Figure 3 .

 3 Figure 3. Toughness estimated from thermodynamic integration for α-cristobalite (left) and CS1000 (right).

Figure 4 .

 4 Figure 4. Fit of the critical stress by LEFM for α-cristobalite (left) and CS1000 (right).

Figure 5 .

 5 Figure 5. Example of loading curve Σ = Σ (L y ) for interface

P

  (δ) dδ = K p.o. (δ p.o ) 2 . Therefore the values of K p.o. and δ p.o are constrained by the ratio K interf ace Ic /K matrix Ic .

Ic

Figure 6 .

 6 Figure 6. Model considered to represent the pull out phenomenon. We display in Figure 7 the toughness ratio K ef f Ic /K matrix Ic in function of the surface fraction of the inclusions and for various values of K interf ace Ic

Figure 7 .

 7 Figure 7. Toughening ratio in function of the various parameters of the problem.

Table 1 .

 1 Summary of the toughnesses estimated by thermodynamic integration, derived from LEFM and measured experimentally

	Estimate / measure	K Ic , MPa.m 1/2 Standard deviation, MPa.m 1/2
	α-cristobalite
	Thermodynamic integration, (001) orientation	1.027	0.053
	Thermodynamic integration, (100) orientation	1.027	0.060
	LEFM, (001) orientation	0.867	0.056
	LEFM, (100) orientation	0.810	0.074
	Experimental measurements	0.818	0.070
	CS1000
	Thermodynamic integration	1.608	0.127
	LEFM	0.797	0.099
	Interface between CS1000 and α-cristobalite
	Thermodynamic integration	0.573	0.036
	LEFM	0.263	0.029


György Hantal
email: ghantal@mit.edu

Laurent Brochard

Benoit Coasne

Franz-Josef Ulm

J.-M Pellenq

CHEMISTRY OF INTERFACES BETWEEN INORGANIC MINERALS AND POROUS CARBONS: IMPLICATIONS FOR THE MECHANICAL PROPERTIES OF GAS SHALE

Introduction

The role organic-rich shale plays in the production of fossil fuel has tremendously increased in the past decade thanks to recent technological developments that allow an efficient exploitation of this unconventional source rock. Gas shale is composed mostly of minerals but also a carbonaceous substance, called kerogen which is primarily responsible for hydrocarbon storage. As kerogen is encapsulated in a nonpermeable mineral matrix, gas and oil recovery operates by means of hydraulic fracturing. One of the many challenges this technique is facing is controlling fracture propagation underground. Understanding how different gas shale components influence the fracture properties of the overall rock is thus of paramount importance for efficient well design. Yet, probing mechanical properties of different gas shale constituents is extremely difficult experimentally due to the very small size of these heterogeneities.

In this study, we make use of an alternative approach, based on molecular simulations, to elucidate fracture properties. We apply our methodology to illite, which is the most prominent clay constituent in gas shale, as well as to a composite of illite and a disordered porous carbon modelling kerogen.

Methodology

The goal of this work is to simulate fracture propagation in our molecular systems and obtain fracture properties. Our approach is based on the energetic description of fracture mechanics. [START_REF] Leblond | Mécanique de la rupture fragile et ductile: Hermès Science[END_REF] Unlike the description based on the evaluation of stress intensity factors at the crack tip, this approach provides an unequivocal description of the system as all quantities involved can be calculated in molecular simulations. [START_REF] Brochard | Fracture toughness from molecular simulation: Application to bulk silica, microporous carbons and their interface[END_REF] In the simulations, a mechanical load is applied to the system (in our case by prescribing the deformation) in every step, and the mechanical answer is measured (the average stress at the boundaries). Upon fracture propagation, the mechanical energy is released by creation of crack surface. The fracture energy, that characterizes the toughness of the system, is obtained by integrating the stress-strain curve of the process. Here we consider two different loading modes: in mode I, a tensile strain is applied normal to the crack orientation, whereas in mode II, a shear strain acts parallel to the crack. In both the pure illite and the composite systems, our focus is mainly on deformations on the 001 surface of illite, as clays are known to cleave primarily along this plane. To perform the molecular simulations, we chose to apply the ReaxFF reactive force field that is able to describe both intermolecular interactions and bond breaking and formation in an effective empirical way. [START_REF] Van Duin | [END_REF]4] To simulate the organic phase in the composite system, we use the CS1000 molecular structure, which is a representative model of porous disordered carbon. [5] Bond formation was stimulated at the organic/inorganic interface by applying an annealing procedure starting form 600K and cooling the system down to 5 K with steps of 50 K, simulating every step for 2 ps.

Results and Discussion 1) Fracture of pure illite Illite was found to exhibit markedly different fracture behaviors depending on the loading mode. In mode I, the material displays a very brittle behavior and a very small fracture energy: it is two orders of magnitude smaller than what was found for silica in our previous study. [START_REF] Brochard | Fracture toughness from molecular simulation: Application to bulk silica, microporous carbons and their interface[END_REF] The fracture energy was found to be twice the surface energy (known as the Griffith criterion), which confirms the very brittle character of illite in this fracture mode.

In mode II, however, no fracture propagation is observed, that is, there is no fracture surface creation, and the mechanical energy is dissipated by the friction induced by clay layers sliding on each other. Technically, this means that the fracture energy (as well as the toughness) is infinity. These two different mechanisms were confirmed by evaluating the crack surface area upon loading by Conolly's method. Figure 1 displays typical stress-strain curves obtained in different loading modes of illite. In mode I, illite exhibits a highly brittle behavior. In mode II, no fracture propagation is observed. Each stress drop in the stress-strain curve of mode II corresponds to one sliding event in the system.

The distinctly different behaviors can be understood by considering the lack of chemical bonds between clay layers where the cohesion is only provided by week electrostatic interactions. This means that none of the releasing mechanisms (fracture propagation or friction) involve chemical bond breaking. This is confirmed by the fact that no initial crack creation is needed in any loading modes. Based of these facts, it is possible to use a non-reactive potential to simulate the different mechanical loadings in the system. To this end, we repeated some simulations with the ClayFF potential. [6] The calculations showed the same behaviors and qualitatively similar fracture energy values in mode I loading. 2) Fracture of illite/kerogen composites In the first step, bond formation at the kerogen/illite interface was studied as resulted from the annealing procedure. To this end, radial distribution functions (RDFs) were calculated to check the existence of possible bonds. A combined analysis of RDFs revealed bond formation between illite and CS1000 via Si-O-C and Al-O-C bonds, but ruled out the possibility of direct chemical bonds between surface Al or Si atoms of illite and C atoms of CS1000. As a consequence of bond formations between the two phases, the surface of illite gets reconstructed and the surface siloxane bonds get reorganized forming lowmembered rings.

In the second part, fracture analysis of the prepared kerogen/illite interface was carried out. In one set of the calculations the system was loaded without any initiated crack, whereas in another set of calculations a crack was created in the interface.

The results showed that fracture propagation can occur both at the interface and the clay phase depending on the relative strength of the interface compared to illite: if the interface is weaker than illite, the fracture will occur in the interface. Figure 2 shows a typical stress-strain curve in the crack free system. The comparison of this figure to Figure 1 reveals that even though the fracture propagation occurs in illite, the material shows significant ductility as well as a much larger surface energy (5 times larger than in the pure case). This observation provides molecular structure based explanation for the apparent toughening of organic rich shale compared to the organic-free case.

Figure 2.

Stress-strain curve of kerogen/illite composite systems without an initiated crack in the system. Even though the fracture propagates in the clay phase, the system displays a much more ductile behavior than pure illite as well as a much larger fracture energy. This is related to the fact that mechanical energy is dissipated in the ductile organic phase upon loading the system. This mechanism can provide explanation for the observed toughening of organicrich shales.

Conclusions

In this work the fracture properties of pure illite and illite/kerogen composites were studied by means of molecular simulations. For pure illite, markedly different fracture behaviors were observed depending on the loading mode. In mode I, illite shows a very brittle behavior with very low toughness, whereas in mode II, no fracture propagation occurs, but the mechanical energy is dissipated by friction induced by clay layers sliding on each other. Technically this means a super plastic behavior with infinite fracture toughness.

In illite/kerogen composites, first the possibility of chemical bond formation induced by an annealing procedure was explored It was found that chemical bond formation can occur between illite and the CS1000 structure used to model kerogen. A detailed fracture study of this composite system showed that fracture propagation can occur either in the kerogen phase or in the interface depending on their relative strength. In cases where the interface is stronger, hence fracture occurs in illite, a very ductile behavior was observed in mode I, as opposed to the very brittle behavior of pure illite. The observed fracture energy is 5 times higher compared to the pure case.

Both the friction based dissipative mechanism of the mode II loading of pure illite and the toughening of illite/kerogen composite systems due to the ductile organic phase can provide explanations for the apparent increased fracability of argillaceous and organic-rich shales.

Figure 1 .

 1 Figure 1. Stress-strain curves obtained from molecular fracture simulations of illite in (a) mode I and (b) mode II.In mode I, illite exhibits a highly brittle behavior. In mode II, no fracture propagation is observed. Each stress drop in the stress-strain curve of mode II corresponds to one sliding event in the system.