# Permeability, Porosity and Klinkenberg Coefficient Determination on Crushed Porous Media

Sandra Profice, Didier Lasseux, Yves Jannot, Naime Jebara, Géraldine Hamon

► To cite this version:

Sandra Profice, Didier Lasseux, Yves Jannot, Naime Jebara, Géraldine Hamon. Permeability, Porosity and Klinkenberg Coefficient Determination on Crushed Porous Media. Petrophysics – The SPWLA Journal of Formation Evaluation and Reservoir Description, 2012, 53 (6), pp.430-438. hal-03830368

HAL Id: hal-03830368

https://hal.science/hal-03830368

Submitted on 2 Dec 2022

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
# Permeability, Porosity and Klinkenberg Coefficient Determination on Crushed Porous Media<fnref n="1" />

Sandra Profice<fnref n="2" />, Didier Lasseux<fnref n="2" />, Yves Jannot<fnref n="3" />, Naime Jebara<fnref n="4" /> and Gérald Hamon<fnref n="4" />

## ABSTRACT

Permeability estimation of poor permeable formations like tight or gas-shale reservoirs using a pulse-decay experiment performed on crushed samples has been shown in earlier works to be an interesting alternative for it is faster and less expensive than traditional transient tests performed on carefully prepared core plugs, although it is restricted to measurement in the absence of overburden pressure. Due to reservoir depletion during production, sample characterization over a wide range of pore-fluid pressure is essential. If the Darcy-Klinkenberg model is thought to be a satisfactory gas-flow model for these tight formations, the full characterization can be achieved by determining both the intrinsic permeability, $k_i$, and Klinkenberg coefficient, $b$.

In this work, the conditions under which reliable estimates of $k_i$, $b$ and porosity, $\phi$ can be expected from this type of measurement are carefully analyzed. Considering a bed of monodisperse-packed spheres and a complete physical model to carry out direct simulations and inversion of the pressure decay, important conclusions are drawn opening wide perspectives for significant operational improvement

of the method. In particular, it is shown that:

i) The particle size of the crushed sample must be well selected for a reliable pressure-decay signal record.

ii) The simultaneous determination of both $k_i$ and $b$ by inversion of the pressure-decay signal is very difficult because the sensitivities of the pressure decay to both coefficients are correlated.

iii) The porosity of the particles can be accurately estimated when the experimental setup has been properly calibrated (volumes of the chambers and of the porous sample). The precision on the estimation of this parameter is however strongly dependent on a bias on the crushed sample volume.

iv) When identification of $k_i$ and $b$ is possible, a very significant error may occur in the determination of the intrinsic permeability due to a bias on the porous sample volume. Errors on the estimated values of $\phi$ and $k_i$ due to a bias on the chamber volume are not very significant. Moreover, $b$ remains insensitive to bias on both the chamber and porous sample volumes.

## INTRODUCTION

Gas shales have become a topic of major interest over the past 10 years due to the large natural gas resource they represent. However, reliable characterization of gas shales is a challenging task since these unconventional reservoirs have very tight pores and therefore specific flow properties, in particular very significant Klinkenberg effects, at least after a sufficiently long period of reservoir depletion. These special features call upon a careful identification of traceable and reliable methods to determine single-phase flow properties that are of major importance to evaluate accumulation and production capabilities of the reservoir. Depletion

of the reservoir may take place over a wide range of pore pressure. Assuming the Darcy-Klinkenberg model is valid for accounting for single-phase gas flow in such a porous medium, a characterization of porosity, $\phi$, intrinsic permeability, $k_i$, and Klinkenberg coefficient, $b$, should provide relevant insights in the reservoir properties.

Since the permeability of shales ranges from tens of microdarcies to nanodarcies, the time needed to reach stationary flow while performing a steady-state experiment on such samples is relatively long (it typically varies as the inverse of $k_i$) even on short plugs. Whereas associated flow rates are small (and difficult to measure), it must be noted that the separate determination of $k_i$ and $b$ requires

<footnote n="1">¹Originally presented at the International Symposium of the Society of Core Analysts, Austin, TX, USA, September 18–21, 2011, Paper SCA2011-32.</footnote>

<footnote n="2">²I2M–TREFLE, Université de Bordeaux, CNRS Esplanade des Arts et Métiers, 33405, Talence Cedex, France;</footnote>

<footnote n="3">³LEMTA, Nancy-Université, CNRS 2, Avenue de la Forêt de Haye, BP 160, 54504, Vandoeuvre Cedex, France;</footnote>

<footnote n="4">⁴TOTAL – CSTJF, Avenue Larribau, 64018, Pau Cedex, France.</footnote>
several different experiments when performed under steady-state conditions. Although this technique is widely used and despite an alternative quasi-steady method proposed elsewhere (Jannot and Lasseux, 2012), methods based on fully unsteady techniques are appealing. In the late 1960s, Brace et al. (1968) elaborated the pulse-decay technique, which consists of applying a pressure pulse at one end of a core plug while recording the response at the opposite end. Since then, numerous studies have sought to improve the investigative setup and associated data interpretation methods. More recently, Luffel et al. (1993) proposed the use of pulse-decay methods on crushed samples. It involves recording the unsteady step of a pycnometry-type experiment using a device schematically represented in Fig. 1.

This option is particularly appealing as it shortens the test duration by increasing the medium exchange area with the penetrating gas while decreasing the penetration depth. In addition, the cost of the test is reduced because a small volume of cuttings recovered during the drilling is sufficient for the analysis. However, experiments under realistic overburden pressures are impossible.

