# Single Phase Compressible Gas Flow in Porous Media: Review and Advances

Rachid Ababou, Mohamed Haythem Bahlouli, Zakaria Saâdi, Israel Canamón Valera

► To cite this version:

Rachid Ababou, Mohamed Haythem Bahlouli, Zakaria Saâdi, Israel Canamón Valera. Single Phase Compressible Gas Flow in Porous Media: Review and Advances. Transport in Porous Media, 2025, 152 (11), pp.89-1007. http://dx.doi.org/10.1007/s11242-025-02226-9

HAL Id: hal-0532226

https://hal.science/hal-0532226v1

Submitted on 20 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Single Phase Compressible Gas Flow in Porous Media: Review and Advances

Rachid Ababou<fnref n="1" /> · Mohamed Haythem Bahlouli<fnref n="2" /> · Zakaria Saâdi<fnref n="3" /> · Israel Canamón Valera<fnref n="4" />

Received: 30 January 2025 / Accepted: 13 August 2025
© The Author(s) 2025

## Abstract

This article focuses on single phase compressible gas flow in porous media, especially hydrogen $H_2$ or other gases like air. It includes a comprehensive literature review on analytical approaches to gas flow, Klinkenberg effect, and other effects like gravitational acceleration (super-gravity cases). The review investigates previous findings for ideal gas flow under isothermal conditions under various conditions – including one-dimensional (1D) permeametric flow conditions – taking into account perfect gas compressibility and the Klinkenberg effect due to gas slippage in fine pores. Usually, gravitational acceleration is neglected in the gas flow literature: this classical assumption is assessed quantitatively, and a new 1D analytical solution is developed at steady state for the case of strong gravitational acceleration, as may arise under centrifugal conditions. On the other hand, new 1D analytical solutions are developed for space-time gas pressure profiles and for mass flux density profiles in the porous column, with or without Klinkenberg effect. These analytical solutions are tested and compared to numerical simulations, both Finite Volume and Finite Element. Both the gas pressure profiles and the mass flux density profiles approach the exact steady state at large times. Furthermore, it is demonstrated that the proposed analytical solution for gas pressure is a fair approximation over a broad range of time scales, from early times up to large times approaching steady state.

**Keyword** Porous media flow · Compressible gas flow · Darcy · Klinkenberg · Analytical solutions · Claystone

<figure><img src="image_2.png" /></figure>

<footnote>Rachid Ababou
rachidababou@gmail.com</footnote>

<footnote n="1">1 IMFT - Institut de Mécanique des Fluides de Toulouse, 31400 Toulouse, France</footnote>

<footnote n="2">2 IRSN - Institut de Radioprotection et de Sûreté Nucléaire, F-92260 Fontenay, France</footnote>

<footnote n="3">3 ASNR - Autorité de Sûreté Nucléaire et de Radioprotection, F-92260 Fontenay, France</footnote>

<footnote n="4">4 ETSIME - Universidad Politécnica de Madrid, Madrid, Spain</footnote>
# 1 Introduction and Summary

Our objectives in this paper are to review the literature on models and solutions for single phase compressible gas flow in porous media, but also, to clarify several effects (gravity, Klinkenberg), and to develop new analytical solutions (steady & transient gas flow). Note that the last section of the paper before conclusion contains a brief recapitulation collecting the various results and analyses on gas flow pressure and flux, on centrifugal acceleration effects, and on the Klinkenberg effect.

This paper is organized as follows.

Section 2 presents an overview and a critical review of the literature on compressible gas flow in porous materials; note that specific aspects of gas flow will be studied in more detail, along with original developments, in later sections of this paper.

Section 3 briefly introduce the physical compressible behavior of ideal gases under isothermal conditions, with a focus on hydrogen ($H_2$) and air.

Section 4 discusses the pressure-dependent gas permeability due to molecular gas slippage in tight pores (Klinkenberg effect), and in particular, discusses estimates of Klinkenberg's coefficient "b".

Section 5 presents the compressible gas flow equations and initiates analytical investigations in the case of a permeametric porous column: assessment of gravitational effect (usually negligible under natural gravity); steady state gas flow under strong acceleration field (e.g. centrifugal); steady state gas flow with negligible gravity with or without the Klinkenberg effect.

Section 6 presents a set of analytical solutions for transient gas flow in a porous column under permeametric pressure conditions. Gas pressure profiles $p(z, t)$ and mass flux density profiles $\phi(z, t)$ are predicted, and these analytical results are applied and compared to a numerical permeametric experiment, implemented using the Finite Volume code TOUGH2 and the Finite Element software platform COMSOL Multiphysics. These novel solutions include the case of Klinkenberg permeability with non zero coefficient b.

Section 7 presents a brief recapitulation of all salient analytical solutions presented in this paper.

Section 8 concludes the paper with a few additional remarks and perspectives regarding analytical studies of gas flow in porous media.

Finally, the paper is completed by three appendices containing details on the analytical solutions, respectively, for steady flow with strong acceleration fields (APPENDIX A, Sec. A), transient gas pressure profiles with Klinkenberg effect (APPENDIX B, Sec. B), and transient flux profiles with Klinkenberg effect (APPENDIX C, Sec. C).

# 2 Literature review on compressible gas flow

In this section, we develop a critical review of the literature on single phase compressible gas flow in porous materials. This overview will be completed in later sections of this paper where some aspects of the literature on gas flow in porous media will be re-analyzed in more detail, and possibly extended. But first, in order to delineate the topic at hand, let us briefly note some relevant features of porous media gas flow (in contrast with liquid flow): gases have a low density and they are significantly compressible; in particular, barotropic gases have pressure-dependent gas density $\rho(p)$; furthermore, at low pressures and in
narrow pores, gases can undergo molecular "slippage" effects; in particular this leads to a pressure-dependent permeability $k(p)$ that increases when pressure decreases (Klinkenberg effect).

Gas flow can occur in porous materials in a variety of cases in natural or man-made environments, for instance in-situ hydrocarbon sequestration and storage (e.g. Oldenburg's Chap.26 in Ho & Webb's 2006 Gas Transport book: Ho and Webb 2006), geysers in geothermal fields, and geological repositories of nuclear waste where hydrogen generation due to corrosion and radiolyzis can occur. Hydrogen generation in nuclear waste repositories has been studied in a number of recent works through numerical simulations and/or gas injection tests (in-situ analogs, often with nitrogen instead of hydrogen as the injected gas): see for instance (Amri et al. 2022; Bahlouli et al. 2023; De La Vaissière et al. 2014; Marschall et al. 2017; Saâdi et al. 2020; Wendling et al. 2019, 2013a, b, c), among others. Furthermore, laboratory experiments involving gas flow in porous media have been developed for various purposes, including for instance with centrifugal methods that enhance the acceleration field to several orders of magnitude beyond natural gravity g. (Note: this is discussed later in Sect. 5.4; a new analytical solution for steady gas flow in the presence of "super-gravity" is then presented in Sect. 5.5; and finally, a sensitivity analysis of the transient effects of strong acceleration fields is also presented in Sect. 6.9).

Let us now consider some specific works in the literature with a focus on analytical solutions for compressible gas flow in porous media.

## 2.1 Historical literature – a brief survey

In Section 5.4.2 of his 1960 book, Scheidegger (1960) developed a historical survey of "general unsteady state flow" of gases in porous media, neglecting gravity. He noted that the resulting equation is a nonlinear diffusion PDE of the form which he proposed to solve numerically by explicit finite differences (under a rather stringent stability condition). In the present paper, after developing the equations for transient compressible gas flow in a one-dimensional column (1D), we will analyze in Sect. 5 a complete solution of these equations under steady state conditions including the effects of the acceleration field (gravitational or centrifugal). We will then return to the transient equation for compressible gases later in Sect. 6. The contemporary literature on transient gas flow is reviewed further below. The reader is referred to Scheidegger 1960 (Scheidegger 1960) for a review of pre-1960 literature on compressible gas flow, although we also devote just below some attention to historical works by Knudsen 1909 (1909) and Klinkenberg 1941 (1941) on the molecular slip model, in addition to a review of more modern literature.

## 2.2 Knudsen, Klinenberg, molecular gas kinetics – overview

In Section 7.3 of his 1960 book, Scheidegger (1960) also developed a historical survey of "gas slippage" effects, starting with Knudsen's 1909 law (1909) which corrects Poiseuille or Darcy laws for molecular "slippage" effects due to molecular mean free path length being longer than tube radius or plate aperture. He then focused on the capillary tube model of molecular slippage in porous media, developed in 1941 by Klinkenberg (1941), a model which is still much in use today. See also the contemporary literature on gas flow, reviewed further below, and reviewed also in Sect. 4 focusing on gas flow with pressure-dependent permeability due to the Klinkenberg effect.
In order to anticipate on the subsequent review of contemporary literature on the Klinkenberg effect, we present here the pressure-dependent Klinkenberg modification of Darcy's permeability as influenced by molecular "slippage" effects in thin pores or tubes at low pressures:

$$k(p) = k_{\infty} \left( 1 + \frac{b}{p} \right) \qquad (2.1)$$

where $k(p)$ [$m^2$] is Klinkenberg's pressure-dependent permeability, $k_{\infty}$ [$m^2$] is the classical Darcy permeability (for large or infinite pressure, i.e., for negligible slippage effects), and $b$ [$Pa$] is the Klinkenberg coefficient.

In their 2018 paper, (Yang and Weigand 2018) investigate the Klinkenberg effect for gas flow in microporous and nanoporous media at the scale of a few grains and pores (typically 12 × 12 pores in 2D) using Bird's "Direct Simulation Monte-Carlo" method of molecular movements. Their underlying DSMC model, which was initiated by G.A. Bird in 1963 for rarefied molecular gas dynamics, is based on a statistical representation of molecular movements as particles undergoing hard-sphere collisions (see the 2013 document Bird 2013 concerning DSMC). The objectives of (Yang and Weigand 2018) were to validate the detailed coefficients proposed by Klinkenberg and others for the pressure-dependent permeability $K(p)$, or equivalently, for the ratio $k(p)/k_{\infty}$, where $k_{\infty}$ is the permeability at high pressure (negligible slippage effects). The ratio $k(p)/k_{\infty}$ increases as pressure decreases due to mean free path effects ("slippage").

The authors (Yang and Weigand 2018), reviewed extensively the various "formulas" proposed in the literature, which are summarized in their Table I in terms of a Knudsen number for comparison with their molecular DSMC results. Some of these "formulas" are of the form $k/k_{\infty} = 1 + a \times Kn$, where the Knudsen number $Kn$ needs to be defined more explicitly depending on molecular mean free path and geometry of pore network (or of single pores). Alternatively, for more realistic porous media, the usual Klinkenberg "regression" is expressed in terms of the macroscopic pore pressure $p$ as follows: $k(p)/k_{\infty} = 1 + b/p$.

We will further discuss these Klinkenberg-type "formulas", and the coefficients $a$ and $b$ defined above, in more detail later in Sect. 4, with a focus on porous materials with tight pores like claystones, and on hydrogen gas $H_2$.

## 2.3 A preview of various applications of gas flow

The 2006 book edited by Ho and Webb (2006) contains 27 chapters on various aspects of Gas Transport in Porous Media. In particular:

- In Chap.5, Webb 2006 (2006) reviews several aspects of gas transport mechanisms, including Darcian advection (neglecting gravity), the Darcy-Brinkman extension to take into account viscous shear in a porous medium near a wall, the Darcy-Forchheimer extension for fast gas flow (neglecting again gravity), and the Klinkenberg effect with pressure-dependent permeability for gas flow "low permeability media". The author (Webb 2006) reports two empirical "correlations" from the literature for the Klinkenberg coefficient $b$ in $k(p)/k_{\infty} = 1 + b/p$, for air flow at 25°C: this will be reviewed in more detail in Sect. 4 devoted to Klinkenberg permeability.
- In Chap.15, Tidwell 2006 (2006) presents steady-state air permeability measurements (gas mini-permeameters, in-situ borehole tests), taking into account ideal gas compressibility and other effects.

- In Chap.26, Oldenburg 2006 (2006) presents in-situ studies of carbon sequestration and storage, where CO2 & CH4 gas density and compressibility are considered for pressures up to 200 bars.

## 2.4 Review of analytical solutions for gas flow in porous media

In their 1998 article, (Wu et al. 1998) develop several analytical solutions for single phase gas flow in porous media involving Klinkenberg effects (as in Eq. 2.1), and they evaluate their validity using both numerical tests and experimental tests. They focus on the following 1D and 2D cases, all neglecting gravity:

- **Steady 1D gas flow** (z) (Wu et al. 1998). The exact pressure solution $p(z)$ is easily obtained in this case. The pressure profile $p(z)$ was compared to a steady state numerical simulation performed with the TOUGH2 code, which matched well. This steady state solution $p(z)$ was also used to illustrate determination of the Klinkenberg coefficient $b$ from permeametric tests with a compressible gas (their eqs. (5.1), (5.2), (5.3)), based on the difference and the arithmetic mean of inlet and outlet boundary pressures. Furthermore, steady-state permeametric gas flow experiments were presented (in their section 5.4): they were used to apply the above method (their eq.5.3) in order to identify both the classical Darcy permeability $k_\infty$ and the Klinkenberg coefficient $b$.

- **Steady 2D cylindrical flow** (r, z) (Wu et al. 1998). Steady state gas flow induced by gas injection from a vertical well: a quasi-analytical solution is proposed in terms of the variable $\Phi = (p + b)^2$ where "$b$ [$Pa$]" is the Klinkenberg coefficient; this solution involves an infinite series with Bessel functions in the arguments.

- **Transient 1D radial gas flow** (r, t) (1st variant) (Wu et al. 1998). This case represents horizontal gas flow towards a well with a given well flux rate $Q$: the problem is expressed in terms of the variable $P_B(r, t) = p(r, t) + b$, and a linearized solution is obtained by "freezing" the nonlinear diffusivity. (note: the variable $P_B = p + b$ was introduced initially in 1953 by Collins and Crawford (1953)).

This approximate solution contains a special function, the Exponential Integral, and it is analogous to Theis' solution of groundwater flow to a well. The solution was compared to a simulation using the TOUGH2 code; the radial pressure profile $p(r, t_{10days})$ matched well with the analytical solution at $t = 10$ days but there was no comparison of profiles versus time (a good fit would have probably been more difficult to obtain given the frozen diffusivity approximation). This approximate linearized solution was also used to illustrate determination of the Klinkenberg coefficient from vertical well tests with imposed flux rate $Q$, by measuring well pressure at two different times (their equations (5.4), (5.5)).

- **Transient 1D radial gas flow** (r, t) (2nd variant) Wu et al. (1998). This case represents horizontal gas flow towards a well with a fixed pressure $P_w$ imposed in the well-bore: the problem is expressed again in terms of the variable $P_B(r, t) = p(r, t) + b$, with the same linearization approximation ("frozen" diffusivity). The resulting linearized solution is expressed as a formal integral containing Bessel functions in the argument (explicit solutions are not furnished).
Various other analytical solutions have been developed in the porous media and the hydrocarbon reservoir literature, some of them similar to the previously reviewed works, but some others using somewhat different combinations of analytical or quasi-analytical approaches.

For instance, in their 2021 paper on gas production from gas wells in high-pressure gas reservoirs, (Zhang et al. 2021) analyze gas flow by taking into account gas compressibility and viscosity-pressure dependence, but also, porous medium "shrinkage" and "expansion" of "irreducible" liquid water. They dub their approach "dynamic material balance". One feature of their mathematical model is the use of a "pseudo-pressure" named $p_P$ and a related "pseudo-time" $t_P$, or similar transformed variables such as as "quasi-pressure" and "quasi-time". These transformed variables have been used instead of pressure and time in the 1960s literature in order to achieve a partial linearization of the transient gas flow equation (which is otherwise nonlinear due to gas compressibility, as will be seen). See (Russell et al. 1966; Al-Hussainy et al. 1966; Al-Hussainy and Ramey 1966) for the initial pseudo-variables concepts, which have also been presented as "quasi-variable" concepts in other more recent literature. The pseudo-pressure $p_P$ used in Zhang et al. (2021) is defined as follows (simplifying here for incompressible porous media and constant viscosity independent of pressure):

