Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability

    Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability Christopher J. Landry, Maˇsa Prodanovi´c, Peter Eichhubl PII: DOI: Reference:

S0166-5162(16)30065-9 doi: 10.1016/j.coal.2016.03.015 COGEL 2611

To appear in:

International Journal of Coal Geology

Received date: Revised date: Accepted date:

7 January 2016 22 March 2016 23 March 2016

Please cite this article as: Landry, Christopher J., Prodanovi´c, Maˇsa, Eichhubl, Peter, Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability, International Journal of Coal Geology (2016), doi: 10.1016/j.coal.2016.03.015

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Revised Main Text to International Journal of Coal Geology March 17, 2016

Authors:

PT

Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability Christopher J. Landry*,a, Maša Prodanovićb, Peter Eichhublc

a

RI

Center for Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX 78712, USA b

SC

Department of Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX 78712, USA c

MA

NU

Bureau of Economic Geology, Jackson School of Geosciences, The University of Texas at Austin, Austin, TX 78757, USA

*

AC CE P

TE

D

Corresponding Author: Christopher J. Landry Center for Petroleum and Geosystems Engineering The University of Texas at Austin 1 University Station C0304 Austin, TX 78712 Email: [email protected] Phone: 248-420-4383 Maša Prodanović

Department of Petroleum and Geosystems Engineering The University of Texas at Austin 200 E. Dean Keeton, Stop C0300 Austin, TX 78712-1585 Email: [email protected] Phone: 512-471-0839 Peter Eichhubl Bureau of Economic Geology The University of Texas at Austin University Station, Box X Austin, Texas 78713-8924 Email: [email protected] Phone: 512-475-8829

1

ACCEPTED MANUSCRIPT Abstract:

D

MA

NU

SC

RI

PT

Mudrocks and shales contain pores within the size range of 2 – 100 nm. Flow of supercritical gas in these pores at reservoir pressure-temperature conditions falls within the slip-flow and early transition-flow regime (0.001 < Kn < 1.0). Currently, the description of supercritical gas flow in these flow regimes is mostly limited to simple tube and channel geometries that are of limited applicability to the sponge-like or platy nanoporous geometry in organic matter or clays. Here, we present a local-effective-viscosity multirelaxation-time lattice Boltzmann model (LEV-LBM) designed to simulate gas flow in the slip- and earlytransition-flow regimes in complex geometries. The LEV-LBM is informed with local effective viscosities at each node to capture the variance of the mean free path of gas molecules in a bounded system. The corrected mean free path for each lattice node is determined using a three-dimensional wall function adaptable to complex pore geometries. To enforce a non-zero slip velocity at solid boundaries, a combined diffusive bounce-back scheme is applied to the pore-walls. The LEV-LBM is first validated in simple tube geometries, where good agreement is found for Knudsen numbers up to 1.0. We then use the LEV-LBM to quantify the finite tube length effect and comment on the implications for pore-network models. We finally demonstrate the utility of the LEV-LBM by simulating pure methane flow in digital reconstructions of nanoporous organic matter at reservoir conditions, and compare the results to bundle of tubes models. We show that the bundle of tubes models overestimate apparent permeability by factors between 1.52 and 153, due to the non-trivial dependence of flow on pore space connectivity and shape.

1. Introduction

TE

Keywords: slip-flow, lattice Boltzmann, mudrock, shale, pore-scale, mesoscopic model

AC CE P

Solid bitumen/pyrobitumen, and organic components of thermally mature mudrocks host nanometer-scale pores (nanopores) that provide pathways for gas flow through mudrocks (Chalmers et al., 2012; Heath et al., 2011; Jiao et al., 2014; Klaver et al., 2015; Loucks et al., 2012; Milliken et al., 2013; Peng et al., 2015; Wood et al., 2015). At reservoir pressure-temperature conditions the mean free paths of supercritical gases, such as methane and carbon dioxide, are at or near the same order in length as the size of the pores. This causes the flow of these gases to significantly deviate from classic Navier-Stokes predictions, resulting in an increase in the apparent permeability of the porous media with decreasing pore pressure. Although this rarefaction effect on, or non-Darcy behavior of, the apparent permeability of porous media is well-known (Klinkenberg, 1941), the prediction of apparent permeability as a function of pore space geometry and gas properties remains elusive. The deviation of gas flow from the classic Navier-Stokes description is qualitatively described by the flow-regime, determined by the Knudsen number, , where is the unbound mean free path of the gas and is a reference length (i.e. the radius of a tube or aperture of a slit). A indicates the onset of non-negligible deviation from Navier-Stokes predictions, with flow regimes being divided up into slip, transition and free-molecular for Knudsen numbers between, 0.001 to 0.1, 0.1 to 10, and greater than 10, respectively. These flow regimes are roughly divided up in this manner to reflect the changing physics as the Knudsen number increases (Colin, 2005; Zhang et al., 2012a). In the slip-flow regime the gas-gas collisions are far greater in number than the gas-wall collisions, however, the gas-wall collisions result in a significant non-zero flow velocity at the walls in comparison to the mean flow velocity, hence the term “slip”. In the transition-flow regime there is a transition from gas-gas collision to gas-wall collision dominance, with both having a significant effect on the overall flow. In the free2

ACCEPTED MANUSCRIPT

MA

NU

SC

RI

PT

molecular-flow regime the gas-wall collisions become increasingly dominant to the extent that the gasgas collisions can be ignored, approaching flow that is described by Knudsen diffusion. It is somewhat common to refer to these flow regimes as non-continuum, however, it is important to note that although these flows deviate from classic Navier-Stokes description, there is no failure of the continuum assumption (Hadjiconstantinou, 2006; Guo et al., 2007). For these flows “the hydrodynamic fields and the associated conservation laws remain well-defined” (Hadjiconstantinou, 2006). The molecular velocities are far greater than the flow velocities of interest here, and subsequently the time scale of individual intermolecular collisions is far smaller than the characteristic time scale of steady-state flow. For steadystate flow the time-averaged microscopic behavior of the fluid in any volume of pore space results in well-defined velocity moments (flow velocity), this suggests that a continuum treatment is valid. Therefore, the deviation of supercritical gas flow in nanoporous media from the classic Navier-Stokes description can be accounted for using an “extended” Navier-Stokes constitution (Guo et al., 2007). Considering the gas flow in nanopores at reservoir pressure-temperature conditions is not actually „rarefied‟, includes flow outside the slip-flow regime, and the term non-continuum is misleading, we refer to this type of flow as mixed viscous-ballistic (MVB).

AC CE P

TE

D

Physical measurement of MVB flow in simple engineered geometries are typically conducted under rarefied conditions, as in, at low pressure, often with large inlet over outlet pressure ratios (e.g. Tison, 1993; Arkilic et al. 1994; Roy et al., 2003; Colin, 2005; Marino, 2009). These pressure conditions are significantly different from those experienced during hydrocarbon production. The prediction of apparent permeability as a function of pore space geometry in complex nanoporous media requires an investigation at the pore-scale. For the investigation presented here this results in a spatial scale of interest on the order of 100-1000 nm (thousands to millions of ~5 nm radius pores). At this length scale the difference between the upper and lower pressure bounds is very small in comparison to typical pore pressures at reservoir conditions. For example, a 100 MPa/m pressure gradient would result in an upper and lower pressure difference for a 100 nm volume of only 10 Pa. The pore-scale flow of supercritical gases in nanoporous media at the spatial scale of interest here is subject to minimal compressibility effects. Direct comparison of pore-scale modeling results, such as the one presented here, to these physical measurements would require an investigation into these compressibility related effects in the model. The flow conditions of interest of this investigation are limited to theoretical work (e.g. Loyalka and Hamoodi, 1990, Sharipov and Seleznev, 1994), due to the difficulty of physically measuring either pressure difference or flow rate while minimizing compressibility effects (Marino, 2009). Predicting the apparent permeability of mudrocks as a function of pore space geometry, pressuretemperature conditions and gas composition for MVB flows has recently received ample attention in the hydrocarbon recovery literature (e.g. Javadpour, 2009; Civan, 2010; Darabi et. al., 2012; Sakhee-Pour and Bryant, 2012; Mehmani et al., 2013; Rahmanian et al., 2013; Yao et al., 2013; Islam and Patzek, 2014; Ma et al., 2014; Singh et al., 2014; Kazemi and Takbiri-Borujeni, 2015; Naraghi and Javadpour, 2015). The over-arching goal is predicting long-term production outcomes from the characterization of mudrock pore space geometry. However, all of these models simplify these pore spaces into a collection of constant cross-section tubes or slits, by use of bundle of tubes models which treat pores as a collection isolated tortuous tubes (e.g. Civan, 2010; Darabi et. al., 2012; Naraghi and Javadpour, 2015) or by use of pore network models which treat pores as a network of constant cross-section conduits (Sakhee-Pour and Bryant, 2012; Mehmani et al., 2013; Ma et al., 2014). Unlike the bundle of tubes models, pore network models allow for the investigation of the dependence of flow on pore space connectivity, including 3