Several variations of this technique were envisaged. For instance, Egermann et al. (2003) proposed a technique based on the flow of a viscous oil within saturated cuttings leading to a compression of the residual-gas content. Luffel et al. (1993) and Lenormand et al. (2010), among others, introduced similar methods relying on the gas desorption of entire plugs subjected to a suction or vacuum. Recently, Wang and Knabe (2010) worked on a pore-pressure oscillation setup based on the generation of a sinusoidal excitation at the upstream side of a sample that is preconditioned at a given pore pressure.

All these transient tests are highly sensitive to the porosity of the medium that governs the transient gas flow. If this parameter is considered as a given datum in the interpretative model, it must be known with extreme accuracy, otherwise, estimates of the remaining parameters will be severely biased. As a consequence, a method allowing the determination of all the three parameters ($k_l$, $b$ and $\phi$) from a single experiment is highly desirable.

The purpose of the present work is to analyze the capabilities of a pulse-decay experiment performed on a crushed sample to achieve this goal. To do so, a physical description on a model configuration is implemented and validated by an approximated analytical solution. Both the direct and inverse procedures are derived, which allows a) investigation of the available information in the pressure-decay record that is strongly related to the radius of the particles of the crushed sample, b) analysis of the possibility of extracting the three parameters $k_l$, $b$ and $\phi$ from the pressure-decay signal using a sensitivity analysis and c) study of the impact of the chamber and sample volume bias on the estimated parameters when inversion is possible.

<figure><img src="image_2.png" /><figcaption>Fig. 1—Schematic experimental setup for a pulse-decay investigation of a crushed sample.</figcaption></figure>

## PROCEDURE

### Modeling

The pulse-decay experiment carried out on the crushed sample must be first described by the flow model.

### Direct Model

The physical model built to describe the pulse-decay experiment on a crushed sample relies on several hypotheses. First of all, the gas is supposed to be ideal, a commonly used approximation that is appropriate for gases like $N_2$ or He at operating pore pressures. Furthermore, the flow is assumed to be isothermal and slow enough to prevent development of inertial (Forchheimer) effects. At the pore level, a first-order slip flow is considered, yielding a Darcy-Klinkenberg flow model at the macroscale. Lastly, the crushed medium is assimilated to a monodisperse-sphere pack, the porous material being rigid, homogeneous isotropic and dry while pressure is assumed to be uniform in the spaces between spherical particles at any time. Under these circumstances, the flow within the pack is one-dimensional along the radial direction of each sphere of radius $R_0$. By combining the mass-conservation equation with the ideal gas law and the momentum conservation equation (Darcy-Klinkenberg), it is described by the following equation on the pressure $P$ at the position $r$ in each particle at time $t$

$$ \phi \mu \frac{\partial P}{\partial t} = \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 k_l \left( 1 + \frac{b}{P} \right) P \frac{\partial P}{\partial r} \right). \quad (1) $$

The associated initial and boundary conditions are given by:

$$ P_0 (t = 0) = P_{0i}, \quad (2) $$

$$ P (r \neq R_0, t = 0) = P_{1i}, \quad (3) $$

$$ \left( \frac{\partial P}{\partial r} \right)_{r=0} = 0, \quad (4) $$
$$-\frac{\mu(V_0 + DV)}{S_T} \left(\frac{\partial P}{\partial t}\right)_{r=R_0} = k_l \left(1 + \frac{b}{P_0}\right) P_0 \left(\frac{\partial P}{\partial r}\right)_{r=R_0} . \quad (5)$$

Equation 3 indicates that the system is initially at equilibrium with the surrounding atmosphere at $P_{li}$ prior to the emission of the pressure pulse at $t=0$ (Eq. 2). The condition (Eq. 4) simply results from the symmetry of the problem while Eq. 5 expresses the balance between the mass flow-rate coming in all the particles of total surface area $S_T$ and the mass flow-rate out of the volume of gas available to the particles, $V_0+DV$. The interest is focused on $P_0(t) = P(r=R_0,t)$ which is the pressure in the reservoir $V_0+DV$ as measured in the experiment.

This initial boundary value problem under this complete form has no analytical solution. It was hence solved by means of a numerical procedure, specifically an explicit finite difference scheme of first-order in time and second-order in space, inspired from previous works on the pulse-decay experiment on core plugs (Jannot et al., 2007). Using $\psi = (P+b)^2$ and the notation $\psi_i^n$ for the nodal value of $\psi$ at $r = (i-1)\Delta r$ and $t = (n-1)\Delta t$ where $\Delta r$ and $\Delta t$ are the grid spacing and time step respectively, the scheme is given by

$$\psi_1^{n+1} = \psi_1^n + 3 \frac{\sqrt{\psi_1^n}}{\alpha} \frac{\Delta t}{\Delta r^2} (\psi_2^n - \psi_1^n), \quad (6)$$

$$\psi_i^{n+1} = \psi_i^n \frac{\sqrt{\psi_i^n}}{\alpha} \frac{\Delta t}{\Delta r} \left[ \left(\frac{1}{r_i} + \frac{1}{\Delta r}\right) \psi_{i+1}^n - \frac{2}{\Delta r} \psi_i^n + \left(\frac{1}{\Delta r} - \frac{1}{r_i}\right) \psi_{i-1}^n \right]; 2 \le i \le m-1, \quad (7)$$

