Shale gas transport behavior considering dynamic changes in effective flow channels

Shale gas transport is influenced by the changes in effective flow channels during production due to an effective stress sensitivity and a sorption layer effect. In this work, shale samples’ core analyses, methane adsorption experiments, and effective stress sensitivity experiments were conducted at in situ conditions. Then, a mathematical modeling was applied to quantitatively determine the real shale gas flow behavior considering the changes in effective flow channels during production. The results show that (1) the gas transport capability in matrix has an obvious tendency of increase in the mid‐late production period; (2) the decline of permeability of inorganic pores/fractures tends to be gentler during production, especially for the smaller flow channels; (3) surface diffusion is dominant in shale matrix throughout production; (4) the contribution of Knudsen diffusion to total gas flux cannot be neglected in flow channels with hydraulic radius of tens of nanometers in the early production period or of tens and hundreds of nanometers in the mid‐late production period; and (5) viscous flow generally dominates the total gas flux in flow channels with hydraulic radius more than 10 nm in the early production period or more than 100 nm in the mid‐late production period. This work is beneficial for an accurate evaluation of shale gas transport during production.

flow of free gas in pores, gas desorption and surface diffusion on pore wall, and configurational diffusion of dissolved gas in organic matter. 12,13 Although the dominate mechanism of gas transport might be different in different sizes of flow channels, the apparent gas permeability of a shale is contributed indeed by all of the above transport mechanisms more or less, which is also demonstrated by numerical simulations of gas flow in shale nanopores with multiple pore geometry. 14,15 Several investigators have presented the unified models on apparent gas permeability of tight rocks. As early as the 1940s, Klinkenberg 16 has put forward a classical mathematical equation to calculate the apparent gas permeability of a rock based on a parameter named slippage factor. After that, the mathematical expressions of slippage factor are determined by some researchers such as Heid et al 17 20 , and Civan. 21 Most of the above studies are based on the fitting results of experimental data, but they might not reflect the gas transport mechanisms. Some researchers, such as Ertekin et al 22 , Thornstenson et al 23 , and Michel et al 24 , believe that the gas transport behavior in rocks is a comprehensive result of Knudsen diffusion nearby pore wall and viscous flow away from the boundary, and establish the more applicable mathematical models to calculate the slippage factor. It is noted that all the above models are based on the classical Klinkenberg model. Some studies, such as Sinha 25 and Kang 26 , indicate that the fitting accuracy of the experimental data is low with a low-pressure or low-permeability conditions using Klinkenberg model. Therefore, some researchers, such as Sakhaee-Pour and Bryan 27 , Civan et al 21 , and Ziarani et al 28 put forward the mathematical models to calculate gas apparent permeability based on the Knudsen number, which is on the basis of Beskok-Karniadakis Equation 29 However, gas transport behavior is actually influenced by a variety of mechanisms, not only one or two mechanisms, especially for low-permeability rocks such as shales, coals, and tight sandstones. Therefore, some researchers further propose the mathematical models coupling several gas transport mechanisms, such as Javadpour et al 30 , Swami et al 31 , and Wu et al. 32 It is necessary to point out that how to comprehensively consider the multiscale pore structure of shale, how to reasonably couple these gas transport mechanisms, and how to consider the influence of microscale effects are still the urgent problems to be overcome. Especially, the microscale effects mentioned above, including effective stress sensitivity, sorption layer effect, and real gas effect, are influenced by the decline of reservoir pressure because not only the permeability of fractures is significantly influenced by formation effective stress but also gas transport behavior in shale nanopores has been demonstrated to be channel-width-dependent. 33 Therefore, it is significant to accurately describe the shale gas transport behavior considering the dynamic changes in effective flow channels sizes.
Multiscale phenomenon is prevalent in various disciplines, which is a hot and difficult issue for current scientific researches. Multigas transport behavior also exists during oil and gas production, but it is not obvious for conventional oil and gas resources. The complexity of shale gas production becomes gradually highlighted as the development of shale gas reservoirs is more and more general in the world, which makes the investigation of multigas transport behavior of shale gas much more important. Therefore, describing the transport behavior of shale gas during production is significant for accurately evaluating the production capacity of a shale gas reservoir. In this work, laboratory measurements and modeling techniques are combined to investigate the multigas transport behavior in a shale gas reservoir during production. The laboratory measurements include HPHT (high pressure and high temperature) methane adsorption experiments and effective stress sensitivity experiments. For F I G U R E 1 Multigas  the modeling techniques, we simulate gas flow behavior in various pore/fracture sizes based on a proposed unified mathematical model, which considers a real gas effect, a sorption layer effect, and an effective stress sensitivity.
Details of the present work are shown as follows: Section 2 describes the shale samples and presents the experimental methods used in this work; Section 3 shows the experimental results; Section 4 establishes the mathematical model of apparent gas permeability considering the dynamic changes in effective flow channels sizes during production; and Section 5 discusses the gas transport behavior during shale gas production based on the calculation results of the proposed mathematical model.

