A Novel Modeling Approach to Stochastically Evaluate the Impact of Pore Network Geometry, Chemistry and Topology on Fluid Transport

Fine-grained sandstones, siltstones, and shales have become increasingly important to satisfy the ever-growing global energy demands. Of particular current interest are shale rocks, which are mudstones made up of organic and inorganic constituents of varying pore sizes. These materials exhibit high heterogeneity, low porosity, varying chemical composition and low pore connectivity. Due to the complexity and the importance of such materials, many experimental, theoretical and computational efforts have attempted to quantify the impact of rock features on fluids diffusivity and ultimately on permeability. In this study, we introduce a stochastic kinetic Monte Carlo approach developed to simulate fluid transport. The features of this approach allow us to discuss the applicability of 2D vs 3D models for the calculation of transport properties. It is found that a successful model should consider realistic 3D pore networks consisting of pore bodies that communicate via pore throats, which however requires a prohibitive amount of computational resources. To overcome current limitations, we present a rigorous protocol to stochastically generate synthetic 3D pore networks in which pore features can be isolated and varied systematically and individually. These synthetic networks do not correspond to real sample scenarios but are crucial to achieve a systematic evaluation of the pore features on the transport properties. Using this protocol, we quantify the contribution of the pore network’s connectivity, porosity, mineralogy, and pore throat width distribution on the diffusivity of supercritical methane. A sensitivity analysis is conducted to rank the significance of the various network features on methane diffusivity. Connectivity is found to be the most important descriptor, followed by pore throat width distribution and porosity. Based on such insights, recommendations are provided on possible technological approaches to enhance fluid transport through shale rocks and equally complex pore networks. The purpose of this work is to identify the significance of various pore network characteristics using a stochastic KMC algorithm to simulate the transport of fluids. Our findings could be relevant for applications that make use of porous media, ranging from catalysis to radioactive waste management, and from environmental remediation to shale gas production.


Introduction
In view of the growing demand for energy, production rocks that were once considered non-reservoirs, such as fine-grained sandstones, siltstones and shales, are becoming increasingly important. Meanwhile, technological advances in horizontal drilling and hydraulic fracturing have permitted dramatic increases in hydrocarbon (gas and oil) production from shale formations (BP 2019; Energy Information Administration (EIA) 2014; Jarvie 2012). It is known that fluid flow and mass transport through rocks are controlled by pore network characteristics such as pore size distribution, pore connectivity, porosity, pore throat width distribution and mineralogy (Bear 1972;Dullien 1992;Chalmers et al. 2012;Kwon et al. 2004). Current research efforts focus on determining which, out of these and other features, control fluid migration. For example, Sahimi described pore connectivity by an equivalent network of pore throats, i.e., narrow passages through which fluids flow, and pore bodies, large voids that meet through the throats (Sahimi 2011). This feature of the porous media, frequently referred to as topology, pore throat connectivity or coordination number, has been widely identified as one of the most important parameters that affect transport in many porous media, including shale rocks (Chen et al. 2003;Vasilyev et al. 2012). According to Rabbani et al. (2016), connectivity can be measured following either forward or backward methods. Forward methods mainly involve 2D or 3D image acquisition of core samples, followed by sophisticated morphological analysis and image processing to extract the network's coordination number (Dong and Blunt 2009;Peng et al. 2012). Backward methods first quantify some macroscopic characteristics of the porous media, e.g., capillary pressure, porosity or relative permeability, and then back-calculate the features of the underlying pore network model (Mata et al. 2001;Tsakiroglou and Payatakes 2000;Davudov et al. 2019;King et al. 2015;Chalmers et al. 2012).
When applied to Barnett and Haynesville shale rocks samples, such approaches revealed pore networks with low connectivity (Hu et al. 2015). Similar results have been reported for UK samples from the East Irish Sea Basin in northern England, and it is now generally assumed that most shale formations consist of poorly connected pore networks (Oluwadebi et al. 2019). Moreover, it is common that during hydraulic fracturing, due to the large pressure gradients applied and in situ stress alterations, part of these poorly connected networks crush and the connectivity further reduces (King et al. 2019;Yu et al. 2018). To support this possibility, Davudov and Moghanloo (2018) identified connectivity loss as the main mechanism leading to permeability reduction in shale samples. From a modeling perspective, the effect of pore network connectivity on permeability can be quantified by modifying selected parameters within the models used. Few examples include Civan's (2001) modified Kozeny-Carman (KC) permeability model, and Pape's et al. (2000) model, which is derived from fractal theory.
In addition to connectivity, other materials properties are important in this quest, for example, the porosity of a medium, defined as the volume fraction of its voids (Thakur 2017). Porosity affects permeability so strongly (Akintunde et al. 2014;Akinlotan 2016) that Archie proposed a porosity-permeability relationship, derived for sandstones, limestones and muddy sands, which predicts a tenfold increase in permeability when porosity increases by merely 3% (Archie 1950). Because it has been suggested that porosity intrinsically depends on the microstructure of porous materials (Shen et al. 2018), we performed a sensitivity analysis on the effect of macro-and microporosity on permeability (Apostolopoulou et al. 2019). In that prior study, we also compared the reliability of our stochastic KMC approach against analytical models to predict permeability.
A plethora of published studies also attempted to correlate the porosity of a formation to its permeability. Magara, e.g., proposed a permeability-porosity relationship to fit laboratory-based porosity-permeability data (Magara 1978). Additionally, several modifications to the classic fractal theory (Mandelbrot 1985) have been proposed to better capture this relationship (Katz and Thompson 1985;Yu and Li 2001;Rieu and Sposito 1991). Zhang et al. (2008) utilized the poroelasticity theory, and Javadpour proposed a dynamic model to describe dynamic porosity and apparent permeability, taking into consideration the poromechanical processes that occur during reservoir stimulation (Sheng et al. 2019). These authors also correlated permeability, porosity and formation depth using neural networks (Jamialahmadi and Javadpour 2000). These efforts are justified by the fact that understanding the effect of porosity on permeability is crucial in developing stimulation strategies for a given formation and for evaluating the profitability of shale formations.
Permeability in low-permeability reservoirs can be affected also by microfractures, which can be naturally occurring or artificially generated via hydraulic fracturing (Cortez-Montalvo et al. 2018). In our previous work, we considered heterogeneous networks containing highly conductive microfractures and anisotropic low-permeability channels (Apostolopoulou et al. 2019;Inyang et al. 2019). In those studies, focused on 2D networks, we provided a detailed description of how to convert heterogeneous matrix fabrics into 2D KMC rates (Apostolopoulou et al. 2019;Inyang et al. 2019). Our stochastic KMC model was validated, in 3D, considering heterogeneous hydrated slit pores (Apostolopoulou et al. 2019). A cross-sectional analysis of the fluid behavior in such pores revealed regions with varying diffusivity, almost by 1 order of magnitude (Apostolopoulou et al. 2019). The algorithm to extract 3D KMC input parameters (i.e., diffusion rates) is described in this prior work (Apostolopoulou et al. 2019).
The sizes of pores and pore throats are also certainly important. In conventional reservoir rocks, pore sizes and throats are relatively large, sufficiently so to store and deliver economic quantities of petroleum. On the contrary, pore throats in unconventional reservoirs are so small that they can limit hydrocarbons passage (Nelson 2009). Nelson reported that the experimentally measured pore throat sizes in samples from Devonian, Jurassic, Cretaceous, Pennsylvanian and Pliocene shale formations range from 5 to 50 nm, although in some samples these values can reach 100 nm (Nelson 2009). Other experimental studies also confirm this assessment (Hu et al. 2015). Zhang et al. (2018) evaluated pore size distribution (PSD) and permeability of a tight reservoir, located in China, and reported that throat sizes are slightly positively correlated with permeability. On the other hand, Katsube and Issler (1998) found that the mode of the PSD significantly correlates with the permeability for various Canadian shale samples. However, even these groups acknowledge that other material properties are likely to affect the relationships reported.
Another material property that is often considered important in determining fluid transport is the chemistry of the porous medium. This is quantified in terms of the mineralogy of the pore throats. It has been observed that in narrow, single pores, the pore throat surface interacts with the fluids, in some cases yielding kinetic barriers that hinder transport (Wang et al. 2018;Franco et al. 2016). This is more pronounced when the throat surfaces are polar and water is present, especially in the case of two-phase systems (Bui et al. 2017;Phan et al. 2016a, b). We recently implemented a stochastic approach to generate digital libraries correlating gas diffusivity to the pore size for five solid substrates (Apostolopoulou et al. 2019). In single slit-shaped pores, it was found that the solid substrate chemistry has little effect on gas diffusivity for pore widths larger than ~ 30 nm, although the effect can be large in smaller pores.
Although the importance of the pore network features has been extensively discussed in the literature, the classification of these features in order of importance remains elusive. The purpose of this work is to identify the significance of various pore network characteristics using a stochastic KMC algorithm to simulate the transport of fluids. A rigorous systematic approach is employed here to generate 3D pore networks in which connectivity, porosity, PSD and mineralogy are treated as variables and are altered individually to quantify their impact on the transport of supercritical methane. Our main findings show connectivity to be the most important parameter that affects the diffusivity of rarefied fluids followed by the pore throat widths and finally porosity. The remainder of this article is organized as follows. In Sect. 2, we describe the protocols implemented to generate the pore networks, the systems investigated and the stochastic 3D KMC model used to model supercritical methane transport. In Sect. 3, we provide the results obtained when the pore network features are altered separately; we also present the sensitivity analysis to reveal the relative importance of the pore networks characteristics. Finally, in Sect. 4, we summarize our observations and provide recommendations. Several computational details and a table summarizing the symbols used throughout the manuscript are provided as Supplementary Material.