$$\psi_m^{n+1} = \psi_m^n + \frac{1}{\tau} \frac{\Delta t}{\Delta r^2} \quad (8)$$

$$\left[-\frac{1}{4}\psi_{m-3}^n + \psi_{m-2}^n + \frac{1}{4}\psi_{m-1}^n - \psi_m^n\right],$$

$$\psi_i^0 = (P_i + b)^2 \quad ; \quad 1 \le i \le m-1, \quad (9)$$

$$\psi_m^0 = (P_{0i} + b)^2 . \quad (10)$$

This scheme, in which we have denoted

$\alpha = \frac{\mu\phi}{k_l}$ and $\tau = \frac{1}{\sqrt{\psi_m^n}} \left(\alpha + \left(\frac{2}{R_0} + \frac{3}{2\Delta r}\right) \frac{\mu V_0}{S_r k_l}\right)$, is subjected to a stability criterion on the time step, classical for a diffusion-like equation, such that $\Delta t < \min_i \left(\frac{\alpha \Delta r^2}{2\sqrt{\psi_i}}\right)$.

## Approximated Analytical Solution

The numerical scheme was validated on a simplified configuration for which an analytical solution can be determined. Such a solution is achieved after linearization of Eqs. 1 to 5. Noting that for perfect gas, $1/P$ represents the fluid compressibility, denoted $\beta_f$ in what follows, and assuming that gas experiences only slight pressure fluctuation during the test, $\beta_f$ remains quasi constant. Under these circumstances, the apparent permeability, $k = k_l (1+b\beta_f)$, is treated as a constant as well and the problem can be solved analytically.

The sequence of derivation of the solution splits in two parts. The set of equations is first Laplace transformed and solved in the Laplace domain. The solution is inverted making use of the Cauchy theorem, as detailed by Hsieh et al. (1980) in a second step. Following this procedure, the analytical solution on $P_0(t)$ is obtained as

$$P_0(t) = \left[ \frac{3\gamma K}{3\delta K R_0^2 + R_0^3} + \frac{2}{R_0} \sum_m \right] \quad (12)$$

$$\frac{\gamma K \exp\left(\frac{-\xi_m^2 K}{R_0^2} t\right) \sin(\xi_m)}{R_0 \left[\delta K \xi_m \cos(\xi_m) + (R_0 + 2\delta K) \sin(\xi_m)\right] + P_{li}^2} \right]^{1/2},$$

$$\text{with:} \qquad \gamma = \delta R_0^2 (P_{0i}^2 - P_{1i}^2), \quad (13)$$

$$\delta = \frac{\mu(V_0 + DV)\beta_f}{S_T k}, \quad (14)$$

$$\text{and} \qquad K = \frac{k}{\Phi \mu \beta_f}. \quad (15)$$

The $\xi_m$ are the roots of:

$$\tan(\xi_m) = \frac{\xi_m}{1 + \frac{\delta K \xi_m^2}{R_0}}. \quad (16)$$

In Fig. 2, we have represented the evolution of $P_0$ obtained from the numerical method along with the evolution predicted by the analytical solution given by Eq. 12. Parameters used for this comparison are $k_l = 10^{-18} m^2$, $b=0 Pa$, $\phi=10\%$, $R_0 = 10^{-2} m$, $V_0 = 20 \times 10^{-6} m^3$, $V_1 = 75.4 \times 10^{-6} m^3$, $DV = 26.5 \times 10^{-6} m^3$, $P_{0i} = 10 \times 10^5 Pa$, $P_{li} = 10^5 Pa$, $\mu = 1.8 \times 10^{-5} Pa.s$, $t_f = 60$ s, $NP = 1200$, $m = 500$ (see nomenclature for the significance of the different parameters).
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2-(a) Numerical and analytical pressure-decay curves, and (b) residues drawn from the comparison of $(P_0)_1$ and $(P_0)_2$ to $(P_0)_3$. $(P_0)_1$, numerical pressure decay keeping a variable $\beta_i$; $(P_0)_2$, numerical pressure decay for a constant $\beta_i$; $(P_0)_3$, analytical pressure decay for a constant $\beta_i$.</figcaption></figure>

The crushed sample volume $V_1$ was fixed to that of the standard core plug and the space trapped between the spheres, i.e. the dead volume (DV), was calculated assuming that the spherical particles are packed according to a perfect face-centered cubic lattice whose compactness reaches 74%. Note that for real crushed samples for which particle shapes may strongly deviate from sphericity, larger compactness, i.e. smaller DV, may be achieved.

In order to estimate the deviation introduced by the hypothesis of a constant compressibility, the numerical procedure was first run without any particular assumption using the direct model Eqs. (1) to (5) and the scheme described above, yielding $P_0(t)$ represented as $(P_0)_1$ in Fig. 2a. In a second step,

the numerical procedure was run using the same parameters with a constant $\beta_f$ equal to $1/P_m$, $P_m$ being the time average over $t_f$ of $(P_0)_1$. The corresponding result is represented as $(P_0)_2$ in Fig. 2a. The analytical solution in this last situation is represented as $(P_0)_3$ in Fig. 2a.

The numerical and analytical results for a constant $\beta_f$ are in perfect agreement since the residue (Fig. 2b) remains globally close to 1 mbar for the chosen number of space nodes. The matching can even be improved when $m$ is increased. This comparison clearly validates our numerical procedure. The comparison between the analytical solution and the numerical one simulated with a variable $\beta_f$ shows a much more significant discrepancy that increases with time as suggested in Wu et al. (1998). Hence, even if the pressure decay remains small during the test, the assumption of constant gas compressibility may introduce a significant error and should be used with great care for further parameter identification since uncontrolled bias may be introduced. The numerical procedure is now used to analyze the impact of the different experimental parameters on the pulse-decay test.

## RESULTS AND DISCUSSION

In this section, the numerical procedure is used to investigate the impact of the particle size of the crushed sample on the pressure decay and to perform a sensitivity analysis of $P_0(t)$ to the parameters $k_i$, $b$ and $\phi$. An inverse procedure based on the direct numerical model is further employed to analyze the bias on $V_0$ and $V_1$ on the estimated values of $k_i$, $b$ and $\phi$ when the estimation is possible.

### Duration of the Experiment

Several numerical tests performed with a particle size $R_0$ of the crushed sample ranging from 1 mm to 1 cm were performed keeping all other parameters the same, namely, $k_i = 10^{-21} m^2$, $b = 68.6 \times 10^5 Pa$, $\phi = 1\%$, $V_0 = 20 \times 10^{-6} m^3$, $V_1 = 75.4 \times 10^{-6} m^3$, $DV = 26.5 \times 10^{-6} m^3$, $P_{0i} = 50 \times 10^5 Pa$, $P_{1i} = 10^5 Pa$, $\mu = 1.8 \times 10^{-5} Pa.s$, $t_f = 200$ s, $NP = 2000$ and $m = 200$. Results of the simulations are reported in Fig. 3.

From these results, it must be emphasized that the period over which the pressure decays and contains useful information about the physical parameters $k_i$, $b$ and $\phi$ is strongly related to the grain size. In fact, this characteristic dimension determines the test duration as it is reduced by a factor of ~10 when the particle radius is decreased from 1 cm to 1 mm. Thus, the smaller the spheres, the faster the decay, an acceleration that can be easily understood as the result of an increase of the exchange surface $S_T$ in conjunction with a reduction of the characteristic dimension to be penetrated by the pressure signal while reducing $R_0$, $V_1$ remaining the same. However, just a few seconds of recording (as for $R_0 = 1$ mm, see Fig. 3) are insufficient to guarantee an accurate post-estimation of the medium properties, taking into account actual pressure fluctuations at short times induced by valve opening and induced
<figure><img src="image_5.png" /><figcaption>Fig. 3—Pressure decay curves for different R₀.</figcaption></figure>

temperature effects resulting from the rapid gas expansion from V₀ to V₀+DV which are not considered in this model. Consequently, for a given rock and predefined experimental conditions, the mean grain size must be carefully selected to provide a usable signal for the estimation of kᵢ, b and φ by signal inversion. This dimension must be considered even more crucial when considering real particle shapes. Indeed, the hypothesis of spherical particles implies that the surface-to-volume ratio is at minimum, hence the time to reach final pressure equilibrium is maximum. With particles of arbitrary shape having the same average dimension the available pressure-decay signal is expected to be even shorter.

## Sensitivity Analysis

We now analyze the sensitivity of P₀(t) to the parameters that are to be estimated, namely kᵢ, b and φ. Several runs were performed with φ = 10%, R₀ = 10⁻² m, V₀ = 20 x 10⁻⁶ m³, V₁ = 75.4 x 10⁻⁶ m³, DV = 26.5 x 10⁻⁶ m³, P₀ᵢ = 50 x 10⁵ Pa, P₁ᵢ = 10⁵ Pa, μ = 1.8 x 10⁵ Pa.s, tᵦ = 12600 s, NP = 12600 and m = 200. The remaining parameters kᵢ and b were those reported in Table 1. Numerical results obtained on P₀(t) for the different (kᵢ, b) pairs are represented in Fig. 4.

Figure 4, in which all pressure signals are superimposed, highlights a key point related to the pulse-decay experiment on crushed core plugs, namely, the difficulty to discriminate the independent role of kᵢ and b on a given pressure decay. Indeed, the direct simulations show that the same signal can be obtained from tests performed with different combinations of the two parameters kᵢ and b that can be varied over an order of magnitude (and even more). The physical explanation of this behavior is as follows. Because the pressure decay remains small, the Klinkenberg term b/P experiences small variations so that kᵢ(1+b/P) remains quasi-constant during the test, making impossible the identification of the separate role of kᵢ and b on P₀(t). In other words, this indicates that kᵢ and b appear

<table><caption>Table 1—(kᵢ, b) Pairs Used for the Simulations of Fig. 4</caption><thead><tr><th>k<sub>i</sub> (m<sup>2</sup>)</th><th>b (Pa)</th></tr></thead><tbody><tr><td>2.26 x 10<sup>-21</sup></td><td>22.4 x 10<sup>5</sup></td></tr><tr><td>2.12 x 10<sup>-21</sup></td><td>24.8 x 10<sup>5</sup></td></tr><tr><td>1.56 x 10<sup>-21</sup></td><td>38.9 x 10<sup>5</sup></td></tr><tr><td>1.00 x 10<sup>-21</sup></td><td>68.6 x 10<sup>5</sup></td></tr><tr><td>7.89 x 10<sup>-22</sup></td><td>90.7 x 10<sup>5</sup></td></tr><tr><td>6.21 x 10<sup>-22</sup></td><td>11.9 x 10<sup>6</sup></td></tr><tr><td>2.84 x 10<sup>-22</sup></td><td>27.8 x 10<sup>6</sup></td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 4—Pressure decay curves for the different (kᵢ, b) pairs reported in Table 1.</figcaption></figure>

