The mechanism of fluid migration in porous networks continues to attract great interest. Darcy’s law (phenomenological continuum theory), which is often used to describe macroscopically fluid flow through a porous material, is thought to fail in nano-channels. Transport through heterogeneous and anisotropic systems, characterized by a broad distribution of pores, occurs via a contribution of different transport mechanisms, all of which need to be accounted for. The situation is likely more complicated when immiscible fluid mixtures are present. To generalize the study of fluid transport through a porous network, we developed a stochastic kinetic Monte Carlo (KMC) model. In our lattice model, the pore network is represented as a set of connected finite volumes (voxels), and transport is simulated as a random walk of molecules, which “hop” from voxel to voxel. We simulated fluid transport along an effectively 1D pore and we compared the results to those expected by solving analytically the diffusion equation. The KMC model was then implemented to quantify the transport of methane through hydrated micropores, in which case atomistic molecular dynamic simulation results were reproduced. The model was then used to study flow through pore networks, where it was able to quantify the effect of the pore length and the effect of the network’s connectivity. The results are consistent with experiments but also provide additional physical insights. Extension of the model will be useful to better understand fluid transport in shale rocks.

## I. INTRODUCTION

The economic success related to shale gas production in the United States has generated great interest worldwide. The combination of horizontal drilling and hydraulic fracturing has provided access to large volumes of unconventional oil and gas, which were previously uneconomic to produce,^{1,2} to the point that shale gas has become one of the most important energy resources for the United States.^{3,4} The term shale play is used to describe a geographic area that contains an organic-rich fine-grained sedimentary rock and exhibits low production rates.^{4,5} Shale rocks are mainly composed of kerogen, quartz, clay, carbonates, and pyrite. Secondary components such as uranium, iron, vanadium, nickel, and molybdenum can be found in the shale matrix.^{4,6} Rock mineralogy is essential for shale evaluation and the design of fracture stimulation, as it greatly affects the mechanical properties of the rock. Shales with abundant quartz are usually brittle, while shales with high clay content are ductile.^{7}

Shale formations seem ubiquitous. For example, the Bowland Shale Formation has been identified as the most promising shale gas play in the UK.^{8–10} However, economically producing shale gas in commercially relevant quantities has proven to be highly challenging. One of several unknowns is which rock features contribute the most to gas transport during the various production steps: the microfractures, the existing pores, and/or the hydraulic fractures. Transport of fluids through shale rocks is complicated because of the chemical heterogeneity of the pores, the low conductivity, the lack of a pore connectivity of significant extent, and the pore width often in the nanometer scale.^{11}

Permeability defines the flow capacity of a porous structure.^{12} The dimensions of permeability are *L*^{2} and relate to the cross-sectional area of the pore throats. The true absolute permeability is considered to be an intrinsic property of the porous material. However, in the presence of low pressure gases, the measured permeability might apparently exceed the true absolute permeability, if gas-slippage phenomena are not considered.^{13,14} Typical shale rock permeabilities are in the range of micro-Darcy ($\mu $D) and nano-Darcy (nD),^{15–17} and experimental data show large variations due to the anisotropy of the rocks, applied pressures, etc.^{18}

To predict the movement of a gas through a porous medium and calculate the matrix permeability, a transport model needs to be defined.^{19} Darcy’s law can describe the laminar viscous flow of one or more phases in a porous matrix, when the rate of the flow is proportional to the pressure gradient. This law for incompressible fluids reads^{20}

where *Q* is the volume flow rate, *A* is the area of porous media normal to the flow, *K* is the hydraulic conductivity (permeability), *h* is the pressure head (pressure divided by the specific weight), *z* is the elevation, and *L* is the length of the flow path. Subscripts 1 and 2 designate the up and downstream positions, respectively.