| Samples' description
Shale samples used in this work are from a Longmaxi Formation on the east of the Sichuan Basin in China, which is a Lower Silurian marine shale gas play. The mineralogical composition includes clay, quartz, K-feldspar, plagioclase, calcite, dolomite, and pyrite. The main minerals are quartz and clay minerals, the contents of which are 45.55% and 33.71%, respectively. The average total organic carbon (TOC) of the samples is 3.09%. The average vitrinite reflectance (Ro) of the shale samples is 2.94%. The geochemical properties measured above indicate that samples in this work belong to the overmatured and organic-rich shales. The samples in the methane adsorption experiments are crushed to a diameter of 180-250 μm (60-80 mesh), and the samples in the effective stress sensitivity experiments are cores with three types of fractures. Figure 2 shows the photographs of shale samples used in the HPHT methane adsorption experiments and effective stress sensitivity experiments, respectively.
Permeability and porosity are two general physical parameters. Shale matrix is quite tight, but natural fractures in shale formation are generally developed. We utilized two methods to measure the permeability of core samples, including the method of CMC300 for samples' permeability higher than 0.01 mD and the method of pressure pulse decay for samples' permeability lower than 0.01 mD. Porosity of the shale samples is measured based on Boyle-Mariotte's law. Permeability and porosity of matrix core samples, which do not show any obvious natural fractures with naked human eyes, are measured under a confining pressure of 3 MPa and the results are shown in Figure 3. Permeability and porosity of the other core samples with obvious natural fractures under a confining pressure of 3 MPa are shown in Table 1. In addition, shale cores with three types of fractures are used in the effective stress sensitivity experiments, including natural-fracture cores, artificial-fracture cores without proppant, and artificial-fracture cores with proppant, which represent three types of fractures developed during shale gas production. Permeability of these cores was measured under in situ effective stress conditions, as shown in Figure 4.
The porosity shown in Figures 3, 4, and Table 1 might not thoroughly consider the nanopores in shale matrix. Therefore, two other methods to characterize samples' porosity are applied in this work, which are low-pressure gas (N 2 ) F I G U R E 2 Shale samples used in the experiments 60-80 mesh adsorption and field emission scanning electron microscopy (FESEM). N 2 adsorption measurement can determine a pore size from 0.35 to 300 nm, and its measuring pore size of the highest accuracy is usually less than 30nm. According to the pore volume measured by the low-pressure N 2 adsorption, shale porosity for diameter from 0.35 to 120 nm is determined, which is shown in Figure 5. It is indicated that 45% of the 29 matrix samples' porosity are between 4% and 5%. Meanwhile, the area porosity of single organic matter in a FESEM photograph is calculated by a binarization processing technique using a software named ImageJ, which is shown in Figure 6. In this work, the area porosity of single organic matter from 11 photographs of FESEM images is calculated and the results are shown in Table 2, indicating that the average area porosity of single organic matter is 8.98%. It is noted that the corresponding pore diameter in the above method is higher than 10nm based on the resolution ratio of FESEM photographs.

| Methane adsorption experimental method
A volumetric method was used in this work to get the methane adsorption isotherms of shale samples at HPHT conditions. The measurement principle of the volumetric method is shown in Figure 7. It is noted that the instrument is in a water temperature bath to maintain the methane adsorption experiments are conducted in a constant temperature condition. The essential technical indexes of the instrument are presented as follows: First, samples were crushed to a diameter of 180-250 μm (60-80 mesh). Second, the crushed samples were put in the sampling cell and the measurement system was put in the water temperature bath at 82℃ that was equal to the in situ reservoir temperature. Third, the measurement system was vacuumized. The sampling cell was kept in a vacuum by closing Valve 6 and Valve 7, and the helium was injected into the referencing cell. Then, the referencing cell was connected to the sampling cell after the helium pressure in referencing cell was stabilizing. The free space volume of the measurement system was determined according to the pressure measured above. Finally, the measurement system was vacuumized again. The sampling cell was kept in a vacuum state until the methane was injected into the referencing cell. Then, the volume of the adsorption phase was determined by measuring the change in methane pressure after the referencing cell was connected to the sampling cell. A methane adsorption isotherm of a shale sample was got through the above method under different methane pressure conditions.

