Study on the Influence of Shale Storage Space Types on Shale Gas Transport.

Shale gas is an important unconventional natural gas resource. Studying the microstructure of shale and the gas transport law is of great significance for the development of shale gas. This paper uses the field emission scanning electron microscope to observe shale samples of the BC shale gas reservoirs in southern China. It is found that there are three types of storage spaces on the micro-nano-scale of shale samples. The storage space can be distributed either in pure organic matter or pure inorganic matter or in both organic matter and inorganic matter. They are called organic storage space, inorganic storage space, and mixed storage space of organic matter and inorganic matter, respectively. According to these types of storage spaces, an ideal conceptual model that reflects various types of storage spaces has been researched and established on the micro-nano-scale. At the same time, the transport mechanisms of slip, diffusion, adsorption, and coupling have been considered, and shale mixed storage space has also been considered in particular. On this basis, a comprehensive equation that can simulate the transport of shale gas in various types of storage spaces is derived. The equation also introduces the proportional parameters of the organic part, fractal characteristics, and water film of the inorganic part in the mixed storage space. Researchers can adjust this parameter to simulate shale gas transportation in different types of storage spaces and then use the finite element method to solve it numerically. This paper analyzes the influence of shale reservoir space types on shale gas transport. The larger the proportion of organic components in the mixed pores, the better the gas transport. The rough fractal dimension of the pores also affects the gas transport. However, when the pore diameter is less than 300 nm, the rough fractal dimension of the pores has a negligible influence on gas transport. For the water film on the inorganic wall surface of mixed pores, the gas transport of the macropore is more sensitive to the change in water film thickness.

Introduction

Shale gas is an important unconventional natural gas resource.[1] With the depletion of traditional energy sources, the industry has begun to pay more attention to unconventional resources.[2] Shale is a fine-grained sedimentary rock (<62.5 μm), which is very dense and usually contains a large amount of organic matter or kerogen.[3] In some organic-rich shale, the content of kerogen or organic matter may account for 40% of the total shale volume. Inorganic minerals and kerogen in shale contain many tiny storage spaces. The small grain size in shale makes the storage spaces in shale extremely dense. Shale gas is produced in source rocks or shale. Because of the low permeability of shale, it is hard for shale gas to be moved. Shale gas is stored in tiny micro–nano reservoir spaces, which are composed of fractures and pores of different sizes. Micro- and nanoscale storage spaces in shale account for the majority. The storage space here refers to the space in shale that can be used to store gas, which includes pores and microfractures. According to the pore size in shale, it can be divided into ultramicropores (pore diameters < 0.7 nm), micropores (0.7–2 nm), mesopores (2–50 nm), and macropores (>50 nm).[4] The pore size is distributed in the range of 10–300 nm. Shale gas is also self-generating and self-storing, often stored in storage spaces in the state of free gas or adsorbed gas. According to the distribution of storage spaces in shale, it can be divided into organic storage spaces and inorganic storage spaces, in which organic storage spaces account for a large part.

For pore characterization, previous studies had used various methods to study the pore characteristics of shale gas, such as scanning electron microscopy (SEM), nano-computed tomography (CT), and nitrogen adsorption.[5−7] Some studies used nano-CT data to construct a three-dimensional pore network; some studies speculated the possible pore structure of shale based on nitrogen adsorption; and some studies investigated the pore size distribution of shale through nuclear magnetic resonance and mercury injection analysis. These studies mainly divided the pores in shale into organic pores and inorganic pores.

At present, there is hardly any literature on organic–inorganic mixed storage spaces.[8] An organic storage space means that the storage space is surrounded by shale organic matter.[9] An inorganic storage space means that the storage space is surrounded by shale inorganic matter. An organic–inorganic mixed storage space means that a part of the storage space is organic and a part is inorganic.

Due to the complex and diverse reservoir spaces of shale gas reservoirs, the state of shale gas is also diverse, which makes the transport mechanism of shale gas very complicated. Understanding the transport law of shale gas in shale reservoir storage spaces is a key scientific issue for accurately evaluating the productivity of gas reservoirs. Many related studies have been conducted on shale gas transport mechanisms. Among them, shale transport mechanisms in shale micro–nano reservoir spaces include viscous flow, slippage, and diffusion.[10] Generally, the Knudsen number (Kn) is used for classifying shale gas transmission types. Kn is defined as the ratio of the average molecular free path to the characteristic length,[11] which represents the severity of the collision of gas molecules with gas molecules and pore walls. When Kn < 0.001, the gas molecular velocity at the pore wall is zero and Darcy’s law is valid. This transport mechanism is also called continuous flow. When 0.001 < Kn < 10, the gas molecular velocity at the pore wall is no longer zero, so the gas flux increases. Darcy’s law is no longer valid, and this transport mechanism is called thin gas transport.[12] Many models have been developed to describe the mechanism of gas transport in shale.[13] Among these models, slip models are used for describing shale gas slippage. This type of model considers the flow mechanism in shale nanopores by adding a slip coefficient to the Darcy formula. The diffusion models are used to characterize the concentration diffusion of shale gas in micro–nano storage spaces, such as the Knudsen diffusion model.[14] Many studies have shown the dissolution and adsorption of shale gas on the surface of organic matter in the shale matrix. The Langmuir adsorption model or the Brunauer–Emmett–Teller (BET) adsorption model is often used to describe the adsorption behavior of shale gas in organic matter.[15]

In the pore channel, the fluid viscous flow rate can be calculated by the Navier–Stokes equation. However, for shale gas, complex mechanisms such as gas volume diffusion, surface diffusion, and adsorption are difficult to be considered in the equation.[16] The modeling of porous media and channels can be well-realized using the Darcy–Brinkman–Stokes (DBS) equation.[17−19] This equation is applicable to both free flow in the channel and seepage in the surrounding porous media. However, this method is difficult to characterize the nonviscous flow in micro–nano pores of shale, and the equations are complicated, so the application range is limited. In addition, the generalized lattice Boltzmann method can be used to simulate gas transport at the microscopic molecular level. However, this method is still under development, and further research is needed to characterize the complex mechanisms of shale gas.[20,21]

