Subcontinuum mass transport of condensed hydrocarbons in nanoporous media

Although hydrocarbon production from unconventional reservoirs, the so-called shale gas, has exploded recently, reliable predictions of resource availability and extraction are missing because conventional tools fail to account for their ultra-low permeability and complexity. Here, we use molecular simulation and statistical mechanics to show that continuum description—Darcy's law—fails to predict transport in shales nanoporous matrix (kerogen). The non-Darcy behaviour arises from strong adsorption in kerogen and the breakdown of hydrodynamics at the nanoscale, which contradict the assumption of viscous flow. Despite this complexity, all permeances collapse on a master curve with an unexpected dependence on alkane length. We rationalize this non-hydrodynamic behaviour using a molecular description capturing the scaling of permeance with alkane length and density. These results, which stress the need for a change of paradigm from classical descriptions to nanofluidic transport, have implications for shale gas but more generally for transport in nanoporous media.


Adsorption isotherms
The adsorption isotherms in Fig. 1c of the article were obtained in CB-GCMC simulations of alkanes in the carbon matrix at different fugacities (equivalent to chemical potential that is µV T -ensemble). To give the adsorption isotherms in terms of pressure instead of fugacity, the fugacity-pressure relation was determined by comparing bulk simulations in the µV T (CB-GCMC) and the N P T (MD) ensemble at the same density according to the method reported in 6 .
The adsorption isotherms can be described by a Langmuir form where ρ ∞ is the maximum possible density of alkanes, which is reached in the limit f → ∞ and κ is connected to the isosteric heat of adsorption Q st . Both ρ ∞ and Q st depend systematically on the alkane length n: • The isosteric heat of adsorption in the low pressures limit Q 0 st scales linearly with the alkane length, because every monomer contributes equally.
• The maximum mass density increases slightly with the alkane length (see Supplementary   Figure 2). In other words, the pore space can be filled a little denser with larger alkanes.
This somewhat counter-intuitive result is due to the relatively high flexibility of the alkanes, and to the fact that the intramolecular distance between two functional groups is smaller than for two groups of different molecules.
The continuous form of the adsorption isotherms shows that no gas-liquid transition (i.e. capillary condensation) occurs. For methane and propane this is expected, since T = 423 K is above the critical temperature for both fluids. This is not the case for the longer alkanes however, where under the same thermodynamic conditions bulk fluid undergoes a phase transition. The dependence of phase behavior on confinement, i.e. in interface dominated systems, is a well-known effect. 9

