György Hantal

Laurent Brochard

M Natália

D S Cordeiro

Franz J Ulm

J.-M Pel- Lenq

Surface chemistry and atomic-scale reconstruction of kerogen-silica composites

The interest in gas shale, a novel source rock of natural gas, has tremendously increased in recent years. Better understanding of the kerogen/rock interaction is of crucial importance for efficient gas extraction and hence for asset management. In this study, we explore the possible chemical bonds between kerogen and silica, one of the most predominant mineral constituents of gas shale by means of quantum chemistry. Energetically favorable bond formation reactions are found between alcoholic hydroxyl, carboxylate and aldehyde groups as well as aliphatic double bonds of kerogen and the silica surface. The performance of a reactive force field is also assessed in a representative set of chemical reactions, which was found satisfactory. The potential impact of bond formation reactions between the two phases on the actual kerogen-silica interface are discussed in function of the kerogen type, maturity and density. Finally, a methodology aiming to reconstruct realistic kerogen-silica interfaces is presented.

Introduction

Gas and oil shale are unconventional deposits of hydrocarbons that represent significant reserves worldwide. Yet, recovering gas and oil from organic-rich shale is much more challenging than from conventional resources. Thanks to recent advances in the recovery techniques, the past 10 years have seen a massive increase in shale gas production: in 2010, the share of shale gas has reached 27% of the total gas production of the U.S. [START_REF] Bonakdarpour | The Economic and Employment Contributions of Shale Gas in the United States[END_REF] Gas shale can be considered a natural composite material as it consists of inorganic minerals as well as two main organic phases, kerogen and bitumen, encapsulated in the mineral matrix. [START_REF] Tissot | Petroleum Formation and Occurrence[END_REF][START_REF] Vandenbroucke | [END_REF] Kerogen, which is by far the major organic component in organic-rich shale, represents only a small fraction, typically around 5 weight percent. However, this organic matter hosts most of the hydrocarbons and therefore plays a central role in the production process. Due to the very low permeability of oil and gas shale rocks, hydraulic fracturing is generally used to propagate cracks in the shale reservoir and thus access those resources. 4 Kerogen, is an amorphous, carbonaceous material with complex hierarchical porosity, which is suspected to have significant effect on the overall rock properties. [START_REF] Vandenbroucke | [END_REF]4 Indeed, laboratory experiments comparing organic-rich shale and normal shale without organic matter have shown that the mechanical behavior during crack propagation is significantly influenced by the presence of organic matter. 5,6 Understanding how kerogen is bonded to the mineral surface and how this bonding depends on the kerogen type, composition, density or maturity is critical for correctly describing the fracture mechanics of oil and gas shale reservoirs and ultimately improving the design of the hydraulic fracturing process. The inorganic matrix in gas shale shows diverse mineralogical compo-sition depending on its sedimentary origin and depositional conditions. The most important inorganic constituents are calcite, clay and silica. In this paper, we focus on silica as it is the major mineral constituent in many shale formations. [START_REF] Vandenbroucke | [END_REF] Despite the large potential in the exploitation of gas shale reserves, the intrinsic complexity of the material has prohibited understanding the nature of cohesion between gas shale components. The present study is the first attempt to describe the chemistry at the kerogen-silica interface in gas shale. Our ultimate goal is to construct a molecular model of such an interface with particular emphasis on the right chemistry. As no experimental information is available regarding bonding at this interface, in the first part of this study we investigate, by means of quantum chemical calculations, bond formation between the silica surface and single organic molecules bearing functional groups typical in kerogen.

In the second part, based on the results from the first part, we investigate the influence of kerogen type, maturity, and density on the repartition of the different chemical bonds at the interface assuming an equilibrium bonding structure. We then propose a realistic molecular reconstruction technique, finally we also present one actual reconstructed model of a kerogensilica interface. Such reconstructions are the basis for a systematic study of the fracture properties at the organicinorganic interface. [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF] The strength of this approach is that it investigates the mechanical behavior of the composite organic-inorganic system based on a comprehensive combination of available experimental data and accurate quantum chemical assessment of the interface chemistry. Although a large body of research has already focused on physisorption of organic molecules on silica surfaces (among others see Refs 8-11), this study is the first, to the best of our knowledge, to characterize chemisorption of organic matter on silica surface.

The paper is organized as follows. In section 2, details of all calculations will be presented. Section 2.1 is dedicated to the studied organic and inorganic interfaces as well as the summary of the applied strategy to examine possible bonds between the organic and inorganic phases. Section 2.2 presents the investigated chemical reactions and the details of the quantum chemical (QM) calculations. Section 2.3 shows the assessment of the performance of an empirical reactive force field with respect to QM calculations to ensure the former is adequate for studying fracture properties. [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF] Section 2.4 introduces the basis of the reconstruction technique. Section 3 discusses the most probable chemical bonds at the interface as well as the implications of these bonds on eight possible kerogen-silica surfaces. Finally, an actual kerogen-silica interface is also reconstructed.

Materials and methods

In this section we review the main steps of our hierarchical approach, as well as we provide the details of our calculations.

