Articles | Volume 27, issue 24
Research article
20 Dec 2023
Research article |  | 20 Dec 2023

A pulse-decay method for low (matrix) permeability analyses of granular rock media

Tao Zhang, Qinhong Hu, Behzad Ghanbarian, Derek Elsworth, and Zhiming Lu

Nanodarcy level permeability measurements of porous media, such as nano-porous mudrocks, are frequently conducted with gas invasion methods into granular-sized samples with short diffusion lengths and thereby reduced experimental duration; however, these methods lack rigorous solutions and standardized experimental procedures. For the first time, we resolve this by providing an integrated technique (termed gas permeability technique, GPT) with coupled theoretical development, experimental procedures, and data interpretation workflow. Three exact mathematical solutions for transient and slightly compressible spherical flow, along with their asymptotic solutions, are developed for early- and late-time responses. Critically, one late-time solution is for an ultra-small gas-invadable volume, important for a wide range of practical usages. Developed to be applicable to different sample characteristics (permeability, porosity, and mass) in relation to the storage capacity of experimental systems, these three solutions are evaluated from essential considerations of error difference between exact and approximate solutions, optimal experimental conditions, and experimental demonstration of mudrocks and molecular-sieve samples. Moreover, a practical workflow of solution selection and data reduction to determine permeability is presented by considering samples with different permeability and porosity under various granular sizes. Overall, this work establishes a rigorous, theory-based, rapid, and versatile gas permeability measurement technique for tight media at sub-nanodarcy levels.

Please read the corrigendum first before continuing.

  • An integrated (both theory and experiments) gas permeability technique (GPT) is presented.

  • Exact and approximate solutions for three cases are developed with error discussion.

  • Conditions of each mathematical solution are highlighted for critical parameters.

  • Essential experimental methodologies and data processing procedures are provided and evaluated.

1 Introduction

Shales, crystalline, and salt rocks with low permeabilities (e.g., <10−17 m2 or approximately 10 microdarcy, µd) are critical components in numerous subsurface studies. Notable examples are the remediation of contaminated sites (Neuzil, 1986; Yang et al., 2015), long-term performance of high-level nuclear waste repositories (Kim et al., 2011; Neuzil, 2013), enhanced geothermal systems (Huenges, 2016; Zhang et al., 2021), efficient development of unconventional oil and gas resources (Hu et al., 2015; Javadpour, 2009), long-term sealing for carbon utilization and storage (Fakher et al., 2020; Khosrokhavar, 2016), and high-volume and effective gas (hydrogen) storage (Liu et al., 2015; Tarkowski, 2019). For fractured rocks, the accurate characterization of rock matrix and its permeability is also critical for evaluating the effectiveness of low-permeability media, particularly when transport is dominated by slow processes like diffusion (Ghanbarian et al., 2016; Hu et al., 2012).

Standard permeability test procedures in both steady-state and pulse-decay methods use consolidated centimeter-sized core-plug samples, which may contain fractures and show dual- or triple-porosity characteristics (Abdassah and Ershaghi, 1986; Bibby, 1981). The overall permeability may therefore be controlled by a few bedding-oriented or cross-cutting fractures, even if experiments are conducted at reservoir pressures (Bock et al., 2010; Gensterblum et al., 2015; Gutierrez et al., 2000; Luffel et al., 1993). Fractures might be naturally or artificially induced (e.g., created during sample processing), which makes a comparison of permeability results among different samples difficult (Heller et al., 2014). Hence, methods for measuring the matrix (non-fractured) permeability in tight media, with the practical necessity to use granular samples, have attracted much attention for their ability to help eliminate the sides effect of fractures (Civan et al., 2013; Egermann et al., 2005; Heller et al., 2014; Wu et al., 2020; Zhang et al., 2020).

A GRI (Gas Research Institute) method was developed by Luffel et al. (1993) and followed by Guidry et al. (1996) to measure the matrix permeability of crushed mudrocks (Guidry et al., 1996; Luffel et al., 1993). This method makes permeability measurement feasible in tight and ultra-tight rocks (with permeability <10−20 m2 or approximately 10 nd), particularly when permeability is close to the detection limit of the pulse-decay approach on core plugs at ∼10 nd (e.g., using the PoroPDP-200 commercial instrument of CoreLab). In the GRI method, helium may be used as the testing fluid to determine permeability in crushed samples at different sample sizes (e.g., within the 10–60 mesh range, which is from 0.67 to 2.03 mm). The limited mesh size of 20–35 (500–841 µm in diameter) was recommended in earlier works, which has led to the colloquial names of “the GRI method/size” in the literature (Cui et al., 2009; Kim et al., 2015; Peng and Loucks, 2016; Profice et al., 2012). However, Guidry et al. (1996) and Luffel et al. (1993) did not document the processing methodologies needed to derive the permeability from experimental data from this GRI method. That is, there are neither standard experimental procedures for interpreting gas pulse-decay data in crushed rock samples nor detailed mathematical solutions available for data processing in the literature (Kim et al., 2015; Peng and Loucks, 2016; Profice et al., 2012). In this work, we seek to (1) develop mathematical solutions to interpret gas pulse-decay data in crushed rock samples without an available published algorithm as this method shares different constitutive phenomena than the traditional pulse-decay method for core-plug samples in Cartesian coordinates and (2) present the associated experimental methodology to measure permeability, reliably, and reproducibly in tight and ultra-tight granular media.

We first derive the constitutive equations for gas transport in granular (unconsolidated or crushed rock) samples. Specifically, we develop three mathematical solutions which cover different experimental situations and sample properties. As each solution has its own pros and cons, we then present in detail the error analyses for the derived exact and approximate solutions and discuss their applicable requirements and parameter recommendation for practical usages. This work aims to fill the knowledge gap in granular rock (matrix) permeability measurement and follow up the literature by establishing an integrated methodology for reproducible measurements of nanodarcy-level permeability in tight rocks for emerging energy and resources subsurface studies.

