# The Impact of Pore Structure on Kerogen Geomechanics

Authors: Saad Alafnan
Journal: Geofluids
 (2021)

---

## Abstract

Production stimulation techniques such as the combination of hydraulic fracturing and lateral drilling have made exploiting unconventional formations economically feasible. Advancements in production aspects are not always in lockstep with our ability to predict and model the extent of a fracturing job. Shale is a clastic sedimentary rock composed of a complex mineralogy of clay, quartz, calcite, and fragments of an organic material known as kerogen. The latter, which consists of large chains of aromatic and aliphatic carbons, is highly elastic, a characteristic that impacts the geomechanics of a shale matrix. Following a molecular simulation approach, the objective of this work is to investigate kerogen’s petrophysics on a molecular level and link it to kerogen’s mechanical properties, considering some range of kerogen structures. Nanoporous kerogen structures across a range of densities were formed from single macromolecule units. Eight units were initially placed in a low-density cell. Then, a molecular dynamic protocol was followed to form a final structure with a density of 1.1 g/cc; the range of density values was consistent with what has been reported in the literature. The structures were subjected to petrophysical assessments including a helium porosity analysis and pore size distribution characterization. Mechanical properties such as Young’s modulus, bulk modulus, and Poisson ratio were calculated. The results revealed strong correlations among kerogen’s mechanical properties and petrophysics. The kerogen with the lowest porosity showed the highest degree of elasticity, followed by other structures that exhibited larger pores. The effect temperature and the fluid occupying the pore volume were also investigated. The results signify the impact of kerogen’s microscale intricacies on its mechanical properties and hence on the shale matrix. This work provides a novel methodology for constructing kerogen structures with different microscale properties that will be useful for delineating fundamental characteristics such as mechanical properties. The findings of this work can be used in a larger scale model for a better description of shale’s geomechanics.


---

## Body

## 1. Introduction

Rock geomechanics is a crucial branch of oil and gas reservoir development. The role of geomechanics is even more influential in unconventional formations, where permeability is mostly attributed to the hydraulic fractures generated. While hydraulic fracturing technologies have advanced, allowing for the efficient stimulation of newly drilled wells and restoration of productivity in old wells, the ability to describe and model the initiation, propagation, and dynamics of induced hydraulic fractures is lacking. Shale is a highly heterogeneous sedimentary rock consisting of a complex mineralogy of clay, quartz, calcite, and fragments of organic matter known as kerogen [1–3]; it is this complexity that causes shale petrophysics and geomechanics to deviate from those of conventional sedimentary rocks.The extent and pattern of hydraulic fractures are determined by the fracturing fluid, as well as the geology of the formation. Highly viscous fluids are used to increase the propagation of the fracture, while low-viscosity fluids are employed to bridge formation discontinuities [4, 5]. A high rate of injectivity is used to enlarge the stimulated reservoir volume SRV, while a lower rate targets interconnecting discontinuities [6]. The formation’s geological characteristics, however, fall beyond predesign control. Hydraulic fractures are likely to progress in the direction of existing natural fractures, where stress accumulation is the highest [7];though, clay minerals and other fragments may act as a barrier to propagation [8, 9]. Fracture geometry is controlled by the depth of the well. Shallow wells tend to form a network of fractures rather than a biwing fracture [10, 11]. The contents of brittle materials, characterizable by lower ranges in the Poisson ratio, correlate positively with the formation of a complex fracture geometry [12, 13]. Generally, four types of fractures are anticipated: simple, fishbone-like, fishbone-like with a fissure opening, and multilateral fishbone-like. The fracture category is determined by several factors, including bedding discontinuity, the presence of natural fractures and faults, and rock mineralogy [14]. The macroscopic mechanical properties of shale depend on the properties of the different constituents on a microscopic level [15, 16]. In general, shale is a composite materials, with clay and quartz serving as the matrix and organic matter present as inclusions [17, 18].The geomechanics of shale has been determined through traditional uniaxial and triaxial testing methods and laboratory scale fracturing experiments [19–24]. Such tests are based on the qualitatively linking of mechanical properties to the mineralogy. Hence, the models developed are vulnerable to a high level of discrepancy when applied to other samples with different mineralogical distributions [25]. Micromechanical testing methods such as quasistatic nanoindentation, mapping by nanoindentation, and atomic force microscope AFM have been employed to study variations in mechanical properties, with shale matrices revealing highly heterogeneous behaviors of feldsbar, carbonate, quarts, and organic matter [26–34]. Young’s modulus of inorganic inclusions is in the range of 50 GPa, while it is less than 10 GPa for organic matter [27–29]. Inorganic constituents exhibit higher anisotropic behavior compared to those that are organic in nature [35–37]. Microscopic mechanical properties vary beyond the elastic region, impacting the development of hydraulic fractures. The extent and complexity of the induced hydraulic fractures are determined by the microstructural constituents of the shale matrix [38–40].Organic matter, characteristically different from inorganic components, is often the target of micromechanical testing. High and low organic matter in shale samples have been tested for tensile failures, unconfined tensile strength, and fracture behavior, revealing some dependency of the geomechanics on the contents of the organic matter [41]. The heterogeneity and complexity of the mineralogy render isolating the impact of organic matter on geomechanics rather challenging. The organic matter consists of kerogen and bitumen. While the latter is soluble in organic solvents, the former represents a solid, best thought of as micro- or mesoporous material [3]. Chemical bonding and intermolecular interactions change with kerogen type and level of maturation, altering the chemical and physical properties (including the geomechanical behavior) and hence impacting the macroscopic geomechanics of shale. Experimental characterization of such kerogen requires extraction, followed by typical petrophysical and geomechanical assessments. Obtaining comprehensive insights into the relationship between the kerogen’s structure and geomechanics means gathering samples containing kerogen that are characteristically different, a task that is both laborious and dubious.The ability to represent naturally occurring or synthetically derived structures on a computational molecular platform has provided an unprecedented degree of flexibility in the study of complex research problems, moving the practice well beyond laboratory-based work. The goal of the present research is to recreate a range of representative kerogen configurations on a molecular level to delineate the relationship between kerogen’s structure and its mechanical behavior.