Because of the unique and complex structure of unconventional reservoirs, flow mechanisms can deviate from Darcy’s law.^{21} Many hypotheses have been proposed to identify the main characteristics responsible for this deviation. One of the prevailing hypotheses suggests the coexistence of different flow regimes due to the various pore sizes in the matrix.^{22,23} Other hypotheses consider the chemical composition of the matrix and the anisotropic permeability to be equally important factors,^{19} and the fact that the behaviour of fluids can vary significantly under confinement.^{24} Indeed, as the distance between the pore walls decreases, the interactions between the molecules and the surface become stronger and the journey one molecule can travel freely shortens. According to the kinetic theory, the average distance a molecule can travel before it collides with another molecule, or with the pore wall, is defined as the mean free path $\lambda $.^{19,22} The mean free path is proportional to the temperature and inversely proportional to the pressure in large pores. In small pores, the mean free path depends on the pore size and shape.^{25}

The Knudsen number, expressed as the ratio of the mean free path $\lambda $ and the characteristic length of the system $\Lambda $, can be used to identify the flow regime.^{26} Four possible transport mechanisms are identified: the Knudsen diffusion, the slip flow, the transitional flow and the viscous flow.^{22} Darcy’s law can be applied in the case of viscous flows, as it inherently assumes laminar flow (Reynolds number = 1).^{27} The transport mechanism also depends on the confinement, as confinement can promote capillary condensation. When the pore pressure exceeds the critical pressure required to fill the pores, the adsorbed molecules condense. This phase change delays the transport of vapour through the condensate-filled pores.^{22,23}

To interpret kinetic data in heterogeneous and anisotropic systems, it is necessary to use complicated models that account for micropore, mesopore, and macropore diffusional resistances.^{28} Different experimental, theoretical, and computational methods have been suggested in the literature. Each approach exhibits advantages and limitations, depending on the time scale of the analysis, and the size of the samples considered.^{29,30} However, by following an experimental approach, no distinction can be made between diffusion and convection and hence, the mechanism of fluid transport cannot be identified with certainty.^{31}

One approach implements empirical equations derived from Darcy’s law. Within this approach, Darcy’s equation is enriched with coefficients designed to match experimental data for systems that, e.g., exhibit slip flow. Most commonly, these models take into consideration effective stress and slip flow contributions, in order to model the apparent permeability. These equations can provide very useful insights, but their accuracy is limited to the experimental conditions used for their derivation.^{32–34}

Lattice Boltzmann (LB) simulations can be applied to obtain transport properties. LB is a mesoscopic approach used to simulate fluid transport, under the assumption of continuum flow. Software packages that implement LB techniques continue to gain attention as computational power increases. However, LB simulations do not account for the effect of micropores and tend to overestimate transport properties.^{35} Another approach involves the use of molecular dynamics (MD) simulations, which can study fluid transport in micropores at the expense of high computational effort, which leads to the so-called scale problem.^{36}

Kinetic Monte Carlo (KMC) simulations employ stochastic descriptions of the dynamical phenomena governing various processes in catalysis, material fabrication, defect evolution in crystals, and diffusion.^{37–41} KMC methods enable access to long time scales (from ms to h) and large spatial scales (from nm to $\mu $m) at comparatively low computational expense.^{42,43} In principle, KMC is able to describe the exact dynamic evolution of a system, by considering that the system performs state-to-state jumps whose statistics follow known distributions. This information determines the timing of the jumps and also the final state (given an initial state). KMC exhibits significant computational savings, compared for instance to MD, as it coarse-grains the vibrational trajectories.^{44} The waiting time before an event takes place depends on the energy barrier the system has to surmount in order to get from one energy basin (state) to another.^{44}

Considering that permeability depends on a number of factors, such as pore characteristics, chemical composition, and transport mechanisms, an approach that accounts for all these factors at low computational cost is required. In this work, we are interested in investigating how the pore characteristics and the pore network connectivity affect the transport properties of light hydrocarbons in hydrated nanopores with different chemical compositions by performing KMC simulations. Methane is selected as it comprises the main component of natural gas. Phan *et al.*^{45} reported the transport properties of methane in three different 1 nm wide slit-shaped pores filled with water. Models for silica, magnesium oxide, and alumina were used as solid substrates, and the diffusion coefficient of methane inside the pores was obtained by performing MD simulations^{45} in the canonical ensemble at 300 K. The three pores filled with water are considered representative of minerals found in the subsurface.^{45} We perform here KMC simulations for the systems considered by Phan *et al.* to validate our model. Then, we use our validated KMC model to investigate the effect of the pore length and the pore network connectivity on methane transport.

