Abstract
Because of the existence of multiscale pores from nano- to macroscale, a multimechanistic shale gas flow process involving the Darcy and Knudsen flows occurs during gas shale well depletion. The respective contribution of the Darcy and Knudsen flows to the permeability is constantly changing with pressure evolution. In this study, laboratory measurements of shale permeability with CO2 injections were carried out under hydrostatic conditions, using the transient pulse-decay method. The “U”-shape permeability curve resulted in both positive and negative effective stress coefficients (Biot’s coefficient) . A permeability turning point was thus created to partition permeability curves into the Darcy and Knudsen sections. The Knudsen effect was proven to be significant at low pressure/late time in the laboratory. Effective stress and sorption-induced deformation have been found to govern the Darcy permeability evolution under the tested experimental conditions. Thus, negative effective stress coefficients, together with the positive ones, should be applied to a nonmonotonic pressure-permeability evolution to explain the concurrent effect of the Darcy flow and Knudsen flow at different pore pressures.
1. Introduction
Gas shale reservoirs have played an important role in natural gas world supply. Hydraulic fracturing stimulation is an effective method to create fracture networks that can connect shale formations to horizontal boreholes [1]. However, low matrix permeability retards gas migration from the tight matrix into the fracture networks. Pore diameters of the shale matrix can be smaller than 2 nm, which is significantly smaller than average pore size in conventional sandstone and carbonate reservoirs [2–5]. The shale matrix, which includes organic matter and clay minerals, serves as the storehouse for gas but limits gas transport due to its low permeability and low diffusivity. The tight structure can involve a variety of flow dynamics processes, making the evaluation of gas migration behaviors intricate and challenging.
Gas transport in shales contains Darcy flow and Knudsen flow [6, 7] and the shrinkage/swelling response of organic components of shale due to gas desorption/adsorption [8–10]. The extremely tight structure of a shale core sample without apparent fractures can limit the gas transmission flux under high effective stress conditions during shale permeability laboratory tests. Since permeability testing time scales with the square of diffusion length, shorter shale core samples were used for permeability measurement which practically made laboratory testing feasible and time effective. In addition, permeability values obtained using classic Darcy’s law may not be valid if non-Darcian flow is taking place [11]. Nanopore and micropore gas flow along with sorption effect needs to be separately characterized and integrated to account for the multimechanistic shale gas flow process [12]. In the above laboratory and modeling researches, they have massively discussed gas flow mechanisms in shale with different stress conditions under various pore scales. However, there is no certain conclusion either on how to reflect those flow types in a gas permeability test or on giving a comprehensive physical understanding to the effective stress in a shale gas flow process. In an actual shale gas flow process, due to the extremely small pore size and/or volume in shale formations and the complex matrix shrinkage/swelling effect, conventional Darcy flow and diffusion models cannot, on their own, describe gas flow transport in shale formations [6]. A multimechanical method includes the analysis and breakdown of the Darcy viscous flow, and the Knudsen flow is thus needed.
In this study, two cut and polished shale samples were used to conduct permeability tests using the pulse-decay method. The Darcy flow and Knudsen flow were separately analyzed to determine their respective contributions to total gas transport. Other important effects—such as effective stress and sorption-induced deformation, which have been frequently reported in unconventional rock studies—were considered to influence gas transport properties. Gas flow mechanisms and effective stress law were fully considered, and their physical meaning was discussed in a laboratory scale through data analysis.
2. Background
2.1. Shale Gas Flow Regimes
During shale gas reservoir production, the Darcy flow permeability is still a prevailing term used to quantify gas deliverability through shale fractures and matrices. However, the Knudsen diffusion is just as important as Darcy’s convective flow in unconventional gas production prediction. Figure 1 shows multiple flow mechanisms during gas shale reservoir depletion according to some existing flow partitioning criteria [13, 14]. The average monthly production for 15 wells in the Eagleford Shale formation shows that peak production occurred late in the second month after the fracturing treatment. This was followed by a rapid decline until the 20th month; then, there was a production plateau at a relatively low production rate. Early in production, fracture (Darcy) flow dominated since free gas molecules were flowing through hydraulic fractures (large-diameter flow channels) and occupying the flow path, driven by pressure gradient. Production reached its peak very quickly due to fracture flow. Following peak production, rapid decline occurred due to quick depletion of compressive fracture gas storage.

