CH4

Article Cite This: Energy Fuels XXXX, XXX, XXX−XXX

pubs.acs.org/EF

Molecular Insights into Kerogen Deformation Induced by CO2/CH4 Sorption: Effect of Maturity and Moisture Liang Huang,*,†,‡ Zhengfu Ning,*,† Qing Wang,† Rongrong Qi,† Zhilin Cheng,† Xiaojun Wu,† Wentong Zhang,† and Huibo Qin§ †

State Key Laboratory of Petroleum Resources and Prospecting and §State Key Laboratory of Heavy Oil Processing, China University of Petroleum (Beijing), Beijing 102249, People’s Republic of China ‡ Department of Chemical and Biomolecular Engineering, University of California, Berkeley, Berkeley, California 94720-1462, United States

Downloaded by UNIV OF SOUTHERN INDIANA at 03:44:35:617 on June 12, 2019 from https://pubs.acs.org/doi/10.1021/acs.energyfuels.9b00409.

S Supporting Information *

ABSTRACT: Kerogen is the source of hydrocarbons in shale. As a soft nanoporous matter, kerogen constantly experiences mechanical deformation induced by subsurface stress environment and interplay with geofluid, significantly affecting the percolation and production of hydrocarbons. However, the associated coupling deformation of kerogen in reservoir conditions remains poorly understood. Here we quantify the kerogen coupling deformation induced by moisture uptake, CH4 and CO2 sorption, and mechanical compression using the combination of molecular simulations and poromechanics theory. The effects of kerogen maturity and preloaded moisture are discussed in detail. Our results show that moisture induced deformation is governed by the dynamic distribution of water molecules and the physicochemical characteristics of kerogen. Moisture induced volumetric strain increases with decreasing kerogen maturity or rising moisture content. The porosity of kerogen skeleton increases upon moisture uptake, although the residual porosity accessible for fluid presents a linearly decreasing trend. The coupling deformation upon gas sorption results from the combined effect of sorption swelling and mechanical compression. The coupling volumetric strain increases with rising kerogen maturity or declining moisture content, contrary to results for moisture uptake. The total deformation, coupling sorption−pressure−moisture effect, rises with increasing kerogen maturity. The evolution the total deformation follows with moisture content depends on the kerogen maturity. Results of this work can help improve the estimation of fluid-in-place in a shale gas reservoir and the understanding of the chemomechanical coupling associated with CO2 sequestration and CO2 fracturing.

1. INTRODUCTION Shale gas is considered as one important substitute for conventional hydrocarbon resources due to its considerable resource abundance, high utilization efficiency, and environmental friendliness.1−3 The share of shale gas in the total production of natural gas in the United States is estimated to be more than 45% in 2035.4 Despite the sharp expansion in production, the recovery factor is still less than 20% for most shale gas reservoirs.5 In addition to depleted development regime, the low recovery can be also attributed to the presence of amorphous kerogen interacting with geofluids, which results in increasing complexity for shale gas production.6−8 Kerogen is a complex heterogeneous matter insoluble in organic solvents, which represents the major structure characteristics of organic matter. It plays a key role in the generation of oil and gas. Also, kerogen has abundant nanopores and large specific surface areas, which makes it the main storage nodules for shale gas.9 The trapped hydrocarbons in kerogen can take up almost half of the total share in gas shale.10 The innate properties of kerogen, including atomic compositions, molecular structure, functional groups, and nanopores, are governed by the organic type and kerogen maturity. The organic type is correlated with the biological source and the depositional environment, while the maturity evolves with the thermal maturation. In addition to physical heterogeneity associated with porous network, © XXXX American Chemical Society

kerogen also possesses chemical heterogeneity related to functional groups,11 leading to the characteristic of mix-wet. Although the kerogen skeleton is widely deemed as hydrophobic,12,13 the carboxylic and hydroxylic functional groups in kerogen are hydrophilic, showing strong affinity with moisture.14,15 Moisture distributed in kerogen can not only occupy gas sorption sites, but also change the interaction between kerogen and gas. The dynamic interaction of various kerogen matrixes and moisture can provide a better understanding of the evaluation of original gas reserves and the fate of residual fracturing fluids.16,17 Extensive efforts have been devoted to the physicochemical properties of kerogen and its interaction with geofluids. Kelemen et al.18 derived the chemical structure parameters of kerogens with different organic types and maturities by using analytical experiments such as solid-state 13C NMR spectroscopy, X-ray photoelectron spectroscopy, and sulfur Xray absorption near edge structure (S-XANES) technique. Based on Kelemen’s work, Ungerer et al.19 constructed the chemical structures of kerogens at atomistic scale and extended these kerogen units to three-dimensional models with molecular modeling. Similarly, we have built the molecular Received: February 8, 2019 Revised: May 29, 2019

