Methane Diffusion in a Flexible Kerogen Matrix | The Journal of

Jun 5, 2019 - This is confirmed by a bulk modulus of ∼1 GPa under dry geological .... So far, we have shown that accounting for matrix flexibility r...
0 downloads 0 Views 3MB Size
Article Cite This: J. Phys. Chem. B 2019, 123, 5635−5640

pubs.acs.org/JPCB

Methane Diffusion in a Flexible Kerogen Matrix Amaël Obliger,*,†,‡ Pierre-Louis Valdenaire,† Franz-Josef Ulm,†,§ Roland J.-M. Pellenq,†,§ and Jean-Marc Leyssale†,∥

Downloaded via UNIV AUTONOMA DE COAHUILA on August 7, 2019 at 09:03:22 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.



MultiScale Materials Science for Energy and Environment (MSE2), The Joint CNRS/MIT/Aix-Marseille University Laboratory, UMI CNRS 3466, Massachusetts Institute of Technology, Cambridge 02139, Massachusetts, United States ‡ Laboratoire des Fluides Complexes et leurs Réservoirs, E2S, UMR 5150, University of Pau and Pays de l’Adour/CNRS/TOTAL, Pau 64000, France § Department of Civil and Environmental Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge 02139, Massachusetts, United States ∥ Institut des Sciences Moléculaires, Université de Bordeaux, CNRS UMR 5255, 351 cours de la libération, Talence 33405, France S Supporting Information *

ABSTRACT: It has been recognized that the microporosity of shale organic matter, especially that of kerogen, strongly affects the hydrocarbon recovery process from unconventional reservoirs. So far, the numerical studies on hydrocarbon transport through the microporous phase of kerogen have neglected the effect of poromechanics, that is, the adsorption-induced deformations, by considering kerogen as a frozen, nondeformable, matrix. Here, we use molecular dynamics simulations to investigate methane diffusion in an immature (i.e., with high H/C ratio) kerogen matrix, while explicitly accounting for adsorption-induced swelling and internal matricial motions, covering both phonons and nonperiodic internal deformations. However, in the usual frozen matrix approximation, diffusivity decreases with increasing fluid loading, as evidenced by a loss of free volume, accounting for adsorption-induced swelling that gives rise to an increase in free volume and, hence, in diffusivity. The obtained trend is further rationalized using a Fujita−Kishimoto free volume theory initially developed in the context of diffusion in swelling polymers. We also quantify the enhancing effect of the matrix internal motions (i.e., at fixed volume) and show that it roughly gives an order of magnitude increase in diffusivity with respect to a frozen matrix, thanks to fluctuations in the pore connectivity. We eventually discuss the possible implications of this work to explain the productivity slowdown of hydrocarbon recovery from shale immature reservoirs.



INTRODUCTION As opposed to conventional reservoirs where they are found as bulk fluids, gas, or oil, hydrocarbons in shale reservoirs are strongly adsorbed within pockets of organic matter (OM) mostly made of kerogen but also, to a lesser extent, resins, asphaltenes, or heavy hydrocarbons.1 The poor connectivity of these pockets at large scales implies using hydrofracking2 to enhance gas/oil flow in shales. As a result, the hydrocarbon flow within shale is an extremely complex process, which is not fully understood.3 Of critical importance is that, despite fracking, the productivity of the hydrocarbon recovery process exhibits strong declines that vary significantly from a shale reservoir to another4−6 and considerably reduce the economic merits of shale hydrocarbons. The fracture network can be partly responsible for such anomalous declines, by introducing a pressure dependence on permeabilities7,8 or interfacial effects.9,10 Other large-scale phenomena, such as wetting, Knudsen diffusion, and so forth,9,11−28 may also be involved. The kerogen itself may also play an important role in the productivity decline. Although kerogen comprises both micropores (pore size ϕ < 2 nm) and mesopores (2 nm < ϕ > 50 nm), it has been recently shown that the mesopores, © 2019 American Chemical Society

isolated among the microporous phase, only plays a minor role on diffusion.29 Conversely, the percolating microporous kerogen is thought to be responsible for the measured ultralow permeabilities.30 However, its specific role regarding wells productivity remains unclear and unsufficiently documented so far. The recent development of realistic molecular models of microporous kerogen31−33 has allowed for several molecular dynamics (MD) studies on hydrocarbon transport properties within these models.34−37 It has been shown that transport is purely diffusive and mainly dictated by strong adsorption effects, resulting in a breakdown of the classical Darcy’s law in such confining situations.34 Common to these studies is the observation that the diffusion coefficient of the considered fluid phase decreases with the increasing fluid concentration.34−37 This result has been recently recasted within the framework of a free volume theory which allows predicting the Received: April 8, 2019 Revised: June 3, 2019 Published: June 5, 2019 5635

