Reduced methane recovery at high pressure due to methane trapping in shale nanopores

By 2050, shale gas production is expected to exceed three-quarters of total US natural gas production. However, current unconventional hydrocarbon gas recovery rates are only around 20%. Maximizing production of this natural resource thus necessitates improved understanding of the fundamental mechanisms underlying hydrocarbon retention within the nanoporous shale matrix. In this study, we integrated molecular simulation with high-pressure small-angle neutron scattering (SANS), an experimental technique uniquely capable of characterizing methane behavior in situ within shale nanopores at elevated pressures. Samples were created using Marcellus shale, a gas-generative formation comprising the largest natural gas field in the United States. Our results demonstrate that, contrary to the conventional wisdom that elevated drawdown pressure increases methane recovery, a higher peak pressure led to the trapping of dense, liquid-like methane in sub-2 nm radius nanopores, which comprise more than 90% of the measured nanopore volume, due to irreversible deformation of the kerogen matrix. These findings have critical implications for pressure management strategies to maximize hydrocarbon recovery, as well as broad implications for fluid behavior under confinement. Liquid-like methane could be trapped in pores of less than 2 nm, rather than recovered, when higher peak pressure is applied during shale gas drilling. This was revealed by integrated molecular simulations and high-pressure small-angle neutron scattering on Marcellus Shale samples.

C urrently, shale reservoirs produce >60% of US natural gas, and are predicted to produce more than three-quarters of natural gas by 2050 1 . However, due to the small porosity and low permeability of tight shale formations, natural gas recovery rates using existing hydraulic fracturing methods are only around 20% 2 . Upon establishment of a tight gas reservoir well, shale gas production generally declines hyperbolically ( Fig. 1) 3 . The initial fast flushing of large, interconnected natural fractures around the well is believed to be responsible for the initial peak in production, while smaller fractures in the damage network, formed by hydraulic fracturing, lead to subsequent sustained production at a lower rate. The intersection of extensive microfracture networks, both natural and formed through hydraulic fracturing, with the nanoporous kerogen matrix is believed to play an important role in production during this stage. For example, Bažant et al. 4 found that the observed methane recovery of 15% was orders of magnitude higher than would be expected without a dense microfracture network. Production then further tapers off over extended time frames, as it becomes limited by the matrix diffusion or desorption of hydrocarbons from the shale matrix.
It is essential to understand how well management operations affect gas transport processes within the matrix in order to raise the production curve at later time periods, and thus maximize recovery of hydrocarbons stored in nanopores. In particular, the establishment of operational pressure controls to improve recovery from shale nanopores is a key strategy that can be applied without significant and costly changes to existing well infrastructure. Although it may take more time for operational pressure changes to propagate into nanopores, the existence of microfracture networks in low permeability shales will enhance propagation into the matrix 5 , and is an impetus for many existing studies on pressure-dependent methane transport in shale nanopores [6][7][8][9] .
Shale formations are sedimentary rocks that have long been considered as source rocks for hydrocarbons. Within these formations, shale nanopores are intrinsically heterogeneous and occur in both organic (e.g., kerogen) and inorganic (e.g., clay) shale components 10 . To better estimate the potential for gaseous hydrocarbon recovery, it is important to characterize these heterogeneous nanopores, as they store the majority of hydrocarbons and constitute the majority of pore volume. For example, mercury porosimetry analysis of Barnett shale found that 80% of pores have a pore size <5 nm 11 .
Owing to the high surface-to-volume ratio, the behavior of fluids confined in smaller nanopores (<100 nm) is dominated by interactions with the pore walls, and can hence differ drastically from behavior in bigger pores (micron-sized or larger), where bulk properties are more prevalent [12][13][14] . For instance, the behavior of water and ice in nanopores was studied extensively and nanoconfinement effects have been found to significantly alter the P-T phase diagram of water/ice 15 . In the case of shale, molecular dynamics simulations predict faster methane flow in nano-sized pores 7 . Nanoconfinement studies have important applications in various fields [16][17][18][19] , particularly nanoscience and nanotechnology [20][21][22] .
It is well documented through a combination of theoretical modeling and experimental observations that nanoscale confinement effects alter the properties of hydrocarbons in nanopores [23][24][25][26][27] . It was recently reported that in Marcellus shale, adsorbed methane forms multilayers in kerogen nanopores, contrary to the common assumption that gas adsorption in shale follows a monolayer Langmuir isotherm 28 . Capillary condensation has also been reported in shale nanopores 29 . According to this mechanism, pore capillary forces lower the vapor pressure of the adsorbate, causing it to condense into a liquid film on the pore walls. Multilayer adsorption and capillary condensation have been widely observed in nanopores [30][31][32] . While theoretical models of these physical processes exist 7,27 , direct experimental validation is challenging, particularly under high-pressure conditions prevailing in hydraulically fractured wells.
Small-angle neutron scattering (SANS) is a powerful technique that can be applied in combination with controlled sample environments to observe in situ the effects of field pressure/ temperature conditions on shale nanopore structure 33 . Neutron interaction with matter is quantified by a nuclear scattering length, the value of which varies considerably among light elements and their isotopes, especially, and relevant to this work, hydrogen and deuterium. Furthermore, neutrons have good penetrability through bulk samples in robust sample environments. These interactions, in addition to the absence of charge interactions, make SANS uniquely suited for the study of geomaterials. SANS has also recently been used to study confinement effects in silica nanopores 34,35 .
One particular benefit to SANS over other methods of measuring porosity is that SANS can measure both open and closed pores and can thus characterize how the matrix accessibility changes with increasing pressure 36 . However, while many previous SANS studies have measured porosity of different types of shale under ambient conditions [37][38][39][40] , only a few have examined pore size distribution and porosity changes in shale as a function of pressure 41,42 . Furthermore, no study has used SANS to examine how cycling pressure impacts the nanopore structure and methane retention.
Pressure management is one of few tools to control the production efficiency that is available to an operator and can be readily adjusted once a well is in operation. In this study, we use SANS to examine how pressure cycling affects methane behavior in Marcellus shale nanopores and, in particular, to link pressure maxima and pore size with methane recovery efficiency. We model the obtained SANS data to better understand the mechanism controlling methane transport and recovery. Our findings indicate that while high pressures are beneficial for methane recovery from larger pores, there is a trade-off involving the trapping of dense, liquid-like methane in smaller shale nanopores. These insights improve our understanding of confined fluid behavior during hydraulic fracturing, and can help us to optimize operational parameters that maximize hydrocarbon Fig. 1 Shale gas production curve. Theoretical shale gas production curve depicting the hypothesized physical mechanisms governing the production. The long tail of low production rates is due to limited transport within the shale matrix.
recovery, particularly at later operation times when recovery is limited by matrix processes.