However, desorption could be involved during the pressure decline stage; only a negligible amount of matrix gas contributed to production because of the slow diffusion process [16–18]. Generally, the Knudsen diffusion flow at this stage does not dominate total gas production compared to the free gas flowing through the fractures. Gas transport mechanism moves from continuum flow to a transition zone between the Darcy flow and Knudsen diffusion. This transition zone includes slip flow () and transition flow (), where neither continuum flow nor free molecular flow is valid [16]. Gas flow in this zone can be considered a combination of the Darcy flow and Knudsen diffusion. After the fractures are depleted, a late-time asymptotic flow takes place, when reservoir fracture pressure becomes very low. At this stage, although fracture permeability is sufficient, gas influx from the matrix towards the fractures will take over the role of fracture permeability to control the well production. This mass influx from the matrix towards the fracture is a multimechanistic micro/nanoscale transport. Methane molecules first desorb from the internal surfaces of the matrix within the nanopore system due to fluid pressure driving force and become free gas [19, 20]. During the desorption process, the distance between the shale gas surface molecule layer and the adjacent molecules decreases, and the effective surface energy on the shale kerogen surface increases [21]. As a result, kerogen shrinkage also occurs in the shale matrix and enhances the nanopore flow path for gas molecules (Figure 2). Because of the wide range of pore sizes within the shale matrix, a combination of multiple flow regimes can be expected. Most gas molecules migrate through the nanopores by means of the Knudsen diffusion flow. The diffusion-dominated flow can remain operational for decades, contributing to the late production flat tail. Gas flow permeability evolution at the late-time production stage is critical and significantly influences the production behaviors [18]. In short, shale gas transport during production includes different flow regimes (Darcy, slip flow, transition flow, and Knudsen diffusion with gas desorption) at every stage, but their respective contributions to the total production change with decreasing pore pressure as the depletion continues.

The breakdown of different flow regimes, such as the Darcy flow, slip flow, and transition flow, can be classified and distinguished by the Knudsen number (), which is a function of pore size and pore pressure [22, 23]: where is the gas mean free path, defined in Nomenclature and Greek Symbols:
The relationship between , pore size, and pressure for CO2 is also plotted in Figure 3, over the experimental pressure range used for the subsequent laboratory study. The temperature was at 296 K (room temperature). For Marcellus Shale reservoirs, the pore diameter is expressed in nanometers and micrometers [24]. Within the pore size range given in Figure 3, gas flow in the shale matrix lies mostly within the transition flow and slip flow region. It can potentially involve the Darcy flow regime if the pressure is sufficiently large or include the Knudsen diffusion when this pressure is low. In order to assess these flow mechanisms, the permeability measurements in this work were all conducted at relatively low pore pressures (<10 MPa) to ensure that 2~3 flow regimes can be covered. Since transition flow and slip flow can be modeled as a weighted combination of the pure Poiseuille flow (Darcy) and Knudsen flow [13, 25], coupled flow regimes were applied to describe gas flow and permeability evolution for gas flow in the shale matrix.

2.2. Permeability Measurement
Brace et al. first introduced the pulse-decay technique as a transient method derived from Darcy’s law and used to measure permeability by applying a pressure difference between two sides of a core sample [26]. Hsieh et al. derived more restrictive analytical solutions of the differential equation, describing decay curves from permeability measurement with compressive storage effect [27]. The exact solution of the differential equation for dimensionless pressure difference and dimensionless time was improved and shown as [28] where and are the ratio of the sample’s storage capacity to that of the upstream reservoir and downstream reservoir and is the th root of the following equation: where and .
To simplify the above method, Jones introduced a factor as follows [29]:
The pulse-decay equation can be described as
Laboratory estimation of shale permeability for unconventional reservoir rocks conducted under hydrostatic conditions with an adsorption effect was reported by Cui et al. An effective adsorption porosity term was introduced to account for the contribution of gas molecule adsorption. The Langmuir model was used to quantify gas adsorption volume as a function of pressure [30] and mathematically described as follows:
So, the sample storage capacity ratio in Cui et al.’s approach becomes where . Recently, Wang et al. pointed out that Cui et al.’s method highlighted the contribution of the sorption effect in permeability calculation, demonstrating how the gas sorption effect can help to predict gas permeability [31]. Since our study involved handling organic-rich shale samples, the pulse-decay permeability calculation method proposed by Cui et al. was chosen to estimate apparent permeability for sorbing gas.
2.3. Effective Stress
To quantify the influence of effective stress on gas permeability, confining and pore pressures should be carefully related to rock and gas properties in pulse-decay tests. Terzaghi effective stress (or simple effective stress), as a function of confining and pore pressure according to Terzaghi’s principle, is defined as [32]
However, for permeability tests, the effective stress law with an appropriate coefficient has been used instead of Terzaghi’s principle [14, 33, 34]:
To experimentally calculate the effective stress coefficient , the ratio of the slope method was presented by Kwon et al.:
The effective stress coefficient represents the sensitivity of permeability to the changes in confining and pore pressures.
3. Experimental Methods
3.1. Sample Procurement and Preparation
With tight structure, equilibrium time for pressure pulse-decay is extremely long for long core samples. This long equilibrium time reduces the efficiency of lab testing, and no pressure drop was observed under very high stress conditions after a few hours of pulse injection. Therefore, shale disks were preferable for permeability measurements; they can significantly shorten equilibrium time and are commonly used for shale permeability measurements [14]. Thus, in this study, Marcellus Shale drilled cores were prepared as thin disks for the gas apparent permeability tests. The prepared thin disks ranged in thickness from 3 to 6 mm and were 25.4 mm in diameter. All trimmed surfaces were polished to enable proper placement in the triaxial cell. Two well-prepared shale samples of different thickness are used in this work, and one of them is shown in Figure 4. Each disk sample has been weighed and placed in an oven for the first 24 hours at 150°C to remove the moisture. The second weighing was conducted before they were put into the oven again. We repeated this procedure until the sample weights remain stable. After the shale sample disks were dried in the oven, they were weighed in a dry and clean plastic sample bag in a lab-use alloy box for 3 hours before testing.