### 1.1. Molecular Simulation of Kerogen

The first attempt to classify kerogens was based on the observation that kerogens contained in coals share a great deal of similarity with the matrix itself (e.g., a graphene-like structure), while kerogens found in other formations contain higher fraction of aliphatic chains [42, 43]. The current classification of kerogens relies on detailed elemental and spectroscopic analyses. In particular, hydrogen to carbon (H/C) and oxygen to carbon (O/C) ratios have been interpreted to determine the degree of aromaticity and oxidation level. Hence, kerogens can be ordered according to their origin (type) and rank (maturation) [44, 45]. Kerogens, generally, are categorized as type I, II, III, or IV, based on their origin. Each kerogen type is subject to a sequence of diagenesis, catagenesis, and mutagenesis processes. Hence, kerogens exist on a spectrum of maturity level [45, 46]. Some references classify type IV as the same as type III, due to their intrinsic similarities [47]. Type I is characterized by relatively high H/C and low O/C ratios, making it more aliphatic in nature [47]. Conversely, type II contains a larger amount of aromatic and cyclic aliphatic fractions, which decreases the H/C ratio [48]. Being originated in an oxygen-rich environment, type III is associated with larger contents of oxygen appearing in the functional groups or connecting cyclic chains of aromatic moieties [48, 49].Through identifying the primary functional groups, complete modeling of the kerogen chemical structures was achieved where two distinct groups of kerogens were recognized. The first was island-like aromatic structures fused by alkoxy or ether while the second was more aliphatic. Large building blocks, also known as macromolecules, were derived based on hypothetical models constrained to match elemental ratios [50–54]. These kerogen macromolecules, when upscaled, failed to match macroscopic experimentally measured parameters such as density, a condition at least partially attributable to the lack of proper description of their three-dimensional (3D) spatial configuration [55]. Nevertheless, representing kerogens from different origins and thermal maturity levels has been achieved for type I and type II at diagenesis, beginning, and end of catagenesis [56, 57]. With the computational power of computer-based techniques, more types have been obtained. Ungerer et al. [58, 59] derived six representative models (i.e., IA, IIA, IIB, IIC, IID, and IIIA) [58, 59]. It is noteworthy to state that the prototypes listed in this review and some others were meant to provide structures that could reproduce the general macroscopic characteristics of kerogens to facilitate thorough investigation of different petrophysical, transport, and other related aspects. In this research, the advancement in the molecular representation of kerogen has been utilized to create representative structures, facilitating a thorough geomechanical analysis.


### 1.2. Atomistic Simulations of Mechanical Properties

Theoretical calculations of the mechanical properties of materials have attracted increasing interest for use in a wide range of engineering applications. Particularly, they have been adopted to study advanced composite materials under a wide range of conditions. The steering principle behind the approach is based on relating structural deformation to the applied stress field, after restoring mechanical equilibrium through energy minimization [60]. Theodorou and Suter [60] conducted a mechanical behavior study of atactic polypropylene at 233 K, considering 15 independent configurations; the authors found Young’s modulus to be within 10% of what was measured experimentally. Commercial software packages have facilitated the creation of various amorphous structures. Fan and Hsu [61] conducted a molecular characterization of aromatic polysulfone elastic moduli within an acceptable range of accuracy and uncertainty. Hutnik et al. [62] considered a more complex system of glassy polycarbonate of 4,4′-isopropylidenediphenol for an analysis of the elastic and plastic properties, observing a high degree of uncertainty attributable to the number of configurations considered. Fan et al. [63] studied a similar system and reduced the uncertainty by considering more configurations. Raaska et al. [64] calculated elastic constants, using both static and dynamic stress-strain concepts. The reported Young’s modulus value was comparable to that which was obtained experimentally. Kerogen, a set of naturally occurring materials consisting of carbon, hydrogen, oxygen, nitrogen, and sulfur and connected chemically through aliphatic, cyclic, and aromatic structures, was found to primarily be similar in chemistry to the aforementioned polymerized compounds. The success achieved in the calculation of the mechanical properties of polymers rendered an atomistic approach a reliable proxy, providing useful insights into kerogen geomechanics and the factors controlling it [65–70]. These studies shed some light on the mechanical properties of kerogen in a general sense. However, a detailed study of how kerogen heterogeneity influences the geomechanics of kerogen is not exhaustive. In this study, the impact of pore structure is linked to the mechanical behavior through an integrated study considering some range of kerogen configurations, temperature, and fluid effects.In the next section, the nanostructure construction of kerogen is detailed, followed by the molecular algorithms utilized to characterize the structure for its relevant petrophysical aspects. The mechanical properties simulation adopted in this research is also presented.


### 1.3. Realistic Kerogen Structures