A

DOI: 10.1021/acs.energyfuels.9b00409 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

kerogen used for gas adsorption in experiments is subject to deformation, the characterization of deformation remains a big challenge due to the stiffness and small volumetric strain of kerogen sample, as well as the difficulty in kerogen extraction from shale. To date, there is no report on kerogen deformation by experimental measurement. Heller and Zoback38 measured the volumetric strain of activated carbon induced by gas adsorption, but only obtained a strain on the magnitude of 10−3. Alternatively, some authors have devoted theoretical efforts to unify fluid sorption and mechanical deformation within porous framework.39−41 Brochard et al.40 developed a poromechanical model to assess the sorption-induced deformation in solids extended to micropores. This model has successfully described the interplay between sorption and deformation in microporous solids such as coal,42 silicon,43 and kerogen.31 Despite previous efforts in the literature, the kerogen deformation associated with gas sorption remains poorly understood. To the best of our knowledge, the effect of preloaded moisture on kerogen deformation induced by gas sorption in kerogens with different maturities has not yet been reported. In this work, we adopt the combination of molecular simulations and poromechanics theory to quantify the kerogen coupling deformation induced by gas sorption considering the effect of preloaded moisture and kerogen maturity. This is the first work to gain deep insights into this topic. The paper is organized as follows. In section 2, we first present the computational methodology including kerogen units, simulation details, and computation theories. In section 3.1, we validate the generated kerogen models by comparing the chemical structure, density, porosity, and excess adsorption with experimental results. In section 3.2, the effects of preloaded moisture and kerogen maturity on kerogen deformation accompanied by changes in nanostructure are discussed based on MD simulations. In section 3.3, we first show the gas sorption isotherms from GCMC simulations, and then adopt the poromechanical model to analyze the coupling deformation of kerogen induced by combined sorption− pressure effect. Next, we explore the total deformation of kerogen taking into account the coupled sorption−pressure− moisture effect. Finally, we draw some conclusions in section 4.

models of kerogens extracted from the Silurian Shale and Cambrian Shale in the Chinese Sichuan Basin by direct experimental characterizations and molecular modeling.20,21 Structural and thermodynamic properties of these kerogen models have also been studied in the work of Ungerer et al.,19 as well as our previous work,20,21 by molecular simulations. Most reported kerogen matrix models in the literature were built following the produce proposed by Ungerer et al.,19 in which the size of nanopores was generally smaller than 10 Å. Some authors have put efforts into improving the porous network of kerogen matrixes at atomistic scale. Tesson and Firoozabadi inserted dummy particles into a kerogen molecular system to expand the pore size.22 Zhou et al.23 inserted cutter clusters to improve both the pore size and heterogeneous pore surface in a kerogen molecular model. In our previous work,20,21 we have attempted to control the micropores in kerogen models by inserting and deleting organic species and gas molecules so as to simulate the realistic paths in reservoir conditions. Kerogen interaction with fluid molecules has been widely investigated by experimental measurements24−26 and molecular simulations.27−34 In reported experimental measurements about CO2 and lighter hydrocarbons adsorption in isolated kerogen,24−26 kerogen samples are allowed to deform upon gas adsorption, which is consistent with the evolutions in realistic reservoir conditions. Kerogen is one kind of soft nanoporous matter exposed to complex stress constraints in subsurface environment. There is a strong coupling between gas sorption and kerogen deformation accompanied by changes in kerogen microstructure. However, most documented work by molecular simulations was based on kerogen models with fixed structures,28−34 ignoring the sorption-induced kerogen deformation. Recently, some investigations have been performed to take into account the flexibility of kerogen upon gas sorption with molecular simulation. Pathak et al.35 obtained the swelling of flexible kerogen by interacting with CH4 and CO2 sorption using molecular dynamic (MD) simulations. In their work, the initially separated systems of kerogen and fixed gas molecules were annealed together to simulate the sorption process at certain pressures, which fail to represent the dynamic deformation of kerogen induced by gas sorption. Tesson and Firoozabadi22 studied CH4 sorption in a flexible kerogen model using a hybrid grand canonical Monte Carlo (GCMC) method based on MD simulation (MD-GCMC). They found that the adsorption in flexible kerogen can be 57% higher than that in rigid structure. The MD in their hybrid MD-GCMC simulation was performed in the NVT ensemble (constant molecular number, volume, and temperature). Although kerogen was flexible, the volume of the kerogen model was fixed in their work. Ho et al.36 also adopted a hybrid MDGCMC method but with NPT ensemble (constant molecular number, pressure, and temperature) for MD to simulate the deformation of kerogen induced by CH4 and CO2 sorption. They observed that kerogen swelling can increase the gas adsorption. Their simulations were conducted at zero effective stress condition to observe the maximum swelling of kerogen, mismatching the stress paths in reservoir conditions. By contrast, Vasileiadis et al.37 compared the adsorption of CH4 in rigid and flexible kerogen structures using GCMC simulations. They concluded that the flexibility of kerogen may not affect CH4 adsorption. Therefore, the kerogen deformation associated with gas sorption using molecular simulations is still in a preliminary stage and remains to be further verified. Although

2. COMPUTATIONAL METHODOLOGY 2.1. Kerogen Structure Units. In this work, we adopt three different mature kerogen units (IIA, IIC, and IID in Figure A.1 in the Supporting Information), covering the maturity from immature to postmature. The kerogen units with type II organic type correspond to the typical sources of conventional and unconventional liquid and gas resources, derived from the organic-rich Duvernay shales. The immature kerogen IIA is a typical source of conventional hydrocarbon resources. The kerogen IIC can represent mature organic matter in a shale reservoir with liquid hydrocarbons, like that in the Bakken shale.44 The kerogen IID is representative of postmature organic matter deposited in a shale gas reservoir, similar to that in the Barnett shale.45 The kerogen units were built by Ungerer et al.19 based on gathered experimental results of 13C NMR spectroscopy, X-ray photoelectron spectroscopy, and S-XANES data.18 The breakdowns of atomic compositions and chemical structural parameters of the kerogen units are listed in Table A.1 in the Supporting Information, which are in good agreement with the analytical results.18 B

DOI: 10.1021/acs.energyfuels.9b00409 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels 2.2. Construction of Kerogen Models. In this work, geometry optimization method and MD simulations were combined to generate the three-dimensional kerogen models with different maturities. The kerogen units were first relaxed with the geometry optimization and annealing dynamics. Geometry optimization is expected to find the relaxed structures with local minimal energy. The smart minimization algorithm with a fine convergence criterion was adopted to optimize the geometry structure of kerogen units. The cutoff distance of the short-range nonbonded interaction was set as 15.5 Å. The annealing dynamics aim at finding the stable structure with global minimal energy. A succession of MD simulations with the NVT ensemble was carried out in each cycle of annealing dynamics simulations. A total of 10 annealing cycles with temperature changing between 300 and 800 K and a total simulation time of 400 ps were executed to obtain the stable structures with the lowest energy. A number of relaxed kerogen units were subsequently incorporated into an amorphous cell to prepare the initial configurations of bulk kerogen models. The initial density of the amorphous cell is targeted at 0.1 g/cm3.19 A relaxation process including a succession of MD simulations19 was thereafter adopted to bring the kerogen models to the equilibrated state. We first performed a structure relaxation with NVT ensemble on these initial configurations at 800 K for 400 ps. Then succeeding MD simulations with NPT ensemble were conducted at 20 MPa with a stepwise decreasing temperature from 800 to 300 K.46 A sufficient simulation time of 400 ps with a time step of 1 fs is adopted for each NPT run. We determined the simulation time by checking the convergence of the density of kerogen models. Furthermore, these configurations of kerogen models were relaxed at 300 K and 0.1 MPa for 2000 ps to obtain the stable structures. 2.3. Computational Details. The kerogen models with various moisture contents were built using the fixed loading task in the Sorption module of Materials Studio. The range of moisture contents (0.7, 1.4, 2.1, and 2.8 wt %) for kerogen is determined so as to be consistent with documented experimental investigations47 and simulation work.33,48 MD simulation with NPT ensemble was then carried out to study the effect of moisture on kerogen deformation. The position of kerogen atoms and the size of simulation box are allowed to change during this relaxation process. The simulation time for each run is 5000 ps with a time step of 1 fs to guarantee the equilibrium of box size. The deformation of kerogen is quantified with the definition of volumetric strain: Vε =

V − V0 V0