$$p_P = p_0 + \frac{1}{\rho_0} \int_{p_0}^{p} \rho(p')dp' \qquad (2.2)$$

where, for an ideal gas, $\rho(p)$ may be deduced from the law of perfect gases. In the isothermal case, $\rho(p)$ is of the form $ap$ (see next section), and the pseudo-pressure becomes equivalent to a rescaled quadratic pressure $p^2$. This is equivalent to the procedure proposed in several works (quoted earlier), i.e., using squared pressure for solving a linearized version of the transient gas flow equation.

Other notable contributions to analytical solutions for compressible gas flow in the porous media literature include the following:

- Singh and Kolditz (2012) investigated gas flow in porous media for different applications like air flow in soils, gas production, and CO2 injection in porous reservoirs underground. They focused on gas compressibility, and considered both isothermal and non-isothermal gas flow, in the latter case with $\rho(p, T)$ and $\mu(p, T)$. They used a few simplified analytical solutions to validate their gas flow code, both isothermal and non-isothermal.

- Ngianga Alalibo and Orukari Mercy (2016) developed approximate solutions for ideal gas flow in porous media considering both the isothermal case ($\rho(p)$ linear), and the adiabatic case ($\rho(p)$ power law). Essentially, they solved a 1D flow problem with boundary conditions on gas density $\rho$, and they used a perturbation approximation, expanding $\rho$ around $\rho_0$.

- Shugard et al. 2012, in a SANDIA report (Shugard and Robinson 2014), developed investigations of isothermal gas flow through powders. In addition to viscous gas flow, which is nonlinear because it involves $\rho(p)$, they also considered an additional "free molecular flow" or Fickian "diffusion along a concentration gradient", which is linear. Based on this, they examined simple steady state solutions for the 1D flow problem, and then considered applications to experimental cases.

- Other authors like (Aadnoy and Finjord 1996) considered liquid flow for a compressible rock, with both liquid and rock behaving elastically, additionally with both perme-
ability k and viscosity µ being pressure-dependent. They used a perturbation approach to solve the resulting nonlinear diffusive flow problem for a transient line sink in an infinite domain oil reservoir, ending up with a solution in terms of the Boltzmann transform $u = (x^2 + y^2)/(4t)$. They also discuss several issues with the details of the perturbation approach in relation to the oil reservoir literature (examining which terms cannot be neglected for consistent results).

# 3 Compressible isothermal behavior of ideal gases

## 3.1 Compressibility law for ideal gases

The compressibility of a fluid is the change in density that will be produced in the fluid by a specified change in pressure. Gases are, in general, highly compressible whereas most liquids have a very low compressibility. The isothermal equation of state for a perfect gas, or ideal gas, can be written (Ababou 2018):

$$pV = nRT \qquad (3.1)$$

where: $n = m/M$ is number of moles (m is total mass and M is molar mass), T is temperature in degrees Kelvin (°K), and R is the Boltzmann constant of perfect gases ($J/(mol.K)$). From the CRC Handbook of Chemistry and Physics (Haynes et al. 2017):

$$R = 8.314460 \text{ J/(Mole.K)}$$

We can also write the law of ideal gases as:

$$p = \rho \frac{RT}{M} \qquad (3.2)$$

where $\rho$ [kg/m³] is the mass density of the gas, and M is the molar mass of the gas molecule in kg/mole. For hydrogen $H_2$ and for dry air, we have:

$$M_{H2} = 2.015 \ 10^{-3} \text{ kg/mole} \quad ; \quad M_{AIR} = 28.963 \ 10^{-3} \text{ kg/mole}$$

The isothermal law of perfect gases can then be expressed as follows:

$$\rho(p) = \alpha_0 p \qquad (3.3)$$

with $\rho$ in [kg/m³], and the compressibility coefficient $\alpha_0$ given in [($kg/m^3$)/Pa] by the following expression:

$$\alpha_0 = \frac{M}{RT} \qquad (3.4)$$

For a given reference pressure $p_0$, the corresponding reference density is $\rho_0 = \alpha_0 p_0$, and the law of ideal gases can then be expressed as:

$$\rho(p) = \rho_0 \frac{p}{p_0} \qquad (3.5)$$

Alternatively, it is convenient to use the reference state ($p_0, \rho_0$) to define the compressibility $\alpha_0$, which does not depend on pressure. Then, since $\alpha = \alpha_0$ is constant (from the law of perfect gases), we can write:
$$ \alpha_0 = \frac{\rho_0}{p_0} \qquad (3.6) $$

where we may take $p_0 = 10^5 Pa = 1 bar$ and $T \approx 293 K$.

For the sake of comparison, let us recall briefly the case of liquids. The mass density of a weakly compressible elastic liquid (e.g., water) follows a very different elastic law:

$$ \rho(p) = \rho_0 \exp\left[\frac{p-p_0}{K}\right] \qquad (3.7) $$

where $K$ is the bulk compressive stiffness of water in Pascals ($K = 2.273 \ 10^9 Pa \approx 23000$ Bars for liquid water). Note that this elastic law holds only for small volume deformations, so that the above law can also be approximated linearly as: $\rho(p) \approx \rho_0 (1 + (p - p_0)/K)$. The inverse stiffness $1/K$ of the liquid plays the role of a compressibility coefficient, but the liquid compressibility law Eq. 3.7 is very different from the law of ideal gases Eq. 3.5.

## 3.2 Compressibility and other properties of H2 and Air

In this section we examine the isothermal compressible behavior of hydrogen ($H_2$) and of Air. The compressibility of hydrogen is important in our analysis because of our interest in studying $H_2$ transport and production by metal corrosion in deep geological disposals of radioactive waste (Saâdi et al. 2020; Saâdi 2024). The compressibility of Air is important in our analysis because of the presence of Air in underground structures like geologic repositories, and also possibly, because of our interest in studying Air (or Nitrogen) injection in claystone rocks (like COx) through borehole experiments to derive the rock hydraulic and transport properties (Amri et al. 2022). However, in this paper, we will focus more on $H_2$ in our analyses.

The molar masses of $H_2$ and of dry Air are:

$$ M_{H2} = 2.01568 \ 10^{-3} \text{ kg/mole} \quad ; \quad M_{AIR} = 28.9647 \ 10^{-3} \text{ kg/mole} $$

Let us consider hydrogen gas $H_2$ as well as dry Air, under standard conditions, that is at $p = 1$ bar and $T \approx 300$ K, and then also at $p = 5$ bars in order to assess the effect of compressibility on their mass density. We can begin with Eq. 3.4 to evaluate the compressibility $\alpha$ for the given gas at the given temperature. Then, we can use Eq. 3.3 to compute the mass density of that gas under those conditions ($p, T$), from which we can deduce a pair of reference values ($p_0, \rho_0$) for instance at $p = 1$ bar and $T \approx 300$ K. We obtain in this fashion the following results:

* For hydrogen, the compressibility $\alpha_{H_2}$ at temperature $T_{300} \approx 300$ K is:

$$ \alpha_{H2} = \frac{M_{H2}}{R T_{300}} = \frac{2.01568 \ 10^{-3}}{8.314460 \times 300} \approx 8.0810 \ 10^{-7} \left[ \frac{kg}{m^3 \times Pa} \right] \qquad (3.8) $$

The value of $\alpha_{H2}$ at $T_0 \approx 293$ K is not much different: $\alpha_{H2} = 8.2754 \times 10^{-7} (kg/m^3)/Pa$.

* For dry air, the compressibility $\alpha_{AIR}$ at standard temperature $T_{300} \approx 300$ K is:
$$ \alpha_{AIR} = \frac{M_{AIR}}{R T_{300}} = \frac{28.9647 \times 10^{-3}}{8.314460 \times 300} \approx 1.16122 \times 10^{-5} \left[ \frac{kg}{m^3 \times Pa} \right] \quad (3.9) $$

The value of $\alpha_{AIR}$ at $T_0 \approx 293 \text{ K}$ is not much different: $\alpha_{AIR} = 1.18896 \times 10^{-5} (\text{kg/m}^3)/\text{Pa}$.

Note that dry air is about 14 times more compressible than $H_2$, and accordingly, the mass densities of $H_2$ and of Air, at the reference state ($p_0, T_0$) = (1 bar, 293 K), are in a ratio of 1/14:

$$ \rho_{H2} = \alpha_{H2} \times p = 8.2754 \times 10^{-7} \times 10^{5} = 0.082754 \text{ kg/m}^{3} $$

$$ \rho_{AIR} = \alpha_{AIR} \times p = 1.18896 \times 10^{-5} \times 10^{5} = 1.18896 \text{ kg/m}^{3} $$

To sum up, under these standard conditions, air is roughly one thousand times "lighter" than water and $H_2$ is about 14 times "lighter" than air. Let us now compute gas density at (5 bars, 293 K), that is, with the pressure increased to 5 bars instead of 1 bar (noting that $1 \text{ bar} = 10^5 \text{ Pa}$):

$$ \rho_{H2} = \alpha_{H2} \times p = 8.2754 \times 10^{-7} \times 5 \times 10^{5} = 0.4137700 \text{ kg/m}^{3} $$

$$ \rho_{AIR} = \alpha_{AIR} \times p = 1.18896 \times 10^{-5} \times 5 \times 10^{5} = 5.9448 \text{ kg/m}^{3} $$

Under increased pressure conditions, with ($p, T$) = (5 bars, 293 K), the gas densities have increased in proportion to pressure, by the linear law $\rho(p) = \rho_0 \frac{p}{p_0}$ given by Eq 3.5. Thus, at 5 bars, $H_2$ is 5 times denser than at 1 bar, because density is a linear function of pressure (at constant temperature). Thus, hydrogen density $\rho_{H2}$ has become $0.4138 \text{ kg/m}^3$ at 5 bars instead of $0.0827 \text{ kg/m}^3$ at 1 bar.

Finally, concerning gas viscosities, the dynamic viscosities $\mu_{H2}$ and $\mu_{AIR}$ are not very sensitive to pressure nor to temperature. Briefly, for hydrogen, we have from various sources:

$$ T \approx 293 \text{ K}; p = 1 \text{ bar} : \mu_{H2} = 8.88 \times 10^{-6} \text{Pa.s} $$

$$ T = 300 \text{ K}; p = 1 \text{ bar} : \mu_{H2} = 8.9 \times 10^{-6} \text{Pa.s} $$

$$ T = 400 \text{ K}; p = 1 \text{ bar} : \mu_{H2} = 10.9 \times 10^{-6} \text{Pa.s} $$

$$ T = 273 \text{ K}; p = 1 \text{ bar} : \mu_{H2} = 8.40 \times 10^{-6} \text{Pa.s} $$

$$ T = 273 \text{ K}; p = 100 \text{ bars} : \mu_{H2} = 8.57 \times 10^{-6} \text{Pa.s} $$

where the first value $\mu_{H2} \approx 8.88 \times 10^{-6}$ was obtained by interpolation from a look-up table in the Tough2 code. Like hydrogen, air viscosity is not very sensitive to pressure and temperature. For air at about 293–300 K and 1 bar, we have $\mu_{AIR} \approx 18.5 \times 10^{-6}$. Thus, air is about twice more viscous than hydrogen (dynamically).

Remark. The compressibility and viscosity properties of $H_2$ gas were also illustrated in Table 11 in the TOUGH2 User Guide v2 (2012) for EOS-5 (Equation Of State 5): the compressibility coefficient $\alpha$, as defined by us, can be inferred from these tabulated values of density and pressure, confirming our evaluation of $\alpha$. Note however, as a warning, that their 2012 Table 11 contained erroneous units in the footnote of the table: the reader should refer instead to the present text for $R$ and for $M_{H2}$.
# 4 Klinkenberg effect: pressure-dependent permeability

In this section we review the Klinkenberg effect due to molecular slippage of gases in tight pores, and we analyze the related coefficient "b" which intervenes in the pressure-dependent Klinkenberg permeability. Gas flow solutions taking into account the Klinkenberg effects are developed in subsequent sections.

## 4.1 Klinkenberg effect: introduction

At sufficiently low pressures and densities, gas molecules in thin pores are affected by a "slip-page" effect which tends to increase the macroscopic velocity of the gas (compared to broader pores). This effect is due to the fact that the molecules' mean free path is smaller (or on the same order) as pores diameters. The consequence for the Darcy flow model is an increased effective permeability compared to classical Darcy permeability. This effect was initially studied by Knudsen in 1909 (1909), and then by Klinkenberg in 1941 (1941). According to the latter (Klinkenberg 1941), effective gas permeability at finite pressure (especially at low pressure) is pressure-dependent and is given by:

$$k_g(P) = k_\infty \left(1 + \frac{b}{P}\right) \qquad (4.1)$$

where $k_\infty$ is the permeability that would obtain at high pressures. That is, $k_\infty$ is the classical Darcy permeability for negligible Klinkenberg effect.

The Klinkenberg effect may have significant impact on gas flow, especially in low permeability media where pore diameters are on the same order (or finer) than the mean free path length of gas molecules. In such cases gas "slippage" increases gas velocity and therefore increases the apparent Darcy permeability. See Klinkenberg 1941 (1941), particularly his equations (5) and (12)-(13) where mean free path length "$\lambda$" and pore radius "$r$" intervene.

## 4.2 Klinkenberg effect: Knudsen number

In their "Table I", Yang & Weigand 2018 (2018) present a compilation of formulas from the literature for the ratio $K_a/k_\infty$, that is the ratio $K_g(P)/k_\infty$ in our notations, expressed in terms of gas Knudsen number $Kn$ and/or gas pressure $p$. The Knudsen number is defined just below by Eq. 4.3. Several of the Knudsen-based formulae were expressed linearly in terms of $Kn$ as follows:

$$K_g(P)/k_\infty = 1 + A Kn \qquad (4.2)$$

where the Knudsen number is defined by:

$$Kn = \frac{\lambda}{R_{PORE}} \qquad (4.3)$$

where $\lambda$ [m] is the gas molecular mean free path length, and $R_{PORE}$ [m] is a typical pore diameter (to be discussed below). The linear formulae $K_g(P)/k_\infty = 1 + A Kn$ were proposed by the following authors: Klinkenberg 1941 (1941) with $A = 4.148$; Sakhaee-Pour and Bryant 2011 (2011) with $A = 63/(4\pi) \approx 5.013$; Mohammadmoradi and Kantzas 2016 (2016) with $A = 9.62$; Zhao et al. 2016 (2016) with $A = 3.20$. Other authors invoked more complex functions of $Kn$, and yet others (including Knudsen 1909) invoked hybrid

<figure><img src="image_3.png" /></figure>
functions depending both on $Kn$ and on pressure $p$. It should be emphasized that the Knudsen number itself depends on pressure, that is $Kn = Kn(p)$, as will be seen shortly below.

## 4.3 Klinkenberg effect: the pressure coefficient $b$ [Pa]

On the other hand, recall that, according to Klinkenberg 1941, $K_g(P)/k_\infty = 1 + b/P$. Several authors have evaluated Klinkenberg's $b$ coefficient based on gas thermodynamics. Thus, Yang & Weigand 2018 (2018) used the gas Knudsen number "Kn" in their equation (8), with the mean free path defined thermodynamically by their equation (9) in terms of temperature, pore diameter, and pressure. They also summarized the various "correlations" obtained in the literature for the ratio $K_g(P)/k_\infty = 1 + b/P$ using the $Kn$ number (see their Table I). In several of the formulas compiled in that table, the authors expressed the ratio $K_g(P)/k_\infty = 1 + b/P$ in the form $1 + A \cdot Kn$, that is:

$$ \frac{K_g(P)}{k_\infty} \equiv 1 + \frac{b}{P} = 1 + A \times Kn \qquad (4.4) $$

...from which some typical expressions and values of our coefficient "b" could be inferred by using:

$$ b \approx A \times P \times Kn \qquad (4.5) $$

Note: we will see that $Kn$ is inversely proportional to $P$ so that the pressure dependence of $b$ disappears.

### Evaluation of Klinkenberg coefficient "b" from gas kinetics

We now proceed to evaluate the Klinkenberg coefficient "b" from Eq. 4.5. The Knudsen number $Kn$ appearing in Eq. 4.5 is defined as $Kn = \lambda/R_{PORE}$, where $\lambda$ [m] is the molecular mean free path length of the gas, and $R_{PORE}$ [m] is a typical pore radius defined from pore size statistics, such as the peak pore diameter, or else the geometric mean pore radius (which gives more weight to fine pores than the arithmetic mean). In Table I of Yang & Weigand 2018 (2018), the value of $A$ is typically between $A \approx 4$ and $A \approx 9$ depending on authors. The mean free path length $\lambda$ depends on the gas molecule and the temperature, and it is given by (for the hard spheres model of gas kinetics):

$$ \lambda = \frac{1}{P} \frac{k_B T}{\sqrt{2} \pi d_{MOLEC}^2} $$

where $d_{MOLEC}$ is some effective gas molecule "diameter", and $k_B$ is the Boltzmann constant related to $R$ via the Avogadro number: we have $k_B = 1.380649 \times 10^{-23} J/K$ (Joules per Kelvin degree), and $R = N_A \times k_B = 8.31446261815324 J \cdot mol^{-1} \cdot K^{-1}$, with $N_A = 6.02214076 \times 10^{23} mol^{-1}$ the Avogadro number. Therefore, in terms of the ideal gas constant $R$, the mean free path length $\lambda$ can be expressed as:

$$ \lambda = \frac{RT}{\sqrt{2} \pi N_A d_{MOLEC}^2} \frac{1}{P} \qquad (4.6) $$

The molecule "diameter" of $H_2$ in the kinetic sense is on the order of 3 Angstrom, that is $d_{MOLEC} \approx 3 \times 10^{-10} m$. This would yield finally, for $H_2$ at $\approx 300$ K:
$$ \lambda_{H2} \approx \frac{8.31446261815324 \times 300}{1.414 \times 3.14159 \times (6.02214076 \times 10^{23}) \times (3 \times 10^{-10})^2} \frac{1}{P} \approx 0.010 \times \frac{1}{P} \quad (4.7) $$

where pressure $P$ is in Pascals. It is important to note, from Eq. 4.6, that the mean free path length $\lambda$ depends inversely on pressure. For instance, for $P = 1 \text{ bar} = 1 \times 10^5 \text{ Pa}$, this yields $\lambda_{H2} \approx 0.010/10^5 \approx 10^{-7} \text{ m} = 0.1 \text{ µm}$.

Concerning pore radii in the COx claystone formation (disturbed or not), the reader is referred to the discussion in Amri et al. (2022). The complexity of accounting for a whole pore size distribution rather than a single pore size should be recognized. The peak value, or the "mode", of the pore radii distribution seems a reasonable choice. In samples of the COx host rock, it is typically about 20 nm ($20 \times 10^{-9}$ meters) for some investigators, or between 50 nm and 90 nm for others, or else, only between 4 nm and 6 nm for others, depending on measurement methods: see detailed discussion in Amri et al. (2022). Opting for a pore radius of 20 nm, that is $R_{PORE} \approx 20 \times 10^{-9} \text{ m}$, and combining with the previous Klinkenberg-type formulas in terms of Knudsen number and mean free path, finally leads to the numerical evaluation of coefficient "b" shown below:

$$ \begin{gathered} Kn = \frac{\lambda}{R_{PORE}}; \lambda_{H2} = 0.010 \times \frac{1}{P}; R_{PORE} \approx 20 \times 10^{-9} \text{ m}; \\ \Rightarrow b = A \times P \times Kn \approx A \times P \times \frac{0.010 \times \frac{1}{P}}{R_{PORE}} \approx A \times \frac{0.010}{R_{PORE}} \end{gathered} $$

With $A \approx 4$, and with $R_{PORE} \approx 20 \times 10^{-9} \text{ m}$ for the COx claystone, this would yield approximately:

$$ b \approx 20 \times 10^5 \text{ Pa} = 20 \text{ bars} $$

According to $K_g(P)/k_\infty = 1 + b/P$, the Klinkenberg gas permeability $K_g(P)$ of COx would be enhanced by a factor 21 for gas at atmospheric pressure (1 bar), according to $1 + b/P = 1 + 20/1 = 21$. On the other hand, $K_g(P)$ would only be increased by a factor 2 for gas at 20 bars pressure. For gas pressures such that $P \gg b = 20 \text{ bars}$ the Klinkenberg effect would become negligible, so that $K_g(P) \approx k_\infty$ will become independent of $P$.

### Evaluation of Klinkenberg coefficient "b" from empirical regressions

Other more empirical evaluations of the Klinkenberg coefficient "b" are available, for instance the empirical regression:

$$ b \approx a \times k_{\infty}^{-0.36} \quad (4.8) $$

is due to Jones 1972 (1972), as mentioned in Wu et al. (1998). Using this regression for $b$ and the few $b$ values provided in Table I of Wu et al. (1998), we infer the constant $a \approx 0.25166$. Applying it to the COx with $k_\infty \approx 10^{-20} \text{ m}^2$ yields $b \approx 40 \text{ bars}$ instead of our previous theoretical evaluation of 20 bars, which is still acceptable given the uncertainties on molecular gas kinetics, pore size distributions, etc. Finally, applying this empirical regression to the COx EDZ (Excavation Disturbed Zone), with permeability $k_\infty \approx 5 \times 10^{-18} \text{ m}^2$ instead of $10^{-20} \text{ m}^2$, would yield $b \approx 4.258 \text{ bars}$ instead of 20 to 40 bars. The smaller value of $b$ for the EDZ corresponds to the predominance of larger pores in the EDZ compared to the intact COx.

In addition, Webb 2006 (2006) reports two other empirical "correlations" from the literature for the Klinkenberg coefficient $b$ in $k(p)/k_\infty = 1 + b/p$ for air flow at 25°C: $b_{air} \approx 0.11 \times k_\infty^{-0.39}$ (from Heid et al. 1950), and $b_{air} \approx 0.98 \times k_\infty^{-0.33}$ (from Jones and Owens
1980), where $k_\infty$ is the classical permeability in $m^2$. The two formulas are in agreement in the range $10^{-14} m^2 < k_\infty < 10^{-17} m^2$, which is probably the range of validity of either formula. The author (Webb 2006) also provided a simple formula to adjust a given reference value $b_0$ to a different ideal gas in terms of its molecular weight ($M$), viscosity $\mu$, and temperature $T$ (the latter in Kelvin degrees):

$$b = b_0 (\mu/\mu_0) (M/M_0)^{-1/2} (T/T_0)^{1/2}.$$

### Evaluation of Klinkenberg coefficient "b" from experiments

The previous discussions on the value of the Klinkenberg "b" coefficient $H_2$ gas flow in the intact COx claystone, seems to be confirmed by direct experimental evidence from (Didier 2012), who studied gas flow and Klinkenberg effect in laboratory experiments on COx samples. Overall, this author found that $b$ varies experimentally between 8.5 bars and 24 bars, depending on the sample used (sample variability). The experimental results also depended on relative air humidity and on temperature, although it was found that $b$ does not vary much with temperature: for instance, $b \approx 18$ bars at $T = 23^\circ C$, and $b \approx 22$ bars at $T = 90^\circ C$.

## 4.4 Closing remarks on gas flow analyses with Klinkenberg effect

In closing, it will be interesting to incorporate the pressure-dependent Klinkenberg permeability $k_g(P)$ of Eq. 4.1 to obtain analytical solutions in terms of flux and pressure profile for compressible gas flow. Thus, in the sequel of this paper, we will insert the Klinkenberg permeability (with $b \neq 0$) in various forms of the 1D gas flow equations, applied to permeametric isothermal gas "flow experiments" under fixed pressures at both ends of the porous column.

The steady flow case is particularly simple, and it has already been studied by several authors including (Dullien 1979; Scheidegger 1960; Shugard and Robinson 2014) [his Sec.3.1.2 & Sec.4.6], and Wu et al. (1998) [their equations (3.1), (3.2) & (3.5)], among others. Nevertheless, for the sake of completion, we will illustrate in Sect. 5.7 the effect of the Klinkenberg coefficient "b" for steady state gas flow (pressure and flux).

Last but not least, note that the analytical solutions of transient gas flow pressure and flux, to be developed in Sects. 6, 6.2, and 6.7, also incorporate the Klinkenberg effect via coefficient "b".

# 5 Compressible gas flow equations & 1D analytical steady flow

Gas flow in porous media behaves very differently from incompressible or weakly compressible liquid flow. For instance, the steady permeametric flow of liquid water under pressure conditions at both ends of the column leads to a linear pressure profile, while compressible gas flow leads to a nonlinear pressure profile, as will be seen in this section shortly.

In this section, we first introduce the basic model equations for transient flow of a compressible gas such as hydrogen and air in a porous medium under isothermal conditions, and we then specialize the analysis to the hydrostatic case (to assess the influence of gravity) and then to the steady state case with pressures imposed at both ends of the porous
column. This leads to explicit expressions for steady state pressure profiles and flux profiles, including for the non-trivial case of super-gravity (e.g. centrifugal acceleration field).

## 5.1 Compressible Darcy's law

For 1D vertical flow with $Oz$ upwards, and for a compressible ideal isothermal gas, Darcy's law can be expressed as:

$$q_z = -\frac{k}{\mu} \left( \frac{dp}{dz} + \rho(p)g \frac{dz}{dz} \right) = -\frac{k}{\mu} \left( \frac{dp}{dz} + \rho(p)g \right) = -\frac{k}{\mu} \left( \frac{dp}{dz} + \alpha g p \right) \quad (5.1)$$

Since the gas is compressible, with variable density $\rho$, it is more convenient for mass balance purposes to express the mass flux by Darcy's law as follows:

$$\vec{\phi} = \rho(p)\vec{q} \rightarrow \phi_Z = \rho(p)q_Z = -\frac{k}{\mu} \alpha p \left( \frac{dp}{dz} + \alpha g p \right) \quad (5.2)$$

where $\vec{\phi}$ in $[(kg/s)/m^2]$ is the Darcy mass flux density vector, and $\phi_Z$ its vertical component. Note that the gravitational term will be neglected in several forthcoming sections, except in one section especially devoted to the gravitational effect.<fnref n="1" />

More generally, to account for the Klinkenberg effect, the classical Darcy permeability $k$ should be replaced with the pressure-dependent Klinkenberg permeability $k_g(p)$ of Eq. 4.1. This yields the more general Darcy's law for compressible gas in terms of mass flux density:

$$\phi_Z = \rho(p) q_Z = -\frac{k_\infty}{\mu} \left(1 + \frac{b}{p}\right) \alpha p \left(\frac{dp}{dz} + \alpha g p\right) \quad (5.3)$$

## 5.2 Compressible mass conservation law

Mass conversation law for 1D vertical flow of compressible isothermal can be expressed as follows (with denoting $\Phi$ the porosity):

$$\Phi \frac{\partial \rho}{\partial t} + \text{div}(\rho(p)\vec{q}) = 0 \rightarrow \Phi \frac{\partial \rho}{\partial t} + \frac{\partial(\rho(p)q_Z)}{\partial z} = S(z, t) = 0 \quad (5.4)$$

In the absence of gas source ($S=0$), and without Klinkenberg effect ($b=0$), inserting Darcy's law Eq. 5.2 for the mass flux $\phi_Z = \rho(p)q_Z$ yields:

$$\Phi \frac{\partial \rho}{\partial t} = \frac{k}{\mu} \left[ \frac{\partial}{\partial z} \left( \rho(p) \frac{\partial p}{\partial z} \right) + \frac{\partial}{\partial z} \left( \rho^2(p) p \right) \right] \quad (5.5)$$

Let us now insert the ideal gas law with $\rho(p) = \alpha p$, with $\alpha$ the constant non-zero compressibility of the gas, and let us divide the entire equation by $\alpha$; this gives:

<footnote n="1">¹ In the sequel, when gravity is neglected, the 1D flow coordinate will be indifferently named "x" or "z" depending on context.</footnote>
$$ \Phi \frac{\partial p}{\partial t} = \frac{k}{\mu} \left[ \frac{\partial}{\partial z} \left( p \frac{\partial p}{\partial z} \right) + \frac{\alpha g L}{L} \frac{\partial}{\partial z} (p^2) \right] \quad (5.6) $$

Finally, we note that the above transient flow equation has the Right-Hand Side operator acting on the squared pressure $p^2$. This can be made more explicit by letting $p \frac{\partial p}{\partial z} = (1/2) \frac{\partial p^2}{\partial z}$, whence the final expression of the nonlinear transient gas flow equation:

$$ \Phi \frac{\partial p}{\partial t} = \frac{k}{\mu} \left[ \frac{1}{2} \frac{\partial^2}{\partial z^2} (p^2) + \frac{\alpha g L}{L} \frac{\partial}{\partial z} (p^2) \right] \quad (5.7) \\ \varphi_Z = \rho(p) q_Z = - \frac{k}{\mu} \alpha p \left( \frac{dp}{dz} + \frac{\alpha g L}{L} p \right) $$

Several remarks are in order concerning Eq. 5.7 (recalling also that it is expressed for the case b=0):

* If the porous medium is compressible (Hydro-Mechanical coupling) then porosity $\Phi$ must be left inside the time derivative $\frac{\partial(\phi_{z(t)})}{\partial t}$; but we neglect this effect here: porosity is constant and only the gas is compressible (not the porous medium).

* In addition we have also neglected any pressure dependence effect on gas viscosity $\mu$.

* We have re-arranged the term $\alpha g$ in the form of the dimensionless gravitational group $\alpha g L$ divided by $L$ (this to facilitate dimensional and scale analyses).

* Alternatively, we can also use the ideal gas law by expressing compressibility $\alpha$ in the form $\alpha_0 = \rho_0/p_0$. Accordingly, we can first choose some reference pressure $p_0$, and secondly, we can define $\rho_0 = \alpha_0 p_0$. Then the previous equation (Eq. 5.7) can be expressed equivalently by substituting $\alpha_0 = \frac{\rho_0}{p_0}$.

Note that the transient gas flow Eq. 5.7 is a nonlinear parabolic PDE (Partial Differential Equation). Its solution, numerically or analytically, would require an initial condition $p(z, 0)$, and two boundary conditions, for instance Dirichlet conditions such as $p(z_{BOT}, t)$ and $p(z_{TOP}, t)$. To be complete, note also that the pressure solution $p(z, t)$ of Eq. 5.7 can be used to obtain also the transient mass flux density $\varphi_Z(z, t)$ $[(kg/s)/m^2]$ from Darcy's mass flux Eq. 5.2.

Next, we will look successively at two specializations of the transient gas flow equations Eq. 5.7 & Eq. 5.2: (i) first the almost trivial case of gas statics in a column ($\varphi_Z(z, t) = 0$); and (ii) secondly, the case of steady state gas flow in a column ($\partial p/\partial t = 0$). The static case is useful to quantify the effects of gravity in a simple way for a static gas column, and the steady case also is of interest in itself, e.g., steady gas flow permeameter (easy to solve exactly with negligible gravity), or steady gas flow systems which appear less trivial to solve exactly if the acceleration field $g$ is not negligible, as will be seen (e.g. centrifuge experiments).

## 5.3 Compressible gas "hydrostatics" and the gravitational group

This subsection aims at evaluating quantitatively how negligible is the effect of natural gravity by looking at the gas pressure profile in a vertical column under static conditions ("gas statics"). This simple analysis will show that the vertical pressure profile of a static gas increases exponentially with depth (instead of linearly for an incompressible liquid). The exponential coefficient is the dimensionless "gravitational group" $\alpha g L$.
As will be seen shortly, this static gas analysis leads to a quantitative evaluation of the error due to neglecting gravity, a common assumption in the literature on gas flow. It is suggested here that in cases like centrifugal gas flow, the acceleration $g *$ in the porous column will be artificially increased to multiples of $g = 9.81 \text{ m}^2/\text{s}$ (e.g. supergravity $g * = 1000 \text{ g}$ or much more): in such cases, the static analysis below indicates that the acceleration $g *$ can no longer be neglected beyond 10 000 g or so.

Note also that, under natural gravity $g = 9.81 \text{ m}^2/\text{s}$, the gravitational effect in a static gas column depends also on the length $L$ of the column. Consider first a vertical porous column of vertical length $L$. For example, with $L = 500 \text{ m}$, the column could represent the vertical well in a High Level Radioactive Waste repository after the end of the back-filling and closure operations. Concerning the possible porosity and permeability of the backfill material, see for instance (Brommundt et al. 2014), (Saâdi et al. 2020), and (Wendling et al. 2019). See for instance Table 1 in Wendling et al. (2019) for properties of the backfill and other materials. We will see that in the case of a 500 m long vertical column, gravity can be indeed completely neglected (as many other authors have assumed).

Let us now express the gas statics equation, simply deduced by inserting zero mass flux in Darcy's mass flux law Eq. 5.2 for an ideal gas. Thus, starting with a known pressure at $z = 0$, the static gas pressure distribution is obtained as follows (without neglecting gravitational acceleration):

$$0 = \varphi_Z = \rho(p)q_Z = -\frac{k}{\mu} \alpha p \left(\frac{dp}{dz} + \alpha g p\right)$$

whence:

$$\frac{dp}{dz} + \alpha g p = 0$$

and:

$$\frac{dp(z)}{p(z)} = -\alpha g \, dz$$

Integrating upwards from $z=0$ to $z=Z$ yields:

$$\int_{p(0)}^{p(Z)} \frac{dp(z)}{p(z)} = -\alpha g (Z - 0)$$

Whence finally, the static gas pressure profile (Z upwards):

$$p(Z) = p(0) \exp\{-\alpha g Z\} \qquad (5.8)$$

Or equivalently, if gas pressure is imposed at the top ($z = L$), for instance $p(L) = p_{ATM} = 1 \text{ bar}$, then:

$$p(z) = p(L) \exp\{+\alpha g(L - z)\} \qquad (5.9)$$

In summary, the static vertical pressure profile of a compressible gas like $H_2$, under a gravitational field $\vec{g} = [0, 0, -g]$, is not linear but exponential.

However, if the dimensionless group $\alpha gL$ is much smaller than unity, that is:
$$ if \alpha gL \ll 1, \qquad (5.10) $$

then it can be shown by a Taylor expansion of the exponential, and by using our definition of compressibility:

$$ \alpha = \rho(p)/p \Rightarrow \alpha = \rho(L)/p(L), \qquad (5.11) $$

...that the pressure profile in the column remains indistinguishable from a linear profile:

$$ p(z) \approx p(L) - \rho(L)gz. \qquad (5.12) $$

For $H_2$ in a porous column of height $L = 500 m$, we have $\alpha gL = 0.396 \times 10^{-2} \approx 0.4\%$ which is indeed much less than unity. If the top pressure is 1 bar, that is $p(L) = p_{ATM} = 1$ bar, the static bottom pressure will not be much different: $p(0) = 1.003969$ bar. It would take a vertical column of height 12 km or more for the discrepancy between bottom and top pressure to reach +10% or more (for reference, note that the deepest human-made drill at the surface of the earth is about 12 km deep).

Thus, for most practical purposes, the effects of natural gravity $g = 9.81 m^2/s$ will be negligible for ideal gases like $H_2$, except perhaps for two specific types of applications: deep drilling to several kilometers depth, and laboratory experiments with centrifuges where the sample can be submitted to accelerations much larger than gravity.

## 5.4 Gas flow in the presence of centrifugal acceleration (discussion)

The case of centrifugal gas flow in porous media may be worth considering. For instance, with a meter long sample ($L = 1 m$), it would take a centrifugal acceleration $g \approx 10 000$ g to induce a notable 10% effect on the pressure profile in a one meter long, gas saturated sample obstructed at both ends (gas statics, no flow).

In the literature, centrifuges have been used for various purposes, e.g. for acting on capillary pressure in two-phase water/air retention and permeability curves measurements, or for porosimetry measurements, or else, for solute separation via pore water extraction (centrifuges are also used to study air flow through quasi-water saturated soil samples). Depending on applications, sample length is on the order of one to several decimeters, and acceleration $g$ ranges from tens of g's up to ten thousand g's and more (with centrifuges rotating at up to tens of thousands rotations per minute): see for instance (Yokoyamaama et al. 2011) for pore water extraction at centrifugal speeds up to 20 000 rpm (rotations per minute, 1 rpm being roughly equivalent to 1 g in laboratory scale experiments), and (Marulanda et al. 2000) concerning an experimental centrifugal model of air sparging for remediation of contaminated aquifers (air flow was observed in a transparent porous medium at varying gravitational accelerations, with recorded images of flow patterns, breakthrough velocities, and plume shapes as a function of g-level and injection pressure). In such cases, acceleration fields on the order of $g \approx 10 000$ g or more should not be neglected in the gas flow equations, whether gas flows as a single phase fluid (the case treated in this paper) or whether it is coupled to liquid flow.

We conclude that there is some interest in accounting for the acceleration field in our focused study of single phase gas flow. We will examine first its effect on steady state gas flow in the next sub-section, and then later, its effect on transient gas flow. (Note: centrifugal acceleration $g$ will be renamed $g$ in the gas flow equations when there is no ambiguity).
## 5.5 Analytic steady-state gas pressure and flux with gravity acceleration

We now summarize the exact analytical solution for pressure $p(z)$ at steady state for compressible gas flow without neglecting gravity, or more generally, taking into account the acceleration field $g$. The detailed calculations are presented in APPENDIX A (Sec. A).

The steady-state version of Eq. 5.7 is:

$$ 
\begin{aligned}
0 &= \frac{k}{\mu} \left[ \frac{1}{2} \frac{\partial^2}{\partial z^2} (p^2) + \frac{\alpha g L}{L} \frac{\partial}{\partial z} (p^2) \right] \\
\varphi_Z &= \rho(p) q_Z = -\frac{k}{\mu} \alpha p \left( \frac{dp}{dz} + \frac{\alpha g L}{L} p \right)
\end{aligned} 
\quad (5.13) $$

where $z$ is upwards, and $g$ is the acceleration field oriented downwards (earth gravity is $g = 9.81 \text{ m}^2/\text{s}$ but we might also consider centrifugal accelerations which are ten thousand times this value or more).

The solution of Eq. 5.13 is obtained after two changes of variables as detailed in APPENDIX A (Sec 8). The steady state pressure profile $p(z)$ between two fixed pressure boundaries at $z_{BOT}$ and $z_{TOP}$ is given by:

$$ p(z) = \sqrt{p^2(z_{BOT}) e^{-2\alpha g(z-z_{BOT})} - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} \left[1 - e^{-2\alpha g(z-z_{BOT})}\right]} \quad (5.14) $$