To study kerogen’s geomechanics, a realistic structure epitomizing kerogen is required to serve as nanoporous media. In this research, a novel molecular modelling approach was followed to construct a nanoporous kerogen model virtually. The extraction of organic matter in sufficient quantities for a petrophysical assessment is arduous. The process generally requires the application of etching agents and/or thermal treatment. Additionally, the desired conditions of pressure and temperature at which different analyses are performed can be experimentally challenging. Attributed to its unprecedented flexibility, molecular modeling framework has the advantage of controlling the petrophysical properties of kerogen during the initiation phase where different kerogen types, bonding strengths, and other chemical properties can be changed in order to independently isolate their impact. The geomechanical assessments carried out in this study relied on the novel work of the kerogen models derived by Ungerer et al. [59]. These units are understood to cover the basic kerogen types (see Figure 1 and Table 1).Figure 1
(a) LAMMPS protocol workflow for creating kerogen nanostructures. (b) from a single kerogen unit to the final kerogen structure [82].
(a)(b)Table 1
Kerogen prototypes presented by Ungerer [58], reproduced from Alafnan [82].
Kerogen typeMaturity levelH/CO/CN/CS/CKIAImmature1.5300.0520.0280.012KIIAImmature1.1700.0950.0240.012KIIBTop of oil window1.1200.0600.0220.009KIICMiddle end of oil window0.9050.0540.0210.008KIIDOvermature0.5800.0510.0230.011KIIIAImmature0.8700.1160.0170.000In this study, kerogen type IIA was used as a representative of the organic matter in shale. Type IIA kerogen macromolecules were recreated on a molecular level, see Figure2 for a graphic representation of the various kerogen types. Polymer-consistent forcefield plus (PCFF+) was used to parameterize the atoms present in kerogen as center forces [71, 72]. Each center is defined and assigned a charge. Lennard Jones 6-9 was employed to describe the interactions between different center forces. The PCFF+ is very well known for its ability to configure organic compounds complementing earlier efforts of forcefields such as with COMPASS and CFF93 [71, 72].Figure 2
Van Krevelen diagram of kerogen maturation [59].An initial low-density structure containing eight units of kerogen IIA was created in relatively large cell volume to avoid any issues caused by larger sized kerogen macromolecules. The final kerogen structure was then formed utilizing a large-scale atomic/molecular massively parallel simulator LAMMPS. The simulation proceeded to a final condition of 20.67 MPa and 350 K which are deemed to represent typical reservoir conditions. The molecular dynamic steps began with initializing the cell by assigning 9.5 cutoff value and 2.0 skin with periodic boundaries, that was followed by a proper configuration of the molecular positions and velocity through energy minimization. Next, two stages of isochoric-isothermal NVT and isobaric-isothermal NPT for 250 ps and 200 ps, respectively, were run at 900 K. Then, the temperature was gradually reduced from 900 K to the final temperature through three consecutive NPT stages. The gradual decrease to the final targeted temperature assured faster convergence of kerogen units. The density of the final kerogen structure was around 1.1 g/cm3, which is close to the measured experimental one for the same type [59]. During the creation of the initial cell, a number of dummy molecules were added to act as a placeholder. The dummy molecules were removed when the final structure was formed, leaving behind a certain number of voids. The same process was repeated with an increasing number of dummy molecules, yielding a total of seven configurations across a range of density values (see Table 2). A description of the molecular construction protocol is given in Figure 1.Table 2
Seven kerogen configurations constructed.
ConfigurationNumber of dummy moleculesDensity (g/cm3)Porosity101.1146.45%251.0968.11%3101.0829.36%4151.05811.43%5201.04113.48%6251.02314.77%7350.96319.68%


### 1.4. Petrophysical Characterization

Helium–kerogen interactions in a Gibbs Monte Carlo simulation was used as a proxy for determining the void space, through the following relationship:(1)Nm=Na−ρaVP.Here, ρa is the density of the helium, Nm is the number of excess helium, and Na is the number of adsorbed molecules of helium. The pore volume VP can be calculated assuming no excess adsorption which is reasonably justified by the confined spaces in kerogen. Once the pore volume is estimated, the porosity could then be calculated when divided by the bulk volume. The calculations were performed on all kerogen structures, revealing increasing void ratios as the number of dummy molecules increased (see Table 1).Porosity, which is an indicator of the total pore volume, can be used to describe the total voids of porous media. A detailed distribution of pore sizes, however, is needed to fully characterize a kerogen structure. For this purpose, a pore size distribution (PSD) analysis was carried out for all seven nanoporous kerogen structures. The insertion of nonoverlapping spheres to fill the void spaces was used to determine the PSD distribution. The process involves predefining a threshold radius (i.e., a sphere was created in a certain location if its radius was greater than or equal to 0.15 nm in this study) [73]. An example of the PSD calculations performed can be found in Figure 3. The PSD and cumulative distribution function are given in Figures 4 and 5, respectively. The analysis revealed an increasing number of larger pores as the number of dummy molecules increased (i.e., larger pores were found to be positively correlated with porosity).Figure 3
Visual illustration of the PSD algorithm for one of the kerogen structures [82].Figure 4
PSD of the kerogen configurations.Figure 5
Cumulative distribution function of the PSD for the seven configurations.


### 1.5. Mechanical Properties