| Effective stress sensitivity experimental method
Effective stress sensitivity experiments were conducted through a specially designed equipment to measure the permeability of cores at different effective stress conditions. A sketch of the effective stress sensitivity experimental setup is shown in Figure 8. In order to model the change in effective stress during production, the confining stress of a core was stabilized as the value of in situ overburden formation pressure. The pore pressure of a core was decreased from the initial formation pressure to the abandonment pressure of the shale gas well.
The experimental procedures are shown as follows: 1. Insert a shale sample into the core holder after checking the airtightness of the equipment. 2. Increase the confining stress and pore pressure simultaneously, and maintain the difference between the above two values equal to 2 MPa. Increase the pore pressure to 38 MPa which is equal to the value of initial formation pressure and then increase the confining stress to 59 MPa which is the value of in situ overburden formation pressure. 3. Maintain the confining stress unchanged. Decrease the backpressure to make a pressure difference between the two ends of core. It should be noted that the pore pressure is still stabilized to be 38 MPa. Then, the core permeability is measured. 4. According to the method shown in the above step, the core permeability at the pore pressure of 36 MPa, 34MPa, 32 MPa, 30 MPa, 27 MPa, 24 MP, 20 MPa, 15 MPa, 10 MPa, and 5 MPa is measured sequentially. 5. Gradually reduce the confining stress and pore pressure until reaching the atmospheric pressure. Then, remove the shale sample from the core holder.

| Methane adsorption isotherms on shales
Five crushed samples were selected in the methane adsorption experiments. One adsorption isotherm of a shale sample is presented in Figure 9. It is shown that the methane adsorption quantity increases as the increase in bulk methane pressure and the adsorption quantity tends to be constant as the pressure increases to a relatively high value. According to the Langmuir theory, the Langmuir volume (V L ) and Langmuir pressure (p L ) values of the five samples are calculated as shown in Table 3. It is indicated that the average values of V L and p L are 2.72 cm 3 /g and 7.57 MPa, respectively.

| Effective stress sensitivity of fracture cores
Permeability of cores with three types of fractures is measured at the different effective stress conditions. Results of three cores representing the three types of fractures are shown in Figure 10. The vertical ordinates in Figure 10 are dimensionless which are the ratio of permeability at a certain effective stress (k i ) to permeability at initial in situ effective stress (k 0 ). It is indicated that the core permeability decreases as the effective stress increases. However, the permeability change behavior is different for the three types of cores. For artificial-fracture cores without proppant, the permeability declines the most. The decline in permeability of artificialfracture cores with proppant is the smallest. Therefore, the effective stress sensitivity of artificial-fracture cores without proppant is generally the strongest, followed by the naturalfracture cores, and the effective stress sensitivity of artificialfracture cores with proppant is the weakest.
Fractures are the dominate flow channels in the effective stress sensitivity experiments. The hydraulic width of fractures during the experiments can be calculated based on the following cubic law of fluid flow in a single fracture.
in which e is the hydraulic width of fractures, k f is the fracture permeability, and D c is the core diameter. Then, the hydraulic width of fractures in the effective stress sensitivity   (1). It is indicated that (1) the hydraulic width of natural fractures at initial effective stress condition is between 481 nm and 878 nm, the average value of which is about 500 nm; (2) for artificial fractures without proppant, its hydraulic width at initial effective stress condition is between 976 nm and 3949 nm and the average value is about 1000 nm; and (3) the hydraulic width of artificial fractures with proppant at initial effective stress is between 9863 nm and 20 421 nm, the average value of which is about 10 000 nm. In order to characterize the change in permeability influenced by effective stress sensitivity, a correction factor of the effective stress sensitivity β p is defined, which is a ratio of the permeability at a certain effective stress to the permeability at initial effective stress. In this work, flow channels are divided into four types including matrix nanopores, natural fractures, artificial fractures without proppant, and artificial fractures with proppant. Different scales of pores or fractures have different β p . For hydraulic radius of organic pores smaller than 100 nm, the corresponding study of Wang et al (2012) 34 is referenced to determine the expression of β p . For the inorganic pores or fractures larger than 100 nm, we choose three representative hydraulic radius values calculated above, including 500 nm, 1000 nm, and 10 000 nm, which correspond to the scale of main flow channels of natural-fracture shale, artificial-fracture shale without proppant, and artificial-fracture shale with proppant, respectively. According to the experimental results of effective stress sensitivity, the mathematical relationship between permeability and effective stress is regressed for shale samples with the three kinds of fractures. Then, β p is expressed as follows.
in which α f is the Biot coefficient for microfracture, α s is the Biot coefficient for matrix, K n is the microfracture normal stiffness, E s is Young's modulus of matrix, b 0 is the initial microfracture aperture, s is the microfracture spacing in matrix, and p 0 is the initial pressure.