-Bond formation energies of reactions between the silica surface and single organic molecules are studied at high quantum chemical level. We consider the hydrated 001 surface of alpha-quartz and monofunctional organic molecules with functional groups characteristic of kerogen. For the same functional group, calculations are repeated with increasing side chain lengths to extrapolate the bond formation energy for large organic compounds (i.e. kerogen).

-We propose a general reconstruction strategy based on results from accurate quantum chemical calculations and available structural information. We discuss the distribution of functional groups at the interface of alpha-quartz and the organic phase in the case of 8 different kerogen compositions corresponding to different types, maturities and densities. Based on the Boltzmann distribution, we calculate probabilities of bond formation between the two phases using the above computed reaction energies.

-We generate an interface between an actual porous disordered carbon model (representing kerogen) and a silica phase based on the information gained from the two previous steps. We use a reactive molecular force field for the reconstruction due to the large number of atoms in the system.

The studied surfaces

The inorganic surface. In the first step, we have to determine what typical functional groups can be expected at the surfaces of the organic and inorganic phases, representative of organic-rich shale. The most stable form of silica at ambient conditions as well as in the upper part of Earth's crust is alphaquartz. The cleavage of silica results in surfaces completely covered with silanol groups as the freshly cleaved pristine silica reacts promptly, if present, with water. As water is common in the subsurface environment, silica can reasonably be expected to be present in its hydrated form. The most stable hydrated crystallographic surface of alpha-quartz is the (001) surface, thus we consider this surface in the first part of our study. 12,13 Note that different crystallographic surfaces of other silica phases differ from the studied one only in the concentration of surface silanol groups. Therefore, silica surfaces are expected to be very similar chemically. 14 The starting point of this study was a structure optimized by Goumans and coworkers. 15 In the first step, a cluster of ca. 130 atoms was cut out near the interface in such a way that only fully coordinated Si atoms were considered. All dangling bonds were saturated with H atoms. This cluster (see Figures 1a and1b) was then used to represent the silica surface throughout the study. The size of this cluster was chosen as a compromise between computational cost and a large enough surface ensuring realistic thickness and width around the chemically bonded organic molecule. To ensure that the cluster retains a structure similar to what it would have as part of an infinite system, some atoms were kept at a constant position in the calculations. To determine the active zone of the cluster, i.e. the group of atoms that are allowed to move in the structure optimization calculations, we applied the following criteria: i) Atoms in OH groups are allowed to move if they are on the top of the cluster (i.e. original surface OH groups, see Figure 1) and they have both H-donor and H-acceptor neighbors in the surface H-bonding system. ii) Silicon atoms directly connected to moving OH groups can move. iii) Every atom that is connected to a moving atom and is at least 3 atoms away from the boundaries of the cluster is also allowed to move. Atoms of the active zone are marked with blue spheres in Figure 1b. In order to probe system size effects, some calculations were repeated on a much larger cluster of 276 atoms (see Figures 1c and1d). The results showed that the obtained energies differ only to a small extent on the two clusters (less than 2%) which is, as we will see, within the error bars of our calculations. This confirms that our smaller cluster is appropriate for this study.

The organic surface. As already discussed, kerogen is a highly heterogeneous substance structurally as well as chemically. [START_REF] Tissot | Petroleum Formation and Occurrence[END_REF][START_REF] Vandenbroucke | [END_REF] Kerogen can be seen as a macromolecule consisting of carbonaceous nuclei cross linked with chain-like bridges. The bridge as well as the nuclei can bear functional groups. These building elements lend a complex 3D-structure to kerogen with nano-and mesoporosity. Elemental analysis of kerogen shows that by far the main atomic constituents are carbon and Oxygen is next in importance with amounts up to 25 mol% in very immature sediments. Sulfur and nitrogen content is secondary, the latter being most probably in a chemically inactive form (heterocycles). [START_REF] Tissot | Petroleum Formation and Occurrence[END_REF] Upon maturation the relative carbon content increases and, depending on the origin of the shale as well as the burial depth, gas or oil or often both are produced. Maturation also promotes the carbonization of kerogen resulting in increasing aromatic content, which is often organized into stacks constituting the nuclei. Functional groups are present inside as well as on the surface of kerogen. Based on infrared spectroscopy, the predominant functional groups were determined. [START_REF] Vandenbroucke | [END_REF]16,[START_REF] Siskin | Detailed structural characterization of the organic material in Rundle Ramsay Crossing and Green River oil shales[END_REF] Here we consider only those that have the possibility to bridge the inorganic matrix: hydroxyl, carboxylate, carbonyl, thiol, aromatic and aliphatic double bonds.

To assess the reactivity of these functional groups on the silica surface, and to determine which functional groups can serve as an anchor for kerogen, the following strategy was applied: i) Reaction energy value of a possible bond formation reaction between the surface and an organic molecule with a specific functional group was determined; ii) calculations were performed first with the smallest possible molecules with a given functional group, then the length of the alkyl side chain was increased until convergence in the reaction energy was reached. This way, bond formation reaction energy can be associated with every important functional group.