ACCEPTED MANUSCRIPT

RI

PT

interconnectivity in heterogeneous networks with bimodal pore size distributions (Mehmani and Prodanović, 2014). The tube/slit models employed by bundle of tubes models and pore network studies were developed for very long tubes, with lengths approximately 20 times larger than the radius of the tube, and often far larger. Unlike for viscous flow, the permeability of a tube for MVB flows is dependent on the length of the tube and pressure gradient as a result of significant end effects (Beskok and Karniadakis, 1999; Colin, 2005; Marino, 2009). None of the above referenced studies considers the effect of finite length on flow at reservoir pressure temperature conditions and we address that issue in this work.

AC CE P

TE

D

MA

NU

SC

Furthermore, the application of tube/slit models for accurate prediction of flow in the complex pore geometries of mudrocks is still an open research question. From scanning electron microscopy (SEM) and transmission electron microscopy (TEM) studies it has been shown that that the pore spaces of mudrocks are best described as “spongey” for organic matter hosted pores or as a complex collection of variable cross-sectional shape for clay hosted pores. It is also evident from imaging studies that the porosity of organic matter varies widely, and is often far greater than the bulk rock porosity (Chalmers et al., 2012; Heath et al., 2011; Jiao et al., 2014; Klaver et al., 2015; Loucks et al., 2012; Milliken et al., 2013; Wood et al., 2015). The tools for simulating supercritical gas flows in nanoporous media in anything but very long tubes, or very long/wide slits has been limited for the most part to molecular dynamics modelling (e.g. Cao et al., 2006; Falk et al., 2015; Firouzi and Wilcox, 2012). The computational requirement for molecular dynamics model prediction of flow increases with the number of molecules being tracked. Molecular dynamics flow simulations are generally limited to very small physical simulation volumes, on the order of 10x10x10 nm3 (Firouzi and Wilcox, 2012; Falk et al., 2015). The complex and often heterogeneous pore spaces of organic matter pore systems in mudrocks may require flow simulations on volumes orders of magnitude larger than this to determine relations between pore space geometry and macroscopic transport properties. However, molecular dynamics simulations of small volumes remain very useful, particularly when surface-gas interaction forces dominate mass flux. Under such scenarios, molecular dynamics simulations offer insight into, and estimations of, surface transport properties (Firouzi and Wilcox 2012, Falk et al. 2015), and can be used in the validation of mesoscopic modeling methods. For MVB flows both local variation of momentum diffusivity and microscopic gas-wall interactions (slip velocity at the wall) must be accounted for, and the mesoscopic approach of lattice Boltzmann (LB) methods has long been considered a natural candidate for the prediction of such flows, at minimal computational cost in comparison to molecular dynamics simulation. Also, LB methods for single phase flows calculate velocity fields explicitly and locally, and are therefore highly scalable on parallel computational resources. A large body of work has been generated attempting to construct LB methods for prediction of MVB flows and great success has been found using a variety of approaches (e.g. Ansumali and Karlin, 2002; Chai et al., 2010; Guo et al., 2006; Guo et al., 2008; Kim et al., 2008; Li et al., 2011; Lim et al., 2002; Liou and Len, 2014; Liu and Guo, 2013; Nie et al., 2002; Niu et al., 2007; Sbargaglia and Succi, 2005; Shen et al., 2004; Succi, 2002; Suga and Ito, 2011; Suga et al., 2010; Verhaege et al., 2009; Zheng et al., 2012). The vast majority of these studies have focused on simple geometries (tubes/slits), with a limited number investigating more complex geometries. Notably, the work of Suga and Ito (2011) showed promising results for three dimensional flows involving obstacles in ducts, and bumpy channels. Similar success was described in the work of Chai et al. (2010) for two-dimensional flow through square arrays of circular cylinders, using a LB model that is partly the basis of the model 4

ACCEPTED MANUSCRIPT

RI

PT

described here. However, in both of these investigations the methods used treated the Knudsen number as a constant. For the porous systems studied by Suga and Ito (2011) and Chai et al. (2010), assuming a constant Knudsen number was acceptable, because the pore spaces, although complex in shape, were similar in size throughout the simulated volumes. Flow in the nanoscale pore geometry of mudrocks presents an issue that has not yet been resolved, construction of a three-dimensional LB model for MVB flows through complex geometries with locally varying pore size, and subsequently a locally varying Knudsen number.

D

MA

NU

SC

Here we adopt the approach of Guo et al. (2007), for which a geometry-dependent local effective viscosity (LEV) is defined to solve the “extended Navier-Stokes constitution”. In Guo et al. (2008) this basis of modelling MVB flows was employed by the multiple-relaxation-time lattice Boltzmann model (MRT-LBM) using Stops (1970) one-dimensional wall function to determine local viscosity for twodimensional parallel-plate flows. Here, we extend this approach to three dimensions for geometries of arbitrary shape. The concept is fairly straightforward, at the wall the mean free path of a gas approaches zero and increases as a function of the distance from the wall until it is equal to the unbound mean free path. By determining the local effective mean free path a local effective kinematic viscosity can be defined, and MVB flow can be modeled using extended Navier-Stokes hydrodynamics in combination with a slip boundary condition. The LEV method can also be implemented using the single-relaxationtime (SRT) LBM, although it should be noted that for a SRT-LBM the LEV method described here may be inaccurate due to discrete effects (Guo and Zheng, 2008; Pan et al., 2006).

AC CE P

TE

Lattice Boltzmann methods have been used to investigate “slip” flows in complex porous media recently (Allan and Mavko 2013; Chen et al. 2015), however these studies employ standard LB methods for pure viscous, no-slip flow to determine the absolute permeability, then apply a slip correction adapted from tube models to determine apparent permeability. Unlike the model here which implicitly includes the ballistic effect in the physics of the simulated fluid, and validates the resulting flow against more rigorous modeling methods at the local level. The investigation presented here has two main objectives. First, to introduce and validate our local effective viscosity lattice Boltzmann model (LEV-LBM). Second, to test the validity of bundle of tubes models for the prediction of apparent permeability as a function of pore geometry at reservoir pressuretemperature conditions. In the following, we review tube models for MVB flows, describe our LEVLBM, and introduce the three-dimensional digital reconstructions of nanoporous organic matter that will act as testing grounds for the bundle of tubes models. We then validate the LEV-LBM by comparison to the Beskok and Karniadakis (1999) tube model, and explore the tube length effect on MVB flows using the LEV-LBM. Finally, we test the utility of bundle of tubes models for the prediction of apparent permeability, by comparing direct MVB flow simulation using the LEV-LBM in three-dimensional nanoporous organic matter reconstructions with the prediction of bundle of tubes models. 2. Methods 2.1. MVB Flow Tube Models Here we critically review three tube models for the prediction of apparent permeability as a function of pore geometry, gas type and pressure-temperature conditions, namely, the dusty-gas model (DGM), the Brown et al. (1946) first-order slip-correction model, and the Javadpour (2009) model. In addition we 5

ACCEPTED MANUSCRIPT review the Beskok and Karniadakis (1999) tube model which is fitted to the linearized Boltzmann solutions Loyalka and Hamoodi (1990), and will be used in the validation of the LEV-LBM. First we review the flow regime of interest in nanoporous media at reservoir pressure-temperature conditions.

(1)

SC

RI

PT

At reservoir pressure-temperature conditions supercritical gas flow in nanoporous media falls within the slip- and early-transition-flow regimes. In Figure 1 we plot the Knudsen number for tubes within the typical size range of nanopores in mudrocks for methane flow at 400 K, and typical reservoir pressures. The unbound mean free path is calculated using the hard-sphere model of the kinetic theory of gases,

TE

D

MA

NU

where is the hard-sphere diameter of methane ( ), and is the number density of methane molecules. At reservoir PT conditions methane phase behavior significantly deviates from that of an ideal gas. Therefore, methane mass density, , is calculated using the Peng-Robinson equation of state (Peng and Robinson, 1976), with , where is Avogadro‟s number, and is molar mass. The majority of the flow within mudrock nanopores falls within the slip- and early-transition-flow regimes, with the exception of nanopores less than 1 nm in radius. For pores of this size the surface-gas interaction potential is expected to dominate gas transport (Wu et al. 2015). This also indicates that at reservoir pressure-temperature conditions, the flow is never free-molecular. While we do not expect MVB flow physics alone to provide a predictive model of transport for pores of this size, we will nevertheless discuss the possibility of inclusion of gas adsorption and surface diffusion in the LEV-LBM framework.

AC CE P

The dusty-gas model (DGM) treats MVB flows as a linear summation of mass flux predicted for viscous flow, , and the mass flux of Knudsen diffusion, , (Mason and Malinauskas 1983; Sakhee-Pour and Bryant 2012), (2)

where is the dynamic viscosity. Knudsen diffusion is free-molecular or fully ballistic flow, and the Knudsen diffusion coefficient, , is a function of the tube radius, , the specific gas constant, , and temperature, , (Roy et al. 2003),

(3)

The apparent permeability for the DGM is then given by, (4) The Brown et al. (1946) tube model applies a slip correction, , to mass flux predicted for viscous flow, . Using the first order slip-correction to the Navier-Stokes equations, the apparent permeability for the Brown et al. (1946) tube model,

6

ACCEPTED MANUSCRIPT (5)

SC

RI

PT

where is a tangential momentum accommodation coefficient, which is based on microscopic arguments regarding gas-surface interactions (Maxwell, 1879) – with representing purely diffusive reflection and for purely specular reflection. Under rarefied conditions, tangential momentum accommodation coefficients as low as 0.2 on polycrystalline and molybdenum surfaces have been observed (Beskok and Karniadakis, 1999; Agrawal, 2008). However, we would argue that the tangential momentum accommodation coefficient is closer to 1 for the organic surfaces of interest here, which are likely to contain a significant amount of molecular-scale structural disorder that would inhibit specular reflection.

NU

The Javadpour (2009) tube model combines the DGM and the slip-flow correction of Brown et al. (1946), , which when expanded gives,

MA

(6)

AC CE P

TE

D

Whereas the DGM and the Brown et al. (1946) tube model were both developed initially from empirical observations, a general tube model developed for prediction in all flow regimes requires a more rigorous treatment of MVB flow, optimally, solving the Boltzmann equation for flow in tubes. Solutions to the Boltzmann equation are limited to numerical methods, and with the advent of greater computational resource availability, Loyalka and Hamoodi (1990) used numerical integration techniques to solve the linearized Boltzmann equation for a hard sphere gas flow in an infinite tube with a fully diffusive tangential momentum accommodation coefficient. The results of Loyalka and Hamoodi (1990) were then used by Beskok and Karniadakis (1999) in the development of a tube model applicable to all MVB flow regimes. The Beskok and Karniadakis (1999) tube model employs a general slip-velocity boundary condition to the Navier-Stokes equations, and a rarefication coefficient, , which describes the decrease in the number of gas-gas collisions compared to gas-wall collisions, (7)

when fitted to Loyalka and Hamoodi (1990), (8) The velocity profile results of Loyalka and Hamoodi (1990) were then used by Beskok and Karniadakis (1999) to verify a local velocity function, (9) where is the local velocity, is the distance from the center axis of the tube, and velocity, which is related to apparent permeability in accordance with Darcy‟s law,

is the mean

7

ACCEPTED MANUSCRIPT (10)

PT

This tube model fits the velocity profiles of the linearized Boltzmann solutions of Loyolka and Hamoodi (1990) with the exception of a small error near the wall at (Beskok and Karniadakis 1999, their Figure 7).

AC CE P

TE

D

MA

NU

SC

RI

To compare the tube models reviewed for typical reservoir pressure-temperature conditions to the more rigorous linearized Boltzmann solution of Loyalka and Hamoodi (1990) and the Beskok and Karniadakis (1999) tube model we present results for pure methane flow in an infinite tube with a radius of 5 nm at a temperature of 400 K (Figure 2). We determined the density of the gas using the Peng-Robinson equation of state, and the unbound mean free path using the hard-sphere gas model. The dynamic viscosity of methane was determined from the formulation of Lee (1966). In Figure 2, pressures of 5 and 100 MPa have Knudsen numbers of 0.344 and 0.0270, respectively. The Brown et al. (1946) tube model using a fully diffusive tangential momentum accommodation coefficient, , agrees with the Beskok and Karniadakis (1999) tube model fairly well, with an increasing underestimation for larger Knudsen numbers (lower pressures). When the tangential momentum accommodation coefficient is set for fully diffusive reflection instead of being treated as a fitting parameter, as it was in Brown et al. (1946) and Javadpour (2009), the Brown et al. (1946) tube model is only applicable to the slip-flow regime, underestimating flow at higher Knudsen numbers (Beskok and Karniadakis, 1999). In contrast, the DGM largely overestimates flow. The DGM assumes Knudsen diffusion is occurring throughout the tube, however, the ballistic effect produces something akin to Knudsen diffusion only within the Knudsen layer, outside of which it is unlikely a molecule will not undergo gas-gas collisions before the next gaswall collision. The Knudsen layer has a thickness on the order of the mean free path. Therefore, in the slip-flow regime, the Knudsen layer would occupy between 0.2% and 20% of the pore space, resulting in the DGM overestimating flow in the slip-flow and early-transition-flow regimes. For the Javadpour (2009) tube model we use a value of 1.0 for , instead of 0.8 as fitted by Javadpour (2009) to the experimental results of Roy et al. (2003). Regardless of the value of , which has a theoretical limit of 1, the Javadpour (2009) model overestimates flow. The Javadpour (2009) tube model inherits the same overestimation associated with the DGM for the slip-flow and early transition-flow regimes, exacerbated by the inclusion of the Brown et al. (1946) slip correction. It is essentially a double correction. Unlike the advection-diffusion equation for solute transport, in which two separate transport mechanisms act to transport the solute, slip-flow and Knudsen diffusion are not separate transport mechanisms that both need be accounted for to describe MVB flow. Slip and Knudsen diffusion are the result of the same transport mechanism – increasing ballistic effect – acting to change the nature of momentum diffusivity in the system with an increasing degree of severity. The variance in momentum diffusivity is most significant near the walls (i.e. the Knudsen layer). This variance in momentum diffusivity has an increasing effect on the overall flow as the Knudsen layer occupies a larger proportion of the system volume. In the free-molecular-flow regime the flow is almost completely ballistic and is well described by diffusion alone, namely, Knudsen diffusion. The differences in the tube model results were previously reported by Yao et al. 2013, however, without the above discussion of the mechanisms responsible. 2.2. MVB Flow Bundle of Tubes Models All of the bundle of tubes models assume an absolute permeability of, 8

ACCEPTED MANUSCRIPT (11)

RI

PT

where is connected porosity, is the mean radius, and is the tortuosity. In Narghi and Javadpour (2015) the tortuosity was given an assumed value of 2, and here we adopt this value. The three models are, the DGM bundle of tubes model, (12)

(13)

NU

SC

the Civan et al. 2010 bundle of tubes model utilizing the Beskok and Karniadakis (1999) tube model,

MA

and the Naraghi and Javadpour (2015) bundle of tubes model utilizing a modified version of the Javadpour (2009) tube model, (14)

(15)

AC CE P

TE

D

where is the fractal dimension of the pore surface (set in Naraghi and Javadpour (2015) to 2), and the tangential momentum accommodation coefficient is treated as a function of the Knudsen number (Agrawal and Prabhu, 2008),

2.3. Local Effective Viscosity Lattice Boltzmann Model In the following we describe our local-effective-viscosity (LEV) LBM for MVB flows. The LEV-LBM uses the same framework as conventional LB methods with the exception of a locally-defined kinematic viscosity and a combined diffusive bounce-back scheme along solid boundaries (Guo and Zheng, 2008; Chai et al., 2010; Suga and Ito, 2011). Details of the conventional MRT-LBM adapted here to MVB flows can be found in the literature (e.g. Qian et al., 1992; d‟Humieres, 2002; Guo and Zheng, 2008; Chai et al., 2010; Suga and Ito, 2011). All simulations were carried out by a parallelized in-house code. Lattice Boltzmann models are mesoscopic models based on microscopic dynamics with locally defined hydrodynamic moments. In short, fluids are simulated as swarms of particles which flow on a discrete lattice, represented by particle distribution functions ( ). Here we use the D3Q19 discrete velocity model (Qian et al., 1992), wherein the particle distribution functions are limited to 19 discrete velocities (18 advecting, 1 rest). Each discrete particle distribution contains the value of the density of the particles moving in its respective direction. During a single time iteration the LB equation is carried out in two steps, streaming of the particle distribution function to the neighboring lattice nodes, and collision of the particle distribution function. In the MRT-LBM the velocity-space particle distribution functions are transformed to moment-space prior to collision (d‟Humieres, 2002). This has a distinct advantage over colliding in velocity space (i.e. SRT-LBM), particularly for LBMs parameterized for MVB flows. In moment space the particle distribution function is represented by hydrodynamic and kinetic moments 9

ACCEPTED MANUSCRIPT which can be relaxed during the collision individually. Therefore, the relaxation time ( ) defining kinematic viscosity ( ),

PT

(16)

SC

RI

where is the square of the speed of propagation (“sound”), and is equal to 1/3 in the D3Q19 lattice used here, can be specified locally without interfering with the other kinetic moments. In the following description of the LEV-LBM the ability to parameterize the viscosity of the simulated fluid on a local level will be key in capturing the physics of MVB flows in complex nanoporous media. ) of a gas is