as correlated parameters for the signal P₀(t), a feature that is more pronounced when kᵢ (and φ) are small, as for instance in the range of interest for tight or gas shale reservoirs. A clear illustration can be made from the sensitivity of P₀(t) to each of the three parameters, kᵢ, b and φ.

The reduced sensitivity (Sᵣ)ᵢ of the measured pressure P₀ to the parameter θᵢ is defined by:

$$ (S_{r})_{i} = \theta_{i} \frac{\partial P_{0}}{\partial \theta_{i}}. \qquad (17) $$

It basically expresses how the evolution of P₀ is modified upon a slight modification of the parameter θᵢ, all other parameters remaining the same. It should be noted that the identification of θᵢ from P₀(t) is possible only if (Sᵣ)ᵢ is larger than the sensitivity of the pressure sensor used to measure P₀(t). An optimal estimation of θᵢ also requires the experiment be designed to maximize (Sᵣ)ᵢ and minimize the θᵢ to θⱼ correlation.
Two different parameters $\theta_i$ and $\theta_j$ are correlated when their respective sensitivities ($S_{ri}$) and ($S_{ri}$) are proportional. In Fig. 5a, we have represented the sensitivity curves of $P_0(t)$ to $k_i$, $b$ and $\phi$ using $k_i = 10^{-21} m^2$, $b = 68.6 \times 10^5 Pa$, $\phi = 1\%$, $R_0 = 10^{-2} m$, $V_0 = 20 \times 10^{-6} m^3$, $V_1 = 75.4 \times 10^{-6} m^3$, $DV = 26.5 \times 10^{-6} m^3$, $P_{1i} = 10^5 Pa$, $\mu = 1.8 \times 10^{-5} Pa.s$, $t_f = 1400$ s, $NP = 14000$ and $m = 200$.

Clearly, the sensitivity of $P_0(t)$ to $\phi$ greatly exceeds those to $k_i$ and $b$, the latter being the lowest one. The large sensitivity to $\phi$ is indeed expected since the pulse-decay test in that case is a pycnometry-type experiment, a result that contrasts with the classical pulse decay test performed on an entire plug. All the three sensitivities remain large compared to the precision of a standard pressure sensor.

Figure 5b points out the correlation of $k_i$ and $b$ for moderate pressure pulses $P_{0i}$ (~50 bars and less) through the quite perfect straight lines of the graph representing the sensitivity to $b$ versus that to $k_i$. Using larger pressure pulses, near 100 bars and more, slightly weakens the correlation in accordance with the fact that the apparent local permeability $k_i$ ($1+b/P$) experiences larger variations over the pressure-decay period. However, with such an experimental design, the correlation remains so significant that the simultaneous identification of the three parameters $k_i$, $b$ and $\phi$ is seriously compromised as further analyzed in the next section.

## Inverse Procedure

The pressure-decay signal inversion is operated using a Levenberg-Marquardt numerical algorithm as in (Jannot et al., 2008). This method which is commonly used to treat least square curve fitting problems returns the vector of parameters $\Theta$ ensuring the best correspondence between a set of measured data ($t_i$, $p_i$) (here they are generated with the direct model) and the estimated values $P(t_i, \Theta)$ derived from the solution to the direct problem. Optimal estimated parameters are found by minimizing the quadratic residues:

$$S(\Theta) = \sum_{i=1}^{NP} [p_i - P(t_i, \Theta)]^2. \quad (18)$$

The inverse procedure was run for two configurations of permeability and porosity. In both cases, the pressure-decay signal was generated using the direct simulation detailed above while superimposing a Gaussian noise which amplitude was fixed according to the pressure-decay magnitude (see Tables 2 and 3) to guarantee the convergence of the algorithm. Inversion was started giving two different initial sets ($k_{lin}$, $b_{in}$, $\phi_{in}$). The ($b_{in}$) were chosen according to the correlation proposed by Jones (1972):

$$(b_{in})_j = 0.189 ((k_{lin})_j)^{-0.36}, \quad (19)$$

whereas the ($\phi_{in}$) were taken equal to the nominal value $\phi$.

Note that since the final equilibrium pressure $P_f$ satisfies:

$$P_f = \frac{P_{0i}V_0 + P_{1i}(DV + \Phi V_1)}{V_0 + DV + \Phi V_1}. \quad (20)$$

$\phi$ can be effectively pre-estimated with accuracy if the various pressures and volumes needed for the calculation are measured with care and the experiment is long enough to reach final equilibrium. Direct simulations and inversions were carried out with the following parameters: $R_0 = 10^{-2} m$, $V_0 = 40 \times 10^{-6} m^3$, $V_1 = 75.4 \times 10^{-6} m^3$, $DV = 40 \times 10^{-6} m^3$, $P_{1i} = 10^5 Pa$, $\mu = 1.8 \times 10^{-5} Pa.s$ and $m = 200$. The results of the inversion are gathered in tables 2 and 3.

As expected, the strong correlation between $k_i$ and $b$ explained above justifies the fact that the estimated values are dependent upon the starting data. Despite the relatively small amplitude of the superimposed noise, errors on $k_i$ and $b$ are quite large. This analysis clearly shows that: a) convergence can only be achieved if initial values of $k_i$ and $b$ are close enough to the actual values, and b) when convergence is achieved, uncertainty of the estimated values of the parameters can be very large.

Although larger pressure pulses could be used to induce larger pressure decays, slightly reducing the correlation