3.2. Experimental Stress Boundary Conditions
Hydrostatic conditions were applied in this series of shale permeability measurements, in which axial stress was equal to confining pressure at a constant value throughout the duration of each experiment. Hydraulic stresses were regulated with computer-controlled syringe pumps using software developed in the LabView environment.
3.3. Experimental Setup and Procedure
The measurements presented in this study are recovered from a standard triaxial apparatus arranged for flow-through or pulse permeability testing as shown in Figure 5. The apparent permeabilities of two Marcellus Shale samples of different thickness were measured under hydrostatic pressure. Using carbon dioxide as the test gas, the permeability of sorbing gas (CO2) was compared by using a high-pressure pump for confining pressure loading and gas pressurization. The experiment process is as follows: (1) apply a certain confining pressure to the sample chamber, and vacuum the fluid pipeline and sample at the same time. (2) In the upstream of the system to the part of valve 1, inject a certain pressure of the test gas, the sample, and the downstream part of the original state. (3) After the injection pressure is stabilized, open valve 1, drive the gas pulse through the shale sample by the pressure difference between the upstream and downstream, read the pressure data after the overall pressure balance of the system, and close valve 1. (4) Repeat steps (1) to (3) until all planned pressure levels are tested.

The balancing time lasts for 30 minutes to ensure that the balancing process is complete. In the experiment, four groups of complete cycle permeability tests were carried out on samples of 5.39 mm and 3.71 mm with CO2 under hydrostatic pressure of 11 and 21 MPa, respectively. Each group of tests was repeated once, and a total of 16 groups of tests were carried out. Finally, based on multiple groups of pressure data obtained in the experiment, equations (6) and (7) were used to calculate sample permeability. An entire CO2 injection cycle on the 5.39 mm sample at 11 MPa is displayed in Figure 6.

4. Results
The measured permeability results of CO2 for two disk samples (two cycles each) are shown in Figure 7. According to Figure 7, the repeatability of the permeability measurements is good. Also, the permeability data for the two samples are extremely close to each other, with very similar trends. In our experiments, the structures of the shale sample matrix were not damaged due to the application of a high confining and/or pore pressure. Our measured permeability values are relatively high among some reported shale permeabilities [35], due to the microfractures present in the shale disk samples. These permeability measurements are still reasonable because Gensterblum reported a collection of permeability data for an organic shale sample case, with permeability ranging from to with a porosity of 2-8%. And the measured porosity of our Marcellus drilled samples ranges from 2% to 5.1%, obtained using the Small-Angle Neutron Scattering (SANS) technique, with a permeability range of 1.5 to . Based on the two sample observations in this study, the consistency of data between samples different in length indicates that the matrix structures for the sample collected from the same source can be very similar. Consistent and accurate results help us justify the utility of collected experimental data and consolidate any hypotheses and conclusions based on them.