DOI: 10.1021/acs.jpcb.9b03266 J. Phys. Chem. B 2019, 123, 5635−5640

Article

The Journal of Physical Chemistry B

maturation of type I/II kerogen from H/C = 1.5 to H/C = 1.1 as investigated here. The choice of using methane instead of nalkanes as a guest molecule was motivated by the requirement of limiting the computation time, methane being lighter, diffuses faster. The kerogen model was prepared using a liquid quench MD simulation and carefully equilibrated for 3−5 ns using MD simulations in the NPT ensemble at various temperatures, pressures, and methane loadings as described in ref 41. For the sake of simplicity, it comprises only carbon and hydrogen; other elements, such as oxygen, sulfur, or nitrogen, in lower concentrations and supposedly playing little to no role on the poroelastic properties have not been considered. Having an H/C ratio of 1.1, the model corresponds to immature kerogen in the early steps of the catagenesis stage.1 Its texture is composed of short aliphatic chains (5 C atoms on average) and small ring clusters (10 C atoms per cluster on average) and has a balanced ring/chain ratio (0.7). From the mechanical point of view, the ring clusters that tend to stiffen the matrix are interconnected by the soft aliphatic chains, giving rise to an overall semiflexible behavior. This is confirmed by a bulk modulus of ∼1 GPa under dry geological conditions, slightly decreasing with increasing T, decreasing P and increasing w. We recall here, as discussed in ref 41, that we treat independently pressure and loading as the presence of hydrocarbons in kerogen results from the chemical maturation of the OM, not from any equilibrium adsorption process. As shown in Figure 1, the immature kerogen matrix exhibits significant swelling, up to 21% volumetric strain, at the maximum considered loading. The pore volume characterized in ref 41 consists in disordered micropores with pore size distributions extending only up to ∼7 Å for the empty matrix. At maximum loading, pores up to ∼13 Å are observed. To better quantify the evolution of porosity with loading, we have calculated the evolution of the accessible porosity φ, defined as the fraction of the volume in which a methane molecule can reside, and the accessible free volume ratio φf, defined as the fraction of the volume that is accessible to methane, yet not occupied by any methane molecule.48 As we can see in Figure 2, adsorption-induced swelling results in a significant increase of the accessible porosity, for instance from 7% in the empty matrix to 23% at 82.4 mg/g under typical geological conditions. More unexpectedly, as shown in Figure 2, increasing loading at constant T and P also leads to a significant increase in the free volume ratio φf, up to ∼16% at the highest loading and geological conditions. We

diffusion coefficients of linear alkanes from the only knowledge of the porosity and the fluid loading.36,37 A common approximation to all these studies is to consider kerogen as a rigid matrix, that is, an undeformable porous solid. Some authors even use dummy particles to mimic neglected compounds in real OMs to control kerogen microporosity in rigid models of kerogen.38−40 Whilst it can be a fair approximation in the limiting case of overmature kerogen,41 it is well known that most kerogens, as other porous carbons (coal, char, etc.), can significantly swell upon gas adsorption.42−44 Adsorption-induced deformation (swelling) has also been evidenced in three recent molecular simulation studies,41,45,46 showing important volume increases (up to 6−8%) under adsorption in typical conditions. An important consequence of these poromechanical couplings is that adsorption can lead to significant increases in porosity and hence possible increases in pore connectivity or free volume. A very recent numerical study accounts for flexibility effects on transport in a kerogen model where the control pressure is applied only in the direction of the flow while restraining some parts of the kerogen with “virtual nails”.47 In the present letter, we investigate the diffusion of methane within a flexible kerogen model and aim at disentangling the effects of (a) adsorption-induced swelling and (b) internal matricial motionscovering both phonons and nonperiodic deformations due, for instance, to fluid displacementon the diffusion coefficients.