The remainder of the paper is organised as follows. Section II provides the theoretical background under which the KMC algorithm is implemented. In Sec. III, we describe the physical systems considered by Phan *et al.* and how a KMC lattice is constructed to represent them, followed by the methodology used to calculate the KMC rates. In Sec. IV, we report flux, permeability, and Mean First Passage Time (MFPT) results through a few single pores. The effect of pore length and network connectivity is also investigated. Finally, Sec. V summarises our conclusions and suggests future studies.

## II. THEORETICAL BACKGROUND

KMC seeks to use state-to-state transition rates to simulate trajectories of the stochastic “wandering” of a system around the state space.^{43} To this end, the KMC implementation requires rate constants that capture the probability per unit time of such state-to-state transitions.^{43} A sequence of such transitions constitutes a sample path or trajectory, whose statistics follow the so-called Master equation (M-equation) that governs the dynamics of the system. The M-equation gives the time course of the probability mass function (p.m.f) of finding the system in a certain state at a given time^{46}

The M-equation describes the system evolution as a balance of probabilities [*P*_{i}(*t*) and *P*_{j}(*t*)] multiplied by the kinetic constants *k*_{ij} and *k*_{ji}, respectively, which quantify the propensity of an event, or, in other words, its probability of occurrence per unit time. Despite its simplicity, the M-equation can hardly be solved. Instead, a KMC algorithm can be employed to simulate sample paths (trajectories) and estimate statistical properties of interest.^{43}

A KMC algorithm can be implemented to address both surface (solid state) and bulk diffusion problems. The mathematical basis of the M-equation is the same but conceptual differences exist. Solid state problems can be described by a state vector *x* and a time coordinate *t* where, *x* represents local minima on a potential energy surface (PES) and the M-equation describes the transition from one such local PES minimum to another.^{47} Diffusion is an activated process and the kinetic constants that describe such transitions depend on the energy barrier between two states. In the case of bulk diffusion problems, the state vector *x* represents the population of species inside a region of a certain energy level. The kinetic constants depend on the diffusivity of the species considered.^{48}

To ensure that the selected transition rates capture correctly the dynamical evolution of the system towards equilibrium, one can consider the Boltzmann relationship

Microscopic reversibility (i.e., detailed balance) needs to be satisfied for any connected pair of states *i* and *j*. This implies the following relationship between rates *k*_{ij} (forward process) and *k*_{ji} (reverse process):^{43}

The features of the KMC algorithm were first described by Young and Elcock in 1966.^{49} Early implementations of the KMC algorithm were reported by Bortz *et al.*^{50} and Gillespie, who utilised a rejection-free method known as the direct method.^{51} We implement here the direct method to describe diffusion as a function of jumps between neighbouring voxels. We discretize the sample space and create finite domains, subsets of the sample space, in which a molecule can be found. One can determine the probability of occupancy inside these subsets, which are referred to as voxels. The movement of molecules from one voxel to others is treated as possible transition pathways. We first identify all possible pathways and then we determine the rates through which the transitions take place. The possible pathways a molecule can follow while it diffuses in a 1D domain with periodic boundaries and *N* square voxels is 2*N*. If the system is closed and both the boundaries are reflective, the possible pathways are reduced to 2*N* − 2 (see Fig. 1).

The algorithm requires the selection of two random, uniformly distributed numbers *u*_{1} and $u2\u2009\u2208\u2009[0,1]$. We use the Mersenne Twister MT19937 uniform random number generator (u.r.n.g).^{52} One of the selected numbers, *u*_{1}, is used to select which event takes place; the second, *u*_{2}, is used to determine the time increment due to the selected transition. A detailed description of the direct method^{51,53} is provided in the supplementary material.

## III. SIMULATION MODELS AND METHODOLOGY

### A. Model system