The kerogen’s molecular structure, initially at equilibrium, was deformed by changing the lattice vectorR=a,b,c to a disturbed cell R′=a′,b′,c′, using a strain matrix e [74–76].
(2)R′=R1+exx12exy12exz12eyx1+eyy12eyz12ezx12ezy1+ezz.The changes imposed led to a shift in the total energy of the system:(3)U=Etot−E0V0=12∑i=16∑j=16Cijeiej,whereU is the stress per unit volume, V0 is the initial volume of the cell, E0 is the energy of prior deformation, Etot is the energy after deformation, and Cij is the element of the stiffness matrix. Both i and j run from one to six {xx, yy, zz, yz, xz, xy}, yielding a stress tensor of 36 elements. Assuming symmetry, the total number of elements was reduced to three for a cubic cell.A predefined straine was defined across one direction (i.e., one component of the lattice). All other strain components were set to zero. The energy change U, defined by Equation (3), reduced as follows [77].
(4)U=C11e22.The above equation was used to calculate element C11. Similarly, the diagonal strains was set as eyz=ezy=1/2e, and all other strain components were set to zero. The energy change equation could then be reduced to calculate C44 [77].
(5)U=C44e22.The bulk modulus, an elastic coefficient representing the response to a uniform applied stress, was calculated by settingexx=eyy=ezz=e, while the shear modulus, a coefficient of the elastic response to shear stress, was estimated by setting ezz=e and exx=eyy=−1/2e [78]. For estimating the elastic moduli, smaller value of strains is used. Elastic moduli are obtained at two values of strains to assure consistency (i.e., 0.02 and 0.04 strains were selected).The aforementioned approach was applied to the different kerogen configurations. The strain implementation and energy calculations were conducted through LAMMPS and assisted by a MedeA interface. The molecular dynamic protocol consisted of initialization with a 3D periodic boundary, 9.5 cutoff distance, and 2.0 skin, followed by energy minimization to optimize the positions of the atoms. Then, strains were applied in predefined directions. The energy was subsequently minimized, and NVT ensemble was run at the specified temperature, allowing for a calculation of the total change in energy and stiffness matrix. The data were then used to estimate the elastic moduli. A summary of the molecular approach is given in Figure6.Figure 6
LAMMPS-based molecular approach used to calculate mechanical properties.The above concept was extended beyond the elastic region by continuously deforming the structure at an incremental strain and minimizing the energy before and after deformation to mimic tensile, compression, and shear failure tests.Kerogen is anticipated to exhibit both elastic and plastic behaviors during the shale stimulation phase and production span. In this study, elastic property and tensile failure tests were carried out to develop a comprehensive understanding of kerogen’s geomechanics. The details are presented in the next section.


## 1.1. Molecular Simulation of Kerogen

The first attempt to classify kerogens was based on the observation that kerogens contained in coals share a great deal of similarity with the matrix itself (e.g., a graphene-like structure), while kerogens found in other formations contain higher fraction of aliphatic chains [42, 43]. The current classification of kerogens relies on detailed elemental and spectroscopic analyses. In particular, hydrogen to carbon (H/C) and oxygen to carbon (O/C) ratios have been interpreted to determine the degree of aromaticity and oxidation level. Hence, kerogens can be ordered according to their origin (type) and rank (maturation) [44, 45]. Kerogens, generally, are categorized as type I, II, III, or IV, based on their origin. Each kerogen type is subject to a sequence of diagenesis, catagenesis, and mutagenesis processes. Hence, kerogens exist on a spectrum of maturity level [45, 46]. Some references classify type IV as the same as type III, due to their intrinsic similarities [47]. Type I is characterized by relatively high H/C and low O/C ratios, making it more aliphatic in nature [47]. Conversely, type II contains a larger amount of aromatic and cyclic aliphatic fractions, which decreases the H/C ratio [48]. Being originated in an oxygen-rich environment, type III is associated with larger contents of oxygen appearing in the functional groups or connecting cyclic chains of aromatic moieties [48, 49].Through identifying the primary functional groups, complete modeling of the kerogen chemical structures was achieved where two distinct groups of kerogens were recognized. The first was island-like aromatic structures fused by alkoxy or ether while the second was more aliphatic. Large building blocks, also known as macromolecules, were derived based on hypothetical models constrained to match elemental ratios [50–54]. These kerogen macromolecules, when upscaled, failed to match macroscopic experimentally measured parameters such as density, a condition at least partially attributable to the lack of proper description of their three-dimensional (3D) spatial configuration [55]. Nevertheless, representing kerogens from different origins and thermal maturity levels has been achieved for type I and type II at diagenesis, beginning, and end of catagenesis [56, 57]. With the computational power of computer-based techniques, more types have been obtained. Ungerer et al. [58, 59] derived six representative models (i.e., IA, IIA, IIB, IIC, IID, and IIIA) [58, 59]. It is noteworthy to state that the prototypes listed in this review and some others were meant to provide structures that could reproduce the general macroscopic characteristics of kerogens to facilitate thorough investigation of different petrophysical, transport, and other related aspects. In this research, the advancement in the molecular representation of kerogen has been utilized to create representative structures, facilitating a thorough geomechanical analysis.


## 1.2. Atomistic Simulations of Mechanical Properties

Theoretical calculations of the mechanical properties of materials have attracted increasing interest for use in a wide range of engineering applications. Particularly, they have been adopted to study advanced composite materials under a wide range of conditions. The steering principle behind the approach is based on relating structural deformation to the applied stress field, after restoring mechanical equilibrium through energy minimization [60]. Theodorou and Suter [60] conducted a mechanical behavior study of atactic polypropylene at 233 K, considering 15 independent configurations; the authors found Young’s modulus to be within 10% of what was measured experimentally. Commercial software packages have facilitated the creation of various amorphous structures. Fan and Hsu [61] conducted a molecular characterization of aromatic polysulfone elastic moduli within an acceptable range of accuracy and uncertainty. Hutnik et al. [62] considered a more complex system of glassy polycarbonate of 4,4′-isopropylidenediphenol for an analysis of the elastic and plastic properties, observing a high degree of uncertainty attributable to the number of configurations considered. Fan et al. [63] studied a similar system and reduced the uncertainty by considering more configurations. Raaska et al. [64] calculated elastic constants, using both static and dynamic stress-strain concepts. The reported Young’s modulus value was comparable to that which was obtained experimentally. Kerogen, a set of naturally occurring materials consisting of carbon, hydrogen, oxygen, nitrogen, and sulfur and connected chemically through aliphatic, cyclic, and aromatic structures, was found to primarily be similar in chemistry to the aforementioned polymerized compounds. The success achieved in the calculation of the mechanical properties of polymers rendered an atomistic approach a reliable proxy, providing useful insights into kerogen geomechanics and the factors controlling it [65–70]. These studies shed some light on the mechanical properties of kerogen in a general sense. However, a detailed study of how kerogen heterogeneity influences the geomechanics of kerogen is not exhaustive. In this study, the impact of pore structure is linked to the mechanical behavior through an integrated study considering some range of kerogen configurations, temperature, and fluid effects.In the next section, the nanostructure construction of kerogen is detailed, followed by the molecular algorithms utilized to characterize the structure for its relevant petrophysical aspects. The mechanical properties simulation adopted in this research is also presented.