respect to box volumes, was computed. The probe insertion method was adopted to detect the free pore volumes. For this geometric determination method, the randomly inserted probe molecule rolls over the van der Waals surfaces of the kerogen atoms. The edge of the probe traces out the Connolly surface,49 while its center traces the solvent surface. In this work, the free pore volumes are identified by the Connolly surfaces rather than the solvent surfaces to represent the realistic physical spaces between kerogen skeletons. The radius of helium (1.3 Å) was assigned for the probe molecule to compute the porosity so as to be consistent with the experiments. Moreover, we can get a number of decreasing free pore volumes by gradually increasing the radius of the probe molecule. Accordingly, the pore size distributions of kerogen models can be computed by differentiating the free pore volumes with respect to different widths of the probe molecule. The deformation of kerogen models induced by gas sorption was investigated using an extended poromechanical model proposed by Brochard et al.40 This extended model takes into account the effect of both free phase and adsorption phase of fluid on the mechanical behaviors of microporous structure. The impact of sorption is represented by adding an integration term on the usual linear poroelastic model: Vε = −

Vw − Vd Vw

∫0

f

nRT df f

(3)

where P is the bulk pressure, K is the bulk modulus, C is the fluid coupling coefficient, n is the total sorption amount of fluid in kerogen models at zero strain, f is fluid fugacity, T is temperature, and R is the universal constant. In this work, the values of K = 2.65 GPa and C = 6.05 for CH4 and 7.06 for CO2 were adopted, as documented in the work of Ottiger et al.50 and Brochard et al.42 To evaluate the influence of gas sorption on kerogen deformation with eq 3, a quantitative estimate of sorption amount is required. The sorption isotherms of CH4 and CO2 were simulated with the GCMC method. The kerogen structures were kept fixed during the sorption process to match the condition of zero strain, as required by the extended poromechanical model. For each pressure point, a total of 6 × 107 Monte Carlo steps are performed, wherein the first 3 × 107 steps are utilized to guarantee the adsorption equilibration, while the remaining steps are utilized to compute the statistical adsorption amount. The GCMC simulation takes into account the attempts for exchanges (insertion and deletion) and MC moves (translation and rotation) of sorbate molecules. The average proportion of attempts for sorbate insertion, deletion, translation, and rotation is 25, 25, 25, and 25%, respectively. The acceptance rate for both translation and rotation is maintained at approximately 50% in the simulation. The method of GCMC simulation outputs the total sorption amount as a function of fugacity rather than pressure. The fluid fugacity is converted from pressure through the Peng− Robinson equation.51 To compare the simulated results with experimental data, excess adsorption is needed, which is converted from the total sorption by ne = n − υρ (4)

(1)

where Vε is the volumetric strain; V and V0 are the box volumes of kerogen models with and without moisture, respectively. Also, the dissolution factor is defined to represent the proportion of loaded water volumes that do not contribute to kerogen deformation, which is expressed as β=

P C + K K

(2)

where β is the dissolution factor, Vw is the volume of loaded water molecules, and Vd is the volume of kerogen deformation due to loading of water molecules. To analyze the deformation of kerogen nanostructure, the porosity, defined as the fraction of free pore volumes with

where ne is the excess adsorption amount, ρ is the bulk gas density, and υ is the free pore volume of kerogen model. All the simulations in this work were performed using the Materials Studio package.52 The COMPASS force field53 was C

DOI: 10.1021/acs.energyfuels.9b00409 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