2 Mathematical solutions for gas permeability of granular samples

For a compressible fluid under unsteady-state conditions, flow in a porous medium can be expressed by the mass conservation equation:

(1a) p t + ρ v = 0

where p is the pressure, t is the time, ρ is the fluid density, and v is the Darcy velocity. In continuity equations derived for gas flow in porous media, permeability can be treated as a function of pressure through the ideal gas law. Constitutive equations are commonly established for a small pressure variation to avoid the non-linearity of gas (the liquid density to be a constant) and to ensure that pressure would be the only unknown parameter (Haskett et al., 1988). For spherical coordinates of fluid flow in porous media, assuming flow along the radial direction of each spherical solid grain, Eq. (1a) becomes

(1b) p t ϕ = 1 c t k μ r 2 r r 2 p r .

The gas compressibility ct is given by

(1c) c t = 1 ρ d ρ d p = 1 p - 1 z d z d p .

In Eqs. (1b) and (1c), ϕ and k are sample porosity and permeability, r is the migration distance of fluid, μ is the fluid viscosity, and z is the gas deviation (compressibility) factor and is constant.

To correct for the non-ideality of the probing gas, we treat gas density as a function of pressure and establish a relationship between the density and the permeability through a pseudo-pressure variable given in Supplement Sect. S1. Detailed steps for deriving mathematical solutions for the GPT, based on heat transfer studies (Carslaw and Jaeger, 1959), can be found in Sect. S2. The Laplace transform is an efficient tool for solving gas transport in granular samples with low permeabilities, as applied in this study. Alternatively, other approaches, such as the Fourier analysis, Sturm–Liouville method, or Volterra integral equation of the second form may be used (Carslaw and Jaeger, 1959; Haggerty and Gorelick, 1995; Ruthven, 1984).

We applied dimensional variables to derive the constitutive equation given in Eq. (S10) for which the initial and boundary conditions are

(2a) 2 U s ξ 2 + s 2 U s = 0 U s = 0 , ξ = 0 ,

(2b) α 2 U s - 1 = 3 K c U s ξ - U s ξ ξ = 1 ,

where Us and ξ represent the dimensionless values of gas density and sample scale and s is the transformed Heaviside operator. α in Eq. (2b) is determined by solving Eq. (S30) for its root. Kc in Eq. (2b) is a critical parameter that represents the volumetric ratio of the total void volume of the sample cell to the pore volume of the porous samples. It is similar to the storage capacity, controlling the acceptable measurement range of permeability and decay time, in the pulse-decay method proposed by Brace et al. (1968).

The fractional gas transfer for the internal (limited Kc value) and external (infinite Kc value) gas transfer of sample is given by

(2c) F f = 1 - 6 n = 1 K c 1 + K c e - α n 2 τ 9 K c + 1 + α n 2 K c 2 ,

(2d) F s = 1 - 6 π 2 n = 1 e - ( n π ) 2 τ n 2 ,

where Ff and Fs represent the uptake rate of gas outside and inside the sample separately as a dimensionless parameter and τ is the Fourier number of dimensionless time. Three approximate solutions of the transport coefficient based on Eqs. (2c) and (2d) for various conditions are presented below.

The late-time solution to Eq. (2c) for a limited Kc value (called LLT hereafter) is

(3a) k = R a 2 μ c t ϕ f s 1 α 1 2 .

The late-time solution to Eq. (2d) when Kc tends to infinity (ILT hereafter) is

(3b) k = R a 2 μ c t ϕ f s 2 π 2 .

The early-time solution to Eq. (2d) when Kc approaches infinity (IET hereafter) is

(3c) k = π R a 2 μ c t ϕ f s 3 36 .

In Eq. (3), Ra is the particle diameter of a sample and s1, s2, and s3 are the three exponents that may be determined from the slopes of data on double logarithmic plots. Table 1 summarizes Eqs. (3a) to (3c) and conditions under which such approximate solutions would be valid.

Table 1Solutions schematic with difference Kc and τ values.

a This defines the volumetric ratio of the total void volume of the sample cell to the pore volume of the porous samples; the classification between the limited and infinity value is proposed as 50 in the following analyses. b The original constitutive equation for different Kc value. c Equations (3a)–(3c) are three approximate solutions of density faction function F.

Download Print Version | Download XLSX

Based on diffusion phenomenology, Cui et al. (2009) presented two mathematical solutions similar to our Eqs. (3a) and (3c). In the work of Cui et al. (2009), however, one of late-time solution is missing, and error analyses are not provided. In addition, the lack of detailed analyses of τ and Kc in the constitutive equations will likely hinder the practical application of Eq. (3b), which is able to cover an experimental condition of small sample mass with a greater τ (further analyzed in Sect. 3). Furthermore, the early-time and late-time solution criteria are not analyzed, and the pioneering work of Cui et al. (2009) does not comprehensively assess practical applications of their two solutions in real cases, which is addressed in this study. Hereafter, we refer to the developed mathematical and experimental gas-permeability-measurement approach holistically as the gas permeability technique (GPT).

3 Practical usages of algorithms for the GPT

As aforementioned, mathematical solutions given in Eqs. (3a) and (3b) were deduced based on different values of Kc and τ, as shown in Sect. S2. This means each solution holds only under specific experimental conditions, which are mostly determined by the permeability, porosity, and mass of samples, as well as gas pressure and void volume of the sample cell. In this section, the influence of parameters Kc and τ on the solution of constitutive equation is analyzed and a specific value of dimensionless time (τ=0.024) is proposed as the criterion required to detect the early-time regime from the late-time one for the first time in the literature. We also demonstrate that the early-time solution of Eq. (3c), which has been less considered for practical applications in previous studies, is also suitable and unique under common situations. In addition, the error in the approximate solution compared to the exact solution and their capabilities are discussed, as it helps to select an appropriate mathematical solution at small τ values. Moreover, we showcase the unique applicability and feasibility of the new solution of Eq. (3b).

