the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A pulsedecay method for low (matrix) permeability analyses of granular rock media
Tao Zhang
Behzad Ghanbarian
Derek Elsworth
Zhiming Lu
Nanodarcy level permeability measurements of porous media, such as nanoporous mudrocks, are frequently conducted with gas invasion methods into granularsized 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 latetime responses. Critically, one latetime solution is for an ultrasmall gasinvadable 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 molecularsieve 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, theorybased, rapid, and versatile gas permeability measurement technique for tight media at subnanodarcy levels.
Please read the corrigendum first before continuing.

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

Article
(4341 KB)
 Corrigendum

Supplement
(838 KB)

The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(4341 KB)  Fulltext XML
 Corrigendum

Supplement
(838 KB)  BibTeX
 EndNote

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.
Shales, crystalline, and salt rocks with low permeabilities (e.g., <10^{−17} m^{2} 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), longterm performance of highlevel 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), longterm sealing for carbon utilization and storage (Fakher et al., 2020; Khosrokhavar, 2016), and highvolume 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 lowpermeability 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 steadystate and pulsedecay methods use consolidated centimetersized coreplug samples, which may contain fractures and show dual or tripleporosity characteristics (Abdassah and Ershaghi, 1986; Bibby, 1981). The overall permeability may therefore be controlled by a few beddingoriented or crosscutting 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 (nonfractured) 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 ultratight rocks (with permeability <10^{−20} m^{2} or approximately 10 nd), particularly when permeability is close to the detection limit of the pulsedecay approach on core plugs at ∼10 nd (e.g., using the PoroPDP200 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 pulsedecay 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 pulsedecay data in crushed rock samples without an available published algorithm as this method shares different constitutive phenomena than the traditional pulsedecay method for coreplug samples in Cartesian coordinates and (2) present the associated experimental methodology to measure permeability, reliably, and reproducibly in tight and ultratight 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 nanodarcylevel permeability in tight rocks for emerging energy and resources subsurface studies.
For a compressible fluid under unsteadystate conditions, flow in a porous medium can be expressed by the mass conservation equation:
where p is the pressure, t is the time, ρ is the fluid density, and $\stackrel{\mathrm{\u203e}}{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 nonlinearity 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
The gas compressibility c_{t} is given by
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 nonideality 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 pseudopressure 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
where U_{s} 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. K_{c} 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 pulsedecay method proposed by Brace et al. (1968).
The fractional gas transfer for the internal (limited K_{c} value) and external (infinite K_{c} value) gas transfer of sample is given by
where F_{f} and F_{s} 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 latetime solution to Eq. (2c) for a limited K_{c} value (called LLT hereafter) is
The latetime solution to Eq. (2d) when K_{c} tends to infinity (ILT hereafter) is
The earlytime solution to Eq. (2d) when K_{c} approaches infinity (IET hereafter) is
In Eq. (3), R_{a} is the particle diameter of a sample and s_{1}, s_{2}, and s_{3} 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.
^{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 K_{c} value. ^{c} Equations (3a)–(3c) are three approximate solutions of density faction function F.
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 latetime solution is missing, and error analyses are not provided. In addition, the lack of detailed analyses of τ and K_{c} 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 earlytime and latetime 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 gaspermeabilitymeasurement approach holistically as the gas permeability technique (GPT).
As aforementioned, mathematical solutions given in Eqs. (3a) and (3b) were deduced based on different values of K_{c} 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 K_{c} 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 earlytime regime from the latetime one for the first time in the literature. We also demonstrate that the earlytime 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 K_{c} value for data quality control
To apply the GPT method, appropriately selecting the parameter K_{c} in Eqs. (3a)–(3c) is crucial, as it is a critical value for data quality control. The dimensionless density outside the sample, U_{f}, is related to K_{c} via Eq. (S33) in Sect. S2. One may simplify Eq. (S33) by replacing the series term with some finite positive value and set
We define ${K}_{\mathrm{f}}={K}_{\mathrm{c}}/(\mathrm{1}+{K}_{\mathrm{c}})$ to interpret the density variance of the system as K_{f} is closely related to the dimensionless density outside the sample, U_{f}.
Equation (4) shows the relationship between U_{f} and K_{c} (Fig. 1). For K_{c} >0, K_{f} falls between 0 and 1. The greater the K_{f} value, the more insensitive to density changes the system is. For K_{c} equal to 50, K_{f} would no longer be sensitive to K_{c} variations as it would have already approached 98 % of the dimensionless density. This means that the U_{f} value needs to be greater than 0.98, and this leaves only 2 % of the fractional value of U_{f} available for capturing gas density change. When K_{c} is 100, the left fractional value of U_{f} 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 K_{c} (less than 10) is recommended, as K_{f} nearly reaches its plateau beyond K_{c}=10 (Fig. 1). When K_{c} is 10, the left fractional value of U_{f} is only as low as 9 %.
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 latetime 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 latetime solutions.
Figure 2a shows the exact solution of F_{s} with their approximate early and latetime solution (Table 1). Two exact solutions of F_{f} where K_{c} equals 10 or 50 are also demonstrated in Fig. 2a. Figure 2b depicts the exact solution from F_{f} for different K_{c} 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 latetime 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 K_{c} value, Eq. (3c) would be more precise than Eqs. (3a) and (3b) and return results close to the exact solution for both F_{f} and F_{s}. Second, the results of Eqs. (3a) and (3b) presented in Fig. 2a are similar; their difference is very small, especially for K_{c}>10. Due to the fact that core samples from deep wells are relatively short in length and their void volume is small (ultralow porosity and permeability such as in mudrocks with k≤ 0.1 nd), in practice, a solution for $\mathrm{10}<{K}_{\mathrm{c}}<\mathrm{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 K_{c} 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.
3.3 Applicability of the earlytime solution
A small K_{c} 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 K_{c} solution of F_{f}, and the infinity K_{c} solution F_{s} is controlled by K_{c}. However, before the selection of K_{c}, the dimensionless time is the basic parameter to be estimated as a priori before the early or latetime solutions are selected.
For pulsedecay methods, the earlytime 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 pulsedecay experiments does not mean that the earlytime 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 ultralowpermeability samples with k≤0.1 nd and small dimensionless times τ<0.024. This situation would only be applicable to earlytime 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 earlytime solution would fit ultralowpermeability 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.
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 F_{f} and F_{s} 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 U_{s} and U_{f} are at least (2n+1)π^{2}. Therefore, the entire series expansion term can be omitted without being influenced by K_{c}. 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 F_{f}, the error differences between the exact and approximate solutions are 3.5 % and 0.37 % for τ=0.05 and 0.1 when K_{c}=10, respectively. When τ≤0.024, the error would be greater than 14.7 %. Figure 2(b) shows that F_{f} can be approximated as F_{s} when K_{c} is greater than 10; the error difference between F_{f} and F_{s} is quite small at this K_{c} value (for K_{c}=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 F_{s}, the error difference is roughly the same as F_{f} 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 K_{c} 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).
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 microfractures 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 X1 and sample X2. 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.
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 laminarflow or diffusive states. This helps obtain a linear trend for ln (1−F_{f}) or ${F}_{\mathrm{s}}^{\mathrm{2}}$ versus time for lowpermeability media. Experimentally, data from 30 to several hundreds of seconds are recommended for tight rocks like shales within the GPT methodology.
In the GPT approach, as mentioned earlier, Eq. (S33) holds for small K_{c} 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 K_{c} would prevent the gas flow from entering into a slippage state as the pressure difference would increase with increasing K_{c}. 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 diffusionlike 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, N_{2}, and Ar, and one sorptive gas, CO_{2} (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 X2 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 <N_{2}<Ar). In terms of pressure drop, argon exhibited the most significant decrease. In a constanttemperature 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 CO_{2}). Argon may provide a wider range of valid penetration zones given its longest decay time (when excluding the adsorbed gas CO_{2}).
Figure 7 shows that the pressure decay curve of the adsorptive gas CO_{2} is different from those of the inert gases used in this study. CO_{2} 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 longterm field observations (Haszeldine et al., 2006; Lu et al., 2009) were carried out to assess the sealing efficiency of mudrocks for CO_{2} storage. In fact, the GPT offers a quick and effective way to identify the adsorption behavior of different mudrocks for both laminarflow and diffusion states.
4.3 Pressure decay behavior for different granular sizes
We compared the pressure drop behavior of gas in the mudrock sample X1 with different granular sizes (averaging from 0.34 to 5.18 mm) using the same sample weight and K_{c}. Results based on the experimental data shown in Fig. 8 indicate that a largersized 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).
^{a} The results are from the SMP200 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 latetime solutions derived from the calculated permeability. For example, the verified dimensionless time would be >0.024 using the earlytime solution result and vice versa.
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 SMP200 instrument. The builtin functions of SMP200 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 curvefit the pressure decay data and determine the permeability. The GRI method built in the SMP200 only suggests the fitting procedure for data processing without publicly available details on the underlying mathematics. The intragranular permeabilities of mudrock sample X1 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 SMP200 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 SMP200 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 ultratight (subnanodarcy levels) samples as there is almost no laminarflow or diffusion state to be captured. The greater pressure drop for largersized 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 microfractures, which complicate the determination of true matrix permeability (Heller et al., 2014), the versatile GPT method can still provide sizedependent permeabilities for a wide range of samples (e.g., from submillimeter to 10 cm diameter fullsize 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 latesolution 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 earlytime 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 latetime 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 K_{c} 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 latetime solutions, as shown in Fig. 9.
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 lowpermeability measurements. The mathematical solutions for the transport coefficient in the GPT were derived for a spherical coordinate system, applicable from laminar flow to slippagecorrected Fickian diffusion. Of the three derived solutions, the earlytime solution is valid when the gas storage capacity K_{c} approaches infinity, while the other two are latetime solutions, valid when K_{c} is either small or tends towards infinity. We evaluated the derived solutions for a systematic measurement of ultralow 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 ultratight rock samples in just tens to hundreds of seconds, with good repeatability.
B_{ij}  Correction parameter for viscosity, constant 
c_{t}  Fluid compressibility, Pa^{−1} 
F_{f}  Uptake rate of gas outside the sample, dimensionless 
F_{s}  Uptake rate in the sample, dimensionless 
f_{1}  Intercept of Eq. (S40), constant 
K_{a}  Apparent transport coefficient defined as Eq. (S9), m^{2} s^{−1} 
K_{c}  Ratio of gas storage capacity of the total void volume of the system to the pore (including adsorptive and 
nonadsorptive transport) volume of the sample, fraction  
K_{f}  Initial density state of the system, fraction 
k  Permeability, m^{2} 
k_{s}  Permeability defined as Eq. (S8), m^{2} (Pa s)^{−1} 
L  Coefficient, unit for certain physical transport phenomenon 
M  Molar mass, kg kmol^{−1} 
M_{m}  Molar mass of the mixed gas, kg kmol^{−1} 
M_{i,j}  Molar mass for gas i or j, kg kmol^{−1} 
M_{s}  Total mass of sample, kg 
N  Particle number, constant 
p  Pressure, Pa 
p_{cm}  Virtual critical pressure of mixed gas, Pa 
p_{p}  Pseudopressure from Eq. (S1), Pa s^{−1} 
R_{a}  Particle diameter of sample, m 
R  Universal gas constant, 8.314 J (mol K)^{−1} 
r  Diameter of sample, m 
s_{1}  Slope of Eq. (S40), constant 
s_{2}  Slope of function Ln(1−F_{s}), constant 
s_{3}  Slope of function ${F}_{\mathrm{s}}^{\mathrm{2}}$, constant 
T  Temperature, K 
T_{cm}  Virtual critical temperature for mixed gas, K 
t  Time, s 
U_{f}  Dimensionless density of gas outside the sample, dimensionless 
U_{s}  Dimensionless density in grain, dimensionless 
U_{∞}  Maximum density defined as Eq. (S37), dimensionless 
V_{1}  Cell volume in upstream of pulsedecay method, m^{3} 
V_{2}  Cell volume in downstream of pulsedecay method, m^{3} 
V_{b}  Bulk volume of sample, m^{3} 
V_{c}  Total system void volume except for sample bulk volume, m^{3} 
$\stackrel{\mathrm{\u203e}}{v}$  Darcian velocity in pore volume of porous media, m s^{−1} 
X  Pressure force, Pa 
y_{i,j}  Molar fraction for gas i or j, fraction 
z  Gas deviation (compressibility) factor, constant 
α_{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}  Pseudodensity from Eq. (S1), kg m^{−3} s^{−1} 
ρ_{s}  Density of gas changing with time in sample, kg m^{−3} s^{−1} 
ρ_{ps}  Pseudodensity of gas changing with time in sample, kg m^{−3} s^{−1} 
ρ_{pf}  Pseudodensity of gas changing with time outside sample, kg m^{−3} s^{−1} 
ρ_{p2}  Initial pseudodensity of gas in sample, kg m^{−3} s^{−1} 
ρ_{p0}  Average pseudodensity 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 pseudodensity 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 (${\mathit{\varphi}}_{\mathrm{f}}={\mathit{\varphi}}_{\mathrm{a}}+{\mathit{\varphi}}_{\mathrm{b}}$) occupied by both free and adsorptive fluids, fraction 
This work did not use any data from previously published sources, and our experimental data and MATLAB processing codes are available at https://doi.org/10.18738/T8/YZJS7Y (Hu, 2020), managed by Mavs Dataverse of the University of Texas at Arlington.
The supplement related to this article is available online at: https://doi.org/10.5194/hess2744532023supplement.
TZ and QHH planned and designed the research, performed the analyses, and wrote the paper with contributions from all coauthors. BG, DK, and ZM participated in the research and edited the paper.
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 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 startup funds to Behzad Ghanbarian.
This paper was edited by Monica Riva and reviewed by Hugh Daigle and one anonymous referee.
Abdassah, D. and Ershaghi, I.: Tripleporosity systems for representing naturally fractured reservoirs, SPE Formation Eval., 1, 113–127, https://doi.org/10.2118/13409PA, 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.: Selfsealing of fractures in argillaceous formations in the context of geological disposal of radioactive waste, ISBN 9789264990951, 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, https://doi.org/10.2118/166473MS, 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 pressurepulse transmission testing of core plugs and crushed samples, Unconventional Resources Technology Conference (URTEC), San Antonio, Texas, 20–22 July 2015, https://doi.org/10.15530/urtec20152154049, 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 ElTonbary, A.: A Review of LongTerm Carbon Dioxide Storage in Shale Reservoirs, in: SPE/AAPG/SEG Unconventional Resources Technology Conference, D023S028R003, https://doi.org/10.15530/urtec20201074, 2020.
Gensterblum, Y., Ghanizadeh, A., Cuss, R. J., AmannHildenbrand, 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, https://doi.org/10.1016/j.advwatres.2022.104123, 2022a.
Ghanbarian, B.: Scale dependence of tortuosity and diffusion: Finitesize scaling analysis, J. Contam. Hydrol., 245, 103953, https://doi.org/10.1016/j.jconhyd.2022.103953, 2022b.
Ghanbarian, B., Hunt, A. G., and Daigle, H.: Fluid flow in porous media with rough poresolid 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, PB96174859/XAB, 1996.
Gutierrez, M., Øino, L. E., and Nygaard, R.: Stressdependent permeability of a demineralised fracture in shale, Mar. Petrol. Geol., 17, 895–907, 2000.
Haggerty, R. and Gorelick, S. M.: Multiplerate mass transfer for modeling diffusion and surface reactions in media with porescale 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 lowpermeability cores, SPE Formation Eval., 3, 651–658, 1988.
Haszeldine, S., Lu, J., Wilkinson, M., and Macleod, G.: Longtimescale interaction of CO_{2} 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 pulsedecay method for low permeability analysis of granular porous media: I. Analytical solutions and experimental methodologies, Texas Data Repository [data set], V1, https://doi.org/10.18738/T8/YZJS7Y, 2020.
Huenges, E.: Enhanced geothermal systems: Review and status of research and development, Elsevier, https://doi.org/10.1016/B9780081003374.000255, 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, https://doi.org/10.3390/en12183559, 2019.
Kamath, J.: Evaluation of accuracy of estimating air permeability from mercuryinjection data, SPE Formation Eval., 7, 304–310, 1992.
Khosrokhavar, R.: Shale gas formations and their potential for carbon storage: opportunities and outlook, Mechanisms for CO_{2} 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.: Longterm performance of a mudrock seal in natural CO_{2} 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, https://doi.org/10.2118/26633MS, 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 lowpermeability environments, Water Resour. Res., 22, 1163–1195, 1986.
Peng, S. and Loucks, B.: Permeability measurements in mudrocks using gasexpansion methods on plug and crushedrock samples, Mar. Petrol. Geol., 73, 299–310, 2016.
Pini, R.: Assessing the adsorption properties of mudrocks for CO_{2} 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 microplug samples for determining the apparent permeability of shale matrix, J Nat. Gas. Sci. Eng., 74, 103126, https://doi.org/10.1016/j.jngse.2019.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 sourcerock samples: Method and application, J. Contam. Hydrol., 103663, https://doi.org/10.1016/j.jconhyd.2020.103663, 2020.
Zhang, T., Hu, Q., Chen, W., Gao, Y., Feng, X., and Wang, G.: Analyses of truetriaxial hydraulic fracturing of granite samples for an enhanced geothermal system, Lithosphere (Special 5), 3889566, https://doi.org/10.2113/2022/3889566, 2021.
 Abstract
 Highlights
 Introduction
 Mathematical solutions for gas permeability of granular samples
 Practical usages of algorithms for the GPT
 Influence of kinetic energy on gas transport behavior
 Conclusions
 Appendix A: Nomenclature
 Appendix B: Greek letters
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
The requested paper has a corresponding corrigendum published. Please read the corrigendum first before downloading the article.
 Article
(4341 KB)  Fulltext XML
 Corrigendum

Supplement
(838 KB)  BibTeX
 EndNote
 Abstract
 Highlights
 Introduction
 Mathematical solutions for gas permeability of granular samples
 Practical usages of algorithms for the GPT
 Influence of kinetic energy on gas transport behavior
 Conclusions
 Appendix A: Nomenclature
 Appendix B: Greek letters
 Data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement