Simulation model

We used the open-source LAMMPS software to carry out molecular dynamics simulations, employing a velocity-Verlet algorithm with a timestep of 2 fs. Figure 2b illustrates the geometry of the simulations containing a porous membrane. We define the coordinate system such that the axis of the pores was aligned with the z axis, with the methane reservoir at the more negative end. The membrane pores are composed of CNTs with a zig-zag configuration, arranged into a triangular lattice with the rows of the lattice aligned in the x direction. Graphene sheets at the ends of the pores block the voids between the CNTs, with the centre of each nanotube aligned with the centre of an graphene ring. The nanotubes ends are cleared by removing atoms in the graphene sheet located within r+1.42 Å of the CNT centre, the extra 1.42 Å equal to the carbon–carbon bond length, ensuring a reasonable spacing between carbon atoms. We apply periodic boundary conditions in the x and y directions, and ‘shrink-wrapped’ boundary conditions in the non-periodic z direction, allowing the system to expand and contract as necessary. The initial state is created by running the simulation at 25 MPa for 1 ns while preventing the movement of methane out or water into the pores using an additional repulsive force field at the pore opening.

We use a Lennard–Jones potential to model the non-electrostatic forces between all particles, with the form

where r ij is the separation between two particles and r c =13.5 Å is the force cutoff distance for interactions not involving piston atoms (which we describe later). The values of the σ ij and ɛ ij parameters for different types of interacting particles are summarized for particles of the same type in Supplementary Table 2. For interactions between unlike particles, we use Lorentz–Berthelot combining rules, such that for particles of type α and β, and . The exception to these combining rules is the methane–piston interaction, which we describe below.

We use a united atom description of methane with the TraPPE force field30. Each methane molecule is represented by a single Lennard–Jones particle, with interaction parameters as listed in Supplementary Table 2. This model accurately describes the liquid–vapour coexistence curve and critical temperature of methane. We include at least 40 methane molecules per square nanometre of membrane cross-sectional area.

We employ the simple point charge (SPC) water model31, with the charges and Lennard–Jones parameters given in Supplementary Table 2. The bond angle is fixed at 109.47° and the hydrogen–oxygen bond length at 1 Å using the SHAKE algorithm32. For the electrostatic interactions, we truncate the force at 9 Å. Owing to the lack of periodicity in the z dimension, Ewald sum methods cannot be used for the long-range electrostatic forces. While the use of a cutoff affects the value of the surface tensions, this simplified water molecule still served the purpose of providing a polar fluid, which is immiscible with methane. We included at least 100 water molecules per square nanometre of membrane lateral area.

For simulations involving CO 2 , we use the ‘elementary physical mode’ described by Harris and Yung33, with interatomic interaction parameters as described in Supplementary Table 2. The carbon–oxygen bond length is fixed at 1.149 Å using the SHAKE algorithm, while constraining the bond angle θ by a harmonic potential with k θ =1,236 kJ mol−1 rad−2.

A Langevin thermostat acting on the methane and oxygen atoms with a friction coefficient of 0.01 fs−1 holds the system temperature at 423 K. The thermostat acts only in the x–y plane so as not to interfere with the transport of the fluid along the axis of the pores.

Opposing pistons with a graphene-like structure apply a prescribed pressure to each side of the membrane. Piston atoms are constrained such that they only move along the z axis. At each time step, the force on each piston atom is set to the average force of all atoms in the piston—causing them to move in unison—plus an additional component corresponding to the external pressure on that piston. We define the interaction between piston and fluid atoms by a Lennard–Jones potential truncated at the minimum energy (r c =21/6σ ij ), and therefore purely repulsive. To prevent methane accumulation at the downstream piston (in contact with the water phase), we use a large σ ij value of 6.8 Å for methane–piston interactions.

Activation time determination

In the simulations used to generate the data shown in Fig. 2c we used a membrane with a pore radius r=0.59 nm and pore spacing D=1.70 nm, and simulation dimensions in the periodic dimensions x and y of 1.704 and 2.951 nm, respectively, such that the system contained two nanotubes. After the initialization process described above, we equilibrated the system with no pressure difference for 2 ns, then linearly decreased the pressure on the left (positive z) side of the membrane containing the water over 10 ps. The simulation was run until escape was observed (and for a short time after). We report results for six pressure differences, for each of which we simulated seven trials with different but equivalent initial conditions. The activation times for the individual trials are shown in Supplementary Fig. 5.

Umbrella-sampling molecular dynamics simulations

We used umbrella sampling to measure the free energy as a function of the amount of extracted methane and gain insight into the origin of the energy barrier. Umbrella sampling is a method commonly used to study the thermodynamics of rare events inhibited by energy barriers34. An order parameter must first be identified describing the process of interest. We used the number of extracted methane molecules per unit area as an order parameter, n ex , which we have approximated using a sigmoidal function to make the variable continuous and differentiable:

where z i is the z coordinate of the ith methane particle, z p =1.8 Å relative to the centres of the carbon atoms in the external graphene surface and λ=0.25 Å. We associate a biasing potential ω(n ex ) with the order parameter to force the system into an ordinarily unlikely state. By measuring the probability distribution of n ex in the biased simulation, PB(n ex ), we deduce the unbiased free energy from the relationship

In practice, rather than attempt to find a bias potential, which allows the entire domain of n ex to be sampled, it is more practical to run many simulation ‘windows’. In each window we applied a simple harmonic bias potential so as to sample a particular range of n ex , with the form

where K and n i determine the strength and centre of the bias for the ith window. To implement such a bias in a molecular dynamics simulation, the force on each methane particle resulting from the bias must be described as a function of the particle position along the z axis:

The need for n ex to be differentiable with respect to the particle position z motivates the use of the logistic form in equation (9), for which the derivative is well known and straightforward to implement:

To calculate the full unbiased probability distribution PU, we used a weighted average of the unbiased probabilities in each window, , according to the weighted histogram analysis method34. The weightings were calculated so as to minimize the statistical error in PU:

where and N i are the measured unbiased probability and number of samples, respectively, in the ith window. The F i terms were calculated according to

We found a self-consistent solution by starting from equation (16) with F i =0 for all windows, then iterating between equation (17) and (16) until convergence34. We considered the solution to have converged when the maximum value of (1−F old /F i )2 for any window was <10−15, where F old is the value of F i in the previous iteration.

After unblocking the pores and applying the bias, we equilibrated the initial simulation window at n i =0 for at least 1 ns. To measure PB(n ex ), we sampled the system every 0.2 ps for at least 1 ns. We generated additional windows by restarting a simulation equilibrated with a similar n ex , equilibrating for at least 0.2 ns. We typically simulated windows separated in n i by 0.4 molecule per nm2, though when close to a free-energy maximum we sometimes required additional windows separated by 0.1 molecule per nm2, and with larger spring constant. We list the membrane parameters for the simulations used to produce the results in Fig. 3a in Supplementary Table 3 (system a).

Disordered membrane

To test the generality of the effects observed for the model CNT membrane, we carried out umbrella-sampling simulations in a system with a 5 × 5 × 5 nm disordered porous carbon membrane obtained using an atom-scale reconstruction technique. The pore size, chemical composition (including sp2/sp3 hybridization), and morphological disorder of this structure are comparable with those of kerogen35,36,37,38.

Heterogeneous membrane

In complex porous materials such as gas shale, a wide variety of geometries and surface chemistries may be present. We have explored the impact of such heterogeneities using a model composite membrane consisting of nanoporous hydrophobic component of width 3.12 nm, similar to the ordered hydrophobic membrane already described, and a hydrophilic component represented by a quartz surface of width 2.16 nm. This membrane is pictured in Fig. 2a(III) and the upper-left inset in Supplementary Fig. 2. The length of the simulation domain in the direction parallel to the stripes was 5.106 nm. The quartz surface was prepared by cutting a bulk quartz lattice39 along the (100) plane, then attaching protons to the resulting dangling oxygen bonds and allowing them to relax over a short NVE simulation. The quartz atoms are frozen during the simulations. The charge and Lennard–Jones parameters for the silicon, oxygen and hydrogen from which the quartz is composed are listed in Supplementary Table 4. The water was completely wetting on the quartz, while showing a contact angle of 133° on the graphene component (see Supplementary Fig. 6 and below in the Methods section), higher than typically observed for graphene because only a single layer was used. The free energy as a function of extracted methane for ΔP=0 and two additional pressure differences as measured using umbrella sampling is shown in Supplementary Fig. 2.

Free-energy barrier dependence on porosity

To test the linear scaling with 1−φ, we used umbrella sampling to determine the free energy as a function of n ex for different combinations of pore diameter and spacing listed in Supplementary Table 3 (systems a, d–h). To calculate the porosity, a correction term c was added to the pore radius to correct for the finite size of the carbon atoms, such that . The results confirm the linear relationship in equation (2), as seen in Fig. 3, with the best fit found for S=16.55±0.41 mJ m−2 and c=0.10±0.03 nm.

Surface tension calculation