## 1.3. Realistic Kerogen Structures

To study kerogen’s geomechanics, a realistic structure epitomizing kerogen is required to serve as nanoporous media. In this research, a novel molecular modelling approach was followed to construct a nanoporous kerogen model virtually. The extraction of organic matter in sufficient quantities for a petrophysical assessment is arduous. The process generally requires the application of etching agents and/or thermal treatment. Additionally, the desired conditions of pressure and temperature at which different analyses are performed can be experimentally challenging. Attributed to its unprecedented flexibility, molecular modeling framework has the advantage of controlling the petrophysical properties of kerogen during the initiation phase where different kerogen types, bonding strengths, and other chemical properties can be changed in order to independently isolate their impact. The geomechanical assessments carried out in this study relied on the novel work of the kerogen models derived by Ungerer et al. [59]. These units are understood to cover the basic kerogen types (see Figure 1 and Table 1).Figure 1
(a) LAMMPS protocol workflow for creating kerogen nanostructures. (b) from a single kerogen unit to the final kerogen structure [82].
(a)(b)Table 1
Kerogen prototypes presented by Ungerer [58], reproduced from Alafnan [82].
Kerogen typeMaturity levelH/CO/CN/CS/CKIAImmature1.5300.0520.0280.012KIIAImmature1.1700.0950.0240.012KIIBTop of oil window1.1200.0600.0220.009KIICMiddle end of oil window0.9050.0540.0210.008KIIDOvermature0.5800.0510.0230.011KIIIAImmature0.8700.1160.0170.000In this study, kerogen type IIA was used as a representative of the organic matter in shale. Type IIA kerogen macromolecules were recreated on a molecular level, see Figure2 for a graphic representation of the various kerogen types. Polymer-consistent forcefield plus (PCFF+) was used to parameterize the atoms present in kerogen as center forces [71, 72]. Each center is defined and assigned a charge. Lennard Jones 6-9 was employed to describe the interactions between different center forces. The PCFF+ is very well known for its ability to configure organic compounds complementing earlier efforts of forcefields such as with COMPASS and CFF93 [71, 72].Figure 2
Van Krevelen diagram of kerogen maturation [59].An initial low-density structure containing eight units of kerogen IIA was created in relatively large cell volume to avoid any issues caused by larger sized kerogen macromolecules. The final kerogen structure was then formed utilizing a large-scale atomic/molecular massively parallel simulator LAMMPS. The simulation proceeded to a final condition of 20.67 MPa and 350 K which are deemed to represent typical reservoir conditions. The molecular dynamic steps began with initializing the cell by assigning 9.5 cutoff value and 2.0 skin with periodic boundaries, that was followed by a proper configuration of the molecular positions and velocity through energy minimization. Next, two stages of isochoric-isothermal NVT and isobaric-isothermal NPT for 250 ps and 200 ps, respectively, were run at 900 K. Then, the temperature was gradually reduced from 900 K to the final temperature through three consecutive NPT stages. The gradual decrease to the final targeted temperature assured faster convergence of kerogen units. The density of the final kerogen structure was around 1.1 g/cm3, which is close to the measured experimental one for the same type [59]. During the creation of the initial cell, a number of dummy molecules were added to act as a placeholder. The dummy molecules were removed when the final structure was formed, leaving behind a certain number of voids. The same process was repeated with an increasing number of dummy molecules, yielding a total of seven configurations across a range of density values (see Table 2). A description of the molecular construction protocol is given in Figure 1.Table 2
Seven kerogen configurations constructed.
ConfigurationNumber of dummy moleculesDensity (g/cm3)Porosity101.1146.45%251.0968.11%3101.0829.36%4151.05811.43%5201.04113.48%6251.02314.77%7350.96319.68%


## 1.4. Petrophysical Characterization

Helium–kerogen interactions in a Gibbs Monte Carlo simulation was used as a proxy for determining the void space, through the following relationship:(1)Nm=Na−ρaVP.Here, ρa is the density of the helium, Nm is the number of excess helium, and Na is the number of adsorbed molecules of helium. The pore volume VP can be calculated assuming no excess adsorption which is reasonably justified by the confined spaces in kerogen. Once the pore volume is estimated, the porosity could then be calculated when divided by the bulk volume. The calculations were performed on all kerogen structures, revealing increasing void ratios as the number of dummy molecules increased (see Table 1).Porosity, which is an indicator of the total pore volume, can be used to describe the total voids of porous media. A detailed distribution of pore sizes, however, is needed to fully characterize a kerogen structure. For this purpose, a pore size distribution (PSD) analysis was carried out for all seven nanoporous kerogen structures. The insertion of nonoverlapping spheres to fill the void spaces was used to determine the PSD distribution. The process involves predefining a threshold radius (i.e., a sphere was created in a certain location if its radius was greater than or equal to 0.15 nm in this study) [73]. An example of the PSD calculations performed can be found in Figure 3. The PSD and cumulative distribution function are given in Figures 4 and 5, respectively. The analysis revealed an increasing number of larger pores as the number of dummy molecules increased (i.e., larger pores were found to be positively correlated with porosity).Figure 3
Visual illustration of the PSD algorithm for one of the kerogen structures [82].Figure 4
PSD of the kerogen configurations.Figure 5
Cumulative distribution function of the PSD for the seven configurations.