| APPARENT PERMEABILIT Y MODEL
Effective stress of shale formation increases during production, which would change the sizes of flow channels including fractures as well as matrix nanopores. It is necessary to propose a unified model to describe apparent gas permeability considering the changes in flow channels sizes. The investigations of shale formation characteristics and gas transport mechanisms show that (1) there are three main types of gas flow channels in shales including organic pores, inorganic pores, and fractures; (2) gas flow mechanisms in shales contain desorption of adsorption gas, configurational diffusion, surface diffusion, Knudsen diffusion, and viscous flow; and (3) shale gas transport behavior during production is influenced by some microscale effects such as a real gas effect, a sorption layer effect, and an effective stress sensitivity. Therefore, gas transport in shales is a multiscale mass transport behavior which happens in a multiple pore/fracture sizes, is influenced by multiple microscale effects, and contains a variety of gas transport mechanisms.

Knudsen diffusion
According to researches of gas diffusion through a porous alumina disk, Roy et al (2003) 35 provided a mathematical expression to describe the relationship among Knudsen diffusion coefficient, pressure gradient, and molar flow rate, which is shown as: in which F is the molar flow rate, D k is the Knudsen diffusion coefficient, A is the cross-sectional area of core, Δp is the pressure drop across the core, R is the gas constant, T is the temperature, and L is the length of core.
Here, we introduce the concept of equivalent permeability and follow the form of Darcy's law to build the relationship between equivalent permeability of Knudsen diffusion and gas flow rate: in which k D is the equivalent permeability of Knudsen diffusion and q is the flow rate.
Knudsen diffusion coefficient in Equation (3) can be given as follows: Comparing the expression shown in Equation (3), (4), and (5), the equivalent permeability of Knudsen diffusion in a shale is expressed as: In order to facilitate the calculation of k D , a real gas state equation (Equation 7) is inserted and Equation (6) can be given as Equation (8).
in which Z is the gas compressibility factor. The Knudsen diffusion is influenced by the changes in flow channels during production, which include the change in adsorption layer thickness and the change in flow channels sizes. In this work, two factors, β p and β s , are used to correct the changes in permeability induced by effective stress sensitivity and adsorption layer effect, respectively. The mathematical expression of β p is shown in Equation 2. β s is a correction factor of the sorption layer effect, which is given as: in which p 0 is an initial pressure, p L is a Langmuir pressure, and ε L is a Langmuir strain.
Therefore, the equivalent permeability of Knudsen diffusion during shale gas production is presented as follows comparing Equation (2), (8), and (9).
in which r H is the hydraulic radius of flow channels.

| Equivalent permeability of viscous flow
Bulk gas transport in shales is a comprehensive result of Knudsen diffusion nearby pore wall and viscous flow away from the boundary. The equivalent permeability of viscous flow is presented as follows considering the sorption layer effect and the effective stress sensitivity, in which β p is expressed in Equation (2) and β s is expressed in Equation (9).