The research method in this paper is based on the mass balance equation and motion equations of the fluid. Most similar models are about organic pore type and inorganic pore type of shale. There are no relevant reports on the mixed pore type. Therefore, it is necessary to develop a novel comprehensive model that simulates shale gas transport in this kind of storage space. The research model mainly focuses on the mixed pore reservoir space in shale. The different wall surfaces, stress sensitivity, and irregular cross sections of pores are considered and analyzed in the model, so it can be more feasible to simulate gas transport in the mixed pore types in shale. The study also considered the influence of the water film attached to the inorganic surface of pores. The model can be used for shale pore research and reservoir numerical simulation, providing strong support for the study of shale gas transport mechanisms and development of shale gas reservoirs.

Analysis of Shale Gas Storage Spaces in BC Gas Reservoirs

Preparation of Shale Samples

The study mainly uses scanning electron microscopy to analyze the shale storage space of shale gas reservoirs. The shale samples come from the shale reservoir BC in Hunan, China. BC is the name of a shale gas reservoir.

The steps to prepare the shale sample for SEM observation are as follows. First, the shale sample is cut to a suitable size and polished carefully with sandpaper. Second, the shale sample slices are pasted on the substrate and then placed in an argon-ion polisher. By bombarding the sample surface with an argon-ion beam, the surface of the shale sample becomes flat. Finally, the sample is fixed on the stage with conductive tape for sputtering, and then, the shale sample for SEM observation is completed. Ion polishing (or milling) technology is the process of removing the top amorphous layer on a sample to reveal the pristine and flat sample surface for high-resolution imaging. The ion milling machine uses a high-energy ion gun to bombard the top surface of the sample. The high-energy ions interact with the loosely bounded surface atoms and remove them to reveal an atomic-level clean surface. The reason for polishing or milling is that the surface of the original rock sample is rough because of its dense structure and small pore size. When SEM observation was performed on unpolished rock samples, the pore morphology at the nanometer level could not be clearly observed. In the experiment, as shown in Figure , a copper shield is often used to cover a part of the rock sample, and another part of the rock sample is bombarded with an argon-ion beam. Eventually, this part of the rock sample will form a slope and an ion-polished surface. The two sides are imaged by SEM; the finally obtained SEM images are clear and reliable, and the pore characteristics are more obvious.

Figure 1

Schematic of the principle of the argon-ion polishing technology.

Determination of Shale Material Composition

The mineral composition of shale samples can be determined by scanning electron microscope observation and X-ray diffraction mineral analysis. The brightness and structure of various minerals can preliminarily distinguish mineral types. In the electron microscope images of shale samples, metallic minerals (such as pyrite) have the highest brightness; organic matter has the lowest brightness; and clay minerals, quartz, calcite, dolomite, and other minerals have moderate brightness. The structure of clay minerals is loose, while the structure of quartz and calcite minerals is dense. Energy-dispersive X-ray spectroscopy (EDS) is a useful analytical technique, which can be used for the elemental analysis or chemical characterization of a shale sample. It relies on the interaction of some source of X-ray excitation and a sample. The study uses EDS technology to identify different components in SEM images by analyzing the spectrum of elements. As shown in Figure , kerogen and quartz materials are identified from the SEM images.

Figure 2

EDS spectrum of shale samples. (a) EDS spectrum of kerogen in a shale SEM image. (b) EDS spectrum of quartz in a shale SEM image.

According to X-ray diffraction mineral analysis, the mineral composition of the shale sample is shown in Figure . In shale samples, the clay content is the highest, followed by the quartz content, and the pyrite content is the lowest. In shale clays, the illite clay content is the highest, followed by the I/S interlayer content, and the kaolinite content is the lowest.

Figure 3

Shale rock mineral content analysis.

Types of Shale Storage Spaces

The shale samples observed are from the Longmaxi formation in the study area, with a burial depth ranging from 1700 to 1900 m. There are four core wells and 38 shale samples in the study area, and a total of 350 SEM photos were obtained. The observation instrument is a ZEISS SIGMA field emission scanning electron microscope (FESEM) (and an energy spectrometer).

As shown in Figure , the surface porosity of the pores can be determined by SEM image analysis technology. The pores in the image can be identified by selecting an appropriate threshold. After identifying the pore area, the surface porosity can be calculated using the ratio of the pore area to the area of the SEM image. For different pore types, mineral analysis can be used for determination.

Figure 4

Identification of pores in SEM images using image analysis technology.

By analyzing the 350 SEM images, the surface porosity of Longmaxi shale samples was found to be between 0.59 and 6.53%, with an average of 3.16%. The surface porosity of organic pores extracted by SEM images was in the range of 0.17–0.46%, with an average of 0.33%. The surface porosity of inorganic pores ranged from 1.12 to 8.03%, with an average of 5.03%. The porosity of mixed pores was in the range of 0.12–1.77%, with an average of 1.53%. The kerogen type of shale in this study area is type I, the vitrinite reflectance Ro is between 1.83 and 2.52%, and the organic matter is in the high-over-mature stage. Some SEM observation and analysis data of the shale samples are shown in Table .

Table 1

Observation and Analysis Data of SEM Images of Shale Samples

sample no.	depth (m)	storage space type	area (μm2)	surface porosity (%)
BC-3944	1755.12	organic	297	0.33
BC-3945	1762.36	mixed	1125	1.25
BC-3947	1771.28	inorganic	2106	2.34
BC-3948	1772.67	organic	387	0.43
BC-3950	1773.85	mixed	873	0.97
BC-3952	1779.84	inorganic	981	1.09
BC-3954	1781.46	mixed	1494	1.66
BC-3955	1783.69	organic	189	0.21
BC-3958	1784.87	mixed	1557	1.73
BC-3961	1785.33	organic	171	0.19
BC-3963	1786.31	inorganic	7128	7.92
BC-3964	1790.28	inorganic	4788	5.32
BC-3969	1800.57	organic	324	0.36

The shale micro–nano storage spaces of the BC shale gas reservoir can be divided into three types, namely, pure inorganic storage spaces, pure organic storage spaces, and inorganic–organic mixed storage spaces. Taking the FESEM image of a shale sample with a sampling depth of 1800 meters as an example to show these three types of features, the specific shale sample image is shown in Figure . In Figure a, nanoscale inorganic pores exist in the sample, and the sizes and shapes of these pores are different. Most of the pores have a circular cross-sectional shape. The diameter of these pores is about tens to hundreds of nanometers. In Figure b, nanoscale natural fractures can be observed in the shale sample. Some of these fractures are inorganic fractures, and some of them are organic–inorganic mixed fractures because they are partly through kerogen organic matters. The white strawberry-like part in the sample image is a pyrite cluster, and usually, inorganic pores can be found between these pyrite particles. In Figure c,d, it can be observed that nanoscale organic pores and organic fractures exist in the kerogen of the example. During hydrocarbon generation, shale gas can flow out through fractures or pores in these kerogens.