## 1.5. Mechanical Properties

The kerogen’s molecular structure, initially at equilibrium, was deformed by changing the lattice vectorR=a,b,c to a disturbed cell R′=a′,b′,c′, using a strain matrix e [74–76].
(2)R′=R1+exx12exy12exz12eyx1+eyy12eyz12ezx12ezy1+ezz.The changes imposed led to a shift in the total energy of the system:(3)U=Etot−E0V0=12∑i=16∑j=16Cijeiej,whereU is the stress per unit volume, V0 is the initial volume of the cell, E0 is the energy of prior deformation, Etot is the energy after deformation, and Cij is the element of the stiffness matrix. Both i and j run from one to six {xx, yy, zz, yz, xz, xy}, yielding a stress tensor of 36 elements. Assuming symmetry, the total number of elements was reduced to three for a cubic cell.A predefined straine was defined across one direction (i.e., one component of the lattice). All other strain components were set to zero. The energy change U, defined by Equation (3), reduced as follows [77].
(4)U=C11e22.The above equation was used to calculate element C11. Similarly, the diagonal strains was set as eyz=ezy=1/2e, and all other strain components were set to zero. The energy change equation could then be reduced to calculate C44 [77].
(5)U=C44e22.The bulk modulus, an elastic coefficient representing the response to a uniform applied stress, was calculated by settingexx=eyy=ezz=e, while the shear modulus, a coefficient of the elastic response to shear stress, was estimated by setting ezz=e and exx=eyy=−1/2e [78]. For estimating the elastic moduli, smaller value of strains is used. Elastic moduli are obtained at two values of strains to assure consistency (i.e., 0.02 and 0.04 strains were selected).The aforementioned approach was applied to the different kerogen configurations. The strain implementation and energy calculations were conducted through LAMMPS and assisted by a MedeA interface. The molecular dynamic protocol consisted of initialization with a 3D periodic boundary, 9.5 cutoff distance, and 2.0 skin, followed by energy minimization to optimize the positions of the atoms. Then, strains were applied in predefined directions. The energy was subsequently minimized, and NVT ensemble was run at the specified temperature, allowing for a calculation of the total change in energy and stiffness matrix. The data were then used to estimate the elastic moduli. A summary of the molecular approach is given in Figure6.Figure 6
LAMMPS-based molecular approach used to calculate mechanical properties.The above concept was extended beyond the elastic region by continuously deforming the structure at an incremental strain and minimizing the energy before and after deformation to mimic tensile, compression, and shear failure tests.Kerogen is anticipated to exhibit both elastic and plastic behaviors during the shale stimulation phase and production span. In this study, elastic property and tensile failure tests were carried out to develop a comprehensive understanding of kerogen’s geomechanics. The details are presented in the next section.


## 2. Results and Discussion

### 2.1. Porosity Effect

The lattice dimension of the seven configurations ranged from 3.6 nm to 3.8 nm, corresponding to a bulk volume of 46.656 nm3 to 54.872 nm3 (i.e., the outcome of applying the construction protocol summarized in Figure 1). To assure that the selected number of kerogen macromolecules was sufficient to yield representative geomechanics, a pragmatic approach of checking the effect of increasing the number of kerogen units was followed. A kerogen structure built from two kerogen units was created, and its density was obtained. The same process was repeated at an increasing number of kerogen units. It was evident that the density converged to the same value for kerogen containing more than four units, while the structure consisting of two units deviated significantly (see Figure 7). The convergence in the density is consistent with estimated mechanical moduli as shown in Figure 8. Similar observations of the number of units needed to form a realistic kerogen structure are reported in the literature [79–84].Figure 7
Density convergence of kerogen structures built using an increasing number of kerogen units.Figure 8
Elastic moduli of kerogen structures, revealing an acceptable degree of reproducibility for kerogen constructed from more than four units.Next, the impact of the kerogen pore structure on its geomechanics was explored in the seven configurations, with porosity values ranging from 6.45% to 19.68% (see Table2). The analysis began with extraction of the elastic moduli. The results showed that kerogen’s elastic moduli decreased as porosity increased, an observation suggesting a higher probability of plastic deformation as the void ratio increases (see Figure 9). The elastic moduli were then utilized to calculate the Poisson ratio.
(6)ν=E2G−1,Figure 9
Elastic moduli as a function of porosity for kerogen containing eight units.whereE is Young’s modulus, G is the shear modulus, and ν is the Poisson ratio.The Poisson ratio is a coefficient used to characterize the tendency of materials to expand in a direction perpendicular to the direction of compression. In geomechanics, it is an indicator of a rock’s likelihood to fracture. For kerogen, the ratio was found to negatively correlate with porosity with some fluctuations at a lower range of porosity (see Figure10). The obtained results of elastic moduli in this computational study are consistent with the experimental work reported in the literature. The experimental approach primarily utilizes nanoindentation coupled with atomic force microscopy to extract Young’s modulus and Poisson’s ratio (see Table 3).Figure 10
Poisson’s ratio for kerogen as a function of porosity.Table 3
Experimental assessments of Young’s modulus and Poisson’s ratio of kerogen samples taken from different shale formations.
FormationPoisson’s ratioYoung’s modulus (GPa)Kumar et al., [85]Woodford shale (USA)0.36.2Kimmeridge (USA)0.35Eliyahu et al., [26]Upper Jurassic Source Rock (USA)0.3Up to 2.5Ahmadov et al., [86]Bazhenov (Russia)0.355.20.454.7Deformation and fracturing occur outside the elastic region. Thus, the stress-strain relationship detailed in Section2.3 was used to explore the mechanisms of mechanical failure. As described previously, an incremental strain was applied along the z-direction to the seven kerogen configurations, with a total strain of 0.4 (see Figure 11).Figure 11
Tensile mechanical deformation of kerogen structures with a porosity range of 6.45% to 19.68%. Three cases were shown to avoid overlapping.The kerogen with the lowest porosity exhibited a prolonged elastic region with an ultimate strength of 0.191 GPa. It showed early signs of failure at a 0.08 strain. The failure progressed through a series of strains, a behavior intrinsically different from what is typically observed for the inorganic constituents of sedimentary rocks. A complete fracture was not observed until a strain value of approximately 0.3. The other kerogen configurations showed similar trends, but with a decreasing order of ultimate strength as porosity increased. The kerogen configurations with the highest porosity experienced the earliest signs of failure, at a 0.05 strain.The geomechanical behavior of the different kerogen configurations intuitively matched the petrophysical analysis conducted in Section2.2. As porosity increased, larger pores were introduced into the structures. These pores tended to serve as points of stress accumulation, enhancing fracture initiation and propagation at smaller magnitudes of applied stresses. A visual evolution of tensile deformation for the two kerogen structures with the lowest and highest porosity values is given in Figure 12. It can be seen that the structure with higher porosity began to deform more severely at a lower strain, in alignment with the larger pores.Figure 12
Visualization of the evolution of tensile deformation for two kerogen structures. The first had a porosity of 6.45% and the second of 19.68%.