| Equivalent permeability of surface diffusion
An organic-rich shale has a very big surface area, which is quite conducive to gas sorption. Therefore, the total gas flux in shales during production is also influenced by surface diffusion besides the above two bulk gas transport mechanisms. In this work, a walk/hopping mechanism is utilized to describe surface diffusion of shale gas. According to the surface diffusion model of Hwang and Kammermeyer (1996) 36 and Chen and Yang (1991), 37 surface diffusion of shale gas is given as: in which D s is surface diffusion coefficient, ∆H is isosteric adsorption heat at the gas coverage of zero that can be obtained by the methane adsorption experiments conducted in this work, θ is the gas coverage on the pore wall, and H(1-κ) is the Heaviside function.
According to the theory of Langmuir, θ is given as: For the Heaviside function in Equation (12), κ is ratio of the blocking-velocity coefficient to the forward-velocity coefficient.
The parameter κ is given as 0.5 based on the study of Wu et al (2016) 32 on shale gas surface diffusion.
According to the theory of Maxwell-Stefan, surface diffusion is driven by a chemical-potential gradient. 38 Meanwhile, since the velocity of surface diffusion is far smaller than that of gas adsorption/desorption, the influence of gas adsorption/desorption on the surface diffusion can be neglected. Therefore, the mass flux of surface diffusion is given as: Here, the concept of equivalent permeability of surface diffusion is still introduced as: in which M is the methane molar mass, and C s is the concentration of adsorption gas that can be given as follows according to the theory of Langmuir.
in which δ is the collision diameter of methane molecule and N A is the Avogadro's constant.
In Equation (16), ξ s is shape correction factor of surface diffusion, which is shown as follows based on the derivation of Wu et al (2016). 32

| Apparent methane permeability
The expressions of equivalent permeability of each gas transport mechanism in shales during production are shown in Equation (10), (11), and (16), respectively. Then, an expression of apparent gas permeability can be given as follows, which is a weighted sum according to the theories of molecule collisions which will be explained in detail later.
in which α and (1-α) are the proportion of Knudsen diffusion and viscous flow in bulk gas flux, respectively. ω is the reduction rate of equivalent permeability of viscous flow caused by the collisions between bulk methane molecules and adsorption methane molecules. β p and β s are expressed in Equation (2) and Equation (9), respectively.
From the perspective of molecule collisions, α in Equation (19) can represent the proportion of collisions between gas molecules and solid molecules of pore wall in the total molecule collisions. And (1-α) represents the proportion of the collisions among gas molecules. According to the studies of Reif (1965) 39 and Thompson and Owens (1975) 40 on gas-gas and gas-solid molecule collisions, the relationship between α and Knudsen number (Kn) is shown as follows.
In Equation (19), ω is the reduction rate of equivalent permeability of viscous flow caused by the collisions between bulk methane molecules and adsorption methane molecules.
As shown in Figure 11, there are two kinds of gas molecule collisions in shales, including collisions among bulk gas molecules and collisions between bulk gas molecules and adsorption gas molecules. The essential of Equation (11) is the collisions among bulk gas molecules. However, Equation (11) cannot reflect the influence of collisions between bulk gas molecules and adsorption gas molecules. From the perspective of equivalent permeability, Equation (11) would overestimate the contribution of collisions between bulk gas molecules and adsorption gas molecules to the equivalent permeability. Therefore, we introduce the parameter ω to make a correction.
For the parameter ω in Equation (19), it is related to properties of pore wall, shale heterogeneity, reservoir temperature, and pressure. Kn is just related to the above factors. Meanwhile, according to the studies of Thompson and Owens (1975) 40 on gas-gas and gas-solid molecule collisions, there should be a relationship between ω and Kn. The experimental data of several steady-state methane flow experiments in shales are used in this work to fit ω and Kn. 26 The results confirm that the relationship between ω and Kn has a fair degree of consistency. Expression of ω is given as follows.
The mathematical model of apparent permeability shown in Equation (19) comprehensively considers the following factors, including a real gas effect, an effective stress sensitivity, a sorption layer effect, porosity, tortuosity, surface roughness, and the weighted contribution of each gas transport mechanism. Moreover, the dynamic changes in effective flow channels sizes, including the change in adsorption layer thickness and the change in pore/ fracture sizes, are thoroughly considered.
In order to validate the model in shale gas flow, the steady-state methane flow experimental data from three shale core samples are used in this work. 26 Comparisons between the mathematical model and experimental data are presented in Figure 12 and show good matches. Therefore, it is demonstrated that the unified model can be used to describe the methane flow experimental results in shale cores. In addition, two other classical and widely used models, including Klinkenberg model 16 and Beskok-Karniadakis model, 29 are compared with the model shown in Equation (19) for validation. First, a correction factor f c is defined, which is the ratio of apparent gas permeability to equivalent liquid permeability. So, three expressions of f c that correspond to the unified model in this work, Klinkenberg model, and Beskok-Karniadakis model are shown as follows, respectively.
The relationship between pressure and f c for different pore size under reservoir temperature (82℃) is shown in Figure  13 based on Equation (22), (23), and (24). It is shown that f c of the three models is smaller and gradually approach 1 with the increase in formation pressure or pore size, indicating an important role of viscous flow in a large pore size or high-pressure condition. However, the difference between the model in this work and two other models increases as the formation pressure or pore size decreases, indicating a non-negligible effect of Knudsen diffusion and surface diffusion in a large Kn condition. Although the Beskok-Karniadakis model introduces the rarefaction coefficient compared with the Klinkenberg model, which improves the calculation accuracy to a certain extent, it still only corrects the gas slip condition of pore wall boundary based on the hydrodynamic continuity model. Meanwhile, it poorly considers the real weighted contribution of each gas transport mechanism. Moreover, it does not consider the collision of gas molecules ( Figure 11) under a small pore size or low-pressure condition. In particular, the Beskok-Karniadakis model and the Klinkenberg model neglect the surface diffusion on the pore wall of organic-rich porous media, and the surface diffusion actually plays an increasingly important role in a smaller pore. Therefore, since the unified model in the work comprehensively considers all the influencing factors besides the dynamic changes in effective flow channels, the model proposed in this work is more robust for shale gas flow.