where $p(z)$ is expressed in terms of the bottom pressure and the steady mass flux density. Squaring this expression on both sides and applying it at $z = z_{TOP}$ yields the constant mass flux density $\varphi_0 \left[ (kg/s)/m^2 \right]$ in terms of both top and bottom pressure, as follows:

$$ \varphi_0 = -\frac{k}{\mu} \alpha^2 g \frac{p^2(z_{TOP}) - p^2(z_{BOT}) e^{-2\alpha g L}}{1 - e^{-2\alpha g L}} \quad (5.15) $$

Finally, inserting Eq. 5.15 in Eq. 5.14 gives a closed form expression for the pressure profile in terms of the top and bottom pressures. This expression, shown below, depends only on the boundary pressures, the compressibility coefficient $\alpha$, and the acceleration field $g$:

$$ p(z) = \sqrt{p^2(z_{BOT}) \left\{ e^{-2\alpha g(z-z_{BOT})} - \Psi(z) e^{-2\alpha g(z_{TOP}-z_{BOT})} \right\} + p^2(z_{TOP}) \Psi(z)} \quad (5.16) $$

with the "scaled pressure profile" $\Psi(z)$ given by Eq. 5.17 (reproduced below from APPENDIX A):

$$ \Psi(z) = \frac{1 - e^{-2\alpha g(z-z_{BOT})}}{1 - e^{-2\alpha g(z_{TOP}-z_{BOT})}} \quad (5.17) $$

which satisfies $\Psi(z_{BOT}) = 0$ and $\Psi(z_{TOP}) = 1$.

## 5.6 Steady state gas flow with negligible or zero gravity

The zero-gravity solution of the previous problem can be easily inferred directly from the steady state equations if one neglects gravity in these equations, as will be shown briefly
below. But we also show in this section that the previous steady flow solution with gravity field $g$ goes asymptotically to the exact zero-gravity solution as the dimensionless group $\alpha gL$ goes to zero: this is achieved by a Taylor development from the gravitational solution given by Eq. 5.15 for the flux, and Eq. 5.16 for the pressure profile.

Let $\epsilon = \alpha g L$, and assume $\epsilon \ll 1$ (negligible acceleration "g"). Note that for natural gravity $g = 9.81 \text{ m/s}^2$, with a column length $L = 0.25 \text{ m}$, and for $H_2$ gas under standard conditions, we have indeed $\epsilon = 2.0285 \times 10^{-6} \ll 1$, so natural gravity is completely negligible (and it remains negligible even for much longer column lengths such as $L = 500 \text{ m}$). Noting that $\alpha^2 g = \alpha \epsilon/L$, the steady mass flux expression Eq. 5.15 can be cast in the form:

$$ \varphi_0 = -(k/\mu)(\epsilon \alpha/L) N_0/D_0 $$

where, noting that $z_{TOP} - z_{BOT} = L$, the numerator $N_0$ and denominator $D_0$ can be expressed as:

$$ N_0 = p_{TOP}^2 - p_{BOT}^2 \exp(-2\epsilon) $$

$$ D_0 = 1 - \exp(-2\epsilon) $$

Using Taylor expansion of the exponential terms yields asymptotically, for $\epsilon \ll 1$:

$$ N_0 \approx p_{TOP}^2 - p_{BOT}^2 (1 - 2\epsilon) $$

$$ D_0 \approx 2\epsilon $$

whence:

$$ \varphi_0 \approx -(k/\mu)(\epsilon \alpha/L) \frac{p_{TOP}^2 - p_{BOT}^2 (1 - 2\epsilon)}{2\epsilon} $$

and finally, for $\epsilon \ll 1$:

$$ \varphi_0 \approx -\frac{k \alpha}{\mu 2L} (p_{TOP}^2 - p_{BOT}^2) \qquad (5.18) $$

This expression of the steady state mass flux density holds for negligible gravity ($\epsilon = \alpha gL \ll 1$). It is has the form of a "quadratic Darcy law", and it is satisfied at any position $z$ within the porous column. This result has been obtained and used since the 1960's by many authors who directly neglected gravity by setting $g = 0$ in the flow equations (e.g. Scheidegger 1960; Wu et al. 1998, etc.).

We now analyze the pressure profile expression $p(z)$ of Eq. 5.16 with the scaled pressure function $\Psi(z)$ of Eq. 5.17, asymptotically for negligible gravity (that is $\epsilon = \alpha gL \ll 1$). First, we insert the previous zero-gravity flux Eq. 5.18. Then, as before, we note that $z_{TOP} - z_{BOT} = L$, we express $\alpha g = (\alpha gL)/L = \epsilon/L$, and we express all exponential terms by their first order Taylor expansion, such as:

$$ \exp(-2 \alpha g (z - z_{BOT})) = \exp(-2 \epsilon (z - z_{BOT})/L) \approx 1 - 2 \epsilon (z - z_{BOT})/L $$