METHODS The molecular model of kerogen used in this work, counting 8032 C and 8832 H atoms, is shown in Figure 1 as equilibrated under typical geological conditions (423 K, 25 MPa), without methane (left) and at a loading w = 82.4 mg of methane (542 molecules) per gram of kerogen, the highest investigated loading. We note that this loading corresponds to a fraction (∼10%) of the weight of n-alkanes produced during the

Figure 1. Snapshots of the kerogen matrix at 423 K and 25 MPa with no fluid adsorbed (left) and 82.4 mg of methane adsorbed per gram of matrix (right). Blue and white sticks represent covalent bonds between C and H atoms of the kerogen; orange spheres stand for methane molecules. Full simulation cell volumes are 188 and 227 nm3 for the empty (left) and loaded (right) systems, respectively.

Figure 2. Evolution of the accessible porosity φ (left) and the free volume ratio φf (right) with methane loading w at 400 K and 0.1 MPa (orange), 400 K and 25 MPa (blue), and 300 K and 25 MPa (red). Dashed and solid lines are linear fits. 5636

DOI: 10.1021/acs.jpcb.9b03266 J. Phys. Chem. B 2019, 123, 5635−5640

Article

The Journal of Physical Chemistry B

up to 20 ns. The obtained diffusion coefficient, 8.92 ± 0.49 × 10−9 m2/s, falls well within the error bars reported in Figure 3 on the basis of the 2.3 ns long simulations. As a consequence of the increase of porosity, free volume, and connectivity of the pore space associated to fluid loading, the self-diffusion coefficient of methane within a flexible kerogen matrix increases (see Figure 3) with increasing loading, which is the exact opposite of the well-accepted behavior obtained when performing the simulations in the rigid matrix approximation (i.e., where larger fluid concentrations imply lower diffusivities because of the lack of free volume to diffuse34−37). The increase of the diffusion coefficient with the fluid concentration can be captured by a free volume theory (eq 3) developed by Fujita and Kishimoto for swelling polymers,52−57 a prototypal example of semiflexible material. Within this model the selfdiffusion coefficient is expressed as ÄÅ É ÅÅ φ (w) − φ ÑÑÑ kBT ÅÅ f 0Ñ ÑÑ D(w) = expÅÅα ÅÅ φ φ (w) ÑÑÑ ξ0 (2) 0 f ÅÇ ÑÖ

notice that the free volume ratio increases linearly with loading according to φf (w) = φ0 + βw

(1)

where φ0 is the porosity of the empty (unswollen) matrix and β is a parameter that quantifies the increase of free volume upon adsorption. From linear fits, we have φ0 = 8.0, 6.1, 4.4 × 10−2 and β = 14.4, 12.9, 8.38 × 10−4 g/mg for conditions 400 K/0.1 MPa, 400 K/25 MPa, and 300 K/25 MPa, respectively. Self-diffusion coefficients of methane are computed from MD simulation in the NVT ensemble performed with an inhouse developed code. As in ref 41, the AIREBO potential49 is used to describe all the interactions (kerogen and all-atom methane molecules). AIREBO was chosen because of its accurate predictions of the elastic properties of graphite50 and, upon slight reparameterization of the Lennard-Jones part, the equation of states of alkanes.41 These features are required to obtain realistic mechanical properties and accurate estimations of the system density that are crucial to the investigation of solid−fluid poromechanical couplings. Temperature is fixed using a three-member Nosé−Hoover chain method as implemented using the reversible scheme of Martyna et al.51 A timestep of 0.25 fs is used, and periodic boundary conditions are applied in all directions. Simulations are run for systems equilibrated in ref 41 for various loadings under geological conditions (400 K, 25 MPa) and at slightly lower T (300 K) or P (0.1 MPa) to investigate the influence of these thermodynamic parameters. The constant volumes imposed during the NVT runs where diffusion is investigated correspond to the values averaged during the NPT equilibration runs of ref 41. Only a few hundreds of picoseconds are needed to reach the linear regime for the evolution of the average mean square displacement of the fluid molecules with respect to time. The diffusion coefficients are fitted in this linear regime (⟨r2⟩ ∝ Dt) up to 2 ns from 2.3 ns long runs.