Note that we do not simulate the bond formation itself, we perform only geometry optimization and compare the energy of the reactants and the supposed products in any realistic reactions that one could expect on the surface of silica. [START_REF] Iler | The chemistry of silica[END_REF]19 Furthermore, we suppose that, in reality, the interface has long enough time to establish equilibrated bonds.

Details of the quantum chemical calculations

Considered reactions. The reactions considered in this study are schematized in a simplified way in Table 1. We considered possible reactions with all the functional groups listed in the previous section. We distinguished condensation reactions (in the case of hydroxyl, carboxylate and thiol groups) and addition reactions on the unsaturated carbonyl C=O and C=C double bonds. Hydroxyl, thiol, carbonyl groups as well as aliphatic carbon double bonds were considered both in normal (i.e. in the end of the chain) and iso position (i.e. non normal), as the position of these functional groups can be expected to influence the bond formation. Note that we did not consider this distinction in the case of the carboxylate group due to its size: the reactive acidic OH group is not directly connected to the chain but through a C=O group. Note that considering a C=O group in normal and iso positions means, respectively, considering ketones and aldehydes.

Note also that we presumed that the formed water molecules in the condensation reactions stay adsorbed on the surface relatively far from the attached organic molecule. The validity of this assumption was tested by exploring the most stable adsorbed positions of a water molecule around anchored organic molecules as well as on pure quartz surface. It was found that it is always more favorable for the water to get adsorbed on the bare quartz surface as this way it is able to form two strong H-bonds with the surface. Figure 2 illustrates one of the studied condensation reactions. Technical details. Quantum chemical calculations were carried out with the hybrid B3LYP DFT functional as it is known to properly describe covalently bonded compounds of main-group elements. [START_REF] Sousa | [END_REF]21 The non-periodic (i.e. gas phase) calculations were performed using localized basis sets: a double zeta basis set was applied for geometry optimization, whereas total energies were obtained on the previously optimized structures with a triple zeta set. To choose the optimal

Table 1. The considered chemical reactions in this study between surface silanol groups and representative organic functional groups of kerogen. Condensation and addition

+

→ + sets, we conducted tests with different basis sets on small carbon and silicon-containing molecules. Optimal geometries and vibrational frequencies were calculated and compared to high-level calculations obtained at the MP2 level of theory with the cc-pvtz basis set. These tests led us to choose the 6-31G** basis for geometry optimization and the 6-311++G** basis for the determination of energies. All calculations were performed using the Gaussian 09 program package. 22 2.3. Assessment of the performance of a reactive force field for application in further work [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF] Although quantum chemistry-based molecular dynamics is, in principle, a very accurate tool to study the dynamics of molecular systems, in practice, such simulations are very limited in terms of accessible time (tens of picoseconds) and system size (several hundreds of atoms). Due to computational inefficiency, we typically use quantum chemistry to gain only a static picture of molecular systems.

However, to consider mechanical cohesion and crack propagation at the kerogen-silica interface, one needs to simulate considerably larger molecular systems containing more than 10,000 atoms. Classical simulations based on an efficient empirical reactive potential can handle systems large enough over a time scale long enough to study the propagation of cracks. In our case, at the same time, such a potential must be able to account for changes of the complex chemistry at the kerogen-silica interface. Empirical reactive potentials rely on a complex mathematical formulation adapted to specific types of bonding situations. As being empirical, these potentials always need to be "trained". The training set might comprise experimental data as well as theoretically determined properties from high-level quantum chemical calculations.

One of the broadest existing reactive force fields, Re-axFF [23][24][25] was chosen for application in a future fracture study in kerogen/silica composites. [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF] Several parameterizations are available in the literature depending on the constituting atoms and the chemical nature of the system. The parameter sets used here were developed by van Duin et al. for small organic molecules 23 and Chenoweth et al. 25 As these parameters were developed for somewhat different systems than ours, it is necessary to make sure that the set of parameters is able to describe our systems with a reasonable degree of accuracy. To this end, 18 gas phase reactions were studied that are representative of processes taking place during crack propagation at the kerogen-silica interface. Optimized geometries and reaction energies were determined at the B3LYP/6-311++G** level of theory, and then compared to results obtained with ReaxFF. The considered reactions as well as the calculated reaction energies are represented in Table 2 and Figure 3.

The results presented in Table 2 and Figure 3 are quite consistent, although, in a few cases, the calculation with the ReaxFF potential differs significantly from the DFT calculation. Except for 3 reactions (italic in Table 2), the average relative error is lower than 20%. ReaxFF seems to somewhat underestimate the strength of single bonds, whereas it slightly overestimates double bond strengths. The trend in bond strength between Si-O and C-O bonds is well captured (the former being stronger, see reactions 6 and 7), which is im-portant for reliable simulations of bond formation and breaking in the systems considered.

Interface reconstruction

In this section, we present the methodology that we developed to propose realistic molecular reconstructions of the kerogen-silica interface. First we will show implications of the different kerogen properties on 8 different example interfaces with alpha-quartz. Our ultimate objective is to reconstruct an actual interface between a crystalline structure of bulk silica and a disordered structure of bulk kerogen by combining the computational performance of the ReaxFF potential and the list of favored reactions at the interface estimated by quantum chemical calculations. In the present work, we apply this methodology to a specific silica mineral and a specific molecular representation of bulk kerogen. Yet, such a methodology can be applied more generally to other silica and kerogen structures.

