Triple-Porosity Modelling for the Simulation of Multiscale Flow Mechanisms in Shale Reservoirs

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China School of Mechanical and Chemical Engineering, The University of Western Australia, 35 Stirling Highway, Perth, WA 6009, Australia Department of Energy and Mineral Engineering, G3 Center and Energy Institute, The Pennsylvania State University, University Park, PA 16802, USA State Key Laboratory of Coal Resources and Safe Mining, School of Safety Engineering, China University of Mining and Technology, Xuzhou 221116, China


Introduction
The organic-rich shales provide the basis of unconventional energy production that has recently come to the forefront of the world's energy discussion.A current area of research is the impact of multiscale pore structure on shale gas transport mechanisms and ultimate gas recovery.Many investigators have been inspired to establish suitable models to characterize gas flow in shale encountering great challenges along the way.
Traditionally, dual-porosity models have been used to model naturally fractured reservoirs, which are composed of a porous matrix surrounded by a larger-scale fracture system.The dual-porosity models assume uniform matrix properties throughout the reservoir which may not capture the relevant physics accurately in shale reservoirs, because the shale matrix is comprised of predominantly clay minerals, quartz, pyrite, and organic matter.Studies of porosity and microstructure on high-quality flat surfaces using SEM indicate that shale exhibits a high degree of microstructure, causing low permeability and low porosity [1,2].These pores ranging from nanometer to micrometer size in shales can be attributed with at least four distinct types: inorganic bulk, organic matter, natural microfractures, and hydraulic fractures [3].Kerogen (organic matter) is finely dispersed within the inorganic matrix, which can adsorb and store free gas simultaneously.In order to describe those complex physics, the triple-porosity model has been used to improve dual-porosity models.Ci-qun [4,5] firstly developed a tripleporosity model for radial flow of slightly compressible fluids through a triple-porosity reservoir under pseudo-steady state interporosity flow.Abdassah and Ershaghi [6] built on this model to consider dual matrix systems with different properties.Wang et al. [3] introduced four types of porous media in productive gas-shale systems.Dehghanpour and Shirdel [7] proposed a triple-porosity model composed of two fracture systems and matrix system without considering desorption, diffusion, and slip flow.Alharthy et al. [8] developed a triple-porosity finite-difference (FD) model which combines a several flow mechanisms in nanopores.Huang et al. [9] presented a new triporosity model for shale gas reservoir with consideration of flow mechanisms from nanoscale to macroscale.Peng et al. [10] investigated the effect of deformation and natural fracture on shale gas recovery rate.Cao et al. [11] developed a coupled multiscale model to analyze the impacts of flow regimes, effective stress, and adsorption on the production rate.However, the geomechanical effects on the gas flow are neglected in all these models which are important to determine shale permeability.
Although many previous studies have contributed to better understanding of the complex flow mechanism from nanoscale to macroscale in shale, little work has been done to simultaneously combine the multiscale flow and coupling effect with geomechanics.For example, permeability evolution in shale under the influence of effective stress is not fully addressed.The flow mechanism of these models also is not able to cover the whole range of shale pore scale.In this work, a triple-porosity model comprising kerogen system, matrix system, and natural fracture system was established.This model was valid over the entire range of flow regimes in shale.The permeability evolution fully coupled with effective stress for matrix and fracture was considered simultaneously.To evaluate the impact of shale deposits and reservoir permeability on the production of shale gas, several case studies will be introduced.

