- Email: [email protected]
Review of multi-scale and multi-physical simulation technologies for shale and tight gas reservoirs
Accepted Manuscript Review of multi-scale and multi-physical simulation technologies for shale and tight gas reservoirs Lei Wang, Shihao Wang, Ronglei Zhang, Cong Wang, Yi Xiong, Xishen Zheng, Shangru Li, Kai Jin, Zhenhua Rui PII:
S1875-5100(16)30887-3
DOI:
10.1016/j.jngse.2016.11.051
Reference:
JNGSE 1961
To appear in:
Journal of Natural Gas Science and Engineering
Received Date: 9 August 2016 Revised Date:
15 November 2016
Accepted Date: 24 November 2016
Please cite this article as: Wang, L., Wang, S., Zhang, R., Wang, C., Xiong, Y., Zheng, X., Li, S., Jin, K., Rui, Z., Review of multi-scale and multi-physical simulation technologies for shale and tight gas reservoirs, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.11.051. 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
Review of Multi-scale and Multi-physical Simulation Technologies for Shale and Tight Gas Reservoirs
Abstract
RI PT
Lei Wanga, Shihao Wanga, Ronglei Zhanga*, Cong Wanga, Yi Xionga, Xishen Zhengb, Shangru Lib, Kai Jinb, Zhenhua Ruic a Energy Modeling Group, Petroleum Engineering Department, Colorado School of Mines, USA b Jiangxi Unconventional Gas Institute, China c Independent Project Analysis, Inc., USA *Corresponding author: Ronglei Zhang, Email: [email protected]
Keyword:
TE D
M AN U
SC
This paper provides a comprehensive literature review on the simulation techniques being developed in recent years for describing unique flow behaviors in shale and tight gas reservoirs. The advances in modeling gas flow and transport mechanisms during the primary and enhanced gas recovery processes are reviewed in detail. The capabilities of reservoir simulation tools are discussed in terms of mathematical treatment (finite difference, finite element, explicit/implicit scheme, sequentially and fully coupled schemes), fluid flow characteristics (Darcy and non-Darcy flow, desorption, Klinkenberg effect and gas slip flow, transitional flow, Knudsen diffusion), reservoir rock properties (pore size distribution, fractures, geomechanics), multidisciplinary coupling scheme (THM, THC, THMC), modeling scale, and computational efficiency. For pore-scale modeling of gas flow and transport in unconventional reservoir rocks, the numerical methods for generating pore network models and the procedure for constructing real 3D digital rocks are first explained. Then, the network modeling methods used for simulating slip and transitional flows in pore networks are introduced and compared. After that, existing lattice-Boltzmann models developed for slip flow simulation are illustrated. Along with the explanation, pros and cons of the models for pore-scale modeling have been identified. Overall, dwelling on the concerns and challenges in shale and tight gas reservoir simulation, current status, progress and bottlenecks of numerical simulation technologies are discussed and perspectives on future development of unconventional gas reservoir simulators are proposed.
AC C
Highlights
EP
Shale and tight gas reservoir, reservoir simulation, pore-scale simulation, adsorption and desorption, Klinkenberg effect and slip flow, transitional flow, Knudsen diffusion, non-linear and non-Darcy flow, geomechanics, THMC
1. Fluid flow and transport mechanisms along with geomechanics in shale and tight gas reservoirs are recapitulated. 2. Advances in macroscopic THMC simulation techniques for shale and tight gas reservoirs are reviewed. 3. Pore-scale modeling methods for fluid flow and transport in unconventional reservoir rocks are illustrated and compared. 4. Challenges from unconventional reservoirs and perspectives on developing new simulation techniques are discussed.
1
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
1. Introduction Unconventional gas reservoirs, such as shale and tight gas reservoirs, are rich resources of natural gas. In unconventional formations, huge amount of gas is trapped in the location where it is generated, without the primary migration process. Featured by extremely low permeability, shale and tight gas reservoirs are different from each other, even with similar porosity, pore pressure, and in-situ stress conditions. As shown in Figure 1 (Wang et al. 2009), permeability values of organic-rich shales are in the range of sub-nanodarcy to tens of microdarcy, and permeability of deep organic-lean mudrocks can be even less. Some correlations can also be observed between the permeability values and the porosity of shales.
Figure 1. The relationships of porosity and permeability for shale gas reservoirs (FWB-Fort Worth Basin, WV #6 well-Marcellus Shale well in West Virginia) in North America obtained from core analysis (Wang et al. 2009).
AC C
EP
The diameter of pore throat in tight gas formations is usually between 0.01µm and 1µm (Holditch 2006, Nelson 2009), while the pore throat diameter of shale formation varies from tens of nanometers to as low as a few nanometers (Nelson 2009). Recent studies even revealed sub-nanometer pores inside shale formations that are beyond the scope of scanning electron microscope (SEM) based on CO2 and N2 adsorption experiments (Clarkson et al. 2013) and atom-scale reconstruction technique (Falk et al. 2015). Due to these extremely small pores, the hydrocarbon transport mechanisms in unconventional reservoirs turn out to be very different from those in conventional reservoirs. Fluid flow behavior in shale and tight gas reservoirs is featured by multiscale single-phase (gas) and/or multiphase fluid (gas, gas condensate and/or brine) flow and transport in ultra-low permeability, highly heterogeneous porous/fractured, and stress-sensitive rocks. Compared with flow in conventional reservoirs, fluid flows in unconventional reservoirs with extremely low permeability are highly non-linear and complex processes, including non-linear adsorption/desorption, non-Darcy flow in the entire range from high flow rate to low flow rate, strong rock-fluid interaction, and rock/organic matter deformation within nano-pores or micro-fractures. Therefore, quantifying flow in unconventional gas reservoirs has become a challenging issue for conventional REV (Representative Elementary Volume)-based Darcy’s law approaches, which may not be suitable and thus must be improved in order to handle these realistic problems in shale and tight gas reservoirs. 2
ACCEPTED MANUSCRIPT
SC
RI PT
2. Gas Flow and Transport Mechanisms Due to the existence of organic matter with strong adsorption capacity in shale gas reservoirs, natural gas components experience desorption from its original adsorbed state as reservoir pressure depletes. In addition, flow in micro- and nano-pores in unconventional reservoirs falls in different regimes from those in conventional reservoirs. Besides, decrease of pore pressure would induce compaction of formations and reduction of the effective pore sizes, resulting in permeability decrease. Also organic matter may shrink due to desorption of adsorbed gas, resulting in dimensional changes of the pores they are residing in. Since shale and tight gas reservoirs are generally completed by hydraulic fracturing, disparities between flow behavior in matrix and fracture networks should be considered as well. All these flow and transport mechanisms are described in the following sections. 2.1 Gas Adsorption and Desorption
AC C
EP
TE D
M AN U
Natural gas presents as free gas and adsorbed gas in shale reservoirs. Significant amount of gas can be adsorbed onto the surface of the organic matter and clays under reservoir conditions (Mengal and Wattenbarger, 2011, EIA, 2011a). In shale matrix, methane molecules are mainly adsorbed onto the organic matter, i.e. kerogen, which is generally measured in terms of total organic carbon (TOC). Adsorbed gas can provide significant portion of gas in place (GIP) as well as affect ultimate recovery rate. In Barnett Shale, when pressure is below 2,000 psi, adsorbed gas accounts for 50 to 80 % of GIP. As pressure goes above 2,000 psi, although the adsorbed gas content tends to level out, it is still significant in percentage, about 30 to 50 % (Mengal and Wattenbarger, 2011). Figure 2 is a case study of Barnett Shale showing typical percentage of adsorbed gas to total gas, and its linear dependence on TOC content (Wang et al. 2009). Gas adsorption capacity also changes with in-situ reservoir pressure. As pressure declines under continuous reservoir depletion, more adsorbed gas is dissociated from adsorbed state to free gas, leading to sustainable gas flow.
Figure 2. Adsorbed and total gas content versus TOC of Barnett Shale (Wang et al. 2009).
2.2 Gas Rarefaction Effects and Flow Regimes Gas rarefaction effects in unconventional reservoirs, including velocity slippage and gas molecular diffusion, have been demonstrated by core flooding experiments (Luffel et al. 1993, Civan et al. 2011, Tinni et al. 2012, Zhao et al., 2015a and 2015b). These two transport phenomena are essentially induced by the gas molecule and 3
ACCEPTED MANUSCRIPT
RI PT
pore wall interactions. In gas sate, molecules keep moving in random directions with thermal velocity. The average distance that one gas molecule can move freely before it collides with any other molecules is defined as the mean free path λ of the gas, which can be calculated as (Bird 2003) k BT (1) λ= 2π Pδ 2 In the above formulation for a single component ideal gas, denotes the Boltzmann constant, T denotes the temperature, P denotes the pressure of the system, δ denotes the collision diameter of the gas. For methane, δ is about 0.40 nm. When the pore diameter d becomes narrow enough to be comparable with the mean free path of the gas molecules in it, the gas molecules will frequently collide with the walls of the flow channel and the molecule-wall interaction cannot be ignored. The gas inside such narrow flow channels becomes rarefied and possesses certain unique transport characteristics. The gas rarefication is quantified by Knudsen number as (Bird 2003)
λ
SC
Kn =
d
(2)
M AN U
where d is the characteristic length, e.g. the diameter of cylindrical pores. The average pore size can either be directly measured from experiments or be estimated from permeability
k ∞ , tortuosity τ, and porosity ϕ of the rock. For example, Civan (2010) used the following correlation to estimate the mean pore radius of shale.
r = 2 2τ
k∞
φ
(3)
AC C
EP
TE D
When a hard-sphere gas molecule collides with the solid wall, its reflection will be either specular or diffuse, depending on the roughness of the wall as well as the molecule type. On one hand, after specular reflection, as shown by the blue arrow in Figure 3, the molecule maintains its previous horizontal velocity component, changing only its moving direction. On the other hand, after diffuse reflection, the molecule will reflect diffusely into a random direction, changing its horizontal and vertical velocity components, as shown by the red arrows in Figure 3. The maintenance of horizontal velocity component results in the ‘slippage’ of gas molecules in the near-wall region of the pores. Meanwhile, the reflected molecules themselves contribute to the diffusion effect in pores with molecular density differences. Therefore, the apparent permeability enhancement is a result of the combined velocity slippage and diffusion effect, both of which come from the molecule-wall interaction. The ratio of the diffuse reflection to the total amount of wall collision is defined as tangential momentum accommodation coefficient (TMAC) σ, ranging from 0 to 1. The velocity slippage as well as the diffusion effect results in the permeability enhancement that the apparent gas permeability of shale and tight gas rocks is larger than their absolute permeability . The relationship between and can be quantified as
ka = f ( Kn ) k∞
(4)
Where f(Kn) is a permeability multiplier, depending on Knudsen number and gas molecule type (Zhang, W.-M. et al. 2012).
4
RI PT
ACCEPTED MANUSCRIPT
Figure 3. Conceptual sketch showing specular reflection and diffuse reflection.
M AN U
SC
To describe gas flow and transport in micro- and nanoscale geometries, the governing equation is chosen upon flow regimes. Table 1 shows the gas flow regimes and corresponding governing equations along with boundary conditions. When is between 0.001 and 0.1, gas flow is in the regime of slip flow, slip boundary condition should be incorporated into Navier-Stokes equation or lattice-Boltzmann method (LBM) to take into account the slippage on the gas solid interface. When is in a higher range of 0.1 and 10, gas flow enters the transitional regime, where neither Navier-Stokes equation nor lattice-Boltzmann model is applicable any more. Then, Burnett equation based on higher order moments of Boltzmann equation should be solved or numerical method of direct simulation Monte Carlo (DSMC) should be used to represent the fluid flow behavior. As goes beyond 10, gas stream is considered as free molecules, and Molecular Dynamics (MD) must be adopted to capture the physics controlling the gas flow.
Continuum model Kn Flow regimes Note
Boltzmann equation
Collisionless Boltzmann equation
Lattice-Boltzmann method
DSMC
Molecular Dynamics
Navier-Stokes equation
Burnett equation
×
(0, 0.001) No slip
(0.001, 0.1) Slip
EP
Particle or molecule model
TE D
Table 1. Knudsen number and flow regimes with applicable mathematical models
(0.1, 10) Transitional Adapted from Bird, 2003.
>10 Free molecular
AC C
Many petrophysical measurements have demonstrated that pore throat radius of shale reservoirs falls in the nanometer range. For instance, Figure 4 shows the pore throat radius distribution of a typical shale formation measured by mercury injection into both crushed and plug core samples (Honarpour et al. 2012). It can be seen that the measurements of the pore size distribution are relatively consistent with each other in that the majority of the pores lie in the range of 1 to 100 nanometers by volumetric fraction. Considering that the major component of unconventional natural gas is methane, more than 80% in general (Schettler Jr. and Parmely 1989), the mean free path of methane gas molecules under reservoir conditions can be estimated from the equation above for ideal gas. Setting the unconventional gas reservoir at a typical temperature of 353.15K or 176ºF (Gong et al. 2013), when the reservoir pressure is 1,000 psi, in Figure 4, the pores bounded by the blue dashed lines, from ~9 to ~900 nanometers, host the slip flow regime defined by 0.001< <0.1. Given the reservoir pressures of 2,000 psi and 3,500 psi, the slip flow regime shifts to smaller pores, as denoted by the pore size ranges between the green dotted lines and red lines, respectively. To the left of these pore size ranges, 5
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
transitional flow occurs, while to the right, it is the familiar viscous flow regime that are accounted for by NavierStokes equation, from which Darcy’s law is derived for fully developed Stokes flow. It is worth noting that under the reservoir temperature and pressure conditions mentioned above, methane is in the supercritical state ( =190.6K, =668.6psi). Experiments for fluid flow through nano-pores demonstrated that diffusive flow mechanisms can be negligible when the fluid is in the supercritical state (Patil et al. 2006).
Figure 4. Measured pore throat radius distribution of a shale reservoir and corresponding flow regimes under different reservoir conditions.
AC C
EP
However, it should be noticed that flow regime is merely a classification, whereas the gas flow behavior is continuously changing with respect to Kn. Moreover, the pore structure in unconventional rocks is essentially irregular, pore sizes are continuous, and thus more than one flow regimes may exist simultaneously in the rock matrix under certain reservoir conditions. Therefore, it is of great importance to develop a unified model that can describe more than one flow regimes in order to accurately capture the commingled gas flow behavior in unconventional reservoir rocks. Figure 5 shows a conceptual model of gas flow along the x direction through a channel, in which a horizontal imaginary plain H and a vertical imaginary plain V are randomly selected for molecule movement analysis. Molecules of Type ① and ③ impinge into H after hitting the upper and lower wall, respectively. Molecules of Type ② and ④ impinge into H without any impact of the wall. The difference between the horizontal flux momentum carried by Type ① and ② molecules and that carried by Type ③ and ④ molecules generates shear stress on plain H. Because of the wall impact on Type ① and ③ molecules, their resultant shear stress along plain H is larger than the shear stress on an imaginary plain that is located inside an infinitely wide flow channel. This mechanism results in the velocity slippage in small channels. The wall impact is a function of the location along the z direction between the two walls. The region that is closer to the wall has high wall impact and the region that is closer to the middle of the flow channel has low or no wall impact. Similarly, the unequal
6
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
molecules passing through the plain V reduces the number density difference across it, resulting in the diffusion effect. The magnitude of the wall impact on diffusion also depends on the pore wall proximity.
Figure 5. Conceptual model of reflected molecules and bulk molecules.
AC C
EP
TE D
Based on the above analysis and for engineering purpose, the flux inside the micro-channel can be approximately divided into two regions: the Knudsen layer and the bulk flow, as shown in Figure 6. The flow inside regions that are close to the upper and lower walls are assumed to be primarily affected by the moleculewall interaction. The regions are therefore called the Knudsen layers. The thickness of the Knudsen layer is usually taken as 2/3 of the mean free path, the derivation of which can be found in (Present, 1958). The flow beyond the Knudsen layer is assumed to be not affected by the wall and is called the bulk flow, which can still be modeled by the Navier-Stokes equation. The transition from the Knudsen layer to bulk layer is not as distinct as shown in Figure 6. Clearly, more advanced techniques are required for accurate simulation of the whole flow field, as discussed in the following sections.
Figure 6. Conceptual velocity profile within a micro-channel. The yellow curve represents the bulk flow part where the wall barely impacts the flow. The ‘straight lines’ inside the Knudsen layer indicate the average velocity of the Knudsen layer. Real flux velocity of gas molecules inside the Knudsen layer is not constant.
7
ACCEPTED MANUSCRIPT
2.3 Klinkenberg Effect
TE D
M AN U
SC
RI PT
In last section, our theoretical analyses clarified that in slip, transitional, and free molecular flow regimes, the apparent gas permeability of rock matrix is higher than the intrinsic permeability. To account for this permeability enhancement, Klinkenberg factor was proposed to lump all the favorable factors together into the Klinkenberg’s correlation for reservoir simulation (Klinkernberg 1941). Practically in conventional gas reservoir engineering, the Klinkenberg or gas-slippage effect is ignored, except when analyzing pressure responses or flow near gas production wells at very low pressure in some cases. In shale and tight gas reservoirs, however, the Klinkenberg effect is expected to be significant, because of the nano-pores in such rocks, even under high pressure conditions. Figure 7 presents Soeder’s data (1988) that were measured under a confining pressure of 3,000 psi with a core sample acquired at 7,449 ft of the West Virgina No. 6 well. It can be seen that the gas permeability is 20 µD under pore pressure of 1,000 psi, and the gas permeability increases to 55 µD under pore pressure of 80 psi, owing to the increases of the mean free path of gas molecules and the enhanced slippage by pore pressure decrease.
Figure 7. The relationship between pore pressure and gas permeability of Marcellus Shale based on core measurements (Soeder 1988).
AC C
EP
In addition, Bryant and Sakhaee-Pour (2011) employed network modeling to investigate the impact of decreasing reservoir pressure on gas permeability. They calculated the ratio of the gas permeability (kg2,in-situ) at low pressures to the permeability (kg1,in-situ) at the initial reservoir pressure of 28 MPa as a function of reservoir pressure (Figure 8) by considering the impact of pressure on both adsorption layer thickness and gas slippage. The network modeling results indicated that the permeability during the initial production stage is strikingly lower than that during late period of production under lower reservoir pressure. It was claimed that the diminishing adsorbed layer is responsible for this impact. Moreover, the relationship between the permeability and the pressure is found to be non-linear. They concluded that unlike the conventional reservoirs, the gas permeability of shale highly depends on the reservoir pressure. These implications are useful for improving reservoir modelling and ultimate recovery prediction and estimation.
8
SC
RI PT
ACCEPTED MANUSCRIPT
2.4 Non-Linear and Non-Darcy Flow
M AN U
Figure 8. The ratio of gas permeability at low pressures, kg2,in-situ, to that at the initial pressure (P = 28 MPa), kg1,insitu, increases as production depletes and pressure drops (Sakhaee-Pour and Bryant, 2011).
EP
TE D
The dynamics of unconventional gas reservoirs are featured by non-linear multiphase flow in ultra-low permeability formation with multiscale heterogeneity, i.e. non-Darcy flow and strong rock-fluid interaction coexisting inside nano-pores and micro-fractures. Hence traditional REV-based Darcy’s law model may not be practical to characterize such flow dynamics. Take the traditional double-porosity model as an example, transient flow regime between fracture and matrix lasts long time period and it also needs long time to reach the pseudo steady-state flow regime. Therefore, it is not applicable to simulate the fluid flow in the fracture-matrix system of unconventional reservoirs. Two comprehensive literature reviews on flow mechanisms in unconventional gas reservoirs performed by Blasingame (2008) and Moridis et al. (2010) pointed out that high-velocity flow might be dominant during the production of shale gas. The non-laminar/non-Darcy flow definition regarding high velocity gas flow cannot be represented by traditional Darcy's law, but Forchheimer’s equation as shown below is probably sufficient for many applications.
−
∂p µ v = + βρ v 2 ∂x k
(5)
AC C
Where p is pressure, µ is viscosity, k is permeability, ρ is density, β is the non-Darcy flow coefficient (Wu 2016). Gas adsorption and desorption on organic matter and solid grains as well as diffusion of gas molecules into flow streams from low permeability rock matrix also contribute to the non-linear flow characteristics in shale and tight gas reservoirs. 2.5 Effect of Geomechanics
In conventional reservoirs, effect of geomechanics on rock deformation and permeability is generally small and is ignored mostly in practice. However, in unconventional reservoirs, such compaction effect can be significant and have a remarkable impact on both fracture and matrix permeability, which has to be taken into account during reservoir simulation. The impact of confining pressure on permeability is owing to the reduction of porosity. Figure 9 (Wang et al. 2009) shows that permeability of the Marcellus Shale is dependent on reservoir 9
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
pressure. It decreases with increasing effective stress, and the same trends can be seen for other shale reservoirs (Bustin et al. 2008). It is obvious that, as the confining pressure increases, the magnitude of permeability reduction in shale is significantly bigger than those of conventional sandstone or carbonate reservoir.
Figure 9. Effect of effective stress on shale gas permeability (Wang et al. 2009).
AC C
EP
TE D
Holditch (2006) compared air permeability values measured at both surface condition and original overburden pressure, as shown in Figure 10, to illustrate the effect of net overburden pressure or total stress on air permeability of Travis Peak core samples drilled from two wells in eastern Texas. For the cores with high permeability of more than 10 md, the permeability at the initial overburden pressure turns out to be relatively smaller than that of the same unstressed core plugs. Nonetheless, the impact of net overburden pressure increases significantly as the permeability of the cores decreases. In terms of the core plugs with unstressed permeability of 0.01 md, the permeability at net overburden pressure becomes ten times lower (0.001 md). The low permeability rocks with smaller pore sizes are much more sensitive to overburden stress than the rock of conventional reservoirs (Holditch 2006). In order to better understand the characteristics of the tight gas formations, such as gas permeability, capillary pressure, formation factor, resistivity index, acoustic velocity, and the rock mechanical properties, special core analyses must be conducted (Soeder and Randolph 1987).
Figure 10. Air permeability values at net overburden pressure versus those at ambient pressure for Travis Peak cores (Holditch 2006) 10
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
The impact of stress changes on fracture permeability in shale is expected to be significant, in sense that the fracture of 2D nature cannot support much stress load in the normal direction of fracture planes. Cipolla (2009) summarized some laboratory studies by approximating the flow conductivity of partially propped and unpropped fractures as a function of closure stress, as shown in Figure 11. The black curve shows the flow conductivity of an unpropped fracture in which the faces of the fracture planes are aligned upon closing. If this is true, very low conductivity would remain for the unpropped fractures under closure stresses over 2,800 psi, which are common closure stress in most of the shale reservoirs. Nonetheless, if the fracture planes are displaced from each other owing to shear offsets or are partially propped, the fracture conductivity as indicated by the blue curve will be enhanced significantly. Furthermore, the orange curve shows that the fracture conductivity can be retained by high-strength proppants, e.g. sintered bauxite, with increasing closure stress. Since the proppant concentration is much less than 0.1 lbm/ft2 in typical fracturing treatments of shale reservoirs, improvement of flow conductivity by partially propping would not be so prominent due to extra proppant crushing. However, shear offsets between fracture planes could happen, thus enhancing unpropped fracture conductivity to several md-ft in unconventional gas reservoirs.
EP
Figure 11. Impacts of closure stress (effective stress) on conductivity of unpropped and partially propped fractures (Cipolla 2009).
AC C
The effect of closure stress on propped fracture conductivity has been investigated for tens of years, however, the flow conductivity of unpropped fractures has not been well studied until recently. It is possible that some of the fracture networks created by hydraulic fracturing could be unpropped or only partially propped (Cipolla et al. 2008 and 2010). As a result, the magnitude of flow conductivity in unpropped or partially propped fracture could restrict the gas production of shale gas reservoirs (Cipolla et al. 2008). The impact of closure stress on unpropped or partially propped fracture conductivity for rocks with larger Young’s modulus was also investigated by experiments. Figure 12 gives the relationship between stress and unpropped fracture conductivity for Barnett Shale under different Young’s moduli. Fracture conductivity dramatically decreases with increasing closure stress and decreasing Young’s modulus (Cipolla et al. 2010), as can be tens of folds in magnitude.
11
ACCEPTED MANUSCRIPT
2
10
6
E=5x10 psi 6
E=4x10 psi 6
1
RI PT
10
0
10
-1
10
0
1000
2000 3000 Stress, psi
SC
Conductivity, mD-ft
E=2x10 psi
4000
5000
M AN U
Figure 12. The relationship between closure stress and flow conductivity of unpropped fracture under different Young’s modulus (Cipolla et al. 2010)
AC C
EP
TE D
As another example, Figure 13 estimates the impact of closure stress on flow conductivity of unpropped fracture in Marcellus shale under Young’s modulus of 2×106 psi based on previous studies (Cipolla et al. 2008) and the history matching results of Barnett Shale (Cipolla et al. 2010). The original flow conductivity of the fracture network tends to be around 2 md-ft. It decreases to 0.02 md-ft as the pressure inside the fractures approaches the flowing bottom hole pressure of 500 psi.
Figure 13. Effect of closure stress on unpropped fracture conductivity in Marcellus Shale (Cipolla et al. 2010) Gas recovery from shale and tight reservoirs is only about 10-30% of GIP, which is much lower than conventional gas reservoirs. Therefore, to promote the development of new technologies or approaches for enhanced gas recovery (EGR), it is necessary to clarify the fundamental mechanisms behind the coupled physical processes occurring in shale and tight gas reservoirs, including slow flow processes, poor connectivity, large immobile zones, high heterogeneity, stress sensitivity of rock, and slow decrease in long-term production rate. However, in-depth laboratory or theoretical studies on gas flow and transport in unconventional reservoirs are still lacking, and how these complicated mechanisms impact gas flow, production, and the ultimate gas recovery rate is still ambiguous. Future research topics in this area may aim at the investigation of surface diffusion inside the 12
ACCEPTED MANUSCRIPT
adsorption layer, the effect of thermal creep on the solid wall, and the effects of tortuosity of porous media on gas flows.
M AN U
SC
RI PT
3. Gas Flow and Transport Modeling Because of lacking quantitative correlations for corresponding gas flow and transport mechanisms inside unconventional gas reservoirs, there are no effective reservoir simulators or efficient modeling approaches available for optimizing shale or tight gas reservoir development (MIT 2011). Traditionally, the modeling approaches for coalbed methane formations have been used for modeling desorption and diffusion effect in shale gas reservoirs (e.g., King and Ertekin 1991, Zhao 1991). Recently, many analyses and discussions have been developed to characterize fluid flow behavior in shale and tight gas reservoirs (e.g. Javadpour et al. 2007, Blasingame 2008, Javadpour 2009, Wu et al. 2009, Cipolla et al. 2010, Moridis et al. 2010, Sun et al., 2010, Wu and Fackahroenphol 2011, Sakhaee-Pour and Bryant 2011, Shabro et al. 2011, Wu et al. 2012, Wu and Wang 2012, Yao et al., 2012, Li et al., 2016). In particular, Blasingame (2008) comprehensively discussed the fundamental technologies in geology, petrophysics, and reservoir engineering that can be applied to evaluate well/reservoir performance and characterize flow behavior in gas reservoirs with low permeability. In this section, relevant findings and conclusions regarding flow mechanisms and physical processes in shale and tight gas reservoirs are discussed. 3.1 Langmuir Isotherm Model
TE D
Many studies (e.g. Mengal and Wattenbarger 2011, Silin and Kneafsey 2011, EIA 2011b, Wu et al. 2012, Li et al., 2016) showed that the mass of adsorbed gas in formation can be correlated to gas pressure according to the Langmuir isotherm (1918). Adsorbed gas content and adsorption isotherm can be measured in laboratory using core samples. By definition, adsorbed gas content is the amount of gas being adsorbed onto the pore surface of the reservoir rock, and adsorption isotherm represents the capacity of the reservoir rock to adsorb gas under a certain pressure while keeping the temperature constant. Then, they are employed to derive the relationship between pressure and gas storage capacity of the reservoir rock.
VE = VL
AC C
EP
In Langmuir isotherm, adsorbed gas volume, VE, is generally measured in cubic feet of gas per ton of net shale (Mengal and Wattenbarger 2011, EIA 2011b),
P P+ PL
(6)
where VL is the Langmuir volume in scf/ton at infinite pressure, P is reservoir gas pressure, and PL is the Langmuir pressure, at which 50% of the gas is desorbed. In general, Langmuir volume is a function of the TOC content and thermal maturity of shale. An example of gas adsorption isotherms measured from low TOC and high TOC Marcellus Shale samples are shown in Figure 14 (EIA 2011b).
13
ACCEPTED MANUSCRIPT
Adsorbed Gas Content: Higher TOC
SC
RI PT
Adsorbed Gas Content: Lower TOC
Figure 14. Marcellus Shale methane adsorption isotherms (EIA, 2011b)
AC C
EP
TE D
M AN U
As discussed above, adsorbed gas contributes significantly to GIP in shale reservoirs and should be accounted for in flow and production performance assessment. Resource concentration of a prospective area in shale basin is evaluated as the combination of free gas in pores and adsorbed gas that can be released at a finite decrease in pressure. Figure 15 illustrates the relative contributions of free gas (porosity) and adsorbed gas (sorbed) to the total GIP as a function of pressure, as adapted from EIA report (2011b). It indicates that gas reserves in shales are dominated by adsorption under low pressure.
Figure 15. Free, adsorbed, and total gas content versus pressure (EIA, 2011b).
3.2 Slip Boundary Condition
To accurately simulate the velocity slippage, both the slip boundary condition on the wall and the velocity distribution inside the Knudsen layer should be considered. The first-order slip boundary condition for gas flow through microscale space is proposed by Maxwell (1879)
uslip = −
2 −σ
σ
14
∂u Kn ∂n s
(7)
ACCEPTED MANUSCRIPT
where σ is the TMAC. uslip is the velocity on the wall, and ( ∂u ∂n ) s is the velocity gradient along the surface normal. As can be seen, the entire flux is accelerated by the velocity slippage on the boundary, which is proportional to the velocity gradient.
RI PT
Since Maxwell’s pioneering work, many slip boundary conditions have been brought up. In general, it is believed that the slip boundary condition depends on Kn and the velocity condition at the wall, as can be expressed in the following form
∂ 2u ∂u uslip = C1 Kn − C2 Kn 2 2 ... ∂n s ∂n s
(8)
M AN U
SC
where C1 and C2 could be either constants or functions of Kn . If C2 is zero, then the slippage velocity reduces to the first-order approximation. The slip velocity can be derived from kinetic theory (Wu and Bogy 2003, Wu 2008), DSMC (direct simulation Monte Carlo) method (Pan et al. 1999), Lattice Boltzmann Method (LBM), and Burnett or Super-Burnett equation with Chapman-Enskog expansion (Loyalka et al. 1975, Lockerby and Reese 2008). Zhang et al. (2012) conducted a thorough literature review of all the existing slip boundary models and their comparison of the predicted values with experimental results concluded that Wu’s model (2008) derived from the gas kinetic theory is the most accurate among all non-empirical correlations. In Wu’s model (2008) 2 2 3 − σ f 3 3 3 (1 − f ) 1 2 , C2 = f 4 + C1 = 1 − f 2 ) − 2 ( 3 σ 2 Kn 4 Kn
where f = min [1 Kn ,1] .
(9)
AC C
EP
TE D
To calculate the velocity profile in the Knudsen layer, there are mainly three models available, namely the wall-function model (O’Hare et al. 2007), the higher-order continuum model (Lockerby et al. 2005) and the power-law model (Aidun and Clausen 2010). The wall-function model is a phenomenological approach. The higher-order continuum model is developed based on flow governing equations, e.g. the Burnett, BGK–Burnett, and super-Burnett equations. According to the power-law model, the velocity inside the Knudsen layer is defined as a power function of distance away from the wall. By extrapolating the results from LBM and DSMC method, Lilley and Sader (2007, 2008) further promoted the power-law model. However, the singularity of the velocity profile predicted by power law is unphysical and brings difficulty in calculating the slip boundary. Fichman and Hetsroni (2005) derived the velocity profile inside the Knudsen layer by expanding the velocity field and calculating the shear stress at each location. Maxwell’s first-order slip boundary condition was adopted on the boundary. However, there are some mistakes in the calculation of the viscosity in the diffuse reflection case. All in all, a unified model that is capable of calculating both the slippage boundary condition and the velocity profile inside the Knudsen layer is still desirable. Theoretical analysis often assumes the flow space as an ideal straight channel or tube. Although it can provide us insights into the physical problems, in view of the complex pore structures in unconventional reservoirs, it cannot fully represent or replace experimental study. In the following section, we present several apparent permeability correction correlations obtained either from theoretical analyses or experiments.
15
ACCEPTED MANUSCRIPT
3.3 Apparent Permeability Correction Due to velocity slippage and gas molecular diffusion, apparent gas permeability of shale and tight gas matrices is higher than the absolute permeability. To facilitate the inclusion of this phenomenon into reservoir simulators, many multipliers have been proposed to correlate the apparent permeability to the absolute permeability in different flow regimes. Representative multipliers are introduced below.
(10)
SC
b ka = 1 + k∞ p
RI PT
3.3.1 Klinkenberg’s correlation for slip flow regime Experimental studies on permeability enhancement effects in tight formations date back to the early 20th century. In 1941, based on gas flooding experiments, Klinkenberg (1941) correlated the apparent gas permeability via ka to the absolute permeability
M AN U
where b is the Klinkenberg factor and p is the average pressure across the core. The Klinkenberg factor is usually obtained by matching experimental data. Compared with the slip boundary conditions discussed in the previous section, Klinkenberg’s correction is equivalently first-order since Kn ∝ 1/ p . Klinkenberg’s correction can be applied in the low Knudsen number range (<0.1), therefore it is widely adopted for simulating low permeability gas reservoirs. The Klinkenberg factor can be expressed as a function of the absolute permeability and the rock porosity. For example, Jones and Owens (1979) formulated it as
b = 12.639 ( k∞ )
−0.33
, Sampath and Keighin (1981) further integrated porosity into it b = 13.851( k∞ / φ )
Recently, square root correlation was proposed by Civan (2010) b = β ( k∞ / φ )
−1/2
, where
−0.53
.
depends on
TE D
molecule type. All the above Klinkenberg factor correlations are obtained by matching experimental data. For unconventional gas reservoirs, since the rock property varies significantly, usage of the above formulations should be cautious.
EP
3.3.2 Apparent permeability for transitional flow regime So far, no satisfactory apparent permeability correction has been developed for the transitional flow regime due to its complexity. Beskok and Karniadakis (1999) brought up a general relationship between Knudsen number and the apparent permeability multiplier
4 Kn f ( Kn ) = (1 + α1 Kn ) 1 + 1 − α 2 Kn
(11)
AC C
where α1 and α 2 are empirical parameters, depending on the cross-section shape of the flow channel. The development of accurate apparent gas permeability in transitional flow regime is still under pursuit. 3.3.3 Apparent permeability for free molecular flow regime Within the free molecular or Knudsen flow regime, the apparent gas permeability can be calculated by considering the diffusivity for Knudsen diffusion from gas kinetics
1 1 8 RT D = du = d 3 3 πMA where u is the thermal velocity of gas molecules, R is the gas constant, and MA is the gas molecular weight.
16
(12)
ACCEPTED MANUSCRIPT
The derivation of Knudsen diffusion coefficient can be found in Knudsen (1909) and Pollard and Present (1948). The apparent (total) gas flux Jt can be calculated by combining the convective flux J c resulted from the absolute permeability and the flux J k resulted from the Knudsen diffusion
k (13) ∇ p = J t = J c + J k = − ∞ ∇ p + D∇ C µg µg By rearranging the above equation, we can get an apparent gas permeability formulation similar to the Klinkbenberg’s correlation (Ertekin et al. 1986, Sun et al. 2015) ka
b ka = 1 + a p is gas viscosity, and
(14)
is gas compressibility factor.
SC
Where ba = pc g µ g D / k ∞ ,
k∞
RI PT
−
3.4 Non-linear and Non-Darcy Flow Model
M AN U
Javadpour et al. (2007) used a series of Gaussian distribution functions to match the experimental data and reproduced the Knudsen diffusion coefficient. However, it should be noted that the Knudsen diffusion coefficient is actually a function of Kn and the Knudsen diffusivity above is only the extreme case when Kn is infinitely large, as pointed out by Pollard and Present (1948).
In the past few years, a series of experimental studies of high velocity non-Darcy flow in propped hydraulic fractures were carried out at Colorado School of Mines (CSM) (e.g. Lopez 2007, Lai 2010, Lai et al. 2012). These studies concluded that the Barree and Conwey model (2004) yields a much better match than Forchheimer’s equation to experimental data of non-Darcy flow at extremely high flow rate, as shown in Figure 16.
TE D
Here Reynolds number is defined as
N Re =
ρv µτ
(15)
EP
where ρ is density (g/cm3), v is superficial velocity (cm/s), µ is viscosity (cp), and τ is the inverse of the characteristic length (100/cm).
AC C
Permeability ratio is defined as
η=
kapp kd
=
kmin 1 − kmin kd + kd (1 + N Re )
(16)
where kapp is apparent rate-dependent permeability, kd is constant Darcy permeability at low flow rate, and kmin is the minimum permeability at high flow rate.
17
ACCEPTED MANUSCRIPT
707.2 m3/D
SC
η
RI PT
283,168.2 m3/D
M AN U
NRe Figure 16. Experimental data of non-Darcy flow through propped fractures (markers) agrees very well with the Barree and Conway model (dashed line) (Lai et al. 2012).
Figure 16 shows that non-Darcy flow happens under high flow rates and the Barree and Conway (2004 and 2007) model is capable of matching the experimental data in the entire range from low to high flow rates. However, it may be insufficient to use Forchheimer’s equation to characterize the observed flow behavior (Lai et al. 2012).
AC C
EP
TE D
Reservoir simulations conducted at CSM, however, found that the Barree and Conway model gives very similar results to Forchheimer’s equation when dealing with single-phase flow using equivalent parameters (AlOtaibi and Wu, 2010 and 2011), as shown in Figure 17. Even for multiphase flow, the two models achieve very similar results if equivalent parameters are used (Wu et al., 2010; Wu et al. 2011a and 2011b; Zhang et al., 2014). Nevertheless, the laboratory studies (Lai et al. 2012) indicated that if extremely high flow rates through hydraulic fractures occur, the Barree and Conway model is better and preferred, while under normally high velocity flow condition, Forchheimer’s equation should be sufficient for non-Darcy flow modeling in fractured shale gas reservoirs.
18
ACCEPTED MANUSCRIPT
Figure 17. Pressure-transient responses of non-Darcy flow in a naturally fractured reservoir using Forchheimer’s equation and the Barree and Conway model (Al-Otaibi and Wu, 2011).
RI PT
Note that most studies on non-linear flow behavior in both conventional and unconventional reservoirs focused on high velocity flow. There are very few studies or little understanding on non-linear flow at low flow rate or low pressure gradient, which has been observed in experiments (e.g. Zeng et al. 2010, Lei et al. 2008). Based on laboratory observations, it is speculated that there may exist a certain threshold value of pressure gradient, below which flow in extra-low permeability rock may not occur. This threshold pressure phenomenon may also exist in shale and tight gas reservoirs or heavy oil reservoirs.
3.5 Reservoir simulation coupling geomechanics
TE D
M AN U
SC
In the past two decades, Thermal-Hydrological-Mechanical-Chemical (THMC) modeling have been widely used in simulating various subsurface engineering problems such as CO2 sequestration, geological disposal of nuclear waste, conventional and unconventional oil/gas production, acid gas injection, etc. Modeling THMC processes are very challenging in consideration of three aspects: 1) it is difficult to develop a THMC code that is fully functional for simulating all the processes, 2) finding constitutive relationships for couplings that are theoretically sound requires interdisciplinary expertise and there is no universal solution, 3) the realistic validation of THMC models is restrained by very rare data sets that can cover all THMC processes. For gas reservoirs without water influx from a contiguous water-bearing zone, rock compaction effect is generally non-negligible. Reservoir pressure could decrease significantly due to production, thus leading to rock deformation, which could be particularly remarkable for over-pressured shale and tight gas reservoirs. For example, Barnett Shale has a higher pore pressure gradient of about 0.54 psi/ft (Bowker 2007, Cipolla et al. 2010) than normal pressure gradient about 0.45 psi/ft. Thus pore pressure decrease induced by production could significantly affect in-situ stress conditions, resulting in rock compaction.
AC C
EP
THM codes have been available since the early 1980s, mostly limited to single phase liquid flow (e.g. Noorishad et al. 1984). A growing number of THM codes for simulating nuclear waste disposal near-field processes were developed in the 1990s. One common approach for simulating THM processes is to link an existing TH code to a mechanical code such as TOUGH-FLAC (Rutqvist 2011) and IPARS-JAS3D (Minkoff et al. 2003). For single phase flow, most of the THM codes solve the coupled equations simultaneously in a so-called fully coupled (monolithic) method, including ROCMAS (Noorishad and Tsang 1996, Rutqvist et al. 2001), COMPASS (Thomas and Sansom 1995), FRACON (Nguyen 1996), THAMES (Ohnishi and Kobayashi 1996), CODE_Bright (Olivella et al. 1994), and Open-Geosys (Wang and Kolditz 2007). For some applications, a multiphase flow approach is required. Exemplary codes that solve the multiphase flow and mechanics using the fully coupled (monolithic) methods include DYNAFLOW (Preisig and Prévost 2011), TOUGH-CSM (Winterfeld and Wu 2012), CODE_Bright (Olivella et al. 1994), and COMPASS (Thomas and Samson 1995). Recently, Kim (2010) advocated a method of sequential coupling of geomechanics and multiphase flow. There are many coupling methods developed for the THM models (Settari et al. 2001, Longuemare et al. 2002, Minkoff et al. 2003, Tran et al. 2004, Samier et al. 2008, Kim et al. 2009, Zhang 2013, Wang et al., 2014, 2015 and 2016), which can be categorized into three kinds: loose coupling, iterative coupling and full coupling. 1. Loose coupling method: For a loose coupling method, the reservoir simulator performs fluid flow calculations at each time step and the flow solutions are passed to the geomechanical model at a selected time step for stress calculations. This approach is also called one-way coupling because only flow solutions are passed for geomechanical calculations while geomechanical solutions do not feedback to 19
ACCEPTED MANUSCRIPT
flow calculations. Basically, two sets of equations for heat and fluid flows as well as geomechanics are solved in the algorithm. The loose coupling method solves the geomechanics after a certain number of time steps of fluid flow. Several THM simulators employ the loosely coupled procedure to solve fluid flow and geomechanics, e.g. ATH2VIS (Longuemare et al. 2002) and IPARS-JAS3D (Minkoff et al. 2003).
M AN U
SC
RI PT
2. Iterative coupling method: The iteratively coupled procedure solves the primary variables for heat and fluid flows and geomechanics sub-systems separately and sequentially in each time step. Usually the fluid flow equation systems are solved first and the intermediate solutions are then passed to the geomechanics system (Xiong et al., 2013; Zhang et al., 2016a and 2016b). The solutions of geomechanics equations then feed back to fluid flow system until the total equation systems converge, e.g. GEOSIM (Settari et al. 2001) and Rocflow (Wang and Kolditz 2007). This approach is a two-way coupling in each time step (Longuemare et al. 2002). One typical example of this iteratively coupled method is explicitly coupled method, which allows one iteration per time step, e.g. FRACTure (Kohl et al. 1997). Another implementation of the sequentially coupled approach is linking an existing TH code to a mechanical code, such as TOUGH_FLAC (Rutqvist and Tsang 2002, Rutqvist et al. 2002). Kim (2010) also studied the sequential coupling method for multiphase flow considering geomechanics.
TE D
3. Full coupling method: As the tightest coupling method, non-linear partial differential equation systems describing heat and fluid flows as well as geomechanics are solved simultaneously in each time step (Minkoff et al. 2003; Zhang et al., 2012a, 2012b and 2012c; Wu et al. 2014; Zhang et al., 2016c). Most of the fully coupled codes are limited to single phase liquid flow (e.g. Noorishad et al. 1984). A growing number of codes have been developed to simulate nuclear waste disposal near-field processes since the 1990s, such as FRACON (Nguyen 1996), ROCMAS (Noorishad and Tsang 1996, Rutqvist et al. 2001), and Open-Geosys (Wang and Kolditz 2007). There are also some codes that fully couple multiphase flow and mechanics, including CODE_Bright (Olivella et al. 1994), COMPASS (Thomas and Sansom 1995), DYNAFLOW (Preisig and Prévost 2011a), and TOUGH-CSM (Winterfeld and Wu 2012).
AC C
EP
Besides the common THM coupling studies abovementioned, some THMC procedures were recently developed to solve the geomechanics coupled flow problems. Yin et al. (2011) developed a fully coupled THMC model for CO2 injection by finite element methods to assess the potential formation damage around the wellbore. To avoid the calculation complexity, some THMC models were developed by linking two existing codes. FLAC3D and TOUGH-REACT has been linked together to investigate the impacts of hydraulic, thermal, and chemical processes on the evolution of enhanced geothermal system and nuclear waste disposal system (Taron et al. 2009, Rutqvist et al. 2014). Kim et al. (2015) presented a sequential implicit algorithm of chemo-thermo-poromechanics for fractured geothermal reservoirs, by linking ROCKMECH to TOUGH-REACT. Thereby, sequential iteration approach is proved sufficient to couple THMC processes together. To evaluate the coupling approaches, Tran et al. (2009) proposed three indexes: accuracy, adaptability, and running speed. Accuracy refers to how close the numerical results are to the real or benchmark solutions. Adaptability, in another word, flexibility, means how easy the existing or mature flow simulators and geomechanics simulators can be coupled without large code change or subsequent maintenance. Running speed is related to computational efficiency and is an important factor for practical full-field simulations. Explicit coupling approach has very good adaptability and high running speed because of loose coupling between two independent simulators but has poor accuracy due to one-way data transfer, usually from flow to 20
ACCEPTED MANUSCRIPT
geomechanics only at selected time steps. Iterative coupling also has quite good adaptability but less running speed than explicit coupling because the geomechanics computations are performed at each time step, on the other hand, its two-way data exchange between flow and geomechanics yields better accuracy than explicit coupling. Full coupling approach has the best accuracy and is unconditionally stable, but it requires much coding work for the tight coupling and does not run as fast as the other two methods.
Change of rock properties
SC
Geomechanical Effect
Change of pore pressure
M AN U
Gas Flow
RI PT
The primary interest to incorporate rock compaction into reservoir simulation lies in its influence on fluid flow, primarily through the change of rock properties, i.e. porosity and permeability (Zhang et al. 2015). Thus it is a two-way interaction between geomechanical effect and gas flow, where gas flow changes pore pressure and insitu stress, and stress-induced changes of porosity and permeability in turn affect fluid flow, as shown in Figure 18.
Figure 18. Two-way fluid flow and geomechanics interaction. Current industry practice in coupling geomechanical effect with reservoir simulation is to correlate the change of rock properties directly with the change of pore pressure instead of stress. For example, the conventional uncoupled reservoir simulator usually approximates the change of porosity as a function of pore pressure through pore volume compressibility by c P−P φ = φ0 e p ( 0 ) ≈ φ0 1 + c p ( P − P0 )
(17)
k
γ ( P ) = TM ( P) *
AC C
EP
TE D
where P0 is the reference pore pressure at which the porosity is 0. This equation approximates the change of porosity as a function of pore pressure with constant pore volume compressibility, which is one simplified method to capture rock deformation effect (Aziz and Settari 1979, Ertekin et al. 2001). In addition to the porosity approximation, conventional uncoupled reservoir simulator also includes rock compaction effect on absolute permeability through correlations with pore pressure. For example, commercial simulator Eclipse (2016) correlates absolute permeability and pore pressure through a transmissibility multiplier ij +
1 Aij 2
dij
(18)
where γ is the pressure dependent transmissibility, TM is a pressure dependent transmissibility multiplier, k
ij +
1 2
is the absolute permeability, Aij is the fluid flow area, and dij is the length for pressure gradient determination. It has been realized that these porosity and permeability approximations are not sufficient to capture the flow behavior in stress-sensitive reservoirs, therefore a variety of methods for improving coupling between fluid flow and geomechanics have been proposed (Dean et al. 2006, Gutierrez et al. 2001, Minkoff et al. 2003, Settari and Walters 2001, Tran et al. 2009, Xiong et al. 2013, Zhang et al. 2015). In geomechanics, rock properties such as porosity and permeability are functions of effective stress rather than pore pressure, as are calculated along with
21
ACCEPTED MANUSCRIPT
flow variables during simulation. In other words, the way to calculate stress may vary for different coupling methods, but in general the following correlations are incorporated into coupled geomechanical simulation k = k (σ ) , φ = φ (σ )
where
(19)
is the stress variable.
RI PT
These correlations are able to produce more accurate simulation results for stress-sensitive gas reservoirs.
M AN U
SC
Based on the illustration of major gas flow and transport mechanisms in shale and tight gas reservoirs, corresponding mathematical correlations are discussed and compared. Though progress has been made for describing an individual mechanism, there still lack of integrated reservoir simulators that can couple most or all of these mechanisms in a comprehensive way. On the other side, involvement of multi-physics problems in pervasive reservoirs significantly increases the computation costs required by serial gas reservoir simulators. To improve the computation efficiency of reservoir simulators, many high performance computing schemes that can simultaneously call thousands of processors have been developed for large-scale reservoir simulations (e.g. Dongarra et al. 1989, Wheeler and Smith 1990, Killough and Bhogeswara 1991, Zhang et al. 2001, Fung and Dogru 2008, Wang et al. 2015). Details of high performance computing in reservoir simulations are not discussed here, due to limited space.
AC C
EP
TE D
4. Pore-Scale Modeling Pore-scale modeling is becoming a commercial routine for evaluating oil and gas reservoirs at micro- and nanoscale, both the imaging techniques and simulation tools have been well developed during the past several decades (Blunt et al. 2013). As a prerequisite for fluid flow simulation, pore imaging techniques generate threedimensional digital porous geometries that are representative of the pore structures in the reservoir rocks. Considering the similarity of digital porous geometries to the real pore structures in rock samples, they can be divided into two major types. The first is pore network models, which statistically represent or extract the typical pore patterns or properties by numerically generating virtual pore networks. The pore network models are generally constructed based on either real 2D rock images (Wu, 2006) or ideal geometrical shapes, such as spheres, cylindrical tubes, and irregular grains (Wang, 2014). The second is real pore networks, which are real pore structures obtained from reconstruction of micro- and nano-Computed Tomography (CT) scanning or Focused Ion Beam-Scanning Electron Microscope (FIB-SEM) images. Pore network models are not limited in size as algorithms are scalable, while real pore networks are constrained by the size and number of the images acquired. More or less, the size of real pore networks is a 1000 times as large as its resolution in magnitude. Pore networks are places inside which fluid flow simulations are carried out, therefore, corresponding to pore network models and real pore networks, simulation tools are categorized into two families as well: network modeling and direct modeling, respectively (Blunt et al. 2013). For single phase flow, the former, due to the geometrical simplicity of the pore network models, mostly applied specific forms of Darcy’s law to account for viscous flow inside regular or irregular capillary tubes interconnecting the big pore spaces. The latter, in contrary, directly solves or indirectly recovers Navier-Stokes equation by discretizing the real pore networks by voxels. It is thus evident that computation cost of direct modeling is much higher than that of the network modeling.
22
ACCEPTED MANUSCRIPT
4.1 Challenges of unconventional reservoirs Unconventional reservoirs, with irregular nano-pores residing in both organic matter and inorganic minerals (Kuila and Prasad 2013), pose great challenges for pore-scale imaging and modeling techniques developed for conventional reservoir characterization.
M AN U
SC
RI PT
On one hand, majority of the shale pores are of nanometers and porosity is quite low, requiring high resolution imaging techniques to provide details of pore spaces for digital core reconstruction. Micro-CT with the resolution of hundreds of nanometers is insufficient to describe the nano-pores, which in general have a median pore size of tens of nanometers (Honarpour et al. 2012, Du and Chu 2012). Nano-CT can achieve a higher resolution of 50 nm (Guo et al. 2015), it is close to most of the reported median pore throat diameters in shale reservoir matrices. Thus, nano-CT scanning would probably snap off most of the interconnecting throats that are narrower than the half of the resolution, resulting in smaller porosity and permeability values. During the last ten years, FIB-SEM imaging has been widely applied to examine the mineralogy, organic matter, matrix texture, porosity, pore size distribution, and pore shapes and types of shale and tight reservoir rocks (e.g. Chalmers et al. 2009, Sisk et al. 2009, Curtis et al. 2010, Lemmens et al. 2010). Extraordinarily, FIB-SEM is capable of going down to about 1 nm to reveal more specific features, unmatchable by nano-CT. However, unlike CT, FIB-SEM is destructive because of layer by layer milling, which also ablates the rock matrix between two slices of 2D images, causing 1-layer thickness information missing. This missed information or discontinuity in-between may either create or lose some connectivity between pores adjacent to each other on neighboring images, yet this probability of occurrence would be quite low due to relatively sparse pore distribution in shale/tight rock matrices. In spite of this deficiency, FIB-SEM imaging is the most precise and reliable technique that is currently available for uncovering details of shale and tight rock matrices.
EP
TE D
On the other hand, for the pore network models, it is becoming increasingly complex to reconstruct representative pore network model for shale, in consideration of the strong heterogeneity and low connectivity of irregular pores. The pores in shale are more irregular than those in the conventional rocks, using simplified geometric shapes as analogs is far from representative of real pore structures. Besides, including micro-fractures, pore size covers a wide range with strong heterogeneity, but for a given pore, there is no clearly distinguishable pore bodies and pore throats, which is challenging to be statistically described. Also, the connectivity among pores is very low, and there are many isolated pores existing in the shale matrix. Moreover, to further classify pore walls as organic or inorganic, which affects the fluid solid interaction, adds another variable into the pore network models. These unconventional features of shale are critical, since they directly affect the pore network models to be built, and further determine the governing equations to be used for fluid flow simulation.
AC C
Based on our flow regime analyses in previous sections, slip boundary condition should be incorporated into Navier-Stokes equation or lattice-Boltzmann model in order to accurately simulate the gas slip flow regime under reservoir conditions. Additionally, in pores of a few nanometers, more physical behavior of transitional flow and even free molecular flow should be considered. Moreover, it is clear that adsorption and desorption play a critical role in organic-rich shale formations (Lu et al. 1995). Boţan et al. (2015) employed grand-canonical Monte Carlo to account for Langmuir adsorption of methane and MD to capture the gas transport without defining flow regimes in nano- and subnano-porous domains, then they incorporated the adsorption-transport results into a lattice model for larger scale simulation of gas transport in 3D real pore networks. To predict the permeability of shale and tight porous media, Singh and Javadpour (2016) proposed a Langmuir slip permeability model for slip flow and transitional flow regimes which uses slip coefficients derived from Langmuir adsorption data and considers pore size enlargement due to desorption.
23
ACCEPTED MANUSCRIPT
Herein, our simplified justification of pore-scale modeling is limited to classification of in-situ flow regimes, without considering adsorption and desorption, and fluid states under confined conditions. Regarding the flow regimes existing in unconventional gas reservoirs, recent progresses in network modeling and direct modeling of single phase gas flow in micro- and nano-porous media are reviewed and discussed below.
RI PT
4.2 Pore network models and real pore networks Network models as an analog for porous media was first reported 60 years ago (Fatt 1956), and have gained significant momentum during the last ten years. For the purpose of simulating gas flow in unconventional reservoirs, efforts have been focused on advancing the representativeness and verisimilitude of pore network models to the real pore structures as well as incorporating more flow mechanisms that are beyond viscous flow into the gas flow models.
AC C
EP
TE D
M AN U
SC
To date, several different classes of methods have been well developed for generating pore network models, essentially these methods try to preserve the pore size distribution and pore connectivity by extracting pore structure information from real images (Okabe and Blunt 2004, Wu et al. 2004) or to generate porous media by simulating the sedimentation and diagenesis processes (Bakke and Øren 1997). Using Markov Chain Monte Carlo simulation (MCMC), Wu et al. (2006) developed a stochastic pore space reconstruction method, which is able to create realistic pore network models based on the morphological information extracted from 2D rock images. With one 2D image, by assuming that the image characteristics are statistically the same in all x, y, and z directions, they reconstructed 3D pore networks for tight mudrock, which has a porosity of 11.73% and a permeability of 0.013mD. In addition, given three 2D images obtained from perpendicular planes, the model is capable of generating a heterogeneous rock cube with qualitatively reproduced morphology and porosity (Wu et al. 2006). Ma et al. (2014) used this MCMC method to extract node-bond pore networks from 2D images of mudstone, the one they obtained has a porosity of 2.9% with more than 80% pores less than 10 nm in radius and the others less than 50 nm. Based on SEM images with resolution of 5 nm taken for shale samples, Chen et al. (2015) constructed four shale pore networks using this MCMC method. The reconstructed pore networks have a permeability range of 0.176-26.8 and pore size between 5 and 80 nm with frequency peak around 24 nm, fairly representing the nano-pores in shale matrix. Dong and Blunt (2009) modified maximal ball algorithm (Silin and Patzek 2006) to extract simplified pore networks from micro-CT images. These pore network models are comprised of pore bodies interconnected by pore throats, commonly represented by spheres and cylinders, respectively. Also, a pore body or a uniform throat can have an angular cross-section (Blunt et al. 2013). The coordination number, and pore body and throat radii can be varied as a distribution during pore network generation. By significantly shrinking these parameters, pore network models generated are being used as analogs for shale matrices (Mehmani et al. 2013, Zhang et al. 2015). Øren and Bakke (2002) developed process-based procedure that directly models the grain sedimentation, compaction, and diagenesis processes to generate pore network models for sandstone. It is expected to be applicable for shale pore network generation, but the computation will be intensive as well. Note that studies on pore network models discussed here are some focusing on unconventional reservoirs among many others. In view of the nano-pores in shale and tight matrices, imaging techniques used for conventional reservoirs, like optical microscope and micro-CT, are incapable of capturing the pore structures due to low resolution. To meet the demand of unconventional reservoir characterization, nano-CT and FIB-SEM have been widely applied to probe the nano-pores in shale and tight matrices. Guo et al. (2015) adopted the non-destructive nano-CT to scan the internal structures of four shale samples with a resolution of 65 nm in pixel or 32.5 nm in pore radius. Using a commercial software, 3D pore network models with an average porosity of 4.6% and a pore radius frequency 24
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
peak at 103 nm have been visualized by stitching up the stacks of 2D images. For the shale samples they tested, nano-CT captured the pore size distribution and pore connectivity, and distinguished the pores in organic and inorganic matter, however, for other shale or tight matrices with much smaller pore radius frequency peak, it might be difficult to achieve representative results. FIB/SEM has demonstrated much higher resolution in visualizing nano-pores in unconventional reservoir rocks (e.g. Chalmers et al. 2009, Sisk et al. 2009, Curtis et al. 2010, Lemmens et al. 2010 and 2011, Zhang et al. 2012). By following procedures similar to those for nano-CT, 3D pore network models can be reconstructed. Figure 19 is a single 2D FIB-SEM image of a stack taken for a core sample from a tight sandstone gas reservoir. The dimension of this image is Y=23.65 µm and Z=16.9 µm with a resolution of 50 nm (25 nm in pore radius). Low resolution was used here in considering that the field of view of SEM is negatively correlated to its resolution. It is obvious that the pores or the micro-fractures are highly irregular, several of which are a few microns in length. By stacking up 200 2D FIB-SEM images in the depth, we reconstructed the 3D digital core with the thickness of 10 µm along X-axis, as shown in Figure 20, where the 2D image in Figure 19 faces front. Through further discretization of the grey-scale images into binary ones, we can extract the pore structures out of these images as a 3D real pore network (Figure 21). This pore network indicates dominant flow pathways along X-axis, while in Y and Z-axes, the pore connectivity is very poor. The porosity of this 3D digital core is 6.5% and the intrinsic permeability measured by lattice-Boltzmann simulation is 1.24 µD. Real pore networks built upon FIB-SEM images authentically describe the pore structures of unconventional reservoir rocks, nonetheless, its application is still limited by the size of the view field and cumbersome laboratory work.
Figure 19. SEM image of a core sample from a tight sandstone gas reservoir.
25
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
Figure 20. 3D digital core reconstructed from FIB-SEM images of a tight sandstone core.
Figure 21. 3D real pore networks extracted from FIB-SEM images of a tight sandstone core.
4.3 Network modeling
Once the pore network models for shale or tight rocks are constructed, mathematical models for simulation of gas flow through the pore networks can be built and implemented. Since Darcy’s law derived from NavierStokes equation is no longer accurate for gas slip flow, many modifications to the intrinsic permeability have been proposed to calculate the apparent permeability . One of the earliest correlations that depict gas slip flow through reservoir rocks is the Klinkenberg’s equation (Klinkenberg 1941)
26
ACCEPTED MANUSCRIPT
λ ka = k 1 + 4c r
(20)
where is a factor slightly less than 1.0. Although it was developed for the purpose of matching core-scale data, the basic assumption is the slippage of gas molecules on solid surface defined by Knudsen number.
RI PT
To simulate gas flow in nano-pores in shale matrix, Javadpour et al. (Javadpour et al. 2007, Javadpour 2009) proposed to calculate the total mass flux J t by simultaneously combining the slip flow J s and Knudsen diffusion
Jk . Jt = J s + J k
(21)
TE D
M AN U
SC
The slip correction used comes from Brown et al. (1946) instead of Klinkenberg’s. Their expression for a cylindrical nano-pore matches experimental data measured from pores of ~200 nm in diameter with an average error of 4.5% (Javadpour 2009). Note that here the simplicity of modifying Darcy’s law is preferred over physics as the gas flow goes beyond = 0.1. Mehmani et al. (2013) extended this model (Javadpour 2009) to simulate gas flows in 3D node-bond pore networks with variable pore sizes. The gas slippage and Knudsen diffusion were only evaluated in cylindrical pore throats. Their results show that the pressure-dependent apparent permeability is affected by the fraction and spatial configuration of the nano-pores (Mehmani et al. 2013). By substituting a gas kinetics expression similar to Klinkenberg’s equation into the slip flux (Javadpour 2009), Ma et al. (2013) simulated gas flow through 3D pore network models statistically constructed upon 2D SEM images of shale samples using MCMC method (Wu 2006). Rather than solely using cylindrical pores, the network model contains throats close to triangle and square shapes, too. Also, in the flow modeling, non-ideal gas compressibility calculated by van der Waal’s equation of state was incorporated into the slip and Knudsen flux for comparison with ideal gas cases, revealing that non-ideal gas effect that favors methane flow can contribute tens of percent to the apparent permeability under reservoir conditions. Superposition of slip and Knudsen flux in total flux is questioned as double counting the contributions from slippage and Knudsen diffusion (Sakhaee-Pour and Bryant 2012, Ma et al. 2013), when compared to the dusty gas
where
AC C
viscous flow J v .
EP
model (DGM) developed by Mason and Malinauskas (1983), which only adds Knudsen diffusion J k to the
J t = J v + J k , J k = − DKn
is the Knudsen diffusivity, ∇ is the pressure gradient,
∇p RT
is the gas constant, and
(22) is temperature.
Sakhaee-Pour and Bryant (2012) chose the first-order slip model in the form of Klinkenberg’s equation over DGM to simulate gas flow in cylindrical nano-pores. Their model, which matches experimental data within 6% deviation, found that as gas production continues with pressure drop in reservoirs, the matrix permeability is expected to increase by 4.5 times due to release of adsorbed layer increasing pore radius and gas slippage. This helps explain the deaccelerating decline in production rate of shale gas wells after the early years (Sakhaee-Pour and Bryant 2012). Pore network modeling simulates the physical behavior of micro- and nanoscale gas flow by modifying the macroscopic correlations, providing a means of directly upscaling microscale physics to field scale applications. 27
ACCEPTED MANUSCRIPT
4.4 Direct modeling
RI PT
Direct modeling solves the gas flow by discretizing pore spaces obtained from SEM or nano-CT images into lattices. As shown in Table 1, either Navier-Stokes equation (Silin 2014) or LBM is capable of capturing the slip flow behavior by incorporating the slip boundary condition. In the transitional flow regime, seeing that Burnett equation is tedious to be solved, DSMC is more commonly used. However, for slow flow as encountered in reservoir matrix, DSMC experiences thermal fluctuations toward convergence, and even simulation of flows through nanoscale digital rock cubes is prohibitively expensive (Li and Sultan 2015). For free molecular flow, MD simulation is too computationally expensive for the representative micrometer scale porous media, only flows in nanotubes have been studied (Feng and Akkutlu 2015). Therefore, we will focus on lattice-Boltzmann method for simulating gas slip flow below.
M AN U
SC
LBM is a numerical method that models the evolution of the single-particle velocity distribution function n ( , t) over spatial lattices according to lattice-Boltzmann equation, which is derived by discretizing Boltzmann equation (He and Luo 1997a) and indirectly recovers NS equation (He and Luo 1997b). As a scalable method suitable for parallel computing (Xiao 2013, Wang 214), it has been widely used to simulate continuum fluid flows in complex geometries (Chen and Doolen 1998, Chen et al. 2014). The velocity distribution function describes the fraction of fluid particles with velocity $ at position and time t in each lattice. Velocity quadrature $ is designed to simplify the continuous molecular velocity distribution in the gas kinetics with an assigned weighting factor %&' , which sums to 1. Evolution of n ( , t) in LBM involves propagation and collision steps, as shown in the left and right equalities below
ni ( r + ci ∆t,t + ∆t ) = ni* ( r ,t ) + fi ( r ,t ) = ni ( r ,t ) + ∆i n ( r ,t ) where the post-collision particle velocity distribution ni* is moved from lattice
(23)
r to the neighbors r + ci ∆t in the
TE D
direction of ci , ∆i is the collision operator, and f i ( r ,t ) is the forcing term. In propagation, populations are required to move to the next lattice in each time step.
AC C
EP
When the fluid particles propagate to a rigid surface, boundary conditions must be applied to model the interactions on the surface. Standard LBM simulates the no-slip boundary condition (Yin et al. 2006) and thus only gives the intrinsic permeability of porous media. By properly treating the particle-surface interactions, many LBM schemes have been developed to simulate the slip boundary conditions. As discrete velocity quadratures approximate the continuum gas velocity distribution, some studies increased the number of velocity quadratures to achieve a high order discretization of the Boltzmann equation (Shan et al. 2006, Zhang et al. 2006). Ansumali et al. (2002) applied the completely diffusive boundary on the wall in LBM to simulate slip flow with ( < 0.03. Higher order velocity set for LBM have been combined with the relaxation time- ( correlation for slip simulation as well (Tang et al. 2008). Chen and Tian (2009) used the Langmuir slip model in the LBM to simulate Couette and Poiseuille flow in the slip flow regime, which matched the analytical solutions. In these abovementioned studies, only gas slippage on a flat surface was calculated, while simulation of slip flow along curved surface in 3D geometries is relatively rare. These considerations are important since the curved surface may result in negativity in slippage contribution to the flow velocity (Einzel et al. 1990, Lockerby et al. 2004, Szalmás 2007). Recently, Cho et al. (2013) conducted slip flow simulation through 3D non-overlapping fibers using multiple relaxation time LBM by relating relaxation time to ( and lattice size, and obtained permeability as a function of (, porosity, and fiber configuration. However, in their simulation, the effect of 28
ACCEPTED MANUSCRIPT
TMAC on slip flow is not investigated (Zhang et al. 2005). Besides, correlation of lattice size to relaxation time renders it complicated to vary the resolution of geometries for more accurate simulations.
SC
RI PT
At present, other than simple channels, cubes, and cylinders, very few LBM schemes with slip boundary condition have been tested for flows through complex porous media, which actually would pose many challenges. For example, pore size in porous media covers a range rather than a unique value, therefore ( is not uniform inside the porous media. However, many models, in essential, defined ( as a specific value but a function of pore size distribution (Zhang et al. 2005, Szalmás 2007, Kim et al. 2008, Tang et al. 2008). Hence in slip flow models, ( should be substituted by mean free path to serve as the input to account for the varying degree of rarefaction. Additionally, application of slip boundary condition in porous media is bounded by the largest ( = 0.1 at the smallest pore throats for the standard LBM, because the assumption of thermodynamic equilibrium cannot be met in Knudsen layer (Tang et al. 2008). To conquer these challenges, Wang (2014) derived a new bounce-back scheme by incorporating the Maxwell’s first-order slip boundary condition into the lattice-Boltzmann equation
2 −σ
σ
τ
(24)
M AN U
us = −
where us is the slip velocity, tangential stress vector τ = n ⋅ π ⋅ ( I − nn ) ,
n is the unit normal vector pointing
Τ from solid wall to fluid, I is the unit tensor. π = µ ∇u + ( ∇u ) is the viscous stress tensor of the local fluid,
is the viscosity, ∇u is the velocity gradient tensor, and
here represents the transpose operation.
Setting a correction term δ i ( r ,t ) in the bounce-back scheme (Ladd and Verberg 2001, Yin et al. 2006), the
δ i ( r ,t ) = −
where
)
ρ ac 1 i
cs2 bl
1 ( ci ⋅ n ) ci ⋅ ( I − nn ) ⋅ j ( r ,t ) ρ ( c i ⋅ n )( n ⋅ u b ) + 1 + 2bl
( ci ⋅ n ) ci ⋅ ub − ( ci ⋅ n )( ub ⋅ n ) −
EP
+
2a ci cs2
TE D
following expression for gas slippage correction is obtained (Wang, 2014)
a ci ci ⋅ ( I − nn ) ⋅ Π neq ,* − Π neq ⋅ ( I − nn ) ⋅ ci 2cs4
(
)
(25)
is the isothermal sound speed in a fluid at equilibrium, ρ is fluid density, *+ is the velocity of the rigid
AC C
surface, slip length bl = ( 2 − σ ) λ / σ ,
j
is the first-order momentum neq
Πneq,*
is the post-collision non-
equilibrium part of the second-order moment Π . In this correction, the boundary velocity is retained, and the mean free path of gas molecules and the TMAC from the Maxwell’s first-order slip boundary condition are both embodied. Also, note that when ,- approaches infinity, this correction term recovers the free-slip boundary condition (Yin et al. 2006). Nevertheless, since ,- is in the denominator, it cannot be assigned very small values to recover the no-slip boundary condition. This slip model has been used to simulate slip flow through three arrays of spheres, demonstrating that the apparent permeability increases with increasing slip length and increases with decreasing porosity, which are formulated into polynomials for upscaling purpose. Also, a local normal calculation method was proposed for this slip model to simulate flows through simple cubic arrays of overlapping spheres (Wang 2014). For complex, digitalized pore surfaces inside shale or tight sandstone matrices where the global normal vector is not available and for situations where the numerical stability of the global normal method is an issue, the method of local normal method offers a realistic approximation. 29
ACCEPTED MANUSCRIPT
RI PT
There are other methods that combined pore network flux models with lattice-Boltzmann simulation. For example, Allan and Mavko (2012) applied LBM to calculate the flux-averaged pore width for Knudsen number, which was further substituted into the unified flow model developed by Karniadakis et al. (2005) and Florence et al. (2007) to compute the apparent permeability. Using DGM, Chen et al. (2015) evaluated the contribution of Knudsen diffusion to the apparent permeability of 3D pore networks reconstructed from SEM images of shale by MCMC method. To calculate the Knudsen diffusivity, lattice-Boltzmann method with single relaxation time was adopted to simulate the diffusion process by taking pore size distribution into consideration. In these hybrid models, the contribution of Knudsen diffusion to the apparent permeability is only counted in the flux model.
SC
5. Conclusions In view of the characteristics of shale and tight gas reservoirs that are different from conventional natural gas reservoirs, a few mechanisms dominating gas flow and transport should be properly addressed in gas reservoir simulation, including gas desorption, rarefaction effect, non-Darcy flow, compaction effect, etc.
M AN U
Though many quantitative correlations for describing gas flow and transport mechanisms in shale and tight gas reservoirs have been proposed, there still lack of generalized and well-validated formulations that can recover the commingled gas behavior in shale and tight gas reservoirs with high heterogeneity. Current large-scale reservoir simulators for shale and tight gas reservoir simulation not only suffer from scarcity of accurate correlations for different flow and transport mechanisms, but also undergo high complexity of coupling them together. High performance reservoir simulations integrating all these predominant factors are still in urgent need.
TE D
Statistical pore network models are readily available for reconstructing variable pore networks from real SEM images of unconventional reservoir rocks by extracting the pore size distribution and pore connectivity information. FIB-SEM imaging and 3D digital core construction are the most authentic means to describe the pore characteristics inside unconventional reservoir rocks, the information completeness and field view are still limited though.
EP
Several pore network modeling methods have been developed for simulating gas flows in slip and transitional flow regimes by successfully capturing the molecular level physics in macroscopic expressions.
AC C
There are a few lattice-Boltzmann models developed for gas slip flow simulation, however, very few of them have been extended to simulate slip flow in 3D complex porous media due to geometrical challenges. To simulate gas flow in transitional flow regimes, DSMC that directly simulates the Boltzmann equation by probabilistically moving and colliding molecular swarms is more suitable.
Acknowledgement
The authors thank Dr. David Wood, the Chief Editor of Journal of Natural Gas Science and Engineering, for his invitation to write this comprehensive review paper.
30
ACCEPTED MANUSCRIPT
Reference Aidun, C. K. and Clausen, J. R. 2010. Lattice-Boltzmann method for complex flows. Annual Review of Fluid Mechanics, 42, 439-472. Allan, A.M. and Mavko, G. November, 2012. The effect of adsorption and diffusion on the gas permeability of kerogen. In 2012 SEG Annual Meeting. Society of Exploration Geophysicists.
RI PT
Al-Otaibi, A. and Y.-S. Wu. An alternative approach to modeling non-Darcy flow for pressure transient analysis in porous and fractured reservoirs, SPE-SAS-984, prepared for presentation at the 2011 SPE Saudi Arabia Section Technical Symposium and Exhibition held in AlKhobar, Saudi Arabia, 15–18 May 2011.
SC
Al-Otaibi, A. and Y.-S. Wu. Transient behavior and analysis of non-darcy flow in porous and fractured reservoirs according to the barree and conway model, SPE 133533, presented at the at the Western North America Regional Meeting, held in Anaheim, California, USA, 26–30 May 2010. Ansumali, S. and Karlin, I. V. 2002. Kinetic boundary condition in the lattice Boltzmann method. Physical Review E 66 (2): 026301.
M AN U
Aziz, Khalid, and Antonin Settari. 1979. Petroleum reservoir simulation. Chapman & Hall. Bakke, S. and Øren, P. E., 1997. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. Spe Journal, 2(02), pp.136-149. Barree, R. D. and M. W. Conway, “Multuiphase Non-Darcy Flow in Proppant Packs,” SPE 109561, presented at the 2007 Annual Technical Conference and Exhibition, Anaheim, CA, 11-14 November, 2007.
TE D
Barree, R.D. and M.W. Conway, “Beyond Beta Factors: A Complete Model for Darcy, Forchheimer and TransForchheimer Flow in Porous Media,” SPE 89325, presented at the 2004 Annual Technical Conference and Exhibition, Houston, Texas 26-29 September, 2004. Beskok, A., Karniadakis, G. E. A model for flows in channels, pipes, and ducts at micro and nano scales. Microscale Thermophys Eng. 3(1), 43–77 (1999). doi:10.1080/108939599199864
EP
Bird, G. A. 2003. Molecular gas dynamics and the direct simulation of gas flows. New York. Oxford University Press Inc.
AC C
Blasingame, T.A., The characteristic flow behavior of low-permeability reservoir systems, SPE 114168, presented at the 2008 SPE Unconventional Reservoirs Conference held in Keystone, Colorado, U.S.A., 10–12 February 2008. Blunt, M. J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A. and Pentland, C., 2013. Pore-scale imaging and modelling. Advances in Water Resources, 51, pp.197-216. Boţan, A., Ulm, F.-J., Pellenq, R. J.-M., and Coasne, B., 2015. Bottom-up model of adsorption and transport in multiscale porous media, Physical Review E, 91, 032133. Bowker, K. A. 2007. Barnett Shale gas production, Fort Worth Basin: issues and discussion. AAPG bulletin 91(4): 523-533. Brown, G. P., DiNardo, A., Cheng, G. K. and Sherwood, T. K. 1946. The flow of gases in pipes at low pressures. Journal of Applied Physics, 17(10), pp.802-813.
31
ACCEPTED MANUSCRIPT
Bustin, R. M., A. M. M. Bustin, X. Cui, D. J. K. Ross, and V. S. Murthy Pathi. 2008. Impact of shale properties on pore structure and storage characteristics, presented at the SPE Gas Production Conference, 16-18 November, Fort Worth, Texas, SPE 119892. Bybee, K. 2008. Horizontal wells in tight gas sands: risk management to maximize success. JPT, 60(10): 61-63.
RI PT
Chalmers, G. R., Bustin, R. M., and Power, I. 2009. A Pore by any other name would be as small: the importance of meso- and microporosity in shale gas capacity. Presented at the AAPG Annual Convention and Exhibition, Denver, Colorado. 7-10 June. Chen, L., Kang, Q., Mu, Y., He, Y.L. and Tao, W.Q., 2014. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. International Journal of Heat and Mass Transfer, 76, pp.210-236.
SC
Chen, L., Zhang, L., Kang, Q., Viswanathan, H. S., Yao, J. and Tao, W., 2015. Nanoscale simulation of shale transport properties using the lattice Boltzmann method: permeability and diffusivity. Scientific reports, 5.
M AN U
Chen, S. and Doolen, G. D. 1998. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics 30 (1): 329-364. Chen, S. and Tian, Z. 2009. Simulation of microchannel flow using the lattice Boltzmann method. Physica A: Statistical Mechanics and its Applications 388 (23): 4803-4810. Chien, M. C. H., Wasserman, M. L., Yardumian, H. E., Chung, E. Y., Nguyen, T., and Larson, J. January, 1987. The use of vectorization and parallel processing for reservoir simulation. In SPE Symposium on Reservoir Simulation. Society of Petroleum Engineers.
TE D
Cho, H., Jeong, N., and Sung, H. 2013. Permeability of microscale fibrous porous media using the lattice Boltzmann method. International Journal of Heat and Fluid Flow 44: 435-443. Cipolla, C. L. Modeling Production and Evaluating Fracture Performance in Unconventional Gas Reservoirs, JPT, Distinguished Author Series, pp.84-90, September 2009.
EP
Cipolla, C. L., E. P. Lolon, J. C. Erdle, and B. Rubin, Reservoir Modeling in Shale-Gas Reservoirs, Reservoir Evaluation and Engineering, August 2010.
AC C
Cipolla, C. L., N. R. Warpinski, M. J. Mayerhofer, E. P. Lolon, and M. C. Vincent, The Relationship Between Fracture Complexity, Reservoir Properties, and Fracture Treatment Design, SPE 115769, presented at the SPE Annual Technical Conference and Exhibition, Denver, 21–24 September 2008. Civan, F., Rai, C. S., and Sondergeld, C. H. 2011. Shale-gas permeability and diffusivity inferred by improved formulation of relevant retention and transport mechanisms. Transport in Porous Media, 86(3), 925-944. doi: 10.1007/s11242-010-9665-x Civan, F. 2010. Effective correlation of apparent gas permeability in tight porous media. Transport in Porous Media, 82(2), 375-384. doi: 10.1007/s11242-009-9432-z Clarkson, C. R., Solano, N., Bustin, R. M., Bustin, A. M. M., Chalmers, G. R. L., He, L., Melnichenko, Y. B., Radliński, A. P., and Blach, T. P. 2013. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel, 103: 606-616.
32
ACCEPTED MANUSCRIPT
Curtis, M. E., Ambrose, R. J., and Sondergeld, C. H. 2010. Structural characterization of gas shales on the microand nano-scales. Presented at the Canadian Unconventional Resources and International Petroleum Conference, Calgary, Alberta. SPE-137693-MS. http://dx.doi.org/10.2118/137693-MS. Darabi, H., Ettehad, A., Javadpour, F., and Sepehrnoori, K. 2012. Gas flow in ultra-tight shale strata. Journal of Fluid Mechanics, 710, 641-658. doi: 10.1017/jfm.2012.424
RI PT
Dean, R. H., X. Gai, C. M. Stone, and S. E. Minkoff. 2006. A comparison of techniques for coupling porous flow and geomechanics. SPE Journal, 11(01): 132-140. DeBaun, D., Byer, T., Childs, P., Chen, J., Saaf, F., Wells, M. and Crumpton, P. January, 2005. An extensible architecture for next generation scalable parallel reservoir simulation. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
SC
Denney, D. 2008. Rock type-understanding productivity in tight gas sands. JPT, 10: 53-56.
M AN U
Dogru, A. H., Fung, L. S., Middya, U., Al-Shaalan, T. M., Pita, J. A., HemanthKumar, K., and Hahn, W. A. October, 2009. A next-generation parallel reservoir simulator for giant reservoirs. In SPE/EAGE Reservoir Characterization and Simulation Conference. Dong, H. and Blunt, M. J. 2009. Pore-network extraction from micro-computerized-tomography images. Physical Review E, 80(3), p.036307. Dongarra, J., Duff, I., Gaffney, P. A. T. R. I. C. K., and McKee, S. 1989. Vector and Parallel Computing. Wiley, New York. Du, L. and Chu, L. 2012, January. Understanding anomalous phase behavior in unconventional oil reservoirs. In SPE Canadian Unconventional Resources Conference. Society of Petroleum Engineers.
TE D
Eclipse, http://www.software.slb.com/products/eclipse.
EIA (U.S. Energy Information Administration). July 2011a. Review of Emerging Resources: U.S. Shale Gas and Shale Oil Plays.
EP
EIA (U.S. Energy Information Administration). April 2011b. World Shale Gas Resources: An Initial Assessment of 14 Regions Outside the United States. Einzel, D., Panzer, P., and Liu, M. 1990. Boundary condition for fluid flow: curved or rough surfaces. Physical Review Letters 64 (19): 2269-2272.
AC C
Ertekin, T., King, G. R., and Schwerer, F. C. 1986. Dynamic gas slippage: a unique dual-mechanism approach to the flow of gas in tight formations. SPE (Society of Petroleum Engineers) Format. Eval., (United States), 1(1). Ertekin, Turgay, Jamal Hussein Abou-Kassem, and Gregory R. King. Basic applied reservoir simulation. Richardson, TX: Society of Petroleum Engineers, 2001. Falk, K., Coasne, B., Pellenq, R., Ulm, F. J., and Bocquet, L. 2015. Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nature Communications, 6. Fatt, I., 1956. The network model of porous media.
33
ACCEPTED MANUSCRIPT
Feng, F. and Akkutlu, I. Y. 2015. November. Flow of Hydrocarbons in Nanocapillary: A Non-Equilibrium Molecular Dynamics Study. In SPE Asia Pacific Unconventional Resources Conference and Exhibition. Society of Petroleum Engineers. Fichman, M. and Hetsroni, G. 2005. Viscosity and slip velocity in gas flow in microchannels. Physics of Fluids (1994-present), 17(12), 123102.
RI PT
Florence, F.A., Rushing, J., Newsham, K.E. and Blasingame, T. A. 2007. January. Improved permeability prediction relations for low permeability sands. In Rocky Mountain Oil and Gas Technology Symposium. Society of Petroleum Engineers. Fung, L. S. and Dogru, A. H. (2008). Parallel unstructured-solver methods for simulation of complex giant reservoirs. SPE Journal, 13(04), 440-446.
SC
Geist, A. 1994. PVM: Parallel virtual machine: a users' guide and tutorial for networked parallel computing. MIT press.
M AN U
Gong, X., Tian, Y., McVay, D. A., Ayers, W. B., and Lee, J. 2013. Assessment of Eagle Ford shale oil and gas resources. Presented at the SPE Unconventional Resources Conference, Calgary, Alberta, 5-7 November. SPE-167241-MS. Guo, X., Shen, Y. and He, S. 2015. Quantitative pore characterization and the relationship between pore distributions and organic matter in shale based on Nano-CT image analysis: A case study for a lacustrine shale reservoir in the Triassic Chang 7 member, Ordos Basin, China. Journal of Natural Gas Science and Engineering, 27, pp.1630-1640.
TE D
Gutierrez, M., R. W. Lewis, and I. Masters. 2001. "Petroleum reservoir simulation coupling fluid flow and geomechanics." SPE Reservoir Evaluation and Engineering 4(03): 164-172.
He, X. and Luo, L.-S. 1997. Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation. Physical Review E 56 (6): 681-6817.
EP
He, X. and Luo, L. S. 1997. Lattice Boltzmann model for the incompressible Navier–Stokes equation. Journal of statistical Physics, 88(3-4), pp.927-944. Holditch, S. A. 2006. Tight gas sands. Journal of Petroleum Technology, 58(06), 86-93.
AC C
Honarpour, M.M., Nagarajan, N.R., Orangi, A., Arasteh, F. and Yao, Z., 2012, January. Characterization of Critical Fluid PVT, Rock, and Rock-Fluid Properties-Impact on Reservoir Performance of Liquid Rich Shales. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Javadpour, F., Fisher, D. and Unsworth, M., 2007. Nanoscale gas flow in shale gas sediments. Journal of Canadian Petroleum Technology, 46(10). Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). Journal of Canadian Petroleum Technology,48(08), pp.16-21. Jones, F. O. and Owens, W. W. 1980. A laboratory study of low-permeability gas sands. Journal of Petroleum Technology, 32(09), 1-631. Karniadakis, G., Beskok, A. and Aluru, N., Microflows and nanoflows: fundamentals and simulation. 2005.
34
ACCEPTED MANUSCRIPT
Killough, J. E., and Bhogeswara, R. 1991. Simulation of compositional reservoir phenomena on a distributedmemory parallel computer. Journal of Petroleum Technology, 43(11), 1-368. Kim, J. 2010. Sequential methods for coupled geomechanics and multiphase flow. Ph.D. Thesis, Department of Energy Resources Engineering, Stanford University, California, 264 pp.
RI PT
Kim, J., Sonnenthal, E., and Rutqvist, J. 2015. A sequential implicit algorithm of chemo-thermo-poro-mechanics for fractured geothermal reservoirs. Computers and Geosciences, 76: 59-71. Kim, S. H., Pitsch, H., and Boyd, I. D. 2008. Slip velocity and Knudsen layer in the lattice Boltzmann method for microscale flows. Physical Review E 77 (2): 026704. King, D. A. 1980. Surface diffusion of adsorbed species: A Review. Journal of Vacuum Science and Technology, 17(1), 241-247.
SC
King, E. G. Thirty Years of Gas Shale Fracturing: What Have We learned? SPE 133456, presented at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19-22 September, 2010.
M AN U
King, G.R. and T. Ertekin, State-of-the-Art Modeling for Unconventional Gas Recovery, SPE Formation Evaluation, March, 1991. Klinkenberg, L. J. 1941. The permeability of porous media to liquids and gases. Drilling and Production Practices API-41-200. Knudsen, M. 1909. Molecular effusion and transpiration. Nature, 80, 491-492.
Knudsen, M. 1995. The laws of molecular flow and of inner friction flow of gases through tubes. Journal of Membrane Science, 100(1), 23-25.
TE D
Kohl, T., Evans, K. F., Hopkirk, J., Jung, R. and Rybach, L., 1997. Observation and simulation of non-Darcian flow transients in fractured rock.Water Resources Research, 33(3), pp.407-418. Kuila, U. and Prasad, M. 2013. Specific surface area and pore‐size distribution in clays and shales. Geophysical Prospecting, 61(2), pp.341-362.
EP
Kurihara M, Ouchi H, Inoue T, Yonezawa T, Masuda Y, Dallimore S. R., Collett T. S. 2005. Analysis of the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate thermal production test through numerical simulation. Geological Survey of Canada, Bulletin 585.
AC C
Ladd, A. J. C. and Verberg, R. 2001. Lattice-Boltzmann simulations of particle-fluid suspensions. Journal of Statistical Physics 104 (516): 1191-1251. Lai, B.T. 2010. Experimental Measurement and Numerical Modeling of High Velocity Non-Darcy Flow Effects in Porous Media, PhD dissertation, Colorado School of Mines, Golden, Colorado. Lai, B. T., J. L. Miskimins, and Y.S. Wu, “Non-Darcy Porous Media Flow According to the Barree and Conway Model: Laboratory and Numerical Modeling Studies,” SPE-122611, SPE Journal, Volume 17, Number 1, pp.70-79, March 2012. Lei, Q. W. Xiong, J. Yuan, and Y.-S. Wu, Behavior of flow through low-permeability reservoirs, SPE 113144, presented at the 2008 SPE Europec/EAGE Annual Conference and Exhibition held in Rome, Italy, 9–12 June 2008.
35
ACCEPTED MANUSCRIPT
Lemmens, H. J., Butcher, A. R., and Botha, P. W. S. K. 2010. FIB/SEM and automated mineralogy for core and cuttings analysis. Presented at the SPE Russian Oil and Gas Conference and Exhibition, Moscow, 26-28 October. SPE-136327-RU. Lemmens, H. J., Butcher, A. R., and Botha, P. W. S. K. 2011. FIB/SEM and SEM/EDX: a new dawn for the SEM in the core lab? Petrophysics 52 (6): 452-456.
RI PT
Li, G., Yao, Y., and Zhang, R., 2016. An improved model for the prediction of liquid loading in gas wells. Journal of Natural Gas Science and Engineering, vol. 32, pp.198-204. Li, J. and Sultan, A.S., 2015, December. Permeability computations of shale gas by the pore-scale monte carlo molecular simulations. In International Petroleum Technology Conference. International Petroleum Technology Conference.
SC
Li, X., Abass, H., Teklu, T. W., and Cui, Q. 2016. A Shale Matrix Imbibition Model-Interplay between Capillary Pressure and Osmotic Pressure. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
M AN U
Lilley, C. R. and Sader, J. E. 2007. Velocity gradient singularity and structure of the velocity profile in the Knudsen layer according to the Boltzmann equation. Physical Review E, 76(2), 026315. Lilley, C. R. and Sader, J. E. August, 2008. Velocity profile in the Knudsen layer according to the Boltzmann equation. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences (Vol. 464, No. 2096, pp. 2015-2035). The Royal Society. Lockerby, D. A., Reese, J. M., Emerson, D. R. et al. 2004. Velocity boundary condition at solid walls in rarefied gas calculations. Physical Review E 70 (1): 017303.
TE D
Lockerby, D. A., Reese, J. M., and Gallis, M. A. 2005. Capturing the Knudsen layer in continuum-fluid models of nonequilibrium gas flows. AIAA journal, 43(6): 1391-1393. Lockerby, D. A. and Reese, J. M. 2008. On the modelling of isothermal gas flows at the microscale. Journal of Fluid Mechanics, 604: 235-261.
EP
Lopez, H. D. 2007. Experimental Analysis and Macroscopic and Pore-level Flow Simulations to Compare NonDarcy Flow Models in Porous Media, PhD dissertation, Colorado School of Mines, Golden, Colorado.
AC C
Loyalka, S. K., Petrellis, N., and Storvick, T. S. 1975. Some numerical results for the BGK model: thermal creep and viscous slip problems with arbitrary accomodation at the surface. Physics of Fluids (19581988), 18(9), 1094-1099. Lu, X.C., Li, F.C. and Watson, A.T. 1995. Adsorption measurements in Devonian shales. Fuel, 74(4), pp.599-603. Luffel, D. L., Hopkins, C. W., and Schettler Jr, P. D. January, 1993. Matrix permeability measurement of gas productive shales. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Ma, J., Sanchez, J.P., Wu, K., Couples, G.D. and Jiang, Z. 2014. A pore network model for simulating non-ideal gas flow in micro-and nano-porous materials. Fuel, 116, pp.498-508. Mason, E.A. and Malinauskas, A.P., 1983. Gas transport in porous media: the dusty-gas model (Vol. 17). Elsevier Science Ltd.
36
ACCEPTED MANUSCRIPT
Masuda Y, Naganawa S, Ando S, Sato K. Numerical calculation of gas-production performance from reservoirs containing natural gas hydrates. paper SPE 38291, Proceedings, Western Regional Meeting, Society of Petroleum Engineers, Long Beach, California, 1997.
RI PT
Masuda Y, Konno Y, Iwama H, Kawamura T, Kurihara M, Ouchi H. Improvement of near wellbore permeability by methanol stimulation in a methane hydrate production well. Paper OTC 19433, Proceedings of the Offshore Technology Conference, Houston, Texas, 2008. Maxwell, J. C. 1879. On stresses in rarified gases arising from inequalities of temperature. Philosophical Transactions of the Royal Society of London, 231-256. Mehmani, A., Prodanović, M. and Javadpour, F. 2013. Multiscale, multiphysics network modeling of shale matrix gas flows. Transport in porous media, 99(2), pp.377-390.
SC
Mengal, S.A. and R.A. Wattenbarger. Accounting for adsorbed gas in shale gas reservoirs, SPE 141085, presented at the SPE Middle East Oil and Gas Show and Conference held in Manama, Bahrain, 25–28 September 2011.
M AN U
Minkoff, S.E., Stone, C.M., Bryant, S., Peszynska, M., Wheeler, M.F., 2003. Coupled fluid flow and geomechanical deformation modeling. Journal of Petroleum Sciences and Engineering 38, 37-56. MIT (Massachusetts Institute of Technology), 2011. The Future of Natural Gas – An Interdisciplinary MIT Study, (Updated from 2010 version). Moridis, G. J., Kowalsky, M. B., and Pruess, K. 2007. Depressurization-induced gas production from class-1 hydrate deposits. SPE Reservoir Evaluation and Engineering, 10(05), 458-481.
TE D
Moridis, G. J., Reagan, M. T., Kim, S. J., Seol, Y., and Zhang, K. 2009. Evaluation of the gas production potential of marine hydrate deposits in the Ulleung Basin of the Korean East Sea. SPE Journal, 14(04), 759-781. Moridis, G.J., T.A. Blasingame, and C.M. Freeman. Analysis of Mechanisms of Flow in Fractured Tight-Gas and Shale-Gas Reservoirs. SPE 139250, presented at the SPE Latin American & Caribbean Petroleum Engineering Conference held in Lima, Peru, 1–3 December 2010.
EP
Moridis, G. J., Silpngarmlert, S., Reagan, M. T., Collett, T., and Zhang, K. 2011. Gas production from a cold, stratigraphically-bounded gas hydrate deposit at the Mount Elbert Gas Hydrate Stratigraphic Test Well, Alaska North Slope: Implications of uncertainties. Marine and Petroleum Geology, 28(2), 517-534.
AC C
Nelson, P. H. 2009. Pore-throat sizes in sandstones, tight sandstones, and shales. AAPG Bulletin, 93(3), 329-340. Nguyen T. S. Description of the computer code FRACON. In: Stephansson, O., Jing, L., Tsang C-F, editors. Coupled thermo-hydro-mechanical processes of fractured media, vol. 79. Elsevier: Developments in Geotechnical Engineering, 1996. p. 539-44. Noorishad J, Tsang C-F, Witherspoon PA. 1984. Coupled thermal-hydraulic-mechanical phenomena in saturated fractured porous rocks: numerical approach. J Geophys Res., 89:10365-73. Noorishad J, Tsang C-F. ROCMAS-simulator: a thermohydro-mechanical computer code. In: Stephansson O, Jing L, Tsang C-F, editors. Coupled Thermo-hydro-mechanical processes of fractured media, vol. 79. Elsevier: Developments in Geotechnical Engineering, 1996. p. 551-8.
37
ACCEPTED MANUSCRIPT
O’Hare, L., Lockerby, D. A., Reese, J. M., and Emerson, D. R. 2007. Near-wall effects in rarefied gas microflows: some modern hydrodynamic approaches. International Journal of Heat and Fluid Flow, 28(1), 3743.
RI PT
Ohnishi, Y. and A. Kobayashi, THAMES. In: Stephanson, O., Jing, L. and Tsang, C. F. editors. Coupled thermohydro-mechanical processes of fractured media: Mathematical and Experimental Studies, vol. 79. Elsevier: Developments in Geotechnical Engineering, 1997, p. 545-549. Ohwada, T. and Xu, K. 2004. The kinetic scheme for the full-Burnett equations. Journal of Computational Physics, 201(1), 315-332. Okabe, H. and Blunt, M. J. 2004. Prediction of permeability for porous media reconstructed using multiple-point statistics. Physical Review E, 70(6), p.066135.
SC
Olivella, S., Carrera, J., Gens, A., Alonso, E.E., 1994. Nonisothermal multiphase flow of brine and gas through saline media. Transport in Porous Media 15, 271-293.
M AN U
Øren, P.E. and Bakke, S. 2002. Process based reconstruction of sandstones and prediction of transport properties. Transport in Porous Media, 46(2-3), pp.311-343. Pan, L. S., Liu, G. R., and Lam, K. Y. 1999. Determination of slip coefficient for rarefied gas flows using direct simulation Monte Carlo. Journal of Micromechanics and Microengineering, 9(1), 89. Patil, V. E., Meeuwissen, J., van den Broeke, L. J., and Keurentjes, J. T. 2006. Permeation of supercritical fluids across polymeric and inorganic membranes. The Journal of Supercritical Fluids, 37(3), pp.367-374. Pollard, W. G. and Present, R. D. 1948. On gaseous self-diffusion in long capillary tubes. Physical Review, 73(7), 762.
TE D
Present, R. D. 1958. Kinatic Theory of Gases.
Rodgers, D. P. June, 1985. Improvements in multiprocessor system design. In ACM SIGARCH Computer Architecture News (Vol. 13, No. 3, pp. 225-231). IEEE Computer Society Press.
EP
Rutqvist, J., 2011. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations.Computers & Geosciences, 37(6), pp.739-750.
AC C
Rutqvist J., Börgesson L., Chijimatsu M., Kobayashi A., Nguyen T. S., Jing L., Noorishad J., and Tsang C.-F. 2001. Thermohydromechanics of partially saturated geological media - Governing equations and formulation of four finite element models. Int. J. Rock Mech. & Min. Sci. 38: 105-127. Rutqvist, J., Zheng, L., Chen, F, Liu H.-H, and Birkholzer, J. 2014. Modeling of Coupled ThermoHydroMechanical Processes with Links to Geochemistry Associated with Bentonite-Backfilled Repository Tunnels in Clay Formations. Rock Mechanics and Rock Engineering, 47: 167-186. Shabro, V., C. Torres-Verdín, and F. Javadpour, Numerical simulation of shale-gas production: from pore-scale modeling of slip-flow, knudsen diffusion, and langmuir desorption to reservoir modeling of compressible fluid, SPE 144355, presented at the SPE North American Unconventional Gas Conference and Exhibition held in The Woodlands, Texas, USA, 14–16 June 2011. Sakhaee-Pour, A. and S. L. Bryant. 2011. Gas permeability of shale. Presented at the SPE Annual Technical Conference and Exhibition held in Denver, Colorado, USA, 30 October–2 November. SPE 146944.
38
ACCEPTED MANUSCRIPT
Sakhaee-Pour, A. and Bryant, S., 2012. Gas permeability of shale. SPE Reservoir Evaluation and Engineering, 15 (04): 401-409. Sampath, K. and Keighin, C. W. 1982. Factors affecting gas slippage in tight sandstones of cretaceous age in the Uinta basin. Journal of Petroleum Technology, 34(11), 2-715.
RI PT
Schettler Jr, P.D. and Parmely, C.R., January, 1991. Contributions to total storage capacity in Devonian shales. In SPE Eastern Regional Meeting. Society of Petroleum Engineers. Settari, Antonin, and Dale A. Walters. 2001. Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction. SPE Journal 6(03): 334-342. Shan, X., Yuan, X., and Chen, H. 2006. Kinetic theory representation of hydrodynamics: A way beyond the Navier-Stokes equation. Journal of Fluid Mechanics 550: 413-441.
SC
Silin, D. and Patzek, T., 2006. Pore space morphology analysis using maximal inscribed spheres. Physica A: Statistical mechanics and its applications, 371(2), pp.336-360.
M AN U
Silin, D. and T. Kneafsey, Gas shale: from nanometer-scale observations to well modeling, CSUG/SPE 149489, presented at the Canadian Unconventional Resources Conference, Calgary, Alberta, Canada, 15-17 November, 2011. Silin, D., 2014. Digital rock studies of tight porous media.
Sisk, C., Diaz, E., Walls, J. et al. 2010. 3D visualization and classification of pore structure and pore filling in gas shales. Presented at the SPE Annual Technical Conference and Exhibition, Florence, 19-22 September. SPE-134582-MS.
TE D
Soeder, D. J. 1988. Porosity and permeability of Eastern Devonian Gas Shale. SPE Formation Evaluation, 3(01): 116-124. Soeder, D. J. and P. L. Randolph, 1987. Porosity, permeability, and pore structure of the tight mesaverde sandstone, Piceance Basin, Colorado, SPEFE 2(6):129–136. SPE 13134, Trans., AIME, 283, 1987.
EP
Sun, Q. Y., Zhang, R. L., Zhang, Y. Z., and Chen, X. Q. 2010. Pressure Transient Analysis of the Low Permeability Composite Reservoir with Threshold Pressure Gradient. In 2010 Asia-Pacific Power and Energy Engineering Conference, pp. 1-4. IEEE.
AC C
Sun, H., Chawathe, A., Hoteit, H., Shi, X., and Li, L. 2015. Understanding Shale Gas Flow Behavior Using Numerical Simulation. SPE Journal, 20(01): 142-154. Szalmas, L. 2007. Slip on curved boundaries in the lattice-Boltzmann model. International Journal of Modern Physics C 18 (1): 15-24. Tang, G., Zhang, Y., Gu, X. et al. 2008. Lattice Boltzmann modeling Knudsen layer effect in non-equilibrium flows. Europhysics Letters 83 (4): 40008. Taron, J., D. Elsworth and K. B. Min, 2009. Numerical simulation of thermal-hydrologicmechanical-chemical processes in deformable, fractured porous media, International Journal of Rock Mechanics Mining Sciences, 46, 842-854. Thomas, H. R., Sansom, M. B., 1995. Fully coupled analaysis of heat, moisture, and air transfer in unsaturated soil. Journal of Engineering Mechanics ASCE 121, 392-405.
39
ACCEPTED MANUSCRIPT
Tinni, A., Fathi, E., Agarwal, R., Sondergeld, C. H., Akkutlu, I. Y., and Rai, C. S. January, 2012. Shale permeability measurements on plugs and crushed samples. In SPE Canadian Unconventional Resources Conference. Society of Petroleum Engineers.
RI PT
Tran, David, Long Nghiem, and Lloyd Buchanan. 2009. Aspects of coupling between petroleum reservoir flow and geomechanics. In 43rd US Rock Mechanics Symposium & 4th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association. Wang, F. P., R. M. Reed, J. A. Jackson, and K. G. Jackson. Pore networks and fluid flow in gas shales, SPE 124253, presented at the 2009 SPE Annual Technical Conference and Exhibition held in New Orleans, Louisiana, USA, 4–7 October 2009.
SC
Wang, K., Liu, H., and Chen, Z. 2015. A scalable parallel black oil simulator on distributed memory parallel computers. Journal of Computational Physics, 301, 19-34. Wang, L., 2014. Simulation of slip flow and phase change in nanopores. PhD dissertation.
M AN U
Wang, S., Xiong, Y., Winterfeld, P., Zhang, K., and Wu, Y. S. 2014. Parallel Simulation of ThermalHydrological-Mechanic (THM) Processes in Geothermal Reservoirs. In Proceeding of 2014 Stanford Geothermal Workshop. Wang, S., Winterfeld, P. H., & Wu, Y. S. 2015. An Efficient Adaptive Nonlinearity Elimination Preconditioned Inexact Newton Method for Parallel Simulation of Thermal-Hydraulic-Mechanical Processes in Fractured Reservoirs. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
TE D
Wang, S., Huang, Z., Wu, Y. S., Winterfeld, P. H., and Zerpa, L. E. 2016. A semi-analytical correlation of thermal-hydraulic-mechanical behavior of fractures and its application to modeling reservoir scale cold water injection problems in enhanced geothermal reservoirs. Geothermics, 64, 81-95.
Wang, W., Kolditz, O., 2007. Object-oriented finite element analysis of thermo-hydro-mechanical (THM) problems in porous media. International Journal Numerical Methods Engineering 69, 162-201.
EP
Wheeler, J. A. and Smith, R. A. 1990. Reservoir simulation on a hypercube. SPE Reservoir Engineering, 5(04), 544-548.
AC C
Wu, K., Nunan, N., Crawford, J.W., Young, I.M. and Ritz, K., 2004. An efficient Markov chain model for the simulation of heterogeneous soil structure. Soil Science Society of America Journal, 68(2), pp.346-351. Wu, K., Van Dijke, M. I., Couples, G. D., Jiang, Z., Ma, J., Sorbie, K. S., Crawford, J., Young, I. and Zhang, X., 2006. 3D stochastic modelling of heterogeneous porous media–applications to reservoir rocks. Transport in Porous Media, 65(3), pp.443-467. Wu, L., and Bogy, D. B. 2003. New first and second order slip models for the compressible Reynolds equation. Journal of Tribology, 125(3): 558-561. Wu, L. 2008. A slip model for rarefied gas flows at arbitrary Knudsen number. Applied Physics Letters, 93(25): 253103. Wu, Y-S., Fakcharoenphol, P., Zhang, R. et al. 2010. Non-Darcy displacement in linear composite and radial flow porous media, in SPE EUROPEC/EAGE Annual Conference and Exhibition, Society of Petroleum Engineers. 40
ACCEPTED MANUSCRIPT
Wu, Y-S., Chen, Z., Kazemi, H., Yin, X., Pruess, K., Oldenburg, C., Winterfeld, P., Zhang, R. 2014. Simulation of Coupled Processes of Flow, Transport, and Storage of CO2 in Saline Aquifers, Tech. rep., Trustees of The Colorado School of Mines. Wu, Y.-S. 2016. Multiphase Fluid Flow in Porous and Fractured Reservoirs. Gulf Professional Publishing.
RI PT
Wu, Y.-S. and P. Fackahroenphol, A unified mathematical model for unconventional reservoir simulation, SPE 142884, submitted for the SPE EUROPEC Conference, 23-26 May 2011 in Vienna, Austria, 2011. Wu, Y.-S., Li, J., Ding, D., Wang, C., and Di, Y. 2014. A generalized framework model for the simulation of gas production in unconventional gas reservoirs. SPE Journal, 19(05), 845-857. Wu, Y.-S. and G. Qin, 2009. A general numerical approach for modeling multiphase flow and transport in fractured porous media,” Communications in Computational Physics, 6(1): 85-108.
SC
Wu, Y.-S. and C. Wang, Transient pressure testing analysis of gas wells in unconventional reservoirs, SPE-SAS312, presented at the 2012 SPE Annual Technical Symposium & Exhibition, Khobar, Saudi Arabia, 811 April 2012.
M AN U
Xiao, F., 2013, Pore-Scale Simulation Frameworks for Flow and Transport in Complex Porous Media. PhD dissertation. Xiong, Y., Fakcharoenphol, P., Winterfeld, P.H., Zhang, R. and Wu, Y-S. 2013. Coupled Geomechanical and Reactive Geochemical Model for Fluid and Heat flow: Application for Enhanced Geothermal Reservoir, in SPE Reservoir Characterization and Simulation Conference and Exhibition, 16-18 September, SPE, Abu Dhabi, UAE.
TE D
Yao, Y., Wu, Y. S., Zhang, R. 2012. The Transient Flow Analysis of Fluid in a Fractal, Double-porosity Reservoir. Transport in porous media, 94(1), 175-187. Yin, S., M. B. Dusseault, and L. Rothenburg. 2011. Coupled THMC modeling of CO2 injection by finite element methods. Journal of Petroleum Science and Engineering 80(1): 53-60. 90.
EP
Yin, X., Koch, D. L., and Verberg, R. 2006. Lattice-Boltzmann method for simulating spherical bubbles with no tangential stress boundary conditions, Physical Review E 73 (2): 026301.
AC C
Zeng, B., L. Cheng, and F. Hao, Experiment and mechanism analysis on threshold pressure gradient with different fluids, SPE 140678, Nigeria Annual International Conference and Exhibition, Tinapa - Calabar, Nigeria, 31 July - 7 August 2010. Zhang, H., Bai, B., Song, K. et al. 2012. Shale gas hydraulic flow unit identification based on SEM-FIB tomography. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 810 October. SPE-160143-MS. Zhang, K. N., Wu, Y.-S., Ding, C., Pruess, K., and Elmroth, E. January, 2001. Parallel computing techniques for large-scale reservoir simulation of multi-component and multiphase fluid flow. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. Zhang, P., Hu, L., Meegoda, J.N. and Gao, S., 2015. Micro/nano-pore network analysis of gas flow in shale matrix. Scientific reports, 5.
41
ACCEPTED MANUSCRIPT
Zhang, R, P. H. Winterfeld, X. Yin, Y. Xiong, and Y-S. Wu. 2015. Sequentially coupled THMC model for CO2 geological sequestration into a 2D heterogeneous saline aquifer. Journal of Natural Gas Science and Engineering 27: 579-615. Zhang, R., Lu, N., and Wu, Y. S., 2012a. Efficiency of a community-scale borehole thermal energy storage technique for solar thermal energy. Proceedings of the GeoCongress.
RI PT
Zhang, R., Yin, X., Winterfeld, P.H. and Wu, Y-S. 2012b. A Fully Coupled Model of Nonisothermal Multiphase Flow, Geomechanics, and Chemistry during CO2 Sequestration in Brine Aquifers, Proceedings of the TOUGH Symposium, pp.838-848.
SC
Zhang, R., Yin, X., Wu, Y-S. and Winterfeld, P.H. 2012c. A Fully Coupled Model of Nonisothermal Multiphase Flow, Solute Transport and Reactive Chemistry in Porous Media, in SPE Annual Technical Conference and Exhibition, 8-10 October, San Antonio, Texas, USA. Zhang, R. 2013. Numerical Simulation of Thermal Hydrological Mechanical Chemical Processes during CO2 Geological Sequestration, PhD thesis, Colorado School of Mines.
M AN U
Zhang, R., Wu, Y-S. and Fakcharoenphol, P. 2014. Non-Darcy Displacement in Linear Composite and Radial Aquifer during CO2 Sequestration, International Journal of Oil, Gas and Coal Technology, Vol. 7, No. 3, pp.244-262. Zhang, R. and Wu, Y.S., 2016a. Sequentially Coupled Model for Multiphase Flow, Mean Stress, Reactive Solute Transport with Kinetic chemical reactions for CO2 Geological Sequestration and EGS System, Journal of Porous Media, Vol. 19 (11), pp. 1101-1021. DOI: 10.1615/JPorMedia.v19.i11.60.
TE D
Zhang, R., Xiong, Y., Winterfeld, P. H., Yin, X. and Wu, Y.-S. 2016b. A Novel Computational Framework for Thermal - Hydrological - Mechanical - Chemical Processes of CO2 Geological Sequestration into a Layered Saline Aquifer and a Naturally Fractured Enhanced Geothermal System. Greenhouse Gas: Science and Technology, Vol. 6 (3), pp. 370-400. doi: 10.1002/ghg.1571.
EP
Zhang, R., Yin, X., Winterfeld, P.H. and Wu, Y-S. 2016c. A Fully Coupled Thermal - Hydrological - Mechanical - Chemical Model for CO2 Geological Sequestration, Journal of Natural Gas Science and Engineering, Vol. 28, pp.280-304. Zhang, R., Shan, X., and Chen, H. 2006. Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation. Physical Review E 74 (4): 046703.
AC C
Zhang, W.-M., Meng, G., and Wei, X. 2012. A review on slip models for gas microflows. Microfluidics and Nanofluidics, 13(6), 845-882. Zhang, Y., Qin, R., and Emerson, D. R. 2005. Lattice Boltzmann simulation of rarefied gas flow in microchannels. Physical Review E 71 (4): 047702. Zhao, X., 1991. An Experimental Study of Methane Diffusion in Coal Using Transient Approach, Ph.D. Dissertation, Department of Mining and Geological Engineering, The University of Arizona. Zhao, X., Rui, Z., Liao, X. and Zhang, R. 2015a. The Qualitative and Quantitative Fracture Evaluation Methodology in Shale Gas Reservoir, Journal of Natural Gas Science and Engineering, Vol. 27, pp.486495.
42
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Zhao, X., Rui, Z., Liao, X. and Zhang, R. 2015b. A Simulation Method for Modified Isochronal Well Testing to Determine Shale Gas Well Productivity, Journal of Natural Gas Science and Engineering, Vol. 27, pp.479-485.
43
ACCEPTED MANUSCRIPT
Highlights 1. Fluid flow and transport mechanisms along with geomechanics in shale and tight gas reservoirs are recapitulated.
RI PT
2. Advances in macroscopic THMC simulation techniques for shale and tight gas reservoirs are reviewed. 3. Pore-scale modeling methods for fluid flow and transport in unconventional reservoir rocks are illustrated and compared. 4. Challenges from unconventional reservoirs and perspectives on developing new simulation
AC C
EP
TE D
M AN U
SC
techniques are discussed.
S1875-5100(16)30887-3
DOI:
10.1016/j.jngse.2016.11.051
Reference:
JNGSE 1961
To appear in:
Journal of Natural Gas Science and Engineering
Received Date: 9 August 2016 Revised Date:
15 November 2016
Accepted Date: 24 November 2016
Please cite this article as: Wang, L., Wang, S., Zhang, R., Wang, C., Xiong, Y., Zheng, X., Li, S., Jin, K., Rui, Z., Review of multi-scale and multi-physical simulation technologies for shale and tight gas reservoirs, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.11.051. 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
Review of Multi-scale and Multi-physical Simulation Technologies for Shale and Tight Gas Reservoirs
Abstract
RI PT
Lei Wanga, Shihao Wanga, Ronglei Zhanga*, Cong Wanga, Yi Xionga, Xishen Zhengb, Shangru Lib, Kai Jinb, Zhenhua Ruic a Energy Modeling Group, Petroleum Engineering Department, Colorado School of Mines, USA b Jiangxi Unconventional Gas Institute, China c Independent Project Analysis, Inc., USA *Corresponding author: Ronglei Zhang, Email: [email protected]
Keyword:
TE D
M AN U
SC
This paper provides a comprehensive literature review on the simulation techniques being developed in recent years for describing unique flow behaviors in shale and tight gas reservoirs. The advances in modeling gas flow and transport mechanisms during the primary and enhanced gas recovery processes are reviewed in detail. The capabilities of reservoir simulation tools are discussed in terms of mathematical treatment (finite difference, finite element, explicit/implicit scheme, sequentially and fully coupled schemes), fluid flow characteristics (Darcy and non-Darcy flow, desorption, Klinkenberg effect and gas slip flow, transitional flow, Knudsen diffusion), reservoir rock properties (pore size distribution, fractures, geomechanics), multidisciplinary coupling scheme (THM, THC, THMC), modeling scale, and computational efficiency. For pore-scale modeling of gas flow and transport in unconventional reservoir rocks, the numerical methods for generating pore network models and the procedure for constructing real 3D digital rocks are first explained. Then, the network modeling methods used for simulating slip and transitional flows in pore networks are introduced and compared. After that, existing lattice-Boltzmann models developed for slip flow simulation are illustrated. Along with the explanation, pros and cons of the models for pore-scale modeling have been identified. Overall, dwelling on the concerns and challenges in shale and tight gas reservoir simulation, current status, progress and bottlenecks of numerical simulation technologies are discussed and perspectives on future development of unconventional gas reservoir simulators are proposed.
AC C
Highlights
EP
Shale and tight gas reservoir, reservoir simulation, pore-scale simulation, adsorption and desorption, Klinkenberg effect and slip flow, transitional flow, Knudsen diffusion, non-linear and non-Darcy flow, geomechanics, THMC
1. Fluid flow and transport mechanisms along with geomechanics in shale and tight gas reservoirs are recapitulated. 2. Advances in macroscopic THMC simulation techniques for shale and tight gas reservoirs are reviewed. 3. Pore-scale modeling methods for fluid flow and transport in unconventional reservoir rocks are illustrated and compared. 4. Challenges from unconventional reservoirs and perspectives on developing new simulation techniques are discussed.
1
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
1. Introduction Unconventional gas reservoirs, such as shale and tight gas reservoirs, are rich resources of natural gas. In unconventional formations, huge amount of gas is trapped in the location where it is generated, without the primary migration process. Featured by extremely low permeability, shale and tight gas reservoirs are different from each other, even with similar porosity, pore pressure, and in-situ stress conditions. As shown in Figure 1 (Wang et al. 2009), permeability values of organic-rich shales are in the range of sub-nanodarcy to tens of microdarcy, and permeability of deep organic-lean mudrocks can be even less. Some correlations can also be observed between the permeability values and the porosity of shales.
Figure 1. The relationships of porosity and permeability for shale gas reservoirs (FWB-Fort Worth Basin, WV #6 well-Marcellus Shale well in West Virginia) in North America obtained from core analysis (Wang et al. 2009).
AC C
EP
The diameter of pore throat in tight gas formations is usually between 0.01µm and 1µm (Holditch 2006, Nelson 2009), while the pore throat diameter of shale formation varies from tens of nanometers to as low as a few nanometers (Nelson 2009). Recent studies even revealed sub-nanometer pores inside shale formations that are beyond the scope of scanning electron microscope (SEM) based on CO2 and N2 adsorption experiments (Clarkson et al. 2013) and atom-scale reconstruction technique (Falk et al. 2015). Due to these extremely small pores, the hydrocarbon transport mechanisms in unconventional reservoirs turn out to be very different from those in conventional reservoirs. Fluid flow behavior in shale and tight gas reservoirs is featured by multiscale single-phase (gas) and/or multiphase fluid (gas, gas condensate and/or brine) flow and transport in ultra-low permeability, highly heterogeneous porous/fractured, and stress-sensitive rocks. Compared with flow in conventional reservoirs, fluid flows in unconventional reservoirs with extremely low permeability are highly non-linear and complex processes, including non-linear adsorption/desorption, non-Darcy flow in the entire range from high flow rate to low flow rate, strong rock-fluid interaction, and rock/organic matter deformation within nano-pores or micro-fractures. Therefore, quantifying flow in unconventional gas reservoirs has become a challenging issue for conventional REV (Representative Elementary Volume)-based Darcy’s law approaches, which may not be suitable and thus must be improved in order to handle these realistic problems in shale and tight gas reservoirs. 2
ACCEPTED MANUSCRIPT
SC
RI PT
2. Gas Flow and Transport Mechanisms Due to the existence of organic matter with strong adsorption capacity in shale gas reservoirs, natural gas components experience desorption from its original adsorbed state as reservoir pressure depletes. In addition, flow in micro- and nano-pores in unconventional reservoirs falls in different regimes from those in conventional reservoirs. Besides, decrease of pore pressure would induce compaction of formations and reduction of the effective pore sizes, resulting in permeability decrease. Also organic matter may shrink due to desorption of adsorbed gas, resulting in dimensional changes of the pores they are residing in. Since shale and tight gas reservoirs are generally completed by hydraulic fracturing, disparities between flow behavior in matrix and fracture networks should be considered as well. All these flow and transport mechanisms are described in the following sections. 2.1 Gas Adsorption and Desorption
AC C
EP
TE D
M AN U
Natural gas presents as free gas and adsorbed gas in shale reservoirs. Significant amount of gas can be adsorbed onto the surface of the organic matter and clays under reservoir conditions (Mengal and Wattenbarger, 2011, EIA, 2011a). In shale matrix, methane molecules are mainly adsorbed onto the organic matter, i.e. kerogen, which is generally measured in terms of total organic carbon (TOC). Adsorbed gas can provide significant portion of gas in place (GIP) as well as affect ultimate recovery rate. In Barnett Shale, when pressure is below 2,000 psi, adsorbed gas accounts for 50 to 80 % of GIP. As pressure goes above 2,000 psi, although the adsorbed gas content tends to level out, it is still significant in percentage, about 30 to 50 % (Mengal and Wattenbarger, 2011). Figure 2 is a case study of Barnett Shale showing typical percentage of adsorbed gas to total gas, and its linear dependence on TOC content (Wang et al. 2009). Gas adsorption capacity also changes with in-situ reservoir pressure. As pressure declines under continuous reservoir depletion, more adsorbed gas is dissociated from adsorbed state to free gas, leading to sustainable gas flow.
Figure 2. Adsorbed and total gas content versus TOC of Barnett Shale (Wang et al. 2009).
2.2 Gas Rarefaction Effects and Flow Regimes Gas rarefaction effects in unconventional reservoirs, including velocity slippage and gas molecular diffusion, have been demonstrated by core flooding experiments (Luffel et al. 1993, Civan et al. 2011, Tinni et al. 2012, Zhao et al., 2015a and 2015b). These two transport phenomena are essentially induced by the gas molecule and 3
ACCEPTED MANUSCRIPT
RI PT
pore wall interactions. In gas sate, molecules keep moving in random directions with thermal velocity. The average distance that one gas molecule can move freely before it collides with any other molecules is defined as the mean free path λ of the gas, which can be calculated as (Bird 2003) k BT (1) λ= 2π Pδ 2 In the above formulation for a single component ideal gas, denotes the Boltzmann constant, T denotes the temperature, P denotes the pressure of the system, δ denotes the collision diameter of the gas. For methane, δ is about 0.40 nm. When the pore diameter d becomes narrow enough to be comparable with the mean free path of the gas molecules in it, the gas molecules will frequently collide with the walls of the flow channel and the molecule-wall interaction cannot be ignored. The gas inside such narrow flow channels becomes rarefied and possesses certain unique transport characteristics. The gas rarefication is quantified by Knudsen number as (Bird 2003)
λ
SC
Kn =
d
(2)
M AN U
where d is the characteristic length, e.g. the diameter of cylindrical pores. The average pore size can either be directly measured from experiments or be estimated from permeability
k ∞ , tortuosity τ, and porosity ϕ of the rock. For example, Civan (2010) used the following correlation to estimate the mean pore radius of shale.
r = 2 2τ
k∞
φ
(3)
AC C
EP
TE D
When a hard-sphere gas molecule collides with the solid wall, its reflection will be either specular or diffuse, depending on the roughness of the wall as well as the molecule type. On one hand, after specular reflection, as shown by the blue arrow in Figure 3, the molecule maintains its previous horizontal velocity component, changing only its moving direction. On the other hand, after diffuse reflection, the molecule will reflect diffusely into a random direction, changing its horizontal and vertical velocity components, as shown by the red arrows in Figure 3. The maintenance of horizontal velocity component results in the ‘slippage’ of gas molecules in the near-wall region of the pores. Meanwhile, the reflected molecules themselves contribute to the diffusion effect in pores with molecular density differences. Therefore, the apparent permeability enhancement is a result of the combined velocity slippage and diffusion effect, both of which come from the molecule-wall interaction. The ratio of the diffuse reflection to the total amount of wall collision is defined as tangential momentum accommodation coefficient (TMAC) σ, ranging from 0 to 1. The velocity slippage as well as the diffusion effect results in the permeability enhancement that the apparent gas permeability of shale and tight gas rocks is larger than their absolute permeability . The relationship between and can be quantified as
ka = f ( Kn ) k∞
(4)
Where f(Kn) is a permeability multiplier, depending on Knudsen number and gas molecule type (Zhang, W.-M. et al. 2012).
4
RI PT
ACCEPTED MANUSCRIPT
Figure 3. Conceptual sketch showing specular reflection and diffuse reflection.
M AN U
SC
To describe gas flow and transport in micro- and nanoscale geometries, the governing equation is chosen upon flow regimes. Table 1 shows the gas flow regimes and corresponding governing equations along with boundary conditions. When is between 0.001 and 0.1, gas flow is in the regime of slip flow, slip boundary condition should be incorporated into Navier-Stokes equation or lattice-Boltzmann method (LBM) to take into account the slippage on the gas solid interface. When is in a higher range of 0.1 and 10, gas flow enters the transitional regime, where neither Navier-Stokes equation nor lattice-Boltzmann model is applicable any more. Then, Burnett equation based on higher order moments of Boltzmann equation should be solved or numerical method of direct simulation Monte Carlo (DSMC) should be used to represent the fluid flow behavior. As goes beyond 10, gas stream is considered as free molecules, and Molecular Dynamics (MD) must be adopted to capture the physics controlling the gas flow.
Continuum model Kn Flow regimes Note
Boltzmann equation
Collisionless Boltzmann equation
Lattice-Boltzmann method
DSMC
Molecular Dynamics
Navier-Stokes equation
Burnett equation
×
(0, 0.001) No slip
(0.001, 0.1) Slip
EP
Particle or molecule model
TE D
Table 1. Knudsen number and flow regimes with applicable mathematical models
(0.1, 10) Transitional Adapted from Bird, 2003.
>10 Free molecular
AC C
Many petrophysical measurements have demonstrated that pore throat radius of shale reservoirs falls in the nanometer range. For instance, Figure 4 shows the pore throat radius distribution of a typical shale formation measured by mercury injection into both crushed and plug core samples (Honarpour et al. 2012). It can be seen that the measurements of the pore size distribution are relatively consistent with each other in that the majority of the pores lie in the range of 1 to 100 nanometers by volumetric fraction. Considering that the major component of unconventional natural gas is methane, more than 80% in general (Schettler Jr. and Parmely 1989), the mean free path of methane gas molecules under reservoir conditions can be estimated from the equation above for ideal gas. Setting the unconventional gas reservoir at a typical temperature of 353.15K or 176ºF (Gong et al. 2013), when the reservoir pressure is 1,000 psi, in Figure 4, the pores bounded by the blue dashed lines, from ~9 to ~900 nanometers, host the slip flow regime defined by 0.001< <0.1. Given the reservoir pressures of 2,000 psi and 3,500 psi, the slip flow regime shifts to smaller pores, as denoted by the pore size ranges between the green dotted lines and red lines, respectively. To the left of these pore size ranges, 5
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
transitional flow occurs, while to the right, it is the familiar viscous flow regime that are accounted for by NavierStokes equation, from which Darcy’s law is derived for fully developed Stokes flow. It is worth noting that under the reservoir temperature and pressure conditions mentioned above, methane is in the supercritical state ( =190.6K, =668.6psi). Experiments for fluid flow through nano-pores demonstrated that diffusive flow mechanisms can be negligible when the fluid is in the supercritical state (Patil et al. 2006).
Figure 4. Measured pore throat radius distribution of a shale reservoir and corresponding flow regimes under different reservoir conditions.
AC C
EP
However, it should be noticed that flow regime is merely a classification, whereas the gas flow behavior is continuously changing with respect to Kn. Moreover, the pore structure in unconventional rocks is essentially irregular, pore sizes are continuous, and thus more than one flow regimes may exist simultaneously in the rock matrix under certain reservoir conditions. Therefore, it is of great importance to develop a unified model that can describe more than one flow regimes in order to accurately capture the commingled gas flow behavior in unconventional reservoir rocks. Figure 5 shows a conceptual model of gas flow along the x direction through a channel, in which a horizontal imaginary plain H and a vertical imaginary plain V are randomly selected for molecule movement analysis. Molecules of Type ① and ③ impinge into H after hitting the upper and lower wall, respectively. Molecules of Type ② and ④ impinge into H without any impact of the wall. The difference between the horizontal flux momentum carried by Type ① and ② molecules and that carried by Type ③ and ④ molecules generates shear stress on plain H. Because of the wall impact on Type ① and ③ molecules, their resultant shear stress along plain H is larger than the shear stress on an imaginary plain that is located inside an infinitely wide flow channel. This mechanism results in the velocity slippage in small channels. The wall impact is a function of the location along the z direction between the two walls. The region that is closer to the wall has high wall impact and the region that is closer to the middle of the flow channel has low or no wall impact. Similarly, the unequal
6
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
molecules passing through the plain V reduces the number density difference across it, resulting in the diffusion effect. The magnitude of the wall impact on diffusion also depends on the pore wall proximity.
Figure 5. Conceptual model of reflected molecules and bulk molecules.
AC C
EP
TE D
Based on the above analysis and for engineering purpose, the flux inside the micro-channel can be approximately divided into two regions: the Knudsen layer and the bulk flow, as shown in Figure 6. The flow inside regions that are close to the upper and lower walls are assumed to be primarily affected by the moleculewall interaction. The regions are therefore called the Knudsen layers. The thickness of the Knudsen layer is usually taken as 2/3 of the mean free path, the derivation of which can be found in (Present, 1958). The flow beyond the Knudsen layer is assumed to be not affected by the wall and is called the bulk flow, which can still be modeled by the Navier-Stokes equation. The transition from the Knudsen layer to bulk layer is not as distinct as shown in Figure 6. Clearly, more advanced techniques are required for accurate simulation of the whole flow field, as discussed in the following sections.
Figure 6. Conceptual velocity profile within a micro-channel. The yellow curve represents the bulk flow part where the wall barely impacts the flow. The ‘straight lines’ inside the Knudsen layer indicate the average velocity of the Knudsen layer. Real flux velocity of gas molecules inside the Knudsen layer is not constant.
7
ACCEPTED MANUSCRIPT
2.3 Klinkenberg Effect
TE D
M AN U
SC
RI PT
In last section, our theoretical analyses clarified that in slip, transitional, and free molecular flow regimes, the apparent gas permeability of rock matrix is higher than the intrinsic permeability. To account for this permeability enhancement, Klinkenberg factor was proposed to lump all the favorable factors together into the Klinkenberg’s correlation for reservoir simulation (Klinkernberg 1941). Practically in conventional gas reservoir engineering, the Klinkenberg or gas-slippage effect is ignored, except when analyzing pressure responses or flow near gas production wells at very low pressure in some cases. In shale and tight gas reservoirs, however, the Klinkenberg effect is expected to be significant, because of the nano-pores in such rocks, even under high pressure conditions. Figure 7 presents Soeder’s data (1988) that were measured under a confining pressure of 3,000 psi with a core sample acquired at 7,449 ft of the West Virgina No. 6 well. It can be seen that the gas permeability is 20 µD under pore pressure of 1,000 psi, and the gas permeability increases to 55 µD under pore pressure of 80 psi, owing to the increases of the mean free path of gas molecules and the enhanced slippage by pore pressure decrease.
Figure 7. The relationship between pore pressure and gas permeability of Marcellus Shale based on core measurements (Soeder 1988).
AC C
EP
In addition, Bryant and Sakhaee-Pour (2011) employed network modeling to investigate the impact of decreasing reservoir pressure on gas permeability. They calculated the ratio of the gas permeability (kg2,in-situ) at low pressures to the permeability (kg1,in-situ) at the initial reservoir pressure of 28 MPa as a function of reservoir pressure (Figure 8) by considering the impact of pressure on both adsorption layer thickness and gas slippage. The network modeling results indicated that the permeability during the initial production stage is strikingly lower than that during late period of production under lower reservoir pressure. It was claimed that the diminishing adsorbed layer is responsible for this impact. Moreover, the relationship between the permeability and the pressure is found to be non-linear. They concluded that unlike the conventional reservoirs, the gas permeability of shale highly depends on the reservoir pressure. These implications are useful for improving reservoir modelling and ultimate recovery prediction and estimation.
8
SC
RI PT
ACCEPTED MANUSCRIPT
2.4 Non-Linear and Non-Darcy Flow
M AN U
Figure 8. The ratio of gas permeability at low pressures, kg2,in-situ, to that at the initial pressure (P = 28 MPa), kg1,insitu, increases as production depletes and pressure drops (Sakhaee-Pour and Bryant, 2011).
EP
TE D
The dynamics of unconventional gas reservoirs are featured by non-linear multiphase flow in ultra-low permeability formation with multiscale heterogeneity, i.e. non-Darcy flow and strong rock-fluid interaction coexisting inside nano-pores and micro-fractures. Hence traditional REV-based Darcy’s law model may not be practical to characterize such flow dynamics. Take the traditional double-porosity model as an example, transient flow regime between fracture and matrix lasts long time period and it also needs long time to reach the pseudo steady-state flow regime. Therefore, it is not applicable to simulate the fluid flow in the fracture-matrix system of unconventional reservoirs. Two comprehensive literature reviews on flow mechanisms in unconventional gas reservoirs performed by Blasingame (2008) and Moridis et al. (2010) pointed out that high-velocity flow might be dominant during the production of shale gas. The non-laminar/non-Darcy flow definition regarding high velocity gas flow cannot be represented by traditional Darcy's law, but Forchheimer’s equation as shown below is probably sufficient for many applications.
−
∂p µ v = + βρ v 2 ∂x k
(5)
AC C
Where p is pressure, µ is viscosity, k is permeability, ρ is density, β is the non-Darcy flow coefficient (Wu 2016). Gas adsorption and desorption on organic matter and solid grains as well as diffusion of gas molecules into flow streams from low permeability rock matrix also contribute to the non-linear flow characteristics in shale and tight gas reservoirs. 2.5 Effect of Geomechanics
In conventional reservoirs, effect of geomechanics on rock deformation and permeability is generally small and is ignored mostly in practice. However, in unconventional reservoirs, such compaction effect can be significant and have a remarkable impact on both fracture and matrix permeability, which has to be taken into account during reservoir simulation. The impact of confining pressure on permeability is owing to the reduction of porosity. Figure 9 (Wang et al. 2009) shows that permeability of the Marcellus Shale is dependent on reservoir 9
ACCEPTED MANUSCRIPT
M AN U
SC
RI PT
pressure. It decreases with increasing effective stress, and the same trends can be seen for other shale reservoirs (Bustin et al. 2008). It is obvious that, as the confining pressure increases, the magnitude of permeability reduction in shale is significantly bigger than those of conventional sandstone or carbonate reservoir.
Figure 9. Effect of effective stress on shale gas permeability (Wang et al. 2009).
AC C
EP
TE D
Holditch (2006) compared air permeability values measured at both surface condition and original overburden pressure, as shown in Figure 10, to illustrate the effect of net overburden pressure or total stress on air permeability of Travis Peak core samples drilled from two wells in eastern Texas. For the cores with high permeability of more than 10 md, the permeability at the initial overburden pressure turns out to be relatively smaller than that of the same unstressed core plugs. Nonetheless, the impact of net overburden pressure increases significantly as the permeability of the cores decreases. In terms of the core plugs with unstressed permeability of 0.01 md, the permeability at net overburden pressure becomes ten times lower (0.001 md). The low permeability rocks with smaller pore sizes are much more sensitive to overburden stress than the rock of conventional reservoirs (Holditch 2006). In order to better understand the characteristics of the tight gas formations, such as gas permeability, capillary pressure, formation factor, resistivity index, acoustic velocity, and the rock mechanical properties, special core analyses must be conducted (Soeder and Randolph 1987).
Figure 10. Air permeability values at net overburden pressure versus those at ambient pressure for Travis Peak cores (Holditch 2006) 10
ACCEPTED MANUSCRIPT
TE D
M AN U
SC
RI PT
The impact of stress changes on fracture permeability in shale is expected to be significant, in sense that the fracture of 2D nature cannot support much stress load in the normal direction of fracture planes. Cipolla (2009) summarized some laboratory studies by approximating the flow conductivity of partially propped and unpropped fractures as a function of closure stress, as shown in Figure 11. The black curve shows the flow conductivity of an unpropped fracture in which the faces of the fracture planes are aligned upon closing. If this is true, very low conductivity would remain for the unpropped fractures under closure stresses over 2,800 psi, which are common closure stress in most of the shale reservoirs. Nonetheless, if the fracture planes are displaced from each other owing to shear offsets or are partially propped, the fracture conductivity as indicated by the blue curve will be enhanced significantly. Furthermore, the orange curve shows that the fracture conductivity can be retained by high-strength proppants, e.g. sintered bauxite, with increasing closure stress. Since the proppant concentration is much less than 0.1 lbm/ft2 in typical fracturing treatments of shale reservoirs, improvement of flow conductivity by partially propping would not be so prominent due to extra proppant crushing. However, shear offsets between fracture planes could happen, thus enhancing unpropped fracture conductivity to several md-ft in unconventional gas reservoirs.
EP
Figure 11. Impacts of closure stress (effective stress) on conductivity of unpropped and partially propped fractures (Cipolla 2009).
AC C
The effect of closure stress on propped fracture conductivity has been investigated for tens of years, however, the flow conductivity of unpropped fractures has not been well studied until recently. It is possible that some of the fracture networks created by hydraulic fracturing could be unpropped or only partially propped (Cipolla et al. 2008 and 2010). As a result, the magnitude of flow conductivity in unpropped or partially propped fracture could restrict the gas production of shale gas reservoirs (Cipolla et al. 2008). The impact of closure stress on unpropped or partially propped fracture conductivity for rocks with larger Young’s modulus was also investigated by experiments. Figure 12 gives the relationship between stress and unpropped fracture conductivity for Barnett Shale under different Young’s moduli. Fracture conductivity dramatically decreases with increasing closure stress and decreasing Young’s modulus (Cipolla et al. 2010), as can be tens of folds in magnitude.
11
ACCEPTED MANUSCRIPT
2
10
6
E=5x10 psi 6
E=4x10 psi 6
1
RI PT
10
0
10
-1
10
0
1000
2000 3000 Stress, psi
SC
Conductivity, mD-ft
E=2x10 psi
4000
5000
M AN U
Figure 12. The relationship between closure stress and flow conductivity of unpropped fracture under different Young’s modulus (Cipolla et al. 2010)
AC C
EP
TE D
As another example, Figure 13 estimates the impact of closure stress on flow conductivity of unpropped fracture in Marcellus shale under Young’s modulus of 2×106 psi based on previous studies (Cipolla et al. 2008) and the history matching results of Barnett Shale (Cipolla et al. 2010). The original flow conductivity of the fracture network tends to be around 2 md-ft. It decreases to 0.02 md-ft as the pressure inside the fractures approaches the flowing bottom hole pressure of 500 psi.
Figure 13. Effect of closure stress on unpropped fracture conductivity in Marcellus Shale (Cipolla et al. 2010) Gas recovery from shale and tight reservoirs is only about 10-30% of GIP, which is much lower than conventional gas reservoirs. Therefore, to promote the development of new technologies or approaches for enhanced gas recovery (EGR), it is necessary to clarify the fundamental mechanisms behind the coupled physical processes occurring in shale and tight gas reservoirs, including slow flow processes, poor connectivity, large immobile zones, high heterogeneity, stress sensitivity of rock, and slow decrease in long-term production rate. However, in-depth laboratory or theoretical studies on gas flow and transport in unconventional reservoirs are still lacking, and how these complicated mechanisms impact gas flow, production, and the ultimate gas recovery rate is still ambiguous. Future research topics in this area may aim at the investigation of surface diffusion inside the 12
ACCEPTED MANUSCRIPT
adsorption layer, the effect of thermal creep on the solid wall, and the effects of tortuosity of porous media on gas flows.
M AN U
SC
RI PT
3. Gas Flow and Transport Modeling Because of lacking quantitative correlations for corresponding gas flow and transport mechanisms inside unconventional gas reservoirs, there are no effective reservoir simulators or efficient modeling approaches available for optimizing shale or tight gas reservoir development (MIT 2011). Traditionally, the modeling approaches for coalbed methane formations have been used for modeling desorption and diffusion effect in shale gas reservoirs (e.g., King and Ertekin 1991, Zhao 1991). Recently, many analyses and discussions have been developed to characterize fluid flow behavior in shale and tight gas reservoirs (e.g. Javadpour et al. 2007, Blasingame 2008, Javadpour 2009, Wu et al. 2009, Cipolla et al. 2010, Moridis et al. 2010, Sun et al., 2010, Wu and Fackahroenphol 2011, Sakhaee-Pour and Bryant 2011, Shabro et al. 2011, Wu et al. 2012, Wu and Wang 2012, Yao et al., 2012, Li et al., 2016). In particular, Blasingame (2008) comprehensively discussed the fundamental technologies in geology, petrophysics, and reservoir engineering that can be applied to evaluate well/reservoir performance and characterize flow behavior in gas reservoirs with low permeability. In this section, relevant findings and conclusions regarding flow mechanisms and physical processes in shale and tight gas reservoirs are discussed. 3.1 Langmuir Isotherm Model
TE D
Many studies (e.g. Mengal and Wattenbarger 2011, Silin and Kneafsey 2011, EIA 2011b, Wu et al. 2012, Li et al., 2016) showed that the mass of adsorbed gas in formation can be correlated to gas pressure according to the Langmuir isotherm (1918). Adsorbed gas content and adsorption isotherm can be measured in laboratory using core samples. By definition, adsorbed gas content is the amount of gas being adsorbed onto the pore surface of the reservoir rock, and adsorption isotherm represents the capacity of the reservoir rock to adsorb gas under a certain pressure while keeping the temperature constant. Then, they are employed to derive the relationship between pressure and gas storage capacity of the reservoir rock.
VE = VL
AC C
EP
In Langmuir isotherm, adsorbed gas volume, VE, is generally measured in cubic feet of gas per ton of net shale (Mengal and Wattenbarger 2011, EIA 2011b),
P P+ PL
(6)
where VL is the Langmuir volume in scf/ton at infinite pressure, P is reservoir gas pressure, and PL is the Langmuir pressure, at which 50% of the gas is desorbed. In general, Langmuir volume is a function of the TOC content and thermal maturity of shale. An example of gas adsorption isotherms measured from low TOC and high TOC Marcellus Shale samples are shown in Figure 14 (EIA 2011b).
13
ACCEPTED MANUSCRIPT
Adsorbed Gas Content: Higher TOC
SC
RI PT
Adsorbed Gas Content: Lower TOC
Figure 14. Marcellus Shale methane adsorption isotherms (EIA, 2011b)
AC C
EP
TE D
M AN U
As discussed above, adsorbed gas contributes significantly to GIP in shale reservoirs and should be accounted for in flow and production performance assessment. Resource concentration of a prospective area in shale basin is evaluated as the combination of free gas in pores and adsorbed gas that can be released at a finite decrease in pressure. Figure 15 illustrates the relative contributions of free gas (porosity) and adsorbed gas (sorbed) to the total GIP as a function of pressure, as adapted from EIA report (2011b). It indicates that gas reserves in shales are dominated by adsorption under low pressure.
Figure 15. Free, adsorbed, and total gas content versus pressure (EIA, 2011b).
3.2 Slip Boundary Condition
To accurately simulate the velocity slippage, both the slip boundary condition on the wall and the velocity distribution inside the Knudsen layer should be considered. The first-order slip boundary condition for gas flow through microscale space is proposed by Maxwell (1879)
uslip = −
2 −σ
σ
14
∂u Kn ∂n s
(7)
ACCEPTED MANUSCRIPT
where σ is the TMAC. uslip is the velocity on the wall, and ( ∂u ∂n ) s is the velocity gradient along the surface normal. As can be seen, the entire flux is accelerated by the velocity slippage on the boundary, which is proportional to the velocity gradient.
RI PT
Since Maxwell’s pioneering work, many slip boundary conditions have been brought up. In general, it is believed that the slip boundary condition depends on Kn and the velocity condition at the wall, as can be expressed in the following form
∂ 2u ∂u uslip = C1 Kn − C2 Kn 2 2 ... ∂n s ∂n s
(8)
M AN U
SC
where C1 and C2 could be either constants or functions of Kn . If C2 is zero, then the slippage velocity reduces to the first-order approximation. The slip velocity can be derived from kinetic theory (Wu and Bogy 2003, Wu 2008), DSMC (direct simulation Monte Carlo) method (Pan et al. 1999), Lattice Boltzmann Method (LBM), and Burnett or Super-Burnett equation with Chapman-Enskog expansion (Loyalka et al. 1975, Lockerby and Reese 2008). Zhang et al. (2012) conducted a thorough literature review of all the existing slip boundary models and their comparison of the predicted values with experimental results concluded that Wu’s model (2008) derived from the gas kinetic theory is the most accurate among all non-empirical correlations. In Wu’s model (2008) 2 2 3 − σ f 3 3 3 (1 − f ) 1 2 , C2 = f 4 + C1 = 1 − f 2 ) − 2 ( 3 σ 2 Kn 4 Kn
where f = min [1 Kn ,1] .
(9)
AC C
EP
TE D
To calculate the velocity profile in the Knudsen layer, there are mainly three models available, namely the wall-function model (O’Hare et al. 2007), the higher-order continuum model (Lockerby et al. 2005) and the power-law model (Aidun and Clausen 2010). The wall-function model is a phenomenological approach. The higher-order continuum model is developed based on flow governing equations, e.g. the Burnett, BGK–Burnett, and super-Burnett equations. According to the power-law model, the velocity inside the Knudsen layer is defined as a power function of distance away from the wall. By extrapolating the results from LBM and DSMC method, Lilley and Sader (2007, 2008) further promoted the power-law model. However, the singularity of the velocity profile predicted by power law is unphysical and brings difficulty in calculating the slip boundary. Fichman and Hetsroni (2005) derived the velocity profile inside the Knudsen layer by expanding the velocity field and calculating the shear stress at each location. Maxwell’s first-order slip boundary condition was adopted on the boundary. However, there are some mistakes in the calculation of the viscosity in the diffuse reflection case. All in all, a unified model that is capable of calculating both the slippage boundary condition and the velocity profile inside the Knudsen layer is still desirable. Theoretical analysis often assumes the flow space as an ideal straight channel or tube. Although it can provide us insights into the physical problems, in view of the complex pore structures in unconventional reservoirs, it cannot fully represent or replace experimental study. In the following section, we present several apparent permeability correction correlations obtained either from theoretical analyses or experiments.
15
ACCEPTED MANUSCRIPT
3.3 Apparent Permeability Correction Due to velocity slippage and gas molecular diffusion, apparent gas permeability of shale and tight gas matrices is higher than the absolute permeability. To facilitate the inclusion of this phenomenon into reservoir simulators, many multipliers have been proposed to correlate the apparent permeability to the absolute permeability in different flow regimes. Representative multipliers are introduced below.
(10)
SC
b ka = 1 + k∞ p
RI PT
3.3.1 Klinkenberg’s correlation for slip flow regime Experimental studies on permeability enhancement effects in tight formations date back to the early 20th century. In 1941, based on gas flooding experiments, Klinkenberg (1941) correlated the apparent gas permeability via ka to the absolute permeability
M AN U
where b is the Klinkenberg factor and p is the average pressure across the core. The Klinkenberg factor is usually obtained by matching experimental data. Compared with the slip boundary conditions discussed in the previous section, Klinkenberg’s correction is equivalently first-order since Kn ∝ 1/ p . Klinkenberg’s correction can be applied in the low Knudsen number range (<0.1), therefore it is widely adopted for simulating low permeability gas reservoirs. The Klinkenberg factor can be expressed as a function of the absolute permeability and the rock porosity. For example, Jones and Owens (1979) formulated it as
b = 12.639 ( k∞ )
−0.33
, Sampath and Keighin (1981) further integrated porosity into it b = 13.851( k∞ / φ )
Recently, square root correlation was proposed by Civan (2010) b = β ( k∞ / φ )
−1/2
, where
−0.53
.
depends on
TE D
molecule type. All the above Klinkenberg factor correlations are obtained by matching experimental data. For unconventional gas reservoirs, since the rock property varies significantly, usage of the above formulations should be cautious.
EP
3.3.2 Apparent permeability for transitional flow regime So far, no satisfactory apparent permeability correction has been developed for the transitional flow regime due to its complexity. Beskok and Karniadakis (1999) brought up a general relationship between Knudsen number and the apparent permeability multiplier
4 Kn f ( Kn ) = (1 + α1 Kn ) 1 + 1 − α 2 Kn
(11)
AC C
where α1 and α 2 are empirical parameters, depending on the cross-section shape of the flow channel. The development of accurate apparent gas permeability in transitional flow regime is still under pursuit. 3.3.3 Apparent permeability for free molecular flow regime Within the free molecular or Knudsen flow regime, the apparent gas permeability can be calculated by considering the diffusivity for Knudsen diffusion from gas kinetics
1 1 8 RT D = du = d 3 3 πMA where u is the thermal velocity of gas molecules, R is the gas constant, and MA is the gas molecular weight.
16
(12)
ACCEPTED MANUSCRIPT
The derivation of Knudsen diffusion coefficient can be found in Knudsen (1909) and Pollard and Present (1948). The apparent (total) gas flux Jt can be calculated by combining the convective flux J c resulted from the absolute permeability and the flux J k resulted from the Knudsen diffusion
k (13) ∇ p = J t = J c + J k = − ∞ ∇ p + D∇ C µg µg By rearranging the above equation, we can get an apparent gas permeability formulation similar to the Klinkbenberg’s correlation (Ertekin et al. 1986, Sun et al. 2015) ka
b ka = 1 + a p is gas viscosity, and
(14)
is gas compressibility factor.
SC
Where ba = pc g µ g D / k ∞ ,
k∞
RI PT
−
3.4 Non-linear and Non-Darcy Flow Model
M AN U
Javadpour et al. (2007) used a series of Gaussian distribution functions to match the experimental data and reproduced the Knudsen diffusion coefficient. However, it should be noted that the Knudsen diffusion coefficient is actually a function of Kn and the Knudsen diffusivity above is only the extreme case when Kn is infinitely large, as pointed out by Pollard and Present (1948).
In the past few years, a series of experimental studies of high velocity non-Darcy flow in propped hydraulic fractures were carried out at Colorado School of Mines (CSM) (e.g. Lopez 2007, Lai 2010, Lai et al. 2012). These studies concluded that the Barree and Conwey model (2004) yields a much better match than Forchheimer’s equation to experimental data of non-Darcy flow at extremely high flow rate, as shown in Figure 16.
TE D
Here Reynolds number is defined as
N Re =
ρv µτ
(15)
EP
where ρ is density (g/cm3), v is superficial velocity (cm/s), µ is viscosity (cp), and τ is the inverse of the characteristic length (100/cm).
AC C
Permeability ratio is defined as
η=
kapp kd
=
kmin 1 − kmin kd + kd (1 + N Re )
(16)
where kapp is apparent rate-dependent permeability, kd is constant Darcy permeability at low flow rate, and kmin is the minimum permeability at high flow rate.
17
ACCEPTED MANUSCRIPT
707.2 m3/D
SC
η
RI PT
283,168.2 m3/D
M AN U
NRe Figure 16. Experimental data of non-Darcy flow through propped fractures (markers) agrees very well with the Barree and Conway model (dashed line) (Lai et al. 2012).
Figure 16 shows that non-Darcy flow happens under high flow rates and the Barree and Conway (2004 and 2007) model is capable of matching the experimental data in the entire range from low to high flow rates. However, it may be insufficient to use Forchheimer’s equation to characterize the observed flow behavior (Lai et al. 2012).
AC C
EP
TE D
Reservoir simulations conducted at CSM, however, found that the Barree and Conway model gives very similar results to Forchheimer’s equation when dealing with single-phase flow using equivalent parameters (AlOtaibi and Wu, 2010 and 2011), as shown in Figure 17. Even for multiphase flow, the two models achieve very similar results if equivalent parameters are used (Wu et al., 2010; Wu et al. 2011a and 2011b; Zhang et al., 2014). Nevertheless, the laboratory studies (Lai et al. 2012) indicated that if extremely high flow rates through hydraulic fractures occur, the Barree and Conway model is better and preferred, while under normally high velocity flow condition, Forchheimer’s equation should be sufficient for non-Darcy flow modeling in fractured shale gas reservoirs.
18
ACCEPTED MANUSCRIPT
Figure 17. Pressure-transient responses of non-Darcy flow in a naturally fractured reservoir using Forchheimer’s equation and the Barree and Conway model (Al-Otaibi and Wu, 2011).
RI PT
Note that most studies on non-linear flow behavior in both conventional and unconventional reservoirs focused on high velocity flow. There are very few studies or little understanding on non-linear flow at low flow rate or low pressure gradient, which has been observed in experiments (e.g. Zeng et al. 2010, Lei et al. 2008). Based on laboratory observations, it is speculated that there may exist a certain threshold value of pressure gradient, below which flow in extra-low permeability rock may not occur. This threshold pressure phenomenon may also exist in shale and tight gas reservoirs or heavy oil reservoirs.
3.5 Reservoir simulation coupling geomechanics
TE D
M AN U
SC
In the past two decades, Thermal-Hydrological-Mechanical-Chemical (THMC) modeling have been widely used in simulating various subsurface engineering problems such as CO2 sequestration, geological disposal of nuclear waste, conventional and unconventional oil/gas production, acid gas injection, etc. Modeling THMC processes are very challenging in consideration of three aspects: 1) it is difficult to develop a THMC code that is fully functional for simulating all the processes, 2) finding constitutive relationships for couplings that are theoretically sound requires interdisciplinary expertise and there is no universal solution, 3) the realistic validation of THMC models is restrained by very rare data sets that can cover all THMC processes. For gas reservoirs without water influx from a contiguous water-bearing zone, rock compaction effect is generally non-negligible. Reservoir pressure could decrease significantly due to production, thus leading to rock deformation, which could be particularly remarkable for over-pressured shale and tight gas reservoirs. For example, Barnett Shale has a higher pore pressure gradient of about 0.54 psi/ft (Bowker 2007, Cipolla et al. 2010) than normal pressure gradient about 0.45 psi/ft. Thus pore pressure decrease induced by production could significantly affect in-situ stress conditions, resulting in rock compaction.
AC C
EP
THM codes have been available since the early 1980s, mostly limited to single phase liquid flow (e.g. Noorishad et al. 1984). A growing number of THM codes for simulating nuclear waste disposal near-field processes were developed in the 1990s. One common approach for simulating THM processes is to link an existing TH code to a mechanical code such as TOUGH-FLAC (Rutqvist 2011) and IPARS-JAS3D (Minkoff et al. 2003). For single phase flow, most of the THM codes solve the coupled equations simultaneously in a so-called fully coupled (monolithic) method, including ROCMAS (Noorishad and Tsang 1996, Rutqvist et al. 2001), COMPASS (Thomas and Sansom 1995), FRACON (Nguyen 1996), THAMES (Ohnishi and Kobayashi 1996), CODE_Bright (Olivella et al. 1994), and Open-Geosys (Wang and Kolditz 2007). For some applications, a multiphase flow approach is required. Exemplary codes that solve the multiphase flow and mechanics using the fully coupled (monolithic) methods include DYNAFLOW (Preisig and Prévost 2011), TOUGH-CSM (Winterfeld and Wu 2012), CODE_Bright (Olivella et al. 1994), and COMPASS (Thomas and Samson 1995). Recently, Kim (2010) advocated a method of sequential coupling of geomechanics and multiphase flow. There are many coupling methods developed for the THM models (Settari et al. 2001, Longuemare et al. 2002, Minkoff et al. 2003, Tran et al. 2004, Samier et al. 2008, Kim et al. 2009, Zhang 2013, Wang et al., 2014, 2015 and 2016), which can be categorized into three kinds: loose coupling, iterative coupling and full coupling. 1. Loose coupling method: For a loose coupling method, the reservoir simulator performs fluid flow calculations at each time step and the flow solutions are passed to the geomechanical model at a selected time step for stress calculations. This approach is also called one-way coupling because only flow solutions are passed for geomechanical calculations while geomechanical solutions do not feedback to 19
ACCEPTED MANUSCRIPT
flow calculations. Basically, two sets of equations for heat and fluid flows as well as geomechanics are solved in the algorithm. The loose coupling method solves the geomechanics after a certain number of time steps of fluid flow. Several THM simulators employ the loosely coupled procedure to solve fluid flow and geomechanics, e.g. ATH2VIS (Longuemare et al. 2002) and IPARS-JAS3D (Minkoff et al. 2003).
M AN U
SC
RI PT
2. Iterative coupling method: The iteratively coupled procedure solves the primary variables for heat and fluid flows and geomechanics sub-systems separately and sequentially in each time step. Usually the fluid flow equation systems are solved first and the intermediate solutions are then passed to the geomechanics system (Xiong et al., 2013; Zhang et al., 2016a and 2016b). The solutions of geomechanics equations then feed back to fluid flow system until the total equation systems converge, e.g. GEOSIM (Settari et al. 2001) and Rocflow (Wang and Kolditz 2007). This approach is a two-way coupling in each time step (Longuemare et al. 2002). One typical example of this iteratively coupled method is explicitly coupled method, which allows one iteration per time step, e.g. FRACTure (Kohl et al. 1997). Another implementation of the sequentially coupled approach is linking an existing TH code to a mechanical code, such as TOUGH_FLAC (Rutqvist and Tsang 2002, Rutqvist et al. 2002). Kim (2010) also studied the sequential coupling method for multiphase flow considering geomechanics.
TE D
3. Full coupling method: As the tightest coupling method, non-linear partial differential equation systems describing heat and fluid flows as well as geomechanics are solved simultaneously in each time step (Minkoff et al. 2003; Zhang et al., 2012a, 2012b and 2012c; Wu et al. 2014; Zhang et al., 2016c). Most of the fully coupled codes are limited to single phase liquid flow (e.g. Noorishad et al. 1984). A growing number of codes have been developed to simulate nuclear waste disposal near-field processes since the 1990s, such as FRACON (Nguyen 1996), ROCMAS (Noorishad and Tsang 1996, Rutqvist et al. 2001), and Open-Geosys (Wang and Kolditz 2007). There are also some codes that fully couple multiphase flow and mechanics, including CODE_Bright (Olivella et al. 1994), COMPASS (Thomas and Sansom 1995), DYNAFLOW (Preisig and Prévost 2011a), and TOUGH-CSM (Winterfeld and Wu 2012).
AC C
EP
Besides the common THM coupling studies abovementioned, some THMC procedures were recently developed to solve the geomechanics coupled flow problems. Yin et al. (2011) developed a fully coupled THMC model for CO2 injection by finite element methods to assess the potential formation damage around the wellbore. To avoid the calculation complexity, some THMC models were developed by linking two existing codes. FLAC3D and TOUGH-REACT has been linked together to investigate the impacts of hydraulic, thermal, and chemical processes on the evolution of enhanced geothermal system and nuclear waste disposal system (Taron et al. 2009, Rutqvist et al. 2014). Kim et al. (2015) presented a sequential implicit algorithm of chemo-thermo-poromechanics for fractured geothermal reservoirs, by linking ROCKMECH to TOUGH-REACT. Thereby, sequential iteration approach is proved sufficient to couple THMC processes together. To evaluate the coupling approaches, Tran et al. (2009) proposed three indexes: accuracy, adaptability, and running speed. Accuracy refers to how close the numerical results are to the real or benchmark solutions. Adaptability, in another word, flexibility, means how easy the existing or mature flow simulators and geomechanics simulators can be coupled without large code change or subsequent maintenance. Running speed is related to computational efficiency and is an important factor for practical full-field simulations. Explicit coupling approach has very good adaptability and high running speed because of loose coupling between two independent simulators but has poor accuracy due to one-way data transfer, usually from flow to 20
ACCEPTED MANUSCRIPT
geomechanics only at selected time steps. Iterative coupling also has quite good adaptability but less running speed than explicit coupling because the geomechanics computations are performed at each time step, on the other hand, its two-way data exchange between flow and geomechanics yields better accuracy than explicit coupling. Full coupling approach has the best accuracy and is unconditionally stable, but it requires much coding work for the tight coupling and does not run as fast as the other two methods.
Change of rock properties
SC
Geomechanical Effect
Change of pore pressure
M AN U
Gas Flow
RI PT
The primary interest to incorporate rock compaction into reservoir simulation lies in its influence on fluid flow, primarily through the change of rock properties, i.e. porosity and permeability (Zhang et al. 2015). Thus it is a two-way interaction between geomechanical effect and gas flow, where gas flow changes pore pressure and insitu stress, and stress-induced changes of porosity and permeability in turn affect fluid flow, as shown in Figure 18.
Figure 18. Two-way fluid flow and geomechanics interaction. Current industry practice in coupling geomechanical effect with reservoir simulation is to correlate the change of rock properties directly with the change of pore pressure instead of stress. For example, the conventional uncoupled reservoir simulator usually approximates the change of porosity as a function of pore pressure through pore volume compressibility by c P−P φ = φ0 e p ( 0 ) ≈ φ0 1 + c p ( P − P0 )
(17)
k
γ ( P ) = TM ( P) *
AC C
EP
TE D
where P0 is the reference pore pressure at which the porosity is 0. This equation approximates the change of porosity as a function of pore pressure with constant pore volume compressibility, which is one simplified method to capture rock deformation effect (Aziz and Settari 1979, Ertekin et al. 2001). In addition to the porosity approximation, conventional uncoupled reservoir simulator also includes rock compaction effect on absolute permeability through correlations with pore pressure. For example, commercial simulator Eclipse (2016) correlates absolute permeability and pore pressure through a transmissibility multiplier ij +
1 Aij 2
dij
(18)
where γ is the pressure dependent transmissibility, TM is a pressure dependent transmissibility multiplier, k
ij +
1 2
is the absolute permeability, Aij is the fluid flow area, and dij is the length for pressure gradient determination. It has been realized that these porosity and permeability approximations are not sufficient to capture the flow behavior in stress-sensitive reservoirs, therefore a variety of methods for improving coupling between fluid flow and geomechanics have been proposed (Dean et al. 2006, Gutierrez et al. 2001, Minkoff et al. 2003, Settari and Walters 2001, Tran et al. 2009, Xiong et al. 2013, Zhang et al. 2015). In geomechanics, rock properties such as porosity and permeability are functions of effective stress rather than pore pressure, as are calculated along with
21
ACCEPTED MANUSCRIPT
flow variables during simulation. In other words, the way to calculate stress may vary for different coupling methods, but in general the following correlations are incorporated into coupled geomechanical simulation k = k (σ ) , φ = φ (σ )
where
(19)
is the stress variable.
RI PT
These correlations are able to produce more accurate simulation results for stress-sensitive gas reservoirs.
M AN U
SC
Based on the illustration of major gas flow and transport mechanisms in shale and tight gas reservoirs, corresponding mathematical correlations are discussed and compared. Though progress has been made for describing an individual mechanism, there still lack of integrated reservoir simulators that can couple most or all of these mechanisms in a comprehensive way. On the other side, involvement of multi-physics problems in pervasive reservoirs significantly increases the computation costs required by serial gas reservoir simulators. To improve the computation efficiency of reservoir simulators, many high performance computing schemes that can simultaneously call thousands of processors have been developed for large-scale reservoir simulations (e.g. Dongarra et al. 1989, Wheeler and Smith 1990, Killough and Bhogeswara 1991, Zhang et al. 2001, Fung and Dogru 2008, Wang et al. 2015). Details of high performance computing in reservoir simulations are not discussed here, due to limited space.
AC C
EP
TE D
4. Pore-Scale Modeling Pore-scale modeling is becoming a commercial routine for evaluating oil and gas reservoirs at micro- and nanoscale, both the imaging techniques and simulation tools have been well developed during the past several decades (Blunt et al. 2013). As a prerequisite for fluid flow simulation, pore imaging techniques generate threedimensional digital porous geometries that are representative of the pore structures in the reservoir rocks. Considering the similarity of digital porous geometries to the real pore structures in rock samples, they can be divided into two major types. The first is pore network models, which statistically represent or extract the typical pore patterns or properties by numerically generating virtual pore networks. The pore network models are generally constructed based on either real 2D rock images (Wu, 2006) or ideal geometrical shapes, such as spheres, cylindrical tubes, and irregular grains (Wang, 2014). The second is real pore networks, which are real pore structures obtained from reconstruction of micro- and nano-Computed Tomography (CT) scanning or Focused Ion Beam-Scanning Electron Microscope (FIB-SEM) images. Pore network models are not limited in size as algorithms are scalable, while real pore networks are constrained by the size and number of the images acquired. More or less, the size of real pore networks is a 1000 times as large as its resolution in magnitude. Pore networks are places inside which fluid flow simulations are carried out, therefore, corresponding to pore network models and real pore networks, simulation tools are categorized into two families as well: network modeling and direct modeling, respectively (Blunt et al. 2013). For single phase flow, the former, due to the geometrical simplicity of the pore network models, mostly applied specific forms of Darcy’s law to account for viscous flow inside regular or irregular capillary tubes interconnecting the big pore spaces. The latter, in contrary, directly solves or indirectly recovers Navier-Stokes equation by discretizing the real pore networks by voxels. It is thus evident that computation cost of direct modeling is much higher than that of the network modeling.
22
ACCEPTED MANUSCRIPT
4.1 Challenges of unconventional reservoirs Unconventional reservoirs, with irregular nano-pores residing in both organic matter and inorganic minerals (Kuila and Prasad 2013), pose great challenges for pore-scale imaging and modeling techniques developed for conventional reservoir characterization.
M AN U
SC
RI PT
On one hand, majority of the shale pores are of nanometers and porosity is quite low, requiring high resolution imaging techniques to provide details of pore spaces for digital core reconstruction. Micro-CT with the resolution of hundreds of nanometers is insufficient to describe the nano-pores, which in general have a median pore size of tens of nanometers (Honarpour et al. 2012, Du and Chu 2012). Nano-CT can achieve a higher resolution of 50 nm (Guo et al. 2015), it is close to most of the reported median pore throat diameters in shale reservoir matrices. Thus, nano-CT scanning would probably snap off most of the interconnecting throats that are narrower than the half of the resolution, resulting in smaller porosity and permeability values. During the last ten years, FIB-SEM imaging has been widely applied to examine the mineralogy, organic matter, matrix texture, porosity, pore size distribution, and pore shapes and types of shale and tight reservoir rocks (e.g. Chalmers et al. 2009, Sisk et al. 2009, Curtis et al. 2010, Lemmens et al. 2010). Extraordinarily, FIB-SEM is capable of going down to about 1 nm to reveal more specific features, unmatchable by nano-CT. However, unlike CT, FIB-SEM is destructive because of layer by layer milling, which also ablates the rock matrix between two slices of 2D images, causing 1-layer thickness information missing. This missed information or discontinuity in-between may either create or lose some connectivity between pores adjacent to each other on neighboring images, yet this probability of occurrence would be quite low due to relatively sparse pore distribution in shale/tight rock matrices. In spite of this deficiency, FIB-SEM imaging is the most precise and reliable technique that is currently available for uncovering details of shale and tight rock matrices.
EP
TE D
On the other hand, for the pore network models, it is becoming increasingly complex to reconstruct representative pore network model for shale, in consideration of the strong heterogeneity and low connectivity of irregular pores. The pores in shale are more irregular than those in the conventional rocks, using simplified geometric shapes as analogs is far from representative of real pore structures. Besides, including micro-fractures, pore size covers a wide range with strong heterogeneity, but for a given pore, there is no clearly distinguishable pore bodies and pore throats, which is challenging to be statistically described. Also, the connectivity among pores is very low, and there are many isolated pores existing in the shale matrix. Moreover, to further classify pore walls as organic or inorganic, which affects the fluid solid interaction, adds another variable into the pore network models. These unconventional features of shale are critical, since they directly affect the pore network models to be built, and further determine the governing equations to be used for fluid flow simulation.
AC C
Based on our flow regime analyses in previous sections, slip boundary condition should be incorporated into Navier-Stokes equation or lattice-Boltzmann model in order to accurately simulate the gas slip flow regime under reservoir conditions. Additionally, in pores of a few nanometers, more physical behavior of transitional flow and even free molecular flow should be considered. Moreover, it is clear that adsorption and desorption play a critical role in organic-rich shale formations (Lu et al. 1995). Boţan et al. (2015) employed grand-canonical Monte Carlo to account for Langmuir adsorption of methane and MD to capture the gas transport without defining flow regimes in nano- and subnano-porous domains, then they incorporated the adsorption-transport results into a lattice model for larger scale simulation of gas transport in 3D real pore networks. To predict the permeability of shale and tight porous media, Singh and Javadpour (2016) proposed a Langmuir slip permeability model for slip flow and transitional flow regimes which uses slip coefficients derived from Langmuir adsorption data and considers pore size enlargement due to desorption.
23
ACCEPTED MANUSCRIPT
Herein, our simplified justification of pore-scale modeling is limited to classification of in-situ flow regimes, without considering adsorption and desorption, and fluid states under confined conditions. Regarding the flow regimes existing in unconventional gas reservoirs, recent progresses in network modeling and direct modeling of single phase gas flow in micro- and nano-porous media are reviewed and discussed below.
RI PT
4.2 Pore network models and real pore networks Network models as an analog for porous media was first reported 60 years ago (Fatt 1956), and have gained significant momentum during the last ten years. For the purpose of simulating gas flow in unconventional reservoirs, efforts have been focused on advancing the representativeness and verisimilitude of pore network models to the real pore structures as well as incorporating more flow mechanisms that are beyond viscous flow into the gas flow models.
AC C
EP
TE D
M AN U
SC
To date, several different classes of methods have been well developed for generating pore network models, essentially these methods try to preserve the pore size distribution and pore connectivity by extracting pore structure information from real images (Okabe and Blunt 2004, Wu et al. 2004) or to generate porous media by simulating the sedimentation and diagenesis processes (Bakke and Øren 1997). Using Markov Chain Monte Carlo simulation (MCMC), Wu et al. (2006) developed a stochastic pore space reconstruction method, which is able to create realistic pore network models based on the morphological information extracted from 2D rock images. With one 2D image, by assuming that the image characteristics are statistically the same in all x, y, and z directions, they reconstructed 3D pore networks for tight mudrock, which has a porosity of 11.73% and a permeability of 0.013mD. In addition, given three 2D images obtained from perpendicular planes, the model is capable of generating a heterogeneous rock cube with qualitatively reproduced morphology and porosity (Wu et al. 2006). Ma et al. (2014) used this MCMC method to extract node-bond pore networks from 2D images of mudstone, the one they obtained has a porosity of 2.9% with more than 80% pores less than 10 nm in radius and the others less than 50 nm. Based on SEM images with resolution of 5 nm taken for shale samples, Chen et al. (2015) constructed four shale pore networks using this MCMC method. The reconstructed pore networks have a permeability range of 0.176-26.8 and pore size between 5 and 80 nm with frequency peak around 24 nm, fairly representing the nano-pores in shale matrix. Dong and Blunt (2009) modified maximal ball algorithm (Silin and Patzek 2006) to extract simplified pore networks from micro-CT images. These pore network models are comprised of pore bodies interconnected by pore throats, commonly represented by spheres and cylinders, respectively. Also, a pore body or a uniform throat can have an angular cross-section (Blunt et al. 2013). The coordination number, and pore body and throat radii can be varied as a distribution during pore network generation. By significantly shrinking these parameters, pore network models generated are being used as analogs for shale matrices (Mehmani et al. 2013, Zhang et al. 2015). Øren and Bakke (2002) developed process-based procedure that directly models the grain sedimentation, compaction, and diagenesis processes to generate pore network models for sandstone. It is expected to be applicable for shale pore network generation, but the computation will be intensive as well. Note that studies on pore network models discussed here are some focusing on unconventional reservoirs among many others. In view of the nano-pores in shale and tight matrices, imaging techniques used for conventional reservoirs, like optical microscope and micro-CT, are incapable of capturing the pore structures due to low resolution. To meet the demand of unconventional reservoir characterization, nano-CT and FIB-SEM have been widely applied to probe the nano-pores in shale and tight matrices. Guo et al. (2015) adopted the non-destructive nano-CT to scan the internal structures of four shale samples with a resolution of 65 nm in pixel or 32.5 nm in pore radius. Using a commercial software, 3D pore network models with an average porosity of 4.6% and a pore radius frequency 24
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
peak at 103 nm have been visualized by stitching up the stacks of 2D images. For the shale samples they tested, nano-CT captured the pore size distribution and pore connectivity, and distinguished the pores in organic and inorganic matter, however, for other shale or tight matrices with much smaller pore radius frequency peak, it might be difficult to achieve representative results. FIB/SEM has demonstrated much higher resolution in visualizing nano-pores in unconventional reservoir rocks (e.g. Chalmers et al. 2009, Sisk et al. 2009, Curtis et al. 2010, Lemmens et al. 2010 and 2011, Zhang et al. 2012). By following procedures similar to those for nano-CT, 3D pore network models can be reconstructed. Figure 19 is a single 2D FIB-SEM image of a stack taken for a core sample from a tight sandstone gas reservoir. The dimension of this image is Y=23.65 µm and Z=16.9 µm with a resolution of 50 nm (25 nm in pore radius). Low resolution was used here in considering that the field of view of SEM is negatively correlated to its resolution. It is obvious that the pores or the micro-fractures are highly irregular, several of which are a few microns in length. By stacking up 200 2D FIB-SEM images in the depth, we reconstructed the 3D digital core with the thickness of 10 µm along X-axis, as shown in Figure 20, where the 2D image in Figure 19 faces front. Through further discretization of the grey-scale images into binary ones, we can extract the pore structures out of these images as a 3D real pore network (Figure 21). This pore network indicates dominant flow pathways along X-axis, while in Y and Z-axes, the pore connectivity is very poor. The porosity of this 3D digital core is 6.5% and the intrinsic permeability measured by lattice-Boltzmann simulation is 1.24 µD. Real pore networks built upon FIB-SEM images authentically describe the pore structures of unconventional reservoir rocks, nonetheless, its application is still limited by the size of the view field and cumbersome laboratory work.
Figure 19. SEM image of a core sample from a tight sandstone gas reservoir.
25
SC
RI PT
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
Figure 20. 3D digital core reconstructed from FIB-SEM images of a tight sandstone core.
Figure 21. 3D real pore networks extracted from FIB-SEM images of a tight sandstone core.
4.3 Network modeling
Once the pore network models for shale or tight rocks are constructed, mathematical models for simulation of gas flow through the pore networks can be built and implemented. Since Darcy’s law derived from NavierStokes equation is no longer accurate for gas slip flow, many modifications to the intrinsic permeability have been proposed to calculate the apparent permeability . One of the earliest correlations that depict gas slip flow through reservoir rocks is the Klinkenberg’s equation (Klinkenberg 1941)
26
ACCEPTED MANUSCRIPT
λ ka = k 1 + 4c r
(20)
where is a factor slightly less than 1.0. Although it was developed for the purpose of matching core-scale data, the basic assumption is the slippage of gas molecules on solid surface defined by Knudsen number.
RI PT
To simulate gas flow in nano-pores in shale matrix, Javadpour et al. (Javadpour et al. 2007, Javadpour 2009) proposed to calculate the total mass flux J t by simultaneously combining the slip flow J s and Knudsen diffusion
Jk . Jt = J s + J k
(21)
TE D
M AN U
SC
The slip correction used comes from Brown et al. (1946) instead of Klinkenberg’s. Their expression for a cylindrical nano-pore matches experimental data measured from pores of ~200 nm in diameter with an average error of 4.5% (Javadpour 2009). Note that here the simplicity of modifying Darcy’s law is preferred over physics as the gas flow goes beyond = 0.1. Mehmani et al. (2013) extended this model (Javadpour 2009) to simulate gas flows in 3D node-bond pore networks with variable pore sizes. The gas slippage and Knudsen diffusion were only evaluated in cylindrical pore throats. Their results show that the pressure-dependent apparent permeability is affected by the fraction and spatial configuration of the nano-pores (Mehmani et al. 2013). By substituting a gas kinetics expression similar to Klinkenberg’s equation into the slip flux (Javadpour 2009), Ma et al. (2013) simulated gas flow through 3D pore network models statistically constructed upon 2D SEM images of shale samples using MCMC method (Wu 2006). Rather than solely using cylindrical pores, the network model contains throats close to triangle and square shapes, too. Also, in the flow modeling, non-ideal gas compressibility calculated by van der Waal’s equation of state was incorporated into the slip and Knudsen flux for comparison with ideal gas cases, revealing that non-ideal gas effect that favors methane flow can contribute tens of percent to the apparent permeability under reservoir conditions. Superposition of slip and Knudsen flux in total flux is questioned as double counting the contributions from slippage and Knudsen diffusion (Sakhaee-Pour and Bryant 2012, Ma et al. 2013), when compared to the dusty gas
where
AC C
viscous flow J v .
EP
model (DGM) developed by Mason and Malinauskas (1983), which only adds Knudsen diffusion J k to the
J t = J v + J k , J k = − DKn
is the Knudsen diffusivity, ∇ is the pressure gradient,
∇p RT
is the gas constant, and
(22) is temperature.
Sakhaee-Pour and Bryant (2012) chose the first-order slip model in the form of Klinkenberg’s equation over DGM to simulate gas flow in cylindrical nano-pores. Their model, which matches experimental data within 6% deviation, found that as gas production continues with pressure drop in reservoirs, the matrix permeability is expected to increase by 4.5 times due to release of adsorbed layer increasing pore radius and gas slippage. This helps explain the deaccelerating decline in production rate of shale gas wells after the early years (Sakhaee-Pour and Bryant 2012). Pore network modeling simulates the physical behavior of micro- and nanoscale gas flow by modifying the macroscopic correlations, providing a means of directly upscaling microscale physics to field scale applications. 27
ACCEPTED MANUSCRIPT
4.4 Direct modeling
RI PT
Direct modeling solves the gas flow by discretizing pore spaces obtained from SEM or nano-CT images into lattices. As shown in Table 1, either Navier-Stokes equation (Silin 2014) or LBM is capable of capturing the slip flow behavior by incorporating the slip boundary condition. In the transitional flow regime, seeing that Burnett equation is tedious to be solved, DSMC is more commonly used. However, for slow flow as encountered in reservoir matrix, DSMC experiences thermal fluctuations toward convergence, and even simulation of flows through nanoscale digital rock cubes is prohibitively expensive (Li and Sultan 2015). For free molecular flow, MD simulation is too computationally expensive for the representative micrometer scale porous media, only flows in nanotubes have been studied (Feng and Akkutlu 2015). Therefore, we will focus on lattice-Boltzmann method for simulating gas slip flow below.
M AN U
SC
LBM is a numerical method that models the evolution of the single-particle velocity distribution function n ( , t) over spatial lattices according to lattice-Boltzmann equation, which is derived by discretizing Boltzmann equation (He and Luo 1997a) and indirectly recovers NS equation (He and Luo 1997b). As a scalable method suitable for parallel computing (Xiao 2013, Wang 214), it has been widely used to simulate continuum fluid flows in complex geometries (Chen and Doolen 1998, Chen et al. 2014). The velocity distribution function describes the fraction of fluid particles with velocity $ at position and time t in each lattice. Velocity quadrature $ is designed to simplify the continuous molecular velocity distribution in the gas kinetics with an assigned weighting factor %&' , which sums to 1. Evolution of n ( , t) in LBM involves propagation and collision steps, as shown in the left and right equalities below
ni ( r + ci ∆t,t + ∆t ) = ni* ( r ,t ) + fi ( r ,t ) = ni ( r ,t ) + ∆i n ( r ,t ) where the post-collision particle velocity distribution ni* is moved from lattice
(23)
r to the neighbors r + ci ∆t in the
TE D
direction of ci , ∆i is the collision operator, and f i ( r ,t ) is the forcing term. In propagation, populations are required to move to the next lattice in each time step.
AC C
EP
When the fluid particles propagate to a rigid surface, boundary conditions must be applied to model the interactions on the surface. Standard LBM simulates the no-slip boundary condition (Yin et al. 2006) and thus only gives the intrinsic permeability of porous media. By properly treating the particle-surface interactions, many LBM schemes have been developed to simulate the slip boundary conditions. As discrete velocity quadratures approximate the continuum gas velocity distribution, some studies increased the number of velocity quadratures to achieve a high order discretization of the Boltzmann equation (Shan et al. 2006, Zhang et al. 2006). Ansumali et al. (2002) applied the completely diffusive boundary on the wall in LBM to simulate slip flow with ( < 0.03. Higher order velocity set for LBM have been combined with the relaxation time- ( correlation for slip simulation as well (Tang et al. 2008). Chen and Tian (2009) used the Langmuir slip model in the LBM to simulate Couette and Poiseuille flow in the slip flow regime, which matched the analytical solutions. In these abovementioned studies, only gas slippage on a flat surface was calculated, while simulation of slip flow along curved surface in 3D geometries is relatively rare. These considerations are important since the curved surface may result in negativity in slippage contribution to the flow velocity (Einzel et al. 1990, Lockerby et al. 2004, Szalmás 2007). Recently, Cho et al. (2013) conducted slip flow simulation through 3D non-overlapping fibers using multiple relaxation time LBM by relating relaxation time to ( and lattice size, and obtained permeability as a function of (, porosity, and fiber configuration. However, in their simulation, the effect of 28
ACCEPTED MANUSCRIPT
TMAC on slip flow is not investigated (Zhang et al. 2005). Besides, correlation of lattice size to relaxation time renders it complicated to vary the resolution of geometries for more accurate simulations.
SC
RI PT
At present, other than simple channels, cubes, and cylinders, very few LBM schemes with slip boundary condition have been tested for flows through complex porous media, which actually would pose many challenges. For example, pore size in porous media covers a range rather than a unique value, therefore ( is not uniform inside the porous media. However, many models, in essential, defined ( as a specific value but a function of pore size distribution (Zhang et al. 2005, Szalmás 2007, Kim et al. 2008, Tang et al. 2008). Hence in slip flow models, ( should be substituted by mean free path to serve as the input to account for the varying degree of rarefaction. Additionally, application of slip boundary condition in porous media is bounded by the largest ( = 0.1 at the smallest pore throats for the standard LBM, because the assumption of thermodynamic equilibrium cannot be met in Knudsen layer (Tang et al. 2008). To conquer these challenges, Wang (2014) derived a new bounce-back scheme by incorporating the Maxwell’s first-order slip boundary condition into the lattice-Boltzmann equation
2 −σ
σ
τ
(24)
M AN U
us = −
where us is the slip velocity, tangential stress vector τ = n ⋅ π ⋅ ( I − nn ) ,
n is the unit normal vector pointing
Τ from solid wall to fluid, I is the unit tensor. π = µ ∇u + ( ∇u ) is the viscous stress tensor of the local fluid,
is the viscosity, ∇u is the velocity gradient tensor, and
here represents the transpose operation.
Setting a correction term δ i ( r ,t ) in the bounce-back scheme (Ladd and Verberg 2001, Yin et al. 2006), the
δ i ( r ,t ) = −
where
)
ρ ac 1 i
cs2 bl
1 ( ci ⋅ n ) ci ⋅ ( I − nn ) ⋅ j ( r ,t ) ρ ( c i ⋅ n )( n ⋅ u b ) + 1 + 2bl
( ci ⋅ n ) ci ⋅ ub − ( ci ⋅ n )( ub ⋅ n ) −
EP
+
2a ci cs2
TE D
following expression for gas slippage correction is obtained (Wang, 2014)
a ci ci ⋅ ( I − nn ) ⋅ Π neq ,* − Π neq ⋅ ( I − nn ) ⋅ ci 2cs4
(
)
(25)
is the isothermal sound speed in a fluid at equilibrium, ρ is fluid density, *+ is the velocity of the rigid
AC C
surface, slip length bl = ( 2 − σ ) λ / σ ,
j
is the first-order momentum neq
Πneq,*
is the post-collision non-
equilibrium part of the second-order moment Π . In this correction, the boundary velocity is retained, and the mean free path of gas molecules and the TMAC from the Maxwell’s first-order slip boundary condition are both embodied. Also, note that when ,- approaches infinity, this correction term recovers the free-slip boundary condition (Yin et al. 2006). Nevertheless, since ,- is in the denominator, it cannot be assigned very small values to recover the no-slip boundary condition. This slip model has been used to simulate slip flow through three arrays of spheres, demonstrating that the apparent permeability increases with increasing slip length and increases with decreasing porosity, which are formulated into polynomials for upscaling purpose. Also, a local normal calculation method was proposed for this slip model to simulate flows through simple cubic arrays of overlapping spheres (Wang 2014). For complex, digitalized pore surfaces inside shale or tight sandstone matrices where the global normal vector is not available and for situations where the numerical stability of the global normal method is an issue, the method of local normal method offers a realistic approximation. 29
ACCEPTED MANUSCRIPT
RI PT
There are other methods that combined pore network flux models with lattice-Boltzmann simulation. For example, Allan and Mavko (2012) applied LBM to calculate the flux-averaged pore width for Knudsen number, which was further substituted into the unified flow model developed by Karniadakis et al. (2005) and Florence et al. (2007) to compute the apparent permeability. Using DGM, Chen et al. (2015) evaluated the contribution of Knudsen diffusion to the apparent permeability of 3D pore networks reconstructed from SEM images of shale by MCMC method. To calculate the Knudsen diffusivity, lattice-Boltzmann method with single relaxation time was adopted to simulate the diffusion process by taking pore size distribution into consideration. In these hybrid models, the contribution of Knudsen diffusion to the apparent permeability is only counted in the flux model.
SC
5. Conclusions In view of the characteristics of shale and tight gas reservoirs that are different from conventional natural gas reservoirs, a few mechanisms dominating gas flow and transport should be properly addressed in gas reservoir simulation, including gas desorption, rarefaction effect, non-Darcy flow, compaction effect, etc.
M AN U
Though many quantitative correlations for describing gas flow and transport mechanisms in shale and tight gas reservoirs have been proposed, there still lack of generalized and well-validated formulations that can recover the commingled gas behavior in shale and tight gas reservoirs with high heterogeneity. Current large-scale reservoir simulators for shale and tight gas reservoir simulation not only suffer from scarcity of accurate correlations for different flow and transport mechanisms, but also undergo high complexity of coupling them together. High performance reservoir simulations integrating all these predominant factors are still in urgent need.
TE D
Statistical pore network models are readily available for reconstructing variable pore networks from real SEM images of unconventional reservoir rocks by extracting the pore size distribution and pore connectivity information. FIB-SEM imaging and 3D digital core construction are the most authentic means to describe the pore characteristics inside unconventional reservoir rocks, the information completeness and field view are still limited though.
EP
Several pore network modeling methods have been developed for simulating gas flows in slip and transitional flow regimes by successfully capturing the molecular level physics in macroscopic expressions.
AC C
There are a few lattice-Boltzmann models developed for gas slip flow simulation, however, very few of them have been extended to simulate slip flow in 3D complex porous media due to geometrical challenges. To simulate gas flow in transitional flow regimes, DSMC that directly simulates the Boltzmann equation by probabilistically moving and colliding molecular swarms is more suitable.
Acknowledgement
The authors thank Dr. David Wood, the Chief Editor of Journal of Natural Gas Science and Engineering, for his invitation to write this comprehensive review paper.
30
ACCEPTED MANUSCRIPT
Reference Aidun, C. K. and Clausen, J. R. 2010. Lattice-Boltzmann method for complex flows. Annual Review of Fluid Mechanics, 42, 439-472. Allan, A.M. and Mavko, G. November, 2012. The effect of adsorption and diffusion on the gas permeability of kerogen. In 2012 SEG Annual Meeting. Society of Exploration Geophysicists.
RI PT
Al-Otaibi, A. and Y.-S. Wu. An alternative approach to modeling non-Darcy flow for pressure transient analysis in porous and fractured reservoirs, SPE-SAS-984, prepared for presentation at the 2011 SPE Saudi Arabia Section Technical Symposium and Exhibition held in AlKhobar, Saudi Arabia, 15–18 May 2011.
SC
Al-Otaibi, A. and Y.-S. Wu. Transient behavior and analysis of non-darcy flow in porous and fractured reservoirs according to the barree and conway model, SPE 133533, presented at the at the Western North America Regional Meeting, held in Anaheim, California, USA, 26–30 May 2010. Ansumali, S. and Karlin, I. V. 2002. Kinetic boundary condition in the lattice Boltzmann method. Physical Review E 66 (2): 026301.
M AN U
Aziz, Khalid, and Antonin Settari. 1979. Petroleum reservoir simulation. Chapman & Hall. Bakke, S. and Øren, P. E., 1997. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks. Spe Journal, 2(02), pp.136-149. Barree, R. D. and M. W. Conway, “Multuiphase Non-Darcy Flow in Proppant Packs,” SPE 109561, presented at the 2007 Annual Technical Conference and Exhibition, Anaheim, CA, 11-14 November, 2007.
TE D
Barree, R.D. and M.W. Conway, “Beyond Beta Factors: A Complete Model for Darcy, Forchheimer and TransForchheimer Flow in Porous Media,” SPE 89325, presented at the 2004 Annual Technical Conference and Exhibition, Houston, Texas 26-29 September, 2004. Beskok, A., Karniadakis, G. E. A model for flows in channels, pipes, and ducts at micro and nano scales. Microscale Thermophys Eng. 3(1), 43–77 (1999). doi:10.1080/108939599199864
EP
Bird, G. A. 2003. Molecular gas dynamics and the direct simulation of gas flows. New York. Oxford University Press Inc.
AC C
Blasingame, T.A., The characteristic flow behavior of low-permeability reservoir systems, SPE 114168, presented at the 2008 SPE Unconventional Reservoirs Conference held in Keystone, Colorado, U.S.A., 10–12 February 2008. Blunt, M. J., Bijeljic, B., Dong, H., Gharbi, O., Iglauer, S., Mostaghimi, P., Paluszny, A. and Pentland, C., 2013. Pore-scale imaging and modelling. Advances in Water Resources, 51, pp.197-216. Boţan, A., Ulm, F.-J., Pellenq, R. J.-M., and Coasne, B., 2015. Bottom-up model of adsorption and transport in multiscale porous media, Physical Review E, 91, 032133. Bowker, K. A. 2007. Barnett Shale gas production, Fort Worth Basin: issues and discussion. AAPG bulletin 91(4): 523-533. Brown, G. P., DiNardo, A., Cheng, G. K. and Sherwood, T. K. 1946. The flow of gases in pipes at low pressures. Journal of Applied Physics, 17(10), pp.802-813.
31
ACCEPTED MANUSCRIPT
Bustin, R. M., A. M. M. Bustin, X. Cui, D. J. K. Ross, and V. S. Murthy Pathi. 2008. Impact of shale properties on pore structure and storage characteristics, presented at the SPE Gas Production Conference, 16-18 November, Fort Worth, Texas, SPE 119892. Bybee, K. 2008. Horizontal wells in tight gas sands: risk management to maximize success. JPT, 60(10): 61-63.
RI PT
Chalmers, G. R., Bustin, R. M., and Power, I. 2009. A Pore by any other name would be as small: the importance of meso- and microporosity in shale gas capacity. Presented at the AAPG Annual Convention and Exhibition, Denver, Colorado. 7-10 June. Chen, L., Kang, Q., Mu, Y., He, Y.L. and Tao, W.Q., 2014. A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications. International Journal of Heat and Mass Transfer, 76, pp.210-236.
SC
Chen, L., Zhang, L., Kang, Q., Viswanathan, H. S., Yao, J. and Tao, W., 2015. Nanoscale simulation of shale transport properties using the lattice Boltzmann method: permeability and diffusivity. Scientific reports, 5.
M AN U
Chen, S. and Doolen, G. D. 1998. Lattice Boltzmann method for fluid flows. Annual Review of Fluid Mechanics 30 (1): 329-364. Chen, S. and Tian, Z. 2009. Simulation of microchannel flow using the lattice Boltzmann method. Physica A: Statistical Mechanics and its Applications 388 (23): 4803-4810. Chien, M. C. H., Wasserman, M. L., Yardumian, H. E., Chung, E. Y., Nguyen, T., and Larson, J. January, 1987. The use of vectorization and parallel processing for reservoir simulation. In SPE Symposium on Reservoir Simulation. Society of Petroleum Engineers.
TE D
Cho, H., Jeong, N., and Sung, H. 2013. Permeability of microscale fibrous porous media using the lattice Boltzmann method. International Journal of Heat and Fluid Flow 44: 435-443. Cipolla, C. L. Modeling Production and Evaluating Fracture Performance in Unconventional Gas Reservoirs, JPT, Distinguished Author Series, pp.84-90, September 2009.
EP
Cipolla, C. L., E. P. Lolon, J. C. Erdle, and B. Rubin, Reservoir Modeling in Shale-Gas Reservoirs, Reservoir Evaluation and Engineering, August 2010.
AC C
Cipolla, C. L., N. R. Warpinski, M. J. Mayerhofer, E. P. Lolon, and M. C. Vincent, The Relationship Between Fracture Complexity, Reservoir Properties, and Fracture Treatment Design, SPE 115769, presented at the SPE Annual Technical Conference and Exhibition, Denver, 21–24 September 2008. Civan, F., Rai, C. S., and Sondergeld, C. H. 2011. Shale-gas permeability and diffusivity inferred by improved formulation of relevant retention and transport mechanisms. Transport in Porous Media, 86(3), 925-944. doi: 10.1007/s11242-010-9665-x Civan, F. 2010. Effective correlation of apparent gas permeability in tight porous media. Transport in Porous Media, 82(2), 375-384. doi: 10.1007/s11242-009-9432-z Clarkson, C. R., Solano, N., Bustin, R. M., Bustin, A. M. M., Chalmers, G. R. L., He, L., Melnichenko, Y. B., Radliński, A. P., and Blach, T. P. 2013. Pore structure characterization of North American shale gas reservoirs using USANS/SANS, gas adsorption, and mercury intrusion. Fuel, 103: 606-616.
32
ACCEPTED MANUSCRIPT
Curtis, M. E., Ambrose, R. J., and Sondergeld, C. H. 2010. Structural characterization of gas shales on the microand nano-scales. Presented at the Canadian Unconventional Resources and International Petroleum Conference, Calgary, Alberta. SPE-137693-MS. http://dx.doi.org/10.2118/137693-MS. Darabi, H., Ettehad, A., Javadpour, F., and Sepehrnoori, K. 2012. Gas flow in ultra-tight shale strata. Journal of Fluid Mechanics, 710, 641-658. doi: 10.1017/jfm.2012.424
RI PT
Dean, R. H., X. Gai, C. M. Stone, and S. E. Minkoff. 2006. A comparison of techniques for coupling porous flow and geomechanics. SPE Journal, 11(01): 132-140. DeBaun, D., Byer, T., Childs, P., Chen, J., Saaf, F., Wells, M. and Crumpton, P. January, 2005. An extensible architecture for next generation scalable parallel reservoir simulation. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
SC
Denney, D. 2008. Rock type-understanding productivity in tight gas sands. JPT, 10: 53-56.
M AN U
Dogru, A. H., Fung, L. S., Middya, U., Al-Shaalan, T. M., Pita, J. A., HemanthKumar, K., and Hahn, W. A. October, 2009. A next-generation parallel reservoir simulator for giant reservoirs. In SPE/EAGE Reservoir Characterization and Simulation Conference. Dong, H. and Blunt, M. J. 2009. Pore-network extraction from micro-computerized-tomography images. Physical Review E, 80(3), p.036307. Dongarra, J., Duff, I., Gaffney, P. A. T. R. I. C. K., and McKee, S. 1989. Vector and Parallel Computing. Wiley, New York. Du, L. and Chu, L. 2012, January. Understanding anomalous phase behavior in unconventional oil reservoirs. In SPE Canadian Unconventional Resources Conference. Society of Petroleum Engineers.
TE D
Eclipse, http://www.software.slb.com/products/eclipse.
EIA (U.S. Energy Information Administration). July 2011a. Review of Emerging Resources: U.S. Shale Gas and Shale Oil Plays.
EP
EIA (U.S. Energy Information Administration). April 2011b. World Shale Gas Resources: An Initial Assessment of 14 Regions Outside the United States. Einzel, D., Panzer, P., and Liu, M. 1990. Boundary condition for fluid flow: curved or rough surfaces. Physical Review Letters 64 (19): 2269-2272.
AC C
Ertekin, T., King, G. R., and Schwerer, F. C. 1986. Dynamic gas slippage: a unique dual-mechanism approach to the flow of gas in tight formations. SPE (Society of Petroleum Engineers) Format. Eval., (United States), 1(1). Ertekin, Turgay, Jamal Hussein Abou-Kassem, and Gregory R. King. Basic applied reservoir simulation. Richardson, TX: Society of Petroleum Engineers, 2001. Falk, K., Coasne, B., Pellenq, R., Ulm, F. J., and Bocquet, L. 2015. Subcontinuum mass transport of condensed hydrocarbons in nanoporous media. Nature Communications, 6. Fatt, I., 1956. The network model of porous media.
33
ACCEPTED MANUSCRIPT
Feng, F. and Akkutlu, I. Y. 2015. November. Flow of Hydrocarbons in Nanocapillary: A Non-Equilibrium Molecular Dynamics Study. In SPE Asia Pacific Unconventional Resources Conference and Exhibition. Society of Petroleum Engineers. Fichman, M. and Hetsroni, G. 2005. Viscosity and slip velocity in gas flow in microchannels. Physics of Fluids (1994-present), 17(12), 123102.
RI PT
Florence, F.A., Rushing, J., Newsham, K.E. and Blasingame, T. A. 2007. January. Improved permeability prediction relations for low permeability sands. In Rocky Mountain Oil and Gas Technology Symposium. Society of Petroleum Engineers. Fung, L. S. and Dogru, A. H. (2008). Parallel unstructured-solver methods for simulation of complex giant reservoirs. SPE Journal, 13(04), 440-446.
SC
Geist, A. 1994. PVM: Parallel virtual machine: a users' guide and tutorial for networked parallel computing. MIT press.
M AN U
Gong, X., Tian, Y., McVay, D. A., Ayers, W. B., and Lee, J. 2013. Assessment of Eagle Ford shale oil and gas resources. Presented at the SPE Unconventional Resources Conference, Calgary, Alberta, 5-7 November. SPE-167241-MS. Guo, X., Shen, Y. and He, S. 2015. Quantitative pore characterization and the relationship between pore distributions and organic matter in shale based on Nano-CT image analysis: A case study for a lacustrine shale reservoir in the Triassic Chang 7 member, Ordos Basin, China. Journal of Natural Gas Science and Engineering, 27, pp.1630-1640.
TE D
Gutierrez, M., R. W. Lewis, and I. Masters. 2001. "Petroleum reservoir simulation coupling fluid flow and geomechanics." SPE Reservoir Evaluation and Engineering 4(03): 164-172.
He, X. and Luo, L.-S. 1997. Theory of the lattice Boltzmann method: from the Boltzmann equation to the lattice Boltzmann equation. Physical Review E 56 (6): 681-6817.
EP
He, X. and Luo, L. S. 1997. Lattice Boltzmann model for the incompressible Navier–Stokes equation. Journal of statistical Physics, 88(3-4), pp.927-944. Holditch, S. A. 2006. Tight gas sands. Journal of Petroleum Technology, 58(06), 86-93.
AC C
Honarpour, M.M., Nagarajan, N.R., Orangi, A., Arasteh, F. and Yao, Z., 2012, January. Characterization of Critical Fluid PVT, Rock, and Rock-Fluid Properties-Impact on Reservoir Performance of Liquid Rich Shales. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Javadpour, F., Fisher, D. and Unsworth, M., 2007. Nanoscale gas flow in shale gas sediments. Journal of Canadian Petroleum Technology, 46(10). Javadpour, F., 2009. Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone). Journal of Canadian Petroleum Technology,48(08), pp.16-21. Jones, F. O. and Owens, W. W. 1980. A laboratory study of low-permeability gas sands. Journal of Petroleum Technology, 32(09), 1-631. Karniadakis, G., Beskok, A. and Aluru, N., Microflows and nanoflows: fundamentals and simulation. 2005.
34
ACCEPTED MANUSCRIPT
Killough, J. E., and Bhogeswara, R. 1991. Simulation of compositional reservoir phenomena on a distributedmemory parallel computer. Journal of Petroleum Technology, 43(11), 1-368. Kim, J. 2010. Sequential methods for coupled geomechanics and multiphase flow. Ph.D. Thesis, Department of Energy Resources Engineering, Stanford University, California, 264 pp.
RI PT
Kim, J., Sonnenthal, E., and Rutqvist, J. 2015. A sequential implicit algorithm of chemo-thermo-poro-mechanics for fractured geothermal reservoirs. Computers and Geosciences, 76: 59-71. Kim, S. H., Pitsch, H., and Boyd, I. D. 2008. Slip velocity and Knudsen layer in the lattice Boltzmann method for microscale flows. Physical Review E 77 (2): 026704. King, D. A. 1980. Surface diffusion of adsorbed species: A Review. Journal of Vacuum Science and Technology, 17(1), 241-247.
SC
King, E. G. Thirty Years of Gas Shale Fracturing: What Have We learned? SPE 133456, presented at the SPE Annual Technical Conference and Exhibition, Florence, Italy, 19-22 September, 2010.
M AN U
King, G.R. and T. Ertekin, State-of-the-Art Modeling for Unconventional Gas Recovery, SPE Formation Evaluation, March, 1991. Klinkenberg, L. J. 1941. The permeability of porous media to liquids and gases. Drilling and Production Practices API-41-200. Knudsen, M. 1909. Molecular effusion and transpiration. Nature, 80, 491-492.
Knudsen, M. 1995. The laws of molecular flow and of inner friction flow of gases through tubes. Journal of Membrane Science, 100(1), 23-25.
TE D
Kohl, T., Evans, K. F., Hopkirk, J., Jung, R. and Rybach, L., 1997. Observation and simulation of non-Darcian flow transients in fractured rock.Water Resources Research, 33(3), pp.407-418. Kuila, U. and Prasad, M. 2013. Specific surface area and pore‐size distribution in clays and shales. Geophysical Prospecting, 61(2), pp.341-362.
EP
Kurihara M, Ouchi H, Inoue T, Yonezawa T, Masuda Y, Dallimore S. R., Collett T. S. 2005. Analysis of the JAPEX/JNOC/GSC et al. Mallik 5L-38 gas hydrate thermal production test through numerical simulation. Geological Survey of Canada, Bulletin 585.
AC C
Ladd, A. J. C. and Verberg, R. 2001. Lattice-Boltzmann simulations of particle-fluid suspensions. Journal of Statistical Physics 104 (516): 1191-1251. Lai, B.T. 2010. Experimental Measurement and Numerical Modeling of High Velocity Non-Darcy Flow Effects in Porous Media, PhD dissertation, Colorado School of Mines, Golden, Colorado. Lai, B. T., J. L. Miskimins, and Y.S. Wu, “Non-Darcy Porous Media Flow According to the Barree and Conway Model: Laboratory and Numerical Modeling Studies,” SPE-122611, SPE Journal, Volume 17, Number 1, pp.70-79, March 2012. Lei, Q. W. Xiong, J. Yuan, and Y.-S. Wu, Behavior of flow through low-permeability reservoirs, SPE 113144, presented at the 2008 SPE Europec/EAGE Annual Conference and Exhibition held in Rome, Italy, 9–12 June 2008.
35
ACCEPTED MANUSCRIPT
Lemmens, H. J., Butcher, A. R., and Botha, P. W. S. K. 2010. FIB/SEM and automated mineralogy for core and cuttings analysis. Presented at the SPE Russian Oil and Gas Conference and Exhibition, Moscow, 26-28 October. SPE-136327-RU. Lemmens, H. J., Butcher, A. R., and Botha, P. W. S. K. 2011. FIB/SEM and SEM/EDX: a new dawn for the SEM in the core lab? Petrophysics 52 (6): 452-456.
RI PT
Li, G., Yao, Y., and Zhang, R., 2016. An improved model for the prediction of liquid loading in gas wells. Journal of Natural Gas Science and Engineering, vol. 32, pp.198-204. Li, J. and Sultan, A.S., 2015, December. Permeability computations of shale gas by the pore-scale monte carlo molecular simulations. In International Petroleum Technology Conference. International Petroleum Technology Conference.
SC
Li, X., Abass, H., Teklu, T. W., and Cui, Q. 2016. A Shale Matrix Imbibition Model-Interplay between Capillary Pressure and Osmotic Pressure. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
M AN U
Lilley, C. R. and Sader, J. E. 2007. Velocity gradient singularity and structure of the velocity profile in the Knudsen layer according to the Boltzmann equation. Physical Review E, 76(2), 026315. Lilley, C. R. and Sader, J. E. August, 2008. Velocity profile in the Knudsen layer according to the Boltzmann equation. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences (Vol. 464, No. 2096, pp. 2015-2035). The Royal Society. Lockerby, D. A., Reese, J. M., Emerson, D. R. et al. 2004. Velocity boundary condition at solid walls in rarefied gas calculations. Physical Review E 70 (1): 017303.
TE D
Lockerby, D. A., Reese, J. M., and Gallis, M. A. 2005. Capturing the Knudsen layer in continuum-fluid models of nonequilibrium gas flows. AIAA journal, 43(6): 1391-1393. Lockerby, D. A. and Reese, J. M. 2008. On the modelling of isothermal gas flows at the microscale. Journal of Fluid Mechanics, 604: 235-261.
EP
Lopez, H. D. 2007. Experimental Analysis and Macroscopic and Pore-level Flow Simulations to Compare NonDarcy Flow Models in Porous Media, PhD dissertation, Colorado School of Mines, Golden, Colorado.
AC C
Loyalka, S. K., Petrellis, N., and Storvick, T. S. 1975. Some numerical results for the BGK model: thermal creep and viscous slip problems with arbitrary accomodation at the surface. Physics of Fluids (19581988), 18(9), 1094-1099. Lu, X.C., Li, F.C. and Watson, A.T. 1995. Adsorption measurements in Devonian shales. Fuel, 74(4), pp.599-603. Luffel, D. L., Hopkins, C. W., and Schettler Jr, P. D. January, 1993. Matrix permeability measurement of gas productive shales. In SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Ma, J., Sanchez, J.P., Wu, K., Couples, G.D. and Jiang, Z. 2014. A pore network model for simulating non-ideal gas flow in micro-and nano-porous materials. Fuel, 116, pp.498-508. Mason, E.A. and Malinauskas, A.P., 1983. Gas transport in porous media: the dusty-gas model (Vol. 17). Elsevier Science Ltd.
36
ACCEPTED MANUSCRIPT
Masuda Y, Naganawa S, Ando S, Sato K. Numerical calculation of gas-production performance from reservoirs containing natural gas hydrates. paper SPE 38291, Proceedings, Western Regional Meeting, Society of Petroleum Engineers, Long Beach, California, 1997.
RI PT
Masuda Y, Konno Y, Iwama H, Kawamura T, Kurihara M, Ouchi H. Improvement of near wellbore permeability by methanol stimulation in a methane hydrate production well. Paper OTC 19433, Proceedings of the Offshore Technology Conference, Houston, Texas, 2008. Maxwell, J. C. 1879. On stresses in rarified gases arising from inequalities of temperature. Philosophical Transactions of the Royal Society of London, 231-256. Mehmani, A., Prodanović, M. and Javadpour, F. 2013. Multiscale, multiphysics network modeling of shale matrix gas flows. Transport in porous media, 99(2), pp.377-390.
SC
Mengal, S.A. and R.A. Wattenbarger. Accounting for adsorbed gas in shale gas reservoirs, SPE 141085, presented at the SPE Middle East Oil and Gas Show and Conference held in Manama, Bahrain, 25–28 September 2011.
M AN U
Minkoff, S.E., Stone, C.M., Bryant, S., Peszynska, M., Wheeler, M.F., 2003. Coupled fluid flow and geomechanical deformation modeling. Journal of Petroleum Sciences and Engineering 38, 37-56. MIT (Massachusetts Institute of Technology), 2011. The Future of Natural Gas – An Interdisciplinary MIT Study, (Updated from 2010 version). Moridis, G. J., Kowalsky, M. B., and Pruess, K. 2007. Depressurization-induced gas production from class-1 hydrate deposits. SPE Reservoir Evaluation and Engineering, 10(05), 458-481.
TE D
Moridis, G. J., Reagan, M. T., Kim, S. J., Seol, Y., and Zhang, K. 2009. Evaluation of the gas production potential of marine hydrate deposits in the Ulleung Basin of the Korean East Sea. SPE Journal, 14(04), 759-781. Moridis, G.J., T.A. Blasingame, and C.M. Freeman. Analysis of Mechanisms of Flow in Fractured Tight-Gas and Shale-Gas Reservoirs. SPE 139250, presented at the SPE Latin American & Caribbean Petroleum Engineering Conference held in Lima, Peru, 1–3 December 2010.
EP
Moridis, G. J., Silpngarmlert, S., Reagan, M. T., Collett, T., and Zhang, K. 2011. Gas production from a cold, stratigraphically-bounded gas hydrate deposit at the Mount Elbert Gas Hydrate Stratigraphic Test Well, Alaska North Slope: Implications of uncertainties. Marine and Petroleum Geology, 28(2), 517-534.
AC C
Nelson, P. H. 2009. Pore-throat sizes in sandstones, tight sandstones, and shales. AAPG Bulletin, 93(3), 329-340. Nguyen T. S. Description of the computer code FRACON. In: Stephansson, O., Jing, L., Tsang C-F, editors. Coupled thermo-hydro-mechanical processes of fractured media, vol. 79. Elsevier: Developments in Geotechnical Engineering, 1996. p. 539-44. Noorishad J, Tsang C-F, Witherspoon PA. 1984. Coupled thermal-hydraulic-mechanical phenomena in saturated fractured porous rocks: numerical approach. J Geophys Res., 89:10365-73. Noorishad J, Tsang C-F. ROCMAS-simulator: a thermohydro-mechanical computer code. In: Stephansson O, Jing L, Tsang C-F, editors. Coupled Thermo-hydro-mechanical processes of fractured media, vol. 79. Elsevier: Developments in Geotechnical Engineering, 1996. p. 551-8.
37
ACCEPTED MANUSCRIPT
O’Hare, L., Lockerby, D. A., Reese, J. M., and Emerson, D. R. 2007. Near-wall effects in rarefied gas microflows: some modern hydrodynamic approaches. International Journal of Heat and Fluid Flow, 28(1), 3743.
RI PT
Ohnishi, Y. and A. Kobayashi, THAMES. In: Stephanson, O., Jing, L. and Tsang, C. F. editors. Coupled thermohydro-mechanical processes of fractured media: Mathematical and Experimental Studies, vol. 79. Elsevier: Developments in Geotechnical Engineering, 1997, p. 545-549. Ohwada, T. and Xu, K. 2004. The kinetic scheme for the full-Burnett equations. Journal of Computational Physics, 201(1), 315-332. Okabe, H. and Blunt, M. J. 2004. Prediction of permeability for porous media reconstructed using multiple-point statistics. Physical Review E, 70(6), p.066135.
SC
Olivella, S., Carrera, J., Gens, A., Alonso, E.E., 1994. Nonisothermal multiphase flow of brine and gas through saline media. Transport in Porous Media 15, 271-293.
M AN U
Øren, P.E. and Bakke, S. 2002. Process based reconstruction of sandstones and prediction of transport properties. Transport in Porous Media, 46(2-3), pp.311-343. Pan, L. S., Liu, G. R., and Lam, K. Y. 1999. Determination of slip coefficient for rarefied gas flows using direct simulation Monte Carlo. Journal of Micromechanics and Microengineering, 9(1), 89. Patil, V. E., Meeuwissen, J., van den Broeke, L. J., and Keurentjes, J. T. 2006. Permeation of supercritical fluids across polymeric and inorganic membranes. The Journal of Supercritical Fluids, 37(3), pp.367-374. Pollard, W. G. and Present, R. D. 1948. On gaseous self-diffusion in long capillary tubes. Physical Review, 73(7), 762.
TE D
Present, R. D. 1958. Kinatic Theory of Gases.
Rodgers, D. P. June, 1985. Improvements in multiprocessor system design. In ACM SIGARCH Computer Architecture News (Vol. 13, No. 3, pp. 225-231). IEEE Computer Society Press.
EP
Rutqvist, J., 2011. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations.Computers & Geosciences, 37(6), pp.739-750.
AC C
Rutqvist J., Börgesson L., Chijimatsu M., Kobayashi A., Nguyen T. S., Jing L., Noorishad J., and Tsang C.-F. 2001. Thermohydromechanics of partially saturated geological media - Governing equations and formulation of four finite element models. Int. J. Rock Mech. & Min. Sci. 38: 105-127. Rutqvist, J., Zheng, L., Chen, F, Liu H.-H, and Birkholzer, J. 2014. Modeling of Coupled ThermoHydroMechanical Processes with Links to Geochemistry Associated with Bentonite-Backfilled Repository Tunnels in Clay Formations. Rock Mechanics and Rock Engineering, 47: 167-186. Shabro, V., C. Torres-Verdín, and F. Javadpour, Numerical simulation of shale-gas production: from pore-scale modeling of slip-flow, knudsen diffusion, and langmuir desorption to reservoir modeling of compressible fluid, SPE 144355, presented at the SPE North American Unconventional Gas Conference and Exhibition held in The Woodlands, Texas, USA, 14–16 June 2011. Sakhaee-Pour, A. and S. L. Bryant. 2011. Gas permeability of shale. Presented at the SPE Annual Technical Conference and Exhibition held in Denver, Colorado, USA, 30 October–2 November. SPE 146944.
38
ACCEPTED MANUSCRIPT
Sakhaee-Pour, A. and Bryant, S., 2012. Gas permeability of shale. SPE Reservoir Evaluation and Engineering, 15 (04): 401-409. Sampath, K. and Keighin, C. W. 1982. Factors affecting gas slippage in tight sandstones of cretaceous age in the Uinta basin. Journal of Petroleum Technology, 34(11), 2-715.
RI PT
Schettler Jr, P.D. and Parmely, C.R., January, 1991. Contributions to total storage capacity in Devonian shales. In SPE Eastern Regional Meeting. Society of Petroleum Engineers. Settari, Antonin, and Dale A. Walters. 2001. Advances in coupled geomechanical and reservoir modeling with applications to reservoir compaction. SPE Journal 6(03): 334-342. Shan, X., Yuan, X., and Chen, H. 2006. Kinetic theory representation of hydrodynamics: A way beyond the Navier-Stokes equation. Journal of Fluid Mechanics 550: 413-441.
SC
Silin, D. and Patzek, T., 2006. Pore space morphology analysis using maximal inscribed spheres. Physica A: Statistical mechanics and its applications, 371(2), pp.336-360.
M AN U
Silin, D. and T. Kneafsey, Gas shale: from nanometer-scale observations to well modeling, CSUG/SPE 149489, presented at the Canadian Unconventional Resources Conference, Calgary, Alberta, Canada, 15-17 November, 2011. Silin, D., 2014. Digital rock studies of tight porous media.
Sisk, C., Diaz, E., Walls, J. et al. 2010. 3D visualization and classification of pore structure and pore filling in gas shales. Presented at the SPE Annual Technical Conference and Exhibition, Florence, 19-22 September. SPE-134582-MS.
TE D
Soeder, D. J. 1988. Porosity and permeability of Eastern Devonian Gas Shale. SPE Formation Evaluation, 3(01): 116-124. Soeder, D. J. and P. L. Randolph, 1987. Porosity, permeability, and pore structure of the tight mesaverde sandstone, Piceance Basin, Colorado, SPEFE 2(6):129–136. SPE 13134, Trans., AIME, 283, 1987.
EP
Sun, Q. Y., Zhang, R. L., Zhang, Y. Z., and Chen, X. Q. 2010. Pressure Transient Analysis of the Low Permeability Composite Reservoir with Threshold Pressure Gradient. In 2010 Asia-Pacific Power and Energy Engineering Conference, pp. 1-4. IEEE.
AC C
Sun, H., Chawathe, A., Hoteit, H., Shi, X., and Li, L. 2015. Understanding Shale Gas Flow Behavior Using Numerical Simulation. SPE Journal, 20(01): 142-154. Szalmas, L. 2007. Slip on curved boundaries in the lattice-Boltzmann model. International Journal of Modern Physics C 18 (1): 15-24. Tang, G., Zhang, Y., Gu, X. et al. 2008. Lattice Boltzmann modeling Knudsen layer effect in non-equilibrium flows. Europhysics Letters 83 (4): 40008. Taron, J., D. Elsworth and K. B. Min, 2009. Numerical simulation of thermal-hydrologicmechanical-chemical processes in deformable, fractured porous media, International Journal of Rock Mechanics Mining Sciences, 46, 842-854. Thomas, H. R., Sansom, M. B., 1995. Fully coupled analaysis of heat, moisture, and air transfer in unsaturated soil. Journal of Engineering Mechanics ASCE 121, 392-405.
39
ACCEPTED MANUSCRIPT
Tinni, A., Fathi, E., Agarwal, R., Sondergeld, C. H., Akkutlu, I. Y., and Rai, C. S. January, 2012. Shale permeability measurements on plugs and crushed samples. In SPE Canadian Unconventional Resources Conference. Society of Petroleum Engineers.
RI PT
Tran, David, Long Nghiem, and Lloyd Buchanan. 2009. Aspects of coupling between petroleum reservoir flow and geomechanics. In 43rd US Rock Mechanics Symposium & 4th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association. Wang, F. P., R. M. Reed, J. A. Jackson, and K. G. Jackson. Pore networks and fluid flow in gas shales, SPE 124253, presented at the 2009 SPE Annual Technical Conference and Exhibition held in New Orleans, Louisiana, USA, 4–7 October 2009.
SC
Wang, K., Liu, H., and Chen, Z. 2015. A scalable parallel black oil simulator on distributed memory parallel computers. Journal of Computational Physics, 301, 19-34. Wang, L., 2014. Simulation of slip flow and phase change in nanopores. PhD dissertation.
M AN U
Wang, S., Xiong, Y., Winterfeld, P., Zhang, K., and Wu, Y. S. 2014. Parallel Simulation of ThermalHydrological-Mechanic (THM) Processes in Geothermal Reservoirs. In Proceeding of 2014 Stanford Geothermal Workshop. Wang, S., Winterfeld, P. H., & Wu, Y. S. 2015. An Efficient Adaptive Nonlinearity Elimination Preconditioned Inexact Newton Method for Parallel Simulation of Thermal-Hydraulic-Mechanical Processes in Fractured Reservoirs. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers.
TE D
Wang, S., Huang, Z., Wu, Y. S., Winterfeld, P. H., and Zerpa, L. E. 2016. A semi-analytical correlation of thermal-hydraulic-mechanical behavior of fractures and its application to modeling reservoir scale cold water injection problems in enhanced geothermal reservoirs. Geothermics, 64, 81-95.
Wang, W., Kolditz, O., 2007. Object-oriented finite element analysis of thermo-hydro-mechanical (THM) problems in porous media. International Journal Numerical Methods Engineering 69, 162-201.
EP
Wheeler, J. A. and Smith, R. A. 1990. Reservoir simulation on a hypercube. SPE Reservoir Engineering, 5(04), 544-548.
AC C
Wu, K., Nunan, N., Crawford, J.W., Young, I.M. and Ritz, K., 2004. An efficient Markov chain model for the simulation of heterogeneous soil structure. Soil Science Society of America Journal, 68(2), pp.346-351. Wu, K., Van Dijke, M. I., Couples, G. D., Jiang, Z., Ma, J., Sorbie, K. S., Crawford, J., Young, I. and Zhang, X., 2006. 3D stochastic modelling of heterogeneous porous media–applications to reservoir rocks. Transport in Porous Media, 65(3), pp.443-467. Wu, L., and Bogy, D. B. 2003. New first and second order slip models for the compressible Reynolds equation. Journal of Tribology, 125(3): 558-561. Wu, L. 2008. A slip model for rarefied gas flows at arbitrary Knudsen number. Applied Physics Letters, 93(25): 253103. Wu, Y-S., Fakcharoenphol, P., Zhang, R. et al. 2010. Non-Darcy displacement in linear composite and radial flow porous media, in SPE EUROPEC/EAGE Annual Conference and Exhibition, Society of Petroleum Engineers. 40
ACCEPTED MANUSCRIPT
Wu, Y-S., Chen, Z., Kazemi, H., Yin, X., Pruess, K., Oldenburg, C., Winterfeld, P., Zhang, R. 2014. Simulation of Coupled Processes of Flow, Transport, and Storage of CO2 in Saline Aquifers, Tech. rep., Trustees of The Colorado School of Mines. Wu, Y.-S. 2016. Multiphase Fluid Flow in Porous and Fractured Reservoirs. Gulf Professional Publishing.
RI PT
Wu, Y.-S. and P. Fackahroenphol, A unified mathematical model for unconventional reservoir simulation, SPE 142884, submitted for the SPE EUROPEC Conference, 23-26 May 2011 in Vienna, Austria, 2011. Wu, Y.-S., Li, J., Ding, D., Wang, C., and Di, Y. 2014. A generalized framework model for the simulation of gas production in unconventional gas reservoirs. SPE Journal, 19(05), 845-857. Wu, Y.-S. and G. Qin, 2009. A general numerical approach for modeling multiphase flow and transport in fractured porous media,” Communications in Computational Physics, 6(1): 85-108.
SC
Wu, Y.-S. and C. Wang, Transient pressure testing analysis of gas wells in unconventional reservoirs, SPE-SAS312, presented at the 2012 SPE Annual Technical Symposium & Exhibition, Khobar, Saudi Arabia, 811 April 2012.
M AN U
Xiao, F., 2013, Pore-Scale Simulation Frameworks for Flow and Transport in Complex Porous Media. PhD dissertation. Xiong, Y., Fakcharoenphol, P., Winterfeld, P.H., Zhang, R. and Wu, Y-S. 2013. Coupled Geomechanical and Reactive Geochemical Model for Fluid and Heat flow: Application for Enhanced Geothermal Reservoir, in SPE Reservoir Characterization and Simulation Conference and Exhibition, 16-18 September, SPE, Abu Dhabi, UAE.
TE D
Yao, Y., Wu, Y. S., Zhang, R. 2012. The Transient Flow Analysis of Fluid in a Fractal, Double-porosity Reservoir. Transport in porous media, 94(1), 175-187. Yin, S., M. B. Dusseault, and L. Rothenburg. 2011. Coupled THMC modeling of CO2 injection by finite element methods. Journal of Petroleum Science and Engineering 80(1): 53-60. 90.
EP
Yin, X., Koch, D. L., and Verberg, R. 2006. Lattice-Boltzmann method for simulating spherical bubbles with no tangential stress boundary conditions, Physical Review E 73 (2): 026301.
AC C
Zeng, B., L. Cheng, and F. Hao, Experiment and mechanism analysis on threshold pressure gradient with different fluids, SPE 140678, Nigeria Annual International Conference and Exhibition, Tinapa - Calabar, Nigeria, 31 July - 7 August 2010. Zhang, H., Bai, B., Song, K. et al. 2012. Shale gas hydraulic flow unit identification based on SEM-FIB tomography. Presented at the SPE Annual Technical Conference and Exhibition, San Antonio, Texas, 810 October. SPE-160143-MS. Zhang, K. N., Wu, Y.-S., Ding, C., Pruess, K., and Elmroth, E. January, 2001. Parallel computing techniques for large-scale reservoir simulation of multi-component and multiphase fluid flow. In SPE Reservoir Simulation Symposium. Society of Petroleum Engineers. Zhang, P., Hu, L., Meegoda, J.N. and Gao, S., 2015. Micro/nano-pore network analysis of gas flow in shale matrix. Scientific reports, 5.
41
ACCEPTED MANUSCRIPT
Zhang, R, P. H. Winterfeld, X. Yin, Y. Xiong, and Y-S. Wu. 2015. Sequentially coupled THMC model for CO2 geological sequestration into a 2D heterogeneous saline aquifer. Journal of Natural Gas Science and Engineering 27: 579-615. Zhang, R., Lu, N., and Wu, Y. S., 2012a. Efficiency of a community-scale borehole thermal energy storage technique for solar thermal energy. Proceedings of the GeoCongress.
RI PT
Zhang, R., Yin, X., Winterfeld, P.H. and Wu, Y-S. 2012b. A Fully Coupled Model of Nonisothermal Multiphase Flow, Geomechanics, and Chemistry during CO2 Sequestration in Brine Aquifers, Proceedings of the TOUGH Symposium, pp.838-848.
SC
Zhang, R., Yin, X., Wu, Y-S. and Winterfeld, P.H. 2012c. A Fully Coupled Model of Nonisothermal Multiphase Flow, Solute Transport and Reactive Chemistry in Porous Media, in SPE Annual Technical Conference and Exhibition, 8-10 October, San Antonio, Texas, USA. Zhang, R. 2013. Numerical Simulation of Thermal Hydrological Mechanical Chemical Processes during CO2 Geological Sequestration, PhD thesis, Colorado School of Mines.
M AN U
Zhang, R., Wu, Y-S. and Fakcharoenphol, P. 2014. Non-Darcy Displacement in Linear Composite and Radial Aquifer during CO2 Sequestration, International Journal of Oil, Gas and Coal Technology, Vol. 7, No. 3, pp.244-262. Zhang, R. and Wu, Y.S., 2016a. Sequentially Coupled Model for Multiphase Flow, Mean Stress, Reactive Solute Transport with Kinetic chemical reactions for CO2 Geological Sequestration and EGS System, Journal of Porous Media, Vol. 19 (11), pp. 1101-1021. DOI: 10.1615/JPorMedia.v19.i11.60.
TE D
Zhang, R., Xiong, Y., Winterfeld, P. H., Yin, X. and Wu, Y.-S. 2016b. A Novel Computational Framework for Thermal - Hydrological - Mechanical - Chemical Processes of CO2 Geological Sequestration into a Layered Saline Aquifer and a Naturally Fractured Enhanced Geothermal System. Greenhouse Gas: Science and Technology, Vol. 6 (3), pp. 370-400. doi: 10.1002/ghg.1571.
EP
Zhang, R., Yin, X., Winterfeld, P.H. and Wu, Y-S. 2016c. A Fully Coupled Thermal - Hydrological - Mechanical - Chemical Model for CO2 Geological Sequestration, Journal of Natural Gas Science and Engineering, Vol. 28, pp.280-304. Zhang, R., Shan, X., and Chen, H. 2006. Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation. Physical Review E 74 (4): 046703.
AC C
Zhang, W.-M., Meng, G., and Wei, X. 2012. A review on slip models for gas microflows. Microfluidics and Nanofluidics, 13(6), 845-882. Zhang, Y., Qin, R., and Emerson, D. R. 2005. Lattice Boltzmann simulation of rarefied gas flow in microchannels. Physical Review E 71 (4): 047702. Zhao, X., 1991. An Experimental Study of Methane Diffusion in Coal Using Transient Approach, Ph.D. Dissertation, Department of Mining and Geological Engineering, The University of Arizona. Zhao, X., Rui, Z., Liao, X. and Zhang, R. 2015a. The Qualitative and Quantitative Fracture Evaluation Methodology in Shale Gas Reservoir, Journal of Natural Gas Science and Engineering, Vol. 27, pp.486495.
42
ACCEPTED MANUSCRIPT
AC C
EP
TE D
M AN U
SC
RI PT
Zhao, X., Rui, Z., Liao, X. and Zhang, R. 2015b. A Simulation Method for Modified Isochronal Well Testing to Determine Shale Gas Well Productivity, Journal of Natural Gas Science and Engineering, Vol. 27, pp.479-485.
43
ACCEPTED MANUSCRIPT
Highlights 1. Fluid flow and transport mechanisms along with geomechanics in shale and tight gas reservoirs are recapitulated.
RI PT
2. Advances in macroscopic THMC simulation techniques for shale and tight gas reservoirs are reviewed. 3. Pore-scale modeling methods for fluid flow and transport in unconventional reservoir rocks are illustrated and compared. 4. Challenges from unconventional reservoirs and perspectives on developing new simulation
AC C
EP
TE D
M AN U
SC
techniques are discussed.