Molecular dynamics simulations
Linear response. A linear relationship was observed for all cases, that is to say for all considered n-alkanes (up to dodecane), pressure gradients (up to ∼ 10 16 Pa/m) and flow directions. As it is typically the case in MD simulations, the pressure gradients were far larger than in realistic systems, excluding non-linear effects in a far larger range than the relevant range of pressure gradients. Concerning the flow direction, the driving force was generally applied along the z-direction of the carbon matrix. But for one dense system of dodecane, flow in all three dimensions was investigated and compared. No differences were observed for the transport properties. This is expected because the structural characteristics of the CS1000a model are isotropic.
Non-Darcy behavior. Based on Darcy's law, the permeance K = −q/∇P is expected to behave like k/η, where η is the fluid viscosity and k a material constant of the matrix -its permeability. In other words, the product K × η should be constant, independently of the fluid type, thermodynamic conditions and applied pressure gradient. Fig. 2 in the article shows that this is not the case here. In order to describe the fluid flow with a Darcy-like expression we would have to assume a fluid and density dependent permeability. One reason for this is the different phase behavior of the confined fluid (see previous discussion of adsorption properties). In the traditional description, adsorption effects are not taken into account. The viscosity of the fluid inside the pores is assumed to be the same as for a bulk fluid under the same thermodynamic conditions.
However, the density of fluids in nanopores can be considerably different than in bulk. Therefore, one might consider taking these adsorption effects into account by using the viscosity of the fluid phase in the pores. However, the definition of viscosity in this disordered and confining pore network is highly non-trivial, 11 which is the reason why we abstain from calculating a confined fluid viscosity. Instead, we make a first estimate from equilibrium simulations of bulk fluid at the same temperature and density as in the pores (NVT ensemble) using the Green-Kubo 10 relation with P ij an off-diagonal element of the stress tensor. This procedure implies that the simulated bulk phase is not strictly at the same pressure as the confined fluid. Interestingly, we find that the temperature and pressure dependence of the viscosity is small compared to the influence of density, meaning that fluids at different (T, P )-conditions, but at the same ρ, will have very similar viscosity -independently of the alkane length. As shown in Supplementary Figure 3, the viscosity can be estimated according to offers some valuable information about the dynamic behavior of a fluid. ¿From the Navier-Stokes equation for a simple incompressible fluid it follows analytically that the long-time behavior of the transverse current autocorrelation function is an exponential decline with a decay parameter that depends quadratically on the mode k, and linearly on the kinematic viscosity ν = η/ρ. Adsorbed inside a solid matrix, we would then expect a relaxation of the form exp [− (γ 0 + k 2 ν)], where γ 0 stems from the friction of the liquid with the solid matrix. This hydrodynamic behavior is found for none of the studied systems. Instead of the simple exponential decay, we find a double exponential of the form and the coefficients a and b do not scale with k 2 , but converge to constant values for large k (see Fig. S4). This is a clear indication that the n-alkane dynamics cannot be described by the classical hydrodynamics equations for a simple viscous fluid. Note that the observed form of the transverse current autocorrelation function does not agree with a visco-elastic behavior either. The deviations from the expected exponential form are more pronounced with increasing alkane length, hinting at the confinement as the origin of the failure of the classical hydrodynamic theory in the studied systems.
Self and collective diffusion. More information about the macroscopic transport properties can be obtained from studying diffusion processes. Indeed the permeance can also be expressed in the form with D 0 a collective diffusion coefficient. The second equality is the fluctuation dissipation theorem for the collective diffusion, relating D 0 to the time correlation function of the velocity fluctuations of the fluid center-of-mass with respect to the (frozen) carbon matrix q(t) = 1 N l v (l) (with v (l) the velocity of molecule l). The factor 3 stems from averaging over three coordinates, which is possible since the transport properties were found to be isotropic. The above relation offers an independent possibility to determine the permeance in equilibrium MD simulations. Fig. S5(a) confirms that results from both methods, equilibrium fluctuations and steady-state flow, agree very well. The collective diffusivity is related to the molecular self diffusivity as follows: It turns out that the remaining cross-correlation part (the integral in Eq. (15)) is often small compared to the self-correlation part D s , as can be seen in Fig. S5(b) and Fig. 4(a) in the article.
The self diffusion of one alkane molecule in a bulk fluid is in agreement with the Stokes-Einstein relation with slip conditions for a particle with effective diameter 2R 0 . Putting in the values for the self diffusion and the viscosity (obtained in equilibrium MD simulations, see discussion above), we find that the effective diameter for the whole alkane chain is about the size of one monomer 2R 0 ≈ 2 1/6 σ CH 2 ≈ 4.2Å (shown in Fig. S6), as is expected 12 . Hence, contrary to the confined case, the length of the alkanes does not play an explicit role for the diffusion in bulk. The origin of the n-dependence of the diffusion in the matrix is the interaction of the alkane chain with the matrix, which scales with the length.

Free volume calculation
In order to estimate the free volume V f ree that appears in Fig. 4b and Eqs. (8)(9)(10) in the article, we performed a numerical calculation of the void space between hard spheres with diameter σ (see Table 1) for representative alkane/CS1000a configurations at different loadings Γ. This calculation of the void space (=free volume) was done by randomly choosing a large number of sample points N (typically N = 50000) inside the simulation box, and checking for each of them if the coordinates lie inside a sphere with σ/2-radius around any solid or fluid atom, or not. In the latter case, the point is counted as lying in the void. The free volume is estimated from the ratio of the number of points in void space to the total number of sample points: V f ree /V = N void /N .