compositions and functional groups (Table A.1 in the Supporting Information) match fairly well the experimental results. These kerogen units were condensed into threedimensional realistic models through a rigorous relaxation process. The physical densities of our kerogen models at ambient conditions are close to those of documented experimental data. The immature kerogen model (IIA) has a density of 1.13 ± 0.02 g/cm3, which is in good agreement with that (1.0−1.15 g/cm3) of immature type II kerogen in New Albany Shale.57 The density (1.17 ± 0.02 g/cm3) of the mature kerogen model (IIC) is close to the type II kerogen densities (1.18−1.25 g/cm3) at a similar maturity.58 The density of the postmature kerogen model (IID) is 1.19 ± 0.02 g/cm3, following the reported trend that kerogen density increases with maturity.19,33,58 The kerogen porosities probed with helium are also comparable to documented experimental results. The porosity (10.9 ± 0.5%) of kerogen IIA is fairly accordant with the porosity (12.5 and 13.8%) for low mature type I−II kerogen of Posidonia Shales by helium pycnometry.59 The porosities of kerogen IIC and IID are 15.3 ± 0.9 and 25.7 ± 1.1%, respectively. The porosities of our kerogen models are coincidental with the reported range (4.45− 22.50%) in organic matter of Barnett mudstone.60 The converted excess adsorption isotherms were compared with the reported experimental data14,26,59,61,62 (Figure A.2 in the Supporting Information). The objective of this comparison is to check the reasonability rather than the accuracy for our kerogen models to match with documented experiments. The magnitude of simulated gas excess adsorption is close to the measured data. The trend that gas excess adsorption first increases to a maximum and then shows a decreasing trend with the increase of pressure is in agreement with the experimental results. Also, the amount of excess adsorption rises with the increase of kerogen maturity, which is consistent with the documented conclusion in the work of Zhang et al.63 by measurement. The difference between simulation and experiment can be mainly attributed to the discrepancy in pore accessibility,64,65 microstructure deformation,66 and sample properties. In Figure A.2, the simulated excess adsorption isotherms in kerogen models with deformation are also compared with the results in rigid kerogen models as documented in our previous work.32 The simulation results without considering the flexibility of kerogen are underestimated compared with the results in deformed kerogen models. The observation is consistent with recent conclusion in the literature,22,36 which further confirms that the flexibility of kerogen upon gas adsorption should be taken into account to accurately estimate the gas reserve in place. The excess adsorption is mainly controlled by the density difference between the adsorption phase and bulk phase according to the definition in eq 4. The density difference increases initially, and then decreases with further increasing pressure. The maximum excess adsorption occurs when the density difference reaches the maximum. If the density of adsorbed phase becomes less than that of bulk phase, a negative excess adsorption can be expected. Besides, the excess isotherm of CO2 reaches the maximum at a lower pressure (∼5.8 MPa) compared with that of CH4 (∼7.8 MPa). This can also be explained by the different change of density difference.26 To show the effect of random error, the standard deviation on excess adsorption is presented in Figure A.2. The uncertainty of excess adsorption for kerogen IIA can be neglected, since the error bars are smaller than the symbols.