MA

NU

For the hard-sphere model of the kinetic theory of gases the unbound mean free path ( directly proportional to the unbound kinematic viscosity,

. From this we define

as a function of the mean free path,

(18)

TE

D

where for the MRT-LBM,

(17)

AC CE P

For the other relaxation times in the LEV-LBM a standard set of values defined for conventional NavierStokes flows are used (d‟Humieres, 2002), with the exception of the relaxation time for the energy-flux moment ( ), which is defined by the second-order slip boundary condition using the parameterization presented by Chai et al. (2010), (19)

where and are the coefficients of the second-order slip boundary condition, and following the work of Chai et al. (2010) are set to 1.11 and 0.61, respectively. At the solid boundaries (pore walls) the LEV-LBM uses a combined diffusive bounce-back scheme (CDBB), (20) where is the location of the solid boundary node, coefficient, also following Chai et al. (2010),

is the discrete velocity and

is the probability

(21)

There are two modes of bounce-back,

represents perfect-reflection or standard bounce-back, 10

ACCEPTED MANUSCRIPT (22)

PT

and represents redistributed or diffusive bounce-back (Ansumali and Karlin 2002) which retains the Maxwellian equilibrium particle distribution ( ) at the pore walls,

RI

(23)

SC

for static walls this reduces to,

, and

, where

NU

where is the sum of the incoming for the solid boundary node at is the sum of the D3Q19 weights, , of the incoming .

(24)

D

MA

In an unbound geometry the mean free path is a constant. In a bound geometry the mean free path is essentially zero at the walls and increases as a function of the distance from the pore walls as illustrated in Figure 3a and 3b (Guo et al., 2007; Zhang et al., 2012a). The Suga and Ito (2011) LBM employed the construction outlined above with a global effective viscosity defined by the infinite parallel wall bound effective mean free path correction ( ) of Guo (2006),

TE

(25)

where

AC CE P

Unlike which is limited to pore geometries with a near consistent , the novelty of the LEV-LBM presented here is the definition of a geometry-dependent local effective viscosity determined by the local effective mean free path, equations 26 to 31. The local effective mean free path, , is determined by the effective mean free path correction function ( ),

is location. Therefore, the effective viscosity also varies locally and

(26) is redefined,

(27)

as well as, . Here is the arithmetic mean of a one-dimensional wall function, wall measured in 18 directions for each lattice node (Figure 3d),

, for

a single planar

(28) where is the distance from the node to the wall in the specified direction. The one-dimensional wall function must meet two boundary conditions, , representing the vanishing free path at the location of the wall, and , the mean free path far away from the wall is equal to the unbound mean free path. These conditions can be met with an inverse tangent function, 11

ACCEPTED MANUSCRIPT

(29)

MA

NU

SC

RI

PT

where and are fitting parameters to be investigated during model validation. A plot of the onedimensional wall function using an illustrative range of possible and values is presented in Figure 3c. In a bound geometry there can also be a significant truncating effect on the molecular free path distribution, which at equilibrium includes a significant probability that the free path of an individual molecule is larger than the geometry. For the discrete wall function presented here, the size of the bounding geometry is equal to . Therefore, a fraction of the molecules, , should have an effective mean free path corrected for the size of the geometry not the unbound mean free path. This fractional construction results in the following “multi-bounce” one-dimensional wall function,

(31)

TE

D

and after some re-arranging,

(30)

2.4 Digital Reconstruction of Nanoporous Media

AC CE P

To demonstrate the utility of the LEV- LBM and the complex relationship between pore space geometry and flow properties, as well as to create a realistic testing ground for the bundle of tubes models, we generate two three-dimensional digital reconstructions of nanoporous organic matter. Nanopores in solid bitumen/pyrobitumen and organic material in mudrocks are generally spherical in shape. Digital reconstruction for this type of pore geometry can be generated as a collection of randomly placed overlapping spheres. Both reconstructions are composed of spheres with a uniform radius distribution between 4.5 and 6 nm. The simulation volumes are 100x100x100 nm3, with a 0.5 nm resolution. For both digital reconstructions spheres are randomly placed in a volume until a target porosity of 0.6 and 0.3 is reached, for reconstruction 1 and reconstruction 2, respectively. Not all of the randomly placed spheres are part of a connected pore network, to determine the connected porosity a percolation algorithm was applied to both reconstructions. Most of the overlapping spheres for the higher porosity target reconstruction 1 are connected, for reconstruction 2 many are not, resulting in a connected of porosity, , of 0.5832 and 0.1963, for reconstruction 1 and 2, respectively (Figure 4). 2.5 Determination of Mean Pore Radius Determining the mean pore radius of nanoporous organic matter can be accomplished by either imagebased methods or bulk measurements, such as mercury intrusion porosimetry (e.g. Sakhaee-Pour and Bryant, 2012; Naraghi and Javadpour, 2015; Peng et al., 2015). For an image-based determination of pore radius distribution we used a maximum inscribed sphere (MIS) algorithm, which in practice could be applied to a three-dimensional serial focused ion-beam SEM image of nanoporous organic matter, and is equivalent to applying a maximum inscribed circle algorithm to a two-dimensional SEM image (Figure 12

ACCEPTED MANUSCRIPT

RI

PT

5). To simulate a measurement made by mercury intrusion porosimetry we employ a full morphology model (Vogel et al., 2005) all-boundary flooding (Figure 6). Mercury intrusion porosimetry and the simulated flooding measure the pore throat distribution, therefore the pore radius distribution will be generally smaller than the pore radius distribution measured using the MIS algorithm. The pore radius distributions for both reconstructions using these two methods are presented in Figure 7. The mean pore radius measured using the MIS algorithm for reconstruction 1 and 2, is 5.084 and 4.821 nm, respectively. The mean pore radius measured using the full morphology model flooding for reconstruction 1 and 2, is 4.223 and 2.585 nm, respectively.

SC

3. Results 3.1. LEV-LBM Validation

AC CE P

TE

D

MA

NU

To validate our LEV-LBM we performed benchmark tube simulations and compared the results to the Beskok and Karniadakis (1999) tube model outlined in Section 2.1. As previously outlined, the Beskok and Karniadakis (1999) tube model is fitted to the linearized Boltzmann solutions of Loyalka and Hamoodi (1990), therefore the validation is the equivalent of a direct comparison to linearized Boltzmann solutions. Note that the LEV-LBM simulation is only informed of the local physics at every node by the effective-mean-free-path correction function, in other words, no individual lattice node is aware of the overall physical geometry of the system, the predicted flow fields emerge from the proper definition of local physics. Simulations were carried out on a tube with periodic boundary conditions (infinite length) with a radius equal to 50 lu, (lu refers to lattice units, 1 lu is equal to the distance between two lattice nodes), and lattice fluid density, , of ~1.0000. A pressure gradient is imposed by applying a body force, , chosen to ensure Stokes flow. Further simulations were then carried out for a resolution study on tubes with radii of 5 and 2.5 lu. The error for bulk flow comparisons between the LEV-LBM and Beskok and Karniadakis (1999) tube model are defined, , where is the mean velocity in the direction of flow predicted by the Beskok and Karniadakis (1999) tube model and is the mean velocity predicted by the LEV-LBM. Four illustrative examples of the fitting values for the one-dimensional wall function, as well as the multi-bounce wall function with and equal to 1 are presented in Figure 8a. In general, the local effective viscosity function with fitting parameters close to 1 are within 6% of the bulk flow predicted by the Beskok and Karniadakis (1999) tube model, with the wall function remaining within 2%. The function with and equal to 1 was then used in the resolution study, the results of which are presented in Figure 8b. As can been seen in Figure 8b the model is very robust, with very similar errors found for tubes as small as 2.5 lu in radius. To ensure that the model is not only reproducing bulk flow rates accurately, but also local velocity, we compare the velocity profiles of the LEV-LBM to the Beskok and Karniadakis (1999) tube model in Figure 8c-f. For the sake of comparison we have also included results using the global effective viscosity method of Suga and Ito (2011), which underestimates the flow in the transition-flow regime. Although it should be noted that this method cannot be readily applied to complex pore systems with significant variance in . The match between the LEV-LBM and the Beskok and Karaniadakis (1999) tube model are excellent in the slip-flow regime for both the single- and multi-bounce wall functions with and . There is however a small under-estimation of the centerline velocity for flows in the early-transition regime, with the multi-bounce wall function producing the most accurate profile for a Knudsen number of 1. The apparent permeability 13

ACCEPTED MANUSCRIPT of these infinite tubes, law,

, determined from the LEV-LBM simulations is calculated using Darcy‟s

where

, where

PT

(32) is the acceleration due to exterior force.

RI

3.2. MVB Flow in Finite Length Tubes

AC CE P

TE

D

MA

NU

SC

In the previous section the tubes were presumed infinite, for flow in tubes of finite length there is an end effect that results in the apparent permeability of the tube being less than that for an infinite tube. Here we demonstrate that the end effect can be attributed to a longer effective mean free path near the ends of the tube than within tube. At the inlet and outlet of a tube connected to reservoirs, where the gas is essentially unbound, there are more molecules originating from the reservoir at the ends than in the central portion of the tube. This difference is quantified by and illustrated in Figure 9a for flow through a tube with a length to radius ratio, , of 4, and . Shown in Figure 9a the effective mean free path at the center of the tube near the end is 0.7948 , while at the central portion of the tube it is 0.7658 . Thus, the effective viscosity at the inlet/outlet of the tube is greater than that at the center of the tube, resulting in a decrease in overall mean velocity, and . To simulate flow in finite tubes with the LEV-LBM a two-layer boundary is used (Figure 9b) with the pressure boundary enforced using the method of Hecht and Harting (2010). From these simulations the apparent permeability was calculated in the same manner as for the infinite tubes, except, . The apparent permeability normalized to infinite tube apparent permeability, , determined by the LEV-LBM as a function is presented in Figure 9c. For pore network models of complex nanoporous media the above simulation set up is analogous to that of a single throat connecting two pore bodies. However, the pore bodies of pore network models have a defined size that is often similar in magnitude to the size of pore throats, which would act to bound the gas in the „reservoirs‟. Therefore, this simulation setup represents the maximum finite length effect for a single throat connecting two pore bodies in a pore network model. The finite length effect on the apparent permeability of MVB flow in finite tubes is fairly small in the slip-flow regime for tubes with a of 2 or larger, with only dipping to 0.928 and 0.806 at of 0.4 for of 0.02 and 0.2 respectively. In the transition flow regime the length effect is much greater, with a of 0.323 at of 0.4 for of 1.0, and retaining a of 0.907 and 0.944 at of 10 for a of 0.5 and 1.0, respectively. Not shown in the plot of Figure 9c, the length effect for of 1.0 was still significant at of 160, with a of 0.963. The length effect can be significant in both the slip-flow and early transition-flow regimes, and should be considered in the construction of pore network models designed to quantify MVB flow. 3.3. Flow in Nanoporous Media To compare the flow results of the LEV-LBM to bundle of tubes models, and to illustrate the shortcomings of such models, we compare LEV-LBM flow simulation in nanoporous organic matter reconstructions 1 and 2 to the three bundle of tubes models reviewed in subsection 2.2 “MVB Flow Bundles of Tubes Model”.

14

ACCEPTED MANUSCRIPT

SC

RI

PT

The simulation set up for flow in nanoporous media is the same as that for flow in finite tubes (Figure 9b) without the inclusion of reservoirs at the inlet and outlet which effect the calculation of . Therefore the nanoporous media are treated as continuous beyond the inlet and outlet, as in, not connected to unbounded reservoirs and not subjected to the finite length effect. A conventional MRT-LBM (viscous flow) – set B of Pan et al. (2006) – using standard bounceback (no-slip) was used to determine , while as a function of pore pressure was determined for pure methane flow at 400 K using the LEVLBM. The determined by a standard MRT-LBM and the bundle of tubes models for both measurements of the mean pore radius ( ) is summarized in table 1. The and apparent permeability normalized to the absolute permeability, , as a function of pressure are presented in Figure 10.

AC CE P

TE

D

MA

NU

The bundle of tubes model greatly overestimates for both reconstructions. For reconstruction 1 the bundle of tubes model predicts a that is 2.18x and 1.50x greater than that predicted by the conventional MRT-LBM, for measured using the MIS algorithm and full morphology model flooding, respectively. For reconstruction 2 the bundle of tubes model predicts a that is 495x and 142x greater than that predicted by the conventional MRT-LBM, for measured using the MIS algorithm and full morphology model flooding, respectively. The overestimation of by the bundle of tubes models for both measurements of the mean pore radius result in the bundle of tubes models overestimating for all of the pressures simulated for both reconstructions. The Civan et al. (2010) bundle of tubes model informed by measured using the full morphology flooding provides the closest estimate of for all pressures and both reconstructions. This model overestimates for reconstruction 1 by 1.96x and 1.52x at pore pressures of 4.06 MPa and 35.5 MPa, respectively, and overestimates for reconstruction 2 by 33.2x and 36.5x at pore pressures of 4.06 MPa and 35.5 MPa, respectively. The Naraghi and Javadpour (2015) bundle of tubes model informed by measured using MIS provides the farthest estimate of for all pressures and both reconstructions. This model overestimates for reconstruction 1 by 2.98x and 4.95x at pore pressures of 4.06 MPa and 35.5 MPa, respectively, and overestimates for reconstruction 2 by 153x and 140x at pore pressures of 4.06 MPa and 35.5 MPa, respectively. The for reconstructions 1 was overestimated by a smaller factor than at all pressures simulated due to the overestimation of (Figure 10d) by all of the tube models. The for reconstructions 2 was overestimated by a larger factor than at all pressures simulated due to the underestimation of (Figure 10d) by all of the tube models. The overestimation of for both reconstructions and the underestimation of for reconstruction 2 by the bundle of tubes models can be attributed to the lack of consideration for the pore shape and pore space connectivity. In Figure 11 and 12 velocity flow fields from simulations using the standard MRT-LBM and LEV-LBM at 4.06 MPa are shown for the entire volume and in cross-section for reconstruction 1 and 2, respectively. From Figure 11a and 11b we can see that the flow in the higher porosity reconstruction is fairly uniform, and from Figure 11c and 11d we observe the largest increase in velocity due to the ballistic effect occurs in the pore throats. The throats occupy a small proportion of the flow volume but have a limiting effect on the overall flow velocity, which leads to the first reason for the failure of the bundle of tubes models, they do not consider the converging diverging nature of the pore space experienced at these bottlenecks, and ultimately overestimate the absolute permeability. From Figure 12a and 12b we can see that the flow in the lower porosity reconstruction is limited to only a handful of flow pathways, and in Figure 12c and 12d we again observe the largest increase in velocity due 15

ACCEPTED MANUSCRIPT

SC

RI

PT

to the ballistic effect occurs in the throats. The limited connectivity of the pore space severely limits the overall flow velocity. This limitation leads to the second reason for the failure of the bundle of tubes models, they do not consider the connectivity of the pore geometry, which will also lead to overestimation of absolute permeability. Moreover, it is for these same reasons that the bundle of tubes models underestimate the increase in for reconstruction 2 with decreasing pressure. For MVB flows, the pore throats – which occupy a small fraction of the pore space – experience the greatest enhancement in flow due to the ballistic effect (Figure 11d), which is exacerbated by the small number of flow paths. For the more homogenous flow of reconstruction 1 the Civan et al. (2010) bundle of tubes model informed by measured using the full morphology flooding accurately captures , while the other bundle of tubes models overestimate at all pressures simulated due to the reasons outlined in subsection 2.1.

NU

4. Discussion

D

MA

The novelty of the LEV-LBM introduced here is the definition of an effective-mean-free-path correction function, for which the arithmetic mean of a heuristic one-dimensional wall function is employed. Within the framework outlined it is possible to use other wall function formulations, of which many have been conceived (Zhang et al. 2012a), and it may be possible to extend this method of mesoscopic modelling to MVB flows with greater ballistic dominance than what was investigated here.

AC CE P

TE

The only input required for simulation of MVB flow using the LEV-LBM is the unbound mean free path, which is a function of gas composition, specifically the gas density, and the effective hard-sphere diameter of the molecules. The gas density can be calculated using an EOS for the pressure and temperature of interest. The effective hard-sphere diameter of gas mixtures can be estimated as the weighted average of the hard-sphere diameters of the individual gas components. Therefore, the LEVLBM can be used to simulate MVB flow for any gas composition. The bundle of tubes models overestimate the apparent permeability of complex pore spaces and underestimate the increase in with decreasing pore pressure. This overestimation of and underestimation of suggests that the apparent permeability of complex nanoporous media cannot be predicted from the pore size distribution alone. A large amount of effort has gone into quantifying the pore size distributions of shales and mudrocks, with less emphasis on quantifying the three-dimensional connectivity and shape of pores. Without a better understanding of this connectivity it is difficult to formulate predictive models of permeability from pore system characterization. As presented here, the LEV-LBM does not consider the effect of surface diffusion, which can be a dominant form of transport in organic matter nanopores, specifically those where a monolayer of adsorbed molecules occupies a significant fraction of the pore space. Recently, Ren et al. (2015) described a method to incorporate the surface diffusion effect on the slip velocity in LB model simulations, while separately tallying the mass flux for the surface diffusion, . Briefly, adsorption of methane to organic matter follows a Langmuir isotherm (Zhang et al., 2012b). The surface concentration, , measured in adsorbed gas mass per solid volume, can be determined from experimental measurements and described as a function of pore pressure in accordance with monolayer adsorbates theory (Ren et al. 2015),

16

ACCEPTED MANUSCRIPT (33)

RI

PT

where is the solid density, is the gas volume per mole at standard pressure and temperature, the Langmuir volume (maximum adsorption), and is the Langmuir pressure (pressure of half monolayer coverage). The mass flux of surface diffusion is described as,

is

(34)

(35)

MA

NU

SC

where is the density of the adsorbed gas, is the mean velocity of the adsorbed gas, is the surface diffusion constant, and is parallel to the direction of the concentration gradient. With experimental or numerical measurements of adsorption behavior, and can be determined and the previously static walls described by equation (23) can now be parameterized with a wall velocity, ,

AC CE P

TE

D

Using this method to model flow in a two-dimensional channel, Ren et al. (2015) found that channels smaller than 20 nm had a significant contribution to total flux from surface diffusion, with surface diffusion contributing >10% of the total flux for channels < 5 nm in width. The contribution to total mass flux from surface diffusion is dependent on , in Ren et al. (2015) a value of 6.82 x 10-4 was used. In Wu et al. (2015) the for organic matter derived for high pressure conditions was found to be on the order -5 of 10 , an order of magnitude smaller than that used by Ren et al. (2015). Therefore, the surface diffusion contribution to total flux may be much smaller. The inclusion of surface diffusion into the LEV-LBM would lead to a model with a more complete description of mass transport in organic matter, but it is presently not clear if the lack of its inclusion in a geometry more complex than considered by Ren et al. (2015) would appreciably alter the total mass flux. Surface diffusion only occurs along walls that run parallel to the pressure gradient, therefore, when a pore wall turns away from being parallel to the pressure gradient, surface diffusion will cease locally. This suggests that pore space complexity could result in the retardation of surface diffusion throughout the pore space. Much like the methods for describing MVB flows in simple geometries were extended for arbitrary geometries by the LEV-LBM, the same would be required for the inclusion of the Ren et al. (2015) method into the LEV-LBM, which will be a part of our future work. However this does not change the conclusions outlined as follows. 5. Conclusions Here we have introduced a local-effective-viscosity multiple-relaxation-time lattice Boltzmann model (LEV-LBM) capable of predicting mixed viscous-ballistic (MVB) flows (slip-flow regime and earlytransition-flow regime) in arbitrary pore geometries. The model is relevant for supercritical gas flow in nanoporous media. We have validated the model against the mean velocity and velocity profiles of Beskok and Karniadakis (1999) tube model fitted to the linearized Boltzmann solutions of Loyolka and Hamoodi (1990). We have also provided a critical review of some common tube models for MVB flows, including the dusty-gas model, the first-order slip correction of Brown et al. 1946, and the Javadpour (2009) tube model. The dusty-gas model generally overestimates flow in the slip- and early-transitionflow regimes, which is attributed to its assumption that Knudsen diffusion occurs throughout the tube regardless of the thickness of the Knudsen layer. The first-order slip correction of Brown et al. (1946) generally agrees with the Beskok and Karniadakis (1999) tube model in the slip-flow regime, but 17

ACCEPTED MANUSCRIPT underestimates flow for higher Knudsen numbers. The Javadpour (2009) tube model overestimates flow in both the slip- and early-transition-flow regimes.

SC

RI

PT

We have also simulated MVB flow using the LEV-LBM in finite length tubes and quantified the effect of the length to radius ratio on the tube apparent permeability. Finite length tubes are shown to have lower apparent permeability than infinite tubes due to the larger mean free path, and subsequently higher effective viscosity, of the gas at the tube ends. This decrease in apparent permeability is shown to be significant for MVB flow in the slip-flow regime for length to radius ratios less than 2, with much larger decreases for length to radius ratios up to 160 for flows in the early transition flow regime. The length of pore throats in pore-network models of MVB flow is non-trivial.

AC CE P

TE

D

MA

NU

To test the utility of bundle of tubes models for prediction of apparent permeability of supercritical gas flow in complex nanoporous media, we created two digital reconstructions of nanoporous organic matter and simulated pure methane flow at 400 K for a range of pore pressures using the LEV-LBM, and compared the simulation results to the predictions of the bundle of tubes models. To inform the bundle of tubes models, the mean radius of the reconstructions was measured using a maximum-inscribed-sphere algorithm and a full morphology model all-boundary flooding. The difference in predicted apparent permeability between the LEV-LBM and the bundle of tubes models are large, with the bundle of tubes models greatly overestimating the apparent permeability by factors between 1.52 and 4.95 for the high porosity reconstruction, and 33.2 and 153 for the low porosity reconstruction, while underestimating the increase in apparent permeability over absolute permeability with decreasing pore pressures for the low porosity reconstruction. The overestimation of absolute permeability that leads to the overestimation of apparent permeability by the bundle of tubes models, and the underestimation of the increase of apparent permeability over absolute permeability with decreasing pressure for the low porosity reconstruction in comparison to LEV-LBM simulation results can be attributed to the non-trivial effect the shape of the pores and the connectivity of the pore space has on gas flow. On the other hand, the overestimation of the increase of apparent permeability over absolute permeability with decreasing pressure for the high porosity reconstruction by the dusty-gas and Naraghi and Javadpour 2015 bundle of tubes model in comparison to LEV-LBM simulation results can be attributed to the overestimation of mass flux by the underlying single tube models in the slip- and early transition-flow regimes. These results indicate that the pore size distribution alone of organic matter pore spaces is an insufficient predictor of apparent permeability. The key to predicting MVB flow using the mesoscopic approach presented here hinges on the accuracy of the effective mean free path correction function, . The formulation for could be further tested and tuned by comparison to molecular dynamics simulations, which will be a part of our future work. Experimental validation will require the development of new MVB flow datasets measured under high pressure conditions in well characterized nanoporous media. However, as it currently stands, the LEVLBM can act as a powerful tool to elucidate the relationship between pore space geometry and apparent permeability for supercritical gas flow in complex nanoporous media.

Acknowledgements 18

ACCEPTED MANUSCRIPT This research was supported by Shell-UT Unconventional Research (SUTUR) collaborative program at The University of Texas at Austin. The authors acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing HPC resources. Publication approved by the Director, Bureau of Economic Geology.

PT

References

RI

Agrawal, A., 2012. A comprehensive review on gas flow in microchannels. Int. J. of Micro-Nano Scale Transport 2(1), 1-40.

SC

Agrawal, A., Prabhu, S.V., 2008. Survey on measurement of tangential momentum accomadation coefficient. J. Vac. Sci. Techol. A. 26(4), 634-635.

NU

Allan, A., Mavko, G., 2013. The effect of adsorption and Knudsen diffusion on the steady-state permeability of microporous rocks. Geophys. 78(2), D75-D83.

MA

Ansumali, S., Karlin, I.V., 2002. Kinetic boundary conditions in the lattice Boltzmann method. Phys. Rev. E 66, 026311.

D

Arkilic, E.B., Schmidt, M.A., Breuer, K.S., 1997. Gaseous slip flow in long microchannels. J. Microelectromech. Syst. 6(2), 167-178.

TE

Beskok, A., Karniadakis, G.E., 1999. A model for flows in channels, pipes, and ducts at micro and nano scales. Microscale Therm. Eng. 3, 43-77.

AC CE P

Brown, G.P., DiNardo, A., Cheng, G.K., Sherwood, T.K., 1946. The flow of gases in pipes at low pressures. J. Appl. Phys. 17, 802-813. Cao, B.Y., Chen, M., Guo, Z.Y., 2006. Effect of surface roughness on gas flow in microchannels by molecular dynamics simulation. Int. J. Eng. Sci. 44, 927-937. Chai, Z., Shi, B., Guo, Z., Lu, J., 2010. Gas flow through square arrays of circular cylinders with Klinkenberger effect: A lattice Boltzmann Study. Commun. Comput. Phys. 8(5), 1052-1073. Chalmers, G.R., Bustin, R.M, Power, I.M., 2012. Characterization of gas shale pore systems by porosimetry, pycnometry, surface area, and field emission scanning electron microscopy/transmission electron microscopy image analyses: Examples from the Barnett, Woodford, Haynesville, Marcellus, and Doig units. AAPG Bull. 96(6), 1099-1119. Chen, L., Zhang, L., Kang, Q., Viswanathan, H.S., Yao, J., Tao, W., 2015. Nanoscale simulation of shale transport properties using the lattice Boltzmann method: permeability and diffusivity. Sci. Rep. 5, 8089. Civan, F., 2010. Effective correlation of apparent gas permeability in tight porous media. Transport Porous Med. 82, 375-384. Colin, S., 2005. Rarefaction and compressibility effects on steady and transient gas flows in microchannels. Microfluid. Nanofluid. 1, 268-279. Darabi, H., Ettehad, A., Javadpour, F., Sepehrnoori, K., 2012. Gas flow in ultra-tight shale strata. J. Fluid Mech. 710, 641-658. 19

ACCEPTED MANUSCRIPT

d‟Humieres, D., Ginzburg, I., Krayczyk, M., Lallemand, P., Luo, L.S., 2002. Multiple-relaxation-time lattice Boltzmann models in three dimensions. Philos. T. Roy. Soc. A 360, 437-451.

PT

Falk, K., Coasne, B., Pellenq, R., Ulm, F.J., Bocquet, L., 2015. Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nat. Comm. 6, 6949.

RI

Firouzi, M., Wilcox, J., 2012. Molecular modeling of carbon dioxide transport and storage in porous carbon-based materials. Micropor. Mesopor. Mat. 158, 195-203.

SC

Guo, Z.L., Shi, B.C., Zheng, C.G., 2007. An extended Navier-Stokes formulation for gas flows in the Knudsen layer near a wall. Earth Planet. Sc. Lett. 80, 24001.

NU

Guo, Z., Zhao, T.S., Shi, Y., 2006. Physical symmetry, spatial accuracy, and relaxation time of the lattice Boltzmann equation for microgas flows. J. Appl. Phys. 99, 074903.

MA

Guo, Z., Zheng, C., 2008. Analysis of lattice Boltzmann equation for microscale gas flows: Relaxation times, boundary conditions and the Knudsen layer. Int. J. Comput. Fluid Dyn. 22(7), 465-473.

D

Guo, Z., Zheng, C., Shi, B., 2008. Lattice Boltzmann equation with multiple effective relaxation times for gaseous microscale flow. Phys. Rev. E 77, 036707.

TE

Hadjiconstantinou, N.G., 2006. The limits of Navier-Stokes theory and kinetic extensions for describing small-scale gaseous hydrodynamics. Phys. Fluids 18, 111301.

AC CE P

Heath, J.E., Dewers, T.A., McPherson, B.J.O.L., Petrusak, R., Chidsey, T.C., Rinehart Jr, A.J., Mozley, P.S., 2011. Pore networks in continental and marine mudstones: Characteristics and controls on sealing behavior. Geosphere 7, 429-454. Hecht, M., Harting, J., 2010. Implementation of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations. J. Stat. Mech.: Theory Exp. 1, P01018. Islam, A., Patzek, T., 2014. Slip in natural gas flow through nanoporous shale reservoirs. J. Unconvent. Oil Gas Resour. 7, 49-54. Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). J. Can. Petrol. Technol. 48(8), 16-21. Jiao, K., Yao, S., Liu, C., Gao, Y., Wu, H., Li, M., Tang, Z., 2014. The characterization and quantitative analysis of nanopores in unconventional gas reservoirs utilizing FESEM-FIB and image processing: An example from the lower Silurian Longmaxi Shale, upper Yangtze region, China. Int. J. Coal Geol. 128129, 1-11. Kazemi, M., Takbiri-Borujeni, A., 2015. An analytical model for shale gas permeability. Int. J. Coal Geol. 146, 188-197. Kim, S.H., Pitsch, H., Boyd, I.D., 2008. Slip velocity and Knudsen layer in the lattice Boltzmann method for microscale flows. Phys. Rev. E 77, 026704.