Each system considered here comprises of three regions (region 1, region 2, and region 3), two fluid species (methane and water), and two fluid phases. Each region is described as a collection of different voxels, as discussed later. In all cases, two bulk areas surround the pore (region 1), two layers of liquid water lay outside both sides of the pore (region 2) and a slit pore is filled with water (region 3). Water molecules are in liquid phase, while methane molecules are in gas phase. As methane migrates through the pore, it first solvates in the liquid water and then it diffuses. Interfaces separate the bulk and the pore space. A schematic is shown in Fig. 2, which corresponds to a silica slit-pore filled with water; methane molecules are present at both bulk regions. In the KMC lattice, each region is composed of voxels. Each voxel is assigned a forward and a backward diffusion rate. In Fig. 2, the transition rates at the interfaces are presented schematically.

In our KMC model, we need to determine the diffusion rates within the three regions and the transition rates at the interfaces between them. To determine the diffusion rates, we considered the diffusion coefficient of methane in the water layer outside the pores and that of methane inside the hydrated pores, enabling us to simulate jumps across voxels inside each region. Potential of mean force (PMF) profiles and methane density profiles was used to assign the transition rates at the interfaces.

The rates that describe the back and forth transitions between neighbouring voxels inside each individual region will be equal due to microscopic reversibility. To calculate such rates, *k* inside the hydrated pore and the water layer, we use the following expression:^{53}

where *D* is the coefficient of methane and l is the voxel size. The diffusion coefficients of methane in the various hydrated pores were obtained by Phan *et al.* and are summarized in Table I. The diffusion coefficient of methane inside the water layer is set to 1.8 $\xd7$ 10^{−9} m^{2}/s.^{45}

Material . | D_{t} ($\xd7$10^{10}, m^{2}/s)
. | $\Delta $F (kcal/mol) . |
---|---|---|

SiO_{2} | 7.82 | 1.2 |

MgO | 5.51 | 1.6 |

Al_{2}O_{3} | 3.26 | 1.3 |

Material . | D_{t} ($\xd7$10^{10}, m^{2}/s)
. | $\Delta $F (kcal/mol) . |
---|---|---|

SiO_{2} | 7.82 | 1.2 |

MgO | 5.51 | 1.6 |

Al_{2}O_{3} | 3.26 | 1.3 |

Information regarding the PMF experienced by the methane molecules across the hydrated pores was obtained by Phan *et al.*^{45} The methane molecules diffuse first through the layer of water outside the pore, where the water molecules are able to move freely. Once they enter the hydrated pores, the adsorbed water molecules can hinder methane transport. To calculate the transition rate at the interfaces, we make use of the energy barriers obtained from the PMF profiles^{45} and then use Boltzmann’s distribution to define the probability of occupancy in regions 2 and 3. The energy barriers considered in our calculations correspond to the maximum height (difference between maxima and minima) of two subsequent PMF curves and are reported in Table I. These energy barriers represent the transport barrier inside the pores. The probability of occupancy in regions 2 and 3, resulting from these energy barriers, is reported in the supplementary material.

From the density profiles of methane and water, we counted the methane population inside the water layer and the gaseous phase and we extracted the probability of occupancy in regions 1 and 2, (reported in the supplementary material). The probabilities obtained from the PMF and density profiles were normalised and used to describe the transition rates at the interfaces. In Table II, we summarise the probability of occupancy at regions 1 (*p*_{1}), 2 (*p*_{2}), and 3 (*p*_{3}) for methane molecules in each system. The transition rates $k1,2int$, $k2,1int$, $k2,3int$, and $k3,2int$ are calculated as

Material . | p_{1}
. | p_{2}
. | p_{3}
. |
---|---|---|---|

SiO_{2} | 0.709 | 0.270 | 0.020 |

MgO | 0.760 | 0.227 | 0.012 |

Al_{2}O_{3} | 0.748 | 0.227 | 0.024 |

Material . | p_{1}
. | p_{2}
. | p_{3}
. |
---|---|---|---|

SiO_{2} | 0.709 | 0.270 | 0.020 |

MgO | 0.760 | 0.227 | 0.012 |

Al_{2}O_{3} | 0.748 | 0.227 | 0.024 |

As a final step, the KMC lattice needs to be constructed. The information available includes methane diffusivities, PMF, and density profiles along the direction parallel to the pore. Because water is considered stagnant, we create a lattice and consider only one type of species, methane, effectively undergoing 1D transport along the direction parallel to the pore. We define the size and the number of the voxels. Each voxel must be bigger than the mean free path of the molecules. Each region contains a different number of voxels. We represent region 1 by a single voxel, assuming the gaseous methane to be well mixed. The methane particles move from one voxel to a neighbouring voxel, located along the x-direction, following a 1D trajectory. An example of the lattice used to perform these calculations (in this regard for the silicon oxide pore) is presented in Fig. 3. Information regarding the lattice implemented for all three substrates is shown in Table III.

Material . | SiO_{2}
. | MgO . | Al_{2}O_{3}
. |
---|---|---|---|

Region 1 dimensions (x,y,z) (nm) | 1.6 $\xd7$ 0.4 $\xd7$ 0.4 | 2.12 $\xd7$ 0.42 $\xd7$ 0.42 | 2.30 $\xd7$ 0.46 $\xd7$ 0.46 |

Region 2 dimensions (x,y,z) (nm) | 1.6 $\xd7$ 0.4 $\xd7$ 0.4 | 2.12 $\xd7$ 0.42 $\xd7$ 0.42 | 2.30 $\xd7$ 0.46 $\xd7$ 0.46 |

Region 3 dimensions (x,y,z) (nm) | 5.2 $\xd7$ 0.4 $\xd7$ 0.4 | 5.09 $\xd7$ 0.42 $\xd7$ 0.42 | 4.60 $\xd7$ 0.46 $\xd7$ 0.46 |

Voxel size | 0.4 $\xd7$ 0.4 $\xd7$ 0.4 | 0.42 $\xd7$ 0.42 $\xd7$ 0.42 | 0.46 $\xd7$ 0.46 $\xd7$ 0.46 |

Number of voxels in region 1 | 1 | 1 | 1 |

Number of voxels in region 2 | 4 | 5 | 5 |

Number of voxels in region 3 | 13 | 12 | 10 |

Number of particles inserted in region 1 | 3 | 8 | 8 |

Material . | SiO_{2}
. | MgO . | Al_{2}O_{3}
. |
---|---|---|---|

Region 1 dimensions (x,y,z) (nm) | 1.6 $\xd7$ 0.4 $\xd7$ 0.4 | 2.12 $\xd7$ 0.42 $\xd7$ 0.42 | 2.30 $\xd7$ 0.46 $\xd7$ 0.46 |

Region 2 dimensions (x,y,z) (nm) | 1.6 $\xd7$ 0.4 $\xd7$ 0.4 | 2.12 $\xd7$ 0.42 $\xd7$ 0.42 | 2.30 $\xd7$ 0.46 $\xd7$ 0.46 |

Region 3 dimensions (x,y,z) (nm) | 5.2 $\xd7$ 0.4 $\xd7$ 0.4 | 5.09 $\xd7$ 0.42 $\xd7$ 0.42 | 4.60 $\xd7$ 0.46 $\xd7$ 0.46 |

Voxel size | 0.4 $\xd7$ 0.4 $\xd7$ 0.4 | 0.42 $\xd7$ 0.42 $\xd7$ 0.42 | 0.46 $\xd7$ 0.46 $\xd7$ 0.46 |

Number of voxels in region 1 | 1 | 1 | 1 |

Number of voxels in region 2 | 4 | 5 | 5 |

Number of voxels in region 3 | 13 | 12 | 10 |

Number of particles inserted in region 1 | 3 | 8 | 8 |

### B. Flux calculations