<table><caption>**Table 2**—Estimates of $k_i$, $b$ and $\phi$ in Case 1 ($P_{0i} = 5 \times 10^5 Pa$, $t_f = 100$ s, $NP = 1000$) (pressure-decay amplitude: 18.7 mbars; -1.5 mbars <noise amplitude < 1.5 mbars)</caption><thead><tr><th></th><th>Real value</th><th>Starting value</th><th>Estimated value</th><th>Error (%)</th></tr></thead><tbody><tr><td>$k_i(m^2)$</td><td>1.0 x 10<sup>-19</sup></td><td>($k_{lin}$)<sub>1</sub> = 10<sup>-18</sup><br/>($k_{lin}$)<sub>2</sub> = 10<sup>-21</sup></td><td>3.17 x 10<sup>-20</sup><br/>1.73 x 10<sup>-20</sup></td><td>68.3<br/>82.7</td></tr><tr><td>$b(Pa)$</td><td>13.1 x 10<sup>5</sup></td><td>($b_{in}$)<sub>1</sub> = 5.7 x 10<sup>5</sup><br/>($b_{in}$)<sub>2</sub> = 68.6 x 10<sup>5</sup></td><td>45.8 x 10<sup>5</sup><br/>86.0 x 10<sup>5</sup></td><td>249.6<br/>556.5</td></tr><tr><td>$\Phi$</td><td>0.01</td><td>($\Phi_{in}$)<sub>1</sub> = 0.05<br/>($\Phi_{in}$)<sub>2</sub> = 0.05</td><td>0.01<br/>0.01</td><td>0<br/>0</td></tr></tbody></table>

<table><caption>**Table 3**—Estimates of $k_i$, $b$ and $\phi$ in Case 2 ($P_{0i} = 5 \times 10^5 Pa$, $t_f = 100$ s, $NP = 1800$) (pressure-decay amplitude: 42 mbars; -3 mbars <noise amplitude < 3 mbars)</caption><thead><tr><th></th><th>Real value</th><th>Starting value</th><th>Estimated value</th><th>Error (%)</th></tr></thead><tbody><tr><td>$k_i(m^2)$</td><td>1.0 x 10<sup>-21</sup></td><td>($k_{lin}$)<sub>1</sub> = 5.0 x 10<sup>-20</sup><br/>($k_{lin}$)<sub>2</sub> = 1.0 x 10<sup>-22</sup></td><td>3.51 x 10<sup>-21</sup><br/>1.53 x 10<sup>-21</sup></td><td>251<br/>53</td></tr><tr><td>$b(Pa)$</td><td>68.6 x 10<sup>5</sup></td><td>($b_{in}$)<sub>1</sub> = 16.8 x 10<sup>5</sup><br/>($b_{in}$)<sub>2</sub> = 157.2 x 10<sup>5</sup></td><td>17.0 x 10<sup>5</sup><br/>44.2 x 10<sup>5</sup></td><td>75.2<br/>35.6</td></tr><tr><td>$\Phi$</td><td>0.01</td><td>($\Phi_{in}$)<sub>1</sub> = 0.05<br/>($\Phi_{in}$)<sub>2</sub> = 0.05</td><td>0.01<br/>0.01</td><td>0<br/>0</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Fig. 5**—(a) Sensitivity curves to $k_l$ and $\phi$ for $P_{0i} = 50 \times 10^5$ Pa. (b) $k_l$-$b$ correlation curves (in absolute value) for different pressure pulses $P_{0i}$.</figcaption></figure>

between $k_l$ and $b$, this may cause experimental artifacts due to pressure fluctuations and thermal effects associated to an abrupt pressure drop at the valve opening as well as mechanical effects on the porous particles. Since the small variation of the apparent permeability leading to correlation between $k_l$ and $b$ is due to a large volume $V_0+DV$ compared to the pore volume $\phi V_1$, an alternative is to considerably reduce $V_0$ and/or $DV$. In the remaining part of this work, $V_0$ is considered to be zero (more precisely, in the above numerical model, $DV$ is set to zero while $V_0$ is taken as the actual value of $DV$) so that $V_s$ is directly connected to the gas supply. The actual experimental situation would be that of Fig. 1, for which the valve between $V_0$ and $V_s$ is opened in order to supply gas and generate the pulse and then closed so that the available gas volume is $DV$ only during the entire test. This configuration is now considered and to reduce $V_0$ (i.e. $DV$) even more, a sample made of a single spherical particle having

the volume of the standard core plug is chosen. Under these circumstances, reliable estimates of $k_l$, $b$ (and $\phi$) are obtained since inversion yields identical results whatever the starting values of the parameters, provided that they lead to convergence. In this situation, the objective is now to investigate the error on the parameters $k_l$, $b$ and $\phi$ estimated by inversion due to bias in the volume of the chamber $V_0$ and in the sample volume $V_1$ that may always occur due to measurement errors in these two parameters.

To do so, pressure-decay signals were generated with the nominal values of $V_0$ and $V_1$ while respective bias of $\pm 2\%$ in $V_0$ and $\pm 6\%$ in $V_1$ were introduced to carry out the inversion. It should be noted that a bias in $V_1$ is unavoidably reflected in $V_0$ ($V_0=DV=V_s-V_1$). Tests were performed with $k_l = 10^{-17} m^2$, $b = 2.49 \times 10^5$ Pa, $\phi = 10\%$, $R_0 = 2.6 \times 10^{-2} m$, $V_0=9.8 \times 10^6$ m$^3$, $V_1 = 75.4 \times 10^6$ m$^3$, $DV=0$ m$^3$, $P_{0i} = 50 \times 10^5$ Pa, $P_{1i} = 10^5$ Pa, $t_f = 30$ s, $NP=300$, $\mu = 1.8 \times 10^{-5}$ Pa.s and $m = 200$.