In this way, after reducing and neglecting terms of order $O(\epsilon)$ in expressions like $1 - 2 \epsilon (z - z_{BOT})/L$, we obtain finally:

$$ p(z) = \sqrt{p_{BOT}^2 + (p_{TOP}^2 - p_{BOT}^2) \frac{z - z_{BOT}}{L}} \qquad (5.19) $$
This steady state pressure profile corresponds exactly to the zero-gravity case $g = 0$, although it was obtained here asymptotically by letting $\epsilon = \alpha gL$) go to zero in the gravitational pressure solution given by Eq. 5.16. Note that here, contrary to the case with gravity acceleration, the steady pressure profile without gravity (or negligible gravity) does not depend on the compressibility coefficient $\alpha$.

*Important remarks.* The above analytical results indicate that, under natural gravity conditions, gravity is negligible, the steady pressure profile does not depend on gas compressibility, and therefore, a steady state permeametric gas flow experiment cannot be used to identify the compressibility coefficient $\alpha$. On the other hand, the permeametric experiment can still differentiate between an ideal gas and a liquid, because the nonlinear steady pressure profile of the compressible gas (Eq. 5.19) is much different from the linear pressure profile that would occur for an incompressible or weakly compressible liquid.

*Numerical example.* Anticipating on the numerical permeametric experiments of $H_2$ gas flow in the next section, let us provide a numerical illustration for the analytical flux (Eq. 5.18) and pressure profile (Eq. 5.19) obtained just above for negligible gravity ($g = 9.81 \text{ m/s}^2$ gives the same results as $g = 0$). The numerical data to be used are:

* Permeability: $k = 5 \times 10^{-18} \text{ m}^2$

* Dynamic viscosity: $\mu = 8.88 \times 10^{-6} \text{ Pa}\cdot\text{s}$

* Gas compressibility: $\alpha = 8.271 \times 10^{-7} (\text{kg/m}^3)/\text{Pa}$

* Length of porous column: $L = 0.25 \text{ m}$

* Boundary pressures: $p_{\text{BOT}} = 1.0 \times 10^5 \text{ Pa}$ and $p_{\text{TOP}} = 4.0 \times 10^5 \text{ Pa}$

The steady mass flux density of $H_2$ gas for these data is then, from previous relations:

$$ \varphi_0 \approx -\frac{5 \times 10^{-18}}{8.88 \times 10^{-6}} \frac{8.271 \times 10^{-7}}{2.0.25} [(4.0 \times 10^5)^2 - (1.0 \times 10^5)^2] \approx -1.397 \times 10^{-7} (\text{kg/s})/\text{m}^2 \quad (5.20) $$

where the minus sign is due to the fact that gas flow is downwards while the z-axis is upwards. This steady state value of analytical mass flux density will be compared in Sect. 6.7 to that obtained from permeametric gas flow simulations performed with the TOUGH2 Finite Volume code, and with the COMSOL Multiphysics Finite Element software: the three steady flow results coincide exactly. That is, mass flux density reaches a spatially constant profile $\varphi = -1.397 \times 10^{-7} (\text{kg/s})/\text{m}^2$ at large times (one day) both with TOUGH2 and COMSOL, as predicted analytically. Finally, the analytical pressure profile for the above input data is given by the following expression (taking here $z_{\text{BOT}} = 0$):

$$ p(z) = 10^5 \times \sqrt{1 + 15 \frac{z}{L}} \quad [\text{Pa}] \quad (5.21) $$

This analytical pressure profile will be plotted and compared to permeametric gas flow simulations with the TOUGH2 code and with COMSOL Multiphysics in Sect. 6.6. Focusing on steady state, the large times pressure profiles obtained with TOUGH2 and with COMSOL at $t \approx 1 \text{ day}$ or more, are indeed indistinguishable from the exact analytical pressure profile given just above.
## 5.7 Steady state gas flow with Klinkenberg permeability

In this section, we calculate in detail the steady state gas flow solution incorporating the Klinkenberg effect (but neglecting gravity, which excludes the case of centrifugal acceleration). This steady solution is not novel: see several authors including (Schei-degger 1960; Dullien 1979) (his Sec.3.1.2 & Sec.4.6), (Wu et al. 1998) (their equations (3.1), (3.2) & (3.5)), (Shugard and Robinson 2014), and others. The (exact) steady gas pressure profile $p(z)$ obtained just below coïncides with the existing literature for that case. Here we provide more details and interpretations, and we include also the analytical expression of the steady state mass flux density $\varphi_0$. In addition, the steady solution is used to illustrate the effect of the Klinkenberg coefficient "b" on pressure profiles and flux (see figures below). Also, finally, this steady solution will be useful for testing full transient solution of permeametric gas flow incorporating the Klinkenberg effect, to be developed further below in Sect. 6.2. Indeed, these transient solutions will be compared at large times with the exact steady solution developed in the present section.

The zero-gravity Klinkenberg flow problem can be expressed as follows by modifying the previous transient flow equations (respectively Darcy Eq. 5.2 and mass conservation Eq. 5.4), first by neglecting gravity ($g \approx 0$), secondly by replacing the classical Darcy permeability $k$ with the pressure-dependent Klinkenberg permeability $k_g(P)$ of Eq. 4.1. This yields the following system of transient flow equations:

$$ 
\begin{aligned}
\varphi_Z &= \rho(p) q_Z = -\frac{k_\infty}{\mu} \left(1 + \frac{b}{p}\right) \alpha p \left(\frac{dp}{dz}\right) \\
\Phi \frac{\partial \rho}{\partial t} + \frac{\partial(\rho(p)q_Z)}{\partial z} &= 0
\end{aligned}
\qquad (5.22) $$

from which we deduce now the corresponding steady-state Klinkenberg flow equation:

$$ 
\begin{aligned}
\varphi_Z &= -\frac{k_\infty}{\mu} \alpha (p + b) \frac{dp}{dz} \\
\frac{d\varphi_Z}{dz} &= 0
\end{aligned}
$$

This equation can then be cast in a simpler form by changing the unknown variable, using the "translated" pressure $P_b = p + b$ instead of $p$. this yields:

$$ 
\begin{aligned}
\varphi_Z &= -\frac{k_\infty}{\mu} \frac{\alpha}{2} \frac{d(P_B^2)}{dz} \\
\frac{d\varphi_Z}{dz} &= 0
\end{aligned}
$$

Finally, naming $\varphi_0 = \varphi_0$ the mass flux density at the inlet boundary (or equivalently at the outlet boundary), we obtain the simple equation in terms of the "translated" pressure $P_b = p + b$:

$$ -\frac{k_\infty}{\mu} \frac{\alpha}{2} \frac{d(P_B^2)}{dz} = \varphi_0 \quad \text{constant in space} \qquad (5.23) $$

Given the prescribed boundary pressure at $z = 0$ (bottom), this yields for $P_B(z)^2$:
$$P_B(z)^2 = -\frac{\varphi_0}{\frac{k_\infty}{\mu} \frac{\alpha}{2}} z + P_B(0)^2$$

Note that $P_B(z) = p(z) + b$ must satisfy $P_B > b \ge 0$. Therefore $P_B$ must be non negative, and one must reject the negative root. The resulting pressure profile $p(z)$ is then:

$$p(z) = -b + \sqrt{(p_{BOT} + b)^2 - \frac{\varphi_0}{\frac{k_\infty}{\mu} \frac{\alpha}{2}} z} \quad (5.24)$$

where we wrote $p(0) = p_{BOT}$ (bottom boundary pressure). Now, to obtain the mass flux density $\varphi_0$ [(kg/s)/m²], we apply the previous pressure profile at the top boundary, $z = L$, whence:

$$\varphi_0 = -\frac{k_\infty \alpha}{\mu 2} \frac{(p_{TOP} + b)^2 - (p_{BOT} + b)^2}{L} \quad (5.25)$$

Note that, in spite of the quadratic terms, the flux $\varphi_0$ depends linearly on the Klinkenberg coefficient "b". Indeed, the above expression can also be formulated equivalently as follows, using the "mean permeametric pressure" $\bar{p}$ defined below, and recalling that isothermal compressibility $\alpha$ is defined as $\alpha = \rho_0/p_0$ for any chosen reference state $(\rho_0, p_0)$:

$$\begin{gathered} \bar{p} = (p_{BOT} + p_{TOP})/2 \quad ; \quad \alpha = \bar{p}/\bar{p} \text{ with } \bar{p} = \rho(\bar{p}) \\ \varphi_0 = -\frac{k_\infty}{\mu} \alpha \frac{\bar{p}(1+b/\bar{p})}{L} \frac{p_{TOP} - p_{BOT}}{L} \quad (5.26) \\ \Rightarrow \varphi_0 = -\bar{p} \frac{k_\infty}{\mu} (1+b/\bar{p}) \frac{p_{TOP} - p_{BOT}}{L} \end{gathered}$$

The last equation in Eq. 5.26 is formulated such as to help interpretation of the flux in permeametric gas experiments in terms of two pressure quantities: (i) the pressure difference, and (ii) the mean pressure. This flux equation Eq. 5.26 is nevertheless exact, and equivalent to the previous Eq. 5.25. It is also similar to equation (3.5) in Wu et al. (1998), where these authors propose to evaluate both $k_\infty$ (Darcy permeability) and $b$ (Klinkenberg coefficient) by proceeding to a set of measurements with varying mean pressure $\bar{p}$ (e.g. varying the pressure at one end while the pressure remains fixed at the other end), measuring the flux $\varphi_0$, and plotting $\varphi_0/(p_{TOP} - p_{BOT})$ as a function of $\bar{p}$.

Indeed, from Eq. 5.26, here is the procedure for interpreting the permeametric gas flow experiment:

$$\begin{align} \text{a. Define } Y(\bar{p}) &= -\varphi_0 \mu \frac{1}{\alpha} \frac{L}{p_{TOP} - p_{BOT}} \tag{5.27} \\ \text{b. Plot the relation } Y(\bar{p}) \text{ experimentally.} \nonumber \\ \text{c. Calculate the linear regression: } \hat{Y}(\bar{p}) &= -k_\infty \bar{p} + k_\infty b \nonumber \end{align}$$

The slope of this linear regression will provide the product $k_\infty \bar{p}$ (whence $k_\infty$), and its intercept will provide $k_\infty b$ (whence the Klinkenberg coefficient $b$).

Finally, inserting the flux Eq. 5.25 or Eq. 5.26 into the pressure profile of Eq. 5.24 yields:
<figure><img src="image_4.png" /><figcaption>Fig. 1 Influence of the Klinkenberg coefficient b on the steady state mass flux density $\phi_0$ (kg/s)/m², with b ranging from 0 to 40 bars (i.e. $40 \times 10^5$ Pa). Note: the number Y = -1.397E-07 marked in the data cursor box at left is the value of $\phi_0$ for b = 0, in (kg/s)/m² units</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2 Influence of the Klinkenberg coefficient b on the steady state pressure profile, with b ranging from 0 to 40 bars (i.e. $40 \times 10^5$ Pa). The pressure profile for b = 0 is represented by a thicker curve</figcaption></figure>

$$p(z) = -b + \sqrt{(p_{\text{BOT}} + b)^2 + z \frac{(p_{\text{TOP}} + b)^2 - (p_{\text{BOT}} + b)^2}{L}} \quad (5.28)$$

As expected, these formulas reduce to the classical Darcy flow of compressible gases under zero gravity when the Klinkenberg coefficient vanishes. Indeed, inserting b = 0 in the previous steady equations, the reader can verify by inspection that these formulas, Eq. 5.25 or Eq. 5.26 for flux, and Eq. 5.28 for pressure, are identical with the equations obtained
previously for zero-gravity steady flow without the Klinkenberg effect (e.g. the expression of $p(z)$ given previously by Eq. 5.19).

The following figures describe the influence of the Klinkenberg coefficient "b" on steady state gas flow in the permeametric experiment with pressure conditions $p = 1 \text{ bar}$ at $z = 0$ and $p = 4 \text{ bars}$ at $z = 0.25 \text{ m}$. Fig. 1 shows the effect of "b" on the steady mass flux density $\phi_0$ (kg/s)/m² calculated exactly from Eq. 5.25 (these results were also obtained numerically with the COMSOL software and were indistinguishable from the exact ones, so they need not be both reproduced here). Fig. 2 shows the effect of "b" on the steady pressure profile $p(z)$ [Pa] calculated exactly from Eq. 5.28. In both figures, we let "b" vary from $b = 0$ to $b = 40 \text{ bars}$. We may focus in particular on the intermediate value $b = 4.258 \text{ bars}$ which could be relevant for the claystone's EDZ at the CIGEO site (Bure, Meuse/Haute-Marne, France), and the higher value $b = 15 \text{ bars}$ which could be relevant for the intact COx claystone at the same site (even higher values of "b" are also possible, as discussed earlier).

Note that the flux increases in absolute value as the Klinkenberg coefficient $b$ increases, and the pressure profiles become less curved (more linear) as $b$ increases. These effects can be explained by the fact that gas permeability $k_g(p)$ increases with $b$, whence the larger flux resulting under fixed pressure conditions. Also, as gas permeability increases (due to $b$), gas compressibility then plays a lesser role and the pressure profiles become more linear as if the gas were incompressible.

# 6 Analytical transient gas flow solutions

## 6.1 Analytical transient gas flow: introduction

Analytical solutions for transient gas flow with or without Klinkenberg effect have been developed, in particular, by Wu et al. (1998) for various flow configurations: their work, and other similar works, were reviewed earlier in the Literature Review section of this paper. Let us focus in particular on the analytical treatment of the 1D transient gas flow equation (neglecting gravity).

We will develop in this section several analytical approximations for transient 1D gas flow with or without Klinkenberg effect, neglecting gravitational acceleration. One of our approximations is "quadratic" (in that the squared pressure is the diffusion variable), analogous to some approximations developed earlier by Wu et al. (1998) for radial flow systems, although we treat here the gas flow problem in a vertical porous column. The other approximation is different, it can be viewed as a "linear" approximation (the pressure is the diffusion variable, not squared pressure). As will be seen, the two approximations perform differently at early times and large times. A third approximation combining the advantages of both is then proposed here, which is better than the existing approximations proposed in the literature.

*Note.* As will be seen, our solution method for solving the *linearized* transient gas flow problems is separation of variables, leading to an infinite eigenfunction expansion (solution of a Sturm-Liouville problem). Alternatively, other techniques like the Laplace transform could be used, perhaps with some difficulties in performing the inverse Laplace transform. At any rate, these techniques (separation of variables, Laplace transform) can only
be applied to linearized PDE's. One key point of our method for obtaining gas pressure profiles $p(z,t)$ is to optimally combine two types of linearizations.

## 6.2 Analytical approximations for transient gas flow

We start by re-considering the previous transient flow system of equations (Eq. 5.7), this time by considering the pressure dependent permeability Eq. 2.1

$$k(p) = k_{\infty} \left( 1 + \frac{b}{p} \right)$$

Inserting $k(p)$ in the Darcy mass flux density $\varphi(z, t)$ [(kg/s)/m²] of Eq. 5.2 yields:

$$\varphi(z,t) = - \frac{k_{\infty} \left(1 + \frac{b}{p}\right)}{\mu} \alpha p \left( \frac{dp}{dz} + \alpha g p \right) \qquad (6.1) \\ \Rightarrow \varphi(z,t) = - \frac{k_{\infty}}{\mu} \alpha (p + b) \left( \frac{dp}{dz} + \alpha g p \right)$$

The form of this Klinkenberg Darcy flux equation (Eq. 6.1) suggests the change of variable $p(z, t) \rightarrow P_B(z, t) = p(z, t) + b$, where $P_B$ can be dubbed the "Klinkenberg pressure" for short. With this, and neglecting gravity ($g = 0$), the flux equation becomes:

$$\varphi(z,t) = - \frac{k_{\infty}}{\mu} \frac{\alpha}{2} \frac{d(P_B^2)}{dz} \qquad (6.2)$$

Inserting this Klinkenberg Darcy flux equation in the mass conservation Eq. 5.4 for zero gravity and zero source term yields finally, in terms of the "Klinkenberg pressure" variable $P_B$, the following "pressure equation":

$$\alpha \Phi \frac{\partial P_B}{\partial t} - \frac{k_{\infty}}{\mu} \frac{\alpha}{2} \frac{\partial^2 (P_B^2)}{\partial z^2} = S(z, t) = 0 \qquad (6.3)$$

Taken together, Eqs. 6.2 and 6.3 form the nonlinear system of flux-pressure equations to be solved for gas flow in the presence of gas compressibility and Klinkenberg permeability. Note that they differ from the previous transient flow Eqs. 5.7 in two ways: first, we do take into account the Klinkenberg effect here; but secondly, we neglect gravity here (no centrifugal acceleration).

To solve this nonlinear transient flow system analytically, we propose at least two different linearization approximations: (i) direct linearization of the mass conservation equation in terms of the variable $P_B = p + b$, and (ii) "quadratic" linearization in terms of the variable $\Psi(z, t) = P_B^2 = (p + b)^2$.

The first linearization consists in approximating $d^2(P_B^2)/dz^2 \approx \bar{P}_B d^2(P_B)/dz^2$ or, alternatively, $d^2(P_B^2)/dz^2 \approx 2 \bar{P}_B d^2(P_B)/dz^2$. Thus:

$$\frac{d^2}{dz^2} (P_B^2) = \frac{d}{dz} \left( 2P_B \frac{dP_B}{dz} \right) \approx \frac{d}{dz} \left( 2\bar{P}_B \frac{dP_B}{dz} \right) \approx 2 \bar{P}_B \frac{d^2(P_B)}{dz^2}$$
where we take the mean Klinkenberg pressure to be the average between top and bottom pressures:

$$ \bar{P}_B = \bar{p} + b \approx \frac{p_{TOP} + p_{BOT}}{2} + b $$

Choosing this alternative for the direct linearization scheme (i) yields for mass conservation the linear diffusion equation (which can be written equivalently in terms of $P_B$ or $p$):

$$ \frac{\partial p}{\partial t} \approx D_1 \frac{\partial^2 p}{\partial z^2} \quad \text{with} \quad D_1 = \frac{1}{\Phi} \frac{k_\infty}{\mu} (\bar{p} + b) \qquad (6.4) $$

The corresponding linearized approximation scheme (i) for the transient flux is:

$$ \varphi(z, t) \approx -\frac{k_\infty}{\mu} \alpha (\bar{p} + b) \frac{dp}{dz} \qquad (6.5) $$

where the pressure $p(z, t)$ must be obtained by solving the pressure-based diffusion Eq. 6.4.

The second approximation scheme (ii) is "quadratic"; it is based on the "quadratic" pressure $\Psi$:

$$ \Psi = P_B^2 = (p + b)^2 \Rightarrow p = -b + \sqrt{\Psi} \qquad (6.6) $$

It consists in expressing all derivatives in terms of this quadratic pressure $\Psi$, and in linearizing the time derivative term as shown below. The $\Psi$-form of the mass conservation equation is the following exact nonlinear diffusion equation:

$$ \frac{\alpha}{2} \frac{\Phi}{\sqrt{\Psi}} \frac{\partial \Psi}{\partial t} = \frac{\alpha}{2} \frac{k_\infty}{\mu} \frac{\partial^2 \Psi}{\partial z^2} $$

The nonlinear term on the left-hand side, $\sqrt{\Psi}$, is then frozen to a constant mean value $\Psi \approx \bar{\Psi}$ defined as follows:

$$ \bar{\Psi} = \frac{\Psi_{BOT} + \Psi_{TOP}}{2} = \frac{(p_{BOT} + b)^2 + (p_{TOP} + b)^2}{2} \qquad (6.7) $$

The resulting mass conservation equation has the form of a linearized diffusion equation in terms of the $\Psi$-variable, as follows:

$$ \frac{\partial \Psi}{\partial t} \approx D_2 \frac{\partial^2 \Psi}{\partial z^2} \quad \text{with} \quad D_2 = \frac{1}{\Phi} \frac{k_\infty}{\mu} (\bar{p} + b) = \frac{1}{\Phi} \frac{k_\infty}{\mu} \sqrt{\bar{\Psi}} \qquad (6.8) $$

Note that the diffusion coefficient $D_2$ for this "quadratic" approximation is the same as for the previous "linear" approximation, and they can both be renamed "$D$". That is:

$$ D_2 = D_1 = D = \frac{1}{\Phi} \frac{k_\infty}{\mu} (\bar{p} + b) \qquad (6.9) $$

Finally, for the quadratic approximation scheme (ii), the transient mass flux density is given by:
$$ \varphi(z, t) \approx -\frac{k_{\infty}}{\mu} \frac{d(p^2)}{dz} \qquad (6.10) $$

This flux can be computed from the quadratic approximation of pressure $p(z, t)$ to be obtained by solving the $\Psi$-based diffusion Eq. 6.8 and using the relation between $\Psi$ and $p$ as defined by Eq. 6.6.

The analytical solutions of the two linearized diffusion equations Eqs. 6.5 and 6.10 are formally similar, although the results in terms of pressure profiles are quite different. We develop a formal infinite series solution for both cases under a spatially constant initial condition, with Dirichlet boundary conditions at $z = 0$ (bottom) and $z = L$ (top). The method is to decompose first the problem into a steady state component and a residual transient component which satisfies a diffusion problem with zero boundary conditions. The latter problem is then solved by the method of separation of variables (i.e., finding the eigenvalues and eigenfunction solutions of a Sturm-Liouville problem). The details are developed in APPENDIX B.

Thus, we reproduce here the final analytical results for both linearization schemes (i) and (ii) in terms of pressure (these will also lead to two distinct approximations (i) and (ii) for flux). Let $\Theta(x, t)$ represent the dependent variable for both pressure linearization schemes (i) and (ii). The infinite series solution obtained for pressure in APPENDIX B, can be summarized as follows:

$$ \Theta(x, t) = \theta(x, t) + \Theta_{\infty}(x) \qquad (6.11) $$

$$ \Theta_{\infty}(x) = \left( \frac{\Theta_2 - \Theta_1}{L} \right) x + \Theta_1 \qquad (6.12) $$

$$ \theta(x, t) = \sum_{n=1}^{\infty} \frac{2}{n\pi} \{\alpha + (\beta - \alpha)(-1)^n\} e^{-\lambda_n^2 Dt} \sin(\lambda_n x), \quad \lambda_n = \frac{n\pi}{L}, \quad n = 1, 2, \dots \qquad (6.13) $$

$$ \alpha = (\Theta_0 - \Theta_1); \quad \beta = (\Theta_2 - \Theta_1) \qquad (6.14) $$

The infinite series $\theta(x, t)$ must be truncated to a finite number of terms: we used 1001 terms in most cases. The "linear" (i) and the "quadratic" (ii) solutions for gas flow are then obtained as follows from Eqs. 6.11, 6.12, 6.13, and 6.14.

1. "Linear" analytical approximation (i): $p(x, t)^{(i)}$

The variable $\Theta(x, t)$ stands for pressure $p$. The space-time pressure $p(x, t)$ is obtained by applying Eqs. 6.11, 6.12, 6.13, and 6.14, with $\Theta(x, t)$ replaced everywhere by pressure $p(x, t)$. If $b = 0$ (no Klinkenberg effect), then the initial and boundary conditions named $\Theta_k$ ($k = 0, 1, 2$) should be replaced by the corresponding initial and boundary pressures $p_k$ ($k = 0, 1, 2$). Otherwise, if $b \neq 0$ (Klinkenberg effect), the initial and boundary conditions named $\Theta_k$ ($k = 0, 1, 2$) should be replaced by $p_k + b$ ($k = 0, 1, 2$).

2. "Quadratic" analytical approximation (ii): $p(x, t)^{(ii)}$

The variable $\Theta(x, t)$ stands for $\Psi = (p+b)^2 = P_B^2$. Eqs. 6.11, 6.12, 6.13, and 6.14, should be applied with $\Theta(x, t)$ replaced everywhere by $\Psi(x, t)$, which is defined in Eq. 6.6. Similarly, the initial and boundary conditions named $\Theta_k$ ($k = 0, 1, 2$) should be replaced by the corresponding initial and boundary values of $\Psi$, that is:
$\Psi_k = (p_k + b)^2$ ($k = 0, 1, 2$). After solution, the $\Psi(x, t)$ infinite series is used to obtain pressure $p(x, t)^{(ii)}$ by letting $p = -b + \sqrt{\Psi}$ according to Eq. 6.6. This procedure is applicable with or without Klinkenberg effect, that is, for $b \neq 0$ as well as for $b = 0$.

## 6.3 Permeametric gas flow parameters used for applications

The permeametric gas flow conditions and parameters are summarized below:

* Compressible hydrogen gas $H_2$, isothermal conditions.

* Gravitational acceleration: $g = 9.81 [m^2/s] \Rightarrow g \approx 0$ (negligible).

* Darcy permeability: $k_\infty = 5 \times 10^{-18} [m^2]$

* Klinkenberg coefficient: $b = 0$ [$Pa$]

* Isothermal gas compressibility ($H_2$): $\alpha = 8.271 \times 10^{-7}$ [($kg/m^3$/Pa)]

* Dynamic gas viscosity ($H_2$): $\mu = 8.88 \times 10^{-6}$ [$Pa \cdot s$]

* Porosity: $\Phi = 0.40 [m^3/m^3]$

* Domain length: $L = 0.25$ [$m$]

* Boundary conditions: $x = 0 : p = p_1 = 1 \times 10^5$ [$Pa$]; $x = L : p = p_2 = 4 \times 10^5$ [$Pa$]; Mean boundary pressure: $\bar{p} = (p_1 + p_2)/2 = 2.5 \times 10^5$ [$Pa$]

* Initial condition (constant pressure profile): $t = 0 : x \in [0, L]$ $p(x, 0) = p_0 = 1.5 \times 10^5$ [$Pa$].

<figure><img src="image_6.png" /><figcaption>**Fig. 3** Analytical pressure profiles $p(x, t)$ (or $p(z, t)$ with negligible gravity) based on the "Linear" approximation option (i). Input data for this permeametric gas flow "experiment" are given in the text. Output times are $t = 0$ (initial condition), and $t = 1$ s, 25 s, 100 s, 1000 s, 5000 s, 30 000 s, and 86 400 s. The dashed curve represents the exact steady state profile, which is different from the large time profile obtained with this "Linear" approximation</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 4 Analytical pressure profiles $p(x, t)$ (or $p(z, t)$ with negligible gravity) based on the "Quadratic" approximation. Input data for this permeametric gas flow "experiment" are described in the text. Output times are $t = 0$ (initial condition), and $t = 1$ s, 25 s, 100 s, 1000 s, 5000 s, 30 000 s, and 86 400 s</figcaption></figure>

## 6.4 Linear & Quadratic approximations of pressure $p(z, t)$ compared with numerical simulations

The analytical $H_2$ gas pressure profiles $p(x, t)$ for the case $b = 0$ (no Klinkenberg effect) are shown in Fig. 3 for the "Linear" option (i), and Fig. 4 for the "Quadratic" option (ii). The infinite series were truncated to 1001 terms in both cases. These analytical pressure profiles will be combined and then compared to the corresponding numerical experiment simulation of permeametric gas flow developed with the Finite Volume TOUGH2 code for the same permeametric gas flow conditions: see pressure profiles shown later in Fig. 5b. Comparing the three pressure profiles $p(x, t)$ (Fig. 3, Fig. 4, Fig. 5b) indicates that the "Linear" approximation (i) is better at early times but becomes inaccurate at large times and yields the "wrong" steady state, while the "Quadratic" approximation (ii) is better at later times and becomes exact at steady state. Based on these preliminary observations, we will then develop, next, an improved analytical solution, before re-considering analytical vs. numerical comparisons.

In summary, comparing the two analytical approximations of $p(x, t)$ with the numerical simulation, it can be clearly seen that the pressure profiles of the first ("Linear") approximation are better at early times but become worse as time grows and tend to the wrong linear pressure profile at large times and at steady state. On the contrary, the pressure profiles of the second ("Quadratic") approximation are inaccurate at early times but become better as time grows and tend to the exact quadratic pressure profile at large times and at steady state. A combination procedure (logistic time average) will be proposed shortly in order to produce a single analytical approximation of $p(x, t)$ that is better than any of the previous two approximations.

In terms of mass flux density $\varphi(x, t)$, the numerical solution obtained with the TOUGH2 code yielded the following value at large times (steady state $\varphi_0$ constant in space): $\varphi_0 = -1.3971 \times 10^{-7}$ [($kg/s$)/$m^2$]. In comparison, the "Quadratic" approximation (ii),
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Comparing analytical (top) and numerical (bottom) gas pressure profiles. The pressure profiles simulated by the TOUGH2 Finite Volumes code compare favorably with the combined linear-quadratic analytical solution</figcaption></figure>

being exact at steady state, yielded the same value of $\varphi_0$ to within the 3rd or 4th digit. More details will be provided shortly concerning the full transient flux solutions $\varphi(z, t)$.

**6.5 Improved analytical pressure profiles for transient gas flow**

We develop here a new procedure to produce a single analytical approximation of $p(x, t)$ that is better than any of the previous two infinite series approximations. The method is to combine the two solutions using a logistically weighted temporal average of the first

<figure><img src="image_10.png" /></figure>
solution (which is good at early times) and of the second solution (which is good at large times and becomes exact as $t \to \infty$). The method requires, first, an estimation of the time scale that separates "early times" from "large times".

### Time scale $T_{SCALE}$ separating early times from large times

Let us designate as $T_{STEADY}$ the diffusion time scale that indicates the approach of a steady state for a diffusion process with coefficient $D$ [$m^2/s$] occurring in a column of length $L$. Note that both the "Linear" and the "Quadratic" diffusion approximations have the same "linearized" diffusion coefficient $D$. The diffusion time scale $T_{STEADY}$ can then be evaluated from:

$$T_{STEADY} = \frac{L^2}{2D} \qquad (6.15)$$

The diffusion coefficient was defined for both analytical approximations by Eq. 6.9. Applying the numerical values of the permeametric gas flow parameters given just above yields:

$$D = \frac{1}{0.40} \frac{5 \times 10^{-18}}{8.88 \times 10^{-6}} (2.5 \times 10^5 + 0) = 3.5191 \times 10^{-7} \quad [m^2/s]$$

As a consequence, the steady state time scale for gas flow is evaluated as:

$$T_{STEADY} = \frac{L^2}{2D} = \frac{(0.25)^2}{2 \times 3.5191 \times 10^{-7}} = 88800 \quad [s]$$

Examining the previous pressure profiles, analytical & numerical, it can be seen that indeed $T_{STEADY} = 88800$ $s$ is a fair approximation of the steady state time scale (about 1 day). Now, in order to obtain a separation time scale ($T_{SCALE}$) between "early" times and "large" times, we use a quarter of the steady time value, that is:

$$T_{SCALE} \approx \frac{T_{STEADY}}{4} = \frac{L^2}{8D} = 22200 \quad [s]$$

### Logistic temporal weighting of "Linear" and "Quadratic" approximations $p(z, t)$

Given the value of $T_{SCALE}$, the idea is to improve the two analytical approximations:

(i) Linear, and

(ii) Quadratic,

by combining the corresponding pressure solutions $p_{(i)}(x, t)$ and $p_{(ii)}(x, t)$ as follows:

$$\hat{p}(x, t) = \frac{p_{(ii)}(x, t)}{1 + \left( \frac{p_{(ii)}(x,t)}{p_{(i)}(x,t)} - 1 \right) \exp(-t/T_{SCALE})} \qquad (6.16)$$

The above improved solution $\hat{p}(x, t)$ can be interpreted as a continuous time weighted average of $p_{(i)}(x, t)$ ("Linear") and $p_{(ii)}(x, t)$ ("Quadratic") based on a logistic curve, using the pre-calculated time scale $T_{SCALE}$ separating early times from large times. This will improve significantly the analytical approximation of the linearized gas flow equation.

## 6.6 Comparing analytical vs. numerical pressure profiles $p(z,t)$

The (combined linear-quadratic) analytical pressure profiles $p(z, t)$ of Fig. 5a compare well with the numerical pressure profiles obtained with the TOUGH2 Finite Volume code for
<figure><img src="image_11.png" /><figcaption>Fig. 7 Analytical mass flux density $\varphi(x, t)$ [or $\varphi(z, t)$ with negligible gravity], based on applying Darcy's law to the "Quadratic" approximation of pressure $p(x, t)$. Input data for this permeametric gas flow "experiment" are given in the text. Output times are $t = 0$ (initial condition), and $t = 1$ s, 25 s, 100 s, 1000 s, 5000 s, 30 000 s, and 86 400 s. The number $X$ given in the box shows the value obtained for the constant profile $\varphi(z, t)$ at $t = 86 400$ s, which coincides with the exact steady state</figcaption></figure>

$$ \varphi(z, t) \approx \varphi(z, t)^{(2)} = -\frac{k_{\infty}}{\mu} \frac{\alpha}{2} \left[ \frac{P_2^2 - P_1^2}{L} + \frac{2}{L} \sum_{n=1}^{\infty} \{\alpha_2 + (\beta_2 - \alpha_2) (-1)^n\} \cos(\lambda_n x) e^{-\lambda_n^2 D t} \right] \quad (6.17) $$

where:

$$ D = \frac{1}{\Phi} \frac{k_{\infty}}{\mu} (\bar{p} + b) : \text{Diffusion coefficient } [m^2/s] $$

$$ \lambda_n = n\pi/L \quad (n = 1, 2, \dots) : \text{Eigenvalues of diffusion PDE} $$

$$ P_0 = p_0 + b, \quad \forall x \in [0, L] : \text{Initial pressure} + "b" $$

$$ P_1 = p(0) + b, \quad P_2 = p(L) + b : \text{Boundary pressures} + "b" $$

$$ \bar{P} = (P_1 + P_2)/2 = \bar{p} + b : \text{"Mean" pressure} + "b" $$

$$ \alpha_2 = P_0^2 - P_1^2 = (p_0 + b)^2 - (p_1 + b)^2 $$

$$ \beta_2 = P_2^2 - P_1^2 = (p_2 + b)^2 - (p_1 + b)^2 $$

More details on analytical mass flux are provided in APPENDIX C, including another approximation $\varphi(z, t)^{(1)}$ which if formally similar, albeit distinct from, Eq. 6.17. Although the transient flux approximations $\varphi(x, t)^{(1)}$ and $\varphi(x, t)^{(2)}$ are distinct in space and time, it was found in APPENDIX C, that the resulting steady state fluxes $\varphi(x, \infty)^{(1)}$ and $\varphi(x, \infty)^{(2)}$ are identical uniform profiles (independent of $z$), and that both are exact and given by:

$$ \varphi(x, \infty) = \varphi_0 = -\frac{k_{\infty}}{\mu} \alpha (\bar{p} + b) \left[ \frac{p_2 - p_1}{L} \right], \quad \forall x \in [0, L] \quad (6.18) $$

For the input data at hand, the corresponding flux value is:
<figure><img src="image_12.png" /><figcaption>Fig. 8 Mass flux density profiles $\varphi(x, t)$ [or $\varphi(z, t)$], comparing analytical flux profiles (solid curves) with Finite Elements simulations implemented with COMSOL Multiphysics (dotted curves). Input data for this permeametric "gas flow experiment" are the same as those given in the text. The initial flux condition is $\varphi(z, 0) = 0, \forall z \in [0, L]$. Output times are: $t = 1$ s for the analytical solution, and $t = 25$ s, 100 s, 1000 s, 5000 s, 30 000 s, 86 400 s for both analytical and COMSOL solutions</figcaption></figure>

$$t \to \infty : \varphi(z, t) \to \varphi_0 \approx -1.383 \ 10^{-7} (kg/s)/m^2 \quad \forall z \in [0, L]$$

By extension, the convergence of the transient flux solutions (1) and (2) towards the same exact steady state suggests that both are reasonably accurate and close to each other also at sufficiently large times. In the sequel, we focus on the quadratic approximation $\varphi(x, t)^{(2)}$ given by Eq. 6.17 above (or Eq. C.4 in APPENDIX C) in order to analyze the transient behavior of mass flux density profiles.

The next figure, Fig. 7, illustrates the space-time behavior of the analytical $\varphi(z, t)$. It can be seen in both cases that the constant steady state gas flux given above ($\varphi_0 \approx -1.383 \times 10^{-7} (kg/s)/m^2$) is reached at large times over the entire porous column (after one day or less). Note also that, at mid-point position, the flux reaches its steady state value as soon as $t \approx 30 000$ s, while at other points of the column the steady state is only reached at larger times $t \approx 86 400$ s (about one day).

Furthermore, the space-time behavior of the analytical flux $\varphi(z, t)$ is validated by a direct numerical solution of the compressible gas flow equations using the Finite Element software platform COMSOL Multiphysics: see Fig. 8, where the analytical solution of Fig. 7 is reprinted together with the numerical flux profiles computed with COMSOL Multiphysics Finite Elements software using mass conservation and Darcy's law for compressible $H_2$ gas (i.e., the PDE's presented earlier). It can be seen that the computed vs. analytical flux profiles $\varphi(z, t)$ are close, and become closer as time grows towards the exact vertical steady state flux profile.

## 6.8 Relative error analyses: analytical-numerical discrepancies

Finally, the discrepancy between analytical and numerical solutions is calculated at each time $t$ in the form of a Mean Relative Error ($\epsilon_{MRE}$). We focus here on the numerical simulation obtained from the COMSOL Multiphysics Finite Element software, and we consider both the
<figure><img src="image_13.png" /><figcaption>Fig. 9 The results of error analysis are plotted here versus time for mass flux density (red curve: $\epsilon_{MRE}(\varphi; t)$) and for pressure (blue curve: $\epsilon_{MRE}(p; t)$). In each case, the ordinate represents the dimensionless Mean Relative Error "MRE" between the approximate analytical solution and the COMSOL numerical simulation, and the abscissa represents time from $t = 0$ to $t \approx 86400$ s (1 day).</figcaption></figure>

pressure variable $p(z, t)$ and the mass flux density $\varphi(z, t)$. The MRE is defined as follows for these two variables:

$$ \epsilon_{MRE}(p; t) = \frac{1}{N_Z} \sum_{j=1}^{j=N_Z} \frac{|p(z_j, t)_{NUMER} - p(z_j, t)_{ANALYT}|}{p_0} \quad (6.19) $$

$$ \epsilon_{MRE}(\varphi; t) = \frac{1}{N_Z} \sum_{j=1}^{j=N_Z} \frac{|\varphi(z_j, t)_{NUMER} - \varphi(z_j, t)_{ANALYT}|}{|\varphi_0|} \quad (6.20) $$

where $N_Z$ is the number of nodes along the porous column in the numerical and analytical solutions: here the length of the column is $L = 0.25$ m and we have taken $N_Z = 100$ nodes. The symbol "ANALYT" designates the analytical solution, and the symbol "NUMER" designates the numerical result (COMSOL simulation).<fnref n="2" /> The parameter $p_0$ designates the

<footnote n="2">² Specifications of COMSOL simulations: Finite Elements with square elements; 2D column of width $L_X = 0.04$ m and height $L_Z = 0.25$ m; the 2D boundary conditions were such that flow is 1D vertical; the square mesh elements were of size $\Delta x = \Delta z = 0.0025$ m; number of elements in each direction: $N_X = 16$ and $N_Z = 100$ (1600 elements in 2D); time stepping: there were 100 output time steps from $t = 0$ to $t = 86400$ s; the numerical time steps were adaptive, variable, and physics controlled.</footnote>
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Log-Log plot of the temporal evolution of mass flux density $|\varphi(z_{MID}, t)|$ at the mid-point position $z_{MID} = L/2$, computed numerically with the Finite Volume code TOUGH2, from time $t = 100$ s until $t = 1\,728\,000$ s, for different values of gravitational acceleration: $g^* = 1 \times g$ (natural gravity), $g^* = 10\,000 \times g$, $g^* = 20\,000 \times g$, and $g^* = 100\,000 \times g$ (strong centrifugal acceleration)</figcaption></figure>

