Effect of Chain Length and Pore Accessibility on Alkane Adsorption in

Article pubs.acs.org/EF

Effect of Chain Length and Pore Accessibility on Alkane Adsorption in Kerogen Kerstin Falk,† Roland Pellenq,†,‡,§ Franz Josef Ulm,† and Benoît Coasne*,†,‡ †

Department of Civil and Environmental Engineering ‡MultiScale Material Science for Energy and Environment UMI:3466 CNRS-MIT, Massachusetts Institute of Technology, Cambridge 02139 Massachusetts United States § CINaM, CNRS/Aix Marseille Université, Campus de Luminy, 13288 Marseille Cedex 09, France S Supporting Information *

ABSTRACT: Configurational-biased Grand Canonical Monte Carlo and Molecular Dynamics simulations were used to investigate the adsorption and structure of n-alkanes in realistic models of kerogen (the organic phase in gas shales). Both the effects of the n-alkane length (from methane to dodecane) and of the adsorbent porosity/density and connectivity were considered. For all n-alkanes, due to the subnano pore size of kerogen, adsorption follows a Langmuir-type adsorption isotherm as the adsorbed amount increases continuously with pressure until maximum loading is reached. While all n-alkanes get adsorbed very efficiently in low density kerogen, size exclusion in high density kerogen lead to low adsorbed amounts for molecules larger than hexane. Using Molecular Dynamics simulations reflecting the setup of actual adsorption experiments, we also probed the effect of pore accessibility on n-alkane adsorption in kerogen. All pores are connected to the surface and, hence, accessible for the low density kerogen. In contrast, most pores are isolated and inaccessible in the kerogen with the high density. We also address the effect of mesoporosity and its connection to nanoporosity on alkane adsorption in kerogen.

1. INTRODUCTION Organic-rich shale formations are heterogeneous multiscale materials made up of tightly packed, fine-grained sedimentary rocks with inclusions of an organic material called kerogen. Due to the very low permeability of this tight rock and a poor understanding of molecular transport in its complex structure, hydrocarbon extraction requires stimulating techniques such as hydraulic fracturing, which raise both environmental and economical concerns.1−3 Researchers are therefore urged to (1) provide tools to assess available resources from data that is amenable with simple experiments and (2) improve our understanding of adsorption and transport in shales in order to optimize hydrocarbon recovery.4,5 Both research efforts require an understanding of the structure, thermodynamics, and dynamics of hydrocarbons trapped in kerogen in organic-rich shale formations at the nanoscopic scale. While kerogen is the key component in shales since it produces hydrocarbons through its decomposition, its physical and chemical properties and the behavior of trapped hydrocarbons remain poorly understood. Kerogen is an intrinsically complex medium as it consists of an amorphous porous material which exhibits significant morphological (pore shape) and topological (pore connectivity) variations.6−8 Moreover, depending on its geographical origin and maturity, kerogen exhibits a broad range of density, chemical composition (both in terms of atomic contents and chemical functions), tortuosity, porosity, etc.9 This paper aims at unraveling the adsorption behavior of hydrocarbons in the nanoporosity of kerogen which consists of irregularly shaped pores of nanometer and subnanometer size, by means of molecular simulations. Recently, Collell et al. performed a molecular simulations study of methane/ethane mixture adsorption in artificially created spherical pores inside very dense kerogen matrices. © XXXX American Chemical Society

Under high pressure (1−20 MPa) they observed confinement effects for subnanometer pore diameters, leading to deviations from the classical adsorption models, i.e., extended Langmuir and ideal adsorbed solution models.10 For larger adsorbed molecules, we expect the influence of the confinement on the thermodynamic and dynamic properties of the adsorbed phase to be even stronger. We investigate both the effects of alkane chain length and pore accessibility by considering various n-alkanes (from methane to dodecane) and kerogen models with different densities (and therefore different pore accessibility). Molecular simulations allow us to study the adsorption process in kerogen directly at the relevant nanoscopic scale. Relevant properties, like density, structure, interaction energies, etc., are directly observable in molecular simulations. For alkanes, we use a well established model that has been parametrized on a wide set of experimentally verified properties, like density, structural data, phase behavior and transport properties.11 Figure 1 shows the nanoporous carbon structures used as kerogen proxies for this study: realistic amorphous carbon structures, namely CS1000 and CS1000a, with different densities ranging from 0.7 to 1.5 g/ cm3.12,13 These disordered structures that have been obtained using an atom-scale reconstruction technique can be seen as reasonable molecular models of mature kerogen as they capture its main features (pore size, density, chemical composition including sp2/sp3 hybridization ratio, morphological disorder).14−17 One difficulty in adsorption studies, experimental as well as numerical, can be the long times needed for equilibration. Configurational-biased Grand Canonical Monte Received: September 9, 2015 Revised: October 26, 2015

A

DOI: 10.1021/acs.energyfuels.5b02015 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 1. Typical molecular configuration of nonane adsorbed in the kerogen proxies CS1000 (left) and CS1000a (right) at T = 423 K and P = 25 MPa. The gray sticks represent the C−C bonds of the CS-matrices18 and the cyan spheres represent CHx-groups of nonane (x = 2 or 3). The size of the simulation box is 5 nm with periodic boundary conditions in all directions.

flexible chain by bonds of fixed length 1.54 Å. Bending and torsion of the chain are described by intramolecular interaction potentials for the angles between three neighboring atoms and for the dihedral angles formed by four neighboring atoms. The bending potential is harmonic

Carlo (CBGCMC) is a simulation technique specifically designed to overcome this limitation: with CBGCMC, systems containing long chain molecules can be equilibrated and sampled with relatively small computational effort. In this work, we use CBGCMC to simulate adsorption of nalkanes (methane, propane, hexane, nonane, and dodecane) in different disordered nanoporous carbons representative of mature kerogen phases at different densities. Specifically, we calculated the adsorption isotherms and the isosteric heat of adsorption in CS1000 and CS1000a, as well as in multiscale systems that combine subnanopores and mesopores. As a first step, we consider a constant temperature T = 423 K and a large pressure range (up to 100 MPa) which are relevant to shale gas/oil recovery from kerogen. From a practical point of view, the relatively high temperature is helpful to reach equilibrium and sampling efficiency.

1 C(θ − θ0)2 2

Ubending(θ ) =

(1)

where θ0 is the equilibrium angle and C/kB = 62500 K (kB is the Boltzmann constant).11 The torsion potential is Utorsion(ϕ) = C1(1 + cos(ϕ)) + C2(1 − cos(2ϕ)) + C3(1 + cos(3ϕ))

(2)

with C1/kB = 355.03 K, C2/kB = −68.19 K and C3/kB = 791.32 K.11 Two functional groups of different molecules, or of the same molecule but separated by more than three bonds, interact with each other via the Lennard-Jones (LJ) potential

2. MOLECULAR MODELS AND SIMULATION METHOD Kerogen Models. The kerogen models considered in this work are amorphous nanoporous carbons referred to as CS1000 and CS1000a (Figure 1). Pikunic, Jain et al. reverse engineered these molecular structures from experimental X-ray diffraction data of pyrolized saccharose (CS1000) and pyrolized, and then activated saccharose (CS1000a) by means of Hybrid Reverse Monte Carlo reconstruction techniques; they are thus realistic representations of the real materials.12,13,18 In particular, they give the correct adsorption isotherms for CO2 and CH4.19 While the reconstructed materials are not rigorously kerogens, they share the main features such as pore size distribution and morphological and topological disorders. The CS1000 and CS1000a matrices shown in Figure 1 are identical to the ones in ref 18 which provides a detailed comparison of their structural properties. Most relevant for the purpose of this study is the difference in density and pore size distribution (Figure 4 in ref 18). They vary by about a factor of 2 between CS1000 and CS1000a, representative for, respectively, very high and low density kerogen. Alkane Model. For the confined fluids, we use a well established united atoms force field to model alkanes.11 This is a coarse grained approach where each functional group CHx is represented by a single interaction site (where x = 2 or 3 depending on the position of the group in the alkyl chain). The pseudoatoms of one n-alkane molecule are connected to a

⎛⎛ σij ⎞12 ⎛ σij ⎞6 ⎞ ULJ(rij) = 4εij⎜⎜ ⎟ − ⎜ ⎟ ⎟ ⎝r ⎠⎠ ⎝⎝ r ⎠

(3)

where rij = |ri − rj| is the distance between the two groups i and j. The interaction parameters εij and σij are calculated with the Lorentz−Berthelot mixing rules:20 εij =

εiεj

σij =

1 (σi + σj) 2

(4)

Interactions between the alkane’s CHx groups and the carbon atoms of the kerogen matrix are also modeled with the LJ potential eq 3, with Steele’s parameters for carbon.21 The interaction strength εi and the kinetic diameter σi for the different group types are summarized in Table 1. All LJ Table 1. Lennard-Jones Parameters for Carbon21 and the Different Functional Groups in the Alkane Chains11a atom or group type

ε/kB (K)

σ (Å)

CH4 CH3 CH2 CCS1000(a)

148 98 46 28

3.73 3.75 3.95 3.36

a

Parameters for interactions between different types are calculated from the Lorentz-Berthelot mixing rules eq 4.

B

DOI: 10.1021/acs.energyfuels.5b02015 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Every step, the position of the next group is chosen according to a procedure introduced by Smit et al.24 To increase the probability of finding a low energy state a position j is generated with probability

interactions are truncated at a distance equal half the simulation box size and tail-corrections are added to the energy.22 The kerogen matrix is frozen, meaning that the position of the atoms remain fixed in the course of the simulations. This implies that possible changes of the kerogen structure upon gas adsorption or desorption, like swelling for example, are neglected. However, these effects were found to be small for methane in nanoporous carbons.19,23 Thus, in order to get a first estimate of the adsorption properties of alkanes, it is reasonable to consider a fixed kerogen matrix. Finally, we neglect all interactions of fluid molecules with the hydrogen atoms in CS1000(a). In doing so, keeping in mind that we use a united atom model for the alkanes, fluid and solid are consistently modeled on the same level of coarse graining, i.e., one single interaction site per functional group. Configurational-Biased Grand Canonical Monte Carlo. In adsorption simulations, we are interested in the number of molecules in the system under certain imposed thermodynamic conditions. From the point of view of statistical physics, this implies considering the system in the grand canonical (μVT) ensemble where the volume V, temperature T, and chemical potential μ are kept constant. The number of molecules N in the system can fluctuate and adapt to the imposed temperature and chemical potential. Note that fixing the temperature and chemical potential automatically fixes the external pressure too. The grand canonical ensemble is thus the natural ensemble of any adsorption experiment and simulation. The generic numerical method to study systems in the μVT ensemble is the Grand Canonical Monte Carlo (GCMC) algorithm. GCMC simulates the equilibrium state of a system by producing and sampling configurations that correspond to the given macroscopic state. The probability to find a specific microscopic configuration is thereby given by the Boltzmann statistics in the Grand Canonical Ensemble. The basic principle of a GCMC simulation is to consider random displacements, insertions, and deletions of molecules and either accept or reject each trial move with a probability that follows the Boltzmann distribution mentioned above. For translations and rotations of a randomly chosen molecule, the acceptance rule is acc(o → n) = min[1, exp(−β(U (n) − U (o)))]

pgenerate (j) = c exp( −βu intra(j))

(6)

which depends on the bending and torsion energy uintra added to the molecule with the new CHx group (c is a normalization constant; an example of an algorithm to implement this random generation can be found in ref 22). In addition, not only one but a set of k hypothetical positions is generated according to eq 6. Finally, one position j is chosen from the k-set with a probability pchoose (j) =

exp( −βuLJ(j)) (7)

w

where uLJ(j) is the energy due to the LJ interaction of a CHx group at trial position j with all other atoms (carbon atoms in the matrix and CHx groups in other fluid molecules) in the system. The normalization factor is k

w=

∑ exp(−βuLJ(j)) (8)

j=1

The new group is added at the chosen position and the procedure is repeated. When the molecule has reached the desired chain length l, it has to be decided whether the whole move is accepted or rejected. At this point, the bias that was introduced during the generation of the molecular configuration has to be compensated for by adjusting the acceptance rule (compared to the unbiased one given in eq 5). In order to ensure the detailed balance, i.e., at equilibrium the transition rate from one microstate to another is the same as for the reverse process, the acceptance rule is ⎛ W (n) ⎞ acc(o → n) = min⎜1, ⎟ ⎝ W (o) ⎠

(9)

with W(n) and W(o) the Rosenbluth factors for the new and old configurations, respectively. For a molecule that is partly regrown starting from group s to the last group l (1 < s < l), the Rosenbluth factor is

(5)

with U(o) and U(n) the energy of the system before (“old” state) and after displacement (“new” state), respectively, and β = (kBT) −1. For a more detailed discussion of the tenets of GCMC simulations, we refer to appropriate textbooks.20,22 In the following, we focus on an advanced version called the Configurational-biased Grand Canonical Monte Carlo (CBGCMC) technique.24,25 That version is an extension of GCMC allowing to sample the phase space for complex molecules efficiently. For flexible molecules, we have to sample the configurational (i.e., conformation of the alkane chain) part of the phase space in addition to translation and rotation. However, a completely random choice of a new molecular configuration most probably results in the rejection of the attempted move according to eq 5. As a consequence, the acceptance rate is too low to sample the phase space efficiently. This is why, in CBGCMC, trial configurations are generated according to a biased procedure that favors low energy states. This is achieved by creating configurations with a probability depending on their intramolecular energy as well as the interaction energy with their neighbors. More precisely, in the case of this study, a new n-alkane configuration is ’grown’ by adding one group CHx (x = 2,3) at a time to form the chain.

l

W=

∏ i=s

w(i) k

(10)

where w(i) is given by eq 8 and has to be calculated for each added group i. Insertion of a molecule consists in growing a complete molecule at a random position in the system. We further combine the configurational bias with an insertion bias, meaning that the position for the first group of the new chain is chosen from a randomly produced set of positions according to eq 7. The overall bias is corrected in the insertion and deletion acceptance rules according to ⎛ VβfW (n) ⎞ acc(N → N + 1) = min⎜1, ⎟ ⎝ N+1 ⎠

(11)

where N is the total number of molecules in the simulation box and f is the fugacity of a fictitious reservoir in thermal equilibrium with the system (f is related to the chemical potential μ: f = f 0 exp(βμ)). Note that due to the bias procedure, this reservoir consists of ideal chain molecules, i.e., C

DOI: 10.1021/acs.energyfuels.5b02015 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels molecules only with bonded interactions (eqs 1 and 2). The chemical potential of this ideal chain reservoir is shifted by a temperature-dependent constant with respect to an ideal gas of nonideal chains (including all interactions), which is usually used as the reference state (see refs 22 and 25 for a more detailed discussion). W(n) is the Rosenbluth factor of the newly inserted molecule: W=

w(1) k insert

l

∏ i=2

w(i) kgrow

Table 2. Characteristics of the Two Different Nanoporous Amorphous Carbon Samples, Named CS1000 and CS1000a, That Were Used As Kerogen Models density

(12)

pore width

amorphous carbon

(g/cm )

(nm)

CS1000 CS1000a

1.58 0.77

≈ 0.3 to 0.7 ≈ 0.3 to 1.2

matrix (CS1000a), the number of adsorbed functional groups CHx (x = 2−4) approaches a saturation value that increases slightly with the alkane length. In other words, the pore space can be filled more effectively with dodecane than with methane, despite its larger molecular size. This somewhat counterintuitive observation can be explained as follows. On the one hand, the distance between two bonded functional groups (1.54 Å) is shorter than the molecular radius of methane (σ = 3.73 Å), meaning that two groups within a molecule take up less volume than two single methane molecules. On the other hand, no size exclusion effect is observed in the lower density kerogen, meaning that the same pore volume is available for dodecane as for methane. Overall, this allows a closer packing of CHx groups with longer alkanes than with methane. Indeed, an analysis of the n-alkane tail−tail distance distribution reveals that, in CS1000a, the molecular structure of the adsorbed hydrocarbons is not altered compared to the bulk reference state (Figure 3). By contrast, size exclusion plays a role in the denser matrix (CS1000) where most pores are narrower than 0.5 nm: in CS1000, the maximum number of functional groups is notably lower for longer n-alkanes than for methane and propane. The effect of size exclusion is also evident in the tail− tail length distributions of hexane, nonane and dodecane: not all conformations fit equally well into the highly confined environment. Consequently, some configurations are favored. This leads to deviations of the tail−tail length distributions for long n-alkanes adsorbed in CS1000 with respect to the bulk distributions. It is worth noting that in regular geometries like sufficiently large graphitic slit pores, conformational changes also happen (see Figure 3); this is not due to size exclusion, but because the alkanes adapt to the atomically smooth flat surface.30 Moreover, we find that the isosteric heat of adsorption in the zero pressure limit qst0 increases linearly with the n-alkane carbon number, meaning that every functional group CHx (x = 2,3) contributes equally to the interaction energy of the molecule with the carbon matrix. As a consequence, adsorption of longer n-alkanes is energetically more favorable. This explains why for a given n-alkane length, the maximum adsorbed amount is reached at a pressure that decreases with increasing chain length. This linear dependence of qst0 on the carbon number n has been observed experimentally for various microporous carbons.31 In Figure 4, we show the isosteric heat

(w(i), i ∈ {1,...,l}, given in eq 8). Inversely, deletion of a randomly chosen molecule with Rosenbluth factor W(o) is accepted with the probability: ⎞ ⎛ N acc(N → N − 1) = min⎜1, ⎟ ⎝ VβfW (o) ⎠

3

(13)

To conclude this rapid survey of the CBGCMC technique, we emphasize again that the basic idea of the Monte Carlo simulation methods is to study systems at thermodynamic equilibrium. It is therefore a priori not suited to reproduce time dependent phenomena. The dynamic behavior can be investigated using molecular dynamics (MD) simulations. Indeed, we also performed an extensive MD study of n-alkane transport in the kerogen nanoporosity, which is discussed in detail in ref 26. In the work presented here, MD was used to address the issue of pore accessibility. Furthermore, the pressure-fugacity relation at constant volume and temperature was obtained numerically by performing both bulk phase GCMC and MD simulations, according to the method reported in ref 27. We simulate bulk fluids at 423 K with CBGCMC for different fugacities to obtain the density as a function of fugacity. Independently, the pressure is obtained as a function of density using equilibrium MD simulations of bulk fluids in the NVT ensemble. For all five considered fluids, we checked that the MD simulations reproduce the phase behavior and density-pressure relation given in ref.28 at 423 K and in the pressure range from 0 up to about 30 MPa (see the Supporting Information, Figures S1 and S2). While methane and propane are supercritical at 423 K, hexane, nonane and dodecane undergo a vapor−liquid phase transition at low pressures at such a temperature (see ref 28 for critical data). In the end, the density-pressure and the density-fugacity relations are mapped to get the pressure as a function of fugacity (see the Supporting Information, Figure S3). All MD simulations were performed with the LAMMPS package29 using a Nosé−Hoover thermostat with a relaxation time of 0.1 ps. Time integration was performed every 1 fs.

3. RESULTS Adsorption of n-Alkanes in Kerogen. We use the CBGCMC method described above to simulate adsorption of n-alkanes in the kerogen analogues from Table 2. We compare the adsorption isotherms (423 K) of pure n-alkanes of different lengths. Instead of plotting the number of adsorbed molecules in the porous material at a given pressure, we report in Figure 2 the number NCHx of adsorbed groups CHx. This allows considering a number of CHx groups almost equivalent in size to methane molecules. All samples adsorb n-alkanes at least up to dodecane. For a given pressure, the amount of adsorbed groups depends systematically on both the alkane length and the material density. For high pressures, in the low density

qst = kBT −

⟨UN ⟩ − ⟨U ⟩⟨N ⟩ ∂⟨U ⟩ = kBT − ∂⟨N ⟩ ⟨N 2⟩ − ⟨N ⟩2

(14)

at low loading (i.e., in the limit p → 0). N is the number of adsorbed alkane molecules and U = Uls + Ull + Uintra is the total interaction energy. Angle brackets ⟨⟩ denote statistical averages. We report both the contribution from the alkane-matrix interactions Uls and the contribution from the intramolecular energy of the adsorbed alkane molecules U intra . The contribution arising from the alkane−alkane interactions Ull is D

DOI: 10.1021/acs.energyfuels.5b02015 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 2. Adsorption isotherms at 423 K for methane, propane, hexane, nonane, and dodecane in two different kerogen proxies: CS1000 [right] and CS1000a [left].

chains. First, it is negative, which reflects the fact that deforming an alkane molecule from its equilibrium structure costs energy. Second, its value increases linearly starting from propane. This means that the alkane molecules deform predominantly through bond torsion. Third, the mean intramolecular energy is independent of the pore space. In fact, under the same pressure and temperature conditions, for a given n-alkane, we find the same values for ⟨Uintra⟩ when the molecules are confined in CS1000 and CS1000a, and for the bulk fluid. This is particularly interesting considering the different tail−tail length distributions shown in Figure 3: although the distributions change for high confinement (CS1000), the adsorbed molecule configurations are still energetically equivalent. Accessibility of pore Space in Kerogen of Different Density. GCMC simulation steps do not represent real molecule trajectories that would result from their dynamics. In particular, molecule insertion/deletion can take place anywhere in the system. For adsorption measurements with the GCMC method, this means that every pore can contain fluid molecules, even if it is not connected to the external surface of the sample. Or looking at the problem the other way around, in a real adsorption experiment, these inaccessible pores would stay empty, leading to deviations from the GCMC predictions. Knowledge about the connectivity of the pore space is therefore important for drawing the right conclusions about the adsorption properties. To get a first idea of the connectivity of the pore space in CS1000 and CS1000a, we look at the 3D density distribution of methane at 25 MPa (Figure 5). At this relatively high pressure, we can assume that all pores that are large enough contain methane, since the adsorption isotherm is close to its maximum value. This density destribution gives us some graphical representation of the pore space. The latter was obtained from the CBGCMC simulations, so that the entire pore space (i.e., including closed and open pores) is seen. For CS1000 (high density kerogen), most pores seem to be isolated. In contrast, the pore space looks well connected for CS1000a (low density kerogen). However, it is not possible to make more detailed assessments, since it is unclear if the pore network in CS1000 percolates, or if some pores are inaccessible in CS1000a. In order to make a more quantitative analysis of the fraction of accessible pores, we performed kinetic adsorption simulations with MD simulations. The basic idea is to reproduce the essentials of a real experimental setup. The empty kerogen matrix is put in contact with a fluid reservoir (Figure 6). For practical reasons, we use two reservoirs on two opposite sides of the sample. This allows using periodic

Figure 3. Distribution of the n-alkane tail−tail distance d (see inset): bulk molecules (dashed lines), molecules adsorbed in CS1000 (full lines), and in CS1000a (symbols) and, for comparison, dodecane in a 1 nm wide slit pore with graphene walls (dotted line), at 423 K and 25 MPa. The conformation with minimal intramolecular potential energy has the length L0 = (n − 1) d0 with d0 = 1.29 Å. The propane distribution is displayed up to approximately one-third of the peak height.

Figure 4. Isosteric heat of adsorption qst = kBT − ∂⟨U⟩/∂⟨N⟩ (see also eq 14; kBT = 3.52 kJ/mol) at low loading, i.e., low pressure p → 0, separated into the different contributions from fluid−solid interactions and from intramolecular fluid interactions. The contributions from intermolecular fluid−fluid interactions (not shown) are negligible in the zero adsorbed amount limit. Simulation data is shown for two different kerogen analogues (see Table 2). For comparison, the black filled circles show experimental data for the qst0 of carbon black, taken from ref 31.

not shown as it is negligible at low adsorbed amount. The intramolecular contribution is an interesting property to look at, as this is not directly accessible in experiments. It provides information about the conformation of the confined alkane E

DOI: 10.1021/acs.energyfuels.5b02015 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

6, we compare the number of adsorbed molecules measured with MD (NVT ensemble) and CBGCMC (μVT ensemble). To exclude interface effects, we restrict the density measurement to the inner part (i.e., a slice about 1.5 nm thick) of the kerogen matrices. In CS1000, pore accessibility is low because most of the pores are not interconnected and only the pores at the surface are open to fluid adsorption, at least on the probed time scale of 1 ns. The amount of adsorbed molecules might still increase with longer simulation times, if the diffusion process is too slow to be captured in 1 ns simulations. In addition, we emphasize that diffusion tends to be underestimated when the flexibility of the kerogen matrix is not taken into account.32,33 In other words, the results in Figure 6 give a lower limit for the pore space accessibility. Nevertheless, in CS1000a, we find the same density with both simulation methods, showing that all pores in CS1000a are accessible on a very short time scale attainable in molecular simulations. The small deviations (