Molecular Models
In order to assess the transport and adsorption properties of n-alkanes in amorphous nanoporous carbons, we performed molecular dynamics (MD) and Grand Canonical Monte Carlo simulations with an insertion and configuration bias (CB-GCMC), respectively. In the following, we describe the models that were used for the alkane molecules and for the carbon matrix.
Porous carbon structure. A realistic model of pyrolized (at 1000 • C) and activated saccharose was used for the carbon structure. The atomistic structure of this model, referred to as imbibition from an outside bulk reservoir. Supplementary Figure 1 shows that the pore size distribution for this numerical sample spans from a fewÅ to ∼ 15Å, which is consistent with the pore sizes probed by CO 2 adsorption in kerogen.
Alkane model. The n-alkanes were modeled with a united atom force field 3 . Bonds between two coarse grained beads and bending between two neighboring bonds are constrained by harmonic where r ij = |r i − r j | is the distance between the two interacting atoms. The parameters ε ij and σ ij are calculated with the Lorentz-Berthelot mixing rules ε ij = √ ε i ε j and σ ij = 1 2 (σ i + σ j ), 4 and the in the CS1000a structure are frozen. Consequently, the results for the alkane flow rates presented in this study are lower limits, because the alkane molecules might deform a flexible structure and thereby increase their mobility. We assume that the effect of flexibility is a systematic one, and that it is small compared to the influence of the alkane size.

Configurational biased Grand Canonical Monte Carlo simulations
Adsorption simulations were performed with the Grand Canonical Monte Carlo technique (µV T ensemble). Monte Carlo steps consisted of translation, rotation, partial regrowth, insertion and deletion of alkane molecules. Due to the relatively large size of the molecules (up to dodecane) and the high density of the system (pore diameters of the order of the molecular size), a biased procedure was used for the insertion and for the regrowth of molecules in order to enhance the probability to find an energetically favorable configuration. Otherwise, the acceptance rate would be impractically low. The applied bias procedure follows the scheme introduced by Smit and can be found in 5 .

Molecular Dynamics
The MD simulations were performed with the LAMMPS package 7 . Integration of the equations of motion was performed with a 1 fs time step. Simulations were performed in the NVT ensemble.
Initial configurations with different alkane densities were taken from the CB-GCMC adsorption simulations. The pressures corresponding to these densities were obtained from a combination of GCMC and MD simulations of bulk fluids. A Nosé-Hoover thermostat with a relaxation time of 0.1 ps was applied to keep the temperature of the fluid at 423 K, a typical value for unconventional reservoirs (From a practical point of view a high temperature has the advantage to accelerate the particle dynamics and thus to equilibrate the system faster). For the non-equilibrium simulations the center of mass velocity was subtracted from the atom velocities before the rescaling. In some simulations with low fluid density the thermostat was applied to a relatively small number of molecules (around 30), which might result in a spurious effect of the thermostat on the flow rate.
We checked that this is not the case by replicating the dodecane system with the lowest density in every dimension (eight times larger volume) and repeating measurements. No influence of the system size on the fluid transport properties was observed.
Flow rates were measured for different driving forces by applying an external gravitational field g of different magnitude. Applying a constant acceleration g is representative for different setups where flow is induced by an external driving force: In particular, it is equivalent to applying a pressure gradient where ρ is the fluid density. Or, instead of a pressure gradient, the driving force could also be a chemical potential gradient -both are related by the Gibbs-Duhem equation ρ N dµ = dP at constant temperature, as was also verified in simulations by Arya et al. 8 In the following, results are presented in terms of pressure gradients, in order to make a direct connection to the extraction of hydrocarbon fluids from an underground reservoir, and to Darcy's law. ¿From the steady-state flow simulations, we obtain a value for the mean fluid flow velocity in the carbon matrix. We studied this mean flow velocity as a function of the driving force (i.e. pressure gradient), the alkane length, and the fluid density. The temperature was kept constant at