### 2.2. Temperature Effect

The mechanical behavior of oncrystalline materials such as kerogen is more vulnerable to be influenced by changes in the temperature. The mechanical assessment of kerogen in the previous section is conducted at 350 K, which is deemed to be a representative of average reservoir temperature. However, reservoir temperature may vary depending on the depth and the location of the reservoir. The adopted framework in this study is capable of addressing the effect of temperature as the construction of the structure, and the quantification of its energy is performed honoring the selected temperature. The stress-strain relationship revealed high sensitivity of the mechanical behavior on the temperature (see Figure13). As the temperature increases, kerogen seems to undergo prolonged plastic deformation before failure.Figure 13
Tensile mechanical deformation of 6.45% porosity kerogen structure at 300 K, 350 K, and 450 K.


### 2.3. Pore Fluid Effect

The reported data in the previous sections were made considering only kerogen with empty pores. However, in the reservoir, pore volume is occupied with fluids in quantities that are determined by the conditions of pressure, temperature, and the available pore volume. To study the effect of fluids, kerogen pore volume was saturated with natural gas (i.e. primarily represented by methane in this study) through running Gibbs Monte Carlo simulation. Two structures (highest and lowest porosity) were subjected to interactions with methane. The output is the number of methane molecules trapped in the structure. The calculations were performed using MedeA interface. Methane fugacity, which is needed for the calculations, was estimated using Peng-Robinson equation of state. In this study, Gibbs Monte Carlo simulation was performed at 6000 psi (i.e., the fugacity is 35.84 MPa at this pressure).Kerogen structures filled with methane were then assessed for their mechanical behavior (see Figures14 and 15). It is evident that the contained methane in pore volume altered the mechanical response of kerogen especially for the high porosity case (i.e., with high porosity, more methane molecules are occupying the structure).Figure 14
Tensile mechanical deformation of 6.45% porosity kerogen structure for empty and filled with methane cases at 6000 psi.Figure 15
Tensile mechanical deformation of 19.68% porosity kerogen structure for empty and filled with methane cases at 6000 psi.


## 2.1. Porosity Effect