20

ACCEPTED MANUSCRIPT Klaver, J., Desbois, G., Littke, R., Urai, J.L. 2015. BIB-SEM characterization of pore space morphology and distribution in postmature to overmature samples from the Haynesville and Bossier Shales. Mar. Pet. Geol. 59, 451-466.

PT

Klinkenberg, L.J. 1941. The permeability of porous media to liquids and gases. API Drill. Prod. Prac, 200-213.

RI

Li, Q., He, Y.L., Tang, G.H., Tao, W.Q., 2011. Lattice Boltzmann modeling of microchannel flows in the transition flow regime. Microfluid. Nanofluid. 10, 607-618.

SC

Lim, C.Y., Shu, C., Niu, X.D., Chew, Y.T., 2002. Application of lattice Boltzmann method to simulate microchannel flows. Phys. Fluids 14, 2299.

NU

Liu, X., Guo, Z., 2013. A lattice Boltzmann study of gas flows in a long micro-channel. Comput. Math. Appl. 65, 186-193.

MA

Liou, T.M., Lin, C.T., 2014. Study on microchannel flows with a sudden contraction-expansion at a wide range of Knudsen number using lattice Boltzmann method. Microfluid Nanofluid 16, 315-327.

D

Loucks, R.G., Reed, R.M., Ruppel, S.C., Hammes, U., 2012. Spectrum of pore types and networks in mudrocks and a descriptive classification for matrix-related mudrock pores. AAPG Bull. 96(6), 10711098.

TE

Loyalka, S.K., Hamoodi, S.A., 1990. Poiseuille flow of a rarefied gas in a cylindrical tube: Solution of linearized Boltzmann equation. Phys. Fluids 2, 2061-2065.

AC CE P

Ma, J., Sanchez, J.P, Wu, K., Couples, G.D., Jiang, Z., 2014. A pore network model for simulating nonideal gas flow in micro- and nano-porous materials. Fuel 116, 498-508. Marino, L., 2009. Experiments on rarefied gas flows through tubes. Microfluid. Nanofluid. 6, 109-119. Mason, E.A., Malinauskas, A.P., 1983. “Gas transport in porous media: The dusty-gas model” Elsevier Science, Amsterdam, Netherlands. Maxwell, J.C., 1879. “The scientific papers of James Clerk Maxwell” (p. 707). Dover Publications, New York, NY. Mehmani, A., Prodanović, M., 2014. The effect of microporosity on transport properties in porous media. Adv. Water Resour. 63, 104-119. Mehmani, A., Prodanović, M., Javadpour, F., 2013. Multiscale, multiphysics network modeling of shale matrix gas flows. Transport Porous Med. 99, 377-390. Milliken, K.L., Rudnicki, M., Awwiller, D.N., Zhang, T., 2013. Organic matter-hosted pore system, Marcellus Formation (Devonian), Pennsylvania. AAPG Bull. 97(2), 177-200. Naraghi, M.E., Javadpour, F., 2015. A stochastic permeability model for the shale-gas systems. Int. J. Coal Geol. 140, 111-124.