Methane molecules in gas phase were inserted in one of the bulk areas (feed), resulting in a pressure rise. The other bulk area (permeate) was kept empty throughout the simulations. Once the methane molecules cross the pore and enter the permeate region, they are deleted and immediately added back into the feed region. Hence, the pressure drop across the membranes was maintained constant. The boundaries of the system were reflective. We calculate the cumulative number of particles that cross the permeate region every 30 ns and report averages obtained from 100 independent simulations. Replicating the simulations as such results in smooth profiles with minimum fluctuations. Starting from *t* = 0, the system was allowed to progress for 30 ns, during which time a counter reports the number of molecules crossing the pore. After 100 simulations, the system was left to diffuse for another 30 ns. This procedure was repeated until a total simulation time of 720 ns was reached. The flux was determined by counting the number of molecules crossing the pore (*Q*_{t}) over time

In Eq. (10), *J* is the molar flux of methane, $\Delta Qt/\Delta t$ is the slope of the fitted line (cumulative number of molecules vs. time), and *A* is the cross-sectional area available for gas permeation, perpendicular to the direction of the diffusion.

The plots of the cumulative number of molecules as a function of time obtained for every system are presented in the supplementary material. A straight line is fitted to this data to calculate the slope of the linear plots. From the slope, Eq. (10) is used to calculate the flux of methane molecules.

### C. Permeability calculation

The substrate permeability *K* was calculated as

where *J* is the molar flux of methane from the KMC calculations, *l*_{p} is the length of the pore, *p*_{1} is the pressure applied in the feed phase, *p*_{2} is the pressure applied in the permeate region (*p*_{2} = 0). In this calculation, the pressure drop remains constant. At most 5% of the molecules in the feed volume escaped to the pore voxels (i.e., region 2).

### D. Mean first passage time (MFPT)

The MFPT reflects the average time required for a particle to reach a specific position inside the pore for the first time. For MFPT calculations, silicon, magnesium, and aluminium oxide pores 40 nm long were considered. A layer of water, approximately 1.5 nm thick, surrounded each side of the hydrated pores.

In the beginning of the KMC simulation, a methane molecule was placed in the middle of the pore (voxel No. 86), as shown in Fig. 4. We performed 1000 independent simulations and monitored the molecule’s trajectory as it exited the pore. The KMC lattice used for these calculations was symmetric and we focused our analysis on the *x*_{a} direction. Starting from voxel No. 86, we recorded the first time at which the molecule was found at each neighbouring voxel sequentially (i.e., voxel No. 85, voxel No. 84, etc.), until it reached voxel No. 1. By averaging these “first times,” we obtained the MFPT as a function of distance from the center of the pore. Whenever the molecule attempted to cross the center of the pore, by moving into voxel No. 87, the simulation was terminated and discarded. The size of the voxels used to simulate the layer of water (region 2) had dimensions 0.25 (nm) × 0.5 (nm) × 0.5 (nm). The voxels used to represent the pore (region 3) had dimensions 0.5 (nm) × 0.5 (nm) × 0.5 (nm).

### E. Heterogeneous pore networks

Our KMC model was used to investigate the transport properties in heterogeneous pore networks. We assessed how the flux observed through the membranes is affected by the length and the connectivity of the pores. To understand how the pore size affects molecular flux, a lattice similar to the one described in Fig. 3 is used. The size of the water layer formed outside of the pores is kept the same. To obtain better statistics, the density of methane particles in the bulk region is increased. The particle flux is expected to gradually decrease, as the length of the pore is increased, while keeping the pressure drop constant.

The effect of the pore connectivity on flux was quantified, starting by a long pore surrounded by water and two bulk regions containing methane (as in Fig. 2). As a next step, two disconnected pores with half the length of the original pore are used in place of the single long pore. These pores are also surrounded by bulk regions of water layers. Afterwards, the length of the pores is reduced by half, yielding four pores, also disconnected. The pore networks studied (type 1, type 2, and type 3) are illustrated in Fig. 5.

## IV. RESULTS AND DISCUSSION

### A. Validation against analytical models

We validated the accuracy of our algorithm by comparing its predictions to those expected by solving the diffusion equation analytically. We consider a finite domain with total length L and reflective boundaries. A non-trivial solution of the 1D diffusion equation on the x-direction is^{54}

where *C*(*x*, *t*) is the spatially and temporally varying concentration, *C*_{0} is the initial concentration of molecules, and $\u2212h\u2264x\u2264+h$ is the region in which the molecules were initially confined. The validation of the KMC algorithm against the 1D diffusion equation was satisfactory, as shown in the supplementary material.