Bulk kerogen can be seen as a disordered porous carbon for which, to the best of our knowledge, no molecular structure has ever been developed that captures the elemental composition, the chemistry, and the 3D nanostructure, i.e., the mechanical, adsorption and transport properties. DFT energy, kcal/mol ReaxFF energy, kcal/mol the review of Vandenbroucke and Largeau, [START_REF] Vandenbroucke | [END_REF] the most advanced molecular representations have generally focused on the elemental composition and chemistry (H/C and O/C ratios, aromaticity, and functional groups), but neglected the threedimensional structure (see Section 8 of Ref 3 and reference therein). Some models were relaxed in three dimensions, but were not validated against experiments characterizing the actual nanostructure such as TEM images, diffraction experiments, or adsorption experiments. Therefore, it is unlikely that those relaxed models would be representative of the mechanical properties of the actual kerogen phase. However, these mechanical properties are crucial in the present work, since the prospective objective of our reconstructions is to enable the simulation of crack propagation at the kerogen-silica interface. [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF] As an alternative to model bulk kerogen, we considered a molecular model of disordered porous carbons developed by Hybrid Reverse Monte Carlo (HRMC): CS1000. 26,27 HRMC is a molecular reconstruction technique based on X-Ray diffraction data validated against adsorption experiments and TEM images. This procedure ensures that the 3D nanostructure is consistent with the actual structure of the material considered. The CS1000 model is not a model of kerogen per se, but of pyrolyzed saccharose. Accordingly, CS1000 is not representative of all the properties of kerogen. In particular, it does not contain oxygen and its hydrogen content is lower than typical values of kerogen. Nonetheless, the three-dimensional structure, which certainly plays a critical role in determining mechanical properties, is typical of disordered porous carbons.

The mechanical properties of such a porous carbon mostly depend on the carbon skeleton, whereas the elemental composition is secondary. In fact, the chosen CS1000 can be regarded a very mature dense kerogen with low oxygen and hydrogen contents. In the reconstruction study, for the silica phase, we considered alpha-cristobalite, a silica polymorph with a tetragonal lattice that can easily adapt to the cubic lattice of CS1000. Therefore, it is possible to consider the interface between alpha-cristobalite and CS1000, as an example of kerogensilica interfaces.

The procedure we followed to reconstruct the interface between CS1000 and alpha-cristobalite has three steps:

1.

The silica surface is covered with reacted organic functional groups and Si-OH terminations. In practice this means replacing the original Si-OH groups with reacted surface terminations (e.g. -C(=O)-O-Si groups) in a random manner but following an equilibrium distribution: The amount of a given organic termination at the surface is determined based on the concentration of reactive functional groups at the surface of kerogen. Reacted functional groups are put on the surface 'manually', by assuming equilibrium bond lengths and angles. The effect of temperature is accounted for through the Boltzmann statistics provided that the system had a long enough time to establish equilibrium bonding structure. The probability that a reactive group at the surface of kerogen actually reacts with the silica surface is:

( ) (
)

RT E RT E P form bond form bond bond . . exp 1 exp - + - = (1)
where E bond form. is the energy of reaction of the functional group considered, R is the universal gas constant and T is the temperature. Note that our reference is the non-bonded case, consequently the Boltzmann factor corresponding to this state of the system is unity. The remaining terminations on the silica surface are the unreacted Si-OH groups. The effect of the concentration of reactive functional groups at the surface of kerogen will be discussed in the next section..

2.

The bulk kerogen model and the silica surface covered with reacted organic functional groups are positioned one in front of the other at a distance of 5Å. To bind the two phases together the empty space between them is filled randomly with atoms representing the same density and elemental composition as the considered kerogen model.

3.

The atoms in the 5 Å-thick layer are relaxed following an annealing procedure with the ReaxFF potential. The atoms outside the 5 Å layer are fixed at their initial position. The atoms initially inside the 5 Å layer cannot move farther than 2.5 Å outside the 5 Å-thick transient region.

Results and Discussion

Bond formation energies

In this section, we present the results for bond formation energies of the typical functional groups at the kerogen surface obtained from quantum chemical calculations. In all cases considered, the calculations were performed for a side chain length of up to 15-20 carbon atoms. Bond formation energies shown in Table 3 were obtained by taking the average of the 5-7 last members of each series. As an example, Figure 4 illustrates this procedure in the case of the carboxylate group. The results show that covalent bonds can be expected between the kerogen and the quartz surface. According to the values in Table 3, bond formation with carboxylate and alcoholic hydroxyl groups (regardless of the position of the OH group in the chain) are very favorable, leading to the formation of -C(=O)-O-Si and >C(H)-O-Si bonds, respectively (see Table 1). Indeed, the esterification reaction with alcohols is one of the most well-known surface modification reactions of silica particles. [START_REF] Iler | The chemistry of silica[END_REF] The condensation reactions of thiols, i.e., the for- Number of carbon atoms Bond formation energy, kcal/mol E bond form. = -(6.44 ± 0.21) kcal/mol mation of Si-S bonds, are not favored: breaking a very strong Si-O bond is not compensated with the formation of a weak Si-S bond and an adsorbed energetic water molecule. This is in accordance with the finding of Ossenkamp et al. who treated silica particles with the solution of alcohols functionalized with mercapto (SH) group. 28 These molecules were found to form bonds to the surface preferentially with their OH termination.