adopted to describe the nonbonded dispersion−repulsion interactions between atomic pairs and the intramolecular interactions in kerogen units. For the COMPASS force field, the van der Waals interaction is computed via the condensedphase optimized Lennard-Jones 6−9 potential, while the electrostatic interaction is described with the unscreened Coulomb potential. Various cross-interaction terms are introduced for the intramolecular interactions to improve its performance. The nonbonded interaction is described as52,53 ÄÅ É 6Ñ ÅÅ i 0 y9 ij rij 0 yz ÑÑÑÑ qiqj Å r j z Å ij u(rij) = u LJ + uC = ∑ EijÅÅÅÅ2jjjj zzzz − 3jjjj zzzz ÑÑÑÑ + ∑ j rij z ÑÑ ÅÅ j rij z rij i>j i>j ÅÅÇ k { k { ÑÑÖ (5)

where u(rij) is the nonbonded interaction, uLJ is the van der Waals interaction, uC is the Coulomb interaction, rij is the distance between atomic sites i and j, qi and qj are the partial charges at atomic sites i and j, and Eij and rij0 are the energetic and size parameters between atomic sites i and j, respectively. The combining LJ parameters Eij and rij0 are derived based on the sixth power mixing rules:52,53 ij (ri 0)6 + (rj 0)6 yz j zz rij = jjj zz j z 2 k {

1/6

0

0 3 0 3 jij (ri ) (rj ) zyz Eij = 2 EiEj jjj 0 6 0 6 zzz j (r ) (r ) z j k i {

(6)

(7)

This force field has been widely used for simulation in kerogen systems.28,30−32,34 It has been proven to provide fairly consistent predictions for adsorption behaviors, microscopic structures and thermodynamic properties of kerogen.28,30−32,34 The cutoff radius for short-range pair interaction is set as 15.5 Å, which is smaller than the half-length of the minimum dimension of the simulation box. The barostat and thermostat methods are Berendsen54 and Nosé−Hoover,55 respectively. The long-range Coulombic interactions are summed by the Ewald method. The fluid molecules including H2O, CH4, and CO2 are represented by rigid full-atom models with point charge on each Lennard-Jones center.56 The source of uncertainty in molecular simulation is mainly derived from the statistical error, which is controlled by the simulation time for MD simulation and the sampling step for GCMC simulation. To reduce the statistical error, sufficient simulation time and sampling step were performed by referring to documented simulation details.19,31,33,46 We checked the system density, sorption amount, pressure, and temperature to guarantee the convergence in the equilibrium stage of MD and GCMC simulations. For MD simulation, the production stage was equally divided into five blocks. For GCMC simulation, five different kerogen configurations were used for sorption. The standard deviations among these blocks for MD simulation or configurations for GCMC simulation were computed to represent the uncertainty, as indicated by error bars. The simulated properties were averaged based on the five blocks or five configurations to show the statistical results.

3. RESULTS AND DISCUSSION 3.1. Validation of Kerogen Models and Methodology. The chemical structures of kerogen units used in this work were built based on analytical experiments. The atomic D

DOI: 10.1021/acs.energyfuels.9b00409 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels Although there are observable deviations for kerogens IIC and IID, the excess isotherms still present a decreasing trend at high pressure considering the deviations. In this work, we study the deformation of kerogen upon fluid uptake using the isotropic volume expansion/contraction assumption. This isotropic deformation assumption has been widely used in recent work in the literature.22,35,36 However, it should be noted that the realistic deformation of subsurface kerogen interacting with fluids may be anisotropic considering the big difference in isothermal compressibility between kerogen and fluids. To justify the isotropic assumption and analyze its possible effect on kerogen deformation, simulation tests for a kerogen/fluid system with isotropic deformation and anisotropic deformation were performed. The MD simulation with NPT ensemble is performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).67 The comparison for cell length and volume between the isotropic and anisotropic scenarios is shown in Figure A.3 in the Supporting Information. Figure A.3a,c shows that the cell lengths along different directions for the anisotropic scenario differentiate from that for the isotropic scenario. The composite system for the anisotropic scenario is observed to have a larger swelling along the y direction. Figure A.3b,d shows that the difference in the volume of the composite system between the isotropic scenario and the anisotropic scenario is small. The volumetric strains for the kerogen/CH4 system are 2.47 and 2.35% for the isotropic and anisotropic scenarios, respectively. The relative error is 5.02% between the two scenarios. The volumetric strains for the kerogen/CO2 system are 5.21 and 4.94% for the isotropic and anisotropic scenarios, respectively. The relative error is 5.26% between the two scenarios. Therefore, we can obtain different linear swelling strains along the three directions using the anisotropic scenario. However, the volumetric swelling strain from this scenario is fairly close to the results from the isotropic scenario. It should be noted that since we are focused on the volumetric swelling of kerogen upon fluid uptake in this work, we believe it is reasonable to adopt this isotropic deformation approximation. It is also worth noting that the importance of anisotropic physics for the kerogen/fluid composite system should not be underestimated, especially for linear deformation analysis. Further investigations are needed to gain more insights into this anisotropic physics of a composite system. 3.2. Effect of Moisture and Maturity on Kerogen Deformation. 3.2.1. Moisture Induced Deformation. Figure 1 shows the volumetric strain of kerogen as a function of moisture content. The volumetric strain is observed to increase with the moisture content, but follow different relationships for different mature kerogen models. Kerogen IIA presents a linear relationship, while kerogens IIC and IID show a quadratic increase of volumetric strain with respect to moisture content. The discrepancy is corelated with the size of kerogen pore. As discussed in our previous work,32 both the porosity and pore size increase with the kerogen maturity. For immature kerogen model with smaller pores (IIA), pore filling is the dominant mechanism for water loading. The moisture induced deformation is governed by the volumes of water molecules loaded into the system. This linear relationship has also been observed for water sorption in flexible nanoporous polymer by molecular simulation.68 The pore size of this polymer ranges between 0.2 and 0.4 nm,68 similar to that of kerogen IIA. For mature and postmature kerogen models (IIC and IID), there are a larger proportion of bigger pores in the system. Water

Figure 1. Volumetric strain of kerogen as a function of moisture content at 30 MPa and 338 K. IID_conden. correspond to kerogen deformation in scenario a, while IIA, IIC, and IID corresponds to kerogen deformation in scenario b. The error bars are estimated using the block averaging method.

loading is governed by the adsorption onto the pore surface. Ho et al.36 has observed the same quadratic relationship between volumetric strain and gas sorption amount on the postmature kerogen IID model through molecular simulation. They also theoretically proved this relationship based on a simple geometric consideration. Figure 1 shows that the volumetric strain decreases as the kerogen maturity rises. Kerogen swelling induced by moisture is insignificant for kerogen IID, and the volumetric strain is only 0.18 ± 0.01% at the moisture content of 2.8 wt %. At the same condition, the volumetric strain can be 2.51 ± 0.20 and 2.78 ± 0.14% for kerogens IIC and IIA, respectively. Although little has been reported on moisture induced kerogen deformation, there are a couple of studies on coal swelling induced by adsorption. The trend of kerogen swelling with respect to maturity is in good agreement with actual measurement on coal. Fry et al.69 concluded by measurement that coal swelling induced by moisture sorption depended on the rank of coal with lower rank coals swelling the most. Day et al.70 experimentally found that the maximum swelling induced by CH4 adsorption in the lowest coal was more than 2 times higher than that in the highest rank sample. Also, the absolute values of volumetric strain on our kerogen models match fairly well with the documented experiments on coal swelling induced by moisture adsorption. Our simulated maximum volumetric strain ranges between 0.18 ± 0.01 and 2.78 ± 0.14% among the three kerogen models. Pan et al.71 measured a maximum volumetric strain of about 1.5% for moisture induced swelling on coal. Fry et al.69 reported a range of 0.5− 5% for volumetric strain induced by moisture adsorption on coal. The trend that moisture induced kerogen swelling weakens with increasing maturity can be attributed to the dynamic distribution of water molecules and the physicochemical characteristics of kerogen models. As concluded in our previous study,28 at low moisture conditions, water molecules tend to occupy oxygen-containing groups in kerogen, while at high moisture conditions water molecules tend to migrate and aggregate into clusters in the middle of enterable pores. Lowmature kerogen contains more abundant oxygen-containing groups and smaller pore sizes, while high-mature kerogen possesses fewer oxygen-containing groups and bigger pore sizes. Therefore, water molecules tend to scatter around polar groups in low-mature kerogen causing a bigger swelling, and E

DOI: 10.1021/acs.energyfuels.9b00409 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 2. Snapshots of water molecule distribution in kerogen IID model at 30 MPa and 338 K (moisture content is 2.8 wt %). Atom representations: black, carbon; red, oxygen; yellow, sulfur; blue, nitrogen. Water molecules are represented by the red and white balls. (a) Scenario a. (b) Scenario b.

Figure 3. (a) Porosity as a function of moisture content. (b) Dissolution factor as a function of moisture content. (c) Porosity differential variation as a function of moisture content. The error bars represent the standard deviations of five kerogen configurations. For the sake of clarity, we only present the uncertainties for kerogen IIC in (c).

locate in the middle of big pores in high-mature kerogen leading to a smaller swelling. In this work, two scenarios were designed to study the effect of moisture loading order on kerogen deformation. In the first scenario, scenario a, kerogen and water molecules were initially placed into the system simultaneously. The system of kerogen and water molecules was then relaxed through the geometry optimization and annealing process. In scenario b, the system of kerogen molecules was first equilibrated, after which water molecules were loaded into the system to study moisture sorption-induced deformation. In Figure 1, we can observe that

kerogen deformation in scenario a is much more significant than that in scenario b. At the moisture content of 2.8 wt %, values of the volumetric strain of kerogen IID are 6.49 ± 0.26 and 0.18 ± 0.01% for scenarios a and b, respectively. The large discrepancy can be attributed to the distribution of water molecules in the system. The configurations of the system after equilibrium for scenarios a and b are shown in Figure 2. Figure 2a shows that water molecules aggregate into a big cluster in the kerogen structure for scenario a. This cluster can expand the kerogen matrix to occupy larger volumes. By contrast, in scenario b (Figure 2b), water molecules are observed to scatter F

DOI: 10.1021/acs.energyfuels.9b00409 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels in the porous skeleton. Consistent with our previous discussion, a small portion of water molecules are located around the heteroatom groups, while dominant water molecules are distributed in the middle of pores in this postmature kerogen model, which results in a smaller structure deformation. Note that scenario a represents the condition in the formation of kerogen, since the sources of kerogen are deposited in water environment. Scenario b represents the conditions in most laboratory measurements for moisture sorption. The real kerogen deformation mode with respect to water in subsurface conditions may be between the two. In this work, kerogen deformation is referred to scenario b without a specific statement, since our focus is on the kerogen deformation induced by moisture and gas sorption. The moisture induced swelling presents a small divergence under different configurations using the block averaging method. The largest uncertainty was observed to be 7.96% for kerogen IIC at the moisture content at 2.8 wt %. 3.2.2. Kerogen Nanostructure Deformation. Associated with the volumetric strain, moisture loaded into the kerogen can induce the deformation of nanostructure. In Figure 3a, we report the kerogen porosity as a function of moisture content at 30 MPa and 338 K. The porosity drops as the kerogen maturity increases among the moisture contents, in agreement with our discussion in section 3.1. Upon moisture uptake, the porosity presents a linearly decreasing trend for the three kerogen models. The decreasing amplitude of porosity among the studied range of moisture contents is in the order kerogen IIA (1.21%) < IIC (2.18%) < IID (3.06%). This can be well explained by the distribution of water molecules in different mature kerogen models.28 Water molecules are inclined to adsorb on polar groups in low-mature kerogen, resulting in smaller pore volume shrinkage, and occupy bigger pores in high-mature kerogen, leading to bigger porosity variation. This also means that not all of the volume of water molecules loaded into the system can contribute to kerogen deformation. As discussed above, water molecules near the kerogen atoms have a bigger effect on kerogen deformation, while water molecules dissolved in the middle of large pores exert a smaller effect. The porosity shows an accepted uncertainty with the standard error smaller than 1.16%. Kerogen IIC presents a larger uncertainty on porosity than kerogens IIA and IID. Figure 3b presents the dissolution factor as a function of moisture content. According to the definition of dissolution factor in section 2.3, it refers to the proportion of water volumes that do not contribute to kerogen deformation. In Figure 3b, the dissolution factor shows a decreasing trend with increasing kerogen maturity. At the moisture content of 0.7 wt %, the dissolution factor is 32.4 ± 1.6, 60.3 ± 3.6, and 94.4 ± 2.8% for kerogens IIA, IIC, and IID, respectively. The discrepancy in the dissolution factor among the three kerogen models matches well with the volumetric strain shown in Figure 1. The dissolution factor decreases as the moisture content rises in the three kerogen models, indicating more water molecules contribute to kerogen deformation. The dissolution factor for kerogens IIA and IID shows a gentle decrease, while that for kerogen IIC shows a sharp decrease after the moisture content of 1.5 wt %. This can be attributed to the dynamic distribution characteristic of water molecules in different mature kerogen models (Figure 4). For immature kerogen IIA (Figure 4a), there are abundant polar groups scattering near many small pores. The newly added water molecules are separately inserted into these small pores. For

Figure 4. Sketch of water molecule distributions in kerogen models. (a) Kerogen IIA. (b) Kerogen IIC. (c) Kerogen IID.

postmature kerogen IID (Figure 4c), there are only limited polar groups located on the surface of big pores. Water molecules aggregate into big clusters near the polar groups with the increase of moisture content. Both of the patterns can lead to a gentle change of dissolution factor. By contrast, mature kerogen IIC has more polar groups on the surface of middle pores (Figure 4b). Water molecules can fill these middle pores in pairs or groups. The strong repulsion between water molecules against the kerogen skeleton can aggravate the deformation and result in a rapid decrease of the dissolution factor. Similar to porosity, the dissolution factor shows a larger uncertainty for kerogen IIC (within 7.86%). The standard error of dissolution factor ranges between 2.34 and 4.58% for kerogen IIC. Note that water molecules loaded into the kerogen have two contrary effects on the pore volumes. On one hand, water molecules can occupy the pore volumes by sorption effect. On the other hand, water molecules can expand the pores by swelling effect. Two scenarios were adopted in this work to determine the contribution of sorption effect and swelling effect to porosity variation. For sorption effect, we performed the simulation of water loading on fixed kerogen structures to eliminate the effect of swelling effect. For swelling effect, we first obtained the total effect of moisture content on kerogen pores based on the flexible kerogen models, and then subtracted the sorption effect to obtain the swelling effect. Figure 3c shows the porosity differential variation (γ) as a function of moisture content. The moisture content is observed to have a neglected effect on γ, indicating a linear relationship between porosity and moisture content. The swelling effect results in positive γ, while the sorption effect and total effect lead to negative γ for the three kerogen models. The results validate the contrary effects of sorption and swelling upon water loading. Also, it indicates that sorption induced pore shrinkage dominates the effect of water uptake. Figure 5 presents the pore size distributions and pore compositions of kerogen models used in the work. In Figure 5a, we check the effect of kerogen maturity on the pore size distributions. The smaller pores (1−3 Å) decline and the bigger pores (5−6 Å) rise as the kerogen maturity increases. Note that we only present part of the pore size distributions of kerogen models in Figure 5a due to the limitation of kerogen size and computational method. There should be wider pore size distributions in these kerogen models, especially for kerogen IID. In Figure 5b, we show the compositions of two types of pores, ineffective pore (diameter