Pore Network Features
As discussed above, we seek to generate pore networks whose features are changed systematically. We describe, in what follows, how pore connectivity, network porosity, pore throat size distribution and pore throat chemistry are changed systematically in our 3D models. In our pore networks, the pores are of the same size and are represented as binary variables; each voxel can either have a single pore or none. Similarly, pore connectivity is represented as a binary variable for a set of adjacent pores, i.e., the neighboring voxels i and j can be either connected via a pore throat or not. The number of connections a voxel has with its neighboring 6 voxels defines the connectivity, which is an integer variable with values between 0 and 6. The pore throats are slit shaped and have the same length and different widths. We implement a stochastic method for generating the pore networks; using distributions of the various pore network features allows the models developed to exhibit local variations of the various properties, as frequently observed in actual samples. Nevertheless, the aim of this study is not to propose a new pore network model, but to use synthetic networks to systematically compare and determine the importance of certain pore network features on the diffusivity of gasses, using a KMC methodology to simulate fluid flow. Examples of pore network models that can be used to realistically represent pore networks include those proposed by Al-Kharusi and Blunt (2007), Jiang et al. (2012) and Raoof and Hassanizadeh (2010), among others. The 3D KMC algorithm proposed in this work can be used to simulate the transport of fluids within any of these pore network models.
Pore Connectivity According to Sahimi (2011), the pores in natural porous media can be divided into two groups: the pore bodies, which have the largest contribution to the porosity, and the pore throats, which are narrow channels connecting pore bodies. To construct our pore networks, we consider both pore bodies and throats. The connectivity level determines the number of connections between neighboring pores and therefore the number of pore throats in the network. The minimum number of throats per pore is 0, when the pore is not connected to any adjacent pores; the maximum is 6, when the pore is surrounded by neighboring pores in all 3 dimensions and it is connected with them all. The maximum number of connections is defined by the geometry of the voxels implemented for the construction of the 3D KMC lattice. Since cubic voxels are used in this 3D KMC model, the maximum number of connections is 6. If a different mesh type was used to generate the 3D KMC lattice (using space-filling polyhedra for the geometry of the voxels), the number of maximum connections would be different.
To render the generated pore networks realistic, we used discretized log-normal distributions to determine the degree of connectivity (see details in Supplementary Material). For each distribution, we specified the values of parameters μ c and σ c . We performed sensitivity analysis on both parameters separately.
The μ c values considered were 2, 3, 4, 5 and 6. In these networks, the σ c value was kept constant at 1 (Fig. 1a). Because connectivity is one of the most important parameters affecting fluid transport, we ensured that our simulations were not biased by a poor choice of the connectivity via the μ c parameter. To this end, we generated two sets of networks, five networks with low connectivity (μ c = 2) and σ c ranging from 0.5 to 2.5 (shown in Fig. 1b) and five networks with high connectivity (μ c = 5) with σ c values also from 0.5 to 2.5 (Fig. 1c). We performed sensitivity analysis on these two sets separately.
The σ c parameter values considered were 0.5, 1, 1.5, 2 and 2.5. In these cases, the μ c value is kept constant and equal to 2 and 5, for the low-and high-connectivity networks, respectively. As σ c increases, the heterogeneity of the system increases because the connectivity distribution becomes broader, and the pores with low connectivity coexist with pores that are highly connected and can easily transport fluids.
Porosity To quantify the effect of porosity on fluid diffusion, we consider pore networks of increasing porosity: 5%, 10%, 15%, 20%, 25%, 30%, 35%, 40%, 45% and 50%. These porosity values, when multiplied by the total number of voxels in the 3D KMC lattice, yield the number of pores present in the network. The pores are randomly inserted in the pore network in an unbiased way (no spatial anisotropy). To ensure that the diffusivity calculated does not depend on the pore distribution, we generate 5 configurations for each porosity. We calculate the diffusion coefficient in each network and report the average, together with the uncertainty quantified using the standard error formula: n] , and n is the number of 3D pore networks considered for the calculation. In Fig. 2, we present an example of the pore networks with porosity 10%, 20%, 30%, 40% and 50%. In Fig. 3, we present an example of three different pore networks with 25% porosity. For all networks presented in Figs. 2 and 3, the connectivity of the networks is high (μ c = 5).
Pore Throat Size Distribution We previously found, for 2D pore networks, that the pore size distribution strongly affects permeability (Apostolopoulou et al. 2019). We perform here a sensitivity analysis on the effect of pore throat size distribution. To be consistent with the literature, we consider log-normal distributions, which characterize shale rock formations (Yang and Aplin 2007). We examine the impact of the distribution's μ t and σ t , within 10 systems. In the first five pore networks, we keep σ t constant and equal to 1, while we change μ t (5, 10, 15, 20 and 25 nm) (see Fig. 4a). These μ t values are relatively low, as shales are characterized by a significant amount of microporosity (Hu et al. 2015;Nelson 2009). In the other five pore networks, we kept μ t constant and equal to 25 nm, while we (1) Error = STD(D(i)) √ n changed σ t (0.5, 1, 1.5, 2 and 2.50 nm). As σ t increases, the distribution becomes broader and the level of heterogeneity in the network increases (see Fig. 4b).
Pore Throat Chemistry We considered systems of varying mineralogy content, namely, silicon oxide (silica), magnesium oxide (MgO), aluminum oxide (alumina), calcite and muscovite. These minerals represent those found in the inorganic part of shale rocks  (Merriman et al. 2003;Zou et al. 2013;Backeberg et al. 2017). Five different compositions, shown in Table 1, were studied. The composition of the pore throats was used, together with their size, to determine the diffusivity of supercritical methane. To assign diffusion coefficient values to the pore throat diameters for each substrate, we used the digital libraries produced previously (see Fig. 5) (Apostolopoulou et al. 2019). Note that to generate these libraries, we considered only supercritical methane, as a pure fluid. If moisture was Fig. 2 Representation of the 3D and 2D pore networks generated while considering porosity values 10-50%. As porosity increases, the pathway a particle can travel, referred to as cluster connectivity, increases  Pore throat size distributions considered when a σ t = 1 nm and remains constant, while μ t changes, and b μ t = 25 nm and remains constant, while σ t changes present, water would fill (partially or fully) some of the pore throats, generating additional barriers to fluid transport (Le et al. 2017). For the pore networks considered in Sects. 2.1.1, 2.1.2 and 2.1.3, the composition of pore throats was 100% silica.