21

ACCEPTED MANUSCRIPT Nie, X., Doolen, G.D., Chen, S., 2002. Lattice-Boltzmann simulations of fluid flows in MEMS. J. Stat. Phys. 107(1/2), 279-289.

PT

Niu, X.D., Hyodo, S.A., Munekata, T., 2007. Kinetic lattice Boltzmann method for microscale gas flows: Issues on boundary condition, relaxation time, and regularization. Phys. Rev. E 76, 036711.

RI

Pan, C., Luo, L.S., Miller, C.T., 2006. An evaluation of lattice Boltzmann schemes for porous medium flow simulation. Comput. Fluids 35, 898-908.

SC

Peng, D.Y., Robinson, D.L., 1976. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 15, 59-64.

NU

Peng, S., Yang, J., Xiao, X., Loucks, R., Ruppel, S.C., Zhang, T., 2015. An integrated method for upscaling pore-network characterization and permeability estimation: Example from the Mississippian Barnett shale. Transp. Porous Med. 109(2), 359-376.

MA

Qian, Y.H., d‟Humieres, D., Lallemand, P., 1992. Lattice BGK models for Navier-Stokes equations. Europhys. Lett. 17(6), 479-484. Rahmanian, M.R., Aguilera, R., Kantzas, A., 2013. A new unified diffusion-viscous-flow model based on pore-level studies of tight gas formations. SPE J. 18(1), 38-49.

TE

D

Ren, J., Guo, P., Guo, Z., Wang, Z., 2015. A lattice Boltzmann model for simulating gas flow in kerogen pores. Transport Porous Med. 106, 285-301.

AC CE P

Roy, S., Raju, R., Chuang, H.F., Cruden, B.A., Meyyappan, M., 2003. Modeling gas flow through microchannels and nanopores. J. Appl. Phys. 93(8), 4870-4879. Sakhaee-Pour, A., Bryant, S.L., 2012. Gas permeability of shale. SPE Reservoir Eval. Eng., 401-409. Sabragaglia, M., Succi, S., 2005. Analytical calculation of slip flow in lattice Boltzmann models with kinetic boundary conditions. Phys. Fluids 17, 093602. Sharipov, F.M., Seleznev, V.D., 1994. Rarefied gas flow through a long tube at any pressure ratio. J. Vac. Sci. Technol. A 12, 2933-2935. Shen, C., Tian, D.B., Xie, C., Fan, J., 2004. Examination of the LBM in simulation of microchannel flow in transitional regime. Microscale Thermophys. Eng. 8, 423-432. Singh, H., Javadpour, F., Ettehadtavakkol, A., Darabi, H., 2014. Nonempirical apparent permeability of shale. SPE Reservoir Eval. Eng. 17(3), 414-424. Stops, D.W., 1970. The mean free path of gas molecules in the transition regime. J. Phys. D 3, 685-696. Succi, S., 2002. Mesoscopic modeling of slip motion at fluid-solid interfaces with heterogeneous catalysis. Phys. Rev. Lett. 89(6), 064502. Suga, K., and Ito, T., 2011.On the multiple-relaxation-time micro-flow lattice Boltzmann method for complex flows. Comput. Model. Eng. Sci. 75(2), 141-172.