Triple-Porosity Systems
2.1.Reservoir Characterization.Understanding the physical multiscale process of gas production from shale reservoirs is important for reservoir evaluation and production optimization.With the scanning electron microscope (SEM) and backscattered electron (BSE) images, many researchers indicate that shales exhibit compositional heterogeneity and variations in pore structures and composition existing on many scales, from the nanoscale to the macroscale [1].Shale gas flows through a network of pores with different diameters ranging from nanometers to micrometers [12].In addition, microstructural features within shales ultimately affect shale's ability to generate, store, and produce gas.
Figure 1 shows a view of shale gas reservoirs in different scales.In the macroscale view, the shale reservoir contains a connected network of large fractures surrounding the large matrix blocks after hydraulic fracturing.The horizontal well is connected to a complex macrofracture network.The shale gas production rate strongly depends on the hydraulic fractures which are highways for gas transfer from shale matrix blocks to the horizontal well.There are also numerous natural fractures existing in the shale matrix in the microscale view.These natural fractures are generally narrow and may enhance permeability locally.These fracture and matrix systems can compose the usual dual-porosity system similar to the conventional natural fracture reservoirs.However, unconventional resources such as shale gas have a complicated pore distribution.The significant difference between pores in conventional and unconventional gas reservoirs is that the number of nanopores is much higher in unconventional shale gas sediments [15] as shown in Figure 1.The SEM image shows that the reservoir matrix is commonly/primarily composed of clay/silica, organic matter/kerogen, and some minerals primarily.Kerogen bulks have a disorderly distribution in matrix and surrounded by inorganic matters.Free gas in shale is stored in the natural fractures and micropores of the matrix system.Moreover, the massive adsorbed gas, which might account for a part of gas storage in gas shale, is stored in the organic nanopores which have a large surface area with strong affinity.Note that the amount of adsorbed gas by the inorganic walls is considered negligible compared to adsorbed gas in kerogen.
Based on the multispace scale of pore distribution, we proposed a triple-porosity system that comprises three contiguous porous media: the kerogen, inorganic matrix, and natural fractures.The organic material in kerogen system mainly consists of micropores (pore lengths less than 2.0 nm) and mesopores (pore lengths between 2 and 50 nm), with an average pore size below 4-5 nm [16].The kerogen system includes organic material and also large interconnected nanopores that provide active sites for gas adsorption.Furthermore, pores in kerogen are on the order of nanometers in size, which is gas transfer channel for the higher mass of gas molecules.Since most kerogen is scattered within inorganic minerals, any other porosity systems in shale should communicate with the kerogen through a few nanopores in the inorganic matrix system.Inorganic slitshaped pore in the matrix is stress sensitive.However, the round-shaped pores in kerogen system can be neglected [14].For the sake of simplicity, we assumed all pores in the matrix are rounded shape.Shale matrix is surrounded by natural fractures which are pathways to connect with hydraulic fractures or the wellbore.The width of natural fracture systems generally is less than 0.05 mm [17].Based on the scale of pores, we assume that mass transferring from kerogen to matrix can be considered as diffusion process and define the gas flow mechanism in the nanopore of the matrix with apparent permeability and viscous flow as the transport mechanism within natural fractures.
The process of gas production in shale gas reservoirs is a combined sequence of gas flow at different length scales.The main gas flows through the natural fractures which feed the hydraulic fractures while they receive flow from the matrix system only.During reservoir depletion, the thermodynamic equilibrium between gas in kerogen and matrix spaces changes.Hence, gas desorbs from the surface of the organic matters and nanopores in the kerogen system.This nonequilibrium process further drives the gas 2 Geofluids molecules to diffuse from the bulk of the organic matters to the surface of the kerogen exposed to the nanopore network [15].The kerogen pore network hydraulically communicates with the inorganic matrix such that mass transport takes place in the one-way sequence during the gas release.Therefore, the gas flow is sequential from one medium to the other in different space scales.In shale gas systems, these multiscale pores form the flow-path network that allows the flow of gas from the kerogen to the wellbore during shale gas production.Each process of gas flow follows its own path and obeys different transport mechanisms at different length scales.The application of the triple-porosity model requires that each porous medium is distributed continuously in space and holds the porous media conditions.