| RESULTS AND DISCUSSION
According to the measured petrophysical properties, geochemical properties, methane adsorption properties, and pore structure as well as the characteristics of shale gas reservoir studied in this work, Table 4 presents the values of parameters in the unified permeability model shown in Equation (19). In Table 4, the value of p 0 equals to the initial in situ pressure of the shale gas reservoir. The value of p final is determined based on a predicted abandonment pressure of a shale gas well. The value of T equals to the in situ temperature of the shale gas reservoir. The seven values of r H are determined based on the pore structure measured in this work as well as the experimental results of effective stress sensitivity. The isosteric adsorption heat ∆H is closely related to the properties of methane-shale system. The value of ∆H in this work is determined as 18.64 kJ/mol based on a corresponding calculation result for the shale gas reservoir studied in this work. 13 The apparent methane permeability for multiple pore/ fracture sizes during production is calculated as shown in Figure 14. It is indicated that the apparent permeability of each representative size of flow channel changes obviously with the depletion of reservoir pressure.
(1) For r H lower than 10 nm, the permeability decreases at first and then increases with the decline of reservoir pressure. Also, it should be noted that when the r H decreases to a very small value, such as 2 nm shown in Figure 14, the permeability improves directly as the reservoir pressure decreases. Since the major peak of the matrix samples is 3-4 nm, 41 the modeling results show that the transport capability of shale matrix nanopores has an obvious tendency of increase in the mid-late production period, which is directly related to contribution of each gas transport mechanism during the decline of pore pressure and increase in formation effective stress that would be discussed in detail in the rest of section 5. Meanwhile, the modeling results are in line with the sensitivity of pore size and pressure for methane gas transport studied by the molecular dynamics simulation of shale matrix pores with a few nanometers. 42 (2) When r H is as high as 50 nm and above, the permeability decreases as the reservoir pressure declines. The results are generally in line with the shale gas well production performance, especially in the early production period, because the stress sensitivity of fractures is usually very obvious. Also, the results show that the decline of permeability for inorganic pores/fractures tends to be more and more gentle as the decrease in reservoir pressure, especially for the smaller flow channels. The modeling results are actually controlled by the response of each gas transport mechanisms to reservoir pressure changes, especially the effective stress changes, which would be explored in detail by the analysis shown in the rest of section 5.
The shale gas transport behavior for different scale of flow channels analyzed above is essentially controlled by the different gas transport mechanisms during production. According to the modeling results of Equation (19), the transport capability of each gas transport mechanism as well as their corresponding contribution to the total gas flux is investigated as follows.
(1) The equivalent permeability of Knudsen diffusion and its contribution to the total gas flux during production are shown in Figure 15. The impacts of Knudsen diffusion during shale gas production include the dynamic changes in flow channels and the collision degree between the gas molecules and the solid molecules on the pore wall. And the dynamic changes in flow channels are determined by the stress sensitivity and the adsorption layer effect. Figure  15 shows that when r H is not larger than 100 nm, the equivalent permeability of Knudsen diffusion decreases slightly as the decline of reservoir pressure. When r H is larger than 100 nm, the equivalent permeability of Knudsen diffusion shows a slight decrease and then increase with the decline of reservoir pressure. According to the impacts of Knudsen diffusion during shale gas production analyzed above, for gas flow channels larger than 100 nm, the mean free path of gas molecules increases with the decrease in reservoir pressure. In addition, the size of flow channel is narrowed induced by the effective stress sensitivity. Therefore, the equivalent permeability of Knudsen diffusion tends to increase for gas flow channels larger than 100 nm during production. For gas flow channels smaller than 100 nm, the decrease in effective size of flow channels due to the decline of reservoir pressure during production might be beyond the tendency of permeability increase due to the increase in the mean free path of gas molecules. Therefore, the equivalent permeability of Knudsen diffusion decreases slightly with the decrease in reservoir pressure. The contribution of Knudsen diffusion to the total gas flux during shale gas production is also presented in Figure 15. It is indicated that the contribution of Knudsen diffusion in the large scales of pores/fractures increases with the decline of reservoir pressure, such as r H of 500 nm, 1000 nm, and 10 000 nm. However, the increase is small, which does not exceed 1% in general. For r H of 2 nm, the contribution of Knudsen diffusion is 5.62% at initial reservoir pressure, and then it gradually decreases to 0.43% with the decline of reservoir pressure. It should be noted that the contribution of Knudsen diffusion changes the most obviously for r H of 10 nm, 50 nm, and 100 nm.
(2) The equivalent permeability of viscous flow and its contribution to the total gas flux during production are shown in Figure 16. It is presented that the equivalent permeability of viscous flow decreases with the decline of reservoir pressure, which is due to the decrease in collision among gas molecules and the dynamic changes in flow channels sizes during production. The decrease in reservoir pressure would make gas desorb. It would weaken the adsorption layer effect and strengthen the effective stress sensitivity, the former of which is conducive and the latter of which is not conducive to the gas transport capacity of viscous flow. According to the modeling results shown in Figure 16, it is indicated that the influence of adsorption layer effect can be ignored. In addition, the shale pore structure influences the gas transport capacity of viscous flow. Although the value of the equivalent permeability of viscous flow is small for nanoscale flow channels, it decreases with the decline of reservoir pressure relatively larger. On the one hand, the smaller the sizes of nanoscale pores are, the smaller the adsorption layer effect weakens as the reservoir pressure declines, and on the other hand, the nanoscale pores are generally weak in strength and easy to deform. Therefore, the effective sizes of nanoscale flow channels in shale matrix have the greater decrease in equivalent permeability of viscous flow during production. For the flow channels sizes with hundreds or thousands of nanometers, they generally belong to inorganic pores, natural fractures, artificial fractures without proppant, and artificial fractures with proppant. Since the above flow channels are relatively large, the collision among gas molecules is strong and the adsorption layer effect is relatively weak. Meanwhile, the ability of flow channels to resist the change in pore/fracture sizes induced by effective stress sensitivity generally increases with the increase in pore/ fracture sizes. Therefore, the larger the pore/fracture sizes are, the smaller the equivalent permeability of viscous flow decreases during production. The contribution of viscous flow to the total gas flux during shale gas production is also shown in Figure 16. It is indicated that the contribution of viscous flow is different in different scales of flow channels. The most obvious changes in contribution F I G U R E 1 4 Apparent permeability of seven representative flow channels during production F I G U R E 1 5 Equivalent permeability of Knudsen diffusion and the contribution to the total gas flux for seven representative flow channels during production of viscous flow are in the flow channels with r H of 10 nm, 50 nm, and 100 nm. The contribution of viscous flow for r H of 2 nm is inherently small, which is not larger than 1% during production. The contribution of viscous flow for r H of 500 nm, 1000 nm, and 10 000 nm is inherently large, which is larger than 95%.
(3) The equivalent permeability of surface diffusion and its contribution to the total gas flux during production are shown in Figure 17. It is indicated that the equivalent permeability of surface diffusion gradually increases with the decline of reservoir pressure. Although the surface diffusion coefficient and the concentration of adsorbed gas decrease due to the decrease in bulk gas pressure that would reduce the gas coverage on the pore wall, the decrease range of adsorption gas concentration is not as large as the pressure decrease in Equation (16). Therefore, the equivalent permeability of surface diffusion increases with the decline of reservoir pressure. In addition, the equivalent permeability of surface diffusion increases obviously with the decrease in pore/fracture sizes. On the one hand, the smaller the flow channels are, the larger the surface area is, which provides a broader space for surface diffusion, and on the other hand, the coverage of adsorbed gas in nanoscale flow channels is much greater compared with the inorganic pores or fractures. Therefore, the equivalent permeability of surface diffusion is modeled to be larger in the smaller pores. The contribution of surface diffusion to the total gas flux during shale gas production is also presented in Figure 17. It is shown that the smaller the flow channels are, the greater the contribution of surface diffusion is. The contribution of surface diffusion is more than 90% in pores with r H of 2 nm during production, and it is less than 1% in pores/fractures larger than 100 nm. In addition, it should be noted that the contribution of surface diffusion changes the most in pores with r H between 2 nm and 10 nm during production. In general, the contribution of surface diffusion in matrix pores below 100 nm cannot be neglected.
According to the above analyses, gas transport behavior in shale formation during production is different in different scales of flow channels. Overall, the contribution of surface diffusion, Knudsen diffusion, and viscous flow in the early production period accounts for a considerable proportion. Especially in the pores with r H of about 2 nm, the surface diffusion absolutely dominates the total gas flux in shales. The role of surface diffusion gradually weakens with the increase in pore/fracture sizes. When r H increases to no more than 100 nm, the contribution of Knudsen diffusion cannot be neglected, but it is generally very small in pores/fractures with r H less than 2 nm and greater than 100 nm. For the contribution of viscous flow, it increases significantly with the increase in pore/fracture sizes, and it is high enough in flow channels with r H larger than 10 nm. In the mid-late production period, surface diffusion still dominates the total gas transport behavior in matrix pores with r H of about 2 nm, but Knudsen diffusion shows a more and more important role in pores/fractures with r H of hundreds of nanometers. For viscous flow, although its transport capability decreases as the reservoir pressure declines, it still dominates in pores/fractures with r H larger than 100 nm. F I G U R E 1 6 Equivalent permeability of viscous flow and the contribution to the total gas flux for seven representative flow channels during production F I G U R E 1 7 Equivalent permeability of surface diffusion and the contribution to the total gas flux for seven representative flow channels during production 6 | CONCLUSIONS This work studied gas transport behavior in shales with multiple pore/fracture sizes at in situ conditions. Methane adsorption experiments and effective stress sensitivity experiments were conducted to determine the changes in effective flow channels sizes during production. Then, a unified model of shale gas apparent permeability was proposed based on the laboratory measurement results and investigation of gas transport mechanisms. Finally, methane transport behavior in shales with multiple pore/fracture sizes was discussed in detail. Conclusions from this work are summarized as follows.
1. Methane adsorption isotherms and permeability of different types of fracture cores were measured to determine the methane adsorption properties and effective stress sensitivity, respectively. The dynamic changes in effective flow channels sizes during production are mathematically characterized according to the adsorption layer effect and effective stress sensitivity. 2. A unified model to calculate the apparent gas permeability of a shale is proposed based on a weighted contribution of each mechanism to the total gas flux. The model comprehensively considers multiple gas transport mechanisms and multiple microscale effects including adsorption layer effect and effective stress sensitivity which make the effective size of flow channels change during production. 3. Apparent permeability of a shale is modeled by the unified model over multiple sizes of flow channels. The change behavior of apparent permeability during shale gas production is clarified. 4. According to the analyses of shale gas transport mechanisms, equivalent permeability of Knudsen diffusion, viscous flow, and surface diffusion are modeled and their contribution to the total gas flux is studied, respectively. The transport behavior of each gas transport mechanism during different production periods is clarified. 2016ZX05061), the Sichuan Province Youth Science and Technology Innovation Team Special Program (No. 2016TD0016), and the National Natural Science Foundation of China (No. 41772150). Langmuir strain, dimensionless θ gas coverage on pore wall, dimensionless μ dynamic viscosity of gas, Pa s ξs shape correction factor of surface diffusion, dimensionless ρ density of gas, kg/m3 τ tortuosity, dimensionless ϕ porosity, % ω reduction rate of kv as collisions between bulk and adsorption gas molecule, dimensionless