4.1. Permeability Evolution with Pore Pressure
For both samples, shale permeabilities showed similar evolution patterns. All the permeability values initially declined when pore pressure ranges from 0.35 to 4 MPa, then started to make a slight or sharp turn, and increased with the decrease in the Terzaghi stress, since the Terzaghi stress is negatively and linearly associated with pore pressure. At low pore pressure, nano-/microscale transport, such as the Knudsen diffusion, dominates and has a decreasing contribution with increasing average pore pressure. With decreasing Terzaghi stress, shale gas permeability still has a decreasing tendency in the pore pressure range 4 to 14 MPa. This is because the nano- and micropores and the low gas pressure limit the contribution of the Darcy flow.
4.2. Permeability Evolution with Effective Stress
In order to further study and address the impact of effective stress on shale permeability, confining pressures () and pore pressures () were linearly plotted with the logarithm of apparent permeability () to determine the effective stress coefficient () and the effective stress (). The result of permeability versus confining pressures at constant pore pressures for the 3.71 mm shale sample is shown in Figure 8(a). Almost the same slopes were obtained from two pore pressure points, which led to equal to -0.0112 when pore pressure is 0.4 MPa and -0.0103 when pore pressure is 6.5 MPa. The negative values indicate that apparent permeabilities decrease with the increase in confining pressure [14, 34].

(a)

(b)
Then, we plotted values when is equal to 11 MPa and 21 MPa for the 3.71 mm sample, respectively, breaking down the pressure at ~4 MPa into two divisions. In Figure 8(b), we can clearly observe that the slopes for versus are always positive at high pore pressure division (>4 MPa) and negative at low pressure division (<4 MPa). These unexpected slopes resulted from the change of the controlling flow regime during pressure depletion. The Knudsen diffusion and slip are more important because the ratio of the mean free path to the pore throat is larger. This is due to lower pressure and narrowing of pore throats (flow paths).
On the other hand, Kwon et al. and Heller et al. obtained only one positive slope when plotting versus . They may have assumed as one property for identical rock samples and only vary with the distribution of quartz and clay minerals. This means that the number of effective stress coefficients obtained depends on the monotony of the permeability data profile, and it may not be valid to apply a unique value to evaluate the relationship between permeability and effective stress.
Using equation (12), we calculated the effective stress coefficients for each individual. Corresponding to the values, became negative at low pressure (<4 MPa) and positive at high pressure (>4 MPa). According to equation (12), we will obtain a positive when is positive and vice versa. Within both low and high pressure sections, permeability decreases with the raise of effective stress based on effective stress law, which is consistent with previous results presented by Kwon et al. and Heller et al. In this study, we applied both positive and negative effective stress coefficients based on “U”-shape permeability evolution trends at varying pressure levels. These two different approaches will allow us to account for the non-Darcian flow components using the traditional effective stress concept.
5. Discussion
5.1. Impact of Flow Mechanisms on Permeability
Our permeability measurements were conducted at a relatively low pressure (<10 MPa), in order to address the diffusion and transition flow effect. Permeability values can reach an extremely high peak due to opening of the transport pore space resulting from confining pressure reduction. So, at the early stage of gas shale well production, low confining pressure (due to high initial reservoir pressure) and highly opened fractures significantly promote the Darcy flow suggested by high apparent permeability. The Darcy flow plays a more important role at this circumstance while the Knudsen flow always happens simultaneously, according to Figure 1. At the late stage, low pore pressure combined with narrow flow paths not only increases the Knudsen number but also increases the contribution of Knudsen flow components [36]. This is one of the main reasons for researchers to integrate multiple flow influxes into an apparent permeability model framework to describe the pore pressure-dependent gas transport property in the shale matrix.
We know that in CO2 capture in coal seams, carbon dioxide acts as a “plasticizer” for coal, lowering the temperature required to cause the transition from a brittle structure to a plastic structure after CO2 is adsorbed by the coal matrix [37]. Permeability can be thus affected by the sequent changes of shale structure due to high stress and matrix swelling/shrinkage. In the shale layer, kerogen softening was observed quite significantly with CO2 injection (supercritical CO2 in many cases) employed as an enhanced recovery process [38]. So, our measured permeability data increase with CO2 injection at higher pressure can be explained, since this enhancement was found when at 295 K when CO2 has encountered phase change—CO2 transforms from gaseous condition to supercritical condition. Supercritical CO2 can result in the decrease in its maximum adsorption amount compared to gaseous condition [39]. With less adsorption, the shale matrix tends to shrink and permeability will thus increase. This observation proves that gas adsorption behavior and sorption-induced matrix deformation in the shale matrix are very critical in describing shale gas apparent permeability evolution. In practice, the matrix swelling strain associated with gas adsorption tends to close the cleat aperture and thus decrease permeability [40]. With the increasing upstream pressure from a very low level, a reduction in permeability can be seen until the Darcian flow starts to take control, and permeability increases again with the gradual loss of effective stress. In conclusion, both sorption-induced deformation and diverse flow behavior between gases explain the CO2 permeability trends. We need to couple the influence of sorption-induced rock deformation and the change of flow regimes to account for shale gas production.
5.2. Effective Stress Coefficient
When pressure is larger than 4 MPa, the effective stress coefficient was found to be positive and less than unity, indicating that the shale sample is experimentally more sensitive to changes in confining pressure than changes in pressure. This may result from very low clay content in the rock sample [34]. On the other hand, the absolute value of is larger than unity at low pressure (<4 MPa), meaning that permeability is more sensitive to pore pressure change. The difference in values and monotony of permeability indicate that the “check”-shape or “U”-shape permeability trends can be explained by certain physical theories rather than by random occurrence, as shown in Figure 9. 4 MPa was found to be a good reference pore pressure for permeability turning point in most of our data. Although this number may not be a guaranteed value for other samples or other conditions, this phenomenon has been captured a lot in theory and practice. The relationship between effective stress and gas apparent permeability contains critical information, not only about tectonic stress change and geological deformation but also about the pore pressure response that reflects how multiple flow regimes influence permeability evolution at varying pressure conditions. All these findings from shale permeability data and effective stresses lead to the conclusion that it is critical to consider the concurrent effect of both geomechanical deformation and macro/microflows.