Results
Small-angle neutron scattering (SANS). The goal of SANS measurements is to examine fluid behavior in shale nanopores in situ during pressure cycling. During this procedure, pores were filled with deuterated methane (CD 4 ), under pressure and the SANS profiles were measured at ambient pressure, then at discrete points with increasing pressure. After making a SANS measurement at a predetermined peak pressure, additional measurements were made at discrete decreased pressures, to remove methane from the pores, and finally again at ambient pressure. As discussed in the methods section, the distribution of methane within different sizes of pores is determined by the change in scattering intensity due to pressure-dependent differences in neutron contrast of pores filled with CD 4 relative to the original pore matrix. The resultant changes in the SANS profiles are used to compute the pressure-dependent apparent pore size probability distributions functions (PDFs) by inverting the intensity according to Eq. (1) in the methods section. The actual number and size distribution of the pores do not necessarily change with pressure. According to Eqs. (1) and (2) in the methods section, the pressure-induced size-dependent changes in the contrast are the main factor governing changes in the apparent PDFs.
The pore size distribution in the sample before pressurization with CD 4 was determined using the Irena "Size Distribution" macro with the maximum entropy method, assuming spherical pore morphology 43 . The SANS intensity profile, shown in Fig. 2a, plots the scattering intensity versus the scattered neutron wave number, Q. The power-law exponent, P, as defined in Eq. (3) in the methods section, is the magnitude of the negative exponent defining the slope of the scattering intensity at low-Q. For the initial sample, the power-law exponent is 3.58, while exposure to methane at the highest peak pressure of 6000 psi decreases the power-law exponent to a minimum value of 3.08. Low-Q powerlaw scattering with an exponent larger than 3, while usually interpreted as scattering from large length scale surface features 44 , may have contributions from volume polydispersity. Over the limited measured range in Q and without further information to discriminate these factors, we do not use the low-Q power-law portion of the data in computing the PDFs. Contributions from the low-Q power-law scattering and the background were instead subtracted from the scattering intensity profile before modeling the pore size PDF. Power law and background components are indicated in Fig. 2a. The PDFs thus cover a domain of sizes where we can reasonably assume that the scattering is due to the size distribution of smaller spherical pores (<10 nm), as expressed in Eq. (1).
A qualitative assessment of the changes in the power law and intensity in the low-Q domain gives information on the relative accessibility to methane over larger length scales (>10 nm). Regardless of the source of the power law, intensity changes are larger in the low-Q region, indicating large pores are more easily accessed by methane than smaller pores. The large change in power-law exponent is indicative of some discrimination in filled versus unfilled pores, or uptake in the local environment over larger length scales. If the power-law exponent change is due to scattering from surfaces, then there will be a transition from surface roughness contributions to contributions from a more convoluted pore filling structure 45 . Calculation of the expected intensity change due to pore filling, if it were solely a result of contrast changes, underestimates the intensity drop at low pressures and overestimates the intensity drop at high pressures.
This supports the interpretation that changes in the power-law exponent occur because methane filling of larger voids impacts both the contrast between the filled void and pore wall surface and the surface scatter from these pores, rather than just filling polydisperse pores. This calculation can be found in Supplementary Fig. 1 in the Supplementary Methods section of the Supplementary Information (SI). Fig. 2a presents the measured ambient SANS spectrum for Marcellus shale and the fitted spectrum, which is calculated using the pore size distribution shown in Fig. 2b. The fitted pore size distribution shows that a significant majority of pore volume in the size range probed by SANS (1 nm-100 nm) belongs to pores that are smaller than 2 nm in radius. SANS measurement of an additional Marcellus shale sample was consistent in both the power law (3.54) and fitted pore size distribution (Supplementary Fig. 2 in the SI). These results are consistent with previous SANS studies of the pore size distribution in Marcellus shale 46 .
To determine the effects of methane pressure cycling, the shale sample was put through two pressure cycles (Fig. 3). For cycle 1 Fig. 2 Marcellus shale pore size distribution. Measured and fitted SANS spectra for Marcellus shale pores under ambient conditions (a) before pressurization with CD 4 . For this graph, blue circles show the raw SANS profile, the black solid line is the fitted SANS profile, and background and power law are represented by the dashed blue and red lines, respectively. The fitted pore size distribution (b) indicates that the majority of these pores are <2 nm in radius. Error bars are included, but in some cases are smaller than the data marker.
(C1, Fig. 3a), SANS profiles were measured at ambient pressure with no methane, at 1500 psi (10.3 MPa), pressurized with methane, and at the peak pressure of 3000 psi (20.7 MPa). SANS was then measured at lower pressures during the drawdown cycle, first at 1500 psi, and again at ambient pressure. For cycle 2 (C2, Fig. 3b) the sample was re-pressurized with methane and measured at 3000 psi, 4500 psi (31.0 MPa), and then at a peak pressure of 6000 psi (41.4 MPa). The sample was then depressurized and measured at 4500, 3000, 1500 psi, and ambient. SANS profiles in Fig. 3 are limited to the peak and ambient pressures. Fig. 3 the SI.