where the prefactor kBT/ξ0 corresponds to the diffusion coefficient at infinite dilution with ξ0 the friction coefficient quantifying the average friction force experienced by an isolated methane molecule. The effect of the fluid concentration is taken into account by the exponential term that represents the average probability for a molecule to find a nearby cavity (free volume) into which it can diffuse. The coefficient α that controls the exponential increase accounts for the ability of the free volume to improve the connectivity of the pore space. Inserting eq 1 into eq 2 leads to ÄÅ ÉÑ ÅÅ ÑÑ kBT βw ÅÅ ÑÑ D(w) = expÅÅα Ñ ÅÅ φ (φ + βw) ÑÑÑ ξ0 (3) ÅÇ 0 0 ÑÖ Fitting this expression to the numerical results (Figures 3 and 4, solid lines) accounting for both the swelling and the internal matricial motions for different conditions, 400 K/0.1 MPa, 400 K/25 MPa, and 300 K/25 MPa gives α = 0.284, 0.270, and 0.134 and ξ0 = 7.82, 21.9 and 21.6 pN s m−1, respectively. We note that data obtained at the same temperature (400 K) show similar values of α, whereas those obtained at the same pressure (25 MPa) show similar ξ0, however, such a possible correlation requires more data to be validated.



RESULTS The diffusion coefficients for the three sets of (P,T) conditions and the various investigated loadings are shown in Figure 3. We note that some error bars at low loadings can be relatively important because of poorer statistics and higher confinement. To verify the convergence of the reported results, the simulations at 400 K, 25 MPa, and 12 mg/g were extended



DISCUSSION So far, we have shown that accounting for matrix flexibility results in a significant increase in diffusivity when increasing fluid loading. From a microscopic perspective, matrix flexibility effects can have two relatively independent origins: (i) a volume change (or swelling), and hence porosity increase, accompanying a change in fluid concentration and (ii) some internal motions at fixed volume, including phonons and nonperiodic deformations that lead to pore size/shape fluctuations and may assist or hinder diffusion. To disentangle the effects of swelling and matrix internal motions on the transport properties, we carried out a second series of simulations at 400 K and 25 MPa in which swelling was accounted for, as in Figure 3, but where the matrix atoms were frozen, as in the usual rigid approximation. The two sets of data are compared in Figure 4A. As can be seen, when only swelling is considered, the same qualitative behavior, namely, that D increases with increasing w

Figure 3. Self-diffusion coefficient of methane (semi-log plot) as a function of loading w for three (T,P) couples. Error bars (95% confidence intervals) are computed from five independent runs. Solid lines are the results of best-fit procedures using eq 3. 5637

DOI: 10.1021/acs.jpcb.9b03266 J. Phys. Chem. B 2019, 123, 5635−5640

Article

The Journal of Physical Chemistry B

Figure 4. (A) Evolution of the self-diffusion coefficient of methane with loading in the kerogen matrix at 400 K and 25 MPa, whilst accounting for both swelling and internal flexibility (s + f), as in Figure 3, or swelling only (s). (B) Evolution with loading of (left axis), the diffusion coefficients ratio Ds+f/Ds, highlighting the diffusion enhancement due to internal matricial motions, and (right axis), the mean pore diameter dm. Error bars in (A) (95% confidence intervals) are computed from five independent runs. All lines but the red dashed line in (B) are based on fits to eq 3. The red dashed line in (B) is a power law fit which serves as a guide to the eye for the evolution of dm with w.

is observed, which confirms that the evolution of the diffusion coefficient is mostly dictated by the change in free volume. However, it is clear from Figure 4A that neglecting internal flexibility leads to considerably underestimated values of the diffusion coefficient. Recalling that the free volumes are the same for the two sets of results, we thus show that internal motions have a strong effect on the fitting parameters α (0.55 to 0.27) and ξ0 (1360−21.9 pN s m−1). From Figure 4B, showing the ratios of the diffusion coefficients with full flexibility D(s+f) over the diffusion coefficients with swelling only D(s), one can clearly see that the diffusivity enhancement due to internal motions increases with decreasing fluid loading, from about 3 to 4 at the highest loading, down to a value of around 60 (i.e., close to 2 orders of magnitude) in the limit of zero loading according to the ratio of free volume fits (dashed blue line in Figure 4B). Also shown in Figure 4B is the evolution of the mean accessible pore diameter (dm), as a function of loading.58 As can be seen, dm increases with increasing w, yet remaining relatively low, especially at low concentrations, with the value for the empty matrix being 4.6 Å. For such small pores, compared to the diameter of methane (3.73 Å), one easily imagine than the deformations (fluctuations) of the porous volume play a key role regarding the mobility of the adsorbed fluid. Such an example of a fluid diffusion event assisted by an internal motion of the matrix is provided in Movie S1 (see the Supporting Information) and summarized in Figure 5 where one can see a methane transit between two initially disconnected pores via the opening of a window by a matricial dihedral rotation. Recently, a continuum theory has been proposed by Marbach et al.59 to describe the diffusion of fluid confined between fluctuating surfaces. By taking the two limiting cases of this theory, namely, the rapid surface fluctuations and the frozen surfaces, we obtain an expression for the ratio of the diffusion coefficients in the flexible and rigid cases as D(s+f)/ D(s) = (1 + 3(dσ/dm)2)/(1 − (dσ/dm)2), where dσ is the standard deviation of the accessible pore size distribution. This model predicts an increase of the ratio from approximately 1.1 at zero loading to 1.6 at 82.4 mg/g. Such a model does fail in capturing both the magnitude of the diffusion enhancement due to flexibility at low loading and its decrease in magnitude with increasing loading. However, we note that the model gives the correct order of magnitude at large loading where the