6. Conclusions
A laboratory study has been conducted to investigate the gas transport in thin shale disk samples through CO2 injections. The influence of effective stress and sorption-induced matrix swelling/shrinkage on permeability was investigated. The following conclusions are made and summarized: (1)The non-Darcy flow effect, basically the Knudsen flow and transition flow, has significant influence on shale gas behavior and overall gas deliverability at low reservoir pressure and in the late production stage(2)The relationship between permeability and pore pressure results in both positive and negative slopes when plotted. The negative slope indicates that the Knudsen flow effect has more significant enhancement on permeability than the reduction from increasing effective stress. Thus, bivalued effective stress coefficients (positive and negative) should be applied to a nonmonotonic pressure-permeability evolution profile(3)For “U”-shape permeability evolution, the effective stress coefficient can be estimated by partitioning permeability into two regions at the turning point. This bivalued effective stress coefficient approach can use the effective stress law concept to explain and model the Knudsen flow at the low-pressure region. In the study, it was found that 4 MPa is the turning point, and we obtained bivalued coefficients that modeled experimental data well
Nomenclature
: | Cross-sectional area of the sample (m2) |
: | Ratio of the sample’s storage capacity to that of the upstream reservoir (dimensionless) |
: | Ratio of the sample’s storage capacity to that of the downstream reservoir (dimensionless) |
: | Gas compressibility (Pa-1) |
: | Molecular diameter (m) |
: | Mass flow correction factor (dimensionless) |
: | Sample length (m) |
: | Gas apparent permeability (m2) |
: | Boltzmann constant () |
: | Pore pressure (Pa) |
: | Confining pressure (Pa) |
: | Pressure of downstream at time (Pa) |
: | Initial pressure of downstream (Pa) |
: | Langmuir pressure (Pa) |
: | Pressure of upstream at time (Pa) |
: | Initial pressure of upstream (Pa) |
: | Pore radius (m) |
: | Temperature (K) |
: | Gas adsorbed volume (cm3·g-1) |
: | Downstream reservoir volume (m3) |
: | Langmuir volume (cm3·g-1) |
: | Effective total sample pore volume (m3) |
: | Mole volume of gas at standard temperature (273.15 K) and pressure (101 325 Pa) () |
: | Upstream reservoir volume (m3). |
: | Slope of the line when plotting the pressure decay against time on a semilog plot (dimensionless) |
: | Dimensionless pressure difference (dimensionless) |
: | The th root of the transcendental equation () (dimensionless) |
: | Gas mean free path (m) |
: | Gas viscosity (Pa·s) |
: | Gas density (kg/m3) |
: | Density of the solid adsorbent (kg/m3) |
: | Effective stress (Pa) |
: | Terzaghi effective stress (Pa) |
: | Matrix porosity (dimensionless) |
: | Effective stress coefficient (dimensionless). |
Data Availability
The data used to support the findings of this study are available from the first author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Acknowledgments
This work is supported in part by the Beijing Municipal Natural Science Foundation (No. 8214063), the Fundamental Research Funds for the Central Universities (No. 2021YQNY08), and the National Natural Science Foundation of China (Nos. 51861145403 and 51874312).