Pore Network Generation Algorithmic Steps
We generated several sets of pore networks, within each of which only one pore feature was changed, while the other parameters remained constant. Table 2 summarizes the systems investigated per pore network connectivity category. In our analysis, we considered three types of pore network connectivity, Type G, which is the general case (see Fig. 1a), Type L and Type H, consisting of poorly and highly connected pores, respectively. The number of voxels in the simulation box was kept constant for the various pore networks. The simulation box is made by inserting 10 voxels in the x and y directions, respectively, and 5 voxels in the z direction. Periodic boundary conditions were applied in all three dimensions. The outcome is a 10 × 10 × 5 matrix, which is presented in Fig. 6.    Table 1 Type H μ c = 5, σ c = 1 Figure 4 5-50 Table 1 To construct the pore networks, we followed the algorithmic steps presented in Fig. 7. The first parameter selected is the matrix porosity, which defines the number of pores inserted in the lattice. To speed up the 3D KMC simulations, in our models all the pores have the same size. Each of these pores can be considered as a well-mixed system, where molecular displacements happen fast in the timescale sampled by the stochastic KMC approach. We assume that pore-to-pore diffusion events via throats are difficult to occur and can be considered as rare events. The positions of the pores are randomly assigned within the lattice, via the uniform random number generator. If a lattice voxel chosen in the routine is already occupied by a pore, a different position is selected for the new pore.
Once all pores are placed in the lattice, the connectivity is required to identify the pore throats. While sampling through a connectivity distribution, a list of connections is generated. To test whether the generated distribution of connectivity matches the target distribution, we analyze μ c and σ c of the two distributions. If the % difference (positive or negative) between these values is less than 5%, we proceed; otherwise, we generate a new distribution. We present the comparison between the two distributions obtained (target and stochastically generated) for some of the networks generated, in Supplementary Material. The number of elements in the produced list of coordination numbers, in other words the number of samples taken from the connectivity distribution, is equal to the number of pores in the network, as defined by the porosity value. At this stage, each pore has a maximum of 6 connections available. Once the list of connections has been generated, each pore is matched randomly with a coordination number. Depending on the coordination number selected, there are up to 6 possible KMC rates that can be determined for a given pore.
These rates define the direction a particle may follow traveling from that pore, toward the voxel on top, bottom, left, right, back and forth. The direction of the jump is randomly selected, so that diffusion is equally probable in all three dimensions. Once the direction is defined, we test to see whether a pore is present in the neighboring voxel. If the voxel does not contain a pore, a different direction is selected, etc. All connections must begin and end in a pore body; otherwise, a new configuration of pore bodies is produced. The process is repeated until all connections are assigned to voxels containing pores and the appropriate number of neighbors surrounding them. However, using this methodology, a pore may be connected to a neighboring pore with zero connectivity. In that case, if a particle ends up in this pore, it remains "stuck" and is unable to transport anywhere else during the simulation. The process of assigning coordination numbers to pores is time-consuming but necessary, especially when generating networks with low porosity and high connectivity.
Next, the width of the pore throats is determined. This parameter determines the kinetic constants for the pore-to-pore transitions. The pore throats are considered to be slit-shaped, and their width is dictated by a selected PSD. Using a stochastic approach, we sample through the target PSD and select the width of each pore throat. To validate the network, we plot the generated PSDs against the reference PSDs presented in Fig. 4. We extract and compare the μ t and σ t values of the generated PSDs. We present the comparison between the two PSDs, for some of the networks generated, in Supplementary Material. If the % difference between the produced and actual PSDs is less than 5%, we proceed to assign pore chemistry, which is an input parameter. According to the % of silica, MgO, alumina, calcite and muscovite in the network, the number of pores made of these materials is specified. The last step of the algorithm is to assign the KMC rates across the network, for which we employ the digital libraries of Fig. 5. Depending on the pore throat chemistry selected, we find the diffusion coefficient associated with the pore diameter selected from step 5 from the corresponding digital library. The diffusion coefficient values are then transformed to KMC rates using (Jansen 2008;Flamm et al. 2009): When all 7 algorithmic steps are completed, we use the generated list of KMC rates to inform our 3D KMC model.