22

ACCEPTED MANUSCRIPT Suga, K., Takenaka, S., Ito, T., Kaneda, M., Kinjo, T., Hyodo, S., 2010. Evaluation of a lattice Boltzmann method in a complex nanoflow. Phys. Rev. E 82, 016701.

PT

Tison, S., 2009. Experimental data and theoretical modeling of gas flows through metal capillary leaks. Vacuum 44, 1171-1175.

RI

Verhaeghe, F., Luo, L.S., Blanpain, B., 2009. Lattice Boltzmann modeling of microchannel flow in slip flow regime. J. of Comput. Phys. 228, 147-157.

SC

Vogel, H.-J., Tolke, J., Schulz, V.P., Krafczyk, M., Roth, K., 2005. Comparison of a lattice-Boltzmann model, a full-morphology model, and a pore network model for determining capillary pressure-saturation relationships. Vadose Zone J. 4, 380-388.

NU

Wood, J.M., Sanei, H., Curtis, M.E., Clarkson, C.R., 2015. Solid bitumen as a determinant of reservoir quality in an unconventional tight gas siltstone play. Int. J. Coal Geol. (In Press).

MA

Wu, K., Li, X., Wang, C., Yu, W., Chen, Z., 2015. Model for surface diffusion of adsorbed gas in nanopores of shale gas reservoirs. Ind. Eng. Chem. Res. 54, 3225-3236. Yao, J., Sun, H., Fan, D., Wang, C., Sun, Z., 2013. Numerical simulation of gas transport mechanisms in tight shale gas reservoirs. Petrol. Sci. 10, 528-537.

TE

D

Zhang, W.M., Meng, G., Wei, X., 2012a. A review on slip models for gas microflows. Microfluid. Nanofluid. 13, 845-882.

AC CE P

Zhang, T., Ellis, S.E., Ruppel, S.C., Milliken, K., Yang, R., 2012b. Effect of organic-matter type and thermal maturity on methane adsorption in shale-gas systems. Org. Geochem. 47, 120-131. Zheng, L., Guo, Z., Shi, B., 2012. Microscale boundary conditions of the lattice Boltzmann equation method for simulating microtube flows. Phys. Rev. E 86, 016712.

23

ACCEPTED MANUSCRIPT Figure Captions Figure 1: Plot of Knudsen numbers as a function of pressure for tubes of varying radius for methane at 400 K. Also included on left are the Knudsen flow regimes.

SC

RI

PT

Figure 2: Comparison of tube models and linearized Boltzmann results of Loyalka and Hamoodi (1990) for pure methane flow at 400 K through a tube with a 5 nm radius. DGM: the dusty-gas tube model, Brown: Brown et al. (1946) tube model, Javadpour: Javadpour (2009) tube model, B-K: Beskok and Karniadakis (1999) tube model, and Linearized Boltzmann: linearized Boltzmann solutions of Loyalka and Hamoodi (1990).

D

MA

NU

Figure 3: Wall function for determination of effective mean free path. (a) Mean free path in an unbound geometry, where the mean free path is constant. (b) Mean free path in a bound geometry, where the mean free path is a function of distance from the wall, with the arrow lengths and colors indicating the relative variance in the length of the mean free path. Smaller mean free paths do not indicate compression of the gas, the density of the gas is constant barring adsorptive effects. (c) One-dimensional wall function ( ) as a function of distance from the wall ( ) for four different choices of the free parameters. (d) Threedimensional 18-direction wall distance measurement schematic with a 2D illustration of the measurement of wall distance for a single lattice node, wall distance is measured in three dimensions for flow simulation.

TE

Figure 4: Images of digital reconstruction 1 (left) and reconstruction 2 (right) of nanoporous organic matter, with connected pore space shown as solid magenta and the solid as transparent.

AC CE P

Figure 5: Images of the maximum inscribed sphere (MIS) algorithm applied to reconstruction 1 (left) and reconstruction 2 (right). Figure 6: Example images of the full morphology model flooding (simulated mercury injection porosimetry) for reconstruction 1 (top) and reconstruction 2 (bottom). The pore spaces are shown as translucent magenta with the drainage front shown in blue, the radii labels indicate the capillary entry pressure radius of the flood for the image below. Figure 7: Pore radius distribution for reconstruction 1 (R1) and 2 (R2) using the maximum inscribed sphere measurement (MIS) and the full morphology model flooding (Flood). Figure 8: Comparison of the LEV-LBM to the Beskok and Karniadakis (1999) tube model. (a) Error (in percent) for the mean velocity for four different choices of the free parameters and the multi-bounce onedimensional wall function ( ) as a function of . (b) Error (in percent) of the mean velocity for four different as a function of lattice resolution using the multi-bounce one-dimensional wall function, including images of the tube at different resolutions. (c-f) Velocity profiles for the Beskok and Karniadakis (1999) tube model and the LEV-LBM using both the effective mean free path correction function and the multi-bounce effective mean free path correction function, and, for comparison, the Suga and Ito (2011) global effective viscosity model for Knudsen numbers (c) 0.02, (d) 0.2, (e) 0.5, and (f) 1. Figure 9: MVB flow in finite length tubes. (a) Illustration of the end effect, showing a cross-sectional image of and the increase in local effective viscosity at the inlet and outlet of the tube, also shown are individual measurements of , with dashed lines indicating infinite measurements resulting in 24

ACCEPTED MANUSCRIPT

PT

(note: measurements are made in three-dimensions, however only those parallel to the plane of the image are shown for illustrative purposes). (b) Pressure bound simulation setup illustration, with layer indicating a single layer of lattice nodes. (c) The dependence of on length to radius ratio, , normalized to the apparent permeability of infinite tubes, .

SC

RI

Figure 10: and as a function of pressure for pure methane flow at 400 K, with reconstruction 1 on the left and reconstruction 2 on the right. CIV: Civan (2010) bundle of tubes model, DGM: dusty-gas bundle of tubes model, and NJ: Naraghi and Javadpour (2015) bundle of tubes model, with MIS: maximum inscribed sphere measurement of mean radius, and Flood: full morphology model flood measurement of mean radius.

NU

Figure 11: Velocity flow fields of reconstruction 1 for viscous no-slip viscous flow (standard LBM) and pure methane flow at 400 K and 4.06 MPa (LEV-LBM), showing the entirety of the simulation domain for (a) standard LBM and (b) LEV-LBM, as well as cross-sections for the (c) standard LBM and (d) LEV-LBM. Note the increase in flow at pore throats for the MVB flow simulated by the LEV-LBM.

TE

D

MA

Figure 12: Velocity flow fields of reconstruction 2 for Navier-Stokes no-slip viscous flow (standard LBM) and pure methane flow at 400 K and 4.06 MPa (LEV-LBM) simulation, showing the entirety of the simulation domain for (a) standard LBM and (b) LEV-LBM, as well as cross-sections for the (c) standard LBM and (d) LEV-LBM. Note the increase in flow at pore throats for the MVB flow simulated by the LEV-LBM and limited number of flow pathways. Tables

AC CE P

Table 1: Summary of mean radius determined using maximum inscribed spheres (MIS) and a full morphology flooding (Flood), bundle of tubes model absolute permeability, and standard LBM determined absolute permeability.

Reconstruction 1

0.5832

MIS [nm] 5.084

Reconstruction 2

0.1963

4.821

MIS [nd] 929.7

Flood [nm] 4.223

Flood [nd] 641.5

LBM [nd] 426.3

281.3

2.585

80.86

0.5826

25

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 1

26

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 2

27

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 3

28

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 4

29

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 5

30

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 6

31

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 7

32

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 8

33

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 9

34

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 10

35

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 11

36

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Figure 12

37

ACCEPTED MANUSCRIPT

AC CE P

TE

D

MA

NU

SC

RI

PT

Graphical abstract

38

ACCEPTED MANUSCRIPT Direct simulation of supercritical gas flow in complex nanoporous media and prediction of apparent permeability Authors:

Christopher J. Landry*,a, Maša Prodanovićb, Peter Eichhublc

a

PT

Center for Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX 78712, USA b

RI

Department of Petroleum and Geosystems Engineering, The University of Texas at Austin, Austin, TX 78712, USA c

SC

Bureau of Economic Geology, Jackson School of Geosciences, The University of Texas at Austin, Austin, TX 78757, USA

TE

D

MA

Developed local effective viscosity lattice Boltzmann model with non-zero slip Apparent permeability of finite length tubes significantly decreases for 0.02
AC CE P

    

NU

Highlights:

39