### B. Flux and permeability through single pores

For validation purposes, we compare the values obtained by our model against those reported by Phan *et al.*^{45} The results are presented in Table IV.

Material . | Flux J (mol/m^{2}s)
. | Permeability K (mol/ms MPa)
. |
---|---|---|

SiO_{2} | 74.76 $\xb1$ 1.42 | (5.38 $\xb1$ 0.11) $\xd7\u200910\u22129$ |

MgO | 71.83 $\xb1$ 1.26 | (3.05 $\xb1$ 0.05) $\xd7\u200910\u22129$ |

Al_{2}O_{3} | 75.91 $\xb1$ 1.27 | (2.04$\xb1$ 0.03) $\xd7\u200910\u22129$ |

Material . | Flux J (mol/m^{2}s)
. | Permeability K (mol/ms MPa)
. |
---|---|---|

SiO_{2} | 74.76 $\xb1$ 1.42 | (5.38 $\xb1$ 0.11) $\xd7\u200910\u22129$ |

MgO | 71.83 $\xb1$ 1.26 | (3.05 $\xb1$ 0.05) $\xd7\u200910\u22129$ |

Al_{2}O_{3} | 75.91 $\xb1$ 1.27 | (2.04$\xb1$ 0.03) $\xd7\u200910\u22129$ |

The flux observed though the three pores is similar. However, the permeability varies. The silicon oxide pore is the most permeable, followed by the magnesium oxide one and the alumina oxide. These trends are in quantitative agreement with the MD results.^{45}

Performing 100 independent KMC simulation runs for the silica, magnesium oxide, and aluminium oxide pores required just 32.9 s, 63.6 s, and 83.1 s, respectively. Comparison between KMC and MD results is presented in Fig. 6. To calculate the error bars in Fig. 6, we perform 100 independent simulations and calculate the standard deviation for flux and permeability. We divide the standard deviation of the mean by the square root of the sample size to obtain the standard error of the mean, which we report here.

### C. Mean first passage times

The MFPT profiles shown in Fig. 7 were obtained for the three substrates by implementing the KMC algorithm. As expected from the diffusion coefficients that describe methane transport inside the slit-shaped pores, transport inside the silicon oxide pores is the fastest. Methane moves slowly inside the alumina oxide pores, for which the MFPT is almost three times slower than that obtained in the silicon oxide pore. The results obtained from the KMC model are in agreement with the ones reported by Phan *et al.*^{45} as shown in the supplementary material. It is worth repeating that the computational effort associated with the KMC approach is significantly smaller than that of the MD approach. The CPU times required were 48.9 s, 48.8 s, and 51.3 s for the simulations performed on the silica, magnesium, and aluminium oxide pores, respectively, when the KMC model was implemented.

### D. Flux versus pore length

Since the KMC model quantitatively reproduces the results obtained from the MD simulations, it can be used to quantify how certain pore characteristics affect permeability. To elucidate the effect of the pore length on the observed flux, we conducted systematic KMC simulations. To improve the statistics, we increased the number of gas molecules inserted in the feed area to 100. The results are shown in Fig. 8 (left panel) where the calculated flux is plotted against the pore length. For these calculations, the pore length was increased from 4 to 25 nm for all three substrates. As expected, as the pore length increases, the flux decreases. We expect the rate of flux decrease to be similar for all substrates. To investigate this hypothesis, we calculate the % of the flux decrease and plot the results obtained against the pore length. The results shown in Fig. 8 (right panel) are indicative of the qualitative validity of our model.

To quantitatively prove the accuracy of our model, we rely on Eq. (13), an analytical expression that describes the relationship between the observed flux and the membrane thickness,^{55}

In Eq. (13), *L* is the membrane thickness, *L*_{m}(*L*) is the mass transport coefficient, *x*_{b} is the mole fraction at the feed-membrane interface, and *x*_{p} is the mole fraction at the membrane-permeate interface. As the pore length increases some of the molecules will occupy positions inside the membrane thus decreasing the driving force for transport. To balance this effect, we monitor the population of methane molecules inside the pore during the simulation and we adjust the number of feed molecules as required to maintain a constant pressure drop. This adjustment keeps the *x*_{b} and *x*_{p} mole fractions in Eq. (13) constant.