3.1 Sensitivity analyses of the Kc value for data quality control

To apply the GPT method, appropriately selecting the parameter Kc in Eqs. (3a)–(3c) is crucial, as it is a critical value for data quality control. The dimensionless density outside the sample, Uf, is related to Kc via Eq. (S33) in Sect. S2. One may simplify Eq. (S33) by replacing the series term with some finite positive value and set

(4) U f - K c 1 + K c > 0 .

We define Kf=Kc/(1+Kc) to interpret the density variance of the system as Kf is closely related to the dimensionless density outside the sample, Uf.

Equation (4) shows the relationship between Uf and Kc (Fig. 1). For Kc>0, Kf falls between 0 and 1. The greater the Kf value, the more insensitive to density changes the system is. For Kc equal to 50, Kf would no longer be sensitive to Kc variations as it would have already approached 98 % of the dimensionless density. This means that the Uf value needs to be greater than 0.98, and this leaves only 2 % of the fractional value of Uf available for capturing gas density change. When Kc is 100, the left fractional value of Uf would be 1 %. This would limit the number of data available (the linear range in Fig. S1) for the permeability calculation, which would complicate the data processing. Thus, for the GPT experiments, a small value of Kc (less than 10) is recommended, as Kf nearly reaches its plateau beyond Kc=10 (Fig. 1). When Kc is 10, the left fractional value of Uf is only as low as 9 %.

Figure 1Dimensionless density, Kf, as a function of dimensionless volume Kc. Major variations in Kf occur for Kc<10, indicating longer gas transmission duration with more pressure decay data available for permeability derivation.


3.2 Recommendation for solution selection

The following three aspects need to be considered before selecting the appropriate solution for permeability calculation: (1) early- or late-time solutions, (2) error between the approximate and exact solutions, and (3) the convenience and applicability of solutions suitable for different experiments. We will first discuss the selection criteria for early- or late-time solutions.

Figure 2a shows the exact solution of Fs with their approximate early- and late-time solution (Table 1). Two exact solutions of Ff where Kc equals 10 or 50 are also demonstrated in Fig. 2a. Figure 2b depicts the exact solution from Ff for different Kc values from 1 to 100 and their corresponding approximate solution for Eq. (3a). The intersection points of the solution Eqs. (3b) and (3c), namely τ=0.024 in Fig. 2a, is used for distinguishing early- and late-time solutions.

Two notable observations can be drawn from Fig. 2b. First, the approximate solution Eq. (3a) would only be applicable at late times when τ is longer than 0.024. For τ< 0.024, regardless of the Kc value, Eq. (3c) would be more precise than Eqs. (3a) and (3b) and return results close to the exact solution for both Ff and Fs. Second, the results of Eqs. (3a) and (3b) presented in Fig. 2a are similar; their difference is very small, especially for Kc>10. Due to the fact that core samples from deep wells are relatively short in length and their void volume is small (ultra-low porosity and permeability such as in mudrocks with k≤ 0.1 nd), in practice, a solution for 10<Kc<100 is the most common outcome, even if the sample cell is loaded as full as possible. Under such circumstances, the newly derived solution, Eq. (3b), becomes practical and convenient: (1) if the Kc and dimensionless time τ have not been evaluated precisely before the GPT experiment, this solution may fit most experimental situations and (2) this solution is suitable for calculation as it does not need the solution from the transcendental equation of Eq. (S30) because the denominator of α has been replaced by π. The data quality control is discussed in Sect. 4.1.

Figure 2Three GPT solutions with different values of τ, Kc; the dashed lines are approximate solutions without a series expansion in Fig. 2b for Ff. Figure modified from Cui et al. (2009).

3.3 Applicability of the early-time solution

A small Kc value can guarantee a sufficient time for gas transfer in samples and provide enough linear data for fitting purposes. We note that the selection of the limited Kc solution of Ff, and the infinity Kc solution Fs is controlled by Kc. However, before the selection of Kc, the dimensionless time is the basic parameter to be estimated as a priori before the early- or late-time solutions are selected.

For pulse-decay methods, the early-time solution has the advantage of capturing the anisotropic information contained in reservoir rocks (Jia et al., 2019; Kamath, 1992). However, it suffers from the shortcoming of uncertainty in data for the initial several seconds, which is why it is not recommended for data processing (Brace et al., 1968; Cui et al., 2009). This is due to (1) the Joule–Thomson effect, which causes a decrease in gas temperature from the expansion, (2) kinetic energy loss during adiabatic expansion, and (3) collision between molecules and the container wall. These uncertainties normally occur in the first 10–30 s, shown in our experiments as a fluctuating period called the “early stage”.

However, the early stage present in pulse-decay experiments does not mean that the early-time solution is not applicable. We demonstrate the relationship between time and dimensionless time in Fig. 3, showing that a short dimensionless time may correspond to a long testing period of hundreds to thousands of seconds in experiments. This is particularly notable for the ultra-low-permeability samples with k≤0.1 nd and small dimensionless times τ<0.024. This situation would only be applicable to early-time solutions but with data available beyond the early stage, and it would provide available data over a long time (hundreds to thousands of seconds). For example, the early-time solution would fit ultra-low-permeability samples in 600 s for 0.1 nd and in at least 1000 s for 0.01 nd, as shown in Fig. 3 in the region below the dark line. Then, using Eq. (3c), the derived permeability would be closer to its exact solution in the earlier testing time (but still after the early stage). The mudrock samples that we tested, with results presented in Sect. 5.3, exhibit low permeabilities, approximately on the order of 0.1 nd.