Figure 5. Snapshots taken from the movie of a MD simulation at 400 K, 25 MPa and low (12.2 mg/g) loading. (a) Methane molecule in the left hand side cannot diffuse toward right as the distance between the two purple atoms is smaller than the required space, roughly σ dmin = σCH4 + 2 2C = 7.09 Å . (b) Thanks to the yellow dihedral rotation, the distance has raised above dmin and the molecule can diffuse. Color code is cyan: kerogen C atoms; white: H atoms; pink: methane C atoms; purple: distance (in Å) between the purple C atoms; yellow: dihedral angle (in degrees) of the three yellow bonds (the full movie is available in the Supporting Information as Movie S1).

approximation made in ref 59, namely, that the fluctuations are small wrt the average pore size, is closer to be valid.



CONCLUSIONS In summary, by simulating methane self-diffusion in a microporous kerogen model while properly accounting for matrix flexibility effects, both adsorption-induced swelling and internal flexibility, we have unambiguously shown that diffusivity increases with increasing loading, as opposed to results obtained in the rigid matrix approximation. Qualitatively, this trend results from an increase in free volume with increasing loading, which can be described by a Fujita− Kishimoto free volume theory. Within this framework, a friction coefficient (ξ0) and a more phenomenological parameter (α) that accounts for the ability of adsorptioninduced swelling to improve the connectivity of the pore space, are needed to capture the evolution of the diffusion coefficient as function of the free volume ratio. We have also shown that these parameters also account for the instantaneous deformations of the matrix at constant volume (phonons and nonperiodic motions), which can tremendously assist methane mobility in micropores. 5638

DOI: 10.1021/acs.jpcb.9b03266 J. Phys. Chem. B 2019, 123, 5635−5640

Article

The Journal of Physical Chemistry B