Figure 5

Field emission scanning electron microscope images of shale samples. (a) Inorganic reservoir space in quartz minerals in shale. (b) Mixed storage space at the boundary of organic matter and pyrite clusters in shale. (c) Organic storage space in an organic matter particle in shale, and the cross section of this storage space is irregular in shape. (d) Organic storage space in an organic matter particle in shale, and the cross section of this storage space is a slit in shape.

Establishment of the Shale Gas Storage Space Geometric Model

Based on the observation and analysis of the shale gas storage space types, the study also combines the shale gas mechanisms, which include shale gas transport mechanisms and gas storage mechanisms, such as shale gas slippage, diffusion, surface diffusion, and desorption/adsorption. A detailed description of the shale gas mechanisms can be found in the fourth section. Based on the above research, a geometric conceptual shale micro–nano (pore diameters in the microscale and nanoscale) storage space model was established (see Figure ).

Figure 6

Schematic diagram of the shale nanoscale pore model.

To facilitate the derivation and numerical solution of the mathematical model, the storage spaces of the shale in this geometric model are simplified as a long cylinder with a length of L and a section diameter of d. This kind of storage space model is prevalent in SEM images of shale samples, which are the circular pores in the images, such as Figure a,c. The storage spaces are composed of two parts, focusing on the organic–inorganic mixed storage spaces. In the storage space model, one part of the storage space is located in the organic matter (see part A in Figure ) and the other part is located in the inorganic part (see part B in Figure ). The inorganic part is mainly composed of clay minerals since the whole-rock analysis of shale samples shows the main component of the inorganic matter in shale samples to be clay. Due to the weak van der Waals force (intermolecular force) between clay mineral molecules and shale gas molecules, it is difficult to adsorb shale gas molecules in this part.[22] The micro–nano storage space model can simulate the transport process of free gas compressed in nanopores, and the adsorbed gas on organic matter pore walls.[23] Under specific temperature and pore pressure conditions, the desorption process and adsorption process of shale gas are in a state of dynamic equilibrium. In the initial state, the shale gas molecules in the pores are in a compressed state. When the pressure at the outlet of the storage spaces decreases, shale gas molecules begin to migrate toward the outlet, and shale gas is produced from the storage spaces in this way. In addition, by changing the ratio of A and B in Figure , the model can simulate all three types of storage spaces of BC shale gas reservoirs. When the ratio of part A is equal to 0, the model can simulate pure organic matter storage spaces. When the ratio of part B is equal to 0, the model can simulate pure inorganic storage spaces.

Establishment of Shale Gas Transport Equations

To study the shale gas production capacity and its influencing factors in the BC shale gas reservoir, a mathematical model describing the transport of shale gas in the storage spaces needs to be derived and established. Aiming at the three types of storage spaces and shale characteristics in the BC shale gas reservoir, the model focuses on the proportion of storage spaces in the organic and inorganic shale parts, along with the impact of free gas in the storage spaces and adsorbed gas on the wall of the storage spaces. In addition, the model also considers the effects of slippage and diffusion of shale gas in the storage spaces.

For this model, the current new contribution is mainly for this type of mixed pores, combining the currently known transport mechanisms reasonably to construct a new comprehensive model, trying to describe the gas transport characteristics in this type of pore.

The model considers various possible influencing factors more comprehensively. In addition to considering the difference in surface composition, the model also considers the influence of shale gas molecules adsorbed on the surface of organic matter on gas transport methods. This includes the coupling of different transport mechanisms. The model also considers the influence of pore roughness and pore deformation.

Model Assumptions

According to the characteristics of shale, the properties of shale gas, and its transport mechanism, the following assumptions are made on the mathematical model: (a) the real gas effect considered; (b) the effect of gravity ignored; (c) gas viscosity changing with pressure; (d) an isothermal flow process; (e) the Langmuir adsorption type of adsorbed gas in shale; (f) gas adsorption on organic matter walls; and (h) single-phase flow.

Mathematical Model of Shale Gas Transport in Shale Storage Space

Establishment of the Motion Equation of Shale Gas

According to literature research, there are both sliding and diffusion effects in shale storage spaces.[24,25] In other words, the total mass flux of shale gas can be written as the sum of the mass flux mslip caused by slippage, the mass flux mdiff caused by diffusion, and the mass flux msurf caused by surface diffusion, which can be expressed by a mathematical formulawhere ωs and ωd are the weighting factors of slippage flow and Knudsen diffusion, respectively. Based on the frequency of molecular collisions, the derived expressions of the weight factor can be written in the following form[26]where Kn reflects the ratio of the mean free path λ of gas to the diameter 2r. The physical expression of Kn is[27]wherewhere T is the temperature, dM is the molecular diameter, and kB is the Boltzmann constant.