Results of the inverse procedure are reported in Table 4. They clearly show that the estimation of the Klinkenberg coefficient is insensitive to bias in both $V_0$ and $V_1$. In addition, a bias in $V_0$ seems to be directly reflected in $k_l$ and $\phi$, which means that the expected error due to this bias is not very significant provided the volume chamber is carefully characterized. However, the results show that the bias in $V_1$ has a very strong impact on both $k_l$ and $\phi$, requiring that the volume of the crushed sample must be very carefully measured if one wishes to extract reliable information on both $k_l$ and $\phi$ from the pulse-decay signal.

## CONCLUSION

A direct model and the related inverse procedure have been proposed to characterize in depth the physical properties of a crushed porous medium using a pulse-decay test. From this analysis, it first appears that the available information in the pressure-decay signal is intimately related to the particle size of the crushed sample. If small particle radii are used, the pulse relaxation time, which is even reduced for increasing permeabilities, may be so small that no interpretation of the pressure-decay signal is possible. This suggests that a pre-estimation (for instance on an entire core plug) of the order of magnitude of both $k_l$ and $\phi$ would provide an advantage in defining the appropriate particle mean radius for a further pulse-decay experiment on the crushed sample.

<table><caption>**Table 4**—Errors in the Estimates of $k_l$, $b$ and $\phi$, Due to a Bias of 2% in $V_0$ and 6% in $V_1$.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Estimated Values</th><th colspan="3">Errors</th></tr><tr><th>k<sub>l</sub> (m<sup>2</sup>)</th><th>b (Pa)</th><th>Φ</th><th>Δk<sub>l</sub>/k<sub>l</sub> (%)</th><th>Δb/b (%)</th><th>ΔΦ/Φ (%)</th></tr></thead><tbody><tr><td>V<sub>0</sub>+ΔV<sub>0</sub></td><td>1.02 x 10<sup>-17</sup></td><td>2.49 x 10<sup>5</sup></td><td>0.102</td><td>2.0</td><td>0</td><td>2.0</td></tr><tr><td>V<sub>0</sub>-ΔV<sub>0</sub></td><td>9.80 x 10<sup>-18</sup></td><td>2.49 x 10<sup>5</sup></td><td>0.098</td><td>2.0</td><td>0</td><td>2.0</td></tr><tr><td>V<sub>1</sub>+ΔV<sub>1</sub></td><td>5.15 x 10<sup>-18</sup></td><td>2.49 x 10<sup>5</sup></td><td>0.050</td><td>48.5</td><td>0</td><td>50.0</td></tr><tr><td>V<sub>1</sub>-ΔV<sub>1</sub></td><td>1.48 x 10<sup>-17</sup></td><td>2.49 x 10<sup>5</sup></td><td>0.154</td><td>48.0</td><td>0</td><td>54.0</td></tr></tbody></table>
A sensitivity analysis performed on the pressure decay indicates that, for a standard experimental design, the identification of both $k_i$ and $b$ is seriously compromised because these two parameters are strongly correlated. The physical origin of this lies in the restricted variation of the apparent gas permeability during the test. Without any pre-estimation (for instance on an entire core plug) of the order of magnitude of both $k_i$ and $b$, their estimation on crushed rock might be off by an order of magnitude. Conversely, and as expected, if volumes of both the chambers $V_0$ (and $DV$) and porous sample ($V_1$) are well characterized, a precise determination of the porosity can be performed as in a classical pycnometry-type experiment. Finally, a quick analysis of the impact of the bias in both $V_0$ and $V_1$ shows that when the inversion can be performed the resulting estimated values of $b$ are unaffected by these bias. The effect of bias in $V_0$ on both $k_i$ and $\phi$ may not be very significant since it seems to be directly reflected as an error in the estimated values of these coefficients. In contrast, a bias in $V_1$ induces very significant errors in both $k_i$ and $\phi$, suggesting that the volume of the crushed porous sample must be precisely measured for reliable identification of these two parameters from the pulse-decay signal.

## ACKNOWLEDGEMENTS

TOTAL is gratefully acknowledged for its support.

## NOMENCLATURE

* b: Klinkenberg coefficient (Pa)

* DV: dead volume (m³)

* k: apparent gas permeability $k = k_i (1+b/P)$ (m²)

* $k_i$: intrinsic permeability (m²)

* m: number of space nodes (-)

* NP: number of measurement points (-)

* $P/p_i$: gas pressure/rth experimental datum (Pa)

* $P_{0i}$: initial pressure pulse (Pa)

* $P_0(t)$: pressure at $r = R_0$ and $t$ (Pa)

* $P_{1i}$: initial steady-state pressure (Pa)

* $P_m$: mean value of $P_0$ over $t_f$ (Pa)

* r: radial coordinate (m)

* $R_0$: particles radius (m)

* $S_T$: particles total exchange area (m²)

* t: time (s)

* $t_i$: date of measurement (s)

* $t_f$: duration of the experiment (s)

* $V_s$: sample cell volume (m³)

* $V_0$: reference cell volume (m³)

* $V_1$: crushed sample volume (m³)

* $\beta_f$: ideal gas compressibility (Pa⁻¹)