Governing Equations of Gas Flow in Kerogen.
Microscale flow of hydrocarbons through porous media involves various distinct transport mechanisms.The conventional Darcy's law cannot be in general applicable to describe the variety of the relevant flow regimes other than the viscous flow regime.For Knudsen numbers less than 0.01, the use of ideal gas constant in Darcy's law and the assumption of continuum flow remain valid.For Knudsen numbers greater than 0.01, an effective permeability must be computed to compensate for the Knudsen diffusion and/or the slippage flow.As for the pore diameter in kerogen on the scale of a nanometer, gas transport is dominated by collisions between gas molecules and the walls, which is called the Knudsen diffusion.Thus, the gas mass flux in kerogen can be expressed by the Knudsen diffusion equation [18,19]: where D K is the Knudsen diffusion constant, M is the molar mass, R (=8.314 J/mol/K) is the gas constant, r k is the pore radius in kerogen, T is the absolute temperature in Kelvin, and P k is the gas pressure in kerogen.Throughout this paper, kerogen, matrix, and fractures are identified with subscripts k, m, and f, respectively.Under the initial condition of shale gas reservoirs, an equilibrium exists between adsorbed gas and bulk gas in nanopores.Langmuir adsorption can be used to calculate the adsorbed gas amount [20].For the sake of simplicity, the basic premise of this work is that kerogen bulk has to be saturated with gas at a particular pressure before it can liberate gas into the pores.During depressurization development of shale gas reservoirs, compared with surface diffusion, adsorption/desorption is a very quick 3 Geofluids physical process.Hence, the adsorbed gas amount can still be calculated using Langmuir adsorption [21]: where ϕ k is porosity of kerogen in shale, ρ s is shale density, ρ ka is gas density at atmospheric pressure in kerogen system, L a represents the Langmuir volume constant, and L b represents the Langmuir pressure.
The mass transfer from kerogen to matrix can be considered as diffusion through nanopores connecting these two systems.Once the gas molecules from the kerogen desorb to the matrix, the difference of gas concentration between two continua at the interface controls gas diffusion behavior.Thus, the mass transfer rate of gas desorption from kerogen into matrix can be described as [22] where a km is the shape factor, ρ m is the gas density in the matrix system, ρ k is the gas density in the kerogen system, and D km is the diffusion coefficient.
Based on the discussions on pore characteristic of kerogen, we assume Knudsen diffusion as the transport mechanism within the kerogen system.Combining ( 1) and ( 4), the mass balance equation of the kerogen system was obtained: 3. Governing Equations of Gas Flow in Matrix.The conventional Darcy equation fails to fully capture the physics of flow in the nanopore structure of shale matrix.To describe gas flow in ultratight natural porous media, apparent permeability function is adopted which accounts for some complexities in the gas flow.Beskok and Karniadakis [23] developed a unified model for gas flow through microtubes that are valid over the entire range of flow regimes.It is a general expression to capture continuum, transition, and Knudsen flow for the apparent gas permeability of tight porous media.Florence et al. [24] derived the following model, which relates the apparent permeability k a and intrinsic permeability k ∞ : where f K n is correction parameter for non-Darcy flow in nanopores.Knudsen number is defined as the ratio of the gas mean free path (λ) and the pore diameter (2r m ): in which B k is the Boltzmann constant, d g is the effective molecular diameter, and P m is gas pressure in matrix.
It is noted that the Knudsen flow relies only on the Knudsen number and the intrinsic permeability of the porous medium.With regard to the permeability, there are many experimental observations suggesting that the change in effective stress varies the intrinsic permeability of the shale matrix [25,26].Therefore, the intrinsic permeability does not remain constant during shale gas extraction.To investigate the effective stress effect, we consider the intrinsic permeability as absolute permeability instead which is an effective strainbased model.According to our previous work on the effective strain-based absolute permeability model [27,28], the multiscale permeability model for shale matrix is described as where the effective strain increment is calculated by where ϕ m0 and k m0 are the initial porosity and permeability of matrix, Δε v is total volumetric strain increment, ΔP m /K s is the changing in compressive strain, K s is the bulk modulus of the shale grains, α is the Biot coefficient, Δε s is the gas sorption-induced volumetric strain increment, and ε L the Langmuir volumetric strain constant representing the volumetric strain at infinite pore pressure.For the sake of simplification, we assume that gas sorption-induced strain in the kerogen system results in the volumetric strain of the matrix system.This permeability model considers the principal controlling factors, including the mechanical deformationinduced pore volume change (first term in the right-hand side of 8), the gas pressure-induced pore volume change (second term), and the sorption-induced pore volume change (third term).The gas mass in the inorganic matrix exists in free phase.The gas mass balance equation in the matrix can be expressed as where u is the gas viscosity.
The mass exchange between the matrix and fracture is captured by a coupling term Q mf which is similar to (4): where a mf is the shape factor and D mf is the diffusion coefficient.