3D KMC Model for Fluid Diffusion
Our KMC approach, applied to 1D, 2D and 3D pore networks, is described by Apostolopoulou et al. (2017Apostolopoulou et al. ( , 2019a. The underlying model of the KMC simulation is the Master equation (Eq. 3), which can be thought of as a "probability balance" (Darby et al. 2016). The Master equation expresses the rate of change for the probability P p (t) of finding the system in state p, at a time t, in terms of the probability influx from other states q, and the probability efflux toward other states q (Van Kampen 2007): The state vectors p and q in Eq. (3) capture the information necessary to describe the location of diffusing fluid particles in the porous network of interest. A state vector stores the number of particles contained in every voxel of the network. This state vector is updated over time. W pq and W qp are the propensities of the p-to-q and q-to-p transitions, respectively, and are calculated by multiplying the correspondent KMC rate constants by the number of particles contained in the correspondent voxels. Generic Master Eq. (3) can be used to describe the diffusion of a particle in a 3D system from voxel A (x, y, z) to voxel B (x, y + 1, z) as follows: in state q, voxel A has n A + 1 particles and voxel B has n B − 1 particles. The probability per unit time (propensity) for the aforementioned diffusion event to happen is given by the KMC rate for the A-to-B transition multiplied by the number of particles in the A voxel, n A + 1. If the transition is accepted and therefore performed, the population in the A voxel will be n A , while the number of particles in the B will be n B , leading to state q.
Conducting a KMC simulation requires the frequent generation of random numbers for the selection of the event at each step and for the calculation of the time required for the transition to happen. To this end, the Mersenne Twister MT19937 uniform random number generator was used (Matsumoto and Nishimura 2000).
To simulate diffusion, we use a single particle randomly inserted in the lattice. To obtain accurate statistics, several independent runs are performed, as summarized in Fig. 8. Because the initial particle position is randomized 75 times, 75 initial configurations are generated. Using each initial configuration, 10 independent KMC runs are performed to compute the average diffusion coefficient.
We monitor the particle trajectory and calculate the mean square displacement (MSD) over time, which, used in Einstein's equation (Eq. (4)), yields the fluid diffusivity.
In Eq. (4), are the MSDs in the XYZ space. The simulation time is 70 ns, and a sample to monitor the particle's position is taken every 0.07 ns. From the 10 independent runs performed, we calculate the standard error using Eq.
(1). It should be noted that the 3D KMC model implemented here was validated previously against results obtained for a series of 3D slit pores (Apostolopoulou et al. 2019).
Our bespoke 3D KMC model was implemented here using the algorithm of Fig. 7 to assign the input KMC rates to the various lattices.
(4) Fig. 8 Algorithm for 3D KMC independent runs. Note that a total of 750 simulations are conducted to obtain the diffusion coefficient for one network

Effect of Network's Connectivity
To quantify the impact of pore throat connectivity, we generated 15 networks with constant porosity (15%), chemical composition (silica) and pore throat size distribution (μ t = 25 nm, and σ t = 1 nm). The only parameter changing among these 15 networks was the connectivity level, referred to as coordination number (see Fig. 1). In the first 5 networks (Type G in Table 2), μ c was increased from 2 to 6, while σ c was kept constant at 1. For the next 5 networks (Type L), connectivity was low (μ c = 2) and σ c was gradually increased from 0.5 to 2.5. For the last 5 networks (Type H), connectivity was high (μ c = 5) and σ c ranged from 0.5 to 2.5. Note that the pore throat widths are described via a size distribution, rather than a single value. Thus, when a pore is connected to other pores, the slit pore throats used as channels have different widths. The smaller the widths, the lower the likelihood of transition of particles, due to the higher resistance to flow. Thus, when the particles are given the option to choose between wide and narrow pore throat widths, it is expected that the fast transitions through the wider pore throats will be more frequent compared to those through narrower throats. Diffusion was estimated via our KMC method as described in Eq. (4). In Fig. 9a, we present the results obtained from the first 5 networks (increasing connectivity). We observe an exponential increase in diffusivity as the connectivity of the networks increases. Indeed, when μ c increases from 2 to 6, the diffusivity increases by almost half an order of magnitude. This is expected for two main reasons. Firstly, as the connectivity increases, more pore-to-pore connections are available, yielding longer pathways. This is captured by the MSD obtained and confirmed by analysis of simulation trajectories (see Supplementary  Material). Secondly, since more pathways are generated, the particles have more choices among pathways, which enables them to choose the "faster routes" over the obstructed ones.
To compare the 15 generated networks presented in Fig. 1, we assume an exponential relationship (see Fig. 9) between the pore network's connectivity and methane diffusivity that can be expressed as: In Eq. (5), y is the diffusion coefficient (in m 2 /s), x is the value of the distribution's parameter (connectivity), K and J are constants, and i = 1, 2 and 3 for the panels A, B and C (5) y = K i e J i x Fig. 9 Effect of the connectivity of the networks on methane diffusivity. a Results obtained when the σ c value of the networks equals 1 and μ c increases. b and c Results for networks Type L and H, respectively, with increasing σ c . The distributions representing the networks in a-c are shown in Fig. 4a-c, respectively presented in Fig. 9, respectively. According to the R 2 ) value, calculated during fitting, the variance in the diffusivity is strongly associated with the variance in the network's connectivity characteristics. Our goal is to assess how the connectivity distribution parameters impact the value of the constant J. The error bars calculated are shown in red. They are relatively small for all networks considered. Visualizing the results obtained, together with the error bars, no overlaps are observed, confirming that the number of independent runs selected for our calculations was appropriate to yield statistically significant diffusivity values.
In the remaining two sets of pore networks (Fig. 9b and c), as the connectivity distributions become broader, pores that are poorly connected co-exist with others that are highly connected. It is not easy to predict how these features affect fluid transport. On the one hand, the increase in the proportion of highly connected pores could promote fast diffusivity. On the other hand, more pores with low connectivity could yield the opposite effect. It is also likely that the results depend on the sensitivity of fluid transport on the existence of poorly vs. highly connected pore bodies.
For both Type L and H systems ( Fig. 1b and c, respectively), diffusivity is found to decrease as σ c increases ( Fig. 9b and c). Considering that the diffusion coefficient reduces when the amount of pores with low connectivity is increased, and that the connectivity distributions considered are log-normal, it is concluded that the diffusivity in our systems is more sensitive to the presence of poorly connected pores. Comparing J 2 and J 3 values, calculated from the exponential trendline fitting (Eq. (5)), it is found that J 2 = -0.07 and J 3 = -0.13, suggesting that when the connectivity is high, the decrease in diffusivity due to increasing σ c is more pronounced. To better understand the reason behind this observation, we investigated the distributions presented in Fig. 1b and c. When the connectivity is low (μ c = 2), the number of highly connected pores increases as σ c increases. On the contrary, when the connectivity is high, an increase in the σ c value results in the appearance of more poorly connected pores. This confirms that low-connectivity pores have a more significant impact on diffusion, compared to high-connectivity ones. In practical terms, this observation suggests that when the aim is to increase diffusivity and permeability for a network characterized by high connectivity, applying techniques to increase the network's connectivity (i.e., generate secondary fractures in a reservoir) will result in modest improvements of the transport properties. On the contrary, the same treatment will have much stronger beneficial effects when applied to a low-connectivity network.
Moreover, the coefficients J 1 (equal to 0.32), J 2 and J 3 are indicative of how sensitive the diffusivity is when the network's connectivity changes. Since J 1 is significantly higher than the absolute values of J 2 and J 3 , it can be deduced that the diffusivity is more sensitive to changes in the network's connectivity (μ c ), rather than the degree of heterogeneity, captured by the increase in the value of σ c .

Effect of Porosity
To investigate the effect of porosity, we consider two sets of pore networks, Type L and H. In each set, we consider 10 porosities, ranging from 5 to 50%. To obtain meaningful statistics, for each porosity we generate 5 equivalent networks, following the protocol described in Fig. 7. In all cases, the pore throats that connect the pores are made of silica and are characterized by pore width distributions described by μ t = 25 nm and σ t = 1 nm.
In Fig. 10, we report the mean values obtained for the diffusion coefficient, as a function of porosity, for the two network families. In panel A we present the results obtained when considering networks with low connectivity, in panel B those with highly connected pore bodies. The results show that, for porosity ranging from 5 to 30%, there is no significant change in the diffusion coefficient calculated. However, a sharp, almost exponential increase is predicted when the porosity exceeds the 35% threshold. This trend is observed in pore networks with both low and high connectivity. The results can be explained as follows.
As the porosity of the system increases, more pore bodies are present within the system. Considering that the connectivity remains unaltered, the new additional pores provide the connections with their adjacent neighbors as described by the connectivity distribution curves. As the number of pores in the system increases, more pores become interconnected (see the trajectory analysis presented in Supplementary Material), generating longer pathways. This is also shown in the 3D pore networks presented in Fig. 2. As this higher effective connectivity due to increased porosity appears, the initial plateau in calculated diffusivity transitions to a rapidly increasing function of porosity, as shown in Fig. 10. Since our 3D lattices are periodic in all 3 dimensions, we do not expect the 35% threshold to be strongly dependent on the representative elementary volume (REV) used to construct our systems. In fact, a similar threshold porosity value (~ 39%) was observed by Prasianakis et al. (2018), who simulated an acidic aqueous mixture diffusing through a calcite mineral matrix, dissolving part of it, thus increasing the system's porosity. An exponential increase in permeability was observed when the 39% porosity threshold was reached. In agreement with those observations, our results suggest that the network's transport behavior can be significantly improved when the porosity is artificially altered. However, the process of modifying the network's porosity utilizing chemical reactions should be carefully planned, as the reaction by-products (mineral precipitation in the case study by Prasianakis et al. 2018) could clog the existing pore bodies, thereby reducing network connectivity, and thus fluid transport.

Effect of Pore Throats Size Distribution
We considered 20 pore networks, in which the pore throat size distribution followed lognormal distributions. In the first 10 pore network sets, we maintained the distributions' σ t equal to 1, varying μ t . Half of the networks generated were Type L (low connectivityμ c ), Fig. 10 Relationship between the porosity of the networks and methane diffusivity for (a) Type L and (b) Type H networks (see Table 2). The standard errors are calculated considering 5 equivalent pore networks for each porosity value and applying Eq. (1) and half were Type H (high connectivityμ c ). In Fig. 11a we present the results obtained for the low-connectivity systems and in Fig. 11b those for the highly connected networks, respectively.
For both panels A and B in Fig. 11, the results show that as the pore throat size increases, the methane diffusion coefficient also increases. This is expected, since wider pore throats allow fluid to transport faster, as shown in Fig. 5. The results in Fig. 11 show that if the pore throats are initially narrow, any further size increase yields moderate increases in diffusivity, while the effect is less pronounced if the pore throats are initially wide (i.e., see the plateau in the datasets). This is a direct consequence of the relationship between diffusion coefficient and pore width observed for single pores (see Fig. 5). The practical implication of this observation is that once the pore throats have reached sizes of ~ 30 nm, widening them further has little effect on the overall diffusion coefficient of supercritical methane.
In Fig. 11, we also report the error bars as estimated by applying Eq. (1). In panel A the error bars are smaller than the symbols, while in panel B they are slightly larger. Analysis of the errors shows no overlaps among results obtained for different pore throat size distributions, which confirms that the results presented are statistically significant. Comparing the results shown in panel A to those in panel B, an almost 2.5 times increase is observed, which is due to the change in the connectivity of the networks. By observing the increasing trend of the diffusivity as a function of the pore throat width, it is observed that the connectivity of the network (low or high) affects the diffusivity in a similar manner, from a qualitative perspective.
The next 10 pore network sets also consisted of 5 Type L and 5 Type H (see μ c in Table 2). The μ t of the PSDs was 25 nm, and σ t ranged from 0.5 up to 2.5 nm. For consistency, the pore throats were all made of silica. The effect of changing σ t in the pore throat size distribution is quantified in Fig. 12. It is expected that as the σ t increases, narrow pore throats coexist with wider ones. We consider pore networks with low and high connectivity (results shown in panel A and B of Fig. 12, respectively). It is observed that the heterogeneity of the pore throat size distributions, i.e., the σ t value, has no impact on the diffusion coefficient of methane. This result is phenomenologically explained by the following observation: when particles have the option of diffusing to neighboring pore bodies traveling across narrow or wide pore throats, they most likely follow the pathways of lowest Fig. 11 Relationship between the pore throats width size distribution and methane diffusivity for (a) Type L and (b) Type H networks (see Table 2). For the μ t of the pore throat widths, we consider the distributions presented in Fig. 4a resistance, i.e., the wider pore throats. The same trend is observed for both low-and highconnectivity pores. However, the self-diffusion coefficients obtained in networks with high connectivity are larger than the others by a factor of almost 2.5. When the networks are highly connected, there are more possible pathways available for the particles to diffuse, and the proportion of the wider pore throats is considerably increased.

Effect of Pore Throat Chemistry
To investigate the effect of the chemistry of the pore throat widths, we generated 2 sets of 5 networks each. The first set contained networks with low connectivity, the second, high connectivity. The pore throats size distribution and the network porosity were kept constant. The composition of the pore networks in each of the 2 sets of 5 networks is shown in Table 1. The results obtained for the self-diffusion of methane, together with the error bars, are presented in Fig. 13. The chemical composition of the pore throats seems to have little effect on the diffusivity calculated. A qualitatively similar observation was reported by lsmail and Zoback, whose permeability measurements indicated that shale mineralogy does not have a strong effect on permeability (Al Ismail and Zoback 2016). Based on our results, this observation could be explained considering the diffusion coefficients Fig. 12 Relationship between the pore throat width size distribution and methane diffusivity for (a) Type L and (b) Type H networks. For the σ t of the pore throat widths, we consider the distributions presented in Fig. 4b   Fig. 13 The effect of pore throat chemistry on the self-diffusion coefficient of supercritical methane. Results for pore networks with low (Type L) and high (Type H) connectivity are shown in panels A and B, respectively obtained in single pores of different widths (see Fig. 5): when the pores are of width ~ 3 nm or wider, the chemistry of the solid substrate has little effect on the diffusion coefficient of supercritical methane. This result could change when the composition of the fluid changes. For example, in hydrated pore networks, water particles could accumulate near the pore throats, depending on the chemistry of the solid, yielding additional kinetic barriers to gas transport that are not explicitly considered in this study (Bui et al. 2017).

Discussion and Sensitivity Analysis
We present here an overview of the results obtained above and perform a sensitivity analysis to determine the impact of connectivity, porosity and pore throat size distribution on the diffusion coefficient of supercritical methane within the various pore networks. To determine the contribution of each network feature on the response variable, i.e., the diffusion coefficient of methane, we plot the % change in this variable as a function of the % change in the network feature studied. An example of how this % change values are calculated is provided as Supplementary Material. The base case for the calculation of the % change in connectivity μ c is μ c = 2, for pore throat connectivity μ t is μ t = 5, and for porosity is 5% (or 0.05). In this analysis we do not consider the effect of the spread for the connectivity and pore throat width sizes, since it was found to have moderate to low significance on methane diffusivity. For the same reason, the effect of the pore throat chemistry is not considered either. Because the analysis presented in Sects. 3.1-3.4 quantifies how changes in pore network features promote or hinder diffusivity, we consider here only the % change (positive or negative) in output. In Fig. 14, we present the results for pore networks with low (panel A) and high (panel B) connectivity.
Comparing the results in panels (A) and (B), it is observed that both the trends and their magnitudes are consistent. It is found that: (1) changes in the pore network connectivity have the greatest impact on methane diffusion; (2) changes in pore throat size distributions yield the second largest impact on methane diffusion, at least when the network feature changes up to 400%; (3) changes in porosity of up to 400% have little effect on methane diffusion. For our analysis, we could consider more PSDs to ensure that the ranking of the porosity and PSD, as presented in Fig. 14, would not change, when the change in variable is greater than 400%. However, we can support this claim by considering the digital libraries presented in Fig. 5. Considering those results (methane diffusion coefficient in isolated pores as a function of pore width and pore chemistry), supercritical methane shows a selfdiffusion coefficient within the pores that is bulk-like when the pore width is of at least 30 nm. To achieve a 700% change in the pore throat size distribution's μ t , μ t should be 40 nm. For this distribution the diffusion coefficient of supercritical methane is expected to be slightly higher than 1.18 × 10 −8 m 2 /s, which corresponds to an almost 40% improvement in diffusivity. However, for 700% change in porosity, the change in diffusivity is almost 45%. Moreover, for higher porosity values the diffusivity increases exponentially, whereas for larger PSD values, the diffusion coefficient should remain constant, as the bulk-like behavior has been reached. Thus, it can be assumed that the porosity and PSD ranking will remain the same.
From a practical perspective, the observations derived from Fig. 14 provide recommendations on how to increase the transport properties within a pore network, when the network's characteristics are known. Since connectivity is the most influential parameter, any strategy aiming to increase it is expected to result in significant improvements in fluid transport. For example, during the hydraulic fracturing process, the creation of more dense secondary fractures could yield significant increases in production rates. This observation is consistent with field measurements during the development of the Wolfcamp and Spraberry fields (Jaripatke et al. 2018). Improvements in the network's porosity are also expected to be beneficial, but only after the 'critical' threshold of ~ 30% porosity is reached. This is consistent with results reported by Prasianakis et al. (2018). The size distribution of the pore throats is an important feature, but, in the case of unconventional reservoirs, it depends on factors that operators have no control over. On the other hand, this feature can be altered in engineering materials used for example in catalysis (e.g., zeolites). Finally, it is important to highlight that the synergistic alteration of both connectivity and porosity is expected to be the most impactful, in qualitative agreement with the results reported by Prasianakis et al. (2018), where chemically active agents dissolved the rock matrix, resulting in a simultaneous increase in porosity and pore network connectivity.
While the factors considered in this study have a significant impact on the diffusivity of supercritical methane, additional parameters will affect permeability. For example, tortuosity can impact the diffusivity of fluids through porous media and therefore permeability. Tortuosity in anisotropic porous media is found to be correlated with porosity and connectivity. In fact, a vast literature explores analytical and empirical relationships between these features (Guo 2012;Shen and Chen 2007). Because homogenous systems were considered in the present study, the impact of tortuosity on diffusivity was considered out of scope. Perhaps future studies could explore how tortuosity, connectivity and porosity affect fluid transport in 3D anisotropic systems.

Comparison Between 2D and 3D Predictions
It is often debated whether using sections of 3D networks (i.e., 2D representations) can yield accurate estimates for a material's permeability. The appeal of using 2D networks over 3D representations is due to the lower computational requirements. For instance, Beckingham et al. (2013) suggested that pore network models informed by analyses of 2D and 3D images at comparable resolutions can produce permeability estimates with relatively good agreement (Beckingham et al. 2013). On the other hand, Begg and King calculated the effective permeability of shale samples and argued that 2D sections can provide valid insights on a limited number of cases, and thus, 3D models should be used instead (Begg and King 1985). Similarly, Trinh et al. (2002) reported that when the percolation threshold of 2D and 3D media is substantially different, the diffusivities calculated are expected to vary significantly even in isotropic media with similar porosities and pore configurations. In our pore networks, the pore bodies and pore throats are randomly assigned to the network, yielding isotropic systems. To investigate whether the 2D pore networks can accurately reproduce the results obtained from the 3D pore networks, we obtained 2 H-type pore networks and sliced them along the z-axis to produce 10 2D pore networks. The 3D pore networks selected have 5% and 45% porosity. In Fig. 15 we present the 5 2D slices obtained for the 5% pore network (panel A) and the 5 slices for the 45% pore network (B).
To calculate the diffusion coefficient, we follow the procedure used in our previous calculations. In Fig. 16, we present the results obtained for the 2 pore networks. We report the individual diffusivities calculated per slice, the average diffusivity obtained from the 5 2D slices and the 3D diffusivity, which is presented in Fig. 10, including the error bars. Our results suggest that although the generated pore networks are isotropic, the diffusivity  Diffusion coefficients obtained when considering 2D and 3D models for pore networks with (a) 5% and (b) 45% porosity calculated from the 2D slices is not equivalent to the one obtained when considering the 3D model. It is also observed that the two predictions deviate the most for the 45% porosity networks. This could be due to the generation of connected clusters in these highly porous networks that result in the appearance of longer connected pathways. This effect is significantly amplified in 3D networks compared to 2D networks resulting in a greater degree of deviation between model predictions. On the contrary, in the 5% porosity network, most of the generated pathways connect 2 or 3 consecutive pores; hence, the addition of the 3 rd dimension has a more limited effect of the diffusivity predictions. Our observations are in agreement with the results reported by Trinh et al. and show that the generation of connected clusters-in other words, the effect of the percolation threshold-can lead to significantly different diffusivity predictions when comparing isotropic 2D and 3D networks with equivalent structural characteristics. If conducting simulations using 3D pore networks is prohibitive due to the features of a given dataset (e.g., inadequate 3D pore resolution), reconstruction of the 3D network from high-resolution 2D thin slices is recommended as an alternative approach, following the methodology reported by Okabe and Blunt (2004).

Conclusions
In this work, we discuss a bespoke 3D stochastic algorithm, which implements the kinetic Monte Carlo (KMC) formalism to predict the diffusivity of fluid molecules through pore network models. In particular, we use the bespoke KMC approach to quantify the impact of four pore network characteristics on the transport of fluids through porous materials. The characteristics considered are the pore network's connectivity, porosity, pore throat widths and pore chemistry. The 3D pore networks simulated consisted of pore bodies, pockets where particles accumulate, connected via slit-shaped pore throats. We have followed a rigorous protocol when generating these synthetic pore networks for our simulations, to ensure that the results reported are both free of bias and statistically significant.
For all the case studies considered, we generated 10 equivalent 3D pore networks and performed 75 independent 3D KMC simulations for each pore network. We report as results the average diffusion coefficient calculated per case study and the corresponding error bars, using a standard error formula. We found that the results obtained together with the reported errors are statistically significant for most of the case studies considered. This confirms that the number of independent runs selected for our stochastic calculations was appropriate to provide results that are statistically valid. The systems that yield overlapping results are those investigating the effect of pore chemistry and the pore throats distributions' variance in the diffusion coefficient of supercritical methane. This is because of the similar transport properties of the substrates considered and the ability of particles to select a low resistance pathway when the networks become highly heterogeneous.
We found the pore network connectivity to have a significant impact on the diffusivity of supercritical methane. This is expected, since the increase in connectivity results in the generation of additional pathways for the particles to diffuse. In fact, as the connectivity increases, a sharp, almost exponential, increase in the diffusivity was observed. This means that if a low-connectivity network was to be treated to improve its connectivity, the improvement in transport properties is expected to be significant. Porosity was also found to yield an exponential increase in methane diffusivity, after a threshold value (~ 35%) is reached. Similar effects have been reported in other computational studies, further substantiating our observations. For low porosity values, the diffusion coefficient exhibited a slow increase with the increase of porosity. We hypothesize that this is because, for low porosity values, the pores may be connected, but the level of cluster connectivity (path length) is low. When considering pore throats with increasing width diameters, the calculated diffusion coefficient increased with a moderate rate, until a plateau was reached. This is probably due to the digital libraries we considered to correlate the pore widths selected and the diffusion coefficient, when considering single pore throats. For the case of two-phase fluids with inherent thermodynamic barriers, it is expected that the corresponding digital libraries would exhibit a different-substrate-sensitive correlation between throat width and diffusivity.
When comparing the diffusivity obtained from 2D to 3D models, we obtained results that were statistically different. In particular, we observed that as porosity increases, the deviation between the prediction of the 2D and 3D models increased significantly.
Our sensitivity analysis revealed the connectivity to be the most important parameter that affects the diffusivity of rarefied fluids (such as supercritical methane) followed by the pore throat widths and finally porosity. However, when considering fluids in different phases, i.e., liquids, this trend may change. The network connectivity and porosity are both characteristics that can be artificially altered. For example, during the hydraulic fracturing process, the creation of denser secondary fractures and the treatment of the rock formation with chemically reactive compounds can increase the network's connectivity and porosity, respectively. We believe that a synergistic approach that combines the improvement of both features would result in the best outcomes, as long as these strategies are carefully planned and sustainable.