Figure 4.

Reaction energies of bond formation reactions between the quartz surface and organic molecules containing a carboxylate group as a function of increasing side chains. The last six values were averaged.

Addition of a surface silanol group on the carbonyl group (C=O) is not favored if the carbonyl group is in iso position (case of ketones), but becomes somewhat favored for normal carbonyls (i.e. aldehydes). This latter functional group leads to the formation of -C(H)(OH)-O-Si bonds (see Table 1). Addition of a surface silanol group on unsaturated carbon-carbon bonds turned out to be favored only if the bond is aliphatic. Indeed, Broge provided experimental proofs for bond formation between heat activated silica surface and both normal and iso-alkenes. [START_REF] Iler | The chemistry of silica[END_REF]29 Note that in this study, cis iso-alkenes were considered in the gas phase. In fact, this isomer is more adapted for the more probable curved conformation of the chain and also provides higher accessibility for surface silanols. Reactions on aromatic systems, as one could expect, are found to be highly unfavorable. Note that we considered here only the first 4 aromatic hydrocarbons. The lowest bond formation energy was found for position 9 of phenanthrene, which is known as the position of the least aromatic character (and therefore associated with the highest reactivity), while the highest value was obtained for position 2 (i.e. beta) of naphthalene. It is reasonable to assume that bond formation energies would not turn negative for larger aromatic moieties, as one cannot expect sites in kerogen with considerably less aromatic character than position 9 in phenantrene. In fact, the typical aromatic sheets were found to be no larger than 4-6 condensed aromatic rings. [START_REF] Vandenbroucke | [END_REF] As expected, in all cases the functional groups in iso position are somewhat less reactive compared to those in normal position due to steric effects. The bond formation reaction energy values given in Table 3 correspond to reactions taking place on the fully hydroxylated (001) alpha-quartz surface. It is worth noting that the density of surface silanol groups on different crystallographic quartz surfaces is different, which could result in different atomic configurations at equilibrium and also different water adsorption energy. This might lead to somewhat different bond formation reaction energies.

Consequences for kerogen-quartz interfaces

With the results of the previous section, it is possible to estimate the approximate concentration and the type of the bonds at the interface between kerogen and silica. Many factors can influence this concentration, such as, among others, kerogen type, maturity, origin, density and burial depth. In this section, we discuss the most important factors and their impact on the concentration and type of bonding at the kerogen-silica interface. We also provide estimations on the number of different bonds depending on those factors.

If a functional group at the surface of kerogen is facing a hydroxylated silica surface, as we have seen, condensation or addition reactions may occur. The approximation that we consider here is that the reaction occurs according to the Boltzmann statistics (Eq. 1). Equation (1) assumes that the reaction occurs independently of the environment. In reality, the reaction of a functional group locally modifies the kerogen structure and therefore, strictly speaking, it is not an independent process. A rigorous approach would be to simulate a whole kerogen molecular structure with surface functional groups over very large time scales, which is not feasible. An argument that would support the validity of the approach we consider is that over very large time scales, like the geologic time of the maturation of the kerogen, the intrinsic creep behavior of the organic matter would relax any local stress generated by the reaction of a functional group.

According to the Boltzmann statistics, the reactions for which the energy in Table 3 is negative are favored; those for which the energy is positive are not favored. At 300K, reactions 1.a, 1.b and 2.b occur with a probability higher than 99%, i.e., almost in all cases. However, reactions 2.a in nor-mal position (i.e. case of aldehyde) occurs with a probability of 91.0%, i.e., a non-negligible proportion of the functional groups does not react. The other functional groups, with positive reaction energies, have a probability to react lower than 0.3%.

In order to translate those probabilities of reaction into the concentration of bonds at the silica-kerogen interface, one needs to know the concentration of functional groups at the surface of kerogen. The surface concentrations of reactive functional groups depend mainly on the elemental composition, the repartition of functional groups in the bulk, the density and the maximum distance between a functional group and the silica surface for the reaction to occur.