(10) Apostolopoulou, M.; Day, R.; Hull, R.; Stamatakis, M.; Striolo, A. A kinetic Monte Carlo approach to study fluid transport in pore networks. J. Chem. Phys. 2017, 147, 134703. (11) Neimark, A. V.; Vishnyakov, A. Phase transitions and criticality in small systems: vapor-liquid transition in nanoscale spherical cavities. J. Phys. Chem. B 2006, 110, 9403−9412. (12) Firouzi, M.; Wilcox, J. Molecular modeling of carbon dioxide transport and storage in porous carbon-based materials. Microporous Mesoporous Mater. 2012, 158, 195−203. (13) Coasne, B.; Galarneau, A.; Pellenq, R. J. M.; Di Renzo, F. Adsorption, intrusion and freezing in porous silica: the view from the nanoscale. Chem. Soc. Rev. 2013, 42, 4141−4171. (14) Malek, K.; Coppens, M.-O. Knudsen self- and Fickian diffusion in rough nanoporous media. J. Chem. Phys. 2003, 119, 2801−2811. (15) June, R. L.; Bell, A. T.; Theodorou, D. N. Molecular dynamics studies of butane and hexane in silicalite. J. Phys. Chem. 1992, 96, 1051−1060. (16) Whitby, M.; Quirke, N. Fluid flow in carbon nanotubes and nanopipes. Nat. Nanotechnol. 2007, 2, 87−94. (17) Liu, Y.-C.; Shen, J.-W.; Gubbins, K. E.; Moore, J. D.; Wu, T.; Wang, Q. Diffusion dynamics of water controlled by topology of potential energy surface inside carbon nanotubes. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 77, 125438. (18) Striolo, A.; Gubbins, K. E.; Gruszkiewicz, M. S.; Cole, D. R.; Simonson, J. M.; Chialvo, A. A.; Cummings, P. T.; Burchell, T. D.; More, K. L. Effect of temperature on the adsorption of water in porous carbons. Langmuir 2005, 21, 9457−9467. (19) Chen, L.; Kang, Q.; Pawar, R.; He, Y.-L.; Tao, W.-Q. Pore-scale prediction of transport properties in reconstructed nanostructures of organic matter in shales. Fuel 2015, 158, 650−658. (20) Chen, L.; Kang, Q.; Dai, Z.; Viswanathan, H. S.; Tao, W. Permeability prediction of shale matrix reconstructed using the elementary building block model. Fuel 2015, 160, 346−356. (21) Chen, L.; Zhang, L.; Kang, Q.; Viswanathan, H. S.; Yao, J.; Tao, W. Nanoscale simulation of shale transport properties using the lattice Boltzmann method: permeability and diffusivity. Sci. Rep. 2015, 5, 8089. (22) Wang, J.; Chen, L.; Kang, Q.; Rahman, S. S. Apparent permeability prediction of organic shale with generalized lattice Boltzmann model considering surface diffusion effect. Fuel 2016, 181, 478−490. (23) Wang, J.; Kang, Q.; Chen, L.; Rahman, S. S. Pore-scale lattice Boltzmann simulation of micro-gaseous flow considering surface diffusion effect. Int. J. Coal Geol. 2017, 169, 62−73. (24) Mehmani, A.; Prodanović, M.; Javadpour, F. Multiscale, multiphysics network modeling of shale matrix gas flows. Transp. Porous Media 2013, 99, 377−390. (25) Clarkson, C. R.; Nobakht, M.; Kaviani, D.; Ertekin, T. Production analysis of tight-gas and shale-gas reservoirs using the dynamic-slippage concept. SPE J. 2012, 17, 230−242. (26) Firouzi, M.; Wilcox, J. Slippage and viscosity predictions in carbon micropores and their influence on CO2 and CH4 transport. J. Chem. Phys. 2013, 138, 064705. (27) Sinha, S.; Braun, E.; Determan, M.; Passey, Q.; Leonardi, S.; Boros, J.; Wood, A., III; Zirkle, T.; Kudva, R. Steady-state permeability measurements on intact shale samples at reservoir conditions-effect of stress, temperature, pressure, and type of gas. SPE Middle East Oil and Gas Show and Conference; Society of Petroleum Engineers: Manama, Bahrain, March 10−13 2013. (28) Ortiz-Young, D.; Chiu, H.-C.; Kim, S.; Voïtchovsky, K.; Riedo, E. The interplay between apparent viscosity and wettability in nanoconfined water. Nat. Commun. 2013, 4, 2482. (29) Berthonneau, J.; Obliger, A.; Valdenaire, P.-L.; Grauby, O.; Ferry, D.; Chaudanson, D.; Levitz, P.; Kim, J. J.; Ulm, F.-J.; Pellenq, R. J.-M. Mesoscale structure, mechanics, and transport properties of source rocks’ organic pore networks. Proc. Natl. Acad. Sci. U.S.A. 2018, 115, 12365−12370.

Although we could not compute collective diffusion coefficients in this work because of the added computational cost of dealing with a fully flexible matrix, we note, however, that noticeable fluid−fluid correlations would arise only at large concentration when the connectivity of the pore space is at its maximum. We can thus reasonably consider that the self and the collective diffusion coefficients are approximately equal, as that already holds in the rigid case.34−37 Finally, coming back to the perspective of hydrocarbon recovery from shale reservoirs, the present results suggest that the hydrocarbon diffusivity should considerably decrease during hydrocarbon expulsion (desorption) from the flexible microporous phase. This could possibly take an important part in the anomalous decline of hydrocarbon production from shale reservoirs, yet this deserves additional investigations.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.9b03266. MD simulation capturing the motion of a fluid molecule diffusing thanks to a dihedral rotation of carbon atoms in the matrix (MPG)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Amaël Obliger: 0000-0002-7765-1691 Roland J.-M. Pellenq: 0000-0001-5559-4190 Jean-Marc Leyssale: 0000-0002-7296-922X Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was supported by the CNRS and the A*MIDEX foundation. REFERENCES