To verify the proportionality with the spreading parameter in equation (2), we measured the surface tensions of the three interfaces of interest in a separate set of molecular dynamics simulations. We prepared three systems: one in which a water phase and methane phase are held in contact under pressure from two graphene pistons, much like the simulations already described, but without any membrane separating the phases; another with only a water phase between the pistons; and another with only methane. The systems are 5-nm wide in the x and y directions, and contain 7,200 water molecules and/or 910 methane molecules. The potentials used to model interactions between the piston atoms and other species are identical to those used for the membrane atoms in umbrella sampling simulations, listed in Supplementary Table 2, with the force truncated at 13.5 Å. We equilibrated the systems for 2 ns before beginning the surface tension measurement. In the methane–water mixed system we initially equilibrated for an additional 1 ns with a planar force separating the phases. We simulated the methane–water systems for 15 ns, water-only for 18.4 ns and methane-only for 10 ns, taking samples of the configuration every 1 ps.

We can express the surface tension between phases A and B using the Kirkwood and Buff approach40, as an integral across the interface of the difference between the pressure components normal and tangential to that interface, p N (z) and p T (z):

where the z axis is normal to the plane containing the interfaces, and z A and z B correspond to points within the bulk of the A and B phases. The normal and tangential pressures can be written in terms of the pressure tensor elements, p N =p zz and .

The pressure tensor components themselves contain a kinetic component, equivalent to that of an ideal gas, plus a second component taking into account the interactions between molecules. The diagonal elements of the pressure tensor at a position z, p αα (z), where α=x, y or z, can be calculated using the expression

where the sums are over each interatomic interaction between each pair of molecules, each molecule being composed of n i atoms. α i and α ia denote the α coordinate of the centre of mass of the ith molecule, and the ath atom of the ith molecule, respectively. r iajb is the separation between atoms ia and jb, ρ(z) is the fluid density and A is the simulation cross-sectional area in the x–y plane. We count the frozen graphene atoms as individual molecules when calculating p zz , but we neglect them in the calculation of p xx and p yy due to the symmetry in these dimensions.

The function depends on the choice of the contour drawn between particles and along which the contribution of that interaction is distributed. There is no unique definition for this contour41. We use the Irving and Kirkwood convention42, which has been shown to agree with other methods43. The contour in this model is a straight line connecting the two interacting molecules, such that ξ(z, z i , z j ) takes the form:

where H(x) is the Heavyside step function. In practice, we divide the system along the z axis into slabs of width 1 Å, and split the contribution to the pressure tensor of each interaction equally between any slabs between or containing the two interacting particles.

At equilibrium, the normal component of the pressure tensor must in principle be constant and equal to the pressure applied by the pistons. In practice the length of time for which the system must be sampled before this limit is reached for the water phase is very long. We have made the assumption that the normal pressure would eventually converge in our calculation of the surface tension.

The tangential and normal pressure profiles calculated using this method are shown in Supplementary Fig. 7. The measurements of γ MA and γ MW in the mixed system are consistent with the more precise measurements in the single-component systems.

Surface energy minimization with Surface Evolver

To overcome the limitation imposed by the small molecular dynamics simulation box size and construct a model for the critical nucleus, we have employed a mesoscale thermodynamic approach. We use the energy minimization algorithm implemented in the open-source Surface Evolver programme22 to determine the nucleus geometry at a range of volume.

Surface Evolver represents surfaces by mesh of triangular facets. In a standard minimization step (default ‘g’ command), the force on each vertex is calculated based on the gradient of the free energy, and then the vertices move according to the strength and direction of that force. Alternatively, the ‘Hessian seek’ command allows the Hessian matrix of second derivatives to be used to find the minimum energy configuration in the direction of motion as determined by the forces. These two techniques can be alternated while minimizing the energy. Further details of the minimization algorithm can be found in the Surface Evolver manual.

In the absence of a pressure difference, the total energy of the system can be described by

where γ base (x,y) is a periodic function of circular domains representing the pores, arranged in a triangular lattice. The function is equal to −γ AW within these circular domains, and γ MA −γ MW everywhere else.

The surface tensions calculated using the molecular dynamics simulations described above were used for the minimization calculations.

During the evolution of a surface in Surface Evolver it is usual to start with a crude approximation to the final geometry, then alternating between moving towards the minimum energy, and refining the surface by splitting existing facet in half. We restricted the refinement of the base of the nucleus such that new vertices are only created on the contact line, since additional vertices within the base perimeter are redundant, experiencing no net force.

The initial geometry consists of 50 vertices arranged on the substrate in a circle to form the contact line. Each of these contact line vertices are connected by an edge to their two neighbours, to a vertex in the centre of the circle, and to a vertex positioned above the centre of the circle at the height required for the target nucleus volume. Vertices located on the contact line are constrained such that they only move in the plane of the surface. Four initial contact radii were tested for each volume, and all but the lowest energy result discarded.