Fourier transform infrared spectroscopy (FTIR) combined with elemental analysis makes it possible to determine the concentration of oxygen-containing functional groups in bulk kerogen. Here, we considered the FTIR analysis of Robin and Rouxhet who made a distinction between the different types of kerogen (see Figure 5). 16 The repartition proposed by Robin and Rouxhet shows significant differences from one type of kerogen to the other. In particular, type I and type II kerogens have a large content of ester (-C(=O)O-) which is a bridging group and as such non-reactive. The other categories are reactive or partially reactive (aldehyde / ketone). Therefore, one may expect the type of the kerogen to be critical factor for the bonding at the kerogen-silica interface. FTIR does not distinguish between the normal and iso positions of the functional groups. In our work, we assumed that half of the functional groups are in normal and the other half in iso positions. As for carbon, 13 C NMR enables estimating the proportion of aromatic carbon. In this work, we considered the analysis proposed by Behar and Vandenbroucke who determined the relative content of aliphatic, naphtenic and aromatic carbons in the function of the kerogen type and maturity (see Figure 5b). 30 In addition to 13 C NMR, Behar and Vandenbroucke also considered mass spectroscopy data to estimate the relative contents of aliphatic and naphtenic carbons. Note that in the previous sections, we did not consider naphtenic carbons as we assumed they are not very different from aliphatic carbons from the point of view of the surface chemistry of kerogen. However, to the best of our knowledge, no quantitative estimate of the relative content of aliphatic sp 2 vs. sp 3 carbons exists in the literature for kerogen; although some techniques such as Raman spectroscopy have this capacity in theory. In the present work, we decided to make a rough estimate: we assumed that one third of the aliphatic C-C bonds are double bonds, half of them in normal configuration, the other half in iso configuration.

The elemental composition, i.e., the atomic ratios O/C and H/C, is well-known for kerogen and serves as a criterion to determine type and maturity. The van Krevelen diagram (see Figure 6a) is a schematic diagram relating the O/C and H/C ratios for the three main types of kerogen at different level of maturation. [START_REF] Tissot | Petroleum Formation and Occurrence[END_REF] Immature type I and II kerogens have typically high H/C ratios (~1.5) and moderate O/C ratios (~0.1 to 0.2), whereas immature type III kerogens have low H/C ratios (~0.8) and large O/C ratios (~0.3). Maturation of kerogen tends to reduce both ratios while the rock is producing oil, gas, water and carbon dioxide. Very mature kerogens have low H/C and O/C ratios (<0.5 and <0.1, respectively). According-ly, the amount of reactive functional groups in kerogen depends significantly on the level of maturity. 16 The percentages represent the proportion of oxygen atoms in each functional group (a). Repartition of aliphatic, naphtenic and aromatic carbons in kerogens of different types and different maturity (left columns immature, right columns mature), from Behar and Vandenbroucke (b). 30 The density of kerogen ranges from 0.8 g/cm 3 to 1.5 g/cm 3 . [START_REF] Vandenbroucke | [END_REF] In addition, density and H/C ratios are usually anticorrelated, i.e., the density would increase with the maturity of the kerogen. [START_REF] Van Krevelen | Coal: Typology -Chemistry -Physics -Constitution[END_REF] However, Vandenbroucke and Largeau have reported high densities even for low maturities. [START_REF] Vandenbroucke | [END_REF] Therefore, the correlation between density and maturity seems valid on a relative scale for a given kerogen field, but not on an absolute scale for all kerogens. Here, we considered the density as an independent parameter. Combining FTIR, 13 C NMR, elemental analysis and density provides an estimation of the bulk concentration of functional groups. To convert bulk concentrations into surface concentrations, one needs to know the thickness of a region close to the surface, in which a functional group is close enough to the silica surface to react. Considering the typical size of the functional groups, the thickness of this region can be reasonably assumed in the range of 2Å to 5Å.

The factors that influence the surface concentration of reactive functional groups in kerogen are numerous. We propose to consider a few situations in order to identify some quantitative trends between the kerogen type, maturity, density and the bonds at the interface with alpha-quartz. We considered case 1 a reference system: it is an immature type II kerogen, of density 1.1 g/cm [START_REF] Vandenbroucke | [END_REF] , assuming a thickness of the reactive zone at the surface of kerogen of 3Å. Case 2 and 3 consider the same parameters but for type III and type I kerogens, respectively; case 4 is a high density kerogen; case 5 is a mature kerogen; case 6 assumes a thick reactive zone. Case 7 and 8 consider two extreme cases with few or many bonds at the interface. For the reference case 1, we estimated that 24% of the silanol groups reacted with the kerogen, producing a majority of >C(H)-O-Si terminations. The proportion of reacted silanol groups increases by 50% for a type III kerogen (case 2: 35%), and decreases moderately for a type I (case 3: 22.9%). The proportion of each type of termination changes slightly from one type of kerogen to the other, but the repartition remains similar with a majority of >C(H)-O-Si type bonds. The concentration of reactive functional groups is proportional to the density of kerogen and the thickness of the reactive zone as is the proportion of reacted silanols (see case 4 and 6, respectively). Accordingly, the proportion doubles between a low density kerogen (~0.8 g/cm 3 ) and a high density one (~1.5 g/cm 3 ). The same trend is valid between a thin reactive zone (2Å) and a thick one (5Å). Maturity reduces the oxygen content and therefore the concentration of oxygencontaining functional groups in kerogen. Accordingly, in case 5 the proportion of reacted silanols (11%) decreased compared to case 1. The last two cases (7 and 8) are extreme cases that illustrate the range of potential bonding between kerogen and alpha-quartz. The range is very large: 6% in case 7 and 80% in case 8. One may discuss the validity of the result in case 8: with such a high concentration of reactive functional groups in kerogen, assuming that each group reacts independently is probably unrealistic. Cases 1 to 6, provide a more realistic range between 11% and 41% of reacted silanols, i.e., 26% ± 15%. These estimations provide an approximate picture of the potential bonding at the interface between kerogen and silica. Since the fracture toughness and the cohesion directly depend on the bonding at the interface, the trends we identified should directly impact the macroscopic properties of silica-rich shales. However, many aspects were neglected in our approach and would need to be addressed for a more realistic study: i) Kerogen and silica were most likely in contact from the early stage of maturation. Therefore the bonding that formed at that time may have evolved, and this evolution may differ from the maturation of bulk kerogen. ii) We considered the same proportion of aliphatic carbon double bonds in all cases. A fine study by Raman spectroscopy would enable us to relate this proportion to the type and maturity of the kerogen. iii) Two close reactive functional groups on the kerogen sur-face may interact when reacting with the silica surface. To understand such interaction, combined reactions should be considered by quantum chemical calculations. iv) A comprehensive HRMC study imposing the right amount of atomic constituents and functional groups in the whole kerogen phase would be needed to more realistically account for the effect of changing kerogen composition at the interface as well as in the bulk structure.

