# Experimental Study on Gas Transport in Shale Matrix with Real Gas and Klinkenberg Effects

Authors: Xiaojun Wu; Qing Wang; Fangtao Lyu; Zhengfu Ning; Zongxiao Ren
Journal: Geofluids
 (2021)

---

## Abstract

Gas transport in shale matrix is complex due to multiple mechanisms and is difficult to be investigated by macroscopic experiment. For Gas Research Institute (GRI) method, which is the most accepted one for gas transport investigation in shale matrix, the apparatus was modified by adding an automatic gas supplement and pressurization (AGSP) system, and a numerical model considering the variation of real gas property and the Klinkenberg effect was established for data interpretation. Then, the intrinsic permeability and Klinkenberg coefficient were effectively obtained by maintaining high expanding speed of gas in apparatus and eliminating the negative effect of low filling degree of sample. By analysis, the ideal gas transports faster than real gas due to the viscosity difference at low pressure and the deviation factor difference at high pressure. For Wufeng-Longmaxi shale matrix, the positive influence of Klinkenberg effect on gas transport would attenuate with increasing pressure and is more powerful than bulk shale sample with fractures. Therefore, the gas transport in real shale matrix could be well known, which is meaningful to production forecast and evaluation in oil and gas fields.


---

## Body

## 1. Introduction

The cognition of gas transport in shale matrix has direct connection to the production of shale gas [1–3], which is significant to the extraction of shale gas reservoir. The ultralow porosity and permeability characteristics of shale matrix [4–6] and the multiple transport mechanisms (real gas effect and Klinkenberg effect [7, 8]) all lead to the complexity of gas transport, which deserves a comprehensive study.By experiment, permeability could be obtained, to reveal the visible gas transport phenomena of shale sample in macroscopic view [9]. Steady-state and non-steady-state experiments were the main experimental methods to investigate gas transport in shale. For steady-state experiment, the permeability of shale core plug was tested at constant pore pressure and overburden stress conditions [10]. Due to the ultralow permeability of shale, much long measure time and persistent gas leakage observation were both required. For non-steady-state experiment, a pressure pulse was applied to the shale pore plug [7, 11–13], or crushed shale rock [11, 14–17], and the permeability of shale could be obtained by analyzing the variable pressure data. Fracture was found to be well developed in tight reservoir [18, 19]. The permeability of core plug is much larger than the particles crushed from the same core plug, owing to the less microcrack of particle [20]. The non-steady-state experiment with crushed particles, called Gas Research Institute (GRI) method, was the most recognized one in studying the gas transport process in shale matrix [13, 15, 21, 22], due to the fact that the particle was small enough to reduce the existence of microcracks greatly.However, the GRI experimental method has three key arguments for gas transport investigation in shale matrix. The first one is the absence of confined pressure [9, 10], for studying the effect of stress sensitivity on gas transport. But in reverse, the other transport mechanisms, e.g., real gas effect and Klinkenberg effect, could be studied obviously without the interference of overburden stress.The second argument is the inaccurate measurement to obtain the experimental pressure drop. The gas expansion from reference chamber (RC) to SC was believed to be an immediate process [11]. After the expansion, the gas would enter the particle slowly, and the target pressure drop could be measured. The filling degree of sample in the sample chamber (SC) should be as high as possible to minimize the void volume in SC [20], for shortening the expanding time from RC to SC. Nevertheless, the volume of RC (VRC) was set to be much smaller than the volume of SC (VSC) [23], which would extend the expanding time from RC to SC. Besides, in the researches with traditional GRI method, the experimental pressures (lower than 1 MPa in general) were much lower than 5 MPa, which is the smallest reservoir pressure of shale gas in North America [24]. The traditional GRI method may have difficulty to simulate the gas transport at reservoir conditions. Lower pressure could cause longer expanding time from RC to SC. Selecting low pressure may be aimed at reducing the pressure drop caused by gas entering the particles, to regard the property of gas as constant [25], although small pressure drop is affected easily by pressure fluctuating [26]. In summary, the comprehensive setting of filling degree, ratio of VRC and VSC, and experimental pressure are difficult to ensure the immediate gas expansion from RC to SC.The last argument is the data analysis approaches. The analytical model was constructed based on Darcy’s law for matrix permeability by data analysis [10, 11, 20]. Analytical model has been further modified to more complicate data analysis cases [13, 15, 17]. The analytical model should be simplified to obtain the analytic solution by the assumption of constant properties of rock and gas and the necessary estimation at specific conditions [11]. The assumption and estimation would affect the accuracy of shale matrix permeability [21, 27]. The numerical approaches [16, 21, 25, 28, 29] considered the variable apparent matrix permeability with Klinkenberg effect and were applied to the data analysis of GRI method for accurate gas transport research. It was pointed out that, in analytical and numerical approaches, the variety of gas properties should not be neglected for accurate permeability [27, 30], especially for high experimental pressure. With the consideration of variable gas properties, the implicit numerical solution is hard to be obtained, while stability analysis is necessary to the explicit numerical solution [23]. However, the mathematical methods above did not consider the variable properties of real gas effect and the Klinkenberg effect, synthetically, to maintain a more accurate gas transport simulation in shale matrix.Real gas effect and Klinkenberg effect both possess significant bearings on gas transport in shale matrix. The molecular interaction of real gas should not be neglected at a certain pressure condition [31]. Great difference would appear between the transports of real gas and ideal gas in single nanopore [7, 31–33]. In GRI method, ideal gas was utilized to analyze the experimental data instead of ideal gas [23]. The impact of real gas property has not caused much attention for determining shale matrix permeability. Klinkenberg effect of gas transport in shale has been studied a lot by experimental methods [25, 26, 34–36]. It was summarized that the Klinkenberg effect is generally essential under 5 MPa [37]. These experiments focused on the shale sample with fractures, rather than the shale matrix. In the traditional GRI researches, extremely low experimental pressure was usually set to obtain the permeability [11, 16, 22]. Few GRI researches studied the permeability of shale matrix with Klinkenberg effect at high pressure. Thereby, the GRI method should consider the real gas effect and Klinkenberg effect at wide pressure range, for an accurate and thorough investigation of these two effects.Herein, the helium was selected as a nonadsorptive gas instead of methane to investigate the gas transport firstly. Utilizing nonadsorptive gas, the transport phenomena with Klinkenberg effect could be better observed with the interference of adsorption mechanism [26, 38–40]. As a crucial assistant experimental method, helium was also commonly applied to investigate the transport of methane in shale sample [41–43]. The studies of nonadsorptive gas transport could be a theoretical foundation to the future studies of adsorptive gas transport at reservoir conditions. For accurate real gas properties of helium, the National Institute of Standards and Technology recommended the equations of compressibility [44], viscosity [45], gas deviation factor, and density [46, 47].Because of the arguments around GRI method, it is necessary to modify the original GRI method to investigate the gas transport in shale matrix with real gas effect and Klinkenberg effect. In this research, with the addition of automatic gas supplement and pressurization (AGSP) system, the experimental apparatus could be utilized for high-pressure experiment easily. The influences of low filling degree and unreasonable volume of RC and SC on gas expanding speed were eliminated. Correspondingly, a numerical model with dynamic gas property and dynamic apparent permeability was established for gas transport with real gas effect and Klinkenberg effect. With the modified apparatus and model, the intrinsic permeability and Klinkenberg coefficient of Wufeng-Longmaxi shale matrix were obtained. Furthermore, in Wufeng-Longmaxi shale matrix, the transports of real and ideal gases were compared and analyzed, and the impact of Klinkenberg effect on gas transport was investigated at different pressure conditions. This research could provide a reliable approach to investigate the shale gas transport property in real shale sample, as a theoretical basis for the development of shale gas reservoir.


## 2. Samples and Experimental Methods

### 2.1. Samples

The three shale samples are from Wufeng-Longmaxi Formation in Changning area of Sichuan Province, China [48]. To attenuate the effect of cracks on permeability and investigate the gas transport in matrix of shale sample, three shale rocks were crushed to 20~30-mesh particle samples [11, 49]. Each particle is assumed to be spherical with the same size [11]. The volume of each particle could be calculated by sphere volume formula. For 20~30 mesh particle, which is 0.6~0.85 mm in diameter size, by taking the middle value, the particle diameter utilized was 0.725 mm for the calculation of the total volume and surface area of particles. Each sample was dried to remove water in a drying oven for 12 hours at 105°C condition before experiment. Then, weigh the dried sample. The total volume of particle samples, Vpar, could be obtained by weight divided by density.
(1)Vpar=Mρs,whereM is the molar mass and ρs is the density of sample. The density was obtained by the standard method from GRI [50]. The number of particles could be obtained by the total volume divided by the single particle volume.
(2)N=3Vpar4πrpar3,whereN is the number of uniform particles and rpar is the radius of particle. Then, the total external surface area of samples could be calculated as follows:
(3)S=t4πNrpar2,whereSt is the total external surface area. The densities, total volumes, and total external surface areas of three particle samples are found in Table 1.Table 1
Densities, volumes, and filling degrees of three shale samples in SC.
SampleDensity (g/ml)Vpar (ml)Vt (ml)St (m2)Filling degree (%)12.57216.63467.0230.550638.36522.56616.63567.0220.550738.36732.63617.05266.6050.564539.331


### 2.2. Experimental Apparatus

As shown in Figure1, the gas cylinder (1), booster pump (2), and regulating valve (3) were in series connection. With the three-way valves (4 and 6), the AGSP system and intake valve (5) were connected in parallel. The AGSP system (Figure 1) is consisted of pneumatic valve (7), solenoid valve (8), pressure controller (9), and pressure transmitter (10). Then, reference chamber (RC) (11), pressure sensor (14), intermediate valve (15), and SC (13) were connected in turn. SC was sealed by flange and sealing ring. The RC, SC, and the gas tank utilized for booster pump were all enclosed in the temperature controlling device (12). The vacuum system and discharging valve (17) were connected in parallel. The vacuum system was consisted by vacuum valve (16) and vacuum (18). Each part was connected by gas pipelines with gas tightness guarantee. At the junction between pipeline and the flange of sample seal, strainer screen is necessary in case of sample entering to the pipeline. Wires were utilized to connect pressure transmitter, pressure controller, and solenoid valve in turn. The AGSP system could effectively reduce measurement error by shortening the time of gas expansion from RC to SC. The detailed discussion would be presented in Section 2.4.Figure 1
The equipment for shale matrix permeability determination with the AGSP system.The volumes of RC and SC could be measured by putting different numbers of stainless-steel balls in SC before the helium expansion from RC to SC [48]. Actually, at the condition that all valves were closed in advance, the volumes of two chambers include the volumes of pipelines connected directly with chambers. Then, with the obtained total volume of particles (Vpar in Table 1), the filling degree in SC, the void volume out of the particles, and the total surface area of particles could be obtained and given in Table 1.


### 2.3. Experimental Procedure

The experimental steps of should be modified for the utilization of AGSP system. Thus, the designed procedure was provided as below.(1)
Keep the intake valve and intermediate valve open and close the discharging valve and vacuum valve. Fill a certain amount of helium to the RC and SC and close the intake valve immediately. The gas tightness of the whole device could be evaluated based on the variation of pressure sensor value in 60 minutes(2)
Place the dry shale particles into the SC. Keep the intake valve and discharging valve closed and open the intermediate valve and vacuum valve to vacuumize the whole device and shale sample for 12 hours. Keep the temperature controlling device running to maintain the temperature at target value. Close all the valves and vacuum after evacuation(3)
Run booster pump and keep the pressure of gas tank with high pressure (larger than 20 MPa). Open the intake valve and the regulating valve to make pressure of RC rise top1. Close the intake valve again and turn on the power of AGSP system. Set the target experimental pressure p1 in the pressure controller and the pressure differential triggering the automatic gas supplement(4)
With the intermediate valve opening, the gas expands from RC to SC. The pressure would drop rapidly and AGSP system starts to work immediately by detecting the pressure differential. As soon as the pressure recovers top1, the AGSP system switches off automatically and the gas supplement would be stopped. Then, cut off the power of AGSP system(5)
With the gas entering to the particles, the computer is utilized to record the pressure dropping with time by pressure sensor. When the pressure is finally stable atp2, the data recording is completed. Then, close the regulating valve, and empty the gas from the whole device by discharging valve


### 2.4. Gas Supplement and Pressurization

In step 4, with the intermediate valve open, the gas in RC expands to SC, as shown in Figure2. For the original GRI method without AGSP system, no gas supplement was presented, and the amount of gas was the same before and after opening the intermediate valve. The pressure in the connected RC and SC would drop dramatically.Figure 2
The function of AGSP system in the experiment.The dramatic drop may cause several concerns. The first concern is that due to the dropping of pressure, the velocity of gas expansion would drop correspondingly. With gas expanding from RC to SC, the gas is also expanding to the particles at the same time. The pressure sensor could only detect the pressure drop caused by gas expanding to the particles after the stopping of the gas expansion from RC to SC. Low velocity of gas expansion from RC to SC would interfere the pressure sensor to detect the pressure drop caused by gas expansion to the particles. This may cause inevitable measurement error. In addition, the low filling degree [16] or unreasonable volumes of RC and SC could also affect the measurement results [23], owing to the caused low expansion velocity from RC to SC. Another concern is that the expansion of fixed amount of gas from RC to SC may induce a low experimental pressure, as shown in Figure 2. The phenomena of shale gas transport are inconsistent at different pressure conditions. Previously, few experiments were conducted with at the pressure higher than 1 MPa. It is necessary to investigate the gas transport at high pressure condition.To avoid these errors, the expansion velocity from RC to SC should be improved and the corresponding expansion process should be short enough. A solution of this issue is to add a system to provide a much higher pressure gradient to ensure the high expansion velocity and a gas supplement to maintain a high experimental pressure.In this research, the AGSP system is added to implement this functionality. In step 4, as soon as opening the intermediate valve, the pressure transmitter detects the pressure drop and the message would be transmitted to the pressure controller. Then, the pressure controller makes the pneumatic valve open by solenoid valve. The gas container of booster pump with much higher pressure starts to supply gas. Thus, the pressure in RC and SC recovers to the experimental pressure in a much short time, shown in Figure2. The pressure transmitter detects the pressure recovery and the pneumatic valve closes again. After that, the slow pressure drop caused by expansion to the particle sample would be well observed. With the AGSP system, the low filling degree and unreasonable volumes of RC and SC would not affect the expansion time from RC and SC. The comparison of experimental results with and without AGSP system would be discussed in detail in Section 4.2.


## 2.1. Samples

The three shale samples are from Wufeng-Longmaxi Formation in Changning area of Sichuan Province, China [48]. To attenuate the effect of cracks on permeability and investigate the gas transport in matrix of shale sample, three shale rocks were crushed to 20~30-mesh particle samples [11, 49]. Each particle is assumed to be spherical with the same size [11]. The volume of each particle could be calculated by sphere volume formula. For 20~30 mesh particle, which is 0.6~0.85 mm in diameter size, by taking the middle value, the particle diameter utilized was 0.725 mm for the calculation of the total volume and surface area of particles. Each sample was dried to remove water in a drying oven for 12 hours at 105°C condition before experiment. Then, weigh the dried sample. The total volume of particle samples, Vpar, could be obtained by weight divided by density.
(1)Vpar=Mρs,whereM is the molar mass and ρs is the density of sample. The density was obtained by the standard method from GRI [50]. The number of particles could be obtained by the total volume divided by the single particle volume.
(2)N=3Vpar4πrpar3,whereN is the number of uniform particles and rpar is the radius of particle. Then, the total external surface area of samples could be calculated as follows:
(3)S=t4πNrpar2,whereSt is the total external surface area. The densities, total volumes, and total external surface areas of three particle samples are found in Table 1.Table 1
Densities, volumes, and filling degrees of three shale samples in SC.
SampleDensity (g/ml)Vpar (ml)Vt (ml)St (m2)Filling degree (%)12.57216.63467.0230.550638.36522.56616.63567.0220.550738.36732.63617.05266.6050.564539.331


## 2.2. Experimental Apparatus

As shown in Figure1, the gas cylinder (1), booster pump (2), and regulating valve (3) were in series connection. With the three-way valves (4 and 6), the AGSP system and intake valve (5) were connected in parallel. The AGSP system (Figure 1) is consisted of pneumatic valve (7), solenoid valve (8), pressure controller (9), and pressure transmitter (10). Then, reference chamber (RC) (11), pressure sensor (14), intermediate valve (15), and SC (13) were connected in turn. SC was sealed by flange and sealing ring. The RC, SC, and the gas tank utilized for booster pump were all enclosed in the temperature controlling device (12). The vacuum system and discharging valve (17) were connected in parallel. The vacuum system was consisted by vacuum valve (16) and vacuum (18). Each part was connected by gas pipelines with gas tightness guarantee. At the junction between pipeline and the flange of sample seal, strainer screen is necessary in case of sample entering to the pipeline. Wires were utilized to connect pressure transmitter, pressure controller, and solenoid valve in turn. The AGSP system could effectively reduce measurement error by shortening the time of gas expansion from RC to SC. The detailed discussion would be presented in Section 2.4.Figure 1
The equipment for shale matrix permeability determination with the AGSP system.The volumes of RC and SC could be measured by putting different numbers of stainless-steel balls in SC before the helium expansion from RC to SC [48]. Actually, at the condition that all valves were closed in advance, the volumes of two chambers include the volumes of pipelines connected directly with chambers. Then, with the obtained total volume of particles (Vpar in Table 1), the filling degree in SC, the void volume out of the particles, and the total surface area of particles could be obtained and given in Table 1.