Governing Equations of Gas Flow in Natural Fracture.
Based on the two-part Hooke model (TPHM) proposed by Liu et al. [29], the fracture aperture b under the condition of compression is defined as where b r is the "hard" part of the fracture aperture or the residual fracture aperture that does not change with stress, b f is the stress-sensitive part, and C f is the fracture compressibility.Then, the porosity and permeability of the fracture system are defined as [30] ∅ Because the most efficient transport mechanism is pressure-driven volume flow, Darcy flow is dominant in fracture networks.The gas mass balance equation in the fracture is given as where ϕ f is the porosity of natural fracture, k f is the permeability of fracture, ϕ f 0 is the initial porosity of natural fracture, k f0 is the initial permeability of natural fracture, P f is the gas pressure in the fracture, ρ f is the gas density in the fracture system.

Governing Equations of Mechanics.
In porous elastic media, such as coal, exists the interaction between the interstitial fluid and coal deformation.The effective stress concept introduced by Terzaghi [31] and refined by Biot [32] points out that the pore pressure helps counteract the mechanical stress carried through grain-to-grain contact.Considering the coal swelling/shrinkage stress σ s induced by gas absorption/adsorption [33], the effective stress equation can be expressed as where σ ij is the total stress (positive in compression), σ ij ′ the effective stress, p the pore pressure, δ ij the Kronecker delta (δ ij = 1 as i = j and 0 in other cases), and σ s the adsorption swelling stress.The parameter α is called the Biot coefficient; it relates the volume of fluid added (or lost) in a porous material element to the volumetric change of the same element.
According to the theory of continuum mechanics, the combination of equilibrium equation with the constitutive equations for the homogeneous, isotropic, and elastic media finds the Navier-type equation [33], where u i is the component of displacement in the i direction, G the shear modulus, ν the Poisson ratio, and f i the component of body force in the i direction.Here, the pore pressure is the gas pressure in the fracture system and δ and β are the Biot coefficient.The Biot coefficient can be written as where K is the bulk modulus, K s is the grain elastic modulus, and a is the uniform spacing between fractures defining the edge dimension of the REV cubic matrix.
The coupling relationship between fluid flow and deformation in shale was built based on the equilibrium equation that describes an equilibrium state among deformation, gas flow-induced friction force, and swelling-induced body force.

Model Validation
This paper validates the numerical model in the context of the process of performing a simulation study.History matching between simulated results and field data for a horizontal well in Barnett Shale is discussed [34].In this case, the well was stimulated by a multistage fracturing with a single, perforated interval for each stage.The derived governing equations for the gas flow in kerogen, matrix, and fracture system are a set of nonlinear partial differential equations (PDE) with the second order in space and the first order in time.Such a complete set of coupled equations is coupled into the interfaces of COMSOL Multiphysics and solved using the powerful desktop computer.The reservoir is assumed to be homogeneous with a volume of 1000 × 500 × 91.4 m.The model contains 28 hydraulic fractures with 30.5 m spacing evenly.The half-length of all fractures is 47.2 m.Detailed reservoir information of the Barnett Shale and parameters used in simulations are all listed in Table 1 [9,22,[35][36][37].
The history matching of field data with the M3 model is presented in Figure 2. It shows a reasonable match between the numerical simulation results and the actual field gas flow data.Compared to the simulated result with single porosity model from Yu and Sepehrnoori [37], the M3 model gives a better agreement with field date over the first 1 year.The simulated result of the well production for 30 years is also shown in Figure 2. It shows a typical production decline curve that production in the initial several years declines hyperbolically and the production decline levels off reflecting an exponential decline rate at some point.
Simulations are performed with the M3 model in which all fractures, inorganic, and organic porosity systems are allowed to flow among themselves and between different porosity types.This is different from conventional single or dual porosity/permeability models which are not sufficient to describe the complex physics of shale gas.A comparison between M3 model and conventional models is shown in Figure 3.The single-porosity model consists of matrix system only, while the dual-porosity model consists of matrix and fracture systems.As can be seen, evidence gap exists between the M3 model and conventional models in the first 1 year and 5 Geofluids late period process of gas production markedly.Again, it supports that the simulated result with the M3 model is in better agreement than conventional models that could account for the multispace and multitime process of shale gas flow.

Application Study
For the results of evaluation of the M3 model, we perform sensitivity study of multispace and multitime process and permeability evolution of the gas flow in shale.This paper focuses on contributions of permeability in different scales which are important factors for shale gas production.

Multispace and Multitime Process of Gas Flow in Shale.
As discussed above, gas production from shale reservoir is a physical multiscale process from nanoscale to macroscale. Figure 4 shows the curves of the gas flow rate in kerogen, matrix and fracture, respectively.The transporting behaviors for shale gas are desorption and diffusion in nanoscale pores, flow and diffusion in the matrix, and flow in fractures.As the fractures exhibit very high permeability but very low porosity, a rapid initial outflow of gas appears in fractures.After a short term, the gas flow rate drops rapidly in the fracture and then keeps at a lower level relatively (stage I).Just when fractures are depleted mostly, the gas in the matrix starts feeding the fractures and then flows into the well through fractures (stage II).The gas flow from matrix-dominated total flow lasts until about 10 days before a dramatical decline of gas flow rate occurs.As the gas pressure drops in the matrix, gas starts desorbing from the kerogen pore walls into the nanopores and feeding the matrix then (stage III).This stage is the major recovery stage, contributing 96.8% of cumulative production.It should be noted that free gas in the matrix system is the main source of the gas flow in the period process of gas production.This phenomenon is more evidence for stack column chart as shown in Figure 5.The contribution of gas in matrix systems to cumulative production is nearly 79%.The gas production from the kerogen system keeps increasing throughout the production time that accounts for the 15% of total gas production.

Permeability Evolution in the Process of Gas Production.
Reservoir permeability is a crucial parameter for shale gas production.The evolution permeability in shale is strongly related to complex geomechanical processes such as the transport of gas, adsorption, desorption, changing horizontal stress, and vertical strain.This study investigated the evolution of matrix and fracture permeabilities.Then, the effect of matrix permeability and fracture permeability on shale gas production was investigated.
During the period of gas production, the evolution of matrix permeability is attributable to a number of mechanisms.The intrinsic permeability decreases slightly with the decrease of gas pressure due to the increasing effective stress.However, the permeability increases with the decrease of reservoir pressure under the influence of slippage effect.Considering both the effects of flow regimes acting and effective stress, the apparent permeability of matrix increases rapidly at the primary stage of gas production, as shown in Figure 6.After then, it increases slightly with the decrease of reservoir pressure.The fracture permeability decreases severely under the effect of effective stress, suggesting that fracture permeability is more sensitive to effective stress than the matrix.The impact of permeability evolution on gas production was studied by carrying out three comparative cases.Although the matrix permeability increases 1.94 times after twelve years of production, it results into an increase by 0.002% in cumulative gas production compared to the The shale permeability varies over several orders of magnitude in the different reservoir.Several cases were further studied to determine the controlling effects of permeability on gas well performance.

The Impact of Different Permeability on Production Rate
Curves. Figure 7 compares the gas production for shales with constant initial fracture permeability, but with initial matrix permeability at 10 nD, 100 nD, 1000 nD and 10000 nD.The cumulative production only increases by 1.05 times as matrix permeability is 1000 times greater.Further increasing the matrix permeability will not result in increased production.This result indicates that matrix permeability has a minor effect on gas production.However, if the natural fracture permeability is much greater than the matrix, then gas is being transported through the natural fracture network at a higher rate than gas flowing through the matrix into the fractures.In this case, the gas production is dependent on the matrix permeability.
The fractures, on the other hand, are the discontinuous natural microfractures with high permeability, which are surrounding shale matrix to communicate with hydraulic fractures and the wellbore.Simulations were conducted for a series of shales with initial fracture permeability varying between 500 nD and 10000 nD.It is apparent in Figure 8 that significantly higher production are achieved for wells producing from shales with higher fracture permeabilities.The    7 Geofluids final cumulative production increases by 1.2, 1.8 and 2.1 times for fracture permeability varying from 500 nD to 1000 nD, 5000 nD and 10000 nD, respectively.Fracture permeability controls the pressure drawn down in the natural fractures, which in turn drives the rate of gas transfer from the matrix.Although free gas in matrix contributes to the main source of gas production, fracture permeability determines the gas production behavior.It indicates that high permeability fracture networks in matrix system accelerate production by providing high conductivity channels for the flow through a reservoir.For unconventional and tight reservoirs with extreme low-matrix permeability, natural fractures have the potential to play a crucial role in the gas production.

Conclusion
The gas production in shale reservoirs is a combined sequence of gas flow mechanisms at different length scales.Therefore, a triple-porosity model (M3 model) was proposed for understanding the complex flow mechanisms occurring   8 Geofluids in these reservoirs.The M3 model capturing multiple pore scales and flow consists of three contiguous porous media: the kerogen, inorganic matrix, and natural fractures.In shale gas systems, these multiscale pores form the flowpath network that allows the flow of gas from the kerogen to wellbore during shale gas production.Each process of gas flow follows its own path and obeys different transport mechanisms at different length scales.We performed history matching of field production data from Barnett Shale.
The results indicate that the M3 model produces better performance than the conventional dual porosity/permeability models.Based on the application study, the following conclusions can be drawn: (1) At an early stage of gas production, free gas in the fracture system contributes to the main source of gas flow.Next, gas in the matrix starts feeding the fractures and flows into the well through fractures.Finally, gas starts desorbing from the pore walls into the nanopores and feeding the matrix.The contribution of gas in matrix systems to cumulative production is nearly 80%.Both gas in kerogen and fracture     (2) Reservoir permeability is another crucial parameter influencing shale gas production.The results show that apparent permeability of the matrix system increases with the decreasing reservoir pressure, by considering both effects of flow regimes acting and effective stress.The fracture permeability decreases severely under the effect of effective stress, suggesting that fracture system is more sensitive to effective stress.Gas production has a strong relationship with natural fracture permeability than matrix permeability, which controls the pressure drawn down in the natural fractures which in turn drives the rate of gas transfer from the matrix.This phenomenon suggests that free gas in the matrix contributes to the main source of gas production, while natural fracture permeability determines the gas production behavior.

Figure 1 :
Figure1: Conceptual triple-porosity model for shale gas reservoir at various scales (SEM image is from Ambrose et al.[13], figure in lower right is from Wasaki and Akkutlu[14]).

Figure 2 :
Figure 2: History matching of Barnett Shale with simulation result.

Figure 3 :
Figure 3: Gas flow rate curves to compare the M3 model with conventional models.

Figure 4 :
Figure 4: Gas flow rate curves in different systems.

Figure 5 :
Figure 5: Cumulative production curves in different systems.

Figure 6 :
Figure 6: Evolution of matrix and fracture permeability.

Figure 7 :
Figure 7: Impact of matrix permeability on gas production.

Figure 8 :
Figure 8: Impact of natural fracture permeability on gas production.

Table 1 :
Values of variables used for the simulation.