Reconstruction of a realistic kerogen-silica interface

We applied the methodology presented in the Materials and methods section to actually reconstruct the molecular structure of kerogen-silica interfaces. In this section, we present the case of one specific situation and assess the performance of the reconstruction methodology. As already discussed, the molecular structure considered for bulk kerogen is CS1000, a porous carbon obtained by the Hybrid Reverse Monte Carlo technique. 22 Note that here the silica surface is represented with the (001) surface of alpha-cristobalite, whose tetrahedral lattice can easily adapt to the cubic lattice of CS1000. Here we assume that the bond formation energies of different functional groups are similar on the considered surface of alpha-cristobalite to those on (001) alpha-quartz as the surface concentration of silanol groups are similar on these surfaces: 8.1 nm -2 and 9.5 nm -2 , respectively. Reacted rganic surface terminations were randomly distributed on the alphacristobalite surface according to the precude presented Sections 2.4 and 3.2. The case we considered is a dense type II kerogen of intermediate maturity: 1.58 g/cm 3 (density of CS1000), H/C = 1.25, O/C = 0.1. We assumed a reactive zone of 3.5 Å and the aliphatic, naphtenic and aromatic carbon content of an immature type II kerogen. The corresponding concentration of reacted groups at the surface of alpha-cristobalite is similar to that obtained for cases 2 and 6 in Figure 6, i.e., a quite high level of bonding at the interface. We randomly distributed those terminations on an alpha-cristobalite surface (Figure 7). Then, as discussed in Section 2.4, we filled randomly the volume between this sur-face and the CS1000 structure with carbon and hydrogen atoms, and performed an annealing procedure (Figure 8). The annealing procedure consists in efficiently driving the system towards a minimum energy state by simulating with ReaxFF while simultaneously decreasing the temperature. Here, only the atoms in the transition layer were moved. We compared several temperature decrease rates (constant or logarithmic), profiles (step by step or continuous), and speeds (slow or rapid). In the procedure that led to the lowest total energy, the starting temperature was 600K which was decreased initially by constant steps of 50 K except at very low temperatures (a slower rate) until reaching the final 5 K. Each step was simulated for 2 ps. Higher initial temperatures, longer simulation steps or smaller temperature steps did not improve the results. We repeated the calculation 10 times considering other random initial configurations between the CS1000 structure and the alpha-cristobalite surface. All these attempts led to energies close to the first calculation with a standard deviation of 4.5 eV for a total of 823 atoms in the transition layer. By comparison, the other annealing procedures we tested led to final energies at least 170 eV higher.

We assessed the performance of the annealing procedure by computing the coordination number and the bond angles for the carbon atoms in the transition layer: for the 11 situations simulated, the carbon coordination number was 2.66 ± 0.03 and the average bond angles were 119.6° ± 0.9°, which are close to the coordination number and average bond angles of the bulk CS1000: 2.77 and 120°, respectively. 26 The coordination number in the transition layer was slightly smaller than in CS1000. We attributed this difference to the influence of the interface which locally affects the atomic environments. The so generated interface is subjected in a further study to determine fracture properties. 7

Conclusions

In this paper, the elements of a reconstruction technique of kerogen-silica interfaces was presented. In the first part of the study, the possible chemical bonds between kerogen and the fully hydroxylated (001) interface of alpha-quartz were studied by means of quantum DFT calculations. The reactivity of different typical functional groups on the surface of kerogen was investigated by calculating reaction energies of possible bond formation reactions. In these reactions, single molecules with a given functional group were considered with an alkyl side chain that was increased until the reaction energy did not show dependency on further increasing the size of the molecule. It was found that alcoholic hydroxyl, carboxylate, and aldehyde groups as well as aliphatic double bonds at the kerogen surface can form energetically favorable chemical bonds with the quartz surface. On the contrary, bond formation of thiols, ketones and aromatic bonds are not favored with surface silanol groups. Although experimental evidence supports bond formation between the functional groups considered here and silica surfaces, further experimental verification in gas shale samples would be desirable.