In eq , the mass flux mslip caused by slippage can be expressed asThis expression ( is based on a modified Darcy formula, where vs is the slip velocity, ks is the slip permeability, and μg is the gas viscosity in shale storage spaces. ks can be written as a product of the intrinsic permeability kp of storage spaces and the slip factor Fs,[28][28] which is ks = kp × Fs, where the slip factor Fs of shale gas in storage spaces can be expressed as followswhere d is the aperture of shale storage spaces, T is temperature, and M is the molecular weight of shale gas. f is a fraction of molecules striking pore walls, which are diffusely reflected.[28,29]kp is the intrinsic permeability of shale storage spaces, and for cylindrical storage spaces, kp can be derived from the Hagen–Poiseuille equation. In the tough pores, the hydraulic conductance g in terms of the average pore radius r in fractal porous media can be written as follows[30]where Δp is the press difference across the pore. It can be seen from the formula that the conductivity g is proportional to the 2(3Ds2) – β power of the pore radius r. q is the volume flow of shale gas, which can be written in the following formwhere v is the transport velocity of shale gas and A is the cross-sectional area of the pore. Based on the principle of equivalent seepage, kp can be regarded as the Darcy permeability of storage spaces, and by substituting eq into eq , the expression of kp can be given bywhere Ds2 is the fractal dimension of the pore cross section and its value ranges from 1 to 2. Ds3 is the fractal dimension of the inner surface of the pore. β is the parameter, and β = (2Ds3 – 3Ds2)/(2Ds3 – 3), which is shown in Figure . Since the cross section of the real shale pore is not perfectly circular, the irregularity of the circular section needs to be considered.[31]C1 is the proportional coefficient, and its unit is s·kg–1·m2−ζ. The proportional variable C1 can be determined by experiments. The fractal parameters Ds2 and Ds3 are determined using image recognition technology to analyze and calculate the pore cross section and wall surface in the SEM image of the shale rock sample. The commonly used algorithms are differential box-counting and the Hausdorff measure method.

Figure 7

Schematic diagram of rough pores in the shale matrix.

Generally, the initial water saturation in a shale formation is extremely low, often lower than its irreducible water saturation. The water adheres to the inorganic walls of the pores in the form of a water film in the shale. Inorganic minerals are hydrophilic, and water molecules can bind to the surface of inorganic minerals such as clay through hydrogen bonds, electrostatic forces, and van der Waals forces. The microscopic force between the water film and the solid surface can be characterized by the separation pressure model, which consists of three parts: the intermolecular force Π1, the electrostatic force Π2, and the structural force Π3. The three parts have corresponding relationships with the water film thickness hw as follows[32,33]where AH represents the Hamaker constant of the gas–liquid–solid three-phase interaction. ε0 is the permittivity of vacuum, ε0 is the relative permittivity of water, ξ1 and ξ2 are electric potentials, k is the coefficient of structural force, λh is the thickness of the hydration layer, Rh is the relative humidity, and σ is the gas–water surface tension.

It can be seen from Figure that the cross section of the mixed pore is composed of an organic part and an inorganic part. The proportion of the organic part is γ, and θ is the coverage of adsorbed gas molecules on the organic wall. The diameter of adsorbed gas molecules is dM. The inorganic wall is covered with a water film whose thickness is hw.

Figure 8

Cross-section model of the mixed reservoir space.

According to this model and considering the influence of shale gas molecules adsorbed on the surface of organic matter and the water film on the inorganic wall on the pore size, the effective radius of the pore can be written as the following formBased on the theory of solid mechanics, the pores in the rock will be deformed under the action of the pressure of the overlying rock, causing the pore size to become smaller. Based on the above theory and mechanism, it can be derived that the effective pore radius considering the force and deformation of the pore is as followswhere E is Young’s modulus, ν is Poisson’s ratio, σr is the overlying rock stress, and χ is the effective stress coefficient. In addition to the above slip mass flux, the diffusion mass flux of shale gas in shale storage spaces is not negligible. The mass flux mdiff caused by slippage can be given by Fick’s law, which is expressed as follows[34−36]where Ds is the diffusion coefficient of gas in shale storage spaces and cg is the shale gas concentration in shale storage spaces. For single gas flow, the shale gas concentration cg can be replaced by the shale gas density ρg in shale storage spaces. The expression of the diffusion coefficient Ds considering the pore roughness can be written by[37−40]where α is a dimensionless probability factor, L is the length of storage spaces, and δ is the ratio of normalized molecular size dm to local average pore diameter, yielding δ = dm/p. Ds3 is the fractal dimension of the pore surface. The calculation of the probability factor α is quite complex and usually requires an understanding of pore geometry. However, for long cylindrical tubes of length L and diameter d, α is given by d/3L, and the derivation method can be found in the literature.[41]

The main component of shale gas is methane, which accounts for 80–100% of shale gas in reservoirs. In research, shale gas is generally considered as pure methane gas. Under the reservoir conditions, the real gas effect is significant and cannot be ignored in the shale storage space. According to the real gas state equation, the density of gas in the shale reservoir can be derived as shown in the following equation[42]where Z, as the compressibility factor of shale gas in the nanopore, is expressed by reduced temperature T and reduced pressure P. The gas compressibility factor Z can be given by[43−45]where Tc and Pc are critical temperature and pressure for a given gas, respectively.

Some shale gas molecules are usually adsorbed on the walls of shale organic matter storage spaces. During the development of shale gas reservoirs, the pressure difference generated along the storage space causes shale gas molecules to move along the wall. The expression for gas surface diffusion can be written in the following formAmong them, msurf represents the mass flux of shale gas surface diffusion. Ds represents the diffusion coefficient of gas molecules on the wall, which is a function of gas type, temperature, and heat, and its unit is m2/s. To describe the surface diffusion of gas in the micro–nano storage space of shale gas reservoirs under high pressure, and to consider the influence of gas coverage on surface diffusion, the surface diffusion coefficient derived based on the kinetic method is shown as follows[46]where Ds0 represents the surface diffusion coefficient of gas molecules when the wall coverage rate is 0 and vk is the speed ratio, which represents the ratio of obstruction speed to movement speed. H is a function of vk, expressed as follows[46]The diffusion coefficient, Csuf, represents the concentration of adsorbed gas molecules on the wall. The surface adsorption gas concentration formula considering the organic matter wall and molecular coverage is as followswhere NA is Afgar’s Rao constant, dm is the molecular diameter, d is the pore size, and M is the molar mass of shale gas.

Establishment of the State Equation of Shale Gas

In this part, the motion equations of shale gas in three kinds of storage states in the shale micro–nano reservoir spaces of BC gas reservoir are established. This study considers the free gas in the pore space and the adsorbed gas on the wall of organic matter.

Under formation conditions, free gas in shale exists in the reservoir spaces in a compressed state. With the development of shale gas reservoirs, the formation pressure decreases, the gas density in the storage spaces decreases, and shale gas is produced from the shale storage spaces. The mathematical expression of the process can be written aswhere Af is the cross-sectional area of the micro–nano reservoir spaces. The shale gas molecules adsorbed on the organic matter surface will affect the cross-sectional area of the micro–nano reservoir spaces. Considering the influence of the diameter of the shale gas molecules on the cross-sectional area, the expression of the cross-sectional area Af can be expressed asShale gas usually adheres to the surface of shale organic matter. This part of shale gas is usually called adsorbed gas. This study uses the Langmuir adsorption model to characterize this phenomenon. According to Langmuir’s isothermal adsorption theory, the desorption and adsorption process of shale gas will reach an equilibrium state. Based on this theory, this article derives the shale gas mass desorbed from the wall per unit of storage space, and its calculation formula is as followsIn eq , ST is the total adsorption site on the wall of organic matter, which can be calculated by the Langmuir adsorption experiment; NA is Avogadro’s constant; and C is the perimeter of the cross section of the storage spaces. For a storage space with a circular cross section, the perimeter is πd, where d is the diameter of the circular cross section. In the above equation, θ is the fraction adsorbed and refers to the fraction of the total surface area occupied by the adsorbed shale gas molecules. It is a function of pressure. As the pressure in the storage spaces decreases, the adsorbed shale gas molecules gradually desorb from the walls of the organic matter.[47] According to Langmuir’s theory, its expression can be written aswhere K is the Langmuir equilibrium constant for the distribution of shale gas molecules between a surface and the gas phase. p is the pressure in the storage spaces. Because the dissolved gas content in organic matter is negligible, this study does not consider the influence of dissolved gas in the model.

Establishment of the Comprehensive Governing Equation

Based on the established geometric model and model assumptions, and based on the principle of conservation of mass, a time-dependent gas flow model in storage spaces is governed by the mass conservation equation, which can be expressed in the following formwhere ρg refers to the density of shale gas in storage spaces, and it is a function of pore pressure. m represents the total mass flux of shale gas through the pore cross section. qdes denotes the mass of gas desorbed from organic walls. γ is the proportion coefficient of shale organic matter occupied by reservoir spaces, and its expression is γ = A/(A + B). When γ is in the range of 0–1, the storage spaces occupy both organic and inorganic matters of the shale, that is, an organic–inorganic hybrid storage space. When γ is equal to 1, all of the storage spaces are pure organic storage spaces. When γ is 0, all of the storage spaces are pure inorganic storage spaces.

By substituting eqs –24 into eq , we can obtain the comprehensive governing equation for shale gas flow in storage spaces as followsFor the convenience of numerical simulation, eq can be simplified to the following formwhere ∇ is the nabla operator, and in this equation, it is ∂/∂x. The gas density ρg and gas viscosity μg are functions of pressure. The above equation is the comprehensive governing equation for the micro–nano shale storage spaces of the BC gas reservoir. This equation can be used to predict the production capacity of pure organic matter storage spaces, pure inorganic matter storage spaces, and organic–inorganic matter mixed storage spaces.

Model Parameter Selection and Numerical Solution

Selection of Model Parameters

There are various types of micro–nano pores in shale. The micro–nano pores in shale refer to pores with pore sizes ranging from several micrometers to a few nanometers. These pores are mainly made up of organic and inorganic minerals. Most of these pores or fractures are composed entirely of organic or inorganic minerals, but some of them are not. They are composed of a certain proportion of organic and inorganic minerals. To study the effect of pore composition on shale gas flow at a micro–nano-scale, our team designed a series of numerical simulation experiments based on the previously established mathematical model and the data of the BC shale gas reservoir. The experiments assume a micro–nano pore in the shale formation, and with the exploitation of the shale gas reservoir, the gas in this pore is gradually produced. The experiments set the diameter of these pores to 10 nm and the pressure at the outlet of the pores to 19 MPa, and the initial pressure in the pore space is 20 MPa. The pore composition is set by the organic proportion parameter γ, and the value of γ ranges from 0 to 1, which indicates the proportion of organic components of pores. When the parameter γ is set to 1, it indicates that the pores are made entirely of organic matter, and when γ is set to 0, it indicates that the pores are made entirely of inorganic mineral.

The parameters required for the mathematical model and numerical simulation come from the shale reservoir in Hunan, China. Some of the data come from actual core sampling tests and experiments, such as the Langmuir pressure constant, the Langmuir volume constant, reservoir temperature, and initial pressure. Some data are from the literature. These parameters are shown in Table .

Table 2

Parameters Used in the Proposed Model for Numerical Simulation

parameters	value
reservoir temperature, T (K)	363.15
molecular weight, M (kg/mol)	1.68 × 10–3
gas diffusion coefficient, Dg (m2/s)	2.18 × 10–5
universal gas constant, R (J/(mol·K))	8.314
nanopore diameter, d (m)	1 × 10–8
molecular diameter, dm (m)	5 × 10–10
Langmuir pressure constant, PL (Pa)	3.45 × 106
Langmuir volume constant, VL (m3/kg)	8.8 × 10–4
initial pressure, pi (Pa)	1.5 × 10–7
outlet pressure, pout (Pa)	1.0 × 107
Avogadro’s contest, NA (1/mol)	6.02 × 1023
total adsorption site, ST (1/m2)	4 × 1018
nanopore length, L (m)	3 × 10–7
organic matter ratio, γ (1)	0–1
Langmuir’s constant, K (1/Pa)	1.5 × 10–7
thickness of the hydration layer, λh (m)[47]	1.5 × 10–9
coefficient of structural force, k (N/m2)[47]	1 × 107
difference of electric potentials, Δξ (V)[47]	0.05
relative permittivity of water, ε (1)	81.5
Hamaker’s constant of solid/water/air, AH (J)	1 × 10–20

Numerical Solution of the Model

Since the derived mathematical model is a nonlinear partial differential equation, it is difficult to find its analytical solution. The numerical simulation method is usually used to solve the equation. This research uses a finite element method to solve the mathematical model because it is a powerful and useful tool that has been widely used in many engineering simulation fields. An important advantage of using the finite element method is that the solution obtained by the finite element method is more accurate than that by the finite difference method. Figure shows the flowchart of the finite element solution of this study. First, the boundary conditions and parameters are determined, and then, the weak formula of the mathematical model is derived by integration.

Figure 9

Flowchart for finite element solution.

To use finite element numerical methods to solve mathematical equations, the geometric model of the storage space needs to be discretized into a space discrete model. In this paper, a three-dimensional micro–nano reservoir space is established, and the reservoir space is discretized into several triangular pyramids. As shown in Figure , the right endpoint of the discrete model represents the outlet of the micro–nano storage space. The walls of the pore model are set as organic walls and inorganic walls as required. Generally, when the discrete elements of the geometric model are more finely divided, the accuracy of the solution will be higher, but the calculation process will consume more computing resources. When the discrete elements are roughly divided, the accuracy of the results will be affected. This study divides the discrete geometric model of the reservoir space into 56 376 elements or units. The gas pressure will change sharply near the outlet of the model, so it is necessary to optimize the area near the outlet locally. In this way, the instability of the numerical solution of the model can be effectively avoided. Then, based on the discrete numerical model combined with the weak formula of the mathematical model, the matrix equation is obtained. Finally, the computer solves the matrix equation to obtain the numerical solution of the problem.

Figure 10

Discretization of the storage space geometry model.

The following explains the discrete format of finite element numerical simulation. The first step is to convert the comprehensive equation into their equivalent weak formulations. The test function φ is introduced, and the partial differential equations are multiplied with these test functions. Then, they are integrated into the modeling domain and partial integration methods are used to handle weak forms. With the weak formulation, the comprehensive equations can be discretized to obtain numerical equations for numerical solution. Using Green’s first identity, the equation derived from eq can be shown as followswhere V is the modeling domain and n is the normal vector at the boundary. Next, we need to find an approximate solution pa, so that the exact solution p can be approximated by pa. The Galerkin method is usually used for numerical discretization, and it is based on the weighted residual method. In the Galerkin finite element method, the shape function and the test function have the same type. The approximate solution pa is a set of linear combinations of basis functions ψ, which is shown as followswhere n is the number of discrete nodes, ψ is the shape function of pressure, and c is the unknown coefficient that needs to be solved. Substituting eq into eq , the resulting spatial discrete equation is as followsFor transient problems, the finite difference method is often used to discretize the time domain, and the equation is transformed into the following formAfter the comprehensive equation is discretized and boundary conditions are applied, a series of equations can be obtained, written in the matrix form as followsAmong them, c is an n-dimensional unknown vector; the coefficient matrix A is an n-dimensional square matrix or stiffness matrix, whose elements are a. The right side of the equation is the n-dimensional vector b. This series of large equations is solved through iterative methods. After the vector c is obtained, the numerical solutions or approximate solutions of the comprehensive equation can be obtained.

Due to the strongly nonlinear characteristics of the model, it is almost impossible to find the analytical solution. In this study, we use COMSOL Multiphysics finite element software to numerically solve partial differential equations and visualize data.

Model Validation

This study uses the experimental data of Roy et al. to verify the comprehensive transport model. Roy et al. used a film with a thickness of 60 μm and an average pore diameter of 200 nm as the porous medium and measured the flow rate of argon through the film under different pressure differences. The specific experimental parameters are shown in Table below.

Table 3

Parameters Used for Model Validation

parameters	value
length, L (μm)	60
diameter, D (nm)	200
outer pressure, Pout (kPa)	4.8
pressure difference, Δp (kPa)	5–120
temperature, T (K)	300
viscosity, μ (Pa·s)	2.22 × 10–5
molecular mass, M (kg/kmol)	39.948

As shown in Figure , the model of this study is compared with experimental data and other models. The model of this study is highly consistent with the experimental data by Roy et al., with a fitting error of 3.02%. Since this model considers more comprehensive transport mechanisms, its fitting effect is better than other models. Compared with the results of other models, this model comprehensively considers factors such as slippage, Knudsen diffusion, transport coupling, real gas effect, and surface diffusion. The comprehensive consideration of transport factors makes the model more in line with the real situation; therefore, the model is more in agreement with the experimental values. The model verification shows the correctness of the model of this study, and it can be used for simulating the transport behavior of shale gas in shale storage space.

Figure 11

Comparison of model computation results and experimental data.

Due to the limitations of manufacturing technology, current physical models can only design relatively simple microscopic pore structures. To verify the complex mixed types of pores, this study uses molecular simulation methods to verify the model. The molecular simulation method is based on the principle of molecular dynamics and the Monte Carlo method, which simulates the fluid flow at the particle level. However, the large amount of calculation is not conducive to large-scale applications. The research first constructed a complex three-dimensional model of rough mixed pores.

The parameters in the model are mainly from references and experiments. The constants and coefficients in this model, such as the general gas constant, the molecular diameter, the total adsorption site, the adsorption constant, and the fitting coefficient, are mainly from references. Reservoir condition parameters (such as reservoir temperature, pore pressure, and stress conditions) are from field data, which can be obtained by collecting logging data. Fluid characteristic parameters (such as viscosity, diffusion coefficient, density, compressibility factor, etc.) are from fluid experiment data, such as PVT experiments and adsorption test data. The microscopic pore structure parameters of pores, such as fractal dimension and organic matter ratio, mixed pore distribution, etc., are mainly from laboratory analysis of electron microscope images. The image recognition program written can extract the topological and color features of the pores in the images and determine the fractal dimension of the pore outline and the fractal dimension of the inner wall surface, the proportion of organic matter, the distribution of mixed pores, and other parameters.

As shown in Figure , the fractal dimension Ds2 of the model section is 1.37, the fractal dimension Ds3 of the pore wall is 2.6, the ratio of the pore organic matter wall is 0.5, and the pressure difference between the two ends of the pore is 1 MPa. The diameter of the microscopic pore model is 20 nm, and the length of the pore model is 80 nm.

Figure 12

Molecular simulation model of rough mixed pores.

As shown in Figure , the calculation results of the model are compared with the molecular simulation results. The research model fits well with the results of molecular simulation. It shows that the model established has certain reliability for the calculation of gas transport in complex mixed pores of shale, and it can be applied to the flow calculation of complex pores in shale.

Figure 13

Validation of the model using molecular simulation methods.

Results and Discussion

Production Dynamics of Different Pore Types

The three-dimensional numerical simulation of the pore is shown in Figure . The pressure of the pore outlet is low, and the pressure of the left end is high. The pressure of the pore gradually decreases from the left end to the outlet end, and the fluid flow direction is approximately parallel to the x-axis direction.

Figure 14

Shale gas transport in the three-dimensional pore.

After performing numerical simulation experiments, the simulation results are shown in Figure a,b. Figure a represents the gas productions from a micro–nano pore under different organic proportion parameters γ. Figure b represents the cumulative gas productions from a micro–nano pore under different organic proportion parameters γ. γ is the proportion of pore organic matter or organic matter wall. In this series of numerical simulation experiments, the parameter γ is set to 1, 0.75, 0.5, 0.25, and 0. In the initial stage, about 0.7 ns before, the shale gas productions from the micro–nano pore under different parameters γ are basically the same. While in the later stage, it can be easily noticed that these gas production curves start to appear different, and the gap between these curves becomes bigger with the passage of time. Since the left end of the pore model is closed and not connected to other pores, the gas production will eventually rapidly decrease to 0. These changes can be seen more clearly in Figure b, in which the cumulative gas production gradually stabilizes over time. That is, the higher the organic matter content in the micro–nano pores, the higher the gas production from these pores, and the higher the cumulative gas production.

Figure 15

Shale gas production process in shale storage spaces.

It can be seen from Figure that the difference between the curves at the beginning stage is small because the gas produced at this time is mainly from compressed free gas. The amount of gas in inorganic pores and organic pores is relatively sufficient. However, as the mining progresses, due to the absence of adsorbed gas in the inorganic pores, the gas production of the pores decreases rapidly. However, the pores of organic matter contain additional adsorbed gas, and the gas production decreases relatively slowly. Thus, the production curves of different storage space types usually have little difference in the early stage, but it becomes obvious in the later stage.

Sensitivity Analysis of the Model

Influence of Pore Diameter and Pressure on the Gas Production Rate

From Figure , it can be seen that pore pressure and pore size have an important influence on the gas production rate of pores. Generally, the greater the gas pressure in the storage space, the greater the gas production rate in the storage space, and as the pore size of the storage space becomes larger, the gas production rate in the storage space gradually increases. When the pressure of the storage space is constant, the larger the pore size of the storage space, the more obvious the influence of pore size on the gas production rate, that is, the curved surface becomes steeper and steeper. When the pore size of the storage space is constant, the larger the pore pressure in the storage space, the smaller the influence of pore size on the gas production rate, that is, the curved surface becomes flatter. The reason for the phenomenon is that the larger pore size allows the gas to pass through a larger cross-sectional area. The greater pore pressure makes gas molecules more closer. Compared with the influence of storage space size and pore pressure on gas transmission, the influence of the ratio of the organic wall surface and the formation temperature is negligible.

Figure 16

Influence of different pore sizes and pressures on the gas production rate. (a) Effect when the organic matter ratio is 0.2 and the reservoir temperature is 90 °C. (b) Effect when the organic matter ratio is 0.4 and the reservoir temperature is 90 °C. (c) Effect when the organic matter ratio is 0.6 and the reservoir temperature is 120 °C. (d) Effect when the organic matter ratio is 0.8 and the reservoir temperature is 150 °C.

Influence of Rough Pores on the Gas Production Rate

Figure shows the effect of rough pores on shale gas production. Two parameters are used to describe the roughness of pores. Ds2 represents the fractal dimension of the pore cross section. The larger the value, the more irregular the pores. Ds3 represents the fractal dimension of the roughness of the inner surface of the pore. The larger the value, the rougher the inner surface of the pore. It can be seen from the figure that the larger the pore roughness fractal dimension, the lower the gas production rate. When the pore size is small, that is, when the pore size is less than 300 nm, the difference between curves is not obvious. The larger the pore size, the more obvious is this difference between curves. It can be seen from the figure that the rough fractal dimension Ds3 of the inner surface of the pore has a greater impact on gas production, which means that the model is more sensitive to this parameter. The reason for this phenomenon may be that rough pore walls change the direction of movement of molecules and increase the frequency of collisions between shale gas molecules, which mainly affect the viscous flow of shale gas.

Figure 17

Influence of fractal dimension of rough pores on the gas production rate.

Influence of Rock Mechanic Parameters on Gas Production

It can be seen from Figure that the rock mechanic parameters of shale have a certain influence on gas production in the reservoir space. It can be seen from the figure that the smaller the Young modulus, the lower the gas production of the storage space. The smaller the Poisson ratio, the lower the gas production of the storage space. When Young’s modulus of shale is low, gas production is more sensitive to the changes in the Poisson ratio. The reason for this phenomenon is that the smaller Young’s modulus causes a greater longitudinal deformation of the reservoir space under the pressure of the overlying rock, which causes the pore diameter and the gas production rate to become smaller. When the longitudinal deformation of the pore is constant, the greater the Poisson ratio of the rock, the greater the lateral deformation of the pore, that is, the greater the lateral elongation will be. Relatively speaking, the pore size will become larger, and the gas production will also increase.

Figure 18

Influence of rock mechanic parameters on gas production.

Contribution Rate of Different Gas Transport Ways

It can be seen from Figure that at different pore sizes, the contribution rates of different gas transport ways are also different. In the smaller pores or storage space, the gas transport way of the surface diffusion is dominant, while in the larger pores or storage space, the gas transport way of the viscous flow is dominant. The contribution rate of bulk diffusion first increases and then decreases with the increase in pore size. The reason for this phenomenon is that when the pores are small, the space occupied by the gas molecules attached to the pore wall is larger, while the proportion of space used for free gas transport is small. When the pore size is large, the proportion of space occupied by the gas molecules attached to the pore wall is relatively small, while the proportion of space used for free gas transport is relatively large. Therefore, when shale gas transports from small pores to large pores, the primary way of its transport can change at any time.

Figure 19

Contribution rate of different transport methods to gas flow.

It can be seen from Figures and 21 that in the mixed pores, the proportion of organic matter walls has an important influence on the transport of shale gas. For the viscous flow, the greater the proportion of organic matter wall in mixed pores, the lower the contribution rate of viscous flow, and the lesser the amount of gas transported by the viscous flow way. Similarly, the greater the proportion of organic matter walls, the lower the contribution rate of bulk diffusion, and the lesser the volume of gas transported by the bulk diffusion way. However, for surface diffusion, the situation is just the opposite. As shown in Figure , the greater the proportion of organic matter walls in mixed pores, the greater the contribution rate of surface diffusion, and the greater the amount of gas transported by surface diffusion. All of these effects are obvious when the pore size is small (pore size less than 20 nm), and the effect cannot be ignored. When the pore size is large, the influence of mixed pores on gas transport can be appropriately ignored.

Figure 20

Influence of the ratio of organic matter walls in mixed pore on the contribution rate of viscous flow.

Figure 21

Influence of the ratio of organic matter walls in mixed pore on the contribution rate of bulk diffusion.

Figure 22

Influence of the ratio of organic matter walls in mixed pore on the contribution rate of surface diffusion.

Effect of Water Film Thickness on Gas Transport

The thickness of the water film in shale pores also has an important effect on gas transport. We designed several sets of simulation schemes to study the influence of water film thickness on gas transport. The diameter of the pore model is set to 50 nanometers (Figure ).

As shown in Figure , as the thickness of the water film increases, the rate of gas production gradually decreases. This is because the increase in the thickness of the water film will reduce the gas flow area. Under the condition of the same water film thickness, the larger the proportion of organic part in pores, the higher the gas production. This is also because the flow area of the pores has increased. Further observation can be found that the smaller the water film thickness, the weaker the influence of the proportion of organic components on the gas production. When the organic components in the mixed pores are relatively small, the gas production is more sensitive to changes in the thickness of the water film. In addition, it can be found that when the pore pressure is smaller, the distance between the curves is closer, which shows that the gas production is not sensitive to the water film thickness and the ratio of organic matter.

Figure 23

Influence of water film thickness on gas transport under different organic matter ratios.

Figure shows the effect of water film thickness on gas transport under different pore diameters. It can be seen from the figure that as the thickness of the water film increases, the gas productions show a decreasing trend. The larger the aperture, the greater the gas production decrease. This shows that gas production in bigger pores is more sensitive to changes in water film thickness.

Figure 24

Influence of water film thickness on gas transport under different pore diameters.

Transport Simulation of Shale Micro–Nano Pore Network

This comprehensive model can also be used to simulate the transport of shale gas in the micro–nano-scale. First of all, as shown in Figure a, based on the SEM image of the shale sample, artificial intelligence algorithms can be used to identify the storage space, organic matter, and inorganic minerals in the image. As shown in Figure b, the pore structure network in the image is extracted, and the organic and inorganic walls of mixed types of pores are also marked. Then, the extracted pore network skeleton is spatially discretized into several units or elements. It should be noted that local densification is required where the pore diameter is relatively small to avoid the problem of nonconvergence of the solution.

Figure 25

Shale gas production process in shale storage spaces. (a) Image recognition of pore network and shale composition in SEM images. (b) Extraction and discretization of the shale pore network skeleton.

Based on the constructed skeleton and combined with the above-mentioned mechanism model, shale gas transport simulation at the SEM scale can be performed. The pore pressure distribution is shown in Figure a. It can be seen from the figure that the pressure in the narrow pores is higher because of the accumulation of gas molecules. From the velocity profile in Figure b, the main transport channel can be seen for gas flow. Moreover, the gas flows faster in the narrow pores or the corners of the channel, which is caused by the gas convergence effect and the larger pressure gradient.

Figure 26

Distribution of pressure and velocity in shale micro–nano pore network.

Production Simulation of Shale Multistage Fractured Horizontal Wells

Applying this method to multistage fractured horizontal wells in shale gas reservoirs can show the impact of different pore types on the production performance of shale gas wells. Table shows the calculation parameters of shale fractured horizontal wells. Some of the data in these parameters are from actual shale gas wellhead production, and some are from the literature. The established model is based on the geometric size of the actual shale gas well, and then, the established model is discrete in geometric space.

Table 4

Main Input Simulation Parameters of Fractured Wells

parameters	value
reservoir thickness (m)	25
specific gravity of shale gas	0.55
molecular weight of shale gas (g/mol)	16
matrix compatibility (Pa–1)	1.3 × 10–7
initial fracture porosity	0.45
fracture compatibility (Pa–1)	1.5 × 10–7
well bore radius (m)	0.076
initial viscosity (cp)	0.0184
fracture width (m)	0.003
reservoir temperature (K)	316
initial reservoir pressure (MPa)	15.8
initial matrix porosity	0.05

The pressure distribution calculated by this model is shown in Figure . It can be seen that the pressure near the artificial fracture is relatively low, and the artificial fracture is the main channel for shale gas transport. The pores in the shale reservoir are set to different pore types. It can be seen from Figure that the model is very consistent with the actual shale gas production data, indicating that the model is reliable and can be used to predict the production performance of shale gas wells. The difference between different pore types is small, but there are certain differences at the corners of the curve. If the pore type in the reservoir is organic pores, the daily gas production of shale gas wells will be relatively higher. If the pore type is inorganic pores, the daily gas production of shale gas wells will be relatively low, and the final cumulative gas production will also be relatively low. If the reservoir pore types are mixed pores, the gas production of shale gas wells is between the organic and inorganic pore types. In fact, the organic pore types, inorganic pore types, and mixed pore types in shale are mixed in a certain proportion. By determining the ratio of each pore type in the core, the production and ultimate recovery of shale gas wells can be predicted more accurately. As shown in Figure , if the average pore size of the shale reservoir is smaller, this difference will be more obvious. The influence of pore types on the gas production of shale gas wells needs to be taken into consideration (Figure ).

Figure 27

Numerical discretization and simulation of the shale fractured horizontal well model.

Figure 28

Influence of different pore types on the production performance of shale gas wells (the average pore diameter of the shale matrix is 800 nm).

Figure 29

Influence of different pore types on the production performance of shale gas wells (the average pore diameter of the shale matrix is 100 nm).

Conclusions

This study proposes a geometric model for characterizing the storage space types of BC shale gas reservoirs based on scanning electron microscope observations of shale reservoir space types. Then, a comprehensive governing equation for shale gas transport is derived for the types of storage spaces. This equation predicts the effect of different storage space types on shale gas transport. The main conclusions of this paper are as follows.

There are three types of storage spaces in BC shale gas reservoirs: pure organic matter storage spaces, pure inorganic matter storage spaces, and mixed storage spaces.

A comprehensive transport equation was established that can be used to simulate the transport of shale gas in various reservoir spaces. This model can be used to simulate the microscopic pore scale and the macroscopic reservoir scale to study the mechanism and production optimization.

The types of shale storage space will affect the transportation of shale gas. The larger the proportion of organic matter in the mixed pores, the more conducive it will be to gas transport.

Larger Poisson’s ratio and Young’s modulus are conducive to gas transport, and the sensitivity of gas production to Young’s modulus is greater than that of the Poisson’s ratio. The rough fractal dimension of the pores also affects the gas transport, but when the pore diameter is less than 300 nm, the rough fractal dimension of the pores has a negligible influence on gas transport. For the water film on the inorganic wall surface of mixed pores, the gas transport of the macropore is more sensitive to the change in water film thickness.