(1) Vandenbroucke, M.; Largeau, C. Kerogen origin, evolution and structure. Org. Geochem. 2007, 38, 719−833. (2) Yethiraj, A.; Striolo, A. Fracking: what can physical chemistry offer? J. Phys. Chem. Lett. 2013, 4, 687−690. (3) Striolo, A.; Cole, D. R. Understanding shale gas: recent progress and remaining challenges. Energy Fuels 2017, 31, 10300−10310. (4) Patzek, T. W.; Male, F.; Marder, M. From the Cover: Cozzarelli Prize Winner: Gas production in the Barnett Shale obeys a simple scaling theory. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 19731−19736. (5) Cueto-Felgueroso, L.; Juanes, R. Forecasting long-term gas production from shale. Proc. Natl. Acad. Sci. U.S.A. 2013, 110, 19660− 19661. (6) Baihly, J. D.; Altman, R. M.; Malpani, R.; Luo, F. Shale gas production decline trend comparison over time and basins. SPE Annual Technical Conference and Exhibition; Society of Petroleum Engineers, 2010. (7) Walsh, J. B. Effect of pore pressure and confining pressure on fracture permeability. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1981, 18, 429−435. (8) Wheaton, R. Dependence of shale permeability on pressure. SPE Reservoir Eval. Eng. 2017, 20, 228−232. (9) Lee, T.; Bocquet, L.; Coasne, B. Activated desorption at heterogeneous interfaces and long-time kinetics of hydrocarbon recovery from nanoporous media. Nat. Commun. 2016, 7, 11890. 5639

DOI: 10.1021/acs.jpcb.9b03266 J. Phys. Chem. B 2019, 123, 5635−5640

Article

The Journal of Physical Chemistry B (30) Aguilera, R. Flow units: From conventional to tight-gas to shale-gas to tight-oil to shale-oil reservoirs. SPE Reservoir Eval. Eng. 2014, 17, 190−208. (31) Collell, J.; Ungerer, P.; Galliero, G.; Yiannourakou, M.; Montel, F.; Pujol, M. Molecular simulation of bulk organic matter in type II shales in the middle of the oil formation window. Energy Fuels 2014, 28, 7457−7466. (32) Bousige, C.; Ghimbeu, C. M.; Vix-Guterl, C.; Pomerantz, A. E.; Suleimenova, A.; Vaughan, G.; Garbarino, G.; Feygenson, M.; Wildgruber, C.; Ulm, F.-J.; et al. Realistic molecular model of kerogen’s nanostructure. Nat. Mater. 2016, 15, 576−582. (33) Atmani, L.; Bichara, C.; Pellenq, R. J.-M.; van Damme, H.; van Duin, A. C. T.; Raza, Z.; Truflandier, L. A.; Obliger, A.; Kralert, P. G.; Ulm, F. J.; et al. From cellulose to kerogen: molecular simulation of a geological process. Chem. Sci. 2017, 8, 8325−8335. (34) Falk, K.; Coasne, B.; Pellenq, R.; Ulm, F.-J.; Bocquet, L. Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nat. Commun. 2015, 6, 6949. (35) Collell, J.; Galliero, G.; Vermorel, R.; Ungerer, P.; Yiannourakou, M.; Montel, F.; Pujol, M. Transport of multicomponent hydrocarbon mixtures in shales organic matter by molecular simulations. J. Phys. Chem. C 2015, 119, 22587−22595. (36) Obliger, A.; Pellenq, R.; Ulm, F.-J.; Coasne, B. Free volume theory of hydrocarbon mixture transport in nanoporous materials. J. Phys. Chem. Lett. 2016, 7, 3712−3717. (37) Obliger, A.; Ulm, F.-J.; Pellenq, R. Impact of nanoporosity on hydrocarbon transport in shales’ organic matter. Nano Lett. 2018, 18, 832−837. (38) Collell, J.; Galliero, G.; Gouth, F.; Montel, F.; Pujol, M.; Ungerer, P.; Yiannourakou, M. Molecular simulation and modelisation of methane/ethane mixtures adsorption onto a microporous molecular model of kerogen under typical reservoir conditions. Microporous Mesoporous Mater. 2014, 197, 194−203. (39) Vasileiadis, M.; Peristeras, L. D.; Papavasileiou, K. D.; Economou, I. G. Modeling of bulk kerogen porosity: methods for control and characterization. Energy Fuels 2017, 31, 6004−6018. (40) Michalec, L.; Lísal, M. Molecular simulation of shale gas adsorption onto overmature type II model kerogen with control microporosity. Mol. Phys. 2017, 115, 1086−1103. (41) Obliger, A.; Valdenaire, P.-L.; Capit, N.; Ulm, F. J.; Pellenq, R. J.-M.; Leyssale, J.-M. Poroelasticity of methane-loaded mature and immature kerogen from molecular simulations. Langmuir 2018, 34, 13766−13780. (42) Bangham, D. H.; Fakhoury, N. The expansion of charcoal accompanying sorption of gases and vapours. Nature 1928, 122, 681− 682. (43) Yakovlev, V. Y.; Fomkin, A. A.; Tvardovski, A. V.; Sinitsyn, V. A.; Pulin, A. L. Adsorption-stimulated deformation of microporous carbon adsorbent. Russ. Chem. Bull. 2003, 52, 354−358. (44) Gor, G. Y.; Huber, P.; Bernstein, N. Adsorption-induced deformation of nanoporous materialsA review. Appl. Phys. Rev. 2017, 4, 011303. (45) Ho, T. A.; Wang, Y.; Criscenti, L. J. Chemo-mechanical coupling in kerogen gas adsorption/desorption. Phys. Chem. Chem. Phys. 2018, 20, 12390−12395. (46) Yu, S.; Bo, J.; Meijun, Q. Molecular Dynamic Simulation of Self- and Transport Diffusion for CO2/CH4/N2 in Low-rank Coal Vitrinite. Energy Fuels 2018, 32, 3085−3096. (47) Wu, T.; Firoozabadi, A. Effect of microstructural flexibility on methane flow in kerogen matrix by molecular dynamics simulations. J. Phys. Chem. C 2019, 123, 10874−10880. (48) The numerical procedure used to compute φ is described in details in ref 41. It consists first in introducing, using random attempts, a large number of spheres having the diameter of methane (σCH4) within the porosity (i.e. excluding overlaps with matrix C or H atoms), and, second, in computing the fraction of the total volume occupied by these spheres, using random drawings again. The accessible free volume ratio is defined as φf = (φV − nvCH4)/V where