## 2.3. Experimental Procedure

The experimental steps of should be modified for the utilization of AGSP system. Thus, the designed procedure was provided as below.(1)
Keep the intake valve and intermediate valve open and close the discharging valve and vacuum valve. Fill a certain amount of helium to the RC and SC and close the intake valve immediately. The gas tightness of the whole device could be evaluated based on the variation of pressure sensor value in 60 minutes(2)
Place the dry shale particles into the SC. Keep the intake valve and discharging valve closed and open the intermediate valve and vacuum valve to vacuumize the whole device and shale sample for 12 hours. Keep the temperature controlling device running to maintain the temperature at target value. Close all the valves and vacuum after evacuation(3)
Run booster pump and keep the pressure of gas tank with high pressure (larger than 20 MPa). Open the intake valve and the regulating valve to make pressure of RC rise top1. Close the intake valve again and turn on the power of AGSP system. Set the target experimental pressure p1 in the pressure controller and the pressure differential triggering the automatic gas supplement(4)
With the intermediate valve opening, the gas expands from RC to SC. The pressure would drop rapidly and AGSP system starts to work immediately by detecting the pressure differential. As soon as the pressure recovers top1, the AGSP system switches off automatically and the gas supplement would be stopped. Then, cut off the power of AGSP system(5)
With the gas entering to the particles, the computer is utilized to record the pressure dropping with time by pressure sensor. When the pressure is finally stable atp2, the data recording is completed. Then, close the regulating valve, and empty the gas from the whole device by discharging valve


## 2.4. Gas Supplement and Pressurization

In step 4, with the intermediate valve open, the gas in RC expands to SC, as shown in Figure2. For the original GRI method without AGSP system, no gas supplement was presented, and the amount of gas was the same before and after opening the intermediate valve. The pressure in the connected RC and SC would drop dramatically.Figure 2
The function of AGSP system in the experiment.The dramatic drop may cause several concerns. The first concern is that due to the dropping of pressure, the velocity of gas expansion would drop correspondingly. With gas expanding from RC to SC, the gas is also expanding to the particles at the same time. The pressure sensor could only detect the pressure drop caused by gas expanding to the particles after the stopping of the gas expansion from RC to SC. Low velocity of gas expansion from RC to SC would interfere the pressure sensor to detect the pressure drop caused by gas expansion to the particles. This may cause inevitable measurement error. In addition, the low filling degree [16] or unreasonable volumes of RC and SC could also affect the measurement results [23], owing to the caused low expansion velocity from RC to SC. Another concern is that the expansion of fixed amount of gas from RC to SC may induce a low experimental pressure, as shown in Figure 2. The phenomena of shale gas transport are inconsistent at different pressure conditions. Previously, few experiments were conducted with at the pressure higher than 1 MPa. It is necessary to investigate the gas transport at high pressure condition.To avoid these errors, the expansion velocity from RC to SC should be improved and the corresponding expansion process should be short enough. A solution of this issue is to add a system to provide a much higher pressure gradient to ensure the high expansion velocity and a gas supplement to maintain a high experimental pressure.In this research, the AGSP system is added to implement this functionality. In step 4, as soon as opening the intermediate valve, the pressure transmitter detects the pressure drop and the message would be transmitted to the pressure controller. Then, the pressure controller makes the pneumatic valve open by solenoid valve. The gas container of booster pump with much higher pressure starts to supply gas. Thus, the pressure in RC and SC recovers to the experimental pressure in a much short time, shown in Figure2. The pressure transmitter detects the pressure recovery and the pneumatic valve closes again. After that, the slow pressure drop caused by expansion to the particle sample would be well observed. With the AGSP system, the low filling degree and unreasonable volumes of RC and SC would not affect the expansion time from RC and SC. The comparison of experimental results with and without AGSP system would be discussed in detail in Section 4.2.


## 3. Data Interpretation Methods

### 3.1. Porosity Determination

Based on the mass conservation before and after gas entering to the particles, the porosity in this experiment could be given by(4)ϕ=pout‐ini−pin‐inipe−pin‐ini−1VRC+VSC−VparVpar,whereϕ is the porosity of sample, pin‐ini and pout‐ini are the initial pressures in and out of the particles at the beginning of pressure drop, and pe is the final pressure at the end of pressure drop. Due to the vacuum treatment, pin‐ini is zero in general.Herein, a dimensionless pressurepD was imported here for better comparing the drop curves of pout at different pin‐ini and pout‐ini pressure conditions. pout refers to the variable pressure out of the particles. The expression of pD is presented as below.
(5)pD=pout−pin‐inipe−pin‐ini.At the beginning of pressure drop,pout in Equation (5) equals to pout‐ini; the values of pD were the same for different pin‐ini and pout‐ini based on Equation (4). At the end of pressure drop, pout in Equation (5) equals to pe, and pD equals to 1 at different conditions. Besides, pD has the same drop tendency with pout.


### 3.2. Gas Transport Models

#### 3.2.1. Analytical Model

Each particle of the shale sample in SC could be assumed to be spherical. For sphere rock, the governing equation could be expressed as [11]
(6)∂ρϕ∂t=1r2∂∂rr2ρkμ∂p∂r,whereρ is the density of gas, t is the transport time, r is the transport distance, p is the pressure, k is the apparent permeability, and μ is the viscosity of fluid.For real gas, the compressibility factor,Cg, and density have the following expressions.
(7)Cg=1ρ∂ρ∂p,(8)ρ=MpZRT,whereZ is the deviation factor and R is the universal gas constant.Based on Equations (7) and (8), Equation (6) could be rewritten as
(9)∂p/Zϕ∂t=1r2∂∂rr2kμCg∂p/Z∂r.With the gas expanding to the particle, the pressure in the particle would rise from zero to the final experimental pressure. In the large scale of pressure changing, it is unreasonable to utilize constant physical parameters,Z, μ, Cg, and k, for different pressures. Consequently, the physical parameters should not be taken out from the partial differential terms, unlike the previous models [11, 21, 23]. If the physical parameters Z, μ, Cg, and k are constant, Equation (9) could be converted to the equation utilized in previous research [11].For real gas, the physical parameters in Equation (9), i.e., Z, μ, and Cg, are all variable with pressure and temperature. For an accurate simulation of real gas at a determined temperature, Z, μ, and Cg at a certain pressure in this research were directly calculated based on the equations recommended by the National Institute of Standards and Technology [44–46], instead of simple functions of pressure in numerical model [14].Due to the Klinkenberg effect, the apparent permeability varies with pressure and has great impact on the gas transport. The apparent permeability with Klinkenberg effect could be expressed as(10)k=kint1+bp,wherekint is the intrinsic permeability and b is the Klinkenberg coefficient.At initial stage, when the gas just stops expansion from RC to SC, the pressure out of the particles equals to the initial experiment pressure, while the pressure in particles could be seemed as zero due to vacuum and short gas expanding time from RC to SC. Then, the initial condition is shown as below.(11)pin=0,pout=pout‐ini,wherepin is the variable pressure in the particles.Owing to the spatial symmetry of sphere, the center of the particle could be regarded as inner sealing boundary. As a result, the inner boundary condition is(12)∂p∂rr=0=0.Based on the law of conservation of mass, the loss of gas out of the particles is equivalent to the amount of gas entering the particles. Thus, the outer boundary condition is(13)Vt∂p/Z∂t=−StkμCg∂p/Z∂r,whereVt is the volume of the external space of particles with the intermediate valve open, which are given in Table 1.


#### 3.2.2. Numerical Model with Dynamic Physical Properties