The lattice dimension of the seven configurations ranged from 3.6 nm to 3.8 nm, corresponding to a bulk volume of 46.656 nm3 to 54.872 nm3 (i.e., the outcome of applying the construction protocol summarized in Figure 1). To assure that the selected number of kerogen macromolecules was sufficient to yield representative geomechanics, a pragmatic approach of checking the effect of increasing the number of kerogen units was followed. A kerogen structure built from two kerogen units was created, and its density was obtained. The same process was repeated at an increasing number of kerogen units. It was evident that the density converged to the same value for kerogen containing more than four units, while the structure consisting of two units deviated significantly (see Figure 7). The convergence in the density is consistent with estimated mechanical moduli as shown in Figure 8. Similar observations of the number of units needed to form a realistic kerogen structure are reported in the literature [79–84].Figure 7
Density convergence of kerogen structures built using an increasing number of kerogen units.Figure 8
Elastic moduli of kerogen structures, revealing an acceptable degree of reproducibility for kerogen constructed from more than four units.Next, the impact of the kerogen pore structure on its geomechanics was explored in the seven configurations, with porosity values ranging from 6.45% to 19.68% (see Table2). The analysis began with extraction of the elastic moduli. The results showed that kerogen’s elastic moduli decreased as porosity increased, an observation suggesting a higher probability of plastic deformation as the void ratio increases (see Figure 9). The elastic moduli were then utilized to calculate the Poisson ratio.
(6)ν=E2G−1,Figure 9
Elastic moduli as a function of porosity for kerogen containing eight units.whereE is Young’s modulus, G is the shear modulus, and ν is the Poisson ratio.The Poisson ratio is a coefficient used to characterize the tendency of materials to expand in a direction perpendicular to the direction of compression. In geomechanics, it is an indicator of a rock’s likelihood to fracture. For kerogen, the ratio was found to negatively correlate with porosity with some fluctuations at a lower range of porosity (see Figure10). The obtained results of elastic moduli in this computational study are consistent with the experimental work reported in the literature. The experimental approach primarily utilizes nanoindentation coupled with atomic force microscopy to extract Young’s modulus and Poisson’s ratio (see Table 3).Figure 10
Poisson’s ratio for kerogen as a function of porosity.Table 3
Experimental assessments of Young’s modulus and Poisson’s ratio of kerogen samples taken from different shale formations.
FormationPoisson’s ratioYoung’s modulus (GPa)Kumar et al., [85]Woodford shale (USA)0.36.2Kimmeridge (USA)0.35Eliyahu et al., [26]Upper Jurassic Source Rock (USA)0.3Up to 2.5Ahmadov et al., [86]Bazhenov (Russia)0.355.20.454.7Deformation and fracturing occur outside the elastic region. Thus, the stress-strain relationship detailed in Section2.3 was used to explore the mechanisms of mechanical failure. As described previously, an incremental strain was applied along the z-direction to the seven kerogen configurations, with a total strain of 0.4 (see Figure 11).Figure 11
Tensile mechanical deformation of kerogen structures with a porosity range of 6.45% to 19.68%. Three cases were shown to avoid overlapping.The kerogen with the lowest porosity exhibited a prolonged elastic region with an ultimate strength of 0.191 GPa. It showed early signs of failure at a 0.08 strain. The failure progressed through a series of strains, a behavior intrinsically different from what is typically observed for the inorganic constituents of sedimentary rocks. A complete fracture was not observed until a strain value of approximately 0.3. The other kerogen configurations showed similar trends, but with a decreasing order of ultimate strength as porosity increased. The kerogen configurations with the highest porosity experienced the earliest signs of failure, at a 0.05 strain.The geomechanical behavior of the different kerogen configurations intuitively matched the petrophysical analysis conducted in Section2.2. As porosity increased, larger pores were introduced into the structures. These pores tended to serve as points of stress accumulation, enhancing fracture initiation and propagation at smaller magnitudes of applied stresses. A visual evolution of tensile deformation for the two kerogen structures with the lowest and highest porosity values is given in Figure 12. It can be seen that the structure with higher porosity began to deform more severely at a lower strain, in alignment with the larger pores.Figure 12
Visualization of the evolution of tensile deformation for two kerogen structures. The first had a porosity of 6.45% and the second of 19.68%.


## 2.2. Temperature Effect

The mechanical behavior of oncrystalline materials such as kerogen is more vulnerable to be influenced by changes in the temperature. The mechanical assessment of kerogen in the previous section is conducted at 350 K, which is deemed to be a representative of average reservoir temperature. However, reservoir temperature may vary depending on the depth and the location of the reservoir. The adopted framework in this study is capable of addressing the effect of temperature as the construction of the structure, and the quantification of its energy is performed honoring the selected temperature. The stress-strain relationship revealed high sensitivity of the mechanical behavior on the temperature (see Figure13). As the temperature increases, kerogen seems to undergo prolonged plastic deformation before failure.Figure 13
Tensile mechanical deformation of 6.45% porosity kerogen structure at 300 K, 350 K, and 450 K.


## 2.3. Pore Fluid Effect

The reported data in the previous sections were made considering only kerogen with empty pores. However, in the reservoir, pore volume is occupied with fluids in quantities that are determined by the conditions of pressure, temperature, and the available pore volume. To study the effect of fluids, kerogen pore volume was saturated with natural gas (i.e. primarily represented by methane in this study) through running Gibbs Monte Carlo simulation. Two structures (highest and lowest porosity) were subjected to interactions with methane. The output is the number of methane molecules trapped in the structure. The calculations were performed using MedeA interface. Methane fugacity, which is needed for the calculations, was estimated using Peng-Robinson equation of state. In this study, Gibbs Monte Carlo simulation was performed at 6000 psi (i.e., the fugacity is 35.84 MPa at this pressure).Kerogen structures filled with methane were then assessed for their mechanical behavior (see Figures14 and 15). It is evident that the contained methane in pore volume altered the mechanical response of kerogen especially for the high porosity case (i.e., with high porosity, more methane molecules are occupying the structure).Figure 14
Tensile mechanical deformation of 6.45% porosity kerogen structure for empty and filled with methane cases at 6000 psi.Figure 15
Tensile mechanical deformation of 19.68% porosity kerogen structure for empty and filled with methane cases at 6000 psi.


## 3. Conclusion

The geomechanical behavior of shale is governed by its individual constituents, including kerogen. The presence of kerogen as finely dispersed fragments renders it inaccessible to traditional mechanical testing methods. In this study, a novel molecular modeling approach was followed to construct, characterize, and mechanically test kerogen configurations across a range of petrophysical properties, in order to determine their geomechanical characteristics. The results revealed that the kerogen exhibited both elastic and plastic deformations. The kerogen with the lowest porosity had the maximum ultimate strength and highest elastic moduli. The higher porosity kerogen showed a higher probability of failure, which was attributable to the presence of larger pores. Kerogen mechanical behavior was found also sensitive to the conditions of temperature and the fluid occupying the pore volume. The findings of this work can be used in a larger scale computational modeling of shale’s geomechanics. It is recommended that this work be validated and the findings of this study experimentally correlated.



---

*Source: 4093895-2021-09-15.xml*