We assume the mass transfer coefficient to linearly decrease as the pore length increases since the conditions under which our simulations take place are kept constant. The relationship between the mass transfer coefficient and the pore length can then be described as,

### E. Effect of pore connectivity

Changing the pore network connectivity, we modify the number of interfaces present in the system (see diagrams in Fig. 5). Assuming that the resistance to diffusion is due to the pore length and the interfaces, we seek to investigate whether these two effects equally hinder the gas transport or whether we reach a rate-limiting step.

The effect of pore connectivity on the overall pore network’s flux was quantified by simulating three 1D pore networks, shown in Fig. 5. The first one consisted of a single long pore (4 nm). This pore was cut into half and the newly generated pores (2 nm) were surrounded by water layers and vapour bulk areas. These two pores were then cut into half to generate the third system. The flux observed for the three pore networks is reported in Fig. 10. The results show that the fluxes in each of the 1D pore networks decrease as the number of interfaces increases, suggesting that the rate-limiting step in these pore networks is due to the entrance of methane into the water-filled pores.

Figures 8 and 10 show that the pore length and the pore network connectivity are equally important in determining the flux. The observed flux decrease is similar whether the pore length increases by 1 nm (from 4 nm to 5 nm and from 5 nm to 6 nm) or the pore network becomes more disconnected (from type 1 to type 2 and from type 2 to type 3). We point out that the percent decrease in flux due to the added interfaces follows a similar trend for all three substrates.

## V. CONCLUSIONS

We have developed a lattice-based KMC model to study fluid transport through slit-shaped pores with different chemical compositions. The substrates analysed here are meant to represent some of the main components of the inorganic material found in shale rocks. The proposed model was used to quantify how the presence of natural fractures and disconnected nanopores affects the permeability of a heterogeneous domain.

The chemistry of the pores affects the transport behaviour of gas methane molecules, as reported in previous studies. The hydrated silicon oxide nanopores exhibit the highest permeability, followed by the permeability observed in hydrated magnesium and aluminium oxide pores. The agreement between the KMC model and MD simulations is quantitative; however, the computational times are significantly reduced when using the KMC model. The model was then used to provide insights regarding the contribution of the pore network characteristics in the transport behaviour. It was found that both the pore length and the network connectivity play a significant role on gas migration.

From our simulations, we observed that the chemical composition of the substrates affects the absolute values of the flux observed for all pore lengths and type of networks, but it does not affect the rate at which fluxes decrease. Hence, one could generate multiple realizations of possible networks for one substrate and extrapolate the findings to approximate the behaviour of substrates with different chemical compositions.

Our model can be considered as a bottom-up approach for mesoscopic studies. Any type of designed or natural network can be simulated, as long as the kinetic (diffusion constants) and thermodynamic (barriers due to the interfaces) properties are provided. The kinetic and thermodynamic properties used for these studies, together with the methodology followed for the calculation of the transition rates, can be applied to any type of pore networks and provide insights regarding the significance of the chemical composition and pore features to the resulting transport properties, provided gas transport occurs via diffusion. Such studies will contribute to a better understanding of the diffusion of multiple species encountered in shale rocks and potentially assist the formulation of strategies to maximise the natural gas or oil recovery through shale rocks.

## SUPPLEMENTARY MATERIAL

See supplementary material for the description of the KMC algorithm and the KMC validation against the 1D diffusion equation and the MD simulations. Density profiles (methane and water) used to calculate the KMC rates at the interfaces are reported together with the probabilities obtained for the three regions.

## ACKNOWLEDGMENTS

We would like to thank financial support provided by Halliburton. The authors acknowledge the use of the UCL Legion High Performance Computing Facility (Legion@UCL), and associated support services, in the completion of this work. This research also received funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 640979 and from the Marie Curie Career Integration Grant No. 2013-CIG-631435 to A.S.