An accurate analytical solution is difficult to be acquired for the above analytical model with dynamic physical properties. By simplification, a simplified analytical solution could be derived for data interpretation [11, 13, 17]. The simplified analytical solution did not consider the significant variation of physical properties with variable conditions [27]. Besides, Cui et al. pointed out that the early-time and late-time analytical solutions could result in different data interpretation results [11]. Therefore, numerical solution is recommended for an accurate data interpretation [14].As a result, due to the complicated variable parameters, i.e.,Z, μ, Cg, and k, it is hard to acquire a solution based on the numerical model with implicit difference. Herein, a numerical model with explicit difference was introduced. Its stability would be discussed in Section 4.1.2 to make sure it is accurate and convincible. As shown in Figure 3, a spherical particle with radius, R, was divided into m spherical shells with identical thickness, Δr. By taking the outer boundary of each shell as the grid point from the inside out, the grid points were denoted from r1 to rm. For a grid point ri, the centers of the shells before and after were ri−1/2 and ri+1/2. The number of time steps is n with regular time interval, Δt.Figure 3
The numerical difference of a single particle.Treating Equation (9) by time difference for the left and spatial difference for the right and arranging the unknown to the left, the discretize governing equation could be given as below
(14)pin+1=ZinΔtϕiΔr2i+122ki+1/2nμi+1/2nCgi+1/2npZi+1n−pZin−i−122ki−1/2nμi−1/2nCgi−1/2npZin−pZi−1n+pin.In Equation (14), the parameters μ, Cg, and k at grid points ri−1/2 and ri+1/2 could be treated by harmonic mean method. The detailed expressions, Equations (A.1)–(A.6), were presented in the appendix.After difference discretization, the initial, inner boundary, and outer boundary conditions could be written as below.(15)pi1=01≤i<m−1,pm1=pout‐ini,p1n=p2n,pmn+1=pmn−StΔtVtΔrZkμCgmnpZmn−pZm−1n.


## 3.1. Porosity Determination

Based on the mass conservation before and after gas entering to the particles, the porosity in this experiment could be given by(4)ϕ=pout‐ini−pin‐inipe−pin‐ini−1VRC+VSC−VparVpar,whereϕ is the porosity of sample, pin‐ini and pout‐ini are the initial pressures in and out of the particles at the beginning of pressure drop, and pe is the final pressure at the end of pressure drop. Due to the vacuum treatment, pin‐ini is zero in general.Herein, a dimensionless pressurepD was imported here for better comparing the drop curves of pout at different pin‐ini and pout‐ini pressure conditions. pout refers to the variable pressure out of the particles. The expression of pD is presented as below.
(5)pD=pout−pin‐inipe−pin‐ini.At the beginning of pressure drop,pout in Equation (5) equals to pout‐ini; the values of pD were the same for different pin‐ini and pout‐ini based on Equation (4). At the end of pressure drop, pout in Equation (5) equals to pe, and pD equals to 1 at different conditions. Besides, pD has the same drop tendency with pout.


## 3.2. Gas Transport Models

### 3.2.1. Analytical Model

Each particle of the shale sample in SC could be assumed to be spherical. For sphere rock, the governing equation could be expressed as [11]
(6)∂ρϕ∂t=1r2∂∂rr2ρkμ∂p∂r,whereρ is the density of gas, t is the transport time, r is the transport distance, p is the pressure, k is the apparent permeability, and μ is the viscosity of fluid.For real gas, the compressibility factor,Cg, and density have the following expressions.
(7)Cg=1ρ∂ρ∂p,(8)ρ=MpZRT,whereZ is the deviation factor and R is the universal gas constant.Based on Equations (7) and (8), Equation (6) could be rewritten as
(9)∂p/Zϕ∂t=1r2∂∂rr2kμCg∂p/Z∂r.With the gas expanding to the particle, the pressure in the particle would rise from zero to the final experimental pressure. In the large scale of pressure changing, it is unreasonable to utilize constant physical parameters,Z, μ, Cg, and k, for different pressures. Consequently, the physical parameters should not be taken out from the partial differential terms, unlike the previous models [11, 21, 23]. If the physical parameters Z, μ, Cg, and k are constant, Equation (9) could be converted to the equation utilized in previous research [11].For real gas, the physical parameters in Equation (9), i.e., Z, μ, and Cg, are all variable with pressure and temperature. For an accurate simulation of real gas at a determined temperature, Z, μ, and Cg at a certain pressure in this research were directly calculated based on the equations recommended by the National Institute of Standards and Technology [44–46], instead of simple functions of pressure in numerical model [14].Due to the Klinkenberg effect, the apparent permeability varies with pressure and has great impact on the gas transport. The apparent permeability with Klinkenberg effect could be expressed as(10)k=kint1+bp,wherekint is the intrinsic permeability and b is the Klinkenberg coefficient.At initial stage, when the gas just stops expansion from RC to SC, the pressure out of the particles equals to the initial experiment pressure, while the pressure in particles could be seemed as zero due to vacuum and short gas expanding time from RC to SC. Then, the initial condition is shown as below.(11)pin=0,pout=pout‐ini,wherepin is the variable pressure in the particles.Owing to the spatial symmetry of sphere, the center of the particle could be regarded as inner sealing boundary. As a result, the inner boundary condition is(12)∂p∂rr=0=0.Based on the law of conservation of mass, the loss of gas out of the particles is equivalent to the amount of gas entering the particles. Thus, the outer boundary condition is(13)Vt∂p/Z∂t=−StkμCg∂p/Z∂r,whereVt is the volume of the external space of particles with the intermediate valve open, which are given in Table 1.


### 3.2.2. Numerical Model with Dynamic Physical Properties

An accurate analytical solution is difficult to be acquired for the above analytical model with dynamic physical properties. By simplification, a simplified analytical solution could be derived for data interpretation [11, 13, 17]. The simplified analytical solution did not consider the significant variation of physical properties with variable conditions [27]. Besides, Cui et al. pointed out that the early-time and late-time analytical solutions could result in different data interpretation results [11]. Therefore, numerical solution is recommended for an accurate data interpretation [14].As a result, due to the complicated variable parameters, i.e.,Z, μ, Cg, and k, it is hard to acquire a solution based on the numerical model with implicit difference. Herein, a numerical model with explicit difference was introduced. Its stability would be discussed in Section 4.1.2 to make sure it is accurate and convincible. As shown in Figure 3, a spherical particle with radius, R, was divided into m spherical shells with identical thickness, Δr. By taking the outer boundary of each shell as the grid point from the inside out, the grid points were denoted from r1 to rm. For a grid point ri, the centers of the shells before and after were ri−1/2 and ri+1/2. The number of time steps is n with regular time interval, Δt.Figure 3
The numerical difference of a single particle.Treating Equation (9) by time difference for the left and spatial difference for the right and arranging the unknown to the left, the discretize governing equation could be given as below
(14)pin+1=ZinΔtϕiΔr2i+122ki+1/2nμi+1/2nCgi+1/2npZi+1n−pZin−i−122ki−1/2nμi−1/2nCgi−1/2npZin−pZi−1n+pin.In Equation (14), the parameters μ, Cg, and k at grid points ri−1/2 and ri+1/2 could be treated by harmonic mean method. The detailed expressions, Equations (A.1)–(A.6), were presented in the appendix.After difference discretization, the initial, inner boundary, and outer boundary conditions could be written as below.(15)pi1=01≤i<m−1,pm1=pout‐ini,p1n=p2n,pmn+1=pmn−StΔtVtΔrZkμCgmnpZmn−pZm−1n.


## 3.2.1. Analytical Model

Each particle of the shale sample in SC could be assumed to be spherical. For sphere rock, the governing equation could be expressed as [11]
(6)∂ρϕ∂t=1r2∂∂rr2ρkμ∂p∂r,whereρ is the density of gas, t is the transport time, r is the transport distance, p is the pressure, k is the apparent permeability, and μ is the viscosity of fluid.For real gas, the compressibility factor,Cg, and density have the following expressions.
(7)Cg=1ρ∂ρ∂p,(8)ρ=MpZRT,whereZ is the deviation factor and R is the universal gas constant.Based on Equations (7) and (8), Equation (6) could be rewritten as
(9)∂p/Zϕ∂t=1r2∂∂rr2kμCg∂p/Z∂r.With the gas expanding to the particle, the pressure in the particle would rise from zero to the final experimental pressure. In the large scale of pressure changing, it is unreasonable to utilize constant physical parameters,Z, μ, Cg, and k, for different pressures. Consequently, the physical parameters should not be taken out from the partial differential terms, unlike the previous models [11, 21, 23]. If the physical parameters Z, μ, Cg, and k are constant, Equation (9) could be converted to the equation utilized in previous research [11].For real gas, the physical parameters in Equation (9), i.e., Z, μ, and Cg, are all variable with pressure and temperature. For an accurate simulation of real gas at a determined temperature, Z, μ, and Cg at a certain pressure in this research were directly calculated based on the equations recommended by the National Institute of Standards and Technology [44–46], instead of simple functions of pressure in numerical model [14].Due to the Klinkenberg effect, the apparent permeability varies with pressure and has great impact on the gas transport. The apparent permeability with Klinkenberg effect could be expressed as(10)k=kint1+bp,wherekint is the intrinsic permeability and b is the Klinkenberg coefficient.At initial stage, when the gas just stops expansion from RC to SC, the pressure out of the particles equals to the initial experiment pressure, while the pressure in particles could be seemed as zero due to vacuum and short gas expanding time from RC to SC. Then, the initial condition is shown as below.(11)pin=0,pout=pout‐ini,wherepin is the variable pressure in the particles.Owing to the spatial symmetry of sphere, the center of the particle could be regarded as inner sealing boundary. As a result, the inner boundary condition is(12)∂p∂rr=0=0.Based on the law of conservation of mass, the loss of gas out of the particles is equivalent to the amount of gas entering the particles. Thus, the outer boundary condition is(13)Vt∂p/Z∂t=−StkμCg∂p/Z∂r,whereVt is the volume of the external space of particles with the intermediate valve open, which are given in Table 1.