As the ReaxFF potential was chosen to be used to study the fracture properties of silica/kerogen interfaces, tests were performed to assess the reliability of this empirical force field in simple reactions representative of bond breaking and formation at the interface. The performance of ReaxFF was found satisfactory which makes it suitable for our approach. [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF] In the second part, combining estimates of functional group surface concentrations with the reaction energies obtained in the first part we presented a way to establish which chemical bonds and how many of them one should expect at the interface depending on the type and maturity of kerogen.

The consequences of different density, origin and maturity on the interfacial bonding were presented and discussed in 8 markedly different cases. Finally, the reconstruction of a realistic interface between a type II kerogen of intermediate maturity and silica was carried out by using the CS1000 model for kerogen and the (001) alpha-cristobalite surface. The annealing procedure to obtain the best optimized interfacial region was also discussed. The presented strategy is the basis for a comprehensive study aiming to determine the impact of the established silica-kerogen interface on the fracture toughness of the composite system. [START_REF] Brochard | Fracture Toughness Calculation by Molecular Simulation: Cases of Silica, Microporous Carbons and their Interface[END_REF]

Figure 1 .

 1 Figure 1. The studied silica clusters (a, b: H 48 O 76 Si 26 and c, d: H 76 O 146 Si 54 ) representing a fully hydroxylated alpha-quartz surface. Atoms whose position can be changed during the optimization are marked in blue (b, d). The yellow, white and red colors correspond to Si, H and O atoms, respectively. Reactions were considered only with OH groups on the top of the clusters. Our tests confirmed that the size of the smaller cluster is adequate for our study.

reactions were distinguished. 1 .Figure 2 .

 12 Figure 2. Illustration of the studied condensation reactions. In the initial state, the reagents are considered infinitely far from each other. After the reaction takes place, the chemical bond is formed between the organic molecule and the quartz surface. The adsorption of the formed water molecule is found to be more favorable on the pure quartz surface than around the organic functional group.

Figure 3 .

 3 Figure 3. Comparison of the reaction energy values obtained at the B3LYP/6-311++G** level of theory and the ReaxFF reactive force field. The straight line shows what would be the case of a perfect match.

Figure 5 .

 5 Figure5. Repartition of oxygen-containing functional groups in kerogens of different types, from Robin and Rouxhet.16 The percentages represent the proportion of oxygen atoms in each functional group (a). Repartition of aliphatic, naphtenic and aromatic carbons in kerogens of different types and different maturity (left columns immature, right columns mature), from Behar and Vandenbroucke (b).30

2 . 4 .*

 24 Type I/II/III mature and density as well as the thickness of the reactive zone. Figure 6a illustrates how our example systems can be distinguished in terms of type and maturity in the van Krevelen diagram. In Figure 6b, we display the proportion of the corresponding surface terminations that we computed considering the alpha-quartz surface with an initial density of silanol groups of 9.56 nm -Table Different cases considered to identify trends between kerogen properties and bonding at the kerogenalpha-quartz interface. Details of the atomic ratios considered: Type II immature: H/C = 1.41, O/C = 0.18 Type III immature: H/C = 0.88, O/C = 0.28 Type I immature: H/C = 1.67, O/C = 0.13 Type I / II / III mature: H/C = 0.5, O/C = 0.06

Figure 6 .

 6 Figure 6. Position of different example kerogens (see Table 4) in the van Krevelen diagram (a). Repartition of the chemical bonds at the α-quartz-kerogen interface for the eight cases listed in Table 4 (b).

Figure 7 .

 7 Figure 7. Alpha-cristobalite surface with organic surface terminations randomly distributed.

Figure 8 .

 8 Figure8. The transition layer between the alpha-cristobalite surface and the CS1000 structure is randomly filled representing the composition and density of CS1000. This transition layer is then relaxed with an annealing procedure.

Table 2 . Test reactions and reaction energies as obtained at the B3LYP/6-311++G** level of theory and the ReaxFF reactive force field.

 2

Table 3 . Bond formation energies obtained for typical functional groups of kerogen on the (001) surface of alpha- quartz. Results are obtained from DFT calculations (see text)

 3

	Functional group	Bond formation energy / kcal/mol
	1.a Hydroxyl	normal	-9.68 ± 0.11
		iso	-5.13 ± 0.75
	1.b Carboxylate		-6.44 ± 0.21
	1.c Thiol	normal	3.52 ± 0.28
		iso	4.52 ± 0.32
	2.a Carbonyl	normal (Aldehyde)	-1.38 ± 0.21
		iso (Ketone)	5.4 ± 1.7
	2.b Aliphatic dou-	normal	-6.49 ± 0.23
	ble bond	iso	-2.77 ± 0.56
	2.c Aromatic bond		16.0 -46.1 a

a The 4 smallest aromatic hydrocarbons were considered

Table 4

 4

	shows eight

ACKNOWLEDGMENT

The authors are grateful to Elixabete Rezabal and Marie-Laure Bocquet for useful discussions. We also thank funding from Royal Dutch Shell and Schlumberger through the MIT X-shale research project.

AUTHOR INFORMATION