uniform initial pressure condition $p_0 = 1.5 \times 10^5$ [Pa] (see permeametric gas flow parameters in Sect. 6.3), and $\varphi_0 \approx -1.383 \times 10^{-7}$ $(kg/s)/m^2$ designates the final steady state mass flux density. Note that the numerical and analytical values of $\varphi_0$ are indistinguishable; the analytical formula of $\varphi_0$ was given earlier in Eq. 6.18 or equivalently in Eq. 5.26 for any value of Klinkenberg's coefficient "b" (we have $b = 0$ in the present error analysis).

The results of error analysis are plotted in Fig. 9, where both $\epsilon_{MRE}(p;t)$ and $\epsilon_{MRE}(\varphi;t)$ are plotted versus time. The mass flux density error ($\epsilon_{MRE}(\varphi;t)$) is generally greater than the pressure error ($\epsilon_{MRE}(p;t)$), but both remain reasonably low except at very short times. The flux error is about 50% at $t = 1$ s (not visible in the plot), but reduces to 17% at $t = 25$ s (visible in the plot), and decreases drastically at later times to reach finally 0.2% at $t = 86400$ s = 1 day. The pressure error is small at short times, reaches a peak of about 5% at $t = 5000$ s, and then decreases down to 0.5% at $t = 86400$ s = 1 day. These good results can be interpreted as a cross-validation between the finite element simulation and our approximate analytical solution.
## 6.9 Transient effects of strong centrifugal acceleration fields (super-gravity)

It was demonstrated earlier that natural gravity is negligible for gas flow, but the case of large centrifugal acceleration $g*$ is different. When $g \gg g$ replaces the natural acceleration $g$ in the equations for compressible gas flow in a porous column, the acceleration may no longer be negligible. Our analytical solution, developed with the acceleration term in APPENDIX A for the steady state case, indicates that the steady pressure profile depends on compressibility $\alpha$ if centrifugal acceleration $g*$ is on the order of tens of thousands of $g$'s or more.

However, there is in fact a gravitational effect especially during early stages of transient gas flow, and the gas flux $\varphi(z, t)$ is also affected in addition to pressure profiles $p(z, t)$. Numerical simulations of transient gas flow indicate that the gravitational effect becomes visible even at moderate accelerations $g \approx g_0$ at sufficiently early times.

We show only one example here: Fig. 10 shows a sensitivity analysis of the time evolution of mass flux density $\varphi(z_{\text{MID}}, t)$ (kg/s)/m$^2$ (in absolute value) at the mid-point section of the column $z_{\text{MID}} = 0.125$ m, for different values of gravitational acceleration: $g*=1$ g, 10 000 g, 20 000 g, and 100 000 g. The simulated $\varphi(z_{\text{MID}}, t)$ appears quite sensitive to $g*$ at early times, although not very sensitive at large times corresponding to the asymptotic steady state.

The steady value of $\varphi(z_{\text{MID}}, t)$ for $g*=100\,000\,g$ in Fig. 10 is $\varphi_0 \approx 2 \times 10^{-7}$ [($kg/s/m^2$)], which is only a bit larger than the value $\varphi_0 \approx 1.4 \times 10^{-7}$ [($kg/s/m^2$)] obtained for natural gravity $g*=1 \times g$. Note that the steady flux $\varphi_0 \approx 1.4 \times 10^{-7}$ obtained numerically for $g*=1 \times g$ is the same as that predicted theoretically for negligible gravity from Eqs. 5.18 and 5.20. This validates our large time numerical flux, at least for the case of natural gravity $g$, and it also confirms that natural gravity $g$ is negligible and equivalent to 0 (as could be expected from our previous analyses of the gravitational group).

Other simulation results (not shown here) indicate that the pressure profiles are only moderately sensitive to $g*$. For instance, considering large times (steady state) pressure profiles $p(z)$, we have $p = 2.5$ bars at elevation $z \approx 9.0$ cm for negligible gravity ($g*=1$ g $\approx$ 0), compared to $z \approx 7.5$ cm for super-gravity $g*=100\,000 \times g$. In summary, the variable most affected by super-gravity is the mass flux density $\varphi(z, t)$ at early or intermediate times.

## 7 Recapitulation of analytical solutions

Finally, this section presents a brief recapitulation of analytical solutions (exact or approximate) discussed and developed in this paper on compressible gas flow in porous media: these are summarized in the list presented below.

1. Steady state gas pressure profile $p(z)$ and mass flux density $\varphi_0$ with or without Klinkenberg effect ($b \neq 0$ or $b = 0$), for negligible gravitational acceleration $g \approx 0$ (natural gravity $g = 9.81 \, m/s^2$ being negligible for any reasonable column length).

- For $b=0$, exact solution $p(z)$ is given by Eq. 5.19; and exact solution for mass flux density $\varphi_0$ is given by Eq. 5.18.

- For $b \neq 0$, exact solution $p(z)$ is given by Eq. 5.28 or Eq. 5.24; and exact solution for mass flux density $\varphi_0$ is given by Eq. 5.25 or equivalently by Eq. 5.26.
These solutions are also available in different forms, more or less explicitly, in several of the references cited in this text.

In addition, the Klinkenberg effect is studied in Fig. 1, which shows the influence of the Klinkenberg coefficient *b* on the steady mass flux density φ₀ (with *b* ranging from 0 to 40 *bars*), and in Fig. 2 showing the steady pressure profiles *p*(*z*) for several values *b* (ranging from 0 to 40 *bars*).

2. Steady state gas pressure profile and mass flux density for $b = 0$ and for strong non negligible gravitational acceleration ("super-gravity"), for instance in the case of centrifugal acceleration ≈ tens of thousands *g*.

- Exact pressure profile $p(z)$ is given by Eq. 5.16 together with Eq. 5.17.

- Exact mass flux density $\varphi_0$ given by Eq. 5.15.

This solution is not available in the literature as far as we know.

3. Transient space-time analytical solutions for gas pressure profiles and mass flux density profiles.

- Two analytical approximations are provided in the form of infinite series for either pressure ("linear" approximation) or quadratic pressure ("quadratic" approximation), both given in Eqs. 6.11, 6.12, 6.13, and 6.14, as explained in the text below these equations. In practice the infinite series were truncated to about 1000 terms in our examples of permeametric gas flow. See Fig. 3 for the "linear" approximation, and Fig. 4 for the "quadratic" approximation. The analytical expressions are valid for $b = 0$ as well as $b \neq 0$ (Klinkenberg).

- A better approximation of transient pressure profiles is then obtained by an optimal weighting of the two previous analytical solutions based on a time scale $T_{SCALE}$ separating short times from large times. This scheme yields the analytical solution described by Eq. 6.16. This analytical expression is also valid for $b = 0$ as well as $b \neq 0$ (Klinkenberg). The resulting pressure profiles are depicted in Fig. 5a.

- Mass flux density $\varphi(z,t)$ is then obtained analytically as an infinite series (Eq. 6.17) by applying Darcy's law for compressible gas to the quadratic pressure approximation. The analytical expression is valid for $b = 0$ as well as $b \neq 0$ (Klinkenberg). The analytical space-time profiles of mass flux density $\varphi(z,t)$ are depicted in Fig. 7. The exact steady state flux $\varphi_0$ ($kg/s/m^2$) is obtained from the analytical transient: see Eq. 6.18.

4. Other findings related to Klinkenberg permeability: estimation of Klinkenberg coefficient *b* [*Pa*] by various methods (molecular gas kinetics, empirical regressions, etc.).

- Eq. 4.1 defines the classical pressure-dependent Klinkenberg permeability.

- Eq. 4.5 and several equations below present molecular gas kinetics estimation of "b", leading to a tentative estimate for $H_2$ in the COx claystone at the CIGEO radioactive waste site in Bure (Meuse-Haute Marne, France).

- Eq. 4.8: empirical regression $b(k_\infty)$ from the literature (where $k_\infty$ denotes the classical Darcy permeability without Klinkenberg effect); application: estimate of "b" for $H_2$ in the COx claystone and its EDZ.

Most of the parameters and variables appearing in the itemized list above were previously defined in this text. However, for the sake of completeness, here is a summary of these parameters and variables:
* $p$ [$Pa$] is pore gas pressure;

* $\alpha$ [$\frac{kg}{m^3}/Pa$] is gas compressibility, as defined earlier in Eq. 3.4 or Eq. 3.6, from the linear $\rho(p)$ relation for ideal isothermal gases;

* $\varphi(z, t)$ or $\varphi_0$ [$\frac{kg}{s}/m^2$] is the mass flux density of the gas in the porous medium;

* $\mu$ [$Pa \cdot s$] is the dynamic viscosity of the gas;

* $k$ [$m^2$] is the Darcy permeability of the porous medium; when discussing the Klinkenberg effect for low gas pressure, the gas permeability becomes the Klinkenberg permeability $k_G(p)$ given by Eq. 4.1, and the classical Darcy permeability is then denoted $k_\infty$ instead of $k$;

* $D$ [$m^2/s$] is the diffusion coefficient of the transient single phase gas flow, defined by casting the transient flow equation in the form of a (linearized) pressure or quadratic pressure diffusion equation. See the original non linear diffusion equation of gas flow (Eq. 6.3), and the resulting linearized version of the diffusion coefficient (Eq. 6.9). Remarkably, this coefficient $D$ is the same for the two distinct "Linear" and "Quadratic" approximations of pressure and flux.

*Note: we use the notations $p, \varphi, \alpha, D, k, k_\infty$ (etc.) in the above list and everywhere in this paper, instead of other notations that may have been used in the cited references.*

## 8 Conclusion

In conclusion, we have reviewed in this paper several aspects of single phase gas flow in porous media, taking into account ideal gas compressibility (isothermal), molecular gas slip-page (the Knudsen effect, with pressure-dependent Klinkenberg permeability), and steady state as well as transient solutions developed in the literature. We have also presented various methods for estimating the "b" coefficient in the pressure-dependent Klinkenberg permeability. Furthermore, we have developed several new analytical solutions focusing on the case of "permeametric" gas flow in a vertical porous column with pressures imposed at both ends, including:

(a) A new analytical solution for steady flow taking into account strong gravity (e.g. centrifugal acceleration); and

(b) Analytical solutions for transient gas flow taking into account the pressure-dependent Klinkenberg permeability.

Some of these analytical solutions are benchmarked with numerical simulations using either Finite Volumes (TOUGH2 code) or Finite Elements (COMSOL Multiphysics). Such analytical benchmarks are very useful for the development and validation of numerical codes.

Concerning topic (a), it is interesting to note that the steady vertical pressure profile $p(z)$ depends on gas compressibility $\alpha$ only for strong gravity, and only in that case, it is possible to devise an experiment to retrieve both the permeability and the compressibility from a steady state permeametric experiment (otherwise, if $g$ is negligible, $p(z)$ does not depend on $\alpha$).

Concerning topic (b), the proposed analytical approximation for transient pressure and flux profiles $p(z, t)$ and $\varphi(z, t)$ in the porous column is novel. It combines two different approximations of the gas flow equation based on a pre-calculated time scale ($T_{SCALE}$). This novel solution combines two optimal features: that of the "Linear" approximation which is preferred at early times, and that of the "Quadratic" approximation which is preferred at larger times.
It is envisioned that this solution technique could also be applied to other transient gas flow configurations, such as transient horizontal gas flow towards a vertical well, in terms of r, t coordinates with $r = \sqrt{x^2 + y^2}$.

In addition, future work could also focus on several effects relevant to gas flow in porous media, including for instance the case of adiabatic compressible gas flow (where $\rho(p)$ is a power law instead of linear), and the case of departure from Darcy's law due to high Reynolds numbers (Darcy-Forchheimer linear-quadratic velocity law), in addition to the Klinkenberg effect due to gas slippage in tight pores (already taken into account in the present work). The analytical approaches outlined in this paper could also serve for considering two-phase liquid-gas flow, for instance in asymptotic cases where the liquid saturation is low, or possibly the reverse case where gas saturation is low, with applications to liquid flow and gas migration in underground radioactive waste repositories and in hydrocarbon oil-gas reservoirs.

## A. APPENDIX A: Steady gas flow solution with acceleration

This APPENDIX presents the detailed calculations that lead to the complete 1D steady state gas flow solution of Sect. 5.5 accounting for both ideal gas compressibility, and a strong acceleration field $g$ [$m^2/s$] (e.g. centrifugal acceleration instead of natural gravity). The strong acceleration was named "g *" in parts of the text discussing centrifugal acceleration effects; it is renamed "g" here for simplicity.

We start with the steady state gas flow Eq. 5.13 of Sect. 5.5, which we repeat in this APPENDIX for convenience: The steady-state version of Eq. 5.7 is (after factoring out $\alpha_0$ in the first equation):

$$ 
\begin{aligned}
0 &= \frac{k}{\mu} \left[ \frac{1}{2} \frac{\partial^2}{\partial z^2} (p^2) + \frac{\alpha g L}{L} \frac{\partial}{\partial z} (p^2) \right] \\
\varphi_Z &= \rho(p) q_Z = -\frac{k}{\mu} \alpha p \left( \frac{dp}{dz} + \frac{\alpha g L}{L} p \right)
\end{aligned} 
\quad (A.1) $$

where $z$ is upwards, and $g$ is the acceleration field oriented downwards (earth gravity is $g = 9.81 \ m^2/s$ but we also consider the possibility of centrifugal accelerations which are ten thousand times this value).

To solve Eq. A.1, we perform a change of variable, introducing the squared pressure $P$:

$$ p(z) \rightarrow P(z) = p^2(z) $$

This yields:

$$ 
\begin{aligned}
0 &= \frac{k}{\mu} \frac{1}{2} \frac{d}{dz} \left[ \frac{dP}{dz} + 2 \frac{\alpha_0 g L}{L} P \right] \\
\varphi_Z &= -\frac{k}{\mu} \frac{\alpha}{2} \left( \frac{dP}{dz} + 2 \frac{\alpha g L}{L} P \right)
\end{aligned} 
\quad (A.2) $$

where we recall that $P = p^2$. For verification, let us perform a brief dimensional analysis of the first equation in Eq. A.2 in terms of the squared pressure variable "P":
$$ \left[ \frac{dP}{dz} \right] = \left[ \frac{Pa^2}{m} \right], \text{ and } \left[ 2 \frac{\alpha g L}{L} P \right] = \left[ \frac{Pa^2}{m} \right], \text{ as it should.} $$

The steady state mass balance flow equation is the first equation in Eq. A.2, expressed in terms of squared pressure $P$. It is completed by the second equation in Eq. A.2 which expresses Darcy's mass flux density in terms of the squared pressure $P = p^2$. It is clear that the steady state mass conservation equation in Eq. A.2 expresses zero mass flux divergence:

$$ 0 = \frac{k}{\mu} \frac{\partial \varphi_z}{\partial z} \rightarrow \varphi_z = \varphi_0 (\text{constant mass flux } \forall z) \qquad (A.3) $$

In particular, with pressure conditions at both boundaries (top and bottom), this implies:

$$ \varphi_{TOP} = \varphi_{BOTTOM} = \varphi_0 [(kg/s)/m^2] $$

which is independent of $z$ in the porous column. Now, integrating the first equation in Eq. A.2 leads to stating that the Darcy mass flux is spatially constant in the column (constant $\varphi_0$). The equation to be solved is now the 1st order ODE:

$$ -\frac{k}{\mu} \frac{\alpha}{2} \left( \frac{dP}{dz} + 2 \frac{\alpha g L}{L} P \right) = \varphi_0 (\text{constant } \forall z) \qquad (A.4) $$

At this stage, a second change of variable is proposed in order to solve the first order differential Eq. A.4 in terms of the squared pressure variable $P = p^2$:

$$ P(z) = W(z) e^{-2\alpha g z} \qquad (A.5) $$

Inserting $P(z)$ of Eq. A.5 into the constant mass flux Eq. A.4 gives (with $W' = \frac{dW}{dz}$):

$$ \begin{aligned} & -\frac{k}{\mu} \frac{\alpha}{2} \left[ \frac{\partial}{\partial z} [W(z) e^{-2\alpha g z}] + 2\alpha g [W(z) e^{-2\alpha g z}] \right] = \varphi_0 \\ & -\frac{k}{\mu} \frac{\alpha}{2} \left[ [W' e^{-2\alpha g z} - W2\alpha g e^{-2\alpha g z}] + 2\alpha g [W(z) e^{-2\alpha g z}] \right] = \varphi_0 \\ & -\frac{k}{\mu} \frac{\alpha}{2} \left[ e^{-2\alpha g z} \frac{\partial W}{\partial z} \right] = \varphi_0 \\ & \implies \frac{\partial W}{\partial z} = -\frac{\varphi_0}{k/\mu} \frac{2}{\alpha} e^{+2\alpha g z} \end{aligned} \qquad (A.6) $$

Integrating Eq. A.6 between the two pressure boundaries $z_{BOT}$ and $z_{TOP}$ yields:

$$ W(z) = W(z_{BOT}) - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} [e^{+2\alpha g z} - e^{+2\alpha g z_{BOT}}] \qquad (A.7) $$

Applying the inverse transform $P(z) = W(z) e^{-2\alpha g z}$ gives the squared pressure $P(z)$:

$$ P(z) = W(z_{BOT}) e^{-2\alpha g z} - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} \left[ 1 - e^{-2\alpha g (z - z_{BOT})} \right] \qquad (A.8) $$

Applying this at $z = z_{BOT}$ gives immediately:

$$ P(z_{BOT}) = W(z_{BOT}) e^{-2\alpha g z_{BOT}} \implies W(z_{BOT}) = P(z_{BOT}) e^{+2\alpha g z_{BOT}} \qquad (A.9) $$
So the squared pressure profile solution $P(z)$ is:

$$P(z) = P(z_{BOT})e^{-2\alpha g(z-z_{BOT})} - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} \left[1 - e^{-2\alpha g(z-z_{BOT})}\right] \quad (A.10)$$

Finally, we back transform the squared pressure $P(z) = p^2(z)$ to obtain the pressure profile $p(z)$. Since pressure $p(z)$ stands for absolute gas pressure, it must be always positive ($p(z) > 0$); therefore:

$$P(z) = p^2(z) \rightarrow p(z) = +\sqrt{P(z)}$$

We also have to account for the pressure boundary conditions at $z = z_{BOT}$ and $z = z_{TOP}$:

$$P(z_{BOT}) = p^2(z_{BOT}); \quad P(z_{TOP}) = p^2(z_{TOP}).$$

So finally, the steady state pressure profile $p(z)$ is given by:

$$p(z) = \sqrt{p^2(z_{BOT}) e^{-2\alpha g(z-z_{BOT})} - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} \left[1 - e^{-2\alpha g(z-z_{BOT})}\right]} \quad (A.11)$$

where $p(z)$ is expressed in terms of the bottom pressure and the steady mass flux. An alternative expression can be obtained by applying Eq. A.11 at the top pressure boundary, and taking into account $L = z_{TOP} - z_{BOT}$:

$$p(z_{TOP}) = \sqrt{p^2(z_{BOT}) e^{-2\alpha gL} - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} \{1 - e^{-2\alpha gL}\}}$$

Squaring this expression on both sides yields the constant mass flux density $\varphi_0 [(kg/s)/m^2]$ in terms of both top and bottom pressure, as follows:

$$p^2(z_{TOP}) = p^2(z_{BOT}) e^{-2\alpha gL} - \frac{\varphi_0}{(k/\mu)} \frac{1}{\alpha^2 g} \left[1 - e^{-2\alpha gL}\right]$$

Whence, re-arranging:

$$\varphi_0 = -\frac{k}{\mu} \alpha^2 g \frac{p^2(z_{TOP}) - p^2(z_{BOT}) e^{-2\alpha gL}}{1 - e^{-2\alpha gL}} \quad (A.12)$$

Finally, inserting Eq. A.12 in Eq. A.11 gives a closed form expression for the pressure profile in terms of the top and bottom pressures. This expression, shown below, depends only on the compressibility coefficient $\alpha$. For convenience, let us first introduce a "scaled pressure profile" function $\Psi(z)$:

$$\Psi(z) = \frac{1 - e^{-2\alpha g(z-z_{BOT})}}{1 - e^{-2\alpha g(z_{TOP}-z_{BOT})}} \quad (A.13)$$

satisfying $\Psi(z_{BOT}) = 0$ and $\Psi(z_{TOP}) = 1$. Dimensional analysis of $\Psi(z)$: $\Psi(z)$ is dimensionless, and both arguments $-2\alpha g(z - z_{BOT})$ and $-2\alpha gL$ inside the exponential terms must be dimensionless. Indeed, given the compressibility $\alpha_0 = \rho_0/p_0$, we have: $[2\alpha g(z - z_{BOT})] = \rho_0 g z / p_0 = [Pa/Pa]$; therefore the arguments of exponentials are dimensionless, as they should.
The closed form pressure profile solution expressed in terms of top and bottom pressures is finally obtained by inserting the mass flux Eq. A.12 into the pressure profile Eq. A.11:

$$p(z)^2 = p^2(z_{\text{BOT}}) \left\{ e^{-2\alpha g(z-z_{\text{BOT}})} - \Psi(z) e^{-2\alpha gL} \right\} + p^2(z_{\text{TOP}}) \Psi(z)$$

...whence:

$$p(z) = \sqrt{p^2(z_{\text{BOT}}) \left\{ e^{-2\alpha g(z-z_{\text{BOT}})} - \Psi(z) e^{-2\alpha gL} \right\} + p^2(z_{\text{TOP}}) \Psi(z)}$$

...and finally:

$$p(z) = \sqrt{p^2(z_{\text{BOT}}) \left\{ e^{-2\alpha g(z-z_{\text{BOT}})} - \Psi(z) e^{-2\alpha gL} \right\} + p^2(z_{\text{TOP}}) \Psi(z)} \quad (\text{A.14})$$

with $\Psi(z)$ given by Eq. A.13. Note that this steady pressure profile satisfies exactly the two boundary conditions at $z_{\text{BOT}}$ and at $z_{\text{TOP}}$. It depends only on these two boundary pressures and on two coefficients (gas compressibility $\alpha$ and acceleration field $g$), but it does not depend on Darcy permeability $k$. On the other hand, the steady mass flux density $\phi_0$, in Eq. A.12, depends on the Darcy permeability $k$, as it should.

## B. APPENDIX B: Pressure solution of linearized transient gas flow

We develop in this APPENDIX infinite series solutions applicable to two different linearized diffusion PDE's for 1D transient gas flow, which appear in the main text in Sect. 6:

* Equation 6.4 for the "Direct Linearization" approach, dubbed "Linear" solution (for short); and

* Equation 6.8 for the "Quadratic linearization" approach, dubbed "Quadratic solution" (for short).

The unknown variable in those diffusion equations is either a "Klinkenberg pressure" named $P_B(z, t)$ with $P_B(z, t) = p(z, t) + b$ (for the "Linear" approximation), or a quadratic pressure named $\Psi(z, t)$ with $\Psi(z, t) = (p(z, t) + b)^2$ (for the "Quadratic" approximation). For convenience, we will name $\Theta(z, t)$ the dependent variable, with either $\Theta = P_B = p + b$ or else $\Theta = \Psi = (p + b)^2$. In both cases, $\Theta(z, t)$ satisfies a classical diffusion PDE with an Initial Condition (IC) $\Theta_0(z)$ which is assumed to be a constant initial profile here ($\Theta_0(z) = \Theta_0$ constant), and Dirichlet Boundary Conditions (BC's) at the two end points ($\Theta_1$ at $z = 0$, $\Theta_2$ at $z = L$).

The $\Theta$-Diffusion equation with its IC and its Non-Homogeneous BC's is:
$$ 
\begin{array}{l}
\displaystyle \frac{\partial \Theta}{\partial t} - D \frac{\partial^2 \Theta}{\partial x^2} = 0 \quad x \in [0, L] \quad t > 0 \\
\\
t = 0 \qquad \Theta(x, 0) = \Theta_0(x) \quad x \in [0, L] \\
\\
x = 0 \qquad \alpha_1 \Theta + \beta_1 \frac{\partial \Theta}{\partial x} = \gamma_1 \quad t > 0 \\
\\
x = L \qquad \alpha_2 \Theta + \beta_2 \frac{\partial \Theta}{\partial x} = \gamma_2 \quad t > 0.
\end{array}
$$

The diffusion coefficient $D$ has units $m^2/s$; it is assumed constant in space and time (this is implied by the linearization approximations that we have applied to the gas flow equation). The boundary parameters ($\alpha_1, \beta_1, \gamma_1$ and $\alpha_2, \beta_2, \gamma_2$) are assumed constant. Obviously, for a well-defined boundary-value problem, $\alpha_1$ and $\beta_1$ cannot be null simultaneously, and similarly, $\alpha_2$ and $\beta_2$ cannot be null simultaneously. In the present case, we will impose Dirichlet Boundary Conditions, with:

$$ \bullet \quad \alpha_1 = 1, \beta_1 = 0, \gamma_1 = \Theta_1 $$

$$ \bullet \quad \alpha_2 = 1, \beta_2 = 0, \gamma_2 = \Theta_2 $$

In general, the initial condition profile $\Theta_0(x)$ is assumed to be square-integrable on $[0, L]$, that is $\Theta_0 \in L^2([0, L])$. In the present case, we are interested in imposing a constant initial profile, that is:

$$ \bullet \quad t = 0: \Theta(x, 0) = \Theta_0(x) = \Theta_0 \text{ constant.} $$

The solution method to solve the above diffusion problem (in both cases: "Linear" & "Quadratic") has two steps. The first step is to decompose the $\Theta$-diffusion problem additively into a steady state component $\Theta_\infty(x)$ plus a residual transient component $\theta(x, t)$ which satisfies a "homogeneous" diffusion problem with zero BC's ("homogeneous" BC's). The second step is to solve the resulting "homogeneous" $\theta$-diffusion problem by the method of separation of variables ($\theta(x, t) = X(x) T(t)$).<fnref n="3" />

Let us then examine first the steady state solution. In the case of two Dirichlet conditions, with $\Theta_1$ imposed at left and $\Theta_2$ imposed at right, the boundary conditions are always compatible, and there exists a unique steady state solution ($t \to \infty$), which is the well known linear conduction profile:

$$ \Theta_{\infty}(x) = \frac{\Theta_2 - \Theta_1}{L}x + \Theta_1 \qquad (B.1) $$

We now study the transient diffusion regime for the $\theta$-diffusion equation with its new initial and boundary conditions, where the dependent variable is defined as: $\theta(x, t) = \Theta(x, t) - \Theta_\infty(x)$. We wish to solve the problem satisfied by $\theta(x, t)$, which represents the "perturbation" of $\Theta(x, t)$ around the steady state solution profile $\Theta_\infty(x)$. Subtracting the steady state solution $\Theta_\infty(x)$ from the transient variable $\Theta(x, t)$, the new "perturbed" variable

$$ \theta(x, t) = \Theta(x, t) - \Theta_{\infty}(x) \qquad (B.2) $$

<footnote n="3">³ Applying the method of separation of variables is mathematically equivalent to finding the spatial eigenvalues and eigenfunctions of the spatial diffusion operator (Sturm-Liouville problem).</footnote>
is found to be governed by the diffusion PDE with homogeneous boundary conditions (notice the lower case θ instead of Θ):

$$ 
\begin{array}{l}
\frac{\partial \theta}{\partial t} - D \frac{\partial^2 \theta}{\partial x^2} = 0 \quad x \in [0, L] \ t > 0 \\
t = 0 \quad \theta(x, 0) = \Theta_0 - \Theta_\infty(x) \quad x \in [0, L] \\
x = 0 \quad \theta(0, t) = 0 \quad t > 0 \\
x = L \quad \theta(L, t) = 0 \quad t > 0.
\end{array}
$$

The next steps of solution are as follows, leading to an infinite series solution for the perturbed variable θ(x, t) (lower case).

1. Separation of variables.

Let formally θ(x, t) = X(x)T(t), and inject this expression in the diffusion PDE, to obtain:

$$ XT' = DX''T \quad \text{whence} \quad \frac{X''}{X} = \frac{T'}{DT}. $$

Note that we have here the equality of two functions, one is a function of space x, the other is a function of time t, this for all values of x and t. These two functions are therefore necessarily constant (independent of both space and time):