## 3.2.2. Numerical Model with Dynamic Physical Properties

An accurate analytical solution is difficult to be acquired for the above analytical model with dynamic physical properties. By simplification, a simplified analytical solution could be derived for data interpretation [11, 13, 17]. The simplified analytical solution did not consider the significant variation of physical properties with variable conditions [27]. Besides, Cui et al. pointed out that the early-time and late-time analytical solutions could result in different data interpretation results [11]. Therefore, numerical solution is recommended for an accurate data interpretation [14].As a result, due to the complicated variable parameters, i.e.,Z, μ, Cg, and k, it is hard to acquire a solution based on the numerical model with implicit difference. Herein, a numerical model with explicit difference was introduced. Its stability would be discussed in Section 4.1.2 to make sure it is accurate and convincible. As shown in Figure 3, a spherical particle with radius, R, was divided into m spherical shells with identical thickness, Δr. By taking the outer boundary of each shell as the grid point from the inside out, the grid points were denoted from r1 to rm. For a grid point ri, the centers of the shells before and after were ri−1/2 and ri+1/2. The number of time steps is n with regular time interval, Δt.Figure 3
The numerical difference of a single particle.Treating Equation (9) by time difference for the left and spatial difference for the right and arranging the unknown to the left, the discretize governing equation could be given as below
(14)pin+1=ZinΔtϕiΔr2i+122ki+1/2nμi+1/2nCgi+1/2npZi+1n−pZin−i−122ki−1/2nμi−1/2nCgi−1/2npZin−pZi−1n+pin.In Equation (14), the parameters μ, Cg, and k at grid points ri−1/2 and ri+1/2 could be treated by harmonic mean method. The detailed expressions, Equations (A.1)–(A.6), were presented in the appendix.After difference discretization, the initial, inner boundary, and outer boundary conditions could be written as below.(15)pi1=01≤i<m−1,pm1=pout‐ini,p1n=p2n,pmn+1=pmn−StΔtVtΔrZkμCgmnpZmn−pZm−1n.


## 4. Results and Discussions

### 4.1. Results

#### 4.1.1. Stability Analysis of Numerical Model

To ensure the stability of numerical model with explicit difference, the following condition should be established for the reasonable time and space intervals [23].
(16)ΔtΔr2≤ϕμCg2k.The permeability should be known in advance to take an effective judgment by Equation (16). Based on the analytical method [11], an estimated permeability could be obtained and help to evaluate the stability of numerical model with explicit difference. Herein, the experimental data of Sample 3 at 8.117 MPa initial pressure and 303.15 K was interpreted by analytical method. As shown in Figure 4, the experimental data was analyzed by late-time (Figure 4(a)) and early-time (Figure 4(b)) methods. It is obvious to find out that the trend lines could not fit well with the experimental data. Based on the slopes of trend lines, the permeabilities from early-time and late-time methods are 0.25×10−23m2 and 1.96×10−23m2 (Table 2). As said by Cui et al., there is great difference between the permeabilities from late-time and early-time methods [11]. It should be noted that, despite the unwell fitting of trend lines, the order of magnitude for calculated permeability is generally the same for different slopes. The unwell fitting maybe due to the ultralow permeability.Figure 4
Curve fitting with (a) late-time and (b) early-time analytical methods at 8.117 MPa and 303.15 K.FR is the mass fraction out of the particle, and FU is the cumulative gas uptake [11].
(a)(b)Table 2
Stability analysis of numerical model based on Sample 3 at 8.177 MPa and 303.15 K conditions.
μ (10-5 Pa·s)Cg (1/MPa)ϕΔt (s)mΔt/Δr2 (108 s/m2)k (early-time)/(10-23 m2)k (late-time)/(10-23 m2)ϕμCg/2k (early-time)/(s/m2)ϕμCg/2k (late-time)/(s/m2)2.030.117890.0380.011007.610.251.961.82×10102.32×109In Equation (16), the values of μ and Cg were calculated at initial pressure condition. The value of porosity of Sample 3 could be calculated by Equation (3). All the values were given in Table 2. By comparison, shown in Table 2, with 0.01 s time interval and 100 grid points, for the permeabilities from late-time and early-time methods, the value in the left of Equation (16) is much smaller than the values in the right. Therefore, 0.01 s time interval and 100 grid points could ensure the stability of numerical model.


#### 4.1.2. Regression Results

As shown in Figure5, the experimental data of 3 and 4 MPa was utilized for the determinations of intrinsic permeability and Klinkenberg coefficient, and the experimental data of 5 MPa was for validation. It should be noted that, in each experiment, it is impossible to set initial value of the pressure drop to be exact 3, 4, or 5 MPa, including the other pressure conditions, owing to the drastic pressure change after the opening of intermediate valve. After regression, the simulation curves fit well with the experimental data.Figure 5
Curve fitting with numerical method at 3 and 4 MPa for Samples (a) 1, (b) 2, and (c) 3 and (d) validation of three samples at 5 MPa. The temperature is 303.15 K.
(a)(b)(c)(d)By regression, the values of porosity, intrinsic permeability, and Klinkenberg coefficient are given in Table3. The porosity was firstly calculated by Equation (4) and adjusted minutely during the curve fitting. The intrinsic permeabilities obtained in this research are close to the values from the original GRI method [22]. The intrinsic permeabilities of Wufeng-Longmaxi shale matrix are appropriate for the GRI method with 20/35-mesh particle size [10]. Moreover, the intrinsic permeabilities obtained by numerical method (Table 3) and analytical method (Table 2) are in the same order of magnitude. The intrinsic permeabilities of Wufeng-Longmaxi shale matrix are much lower than the permeability of shale sample with fractures [43], due to that the permeability of bulk shale was believed to be dominated by fractures [11].Table 3
Regressed parameters for Wufeng-Longmaxi samples.
Sampleϕkint (10-23 m2)b (MPa)10.033.42.320.0211.741.530.0383.12


### 4.2. The Effect of AGSP System

With and without the AGSP system, there is great difference for the pressure drop curves. Keep Sample 3 in SC and set the pressure in RC to be 9 MPa before the opening of intermediate valve. As shown in Figure6(a), as soon as opening the intermediate valve, the pressure would drop instantly. The pressure would take 9 seconds (Δt2 in Figure 6(a)) to balance before drop. In the 9 seconds, a certain amount of gas has entered the vacuumized particle, which affected the accuracy of pressure drop data. On the contrary, with the AGSP system, the gas was supplemented, and pressure came back to 9 MPa in 2 seconds (Δt1 in Figure 6(a)), which ensured the accuracy of pressure drop data. Besides, if the gas supplement is fast enough, the time interval, Δt1, is easy to be undetected.Figure 6
The (a) initial stage and the (b) whole process of pressure drops with and without AGSP system at 303.15 K, with Sample 3.
(a)(b)The pressure drops with and without the AGSP system could be compared effectively by the dimensionless pressure,pD (Figure 6(b)). The initial pD without the AGSP system is much smaller than the one with the AGSP system. This could be explained by the fact that, without AGSP system, large amount of gas has entered the particles before the pressure’s drop. Therefore, the regressed parameters without AGSP would appear large errors. Actually, the balance time without AGSP system has connection with the volume ratio of RC and SC. Small VRC and large VSC could bring large pressure drop, due to the larger capacity to samples, just as the treatment in the previous research [23]. However, this approach could extend the time gas expanding from RC to SC and larger Δt2 would be obtained. With the AGSP system, there is no requirement to the ratio of VRC and VSC, and the large pressure drop could be still realized.The filling degree of samples could affect the dead volume of SC, which determines the gas expanding time. Low filling degree may result in a long expanding time from RC to SC, influencing accuracy of the following pressure drop caused by gas entering the particle. The filling degree should be adjusted to minimize the void volume [20]. At least 75% filling degree was recommended for the minimization of void volume [16]. However, in our experiment, the filling degree of three 20~30-mesh samples is difficult to be larger than 50% and is only 40% around (Table 1). Thus, the low filling degree in our experiment led to a long expanding time from RC to SC, without the help of the AGSP system (Figure 6(a)). The effect of filling degree on the experiment results should be reevaluated with the help of the AGSP system. For 30.83% and 20.66% filling degrees, the experimental and simulative results are consistent with each other (Figure 7). The simulation utilized the regressed parameters of Sample 3 in Table 2. In the low filling degree situation, the gas still expands rapidly to the SC due to the high-pressure gas supplement conducted by the AGSP system. Then, the following pressure drop would be monitored accurately, without the interference of low filling degree.Figure 7
The experimental and simulative results for Sample 3 with 30.83% and 20.66% filling degrees, at 5 MPa and 303.15 K.


### 4.3. Real and Ideal Gases

For real gas, the properties, i.e.,Z, μ, and Cg, vary differently with ideal gas, resulting in the transport difference of real and ideal gas. The parameters of ideal gas in previous research [23] were applied here for a comparison with real gas transport. Same as the previous research, the simulations of real and ideal gases were conducted at 5 MPa, shown in Figure 8. The pressure of ideal gas is lower than that of the real gas at the same time, indicating that the ideal gas transports faster in the shale matrix than real gas. The distinction of real and ideal gas transports has been found in the study of ideal circular nanotubes [31, 33]. For the gas transport in the real shale sample, there is obvious and nonnegligible distinction for real and ideal gas transports in shale matrix, too.Figure 8
The simulation of real and ideal gases and the real gas simulations with ideal parameters, i.e.,Z, μ, and Cg, for Sample 3, at 5.003 MPa and 303.15 K, with Sample 3. The time interval is 50~250 s.In further comparative simulation, one of the real gas parameters,Z, μ or Cg, was replaced by the corresponding ideal parameter, as shown in Figure 8. The ideal parameters were from the previous research [23] and presented in Figure 8. In addition, at the center of particle, based on the variation of real gas pressure with time (Figure 9(a)), the corresponding variations of real and ideal compressibility factors (Figure 9(a)), real deviation factor, and real viscosity (Figure 9(b)) were provided. From Figure 9(a), the real Cg is smaller than the ideal Cg and the difference grows with pressure. In the previous research [23], the ideal Z and μ were constant. But for real gas, Z and μ both increase with pressure and have the similar tendencies with pressure (Figure 9(b)). As a result, the real Z and μ are larger than ideal Z and μ. Comparing the real gas pressure drop curves with and without the ideal parameters (Figure 8), it could be concluded that smaller Cg, Z, and μ could improve the gas transport ability in shale matrix.Figure 9
The variation of real gas pressure with time (a), and the corresponding variations of real and ideal compressibility factors (a), real deviation factor, and real viscosity (b) at 5.003 MPa and 303.15 K, with Sample 3.
(a)(b)With the pressure at the center growing rapidly after 100 s, the real and ideal compressibility factors both decrease quickly, and more obvious deviation appears (Figure9(a)). However, the simulation of real gas with ideal Cg does not deviate a lot from the simulation of real gas (Figure 8). The deviation of real and ideal compressibility factors is not large enough to lead to the great difference of real and ideal gas transports. From Figure 8, the curve of real gas with ideal μ is close to the ideal gas curve firstly and close to the real gas curve later. This is due to the fact that, with the increasing of real μ, the disparity of ideal and real μ is also growing. Conversely, the curve of real gas with ideal Z is close to the curve of ideal gas firstly and close to the curve of real gas later. It agrees with the theory that the interaction of gas molecules would be improved at high pressure and the deviation factor should not be neglected [31, 51]. It indicates that the viscosity and deviation factor, μ and Z, dominate the distinction of real and ideal gas transports at relatively low and high pressure, separately. Therefore, the real gas transport varies with ideal gas transport mainly owing to the effect of Z and μ at different pressure conditions.


### 4.4. Klinkenberg Effect

Klinkenberg effect has critical impact on gas transport, deserving a deeper exploration. As shown in Figure10(a), different initial pressures in and out of the particle (pin‐ini and pout‐ini) and the same differential pressure (3 MPa) were set for the gas transport simulation in shale matrix with and without Klinkenberg effect (b=2.3MPa and 0 MPa). This kind of setting was aimed at comparing effectively at different pressure conditions. For different pin‐ini and pout‐ini, without Klinkenberg effect, the pressure drops are not the same. The higher the initial pressure, the faster the pressure drop. It is caused by the properties’ variations with pressure, which was discussed in Section 4.3. Considering the Klinkenberg effect, the pressure drop of each initial pressure was accelerated. Thus, it could be concluded that, at a certain pressure condition, the variation of real gas property and the Klinkenberg effect are both beneficial to the gas transport in shale matrix.Figure 10
The impact of Klinkenberg effect on gas transport in shale matrix at different pressure conditions. Differentpin‐ini and pout‐ini and same differential pressure with and without Klinkenberg effect (a). pin‐ini=0 and different pout‐ini, i.e., (b) 3.01 MPa, (c) 5.003 MPa, and (d) 9.184 MPa, with 1 MPa, 2 MPa, and 3 MPa Klinkenberg coefficients.
(a)(b)(c)(d)In Figure10(a), the higher the pressure is, the less difference of pressure drop curves with and without Klinkenberg effect. For further evaluating the influence of Klinkenberg effect on gas transport at a certain pressure condition, a series of simulations were conducted with different pout‐ini and Klinkenberg coefficients b (Figures 10(b)–10(d)). The initial pressures in the particle were zero. Same as Figure 10(a), the larger the Klinkenberg coefficient is, the faster the pressure drop would be. Based on Equation (7), the apparent permeability possesses positive relationship with Klinkenberg coefficient. Stronger Klinkenberg effect could bring higher apparent permeability, improving the gas transport in shale matrix. But, with higher pout‐ini (Figures 10(b)–10(d)), the discrepancy of pressure drop curves with b=1, 2, and 3 MPa tends to be less obvious. At 0 MPa pin‐ini and 9.184 MPa pout‐ini condition (Figure 10(d)), there is little difference for different Klinkenberg coefficients. It indicates that the influence of Klinkenberg effect on gas transport would be weakened with pressure growing. The Klinkenberg effect was believed to have important impact on the gas transport under the threshold pressure (5 MPa) [37]. However, for the shale sample in our research, the Klinkenberg effect still obviously influences the gas transport in matrix with 5 MPa pin‐ini and 8 MPa pout‐ini (Figure 10(a)). The Klinkenberg effect could be evaluated by Klinkenberg coefficient. Investigated by Yang et al. [52], the Klinkenberg coefficient of shale core for helium is among 0.43~0.95 MPa, which is much lower than the value obtained for the shale matrix in this research (1.5~2.3 MPa). This could be explained by the negative relationship of Klinkenberg coefficient and intrinsic permeability, found by Yang et al. [52]. The intrinsic permeability of shale matrix is much lower than shale core with fractures, due to the domination of fracture on shale core’s permeability [16]. Thus, the Klinkenberg coefficient of shale matrix would be larger than the shale core. As a result, the Klinkenberg effect in shale matrix maybe stronger than the bulk shale sample with fractures.


## 4.1. Results

### 4.1.1. Stability Analysis of Numerical Model

To ensure the stability of numerical model with explicit difference, the following condition should be established for the reasonable time and space intervals [23].
(16)ΔtΔr2≤ϕμCg2k.The permeability should be known in advance to take an effective judgment by Equation (16). Based on the analytical method [11], an estimated permeability could be obtained and help to evaluate the stability of numerical model with explicit difference. Herein, the experimental data of Sample 3 at 8.117 MPa initial pressure and 303.15 K was interpreted by analytical method. As shown in Figure 4, the experimental data was analyzed by late-time (Figure 4(a)) and early-time (Figure 4(b)) methods. It is obvious to find out that the trend lines could not fit well with the experimental data. Based on the slopes of trend lines, the permeabilities from early-time and late-time methods are 0.25×10−23m2 and 1.96×10−23m2 (Table 2). As said by Cui et al., there is great difference between the permeabilities from late-time and early-time methods [11]. It should be noted that, despite the unwell fitting of trend lines, the order of magnitude for calculated permeability is generally the same for different slopes. The unwell fitting maybe due to the ultralow permeability.Figure 4
Curve fitting with (a) late-time and (b) early-time analytical methods at 8.117 MPa and 303.15 K.FR is the mass fraction out of the particle, and FU is the cumulative gas uptake [11].
(a)(b)Table 2
Stability analysis of numerical model based on Sample 3 at 8.177 MPa and 303.15 K conditions.
μ (10-5 Pa·s)Cg (1/MPa)ϕΔt (s)mΔt/Δr2 (108 s/m2)k (early-time)/(10-23 m2)k (late-time)/(10-23 m2)ϕμCg/2k (early-time)/(s/m2)ϕμCg/2k (late-time)/(s/m2)2.030.117890.0380.011007.610.251.961.82×10102.32×109In Equation (16), the values of μ and Cg were calculated at initial pressure condition. The value of porosity of Sample 3 could be calculated by Equation (3). All the values were given in Table 2. By comparison, shown in Table 2, with 0.01 s time interval and 100 grid points, for the permeabilities from late-time and early-time methods, the value in the left of Equation (16) is much smaller than the values in the right. Therefore, 0.01 s time interval and 100 grid points could ensure the stability of numerical model.