V is the total volume, n is the number of methane molecules and vCH4 is the molecular volume of a methane molecule, considered as a sphere of diameter σCH4; the van der Waals diameters for methane and matrix atoms are σCH4 = 3.73 Å, σC = 3.36 Å and σH = 2.42 Å. (49) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions. J. Chem. Phys. 2000, 112, 6472−6486. (50) Farbos, B.; Da Costa, J.-P.; Vignoles, G. L.; Leyssale, J.-M. Nanoscale elasticity of highly anisotropic pyrocarbons. Carbon 2015, 94, 285−294. (51) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit reversible integrators for extended systems dynamics. Mol. Phys. 1996, 87, 1117−1157. (52) Fujita, H.; Kishimoto, A. Diffusion-controlled stress relaxation in polymers. II. Stress relaxation in swollen polymers. J. Polym. Sci. 1958, 28, 547−567. (53) Kishimoto, A.; Fujita, H. Diffusion-controlled stress relaxation in polymers. III. Stress relaxation in a swelling polymer. J. Polym. Sci. 1958, 28, 569−585. (54) Fujita, H. Notes on free volume theories. Polym. J. 1991, 23, 1499−1506. (55) Fujita, H. Free volume interpretation of the polymer effect on solvent dynamics. Macromolecules 1993, 26, 643−646. (56) Korsmeyer, R. W.; Lustig, S. R.; Peppas, N. A. Solute and penetrant diffusion in swellable polymers. I. Mathematical modeling. J. Polym. Sci., Part B: Polym. Phys. 1986, 24, 395−408. (57) Brazel, C. S.; Peppas, N. A. Modeling of drug release from swellable polymers. Eur. J. Pharm. Biopharm. 2000, 49, 47−58. (58) dm was calculated by averaging over the accessible part of the pore size distribution (i.e., the part with pore diameters larger than σCH4). (59) Marbach, S.; Dean, D. S.; Bocquet, L. Transport and dispersion across wiggling nanopores. Nat. Phys. 2018, 14, 1108.

5640

DOI: 10.1021/acs.jpcb.9b03266 J. Phys. Chem. B 2019, 123, 5635−5640