Figure 3Dimensionless time τ versus actual time for different permeability values through Eq. (S14) using He gas, sample porosity of 5 %, and sample diameter of 2 mm, where 1 nd is nearly equal to 10−21 m2.


3.4 Error analyses between exact and approximate solutions

It is unpractical to use the exact solutions with their series part to do the permeability calculation; thus, only the approximate solutions are used and the error difference between the exact and approximate solutions is discussed here. The original mathematical solutions, Eqs. (S39) and (S49), are based on series expansion. For dimensionless densities Ff and Fs in Eqs. (S39) and (S49), their series expansion terms should converge. However, the rate of convergence is closely related to the value of τ. For example, from Eq. (S30), when τ≥1, the exponent parts of Us and Uf are at least (2n+1)π2. Therefore, the entire series expansion term can be omitted without being influenced by Kc. In practical applications, the solutions given in Eqs. (3a)–(3c) are approximates without series expansion. In this study, we provide diagrams of the change in errors with dimensionless time in the presence of adsorption (Fig. 4).

For Ff, the error differences between the exact and approximate solutions are 3.5 % and 0.37 % for τ=0.05 and 0.1 when Kc=10, respectively. When τ≤0.024, the error would be greater than 14.7 %. Figure 2(b) shows that Ff can be approximated as Fs when Kc is greater than 10; the error difference between Ff and Fs is quite small at this Kc value (for Kc=10, 6.6 % is the maximum error when τ=0.01; 4.4 % when τ=0.05; and 2.9 % when τ=0.1), as shown in Fig. 4.

For Fs, the error difference is roughly the same as Ff and equal to 3.6 % for τ=0.05 and 0.38 % for τ=0.1. This verifies that newly derived Eq. (3b) is equivalent to Eq. (3a) when Kc is greater than 10. As for the evaluation of Eq. (3c), the error difference with the exact solution will increase with dimensionless time (5.1 % for τ=0.003, 9.7 % for τ=0.01, and 16 % for τ=0.024).

Figure 4Error analyses of Ff and Fs for their exact and approximate solutions.


4 Influence of kinetic energy on gas transport behavior

4.1 Flow state of gas in granular samples

In the following, we apply the approximate solutions, Eqs. (3a)–(3c), to some detailed experimental data and determine permeability in several mudrock samples practically compatible with sample size, gases, and molecular dynamics analyses.

During the GPT, with the boundary conditions described in Sect. S2, the pressure variation is captured after gas starts to permeate the sample from the edge, and the model does not take into account the gas transport between particles or into any micro-fractures if present. Thus, the transport that conforms to the “unipore” model and occurs after the early stage (defined in Sect. 3.3) or in the “penetration zone” (the area between the two vertical lines in Fig. 5) should be used to determine the slope. Figure A2 shows how to obtain the permeability result using the applicable mathematical solutions (Eqs. 3a–3c). Figure 5 shows the pressure variance with time during the experiment using sample size from 0.34 to 5.18 mm for sample X-1 and sample X-2. From Fig. 5, the time needed to reach pressure equilibrium after the initial fluctuation stage is 20–100 s, and the penetration zone decreases with decreasing grain size over this time period.

Figure 5Fitting region (the penetration zone in the shadowed area) for mudrock sample X-1 with different granular sizes; the penetration zone illustrating the pressure gradient mainly happens at 20 to 200 s for this sample.


In fact, the penetration zone, as an empirical period, is evaluated by the pressure change over a unit of time before gas is completely transported into the inner central part of the sample to reach the final pressure. Owing to the sample size limitation, a decreasing pressure could cause multiple flow states (based on the Knudsen number) to exist in the experiment. The pressure during the GPT experiment varies between 50 and 200 psi (0.345 to 1.38 MPa). Figure 6 shows the Knudsen number calculated from different pressure conditions and pore diameters together with their potential flow state. Based on Fig. 6, the flow state of gas in the GPT experiments is mainly dominated by Fickian and transition diffusion. Essentially, the flow state change with pressure should be strictly evaluated through the Knudsen number in Fig. 6 to guarantee that the data in the penetration zone are always fitted with the GPT's constitutive equation for laminar-flow or diffusive states. This helps obtain a linear trend for ln (1−Ff) or Fs2 versus time for low-permeability media. Experimentally, data from 30 to several hundreds of seconds are recommended for tight rocks like shales within the GPT methodology.

Figure 6Flow state of gas under different testing pressures, with 145 psi being approximately equal to 1 MPa; modified from Chen and Pfender, 1983, and Roy et al., 2003.

In the GPT approach, as mentioned earlier, Eq. (S33) holds for small Kc values (e.g., <10) so that the approximately equivalent void volume in the sample cell and sample pore volume would allow for sufficient pressure drop. It also gives time and allows the probing gas to expand into the matrix pores to have a valid penetration zone and to determine the permeability. Greater values of Kc would prevent the gas flow from entering into a slippage state as the pressure difference would increase with increasing Kc. However, large pressure changes would result in a turbulent flow (Fig. 6), which would cause the flow state of gas to be no longer valid for the constitutive equation of the GPT. Overall, the GPT solutions would be applicable to gas permeability measurement, based on the diffusion-like process, from laminar flow to Fickian diffusion, after the correction of the slippage effect.

4.2 Pressure decay behavior of four different probing gases