### 4.1.2. Regression Results

As shown in Figure5, the experimental data of 3 and 4 MPa was utilized for the determinations of intrinsic permeability and Klinkenberg coefficient, and the experimental data of 5 MPa was for validation. It should be noted that, in each experiment, it is impossible to set initial value of the pressure drop to be exact 3, 4, or 5 MPa, including the other pressure conditions, owing to the drastic pressure change after the opening of intermediate valve. After regression, the simulation curves fit well with the experimental data.Figure 5
Curve fitting with numerical method at 3 and 4 MPa for Samples (a) 1, (b) 2, and (c) 3 and (d) validation of three samples at 5 MPa. The temperature is 303.15 K.
(a)(b)(c)(d)By regression, the values of porosity, intrinsic permeability, and Klinkenberg coefficient are given in Table3. The porosity was firstly calculated by Equation (4) and adjusted minutely during the curve fitting. The intrinsic permeabilities obtained in this research are close to the values from the original GRI method [22]. The intrinsic permeabilities of Wufeng-Longmaxi shale matrix are appropriate for the GRI method with 20/35-mesh particle size [10]. Moreover, the intrinsic permeabilities obtained by numerical method (Table 3) and analytical method (Table 2) are in the same order of magnitude. The intrinsic permeabilities of Wufeng-Longmaxi shale matrix are much lower than the permeability of shale sample with fractures [43], due to that the permeability of bulk shale was believed to be dominated by fractures [11].Table 3
Regressed parameters for Wufeng-Longmaxi samples.
Sampleϕkint (10-23 m2)b (MPa)10.033.42.320.0211.741.530.0383.12


## 4.1.1. Stability Analysis of Numerical Model

To ensure the stability of numerical model with explicit difference, the following condition should be established for the reasonable time and space intervals [23].
(16)ΔtΔr2≤ϕμCg2k.The permeability should be known in advance to take an effective judgment by Equation (16). Based on the analytical method [11], an estimated permeability could be obtained and help to evaluate the stability of numerical model with explicit difference. Herein, the experimental data of Sample 3 at 8.117 MPa initial pressure and 303.15 K was interpreted by analytical method. As shown in Figure 4, the experimental data was analyzed by late-time (Figure 4(a)) and early-time (Figure 4(b)) methods. It is obvious to find out that the trend lines could not fit well with the experimental data. Based on the slopes of trend lines, the permeabilities from early-time and late-time methods are 0.25×10−23m2 and 1.96×10−23m2 (Table 2). As said by Cui et al., there is great difference between the permeabilities from late-time and early-time methods [11]. It should be noted that, despite the unwell fitting of trend lines, the order of magnitude for calculated permeability is generally the same for different slopes. The unwell fitting maybe due to the ultralow permeability.Figure 4
Curve fitting with (a) late-time and (b) early-time analytical methods at 8.117 MPa and 303.15 K.FR is the mass fraction out of the particle, and FU is the cumulative gas uptake [11].
(a)(b)Table 2
Stability analysis of numerical model based on Sample 3 at 8.177 MPa and 303.15 K conditions.
μ (10-5 Pa·s)Cg (1/MPa)ϕΔt (s)mΔt/Δr2 (108 s/m2)k (early-time)/(10-23 m2)k (late-time)/(10-23 m2)ϕμCg/2k (early-time)/(s/m2)ϕμCg/2k (late-time)/(s/m2)2.030.117890.0380.011007.610.251.961.82×10102.32×109In Equation (16), the values of μ and Cg were calculated at initial pressure condition. The value of porosity of Sample 3 could be calculated by Equation (3). All the values were given in Table 2. By comparison, shown in Table 2, with 0.01 s time interval and 100 grid points, for the permeabilities from late-time and early-time methods, the value in the left of Equation (16) is much smaller than the values in the right. Therefore, 0.01 s time interval and 100 grid points could ensure the stability of numerical model.


## 4.1.2. Regression Results

As shown in Figure5, the experimental data of 3 and 4 MPa was utilized for the determinations of intrinsic permeability and Klinkenberg coefficient, and the experimental data of 5 MPa was for validation. It should be noted that, in each experiment, it is impossible to set initial value of the pressure drop to be exact 3, 4, or 5 MPa, including the other pressure conditions, owing to the drastic pressure change after the opening of intermediate valve. After regression, the simulation curves fit well with the experimental data.Figure 5
Curve fitting with numerical method at 3 and 4 MPa for Samples (a) 1, (b) 2, and (c) 3 and (d) validation of three samples at 5 MPa. The temperature is 303.15 K.
(a)(b)(c)(d)By regression, the values of porosity, intrinsic permeability, and Klinkenberg coefficient are given in Table3. The porosity was firstly calculated by Equation (4) and adjusted minutely during the curve fitting. The intrinsic permeabilities obtained in this research are close to the values from the original GRI method [22]. The intrinsic permeabilities of Wufeng-Longmaxi shale matrix are appropriate for the GRI method with 20/35-mesh particle size [10]. Moreover, the intrinsic permeabilities obtained by numerical method (Table 3) and analytical method (Table 2) are in the same order of magnitude. The intrinsic permeabilities of Wufeng-Longmaxi shale matrix are much lower than the permeability of shale sample with fractures [43], due to that the permeability of bulk shale was believed to be dominated by fractures [11].Table 3
Regressed parameters for Wufeng-Longmaxi samples.
Sampleϕkint (10-23 m2)b (MPa)10.033.42.320.0211.741.530.0383.12


## 4.2. The Effect of AGSP System

With and without the AGSP system, there is great difference for the pressure drop curves. Keep Sample 3 in SC and set the pressure in RC to be 9 MPa before the opening of intermediate valve. As shown in Figure6(a), as soon as opening the intermediate valve, the pressure would drop instantly. The pressure would take 9 seconds (Δt2 in Figure 6(a)) to balance before drop. In the 9 seconds, a certain amount of gas has entered the vacuumized particle, which affected the accuracy of pressure drop data. On the contrary, with the AGSP system, the gas was supplemented, and pressure came back to 9 MPa in 2 seconds (Δt1 in Figure 6(a)), which ensured the accuracy of pressure drop data. Besides, if the gas supplement is fast enough, the time interval, Δt1, is easy to be undetected.Figure 6
The (a) initial stage and the (b) whole process of pressure drops with and without AGSP system at 303.15 K, with Sample 3.
(a)(b)The pressure drops with and without the AGSP system could be compared effectively by the dimensionless pressure,pD (Figure 6(b)). The initial pD without the AGSP system is much smaller than the one with the AGSP system. This could be explained by the fact that, without AGSP system, large amount of gas has entered the particles before the pressure’s drop. Therefore, the regressed parameters without AGSP would appear large errors. Actually, the balance time without AGSP system has connection with the volume ratio of RC and SC. Small VRC and large VSC could bring large pressure drop, due to the larger capacity to samples, just as the treatment in the previous research [23]. However, this approach could extend the time gas expanding from RC to SC and larger Δt2 would be obtained. With the AGSP system, there is no requirement to the ratio of VRC and VSC, and the large pressure drop could be still realized.The filling degree of samples could affect the dead volume of SC, which determines the gas expanding time. Low filling degree may result in a long expanding time from RC to SC, influencing accuracy of the following pressure drop caused by gas entering the particle. The filling degree should be adjusted to minimize the void volume [20]. At least 75% filling degree was recommended for the minimization of void volume [16]. However, in our experiment, the filling degree of three 20~30-mesh samples is difficult to be larger than 50% and is only 40% around (Table 1). Thus, the low filling degree in our experiment led to a long expanding time from RC to SC, without the help of the AGSP system (Figure 6(a)). The effect of filling degree on the experiment results should be reevaluated with the help of the AGSP system. For 30.83% and 20.66% filling degrees, the experimental and simulative results are consistent with each other (Figure 7). The simulation utilized the regressed parameters of Sample 3 in Table 2. In the low filling degree situation, the gas still expands rapidly to the SC due to the high-pressure gas supplement conducted by the AGSP system. Then, the following pressure drop would be monitored accurately, without the interference of low filling degree.Figure 7
The experimental and simulative results for Sample 3 with 30.83% and 20.66% filling degrees, at 5 MPa and 303.15 K.