Discussion
There are two significant differences between the SANS data obtained at the two peak pressures, 3000 and 6000 psi. The first key difference is in the scattering intensities of the ambient measurements before and after pressure cycling (Fig. 3c). The ambient measurement post C1 has a uniformly lower intensity over the entire measured Q range than the initial measurement at Small-angle neutron scattering log-log profiles for the first (a) and second (b) pressure cycles. For these graphs, red circles show initial ambient profiles, blue squares show peak pressure profiles, and purple triangles show the post-pressurization ambient profiles. At low Q, intensity decreases with increasing pressure due to pore filling by methane. Zoomed-in graphs on the right show unique methane behavior in small nanopores, which is reversible for the lower peak pressure cycle but irreversible for the higher peak pressure cycle. Comparison of ambient SANS (c) shows decreased methane retention at low Q for C2. For these graphs, blue squares represent the initial ambient profile, red circles represent the ambient measurement after the 3000 psi cycle, and purple triangles represent the ambient measurement after the 6000 psi pressure cycle. Error bars are included, but in some cases are smaller than the data marker. These error bars indicate measured error in scattering intensity at each Q. The confidence interval is the standard deviation. ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-020-00047-w ambient pressure before introduction of CD 4 , meaning that the contrast between pores and the surrounding matrix material is lower over all the length scales after the C1 pressure cycle or that the number of scatterers, i.e., pores, is decreasing. This decrease is most likely indicative of residual methane in the shale matrix rather than pore loss, since some of the scattering intensity is regained after the second pressure cycle. After the second, C2, cycle the ambient scattering was slightly higher than the C1 scattering for Q < 0.1 Å −1 , but still lower than the initial ambient scattering for that range. This indicates that while residual methane is also present after the 6000 psi peak pressure cycle, more methane was removed during the drawdown in C2 relative to the lower pressure, C1, cycle.
Molecular simulation modeling is used to shed light on the mechanism behind higher methane recovery from larger pores, as observed in the SANS measurements during drawdown in C2, the higher pressure cycle. We simulated methane recovery from a slit pore with varying widths where periodic boundary conditions were applied. Additional details can be found in the methods section. Fig. 4a shows the remaining numbers of molecules after drawdown from peak pressures of 3000 and 6000 psi at 300 K for different slit widths, the modeling equivalent to the pore channel diameter. These numbers are normalized to the density of methane at each of the peak pressures. While the channel length may impact retention extent, trends will be preserved. The model predicts less retention following drawdown from higher peak pressures. This trend matches experimental observation of SANS intensity changes after high and low peak pressure drawdowns, although this comparison is qualitative rather than quantitative (Fig. 4b). The greater decrease in intensity for the 3000 psi drawdown cycle indicates that there is more methane remaining in the matrix pore spaces. This finding is consistent with prior finding of pressure drawdown effects on methane production 47 . Nevertheless, this observation only holds for scattering at Q < 0.1 Å −1 , while unique trends are observed in the high Q regime (Q > 0.1 Å −1 ) between C1 and C2, as described below. Self-diffusion coefficients were also calculated for methane in these nanopores during the equilibrium state between the pore surface and methane (Fig. 4c). We observed a reduction in the diffusion coefficient as the pore width decreases or the pressure increases.
The second key difference is the high Q behavior, as highlighted in the zoomed-in graphs on the right-hand side of Fig. 3a and 3b. While the intensity at Q < 0.1 Å −1 decreases with methane pressurization relative to ambient scattering, the opposite trend is noted for the high Q region (Q > 0.1 Å −1 ). For the lower peak pressure cycle, C1, intensity increases relative to ambient scattering in the high Q region at 3000 psi, but decreases upon depressurization to a value that is only slightly below the ambient scattering intensity. For the higher peak pressure cycle, C2, a similar intensity increase is observed at high Q, but the intensity in that region remains elevated after depressurization. Thus, the mechanism responsible for this change appears to be reversible at a peak pressure of 3000 psi but irreversible after the higher peak pressure of 6000 psi. The difference in intensity caused by this increase is well outside the confidence interval of measurements, as indicated by the error bars in Fig. 3b. Thus, we are confident that this is a real observation.
Fitting the scatterer size distribution for the 6000 psi peak pressure measurement shows that this irreversible intensity increase is caused by increased scattering on a length scale <1 nm in radius, comparable with empty pore scattering from the ambient pore size distribution (Fig. 5). While this fitting gives an accurate picture of length scale where increased scattering occurs, it is not a measurement of increased quantity of scatterers (i.e., pores) because it does not account for differences in contrast caused by methane.
In a previous study, Ruppert et al. 41 observed similar increases in high Q scattering upon the introduction of pressurized methane to shale pore spaces. They hypothesize that this phenomenon could not be caused by the uniform condensation of CD 4 on pore walls, as this would require the density of the condensed methane to be~1 g/cm 3 , higher than that of liquid CD 4 (∼0.5 g/cm 3 ). Scattering intensity increases were thus attributed to the formation of clusters of CD 4 molecules in small nanopores, introducing pressure-dependent scattering length density (SLD) variations. The presence of these clusters in nanopores increases the number of scattering interfaces, and as a result the overall scattering intensity increases. However, as that study only collected data upon increasing methane pressure, it did not reveal the trapping of CD 4 clusters during methane pressure drawdown, reported here. Our observation has important implications for pressure management strategies to maximize recovery from the matrix at later operation times, particularly since the vast majority of the nanopore volume falls in the range where this trapping effect is observed (Fig. 2). This clustering behavior is supported by observations in the literature. Gidcumb 48 found that high effective pressures in nanopores leads to greater than expected densities of adsorbed methane, while Sun et al. 49 recently reported that methane will be inhomogeneously distributed within kerogen nanopores.
Close examination of previously reported results allows us to establish a mechanism for the observed phenomenon. Although it is well documented that increasing pressure leads to more methane uptake by kerogen in shale 50 , it is not well understood how this uptake process affects the kerogen pore structure. Recently, Tesson and Firoozabadi 51 used Molecular Dynamics-Grand Canonical Monte Carlo (MD-GCMC) simulations to shed light on this process. They found that up to 3000 psi, uptake of methane corresponds with kerogen swelling of up to~2.5%. When pressure is increased further to 6000 psi, the swelling percentage decreases slightly, despite a continued increase in methane uptake. This was attributed to mechanical compression, which leads to the deformation of the flexible kerogen matrix, changing the structure of kerogen pores, including their size, shape and connectivity.
For the first time, we present experimental evidence that this nanoscale deformation exists, and significantly impacts methane recovery from shale nanopores during depressurization. Combining our experimental observations with the previously modeled kerogen mechanical behavior at high pressures, we propose the following mechanism for methane trapping in smaller nanopores in kerogen during the high-pressure drawdown scenario. According to this mechanism, methane incorporation into kerogen results in swelling of the kerogen pore matrix. Up to 3000 psi, this swelling is reversible due to the high mechanical flexibility of the relatively "soft" kerogen, and we call this the "elastic" regime. However, continued methane infiltration beyond this pressure causes permanent kerogen matrix deformation, dubbed the "plastic" regime, closing pore throats and irreversibly trapping methane clusters in the kerogen pore space. This proposed mechanism is illustrated in Fig. 6.
Our observation of methane trapping due to kerogen deformation at high pressure has significant implications for the operation of wells in gas reservoirs. Pore size distribution indicates that the majority of nanopores in this Marcellus shale sample are <2 nm in radius, within the size range where the methane trapping behavior was observed. For shales within the gas window, as is the case for Marcellus shale, the predominant type of porosity is intraparticle kerogen porosity 52 . Hence, this mechanism can contribute significantly to methane retention, and operation changes such as lowering the peak pressure during cycling are recommended to prevent trapping in order to maximize recovery from the matrix at later times. Additional factors such as shale mineralogy, frac fluid chemistry, reservoir temperature, and overburden stress may also affect the trapping behavior. Further studies are needed to delineate how these factors enhance or diminish methane recovery from shale nanopores, with a goal of identifying the optimal parameters for maximizing hydrocarbon recovery at later time scales, when recovery is primarily from the shale matrix. This includes determination of the optimal pressure condition, which we anticipate to fall between 3000 and 6000 psi, that will prevent trapping while taking advantage of increased recoveries observed at high pressures. Additional pressure measurements will also allow us to quantify the relationship between pressure and kerogen elasticity in order to determine the pressure cutoff for plastic deformation. By relating this with kerogen properties such as thermal maturity, observations can be expanded to other shale formations.