$$ \frac{X''}{X} = \frac{T'}{DT} = \text{constant}. $$

The temporal part T(t) can be easily integrated into an exponential function of time. This function T(t) can be assumed not to diverge as time grows; in fact it will later be seen (proved) that it does not diverge with time.

Let us impose that the constant in the previous equation be negative, that is, we let constant = -λ². Then X(x) and T(t) must satisfy the two ODE's:

$$ 
\begin{array}{l}
X'' + \lambda^2 X = 0 \\
T' + D\lambda^2 T = 0.
\end{array}
$$

2. Elementary solutions for the spatial eigen-functions $X_n(x)$

$X_n(x)$ is solution of the following 2nd order ODE problem with 2 zero order boundary conditions:

$$ X'' + \lambda^2 X = 0 \quad \text{and} \quad X(0) = X(L) = 0. $$

The general solution is of the form

$$ X(x) = A e^{i \lambda x} + B e^{-i \lambda x}. $$

The two boundary conditions at $x = 0$ and $x = L$ then impose that A and B satisfy the system:

$$ \begin{bmatrix} 1 & 1 \\ e^{i \lambda L} & e^{-i \lambda L} \end{bmatrix} \begin{bmatrix} A \\ B \end{bmatrix} = 0. $$
But for (A, B) to be different from the trivial solution (0, 0), the determinant of the system must be zero; this yields the auxiliary condition :

$$ \begin{vmatrix} 1 & 1 \\ e^{i\lambda L} & e^{-i\lambda L} \end{vmatrix} = e^{-i\lambda L} - e^{i\lambda L} = -2i \sin \lambda L = 0. $$

There exists a countable infinity of solutions:

$$ \lambda_1 = \frac{\pi}{L} \quad \lambda_2 = \frac{2\pi}{L} \quad \dots \quad \lambda_n = \frac{n\pi}{L} \quad \dots $$

For each $\lambda_n$, there exists an infinite set of functions that differ only by a multiplicative factor, and such that $A + B = 0$. Imposing also that $||X_n||_2 = \int_0^L |X_n(x)|^2 dx = 1$, one can choose a (normalized) family of functions which generates all these solutions:

$$ X_n(x) = \sqrt{\frac{2}{L}} \sin(\lambda_n x) $$

Any solution $X(x)$ is then of the form $C X_n(x)$, where C is a constant coefficient.

### 3. Orthogonality of the eigenfunctions $X_n(x)$

It can then be verified that

$$ \int_{0}^{L} X_{n}(x)X_{m}(x) dx = \delta_{nm}, \quad n, m = 1, 2, \dots, $$

where $\delta_{nm}$ is the Kronecker symbol (1 if $n = m$, 0 otherwise). This equality means that $(X_n)_{n=1,2,...}$ is a set of orthonormal functions of the Hilbert space $L^2_{[0,L]}$.

### 4. Infinite orthogonal expansion for $\theta(x, t)$ in space

From Sturm-Liouville theory, we may assume that any function $f(x)$ of $L^2_{[0,L]}$ can be decomposed on the family $(X_n)_{n=1,2,...}$ which (therefore) forms an orthogonal basis of the space $L^2_{[0,L]}$. One can then find a unique family of coefficients $(c_n)_{n=1,2,...}$ such that:

$$ f(x) = \sum_{n=1}^{\infty} c_n X_n(x). $$

Each coefficient $c_n$ can be expressed by projecting $f$ on the elements of the basis, that is:

$$ \int_{0}^{L} f(x) X_{m}(x) dx = \int_{0}^{L} X_{m}(x) \sum_{n=1}^{\infty} c_{n} X_{n}(x) dx. $$

Assuming that the summations can be inverted, this yields:

$$ \int_{0}^{L} f(x) X_{m}(x) dx = \sum_{n=1}^{\infty} c_{n} \int_{0}^{L} X_{m}(x) X_{n}(x) dx = c_{m} \int_{0}^{L} X_{m}(x)^{2} dx. $$

whence:

$$ c_m = \frac{\int_0^L f(x) X_m(x) dx}{||X_m||_2^2}. $$

Note that choosing $||X_m||_2 = 1$ yields equivalently: $c_m = \int_0^L f(x) X_m(x) dx$.
5. Temporal functions $T_n(t)$

One can now solve the ODE governing $T(t)$ in order to determine the $T_n(t)$'s. The result is:

$$T_n(t) = A_n e^{-D\lambda_n^2 t}.$$

It is seen here, as announced before, that the temporal functions $T_n(t)$ do not diverge indeed as time grows, since $D > 0$ and $\lambda_n^2 > 0$.

6. Infinite expansion for the space-time solution $\theta(x, t)$

At this point, we can alternatively express the decomposition of $\theta(x, t)$ for any fixed time $t$, over the spatial functions $X_n(x)$, with time-dependent coefficients as follows:

$$\forall t > 0 \ \exists(a_n(t))_{n=1,\infty} \quad \theta(x, t) = \sum_{n=1}^{\infty} a_n(t)X_n(x).$$

Injecting this expression in the PDE yields:

$$\sum_{n=1}^{\infty} (a'_{n}(t)X_{n}(x) - D a_{n}(t)X''_{n}(x)) = 0$$

whence

$$\sum_{n=1}^{\infty} [a'_{n}(t) + \lambda_{n}^{2} D a_{n}(t)] X_{n}(x) = 0$$

Therefore, we obtain:

$$\forall n > 0 \quad a'_{n} = -\lambda_{n}^{2} D a_{n}$$

Time integration of this simple ODE yields finally:

$$a_n(t) = A_n e^{-\lambda_n^2 D t}.$$

This confirms that the general space-time solution of the diffusion problem is of the form:

$$\theta(x, t) = \sum_{n=1}^{\infty} A_n e^{-\lambda_n^2 D t} \sin(\lambda_n x), \quad \lambda_n = \frac{n\pi}{L}, n = 1, 2, \dots \qquad (B.3)$$

where (recall) the "lower case" variable $\theta(x, t)$ (lower case) is the perturbation of the original variable $\Theta(x, t)$ (upper case) around the steady state solution $\Theta_{\infty}(x)$. Finally, from Eq. B.3, the initial condition in terms of the $\theta$-variable can also be expressed as an infinite series:

$$t = 0: \quad \theta(x, 0) = \sum_{n=1}^{\infty} A_n \sin(\lambda_n x), \quad \lambda_n = \frac{n\pi}{L}, n = 1, 2, \dots \qquad (B.4)$$

where the coefficients $A_n$ remain unknown so far.

7. Using the initial condition to determine the $A_n$ coefficients of the infinite series $\theta(x, t)$
We now take into account the known initial condition $\theta(x, 0)$ of the $\theta$-diffusion Eq., and we use it in order to determine the $A_n$'s in Eq. B.4. The initial condition for the $\theta$-variable is:

$$t = 0 : \quad \theta(x, 0) = \Theta_0 - \Theta_{\infty}(x) = \Theta_0 - \left( \frac{\Theta_2 - \Theta_1}{L} x + \Theta_1 \right), \quad x \in [0, L]$$

The above expression is simplified for convenience by introducing parameters $(\alpha, \beta)$, whence:

$$\begin{gathered} t = 0, \quad x \in [0, L] \quad \theta(x, 0) = (\Theta_0 - \Theta_1) - (\Theta_2 - \Theta_1) \frac{x}{L} = \alpha - \beta \frac{x}{L} \tag{B.5} \\ \text{with:} \quad \alpha = (\Theta_0 - \Theta_1); \quad \beta = (\Theta_2 - \Theta_1) \end{gathered}$$

The initial condition $\theta_0(x) = \alpha - \beta x/L$ of Eq. B.5 is now injected into the infinite expansion solution $\theta(x, t)$ using the following projection method: we obtain the coefficients $A_n$ by forming the scalar product of the infinite series initial function $\theta(x, 0)$ with the sinusoid $\sin(m\pi x/L)$ divided by $L$:

$$\begin{aligned} m &= 1, 2, \dots : \quad \frac{1}{L} \int_0^L \theta_0(x) \sin(m\pi x/L) dx \\ &= \frac{1}{L} \int_0^L \sum_{n=1}^\infty A_n \sin(n\pi x/L) \sin(m\pi x/L) dx \\ &= -\frac{1}{2L} \sum_{n=1}^\infty A_n (\delta(n+m, 0) - \delta(n-m, 0)) \\ &= -\frac{1}{2L} (A_{(-m)} - A_{(+m)}) \\ &= \frac{1}{2L} A_m \end{aligned}$$

where we used the identity $\cos(a+b) - \cos(a-b) = -2\sin(a)\sin(b)$, and $\delta(i,j)$ designates the Kronecker symbol: $\delta(i,j) = 1$ if $i = j$, else $\delta(i,j) = 0$ if $i \neq j$. The final result above is obtained by letting $A_{(-m)} = 0$ for $m \ge 1$ since the $A_m$'s are not defined for $m \le 0$. Thus, we have:

$$m = 1, 2, \dots : \quad A_m = \frac{2}{L} \int_0^L \theta_0(x) \sin(m\pi x/L) dx \qquad (B.6)$$

The last step for computing the $A_m$'s is now to insert the initial condition $\theta_0(x)$ of Eq. B.5 into Eq. B.6. This yields, for $m = 1, 2, \dots$:
$$
\begin{align*}
A_m &= \frac{2}{L} \int_0^L \theta_0(x) \sin(m\pi x/L) \, dx \\
&= \frac{2}{L} \int_0^L \left(\alpha - \beta \frac{x}{L}\right) \sin(m\pi x/L) \, dx \\
&= 2\alpha \int_0^1 \sin(m\pi s) \, ds - 2\beta \int_0^1 s \sin(m\pi s) \, ds \\
&= 2\alpha \frac{1 - (-1)^m}{m\pi} - 2\beta \frac{(-1)^m}{m\pi}
\end{align*}
$$

where we used integration by parts for the second term, and we used the identity $cos(m\pi) = (-1)^m$. The result for the coefficients $A_m$ is finally:

$$m = 1, 2, \dots : \quad A_m = \frac{2}{m\pi} \{\alpha + (\beta - \alpha)(-1)^m\} \qquad (B.7)$$

where the parameters contain the initial and boundary conditions: $\alpha = (\Theta_0 - \Theta_1)$; $\beta = (\Theta_2 - \Theta_1)$.

8. **Final form of the solution:**

The final explicit form of the infinite series solution of Eq. B.3 is then, in terms of the $\theta$-variable:

$$\theta(x, t) = \sum_{n=1}^{\infty} \frac{2}{n\pi} \{\alpha + (\beta - \alpha)(-1)^n\} e^{-\lambda_n^2 D t} \sin(\lambda_n x), \quad \lambda_n = \frac{n\pi}{L}, \quad n = 1, 2, \dots \quad (B.8)$$

The explicit solution in terms of the original $\Theta$-variable is then obtained by adding the steady state solution $\Theta_{infty}(x)$ to the above infinite series solution $\theta(x, t)$ of Eq. B.8, that is:

$$\Theta(x, t) = \theta(x, t) + \Theta_{\infty}(x) = \theta(x, t) + \left( \frac{\Theta_2 - \Theta_1}{L} x + \Theta_1 \right) \qquad (B.9)$$

In practice, the infinite series $\theta(x, t)$ must be truncated to a finite number of terms: we used 1001 terms in most of our applications in the text (which range from early times to large steady state times), although a much smaller number of terms could be adequate for most cases depending on the time scale of interest.

## C APPENDIX C: Flux solution of linearized transient gas flow

This APPENDIX presents, briefly, two linearized approximations for mass flux density $\phi(x, t)$ ($kg/s/m^2$, which are related to the pressure approximations (i) and (ii), [or equivalently (1) and (2)], developed in the previous APPENDIX B (Sec. B), and presented in the text in Eqs. 6.11, 6.12, 6.13, and 6.14. These equations were completed by a description of the procedures to obtain, respectively, the "Linear" approximation $p^{(1)}(x, t)$ and the "Quadratic" approximation $p^{(2)}(x, t)$.

The resulting flux approximations, named $\phi^{(1)}(x, t)$ and $\phi^{(2)}(x, t)$, are deduced from the infinite series pressure approximations $p^{(1)}(x, t)$ and $p^{(2)}(x, t)$ by applying Darcy's law for compressible gas flux with Klinkenberg effect (Sect. 5.1, Eq. 5.3). This Darcy equation is reformatted and reproduced here after neglecting the gravity term in $(dp/dz + \alpha g p)$:
$$ \varphi_Z = \rho(p) q_Z = -\frac{k_\infty}{\mu} \alpha (p+b) \left(\frac{dp}{dz}\right) \quad (C.1) $$

### C.1 "Linear" approximation of mass flux density: $\varphi(x, t)^{(1)}$

Inserting the first infinite series approximation of pressure, $p^{(1)}(x, t)$, in Darcy's law yields, after some calculations:

$$ \varphi(x, t)^{(1)} = -\frac{k_\infty}{\mu} \alpha (\bar{p} + b) \left[ \frac{p_2 - p_1}{L} + \frac{2}{L} \sum_{n=1}^{\infty} \{\alpha_1 + (\beta_1 - \alpha_1) (-1)^n\} \cos(\lambda_n x) e^{-\lambda_n^2 D t} \right] \quad (C.2) $$

where:

$D = \frac{1}{\Phi} \frac{k_\infty}{\mu} (\bar{p} + b) : \text{Diffusion coefficient } [m^2/s]$

$\lambda_n = n\pi/L (n = 1, 2, \dots) : \text{Eigenvalues of diffusion PDE}$

$p_0 = p(x, 0) \forall x \in [0, L] : \text{Initial pressure (uniform in space)}$

$p_1 = p(0), p_2 = p(L) : \text{Boundary pressures at } x=0 \text{ and } x=L$

$\bar{p} = (p_1 + p_2)/2 : \text{"Mean" pressure}$

$\alpha_1 = p_0 - p_1$

$\beta_1 = p_2 - p_1$

### Remark on the diffusion coefficient D.

The linearized diffusion coefficient $D$ [$m^2/s$] was given in the text by Eq. 6.9, also reproduced just above. It was noted in the text that the two linearized approximations ("linear" and "quadratic") are governed by the same diffusion coefficient D.

### Steady state flux.

By letting $t \to \infty$ in Eq. C.2, the steady state flux $\varphi(x, \infty)^{(1)}$ is retrieved:

$$ \varphi(x, \infty)^{(1)} = -\frac{k_\infty}{\mu} \alpha (\bar{p} + b) \left[ \frac{p_2 - p_1}{L} \right], \quad \forall x \in [0, L] \quad (C.3) $$

Remarkably, this steady state flux obtained from the "Linear" approximation (1) is exact. This can be seen by comparing it to the exact expression obtained in the text, Eq. 5.25, or even more clearly, to Eq. 5.26 expressed in terms of $\bar{p}$. Interestingly, the exact steady flux can also be expressed as follows in terms of Klinkenberg gas permeability $k_G(p)$ and gas density $\rho(p)$ at $p = \bar{p}$:

$$ \varphi(x, \infty)^{(1)} = -\rho(\bar{p}) \frac{K_G(\bar{p})}{\mu} \left[ \frac{p_2 - p_1}{L} \right], \quad \forall x \in [0, L] $$

### Remark on notations.

The boundary pressures $p_1 = p(0)$, $p_2 = p(L)$ are also named $p_{BOT}$ and $p_{TOP}$ in the text, and the 1D coordinate $x \in [0, L]$ is also named $z \in [0, L]$ in the text. This makes no difference, since gravity is neglected here, so that the coordinate axis can be horizontal or vertical equivalently.
### C.2 "Quadratic" approximation of mass flux density: $\varphi(x, t)^{(2)}$

Inserting the second infinite series approximation of pressure, $p^{(2)}(x, t)$, in Darcy's law yields, again after some calculations:

$$ \varphi(x, t)^{(2)} = -\frac{k_{\infty}}{\mu} \frac{\alpha}{2} \left[ \frac{P_2^2 - P_1^2}{L} + \frac{2}{L} \sum_{n=1}^{\infty} \{\alpha_2 + (\beta_2 - \alpha_2) (-1)^n\} \cos(\lambda_n x) e^{-\lambda_n^2 D t} \right] \quad (C.4) $$

where:

$$ D = \frac{1}{\Phi} \frac{k_{\infty}}{\mu} (\bar{p} + b) : \text{Diffusion coefficient } [m^2/s] $$

$$ \lambda_n = n\pi/L (n = 1, 2, \dots) : \text{Eigenvalues of diffusion PDE} $$

$$ P_0 = p_0 + b, \quad \forall x \in [0, L] : \text{Initial pressure + "b"} $$

$$ P_1 = p(0) + b, \quad P_2 = p(L) + b : \text{Boundary pressures + "b"} $$

$$ \bar{P} = (P_1 + P_2)/2 = \bar{p} + b : \text{"Mean" pressure + "b"} $$

$$ \alpha_2 = P_0^2 - P_1^2 = (p_0 + b)^2 - (p_1 + b)^2 $$

$$ \beta_2 = P_2^2 - P_1^2 = (p_2 + b)^2 - (p_1 + b)^2 $$

This second, "quadratic" flux approximation, is formally similar to the previous "linear" flux approximation: the infinite expansion is formally the same, but the coefficients are somewhat different; in particular $(\alpha_2, \beta_2)$ are different from the previous $(\alpha_1, \beta_1)$.

**Steady state flux.**

By letting $t \to \infty$ in Eq. C.4, the steady state flux $\varphi(x, \infty)^{(2)}$ is retrieved:

$$ \varphi(x, \infty)^{(2)} = -\frac{k_{\infty}}{\mu} \frac{\alpha}{2} \left[ \frac{(p_2 + b)^2 - (p_1 + b)^2}{L} \right], \quad \forall x \in [0, L] \quad (C.5) $$

Remarkably, this steady state flux obtained from the present "Quadratic" approximation (2) is again identical with the exact relation given in the text by Eq. 5.25 or Eq. 5.26, same as in the previous paragraph. This can be seen by developing the quadratic terms in $\varphi(x, \infty)^{(2)}$ using the identity $a^2 - b^2 = (a - b)(a + b)$. This leads finally to the same exact relation:

$$ \varphi(x, \infty)^{(2)} = -\frac{k_{\infty}}{\mu} \alpha (\bar{p} + b) \left[ \frac{p_2 - p_1}{L} \right], \quad \forall x \in [0, L] \quad (C.6) $$

### C.3 Conclusion on "Linear" and "Quadratic" approximations of flux $\varphi(x, t)$

In conclusion, while the transient flux approximations $\varphi(x, t)^{(1)}$ and $\varphi(x, t)^{(2)}$ are distinct in space and time, it is found that the resulting steady state fluxes $\varphi(x, \infty)^{(1)}$ and $\varphi(x, \infty)^{(2)}$ are identical, and both are exact. By extension, this suggests that the transient flux solutions (1) and (2) are reasonably accurate and close to each other at sufficiently large times. In the main text, we choose to focus on the quadratic approximation $\varphi(x, t)^{(2)}$ given by Eq. C.4 to analyze and illustrate the transient behavior of the mass flux density profiles (i.e., we exploit the analytical approximation $\varphi(x, t) \approx \varphi(x, t)^{(2)}$ over the full range of time scales).
**Acknowledgements** The first author (RA) acknowledges logistic support of the IRSN (Institut de Radioprotection et de Sûreté Nucléaire, France) per contract between IRSN and the Institut National Polytechnique de Toulouse (INPT): contract IRSN 2022-0115 INP 06247. The second author (MHB) was supported by a doctoral grant from the IRSN while he was a PhD student at the INPT. Finally, the authors thank the anonymous reviewers for their helpful comments and suggestions.

**Author Contributions** Authors contributions can be summarized as follows: article writing, literature review, and analytical solutions were performed by R. Ababou; numerical solutions were performed by M. H. Bahlouli, Z. Saâdi, and I. Canamón Valera; all authors reviewed and revised the intermediate versions and the final version of the manuscript, and all authors read and approved the final manuscript.

**Funding** Open access funding provided by Institut National Polytechnique de Toulouse. The first author (R. Ababou) was granted logistic support by the IRSN (Institut de Radioprotection et de Sûreté Nucléaire, France) per contract between IRSN and the Institut National Polytechnique de Toulouse (INPT): contract IRSN 2022-0115 INP 06247, and the second author (M. H. Bahlouli) was supported by a doctoral grant from the IRSN while he was a PhD student at the INPT. Besides, the authors have no other competing interests to declare that are relevant to the content of this article. The authors are solely responsible for the statements and results presented in this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Aadnoy, B.S., Finjord, J.: Analytical solution of the Boltzmann transient line sink for an oil reservoir with pressure-dependent formation properties. J. Petrol. Sci. Eng. 15(24), 343–360 (1996). https://doi.org/10.1016/0920-4105(95)00078-X

Ababou, R.: Capillary Flows in Heterogeneous and Random Porous Media, vol. 1, pp. 1–371. John Wiley & Sons, New Jersey (2018)

Amri, A., Saâdi, Z., Ababou, R.: Parametric sensitivity to capillary entry pressure in two-phase water-gas flow models: deep geologic disposal of radioactive waste. Transp. Porous Media 145(1), 13–43 (2022)

Bahlouli, M.H., Saâdi, Z., Ababou, R.: Validation of iTOUGH2-EOS5 against analytical solution of compressible gas flow in a porous medium. In: Abstract in 10th TOUGH Symposium 2023, Lawrence Berkeley National Laboratory (LBNL), Sep. 2023, Berkeley (CA), USA (2023). https://ut3-toulouseinp.hal.science/hal-04483193v1

Bird, G.A.: The DSMC Method. [Direct Simulation Monte Carlo]. Vol. 1. CreateSpace Independent Publishing Platform, (2nd ed.) (2013)

Brommundt, J., Kaempfer, T.U., Enssle, C.P., Mayer, G., Wendling, J.: Full-scale 3D modelling of a nuclear waste repository in the Callovo–Oxfordian clay. Part 1: thermo-hydraulic two-phase transport of water and hydrogen. Geol. Soc. London Special Publ. 400, 443–467 (2014)

Collins, R.E., Crawford, P.B.: Calculations of unsteady-state gas flow through porous media, corrected for Klinkenberg effects. AIME Trans. 198, 339–340 (1953)

De La Vaissière, R., Gerard, P., Radu, J.-P., Charlier, R., Collin, F., Granet, S., Talandier, J., Piedevache, M., Helmlinger, B.: Gas injection test in the Callovo–Oxfordian claystone: Data analysis and numerical modelling. Geol. Soc. London Special Publ. 400(1), 427–441 (2014)

Didier, M.: Etude du transfert réactif de l'hydrogène au sein de l'argiletexte. Theses. Université de Grenoble, France (Oct. 2012). https://theses.hal.science/tel-00770176

Dullien, F.A.L.: Porous Media: Fluid Transport and Pore Structure, pp. 1–396. Academic Press Inc, Amsterdam (1979)

Haynes, W.M., Lide, D.R., Bruno, T.J. (eds). CRC Handbook of Chemistry and Physics (97th Edition). CRC Press, Taylor & Francis Group, pp. 1–2643 (2017)

<figure><img src="image_15.png" /></figure>
Heid, J.G., McMahon, J.J., Nielson, R.F., Yuster, S.T.: Study of the Permeability of Rocks to Homogeneous Fluids. API Drilling and Production Practice. pp. 230–244 (1950)

Ho, C.K., Webb, S.W. (eds.): Gas Transport in Porous Media, vol. 20, pp. 1–444. Theory and Applications of Transport in Porous Media. Springer (2006)

Al-Hussainy, R., Ramey, H.J.: Application of real gas flow theory to well testing and deliverability forecasting. J. Petrol. Technol. **18**(5), 637–642 (1966)

Al-Hussainy, R., Ramey, H.J., Crawford, P.B.: The flow of real gases through porous media. J. Petrol. Technol. **18**(5), 624–636 (1966)

Jones, F.O., Owens, W.W.: A Laboratory Study of Low-Permeability Gas Sands. J. Petroleum Technology 1631–1640 (Sept. 1980)

Jones, S.C.: A rapid accurate unsteady-state Klinkenberg parameter. SPE J. 383–397 (1972)

Klinkenberg, L.J.: The Permeability of Porous Media To Liquids and Gases. Drill. Prod. Practice 200–213 (1941)

Knuksen, M.: Die Gesetze der Molekularstromung und der inneren Reibungsstromung der Gase durch Rohren. Ann. Phys. **28**, 75–130 (1909)

Marschall, P., Silvio Giger, S., De La Vaissière, R., Shao, H., Leung, H., Nussbaum, C., Trick, T., Lanyon, B., Senger, R., Lisjak, A., Alcolea, A.: Hydro-mechanical evolution of the EDZ as transport path for radionuclides and gas: insights from the Mont Terri rock laboratory (Switzerland). Swiss J. Geosci. **110**, 173–194 (2017). https://doi.org/10.1007/s00015-016-0246-z

Marulanda, C., Culligan, P.J., Germaine, J.T.: Centrifuge modeling of air sparging - a study of air flow through saturated porous media. J. Hazard. Mater. **72**(2–3), 179–215 (2000)

Mohammadmoradi, P., Kantzas, A.: Pore-scale permeability calculation using CFD and DSMC techniques. J. Petr. Sci. Eng. **146**, 515 (2016)

Ngianga Alalibo, T., Orukari Mercy, A.: Seepage of isothermal and adiabatic fluids through granular media. Int. J. Eng. Sci. Math. **5**, 16–24 (2016)

Oldenburg, C.M.: Geologic Carbon Sequestration: CO2 Transport in Depleted Gas Reservoirs. In: Ho, C.K., Webb, S.W. (eds). Gas Transport in Porous Media, Vol. 20. Theory and Applications of Transport in Porous Media. Springer, Chap. 26, pp. 419–425 (2006)

Russell, D.G., Goodrich, J.H., Perry, G.E., Bruskotter, J.F.: Methods for predicting gas well performance. J. Petrol. Technol. **18**(1), 99–108 (1966)

Saâdi, Z.: Gas-entry pressure impact on the evaluation of hydrogen migration at different scales of a deep geological disposal of radioactive waste. Sci. Rep. **14**, 6221 (2024). https://doi.org/10.1038/s41598-024-56454-y. (Tech. rep. Nature)

Saâdi, Z., Dymitrowska, M., Deleruyelle, F., Marsal, F.: An evaluation model of the impact of hydrogen ‘piston effect’ on water displacement in a deep geological repository for radioactive waste. Environ. Earth Sci. **79**, 434 (2020)

Sakhae-Pour, A., Bryant, S.L.: Gas permeability of shale. In: Proceedings of the SPE Annual Technical Conference and Exhibition Society of Petroleum Engineers. pp. 401–409 (2011)

Scheidegger, A.E.: The Physics of Flow Through Porous Media, revised, pp. 1–313. University of Toronto Press, Toronto (1960)

Shugard, A.D., Robinson, D.B.: A Simple Model of Gas Flow in a Porous Powder Compact. SANDIA Report SAND2014-2858. Tech. rep. Sandia National Laboratories, pp. 1–33 (Apr. 2014)

Singh, A., Kolditz, O.: Chap.8: Gas Flow. In: Kolditz, O. et al. (eds.) Thermo-Hydro-Mechanical-Chemical Processes in Fractured Porous Media. Springer-Verlag Berlin Heidelberg (2012). https://doi.org/10.1007/978-3-642-27177-9

Tidwell, V.C.: Air Permeability Measurements in Porous Media. In: Ho, C.K., Webb, S.W. (eds). Gas Transport in Porous Media. Vol. 20. Theory and Applications of Transport in Porous Media. Springer, Chap. 15, pp. 273–278 (2006)

Webb, S.W.: Gas Transport Mechanisms. In: Ho, C.K., Webb, S.W. (eds.) Gas Transport in Porous Media. Theory and Applications of Transport in Porous Media, vol. 20, pp. 5–26. Springer, New York (2006). (Chap. 2)

Wendling, J., Justinavicius, D., Sentis, M., Amaziane, B., Bond, A., Calder, N.J.: Gas transport modelling at different spatial scales of a geological repository in clay host rock. Environmental Earth Sciences **78**(221), 1–16 (2019). https://doi.org/10.1007/s12665-019-8230-3

Wendling, J., Treille, E., Dymitrowska, M., Pellegrini, D., Ahusborde, E., Jurak, M., Amaziane, B., Caro, F., Genty, A., Poskas, P.: Final report on benchmark studies on repository-scale numerical simulations of gas migration. Part 1: cell scale benchmark. FORGE Report D1.6. Tech. rep. European Commission (2013a)

Wendling, J., Treille, E., Dymitrowska, M., Pellegrini, D., Ahusborde, E., Jurak, M., Amaziane, B., Caro, F., Genty, A., Poskas, P.: Final report on benchmark studies on repository-scale numerical simulations of
gas migration, Part 2: module scale benchmark. FORGE Report D1.6. Tech. rep. European Commission (2013b)

Wendling, J., Treille, E., Dymitrowska, M., Pellegrini, D., Ahusborde, E., Jurak, M., Amaziane, B., Caro, F., Genty, A., Poskas, P.: Final report on benchmark studies on repository-scale numerical simulations of gas migration, Part 3: repository scale benchmark. FORGE Report D1.6. Tech. rep. European Commission (2013c)

Wu, Y.-S., Pruess, K., Persoff, P.: Gas Flow in Porous Media with Klinkenberg Effects. Transp. Porous Media **32**, 117–137 (1998). https://doi.org/10.1023/A:1006535211684

Yang, G., Weigand, B.: Investigation of the Klinkenberg effect in a micro/nanoporous medium by direct simulation Monte Carlo method. Phys. Rev. Fluids **3**(044201), 1–17 (2018)

Yokoyamaama, T., Nakashima, S., Murakami, T., Mercury, L., Kirino, Y.: Solute distribution in porous rhyolite as evaluated by sequential centrifugation. Appl. Geochem. **26**(8), 1524–1534 (2011). https://doi.org/10.1016/j.apgeochem.2011.06.017

Zhang, L., He, Y., Guo, C., Yu, Y.: Dynamic Material Balance Method for Estimating Gas in Place of Abnormally High-Pressure Gas Reservoirs. GeoScienceWorld: Lithosphere **2021**(6669012), 1–19 (2021). https://doi.org/10.2113/2021/6669012

Zhao, J., Yao, J., Zhang, M., Zhang, L., Yang, Y., Sun, H., et al.: Study of gas flow characteristics in tight porous media with a microscale lattice Bolzmann model. **6**(32393) (2016). https://www.nature.com/articles/srep32393

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.