## 4.3. Real and Ideal Gases

For real gas, the properties, i.e.,Z, μ, and Cg, vary differently with ideal gas, resulting in the transport difference of real and ideal gas. The parameters of ideal gas in previous research [23] were applied here for a comparison with real gas transport. Same as the previous research, the simulations of real and ideal gases were conducted at 5 MPa, shown in Figure 8. The pressure of ideal gas is lower than that of the real gas at the same time, indicating that the ideal gas transports faster in the shale matrix than real gas. The distinction of real and ideal gas transports has been found in the study of ideal circular nanotubes [31, 33]. For the gas transport in the real shale sample, there is obvious and nonnegligible distinction for real and ideal gas transports in shale matrix, too.Figure 8
The simulation of real and ideal gases and the real gas simulations with ideal parameters, i.e.,Z, μ, and Cg, for Sample 3, at 5.003 MPa and 303.15 K, with Sample 3. The time interval is 50~250 s.In further comparative simulation, one of the real gas parameters,Z, μ or Cg, was replaced by the corresponding ideal parameter, as shown in Figure 8. The ideal parameters were from the previous research [23] and presented in Figure 8. In addition, at the center of particle, based on the variation of real gas pressure with time (Figure 9(a)), the corresponding variations of real and ideal compressibility factors (Figure 9(a)), real deviation factor, and real viscosity (Figure 9(b)) were provided. From Figure 9(a), the real Cg is smaller than the ideal Cg and the difference grows with pressure. In the previous research [23], the ideal Z and μ were constant. But for real gas, Z and μ both increase with pressure and have the similar tendencies with pressure (Figure 9(b)). As a result, the real Z and μ are larger than ideal Z and μ. Comparing the real gas pressure drop curves with and without the ideal parameters (Figure 8), it could be concluded that smaller Cg, Z, and μ could improve the gas transport ability in shale matrix.Figure 9
The variation of real gas pressure with time (a), and the corresponding variations of real and ideal compressibility factors (a), real deviation factor, and real viscosity (b) at 5.003 MPa and 303.15 K, with Sample 3.
(a)(b)With the pressure at the center growing rapidly after 100 s, the real and ideal compressibility factors both decrease quickly, and more obvious deviation appears (Figure9(a)). However, the simulation of real gas with ideal Cg does not deviate a lot from the simulation of real gas (Figure 8). The deviation of real and ideal compressibility factors is not large enough to lead to the great difference of real and ideal gas transports. From Figure 8, the curve of real gas with ideal μ is close to the ideal gas curve firstly and close to the real gas curve later. This is due to the fact that, with the increasing of real μ, the disparity of ideal and real μ is also growing. Conversely, the curve of real gas with ideal Z is close to the curve of ideal gas firstly and close to the curve of real gas later. It agrees with the theory that the interaction of gas molecules would be improved at high pressure and the deviation factor should not be neglected [31, 51]. It indicates that the viscosity and deviation factor, μ and Z, dominate the distinction of real and ideal gas transports at relatively low and high pressure, separately. Therefore, the real gas transport varies with ideal gas transport mainly owing to the effect of Z and μ at different pressure conditions.


## 4.4. Klinkenberg Effect

Klinkenberg effect has critical impact on gas transport, deserving a deeper exploration. As shown in Figure10(a), different initial pressures in and out of the particle (pin‐ini and pout‐ini) and the same differential pressure (3 MPa) were set for the gas transport simulation in shale matrix with and without Klinkenberg effect (b=2.3MPa and 0 MPa). This kind of setting was aimed at comparing effectively at different pressure conditions. For different pin‐ini and pout‐ini, without Klinkenberg effect, the pressure drops are not the same. The higher the initial pressure, the faster the pressure drop. It is caused by the properties’ variations with pressure, which was discussed in Section 4.3. Considering the Klinkenberg effect, the pressure drop of each initial pressure was accelerated. Thus, it could be concluded that, at a certain pressure condition, the variation of real gas property and the Klinkenberg effect are both beneficial to the gas transport in shale matrix.Figure 10
The impact of Klinkenberg effect on gas transport in shale matrix at different pressure conditions. Differentpin‐ini and pout‐ini and same differential pressure with and without Klinkenberg effect (a). pin‐ini=0 and different pout‐ini, i.e., (b) 3.01 MPa, (c) 5.003 MPa, and (d) 9.184 MPa, with 1 MPa, 2 MPa, and 3 MPa Klinkenberg coefficients.
(a)(b)(c)(d)In Figure10(a), the higher the pressure is, the less difference of pressure drop curves with and without Klinkenberg effect. For further evaluating the influence of Klinkenberg effect on gas transport at a certain pressure condition, a series of simulations were conducted with different pout‐ini and Klinkenberg coefficients b (Figures 10(b)–10(d)). The initial pressures in the particle were zero. Same as Figure 10(a), the larger the Klinkenberg coefficient is, the faster the pressure drop would be. Based on Equation (7), the apparent permeability possesses positive relationship with Klinkenberg coefficient. Stronger Klinkenberg effect could bring higher apparent permeability, improving the gas transport in shale matrix. But, with higher pout‐ini (Figures 10(b)–10(d)), the discrepancy of pressure drop curves with b=1, 2, and 3 MPa tends to be less obvious. At 0 MPa pin‐ini and 9.184 MPa pout‐ini condition (Figure 10(d)), there is little difference for different Klinkenberg coefficients. It indicates that the influence of Klinkenberg effect on gas transport would be weakened with pressure growing. The Klinkenberg effect was believed to have important impact on the gas transport under the threshold pressure (5 MPa) [37]. However, for the shale sample in our research, the Klinkenberg effect still obviously influences the gas transport in matrix with 5 MPa pin‐ini and 8 MPa pout‐ini (Figure 10(a)). The Klinkenberg effect could be evaluated by Klinkenberg coefficient. Investigated by Yang et al. [52], the Klinkenberg coefficient of shale core for helium is among 0.43~0.95 MPa, which is much lower than the value obtained for the shale matrix in this research (1.5~2.3 MPa). This could be explained by the negative relationship of Klinkenberg coefficient and intrinsic permeability, found by Yang et al. [52]. The intrinsic permeability of shale matrix is much lower than shale core with fractures, due to the domination of fracture on shale core’s permeability [16]. Thus, the Klinkenberg coefficient of shale matrix would be larger than the shale core. As a result, the Klinkenberg effect in shale matrix maybe stronger than the bulk shale sample with fractures.


## 5. Conclusions

In this research, based on GRI method, the experimental apparatus and the data interpretation model were modified for an accurate and deep study on gas transport with real gas effect and Klinkenberg effect in Wufeng-Longmaxi shale matrix. The detailed conclusions of this research are shown as below.(1)
An AGSP system was added to the original apparatus for GRI method. The AGSP system helps to carry out experiment at high pressures. By experimental evaluation of AGSP system, the gas expands from RC to SC much faster than the original apparatus, and the negative effect of low filling degree on gas expansion would be eliminated. Accordingly, the reliability of the obtained experimental data could be ensured(2)
Considering the real gas effect and Klinkenberg effect, an explicit numerical model was established with variable gas properties and apparent permeability. The stability of numerical model was also conducted. By experimental data interpretation with the model, the intrinsic permeability and Klinkenberg coefficient were obtained. The intrinsic permeability of Wufeng-Longmaxi shale matrix is much smaller than the shale core plug with developed fractures(3)
There is obvious discrepancy for real and ideal gas transports in shale matrix. SmallerCg, Z, and μ of gas could improve the transport ability. The ideal gas transports faster than real gas owing to the dominance of smaller μ at low pressure condition and the dominance of smaller Z at high pressure condition(4)
At a certain pressure condition, the variable real gas property and the Klinkenberg effect are both beneficial to the gas transport in shale matrix. The higher the pressure, the less impact of Klinkenberg effect on transport would be. The impact of Klinkenberg effect is still obvious at 5 MPa and could be ignored at around 9 MPa



---

*Source: 5579307-2021-03-19.xml*