Evidence of Facilitated Transport in Crowded Nanopores

Fluid transport in nature often occurs through crowded nanopores, where a number of phenomena can affect it, because of fluid–fluid and fluid–solid interactions, as well as the presence of organic compounds filling the pores and their structural fluctuations. Employing molecular dynamics, we probe here the transport of fluid mixtures (CO2–CH4 and H2S–CH4) through silica nanopores filled with benzene. Both CO2 and H2S are strongly adsorbed within the organic-filled pore, partially displacing benzene. Unexpectedly, CO2/H2S adsorption facilitates CH4 transport. Analysis of the trajectories suggests that both CO2 and H2S act as vehicle-like carriers and might swell benzene, generating preferential transport pathways within the crowded pore. The results are useful for identifying unexpected transport mechanisms and for developing engineering approaches that could lead to storage of CO2 in caprocks.


Methods and Algorithms
Simulation Setup. Organic-rich shale formations are source rocks as well as the reservoir basement and caprocks that trap oil and gas. To build nano-pores representative of those found in model organic-rich shales, we considered an 13.4´13.4´15 Å 3 hydroxylated amorphous silica substrate with 4.5 OH nm -2 on its surface, taken from those published by Ugliengo et al., 1 as a unit cell. We duplicated this unit-cell to create two parallel silica slabs placed at a distance such to yield a slit-shaped nano-pore of width ~20 Å.
The X and Y dimensions of this substrate were Lx,p = 53.55 Å and Ly,p = 53.15 Å, respectively. We saturated the pore edge surfaces with OH groups to achieve a surface density of 4.5 OH nm -2 , which is consistent with experiments. 2 The simulation box is periodic in the three directions. Its Y dimension (53.15 Å) reflects the periodicity of the silica substrate; the X and Z dimensions were set to 193.55 and 31.71 Å, respectively.
Due to periodic boundary conditions, the model pore is effectively infinite along the Y direction.
Conversely, the pore is finite along the X direction, along which it is exposed to the fluid reservoirs, as illustrated in Figure S1, panel A. The pore obtained following this procedure could represent nano-pores existing in some shale gas plays with tectosilicates (> 30 wt% of quartz) as main components. [3][4] To model organic content trapped within the amorphous silica pore, we filled the slit-shaped silica pore with 400 benzene molecules, resembling aromatic oil trapped in sedimentary rocks. This model is meant to represent organic-rich shale caprocks that contain significant amount of organic carbon (> 11.7 wt %). 4 The number of benzene molecules introduced were sufficient to fill the pore volume and form thin layers on the solid substrate outside the pore (see Figure S1, panel A).
Force Fields. The CLAYFF force field was implemented to simulate the SiO2 substrates. 5 Silicon and oxygen atoms were tethered to their initial positions by applying a harmonic restraint force with a spring constant of 100 kcal/mol.Å. The surface hydroxyl hydrogen atoms were allowed to vibrate. The transferable potentials for phase equilibria (TraPPE) force field was implemented to model CO2, H2S, and CH4 molecules. [6][7][8] The TraPPE models for CO2, H2S, and CH4 were found to reproduce experimental results of vapour-liquid phase equilibria for the pure components as well as CO2-CH4 and H2S-CH4 binary mixtures, achieving high accuracy. [6][7] Benzene molecules were modelled using the second generation of general AMBER force field (GAFF2), 9 which exhibits a moderate improvement with respect to GAFF1 in its ability to reproduce intermolecular energy, liquid density, heat of vaporization, and hydration free energy for organic molecules. 10 Non-bonded molecular interactions were modelled by means of dispersive and electrostatic forces. The electrostatic interactions were modelled by the Coulomb potential, with long-range corrections treated using the particle-particle particle-mesh method (PPPM). 11 Dispersive interactions were modelled by 12-6 Lennard-Jones (LJ) potentials. The LJ parameters for unlike interactions were determined by Lorentz-Berthelot combining rules from the values of like components. 12 The cut-off distance for all interactions was set to 12 Å. No corrections were applied for long-range LJ interactions because it was found that such corrections do not affect significantly the amount of fluid adsorbed in the pore used here.
Algorithms. We conducted both equilibrium and boundary-driven non equilibrium MD (BD-NEMD) 13 for studying CO2-CH4 and H2S-CH4 binary mixtures passing through the benzene-filled pore. Initially, we equilibrated the system containing benzene in the presence of (a) pure CH4, and (b) CO2-CH4 and H2S-CH4 binary mixtures. The equilibration simulations lasted for 50-100 ns, with no external force applied. The simulations conducted for pure CH4 served as a base case against which to quantify the effect of CO2 or H2S on CH4 transport. Several mixtures were considered, as listed in Figure S1, panel B. The numbers of CO2 and H2S molecules were varied from 50 to 800. The number of CH4 molecules in each binary system was also changed, to maintain the pressure outside the pore for all systems at ~13.9 MPa at equilibrium.
The amount of CH4 needed for each system was calculated using the Peng-Robinson equation of states (EOS), based on the molecular density of CO2, H2S, and CH4 outside of the pore, and the customary mixing rules for multi-component mixtures. [14][15][16] The BD-NEMD simulations were conducted by applying a force along the X-axis to all CO2, H2S, and CH4 molecules located in a thin slab of width dext = 20 Å outside of the pore (within the permeate) to establish and maintain a constant pressure difference across the pore. In each simulation, 100 ns of simulations were required to reach steady state. Equilibration was confirmed by examining the variations in the energy and temperature as a function of simulation time, as well as CO2, H2S, and CH4 density profiles along the pore. Once steady state was achieved, production simulations were conducted for 50-150 ns.
Implementation. Equilibrium and BD-NEMD simulations were conducted using the package LAMMPS. 17 Simulations were performed in the canonical ensemble (NVT), where the number of particles, volume, and temperature are kept constant. The equations of motion were solved by implementing the leapfrog algorithm 18 with 1.0 fs time steps. The simulated temperature was maintained at 300 K by Nosé-Hoover thermostats 19-20 with a relaxation time of 100 fs. We applied separate thermostats to fluid molecules (CO2, H2S, and CH4), solvent molecules (benzene), and to the atoms in the silica substrates. [21][22] This allowed us to reduce the perturbations on the dynamics of the system due to the application of the external force, which causes a constant energy input on fluid molecules. mixtures inside the benzene-filled pore at 300K using equilibrium MD simulations. Assuming that in the range of compositions considered the Henry's law is applicable, the solubility coefficient might be obtained as the slope of the adsorption isotherms shown in Figure S2 at the limit of zero mole fraction of CO2 and H2S (or the mole fraction of CH4 approaches 1). From the simulation results we estimate that the solubility coefficient for methane in the crowded pore is lower in CO2-CH4 (~3.8´10 -3 Å -3 ) and H2S-CH4 (~4.6´10 -3 Å -3 ) mixtures, compared to that of CO2 (11.6´10 -3 Å -3 ) and H2S (32.2´10 -3 Å -3 ). Note that solubility coefficients are calculated based on the volume available to the fluids in the pore.
To assess whether all simulated processes are reversible and that the results presented are representative of equilibrated systems (see Figure 1A and 1B in the main text), we carried out CO2/H2S adsorptiondesorption cycles. By employing the 'evaporate' and 'deposit' procedures available in the software package, 17 we extract CO2/H2S molecules and simultaneously insert CH4 molecules randomly into the benzene-filled pore as well as the bulk reservoirs starting from the configuration of system 8C/7H (see Figure S1, panel B for system composition), respectively, to slowly retrieve the compositions of the other binary mixtures. The results for solubility of CO2, H2S, and CH4 and adsorption and re-adsorption of benzene during CO2/H2S adsorption-desorption cycles are shown in Figure 1A and 1B in the main text. Benzene Displacement. In Figure S3, we report the number of benzene molecules displaced out of the SiO2 nano-pore upon changing the composition of the binary mixtures simulated. Focus is on the number of CO2 (blue) and H2S (yellow) molecules adsorbed within the pore. The results suggest that H2S is more effective than CO2 in displacing benzene. Our results suggest that ten H2S molecules displace four C6H6 molecules from the pore, while ten CO2 molecules displace only three C6H6 molecules.  Figure S4, the density profiles of CH4  In Figure S5 we show one representative simulation snapshot (top panel) and the averaged density profiles (bottom panel) for CH4 and CO2 molecules distributed along the length of the simulation box along the X direction when the applied external force is 0.01 kcal/(mol.Å). The results are shown for the system 3C (see Figure S1, panel B, for its composition). The results show that the external force yields a pressurised zone on the right of the porous media, within which an increase in CH4 and CO2 densities is observed (Figure S5, bottom panel). This is the 'retentate' volume. The external field is considered equivalent to imposing a density (or a pressure) gradient, leading to a macroscopic flux in the direction of the arrow in Figure S5, top panel. Once a density gradient is imposed, CH4 and CO2 molecules diffuse from the retentate to the 'eluate' volume, located on the left of the pore, where the fluid density is lower.