Surface Evolver minimizations of drop or bubble geometries on patterned surface can be challenging due to a tendency for the system to become stuck in local minima. If care is not taken a situation can arise in which even the liquid–vapour interface obviously fails to adopt a physically reasonable shape. We have found an effective strategy to avoid such problems is to alternately evolve the surface with the contact line vertices fixed in place or allowed to move freely. We used a combination of the regular linear gradient decent method (the default ‘g’ command) and the Hessian seek method (‘Hessian_seek’ command), described above. Motion of the vertices in both cases is multiplied by a ‘scale factor’ (<1), automatically calculated by the software, which dampens the motion as an energy minimum is approached. If the scale factor approaches zero, the surface stops evolving, while not necessarily having found the minimum energy configuration.

The algorithm described in Supplementary Note 1 was used within Surface Evolver to minimize the interfacial free energy.

Cassie–Baxter model of the nucleus

Here we describe the derivation of the thermodynamic model to predict the size and energy of the critical nucleus, which we compare with the Surface Evolver results in Fig. 4.

Consider a methane nucleus with volume V act on a porous surface. We make the simplifying assumption that the nucleus adopts an idealized spherical cap geometry with base radius R and an effective contact angle θ eff . The base radius R, contact angle, volume and curved surface area A AW (methane–water interface) are related by

where and . In the limit of low θ eff , β v (θ eff )=θ eff /4 and β a =1.

Assuming the length scale of the pattern is significantly smaller than the radius of the nucleus, the total free energy can be written

In the case of the simple geometry used in the molecular dynamics and surface evolver calculations, . The free energy can be rewritten in terms of the contact angle of the methane on the solid phase:

into which we may substitute the expressions for the area, volume and the Cassie–Baxter effective contact angle, to produce an expression for the free energy in terms of γ AW , V act , θ eff and ΔP:

where . In the limit of low θ eff , . It is straightforward then to determine the maximum in the free energy:

where we introduce the Kelvin radius R K =γ AW /|ΔP| and a geometric factor . In the small angle limit .

Measuring hydrophobicity of the membrane

To confirm that the model membrane is hydrophobic, we carried out a molecular dynamics simulation of a water drop on a graphene substrate. We simulated 700 water molecules represented by the SPC model described above section in an initially cubic lattice. During the initial 0.5 ns, a velocity rescaling thermostat ramped the temperature linearly from 400 to 300 K—a process we find allows the drop to relax more rapidly. The drop momentum in the xy plane was reset to zero every 20 ps during this first 0.5 ns. After the initial relaxation, the temperature was held at 300 K using Langevin thermostat. After equilibrating for a further 3.5 ns the system was sampled every 1 fs over a 2-ns period so as to measure the contact angle.

The radial density profile of the drop, measured relative to the centre of mass, is shown in Supplementary Fig. 6. To obtain the position of the surface as a function of z and , the system is divided into slabs parallel to the substrate and 0.5 Å thick in z. In each slab we fit the radially averaged density to the function via , the position of the interface; and , the bulk density, fixing nm. To measure the contact angle, a circle is fitted to r s (z) in the central region of the drop, avoiding deviations due to contact line tension near the base and low density regions near the top.

Statistical model of long-time recovery

Here we describe with additional detail the derivation of the statistical model for long-time recovery kinetics. We begin by considering the shale as containing a large number of trapped gas pockets, A, each containing a volume of gas V 0 , and each of which must overcome some energy barrier ΔG* before the gas within may be recovered. This scenario is illustrated schematically in Fig. 1. These pockets have a wide distribution of energy barriers as a result of variations in the local geometry and surface chemistry of the shale. The variation in energy barriers leads in turn to variations in activation time τ act , which could potentially span many orders of magnitude.

The probability of a pocket with a particular energy barrier being overcome at a particular time is given by . Once the energy barrier associated with a pocket has been overcome, let the volume recovered from that pocket vary with time according to =V 0 Ф(t), which increases from zero to a maximum of, V 0 over a time τ h .

We can write the total recovered volume at a given time in terms of an ensemble average over all the gas pockets of the probable amount of gas extracted from a given pocket at that time, given by equation (5).

Given a distribution of activation times, p τ , we can write equation (5) as an integral:

To explore the scaling of gas recovery with time according to equation (28), we assume an exponential distribution of energy barriers:

This simply represents one possible distribution, and the precise form is not critical to our analysis. The corresponding distribution of activation times is

where we define ; θ 0 and γ 0 fix the typical range spanned by these parameters over the reservoirs.

Taking this distribution of barriers and substituting into equation (28), we find the expression

where Γ(α,t) is the incomplete gamma function. The solution to this expression is obtained for the two regimes t<τ h and t>τ h by using the Laplace transform for the convolution to get the asymptotics, with the results and , respectively, as discussed above.

Data availability

The data that support the findings of this study are available from the corresponding author on request.