Methods
Marcellus shale samples. Samples of Marcellus shale were provided by the Marcellus Shale Energy and Environmental Laboratory (MSEEL), a multidisciplinary and multi-institutional effort focused on improving unconventional resource recovery. MSEEL consists of two parallel legacy horizontal shale gas production wells in Morgantown, West Virginia. These wells are operated by North Northeast Energy and West Virginia University as a hydraulically fractured field site for the U.S. Department of Energy's National Energy Technology Laboratory (NETL). The Marcellus shale formation is the largest natural gas field in the U.S. in terms of its cubic feet of reserves. The initial reservoir pressure for Marcellus shale can vary but is generally reported to be between 3000 and 5000 psi [53][54][55] .
As part of MSEEL's initiative, the geology of this site has been wellcharacterized [56][57][58] . Marcellus shale at the MSEEL site primarily consists of siliciclastic mudstones with interbedded carbonates 57 . SANS samples were created from a Marcellus shale core, which was taken from a depth of 7544. 15-7546.80 feet. Powder from this core was analyzed using X-ray diffraction (XRD) and simultaneous thermogravimetric analysis and differential scanning calorimetry (TGA-DSC) (Fig. 7). XRD results (Fig. 7a) show that the sample consists primarily of quartz, illite, pyrite, calcite, and dolomite. This is consistent with previous observations of Marcellus clay being primarily illitic 57 . TGA-DSC results (Fig. 7b) show a total mass loss of 13.1%, with about 0.2% occurring before 200°C, indicative of water loss. At higher temperatures, there are three distinct endothermic heat events with peaks at 536°C, 735°C, and 876°C. The peak at 536°C indicates kerogen decomposition, confirming the presence of a significant amount of kerogen in the shale sample. The peak at 735°C is likely carbonate mineral decomposition, while the smaller peak at 876°C can be attributed to continued kerogen pyrolysis or the decomposition of inorganic minerals present at a lesser concentration than carbonates.
To prepare the MSEEL core for SANS analysis, rock samples were cut into 3/8" (9.525 mm) diameter, 300 µm thin sections. A picture of a sample can be found in Fig. 8a. This shale thickness was previously established to mitigate multiple neutron scattering that, if present, would complicate data analysis. The shale sections were mounted onto 1 mm thick quartz glass of the same diameter using epoxy on just the edges of the sample, so as not to interact with the ¼" (6.35 mm) diameter neutron beam. The sample size was chosen to be compatible with high-pressure titanium cells used for SANS gas experiments.
SANS methodology. Measurements were conducted on the NG7 30 meter SANS instrument at the NIST (National Institute of Standards and Technology) Center for Neutron Research (NCNR) 59 . Incident neutron λ was 6 Å with a resolution, Δλ⧸λ, of 13%. The sample-to-detector distances were set at 1, 4, and 13 m to give a Q (scattered neutron wave number) measurement domain of 0.003 to 0.43 Å −1 . The intensity distribution of the scattered neutrons as a function of their scattering angle, 2θ, was recorded by a 2D detector. These intensities were radially averaged for a given scattering angle using the NIST SANS/USANS data reduction package 60 to convert the data to differential cross sections per steradian per unit volume (cm −1 ) as a function of the scattered neutron wave number, Q = 4π/λsin θ(Å −1 ). Radially averaged SANS data are reported as absolute scale of differential crosssection per steradian per unit volume: Intensity, I (cm −1 ), as a function of scattering wave number, Q (Å −1 ). Q is inversely proportional to the pore radius, and intensity correlates positively with the pore quantity. Thus, the shape of the reduced scattering curve provides a powerful means of determining where methane is distributed among differently sized pores.
In SANS measurements, the neutron elastic scattering intensity is proportional to the average spatial fluctuations in contrast squared, Δρ 2 , over the volume sampled by the neutron beam. Here, Δρ ¼ ρ p À ρ m , where ρ p refers to the SLD of the pore and ρ m refers to the SLD of the surrounding matrix. The SLD of air in empty pores is approximately zero. When pressurized methane is introduced to shale, it enters accessible pores, changing the pore SLD and thus the scattered intensity. For SANS measurements, deuterated methane (CD 4 ) was used because of its positive bound coherent scattering length (3.33 × 10 −4 Å), as compared to that of hydrogenated methane, CH 4 (−8.31 × 10 −5 Å). By convention, the positive sign indicates scattering that is in phase with the incident beam wave function, whereas the negative sign indicates scattering that is 180°out of phase. Thus, by using CD 4 , the pore SLD is increased, and there is a decrease in the SANS intensity that is very sensitive pore filling. Furthermore, when the pressure of CD 4 is decreased during the drawdown phase of the pressure cycle, the extent to which the scattering intensity does not return to the original scattering is a sensitive measure of remnant CD 4 in the pores. Supplementary Table 1 in the SI contains the calculated matrix SLD and SLD of methane in the pores at each system pressure.
A general term for SANS intensity, I(Q), is given as: Here, N T is the number density of voids contributing to the scatter (cm −3 ); P(R) is the PDF of voids with size, R; and I(Q, R) is the scattering intensity from a single void depending on these parameters. The intensity from a single void is where V is the void volume and F is the form factor, which is the scatter from a single void, normalized for its volume and contrast. The contribution of the  contrast term for CD 4 filled voids will depend on the density and state of the deuterated methane and the composition of the surrounding matrix. However, the largest contribution to the scatter will be from the unfilled pores, which will have the largest contrast relative to the shale matrix due to the lower contrast between the matrix and deuterated methane filled pores. In Eq. (1), we ignore the inter-void connectivity (e.g., correlations in the arrangement of pores) and large-scale surface structure that contribute to the scatter at low-Q, as this contribution is factored out in the analysis. This low Q scattering from shale obeys a power law: For values of P between 3 and 4, scattering can be interpreted as resulting from surface geometry or from pore size polydispersity 61 . Values of P < 3 indicate particle and particle aggregate geometry, with non-integer values interpreted as being characteristic of fractal-like hierarchical structure 62 .
Pressure system. The gas pressure cell was machined out of titanium, with titanium windows. Titanium is relatively transparent to neutrons and has limited background scattering 63 . The shale sample was placed on a quartz glass spacer. The pressure system was sealed with 'o'-ring gaskets and the Ti windows were compressed against the spacer with a screw-on retaining ring. A SANS background for data reduction was taken using the assembled cell without a sample. The titanium cell has one inlet/outlet attached to the methane gas pressurization system. The pressure system consists of two ISCO 100HLf continuous flow pumps (Teledyne/ ISCO. Lincoln, NE, USA) (Note: The identification of any commercial product or trade name does not imply endorsement or recommendation by the National Institute of Standards and Technology). The two pumps were used in series to fill the pressure cell from a lecture bottle of CD 4 gas (Cambridge Isotopes Laboratory Inc, Andover MA, USA) and regulate pressure to set values from 1500 to 6000 psi. A schematic of the pressure system can be found in Fig. 8b. The pressurization process was controlled and automated using a program provided by NCNR that utilizes the LabView (National Instruments, Dallas, TX, USA) development platform. For SANS data collection at ambient pressures after methane pressurization, the cell was detached from the system and its pressure valve was opened.
Each SANS measurement took approximately two hours, and pressure was increased between measurements at a ramp rate of 300 psi/min with an equilibration time of 20 min. Previous high-pressure methane adsorption measurements on shale have reported an equilibrium time of <1 h, indicating that samples are close to equilibrium in our system 64 . As the ramp up and sample measurement times are identical for each pressure measurement, comparisons between the two cycles are valid even if the sample is not fully equilibrated with methane.
Modeling approach. Molecular simulation is a microscopic simulation approach, which considers the molecular structure of the matter, with the force field parameters describing the intermolecular and intramolecular interactions. The force field is represented by a set of equations optimized to maintain the molecular entity and reproduce its macroscopic properties. For atomistic simulation, the classical force field includes terms for bonded and non-bonded atoms 65 . They include the bond interactions between two atoms, the angle interactions between three atoms and dihedral interactions between four atoms.
There are two distinct techniques for system evolution in molecular simulations: Molecular Dynamics (MD) and Monte Carlo (MC) approaches. While the system evolution in MD simulations is a time dependent process controlled by Newton's equation of motion, MC simulation is optimized through the Markov chain to achieve more than a 50% acceptance rate for the system configuration states 66 . Therefore, Monte Carlo is more suited to study the equilibrium properties of the system and MD is more suited to study the non-equilibrium ones.
MD simulations have been used previously to quantify methane transport in shale nanopores as a function of slit width and pressure. Although the systems are not identical, these studies can provide some insight into shared mechanisms governing transport. For example, Yu et al. 9 found that with decreasing system pressure, interactions between the gas and wall increased, and this phenomenon is especially prominent in smaller nanopores. Xiong et al. 67 also modeled methane behavior in organic nanopores and found that with increasing pressure, the proportion of adsorbed methane gas decreased and most gas was instead in the free state.
To model methane production for comparison to SANS results, we started with an empty graphite pore and used Monte Carlo simulation to achieve the equilibrium state of the pore and methane. It is not possible to perform direct molecular simulation of the experimental data, which would require inclusion of heterogeneity of the porous media (both organic and inorganic materials) and sample size. Instead, we designed computational experiments to investigate salient aspects of the SANS measurements, reconciling the simulations with the modeling when possible. According to previous studies, graphite has a surface energy of 63 ± 7 mJ/m 2 68 , while the surface energy for shale ranges between 40-50 mJ/m 2 69 . Although these surface energies are somewhat different, graphite is still expected to capture confinement and pressure effects.
While shale nanopores are both rough and irregular, smooth graphite slit pore models have been widely used as a simplified corollary for kerogen nanopores in shale [70][71][72] . Surface roughness affects both the flow dynamics and the adsorption capacity of methane in pores [73][74][75][76] . However, given that our primary goal is to qualitatively compare retention in our model with neutron scattering observations, it warrants the use of smooth surface. Note that while the inclusion of surface roughness may change absolute values of methane retention, we expect that pressure and size trends hold as the roughness is present in all systems.
All the simulation were performed using LAMMPS 77 . We used TraPPE_UA force field to model the methane 78 . Fig. 9 presents model images and describes the workflow for MD simulations. We started with an empty graphite pore and used  Workflow of the molecular simulation starts with an empty silt pore. Monte Carlo simulation is used to achieve the desired equilibrium between pore and fluid. After that, the simulation system is extended to include a production region where the molecular dynamics describe the fluid transport from the pore to the production region.
Grand Canonical Monte Carlo simulation of the system for 1 million steps to populate the pore with methane at the specified physical conditions, achieving the equilibrium state between the pore surface and methane. We simulated the system for 5 ns, when the mean square displacement (MSD) is collected to estimate the self-diffusion coefficients using the Einstein relation 79 . Then, we conducted molecular dynamics of the system for another 1 ns using the canonical ensemble. After establishing equilibrium, we extended the simulation box in the Y-direction to define a production region, where methane can be extracted, which was controlled using the "Fix Evaporate" command. Modeling was done in three dimensions. The X and Y dimensions were set at 3.5 × 25.5 nm for all systems. The Z dimension was the slit width, a corollary to the pore diameter. Slits widths were set at 58, 12, 6, and 3 nm to attain the data represented in Fig. 4a.

Data availability
All SANS and source data shown in Figs. 2-5, and 7 can be found as Supplementary Data 1 in the Supplementary Information. Irena, which was used for all SANS data processing and analysis, can be downloaded for free here: https://usaxs.xray.aps.anl.gov/ software/irena.

Code availability
All the simulation were performed using LAMMPS with a TraPPE_UA force field model for methane. LAMMPS can be downloaded for free at: https://lammps.sandia.gov/. TraPPE_UA models can be downloaded for free at: http://trappe.oit.umn.edu/.