* $\Delta r$: node spacing (m)

* $\Delta t$: time step (s)

* $\Delta X/X$: relative bias on X (-)

* $\phi$: porosity (-)

* $\Theta$: vector of parameters

* $\theta_i$: iᵗʰ parameter

## REFERENCES

* Brace, W.F., Walsh, J.B., and Frangos, W.T., 1968, Permeability of Granite Under High Pressure, *Journal of Geophysical Research*, **73**(6), 2225-2236.

* Egermann, P., Lenormand, R., Longeron, D., and Zarcone, C., 2003, A Fast and Direct Method of Permeability Measurement on Drill Cuttings, *Petrophysics*, **44**(4), 243-252.

* Hsieh, P.A., Tracy, J.V., Neuzil, C.E., Bredehoeft, J.D., and Silliman, S.E., 1980, A Transient Laboratory Method for Determining the Hydraulic Properties of Tight Rocks-1, Theory, *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, **18**(3), 245-252.

* Jannot, Y., and Lasseux, D., 2012, A New Method to Measure Gas Permeability of Weakly Permeable Porous Media, *Review of Scientific Instruments*, **83**(1), Paper 015113.

* Jannot, Y., Lasseux, D., Delottier, L. and Hamon, G., 2008, A Simultaneous Determination of Permeability and Klinkenberg Coefficient from an Unsteady-State Pulse-Decay Experiment, Paper SCA2008-09, presented at the SCA International Symposium, Abu Dhabi, UAE, October 29 - November 2.

* Jannot, Y., Lasseux, D., Vizé, G., and Hamon, G., 2007, A Detailed Analysis of Permeability and Klinkenberg Coefficient Estimation from Unsteady-State Pulse-Decay or Draw-Down Experiments, Paper SCA2007-08, presented at the SCA International Symposium, Calgary, Canada, September 10-12.

* Jones S.C., 1972, A Rapid and Accurate Unsteady-State Klinkenberg Permeameter, Paper SPE-3535, *SPE Journal*, **12**(5), 383-397.

* Lenormand, R., Bauget, F., and Ringot, G., 2010, Permeability Measurement on Small Rock Samples, Paper SCA2010-32, presented at the SCA International Symposium Halifax, Nova Scotia, Canada, September 18-21.

* Luffel, D.L., Hopkins, C.W., and Shettler, P.D., 1993, Matrix Permeability Measurement of Gas Productive Shales, Paper SPE-26633, presented at the SPE Annual Technical Conference and Exhibition, Houston, Texas, October 3-6.

* Wang, Y., and Knabe, R. J., 2010, Permeability Characterization on Tight Gas Samples Using Pore Pressure Oscillation Method, Paper SCA2010-30, presented at the SCA International Symposium, Halifax, Nova Scotia, Canada, September 18-21.

* Wu, Y., S., Pruess, K., and Persoff, P., 1998, Gas Flow in Porous Media with Klinkenberg Effects, *Transport in Porous Media*, **32**(1), 117-137.
## ABOUT THE AUTHORS

**Sandra Profice** is a PhD student at the University of Bordeaux (Institute of Mechanics and Engineering). She obtained an engineering degree from the National School of Electricity and Mechanics (Nancy), in Fluid Mechanics in 2010. Her research, carried out in collaboration with TOTAL, deals with the characterization of the physical properties of low permeable porous media, such as tight gas and gas shale reservoirs.

**Didier Lasseux** is a senior scientist at the Centre National de la Recherche Scientifique (Institute of Mechanics and Engineering in Bordeaux). He received an engineering degree from the Ecole Nationale Supérieure des Arts et Métiers and a PhD from the University of Bordeaux. He has been working on transport in porous media, more specifically on single and multiphase flow, for about 20 years. His recent research deals with the modeling of non Darcy flow and characterization of porous media properties. He has published more than 90 papers in journals and conference proceedings including SPE and SCA and has been granted 3 patents. He has co-founded the International Society for Porous Media (INTERPORE).

**Yves Jannot** is a senior scientist at the Centre National de la Recherche Scientifique (Laboratory of Theoretical and Applied Energetic and Mechanics). He received an engineering degree from the Ecole des Mines de Nancy and a PhD from the University of Nancy. He has been working on characterization of thermal properties of porous media for about 15 years. His recent research deals with the adaptation of methods for the characterization of thermal properties to the determination of the parameters related to fluid flow in porous media. He has published more than 60 papers in journals and conference proceedings and has been granted one patent.

**Naime Jebara** is a reservoir engineer at “Total SA” main technical and scientific research center in Pau, France. He is in charge of experimental pertrophysical characterization of conventional and unconventional special core analysis. He graduated from the Institut Français du Pétrole “IFP-School” with a MSc in petroleum engineering in 2009 and earned a degree in mechanical engineering from Institut National des Sciences Appliquées de Lyon (France) in 2008.

**Gerald Hamon** is currently an expert for Petrophysics with TOTAL. He earned his PhD in Fluid Mechanics in 1980 from the Institut National Polytechnique de Grenoble (France). He joined the oil industry in 1980 and worked in a variety of senior technical positions, including core analysis, field development in the North Sea and special expertise in the evaluation and numerical simulation of naturally fractured reservoirs. He authored or co-authored more than 60 papers for symposiums/journals of the Society of Core Analysts, the Society of Petroleum Engineers and the Society of Professional Well Log Analysts. In 2007, he received the Darcy Award from the Society of Core Analysts for lifetime technical achievement.