We used three inert gases, He, N2, and Ar, and one sorptive gas, CO2 (Busch et al., 2008), to compare the pressure drop behavior for the sample size with an average granular diameter of 0.675 mm. Results for the mudrock sample X-2 are presented in Fig. 7. Among the three inert gases, helium and argon required the shortest and longest time to reach pressure equilibrium (i.e., He <N2<Ar). In terms of pressure drop, argon exhibited the most significant decrease. In a constant-temperature system, the speed (or rate) at which gas molecules move is inversely proportional to the square root of their molar masses. Hence, it is reasonable that helium (with the smallest kinetic diameter of 0.21 nm) has the shortest equilibrium time. However, the pressure drop is more critical than the time needed to reach equilibrium for the GPT, as the equilibrium time does not differ much (basically within 10 s for a given sample weight, except for the adsorptive CO2). Argon may provide a wider range of valid penetration zones given its longest decay time (when excluding the adsorbed gas CO2).

Figure 7Measured pressure decay curves from mudrock sample X-2 for gases of different molecular diameters σ and molecular weights M (g mol−1); decay happened at around 50 psi (approximately 0.345 MPa)


Figure 7 shows that the pressure decay curve of the adsorptive gas CO2 is different from those of the inert gases used in this study. CO2 has a slow equilibrium process due to its large molar mass and the greatest pressure drop among the four gases due to its adsorption effect. This additional flux needs to be taken into account to obtain an accurate transport coefficient. Accordingly, multiple studies including laboratory experiments (Pini, 2014) and long-term field observations (Haszeldine et al., 2006; Lu et al., 2009) were carried out to assess the sealing efficiency of mudrocks for CO2 storage. In fact, the GPT offers a quick and effective way to identify the adsorption behavior of different mudrocks for both laminar-flow and diffusion states.

4.3 Pressure decay behavior for different granular sizes

We compared the pressure drop behavior of gas in the mudrock sample X-1 with different granular sizes (averaging from 0.34 to 5.18 mm) using the same sample weight and Kc. Results based on the experimental data shown in Fig. 8 indicate that a larger-sized sample would provide more data to be analyzed for determining the permeability. This is because the larger the granular size, (1) the larger the pressure drop and (2) the greater the decay time, as Fig. 8 demonstrates. This is consistent with the simulated results reported by Profice et al. (2012).

Figure 8Pressure decay curves (50 psi equals approximately 0.345 MPa) measured by helium on sample X-1 with five different granular sizes. The intra-granular porosity was 5.8 %, independently measured by mercury intrusion porosimetry.


Table 2Permeability results from the methods of GPT and SMP-200 for X-1.

a The results are from the SMP-200 using the GRI default method. b The results are from the GPT method we proposed. c The result which failed the criteria of dimensionless time. d CR means conflicting results, as the verified dimensionless time fails to corroborate the early- or late-time solutions derived from the calculated permeability. For example, the verified dimensionless time would be >0.024 using the early-time solution result and vice versa.

Download Print Version | Download XLSX

As reported in Table 2, the permeability values measured by the GPT method are 1 or 2 orders of magnitude greater than those measured by the SMP-200 instrument. The built-in functions of SMP-200 can only be used for two default granular sizes (500–841 µm for GRI and 1.70–2.38 mm for what we call GRI+) to manually curve-fit the pressure decay data and determine the permeability. The GRI method built in the SMP-200 only suggests the fitting procedure for data processing without publicly available details on the underlying mathematics. The intra-granular permeabilities of mudrock sample X-1 vary from 0.02 to 1.17 nd for five different granular sizes using the GPT. With the same pressure decay data selected from 30 to 200 sec, the permeability results for GRI and GRI+ sample sizes from the SMP-200 fitting are 0.65 and 14.2 nd as compared to 0.06 and 0.43 nd determined by the GPT using the same mean granular size. Our results are consistent with those reported by Peng and Loucks (2016), who found differences of 2 to 3 orders of magnitude between the GPT and SMP-200 methods (Peng and Loucks, 2016).

There are several issues associated with granular samples with diameters smaller than, on average, 1.27 mm. First, the testing duration is short, and second, there is not sufficient pressure variation analyzed in Fig. 8. Both may cause significant uncertainties in the permeability calculation and, therefore, make samples with diameters smaller than 1.27 mm unfavorable for the GPT method, particularly ultra-tight (sub-nanodarcy levels) samples as there is almost no laminar-flow or diffusion state to be captured. The greater pressure drop for larger-sized granular samples would result in greater pressure variation and wider data region compared to smaller granular sizes (see Figs. 6 and 9). Although samples of large granular sizes may potentially contain micro-fractures, which complicate the determination of true matrix permeability (Heller et al., 2014), the versatile GPT method can still provide size-dependent permeabilities for a wide range of samples (e.g., from sub-millimeter to 10 cm diameter full-size cores) (Ghanbarian, 2022a, b). In addition, the surface roughness of large grains may also complicate the determination of permeability, which requires attention (Devegowda, 2015; Rasmuson, 1985; Ruthven and Loughlin, 1971). Overall, our results demonstrated that sample diameters larger than 2 mm are recommended for the GPT to determine the nanodarcy permeability of tight mudrocks, while smaller sample sizes may produce uncertain results.

4.4 Practical recommendations for the holistic GPT

Here, we evaluate the potential approximate solution for tight rock samples using frequently applied experimental settings by considering the critical parameters, such as sample mass, porosity, and estimated permeability (as compiled in Fig. 9 showing the dimensionless time versus porosity). Based on the results presented in Figs. 3 and 6, only t<200 s is dominant and critical for the analyses of dimensionless time and penetration zone. Thus, we take 200 s and use helium to calculate the dimensionless time. Another critical parameter to ensuring sufficient decay data is a sample diameter of greater than 2 mm. Thus, we only show the dimensionless time versus porosity for sample diameters greater than the criteria of 2 mm.