Fluid Density Profiles at Equilibrium and During Flow. In
In Figure S6, we report the in-plane density distributions of benzene carbon atoms within first, second, third and fourth layers at equilibrium (left) and during flow (right) for the system shown in Figure S5, top panel. The high-density areas (red-yellow spots) of the contour plots indicate positions where the benzene molecules preferentially reside. The contour plots for the first and fourth layers at equilibrium and during flow suggest that some benzene molecules in these layers are parallel to the surfaces, as illustrated by the presence of planar hexagonal rings reminiscent of the benzene carbon atoms, while others orient themselves at an angle with respect to the surfaces. The results in Figure S6 show that the benzene molecules in the middle of the pore manifest a somewhat less ordered structure than those close to the pore surfaces. Our analysis reveals that the positions within the pore where benzene molecules occupy once the external forces are applied are similar to those in which they accumulate at equilibrium, as long as the applied pressures are moderate. The same qualitative results are also observed for the H2S-CH4 binary mixtures. This suggests that the external force of 0.01 kcal/(mol.Å) is not large enough to significantly perturb the structural properties of confined benzene in the presence of CH4 with either CO2 or H2S. Note that the transport diffusion coefficient Dt is quantified in the limit of the eternal force to zero, in which case the structure of the benzene-filled pore remains unchanged; hence it is essential to apply such a small external force, which would guarantee that the morphology of pore saturated with benzene is uninterrupted. Upon applying the external force, the density profiles show (see Figure S5, bottom panel) that once the flow is induced, the molecular density profiles show a decreasing gradient from the retentate (right) to the eluate regions (left).
For the density profiles, the centres of mass of the various molecules are considered.
In Eq. (1), Ni + and Nidenote the number of molecule species i that have passed the plane from right to left and left to right, respectively, in the representation of Figure S5.
Once the flux is known, permeability, K, can be determined based on Darcy's law, 23 which shows a linear relationship between molar flux and pressure drop across a porous medium.
The pressure drop, dP, for each species i (CH4, CO2, or H2S) was calculated using the Peng-Robinson equation of states using the molecular density of each specie in the feed and permeate regions as inputs. 24 (3) Note that the cross-sectional area of one plane within the pore A(x) varies along the X direction of the pore because of the roughness of the amorphous silica substrates as well as the non-uniform distribution of benzene inside the pore. Thus, A(x) should be calculated based on the volume available for the various species CH4, CO2, and H2S to occupy.
The results for A(x) as a function of system composition are reported in Figure S7. Knowing the permeability, the transport diffusivity of CH4, CO2, and H2S can be extracted by 26 where Si is the solubility of species i within the pore. The solubility of species i is calculated as the ratio of its density in the pore, ri, and its partial pressure in the bulk feed side, pi, using the slope of the simulated adsorption isotherms (see Figure S2) obtained in the linear regime, at high CH4 mole fraction: We report the solubility for CH4, CO2 and H2S in the pore filled with benzene in Table S1 in the SI. We also computed the self-diffusivity, which describes the motion of individual species within the pores filled with benzene, using the Green-Kubo formulation: 27 In Eq. (6), vk i is the centre of mass velocity of molecule k of species i.
The results of permeability, transport-and self-diffusivity of CH4, CO2 and H2S are reported in the main text. Figure S8. In-plane surface density distributions of CO2 (top) and H2S (bottom) molecules found within first, second, third and fourth layers for the representative 100CO2-1060CH4 and 100H2S-1080CH4 mixtures in the pore filled with benzene at 300 K, respectively.

Molecular Clustering.
To gain insights into the molecular structure of CO2 and H2S molecules inside the pore filled with benzene when injecting CO2-CH4 and H2S-CH4 mixtures, we calculated the in-plane surface density distributions of CO2 and H2S molecules found in the first, second, third and fourth layers formed in the SiO2 pore. The projections are obtained along the X-Y plane. The results are shown for CO2 (top panels) and H2S (bottom) in the representative 100CO2-1060CH4 and 100H2S-1080CH4 mixtures in the pore filled with benzene at 300 K, respectively, in Figure S8. Details regarding the positions of the layers within the nano-pore are presented in Figure S4. The first and fourth layers of CO2/H2S are near the bottom and top silica slabs, respectively, while the second and third layers are located approximately in the middle of the pore. The high-density areas (shown as red-yellow spots) of the contour plots indicate positions where the CO2/H2S molecules preferentially reside. The results show the existence of CO2/H2S molecular clusters in the middle of the pore. These clusters appear to be more pronounced in the presence of H2S within the pores rather than CO2.