Figure 9 demonstrates that the sample permeability has dominant control on early- or late-solution selection, and we decipher a concise criterion for the selection of the three solutions. We classify the dimensionless time versus porosity relationship into three cases. First, among the curves shown in Fig. 9, only that corresponding to k=0.1 nd and sample diameter of 2 mm stays below the dashed line representing τ=0.024. Therefore, the early-time solution is appropriate for tight samples with permeabilities less than 0.1 nd (as shown in the analyses of Sect. 4.3, which also conforms to the situation of the molecular sieve sample that we tested in Sect. S3). Second, for permeabilities greater than 10 nd (the curve is above the line of τ=0.024), the newly derived late-time solution, Eq. (3b), is recommended as it is more convenient for mathematical calculation without the consideration of transcendental functions. The reason is that the sample cell can be filled as much as possible (∼90 % of the volume) with samples and solid objects. However, as the tight rock hardly presents a large value of porosity, it is difficult to achieve a small Kc value, with an inconsequential influence between Eqs. (3b) and (3a). Lastly, in the case of permeability around 1 nd, the value of porosity would be critical in the selection of the early- or late-time solutions, as shown in Fig. 9.

Figure 9Holistic GPT used to explore the appropriate solution based on diameter, permeability, and porosity of samples. The legend shows the diameter of granular sample and permeability, along with a dashed line for a dimensionless time of 0.024, while regions above and below this value fit for the late- and early-time solutions, respectively (1 nd is nearly equal to 10−21 m2).


5 Conclusions

In the present work, we solved a fluid flow state equation in granular porous media and provided three exact mathematical solutions along with their approximate ones for practical applications of low-permeability measurements. The mathematical solutions for the transport coefficient in the GPT were derived for a spherical coordinate system, applicable from laminar flow to slippage-corrected Fickian diffusion. Of the three derived solutions, the early-time solution is valid when the gas storage capacity Kc approaches infinity, while the other two are late-time solutions, valid when Kc is either small or tends towards infinity. We evaluated the derived solutions for a systematic measurement of ultra-low permeabilities in granular media and crushed rocks using experimental methodologies with the data processing procedures. We determined the error for each solution by comparing with the exact solutions presented in the Supplement. The applicable conditions for such solutions of the GPT were investigated, and we provided the selection strategies for three approximate solutions based on the range of sample permeability. In addition, a detailed utilization of GPT was given to increase confidence in the GPT method through the molecular sieve sample, as it allows for a rapid permeability test for ultra-tight rock samples in just tens to hundreds of seconds, with good repeatability.

Appendix A: Nomenclature
Bij Correction parameter for viscosity, constant
ct Fluid compressibility, Pa−1
Ff Uptake rate of gas outside the sample, dimensionless
Fs Uptake rate in the sample, dimensionless
f1 Intercept of Eq. (S40), constant
Ka Apparent transport coefficient defined as Eq. (S9), m2 s−1
Kc Ratio of gas storage capacity of the total void volume of the system to the pore (including adsorptive and
non-adsorptive transport) volume of the sample, fraction
Kf Initial density state of the system, fraction
k Permeability, m2
ks Permeability defined as Eq. (S8), m2 (Pa s)−1
L Coefficient, unit for certain physical transport phenomenon
M Molar mass, kg kmol−1
Mm Molar mass of the mixed gas, kg kmol−1
Mi,j Molar mass for gas i or j, kg kmol−1
Ms Total mass of sample, kg
N Particle number, constant
p Pressure, Pa
pcm Virtual critical pressure of mixed gas, Pa
pp Pseudo-pressure from Eq. (S1), Pa s−1
Ra Particle diameter of sample, m
R Universal gas constant, 8.314 J (mol K)−1
r Diameter of sample, m
s1 Slope of Eq. (S40), constant
s2 Slope of function Ln(1−Fs), constant
s3 Slope of function Fs2, constant
T Temperature, K
Tcm Virtual critical temperature for mixed gas, K
t Time, s
Uf Dimensionless density of gas outside the sample, dimensionless
Us Dimensionless density in grain, dimensionless
U Maximum density defined as Eq. (S37), dimensionless
V1 Cell volume in upstream of pulse-decay method, m3
V2 Cell volume in downstream of pulse-decay method, m3
Vb Bulk volume of sample, m3
Vc Total system void volume except for sample bulk volume, m3
v Darcian velocity in pore volume of porous media, m s−1
X Pressure force, Pa
yi,j Molar fraction for gas i or j, fraction
z Gas deviation (compressibility) factor, constant
Appendix B: Greek letters
αn nth root of Eq. (S30), constant
μ Dynamic viscosity, Pa s or N s m−2
μi,j Dynamic viscosity for gas i or j, Pa s or N s m−2
μmix Dynamic viscosity of mixture gas, Pa s or N s m−2
μp Correction term for the viscosity with pressure, Pa s or N s m−2
ξ Dimensionless radius of sample, dimensionless
ρ Density of fluid, kg m−3
ρ0 Average gas density on the periphery of sample, kg m−3
ρ1 Gas density in reference cell, kg m−3
ρ2 Gas density in sample cell, kg m−3
ρb Average bulk density for each particle, kg m−3
ρf Density of gas changing with time outside sample, kg m−3 s−1
ρf Maximum value of ρf defined as Eq. (S38), kg m−3 s−1
ρp Pseudo-density from Eq. (S1), kg m−3 s−1
ρs Density of gas changing with time in sample, kg m−3 s−1
ρps Pseudo-density of gas changing with time in sample, kg m−3 s−1
ρpf Pseudo-density of gas changing with time outside sample, kg m−3 s−1
ρp2 Initial pseudo-density of gas in sample, kg m−3 s−1
ρp0 Average pseudo-density of gas on sample periphery, kg m−3 s−1
ρrm Relative density to the mixed gas, kg m−3 s−1
ρsav Average value of ρsr defined as Eq. (S47), kg m−3 s−1
ρsr Average value of pseudo-density of sample changing with diameter, kg m−3 s−1
ρs Maximum value of ρsr defined as Eq. (S46), kg m−3 s−1
τ Dimensionless time, dimensionless
ϕ Sample porosity, fraction
ϕf Total porosity (ϕf=ϕa+ϕb) occupied by both free and adsorptive fluids, fraction
Data availability

This work did not use any data from previously published sources, and our experimental data and MATLAB processing codes are available at (Hu, 2020), managed by Mavs Dataverse of the University of Texas at Arlington.


The supplement related to this article is available online at:

Author contributions

TZ and QHH planned and designed the research, performed the analyses, and wrote the paper with contributions from all co-authors. BG, DK, and ZM participated in the research and edited the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We extend our deepest appreciation to handling editor Monica Riva for their diligent assistance and to Hugh Daigle and one anonymous referee for their insightful comments on this paper.

Financial support

Financial assistance for this work was provided by the National Natural Science Foundation of China (grant nos. 41830431 and 41821002), PetroChina International Cooperation Project (grant no. 2023DQ0422), Shandong Provincial Major Type Grant for Research and Development from the Department of Science & Technology of Shandong Province (grant no. 2020ZLYS08), Maverick Science Graduate Research Fellowship for 2022–2023, the AAPG and the West Texas Geological Foundation Adams Scholarships, and Kansas State University through faculty start-up funds to Behzad Ghanbarian.

Review statement

This paper was edited by Monica Riva and reviewed by Hugh Daigle and one anonymous referee.


Abdassah, D. and Ershaghi, I.: Triple-porosity systems for representing naturally fractured reservoirs, SPE Formation Eval., 1, 113–127,, 1986. 

Bibby, R.: Mass transport of solutes in dual – porosity media, Water Resour. Res., 17, 1075–1081, 1981. 

Bock, H., Dehandschutter, B., Martin, C. D., Mazurek, M., De Haller, A., Skoczylas, F., and Davy, C.: Self-sealing of fractures in argillaceous formations in the context of geological disposal of radioactive waste, ISBN 978-92-64-99095-1, 2010. 

Brace, W. F., Walsh, J. B., and Frangos, W. T.: Permeability of granite under high pressure, J. Geophys. Res., 73, 2225–2236, 1968. 

Busch, A., Alles, S., Gensterblum, Y., Prinz, D., Dewhurst, D. N., Raven, M. D., Stanjek, H., and Krooss, B. M.: Carbon dioxide storage potential of shales, Int. J. Greenh. Gas Con., 2, 297–308, 2008. 

Carslaw, H. S. and Jaeger, J. C.: Conduction of Heat in Solids, ISBN 0198533683, 1959. 

Chen, X. and Pfender, E.: Effect of the Knudsen number on heat transfer to a particle immersed into a thermal plasma, Plasma Chem. Plasma P., 3, 97–113, 1983. 

Civan, F., Devegowda, D., and Sigal, R. F.: Critical evaluation and improvement of methods for determination of matrix permeability of shale, Society of Petroleum Engineers,, 2013. 

Cui, X., Bustin, A., and Bustin, R. M.: Measurements of gas permeability and diffusivity of tight reservoir rocks: different approaches and their applications, Geofluids, 9, 208–223, 2009. 

Devegowda, D.: Comparison of shale permeability to gas determined by pressure-pulse transmission testing of core plugs and crushed samples, Unconventional Resources Technology Conference (URTEC), San Antonio, Texas, 20–22 July 2015,, 2015. 

Egermann, P., Lenormand, R., Longeron, D. G., and Zarcone, C.: A fast and direct method of permeability measurements on drill cuttings, SPE Reserv. Eval. Eng., 8, 269–275, 2005. 

Fakher, S., Abdelaal, H., Elgahawy, Y., and El-Tonbary, A.: A Review of Long-Term Carbon Dioxide Storage in Shale Reservoirs, in: SPE/AAPG/SEG Unconventional Resources Technology Conference, D023S028R003,, 2020. 

Gensterblum, Y., Ghanizadeh, A., Cuss, R. J., Amann-Hildenbrand, A., Krooss, B. M., Clarkson, C. R., Harrington, J. F., and Zoback, M. D.: Gas transport and storage capacity in shale gas reservoirs – A review. Part A: Transport processes, Journal of Unconventional Oil and Gas Resources, 12, 87–122, 2015. 

Ghanbarian, B.: Estimating the scale dependence of permeability at pore and core scales: Incorporating effects of porosity and finite size, Adv. Water Resour., 161, 104123,, 2022a. 

Ghanbarian, B.: Scale dependence of tortuosity and diffusion: Finite-size scaling analysis, J. Contam. Hydrol., 245, 103953,, 2022b. 

Ghanbarian, B., Hunt, A. G., and Daigle, H.: Fluid flow in porous media with rough pore-solid interface, Water Resour. Res., 52, 2045–2058, 2016. 

Guidry, K., Luffel, D., and Curtis, J.: Development of laboratory and petrophysical techniques for evaluating shale reservoirs. Final technical report, October 1986–September 1993, ResTech Houston, Inc., TX, United States, PB-96-174859/XAB, 1996. 

Gutierrez, M., Øino, L. E., and Nygaard, R.: Stress-dependent permeability of a de-mineralised fracture in shale, Mar. Petrol. Geol., 17, 895–907, 2000. 

Haggerty, R. and Gorelick, S. M.: Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity, Water Resour. Res., 31, 2383–2400, 1995. 

Haskett, S. E., Narahara, G. M., and Holditch, S. A.: A method for simultaneous determination of permeability and porosity in low-permeability cores, SPE Formation Eval., 3, 651–658, 1988. 

Haszeldine, S., Lu, J., Wilkinson, M., and Macleod, G.: Long-timescale interaction of CO2 storage with reservoir and seal: Miller and Brae natural analogue fields North Sea, Greenhouse Gas Control Technology, 8, 2006. 

Heller, R., Vermylen, J., and Zoback, M.: Experimental investigation of matrix permeability of gas shalesExperimental Investigation of Matrix Permeability of Gas Shales, AAPG Bull., 98, 975–995, 2014. 

Hu, Q., Ewing, R. P., and Dultz, S.: Low pore connectivity in natural rock, J. Contam. Hydrol., 133, 76–83, 2012. 

Hu, Q., Ewing, R. P., and Rowe, H. D.: Low nanopore connectivity limits gas production in Barnett formation, J. Geophys. Res.-Sol. Ea., 120, 8073–8087, 2015. 

Hu, Q.: A pulse-decay method for low permeability analysis of granular porous media: I. Analytical solutions and experimental methodologies, Texas Data Repository [data set], V1,, 2020. 

Huenges, E.: Enhanced geothermal systems: Review and status of research and development, Elsevier,, 2016. 

Javadpour, F.: Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone), J. Can. Petrol. Technol., 48, 16–21, 2009. 

Jia, B., Tsau, J., Barati, R., and Zhang, F.: Impact of heterogeneity on the transient gas flow process in tight rock, Energies, 12, 3559,, 2019. 

Kamath, J.: Evaluation of accuracy of estimating air permeability from mercury-injection data, SPE Formation Eval., 7, 304–310, 1992. 

Khosrokhavar, R.: Shale gas formations and their potential for carbon storage: opportunities and outlook, Mechanisms for CO2 Sequestration in Geological Formations and Enhanced Gas Recovery, Environmental Processes, ISBN 9783319230870, 3319230875, 67–86, 2016. 

Kim, C., Jang, H., and Lee, J.: Experimental investigation on the characteristics of gas diffusion in shale gas reservoir using porosity and permeability of nanopore scale, J. Petrol. Sci. Eng., 133, 226–237, 2015. 

Kim, J., Kwon, S., Sanchez, M., and Cho, G.: Geological storage of high level nuclear waste, KSCE J. Civ. Eng., 15, 721–737, 2011. 

Liu, W., Li, Y., Yang, C., Daemen, J. J., Yang, Y., and Zhang, G.: Permeability characteristics of mudstone cap rock and interlayers in bedded salt formations and tightness assessment for underground gas storage caverns, Eng. Geol., 193, 212–223, 2015. 

Lu, J., Wilkinson, M., Haszeldine, R. S., and Fallick, A. E.: Long-term performance of a mudrock seal in natural CO2 storage, Geology, 37, 35–38, 2009. 

Luffel, D. L., Hopkins, C. W., and Schettler Jr., P. D.: Matrix permeability measurement of gas productive shales, Society of Petroleum Engineers, SPE,, 1993. 

Neuzil, C. E.: Can shale safely host US nuclear waste? Eos, Transactions American Geophysical Union, 94, 261–262, 2013. 

Neuzil, C. E.: Groundwater flow in low-permeability environments, Water Resour. Res., 22, 1163–1195, 1986. 

Peng, S. and Loucks, B.: Permeability measurements in mudrocks using gas-expansion methods on plug and crushed-rock samples, Mar. Petrol. Geol., 73, 299–310, 2016. 

Pini, R.: Assessing the adsorption properties of mudrocks for CO2 sequestration, Enrgy Proced., 63, 5556–5561, 2014. 

Profice, S., Lasseux, D., Jannot, Y., Jebara, N., and Hamon, G.: Permeability, porosity and klinkenberg coefficien determination on crushed porous media, Petrophysics, 53, 430–438, 2012. 

Rasmuson, A.: The effect of particles of variable size, shape and properties on the dynamics of fixed beds, Chem. Eng. Sci., 40, 621–629, 1985. 

Roy, S., Raju, R., Chuang, H. F., Cruden, B. A., and Meyyappan, M.: Modeling gas flow through microchannels and nanopores, J. Appl. Phys., 93, 4870–4879, 2003.  

Ruthven, D. M.: Principles of Adsorption and Adsorption Processes, John Wiley & Sons, ISBN 9780471866060, 0471866067, 1984. 

Ruthven, D. M. and Loughlin, K. F.: The effect of crystallite shape and size distribution on diffusion measurements in molecular sieves, Chem. Eng. Sci., 26, 577–584, 1971. 

Tarkowski, R.: Underground hydrogen storage: Characteristics and prospects, Renew. Sust. Energ. Rev., 105, 86–94, 2019. 

Wu, T., Zhang, D., and Li, X.: A radial differential pressure decay method with micro-plug samples for determining the apparent permeability of shale matrix, J Nat. Gas. Sci. Eng., 74, 103126,, 2020. 

Yang, M., Annable, M. D., and Jawitz, J. W.: Back diffusion from thin low permeability zones, Environ. Sci. Technol., 49, 415–422, 2015. 

Zhang, J. J., Liu, H., and Boudjatit, M.: Matrix permeability measurement from fractured unconventional source-rock samples: Method and application, J. Contam. Hydrol., 103663,, 2020. 

Zhang, T., Hu, Q., Chen, W., Gao, Y., Feng, X., and Wang, G.: Analyses of true-triaxial hydraulic fracturing of granite samples for an enhanced geothermal system, Lithosphere (Special 5), 3889566,, 2021. 


The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.

Short summary
Tight rock is essential to various emerging fields of energy geosciences such as EGS and CCUS, but its ultra-low permeability is not easily measurable as a rigorous and rapid theory-based measurement technique for sub-nanodarcy levels is lacking. For the first time, we resolve this by providing an integrated technique (termed gas permeability technique) with coupled theoretical development, experimental procedures, and a data interpretation workflow.