High-velocity compressible gas flow modelling through porous media considering the quadratic pressure gradient: a generalised analytical solution

Abstract

Transient compressible gas flow through porous media plays a pivotal role in energy production and storage. Understanding gas flow dynamics is crucial for interpreting subsurface reservoir complexities. Underground formations hold significant potential for natural gas production, or gas storage including CO₂ and hydrogen. Traditional steady-state core analysis methods face limitations in tight formations, necessitating the development and utilisation of transient methods to study flow behaviours. This paper introduces two analytical solutions addressing transient highly compressible gas flow through porous media. By incorporating quadratic pressure gradient and variable compressibility terms into the governing equations, more realistic models are proposed, applicable across a wide spectrum of high-pressure gas flow scenarios. These terms transform the linear diffusivity equation into a nonlinear Partial Differential Equation (PDE), which is solved analytically. The research evaluates the developed models, which are required for permeability quantifications. Comparative analysis against more simplified available models reveals discrepancies between the two, particularly in the high-pressure range, demonstrating the unreliability of available models for accurate representation of transient compressible gas flows.

1 Introduction

In the past few decades, the rising global energy demand has driven the exploration and extraction of fossil fuels to more challenging and complex environments (Nabizadeh et al., 2024). The drilling of production wells is now more frequent in geologically intricate high-pressure and high-temperature formations. This could be a more challenging process where the operating fluid in the subsurface porous media contains a compressible gas (Ababou et al., 2025; Scibek et al., 2025). The presence of a compressible gas under high pressure (or high-pressure gradient) can be problematic under certain circumstances (Nabizadeh, 2024; Nabizadeh and Jamiolahmady, 2025).

Transient compressible gas flow through porous media has various applications in both energy production and storage sectors. Accurate evaluation of gas flow is a critical step in achieving a deep understanding of subsurface reservoirs. Tight and ultra-tight reservoirs have received significant attention during the past decades due to their considerable potential for natural gas production and gas storage (CO₂, hydrogen, and air) (Jadhawar and Doryani, 2025; Larki et al., 2025; Larki et al., 2026; Nasieef et al., 2021; Seteyeobot et al., 2022; Shojaee et al., 2026; Shojaee et al., 2025). Conventional steady-state core analysis cannot be used for tight cores for samples with <10 mD permeability because of the long time required to reach steady-state. Instead, various transient methods have been proposed to study flow in such rocks.

Transient gas-flow behaviour is also relevant to subsurface CO₂ storage, where CO₂ is injected into porous formations under high-pressure conditions (Zheng et al., 2025). In particular, the nonlinear pressure response and pressure-dependent compressibility considered in this study can influence the interpretation of transient injection and pressure behaviour in low-permeability storage formations. The generalised analytical framework presented here can therefore be extended to CO₂ storage systems beyond conventional hydrocarbon reservoirs.

In recent years, the Pressure Pulse Decay (PPD) approach has become a widely used transient method for the permeability measurement of low-permeability samples. This method is initially presented by (Brace et al., 1968). The experimental setup of this method is straightforward. This method comprises an upstream tank connected to a downstream tank through a core sample in the middle. A pressure pulse is applied at the upstream tank, and the gas moves from the high-pressure tank (i.e., upstream tank) to the low-pressure one (i.e., downstream tank). The schematic of the PPD setup is presented in Figure 1. At the initial condition, the rock sample is in the equilibrium condition. Opening the upstream valve leads to gas flow through the core. The test ends as soon as the pressure difference between the upstream and downstream containers becomes negligible. Analysis of the differential pressure between those two tanks determines the permeability of the sample.

FIGURE 1

Feng et al., 2017b) (b) unsteady-state pulse-decay arrangement for core permeability measurement, illustrating the upstream and downstream pressure measurements used to monitor the pressure response during the test (Nazari Moghaddam and Jamiolahmady, 2016)

After pressure data of a PPD test are recorded, they need to be analysed to obtain permeability. There are analytically based methods for such analysis, albeit with simplifying assumptions (Brace et al., 1968). presented a simplified analytical solution for pressure transient data analysis that yields an estimation of core permeability. Several simplifying assumptions are involved in their method. First, a constant linear pressure gradient across the core is assumed. Moreover, gas compressibility variation with pressure is ignored in their approach, which could produce a considerable error. Additionally, the developed analytical solution only applies to non-adsorptive gases. Also, the Brace approach is known as the zero-porosity model since the rock’s compressive capacity (i.e., its pore volume) is neglected.

Several authors tried to increase the accuracy and reliability of the Brace method by reducing the number of simplifying assumptions mentioned above (Dicker and Smits, 1988; Jones, 1997). modified the method by including the compressive storage effect (i.e., pore volume of the sample) is considered. This makes the data analysis section more robust and reliable. The gas sorption effect is another unique feature of tight samples, which was neglected in previous models (Cui et al., 2009). combined the Langmuir adsorption isotherms with the Dicker and Smits method to extend the application of the PPD analytical solution when adsorption is important. However, the abovementioned models do not account for the compressibility effect (Haskett et al., 1988). suggested that pseudo-pressure and pseudo-time parameters should be considered for the PPD data analysis. Several authors used pseudo-time and pseudo-pressure definitions to consider pressure-dependent parameters of highly compressible gas and linearise the nonlinear diffusivity equation for compressible flows (Abdelmalek et al., 2018; Anderson and Mattar, 2007; Haskett et al., 1988; He et al., 2020; Lee and Holditch, 1982; Liu and Emami-Meybodi, 2021; Zhao et al., 2020).

Although the definition of pseudo-pressure and pseudo-time consider the compressibility effect, all the mentioned references only consider the linear pressure gradient term within the diffusivity equation (referred to as the linear form of diffusivity equation). In other words, they ignore the quadratic pressure gradient term in the diffusivity equation (referred to as the non-linear form of the diffusivity equation) (Liang et al., 2001; Odeh and Babu, 1988) presented an analytical solution of the nonlinear diffusivity equation with the quadratic pressure term, assuming several boundary conditions. They observed that the error between the linearised and nonlinear methods was negligible for cases with a low-pressure gradient, but it deviated as the pressure gradient increased. Several authors attempted to consider the quadratic pressure term for steady-state and pseudo-steady-state flows (Braeuning et al., 1998; Finjord and Aadnoy, 1989; Jelmert and Vik, 1996). Those references concluded that the quadratic pressure term should be included in the diffusivity equation for applications with a large pressure gradient.

The importance of the quadratic pressure gradient in dual-porosity porous media is highlighted by (Nie and Ding, 2010). They numerically solved nonlinear spherical flow equations for a fractured reservoir and compared the results with the simplified model without the quadratic pressure gradient term. Some distinct differences were observed, and it was suggested that quadratic pressure gradient should not be neglected, especially for tight and heavy oil reservoirs. Several other authors also focused on the importance of the quadratic pressure gradient term in dual-porosity models using numerical, semi-analytical, and analytical methods (Abbasi et al., 2017; Li et al., 2015; Liu et al., 2016).

As mentioned above, appropriate data analysis methods for PPD experiments with limited simplifying assumptions are required to achieve the best possible outcome from recorded experimental data. It has been highlighted that several publications have been focused on the importance of quadratic pressure gradient effect that exists in PPD experiments.

Another aspect of compressible flow modelling through porous media, which has received significant attention in recent years, is related to Forchheimer and inertial effects. Generally, inertia force implies the force due to the fluid’s momentum and it is added to the momentum equation by the term , which is a generalised term that encompasses the broader effects of fluid motion. The Forchheimer equation, on the other hand, only considers the additional resistance to fluid motion in porous media beyond the prediction of Darcy’s law (Jamiolahmady, 2003). Inertia in porous media generally stems from the change in flow direction rather than chaotic motion of particles that happens in flow through pipes. In any case, it is believed that such non-linear flow behaviour occurs at high velocities typically occurring in high permeability rocks rather than low permeability rocks studied through the PDD method, therefore, its effect is often ignored in the interpretation of data obtained by this technique for very low permeability rocks. Accordingly, in the derivation of equations presented here, the effect of inertia expressed by the Forchheimer equation in the porous media has been ignored here too.

To the best of the authors’ knowledge, an exact analytical solution for highly compressible flow through porous media considering the quadratic pressure gradient term for PPD setup, with time-dependent boundary conditions, has not yet been developed.

This study presents two analytical solutions to the transient nonlinear diffusivity equation: constant compressibility and variable compressibility. The solutions proposed in this research apply to both conventional linear (without quadratic pressure gradients) and nonlinear scenarios (with quadratic pressure gradients). These derivations are analytical solutions in the time domain, with no numerical assumptions during their derivations. The present paper constitutes Part I of this research and focuses on the analytical derivation and theoretical assessment of these solutions, while Part II presents complementary experimental investigations and core-scale numerical simulations to further evaluate and consolidate the findings of the present analytical study (Nabizadeh et al., 2024). The one-dimensional diffusivity equation and corresponding boundary conditions, together with the proposed solution technique, are discussed in Section 2. Then, the results of the analytical solutions are compared with previously published data analysis methods in Section 3. Finally, the main conclusions of this study are presented in Section 4.

2 Problem statement and analytical derivations

The general formulation of one-dimensional diffusivity equation for compressible gas flow is derived by combining the continuity equation and Darcy’s law. It is known that combining these equations produce a nonlinear PDE. In the literature, the quadratic pressure gradient term is usually ignored for simplicity. However, as mentioned in the introduction section, this simplifying assumption is not valid for high-velocity compressible gas flow during PPD experiments. This section aims to derive the governing equations without neglecting the quadratic pressure gradient term, which produce a nonlinear PDE.

Figure 1 presented a conceptual schematic view of the laboratory equipment, where the core sample is flooded with the operating fluid (i.e., gas) from the upstream to the downstream container. Similar to other relevant studies on the public domain, the flow paths are often assumed to be straight lines, making Cartesian coordinates the best option to represent flow behaviours. Therefore, continuity and Darcy equations in Cartesian coordinates are provided as follows:

The general 1D continuity equation in is (Brace et al., 1968):

The general form of Darcy’s law for 1D flow is:

The combination of these two equations and the definition of compressibility produce the diffusivity equation as follows:

Where P is pressure, is velocity, is density, is porosity, is permeability, and are fluid and total compressibility, respectively. As mentioned in the introduction section, this study presents two different solutions for the constant compressibility and variable compressibility cases in Section 2.1, and Section 2.2, respectively.

2.1 Nonlinear diffusivity equation considering quadratic pressure gradient and constant compressibility

For boundary conditions, the PPD setup is similar to the standard experiment mentioned in the literature (Brace et al., 1968). At initial condition, a high-pressure upstream container connects to a core and a downstream container, equilibrated at a lower pressure. Due to the unsteady-state nature of PPD experiment, boundary conditions are time dependent equations. Equation 7 represents the mass balance equations at and .

The corresponding initial and boundary conditions are:

Where is sample (core) pore volume and / are upstream/downstream volumes. These time-dependent boundary conditions specify the rate at which gas flows due to pressure gradient () through the cross-sectional area. Equations 3–7 can be subjected to a variable change to replace P(x,t) with a new dimensionless parameter as follows:

Applying this transformation yields:

The new diffusivity equation using this transformation is:

The modified boundary and initial conditions are:

In this format, is a dimensionless parameter. The following additional dimensionless parameters are employed to convert the diffusivity equation (Equation 10) and its associated boundary conditions (Equations 11 and 12a) into the dimensionless format.

Then:

This nonlinear PDE, Equation 14, can be transformed into a linear one using the following expression:

A new parameter ' ' is defined as:to give

Using the above transformation (Equation 19), one can convert the nonlinear PDE (Equation 14) into the following linear PDE by inserting Equation 19 into Equation 14:

The boundary and initial conditions are re-written as:

Now, the linearised PDE and its boundary and initial conditions are solved. This equation can be solved by the Laplace transform. Applying the Laplace transform (Nabizadeh et al., 2022; Nabizadeh et al., 2020; Widder, 2015) to Equations 20–22a yields:

Introducing two new parameters:

Then the new boundary conditions in the Laplace domain are:

In summary, using the Laplace transform method, a PDE at the time domain is converted into an ODE at the Laplace domain. The solution to this ODE is:

There are two unknown parameters (b and d) in Equation 26, and there are two boundary conditions. The solution of Equation 26 using these two boundary conditions (Equation 25) is:

Assuming and :

Equation 28 is the final solution at the Laplace domain, and a Laplace inverse transformation is needed to shift back from the Laplace domain to the time domain.

In the above equation, the Laplace parameter () is assumed to be complex, and the constant parameter (c) is a real positive that is large enough to cover all the singularities of the function .

To calculate the above integral, Cauchy’s residue theorem (Mitrinovic and Keckic, 1984) can be used as follows:

Where are the poles and are the associated residues. Now poles of the integrand should be examined. There is a simple pole at s = 0. The remaining poles can be determined using Equation 27:

It seems a complicated equation, but with a deeper look, one can conclude that the only solution to this equation occurs if is imaginary:

Where is a real number. Then the following equation can be derived:

Then:

The above equation has an infinite number of roots; therefore, the integrand has unlimited poles in addition to the poles at the origin. First, the residue of the pole at s = 0 should be determined:

Now, the residues of the remaining poles should be determined. Let’s assume:

It is known that is non-zero; therefore, all other remaining poles are simple poles, and their residue can be calculated by:

Then:

Now, all the residues are calculated, and one can compute the Laplace inverse transform using Equation 30:

At this stage, the equation can be transformed back to its original parameters form:

Introducing three new parameters to compress the final solution:

The final solution is obtained as:

In this form, the final pressure decay () is presented as a function of dimensionless time (), θ, , , and .

2.2 Nonlinear diffusivity equation considering quadratic pressure gradient and variable compressibility

In the previous section, it was assumed that compressibility was constant. This section aims to eliminate that simplifying assumption and make the solution more realistic, which could be applied to a broad range of applications with varying compressibility. Assuming gas compressibility is much higher than rock compressibility (), and the gas is an ideal gas (), Equation 3 would be converted to the following equation:

With the initial and boundary conditions:

To improve the readability of the article, detailed derivation has been provided in Supplementary Appendix A. Following the procedure mentioned in the appendix yields:

Where:

The results of the two different solutions presented in this and previous sections are compared and discussed in the next section.

3 Result and discussion

This study presents two analytical solutions to the generalised diffusivity equation for transient compressible flow through one-dimensional porous media, considering quadratic pressure term and with constant and variable compressibility assumptions. Since these two analytical solutions consist of an infinite series, two sets of MATLAB codes have been developed for the post-processing and data analysis. Results are validated against conventional diffusivity equation solution methods to determine the deviations between different models. In particular and for the verification section, the results are also compared with those of the (Dicker and Smits, 1988), the most common method in the literature for analysing PPD pressure data.

It must be noted that complementary experimental and numerical investigations that further verify and consolidate the findings of this research are available in a separate publication, as Part 2 of this research (Nabizadeh et al., 2024).

The exciting feature of the presented solution is that it depends on the initial absolute pressure of the containers and the pressure difference value. All previously presented analytical solutions are independent of initial absolute pressure conditions. It should be noted that the solutions provided here apply to continuum flows, where slip flows and net stress effects are negligible.

3.1 Solution 1: constant compressibility

In this section, the results of a parametric study to observe the impacts of different input parameters on the flow behaviour for the first solution (presented in Section 2.1) is discussed. This solution requires three input parameters: 1) gas compressibility, 2) initial upstream and downstream pressure, and 3) pore/container volume ratios. The results of different cases with different operational pressure conditions all with constant gas compressibility have been plotted in Figures 2–4 to investigate the importance of operating pressure on several parameters, including the pressure decay (), logarithm of pressure decay (), and the slope of logarithmic pressure decay versus time (i.e., α), used for analysing the PPD results. It has been assumed that both containers are identical and ten times bigger than core pore volume; therefore, . The definition of Lambda is presented in Equation 24. In this section, the arithmetic average of upstream and downstream containers pressure has been used to calculate the compressibility of an ideal gas expressed by .

FIGURE 2

FIGURE 3

FIGURE 4

Figure 2 presents pressure decay () versus dimensionless time (Equation 13). First, for a case with infinitesimal pressure difference (between upstream and downstream containers), the solution of this research is compared with that of the conventional Dicker and Smiths method. For this case, the quadratic pressure gradient term is negligible; therefore, the solution proposed here and that of the Dicker and Smits method should be similar, as demonstrated in this figure. These results confirm the accuracy of the developed generalised analytical solution producing the same results as those of the Dicker and Smits method. There is only a minor difference in the early times, caused by simplifying an infinite series to only the first parameter (Feng, 2017; Yang et al., 2015; Zhao et al., 2019). As mentioned in Equation 45 (i.e., final equation for solution #1), the final solutions consist of an infinite series. In the literature, such series are simplified by only considering the first term, whilst here, all terms of the infinite series are considered, which is more accurate than the previously published methods. This could explain the abovementioned minor error observed at the early times.

After this verification step, three more cases were considered with larger pressure difference values but the same average pressure value (to ensure that compressibility factor based on average pressure is constant). In this figure, the legend 3–1 refers to inlet and outlet pressure values of 3 and 1 bars, respectively. The same goes for other legends and the rest of this research. It can be observed that increasing the pressure difference leads to deviations from the first two plots, which demonstrates the impact of the initial pressure difference and the quadratic pressure term on the gas flow.

In the literature, the slope of logarithmic pressure decay () versus time/dimensionless time is used to determine the core sample properties (Feng et al., 2017b; Jia et al., 2018; Wang et al., 2019; Zhang et al., 2020; Zhao et al., 2019). In Figure 3, logarithmic pressure decay versus dimensionless time of data of Figure 2 is plotted, and the slope of this graph (i.e., referred to as α) is plotted in Figure 4. As mentioned above, the higher the initial pressure difference, the higher the deviation from previous conventional methods.

As shown in Figure 4, Dicker and Smits method produces a constant α (slope of logarithmic pressure decay versus pseudo time) versus time which means the permeability is constant with time. The case generated for verification purposes of the proposed solution also generated the same trendline, which means that the analytical solution proposed here produces a constant permeability value for scenarios with infinitesimal pressure differences. However, a deviation from the flat line is observed for other cases, suggesting the permeability is not constant with time, and it converges to that flat line at later times when the pressure difference is diminished. There is not a huge difference between the 2.4–1.6 bar case and the flat line because the pressure difference between the two containers is small enough to neglect the quadratic pressure gradient term. On the contrary, for the 3–1 bar scenario, the pressure difference is comparatively high; therefore, the deviation from the flat line at the early time is considerable, and at later times (dimensionless time = 5), it converges to that of the flat line. Other investigators (e.g. (Wang et al., 2019)) reported that assuming a constant spatial pressure gradient may not apply to those cases with high-pressure differences and may lead to significant permeability error. From the experimental point of view (Yang et al., 2015), performed several PPD experiments under different operating pressure conditions and observed that measured permeability increases as the gas pressure increases, although the net confining pressure remains constant (Abdelmalek et al., 2018). performed several PPD experiments and noted that experimental results show sensitivity to pore pressure. It can be concluded that the proposed analytical solution is compatible with experimental observations and previous theoretical understandings of compressible gas flow behaviours.

As mentioned before, gas compressibility is an input parameter and should be defined accordingly. Figure 5 depicts α versus dimensionless time for 3–1 bar case with three different pressure values used to determine compressibility factor according to (). These cases are: 1) Arithmetic average of upstream and downstream pressure (i.e., (1 + 3)/2 = 2 bar) 2) Upstream pressure (i.e., 3 bar), 3) Downstream pressure (i.e., 1 bar). These three cases are compared with the Dicker and Smits method results, as shown in Figure 5. This figure shows that the α value varies significantly between cases with different compressibility values, especially at early times. This suggests that the magnitude of α strongly depends on the compressibility factor, especially at earlier times, when the pressure difference is considerably high. Figure 5 also highlights the importance of correctly estimating gas compressibility for PPD data analysis. In this solution and for simplification purposes, the compressibility factor is constant. However, the results of Figure 5 highlight the importance of accurate compressibility factor estimation.

FIGURE 5

In the literature, it is common to apply core final equilibrium pressure (i.e., average core pressure) for further data analysis calculation (Abdelmalek et al., 2018; Feng et al., 2017a). To further investigate the impact of such an assumption, additional cases have been considered with gas compressibility calculated based on average pressure, with results shown in Figure 6. It is known that the higher the pressure, the lower the compressibility factor and its impact on the results. Figure 6, in particular, investigates the high-pressure cases in which the magnitude of compressibility factor and its impact are less.

FIGURE 6

It can be concluded from Figure 6 that although the magnitude of the impact of compressibility factor reduces at higher pressure values, the impact of the quadratic pressure term is not negligible, i.e., the same early time deviation noted in Figure 5 is also observed here. These results suggest that ignoring this parameter leads to a considerable error, even at high operational pressure conditions, with several applications in various geological environments.

3.2 Solution 2: variable compressibility

The results of the constant compressibility analytical solution presented in Section 2.1 were discussed in the previous section. The conducted sensitivity study demonstrated the considerable impact of the compressibility factor on gas flow behaviour. As described before, Section 2.2 relaxed the constant compressibility assumption and presented a new analytical solution considering variable compressibility. Before showing the corresponding results of this solution, the difference between these two solutions should be highlighted. For Solution #1, dimensionless time () is used, while Solution #2 utilises the pseudo dimensionless time () to capture the compressibility effect. As explained before, the definition of pseudo-time is not physical, since it is a function of both time and space. Without an additional simplifying assumption, the results of those two solutions are not comparable in a 2D plot due to the physical nature of pseudo-time. However, if is substituted with a constant (e.g., average core pressure), the pseudo-time and dimensionless time are similar.

The exact average pore pressure within the core evolves throughout a pressure pulse decay (PPD) test and cannot be measured directly using conventional experimental setups, where only the upstream and downstream reservoir pressures are recorded. Consequently, the average core pressure is commonly approximated by assuming a linear pressure distribution between the two ends of the core, which leads to a constant average pressure equal to the arithmetic mean of the boundary pressures. This assumption is widely adopted in pseudo-pressure and pseudo-time formulations (Wang et al., 2021), resulting in diffusion-dominated pressure propagation and a pressure profile that remains close to linear over much of the transient response. Under these conditions, the higher-order spatial variations in pressure are small, and the associated error in the average pressure is negligible. However, the assumption may become less accurate for ultra-tight rocks, highly heterogeneous formations, or fractured cores, where significant nonlinear pressure gradients may develop within the sample.

It should be noted that this does not mean the results of the first solution and those of the second one with the higher accuracy are similar. They have different simplifying assumptions with different derivations, and the second solution is more accurate.

Figure 7 depicts pressure decay for both solutions versus pseudo time at high operating pressure conditions. As expected, the higher the pressure difference, the higher the deviation from the Dicker and Smits method, which as mentioned before, is an analytical solution that ignores the quadratic term of diffusivity equation. Moreover, one can detect a considerable difference between Solution #1 (i.e., constant compressibility) and Solution #2 (i.e., variable compressibility) for the same operating pressure. It proves the importance of considering variable compressibility and the error arising from the constant gas compressibility assumption.

FIGURE 7

Figure 8 presents the slope of logarithmic pressure decay versus pseudo time (i.e., α) for the cases noted in Figure 7. As mentioned before, this parameter is used for further data analysis to obtain permeability. Figure 8 reveals a considerable gap between variable compressibility and constant compressibility models, leading to a substantial error during permeability calculation. For example, for 80–10 bar scenario, the predicted α from Solution #2 is ∼30% higher than the calculated one from Solution #1. This error could incorrectly be attributed to other physical phenomena (e.g., pressure-dependent permeability or sorption effect) for real experimental data if an unsuitable data analysis method is implemented.

FIGURE 8

The results presented in Figures 7, 8 demonstrate that the compressibility effect at high operating pressure conditions plays a crucial role and significantly affects the flow behaviours. At this stage, its impact under low operating pressure conditions is examined. To this end, three low-pressure scenarios have been designed, and the results of both constant compressibility and variable compressibility methods have been compared.

Figure 9 depicts pressure decay versus pseudo time for those low operating pressure condition cases. Although the gap between different cases is smaller compared to those in Figure 7, it cannot be neglected. α versus pseudo time data is presented in Figure 10. It is noted that here the error arising from neglecting variable compressibility is around 15%, which is again lower than 30% noted in Figure 8 between the constant and variable compressibility cases. Thus, the reported percentage should be interpreted as the maximum deviation over the simulated transient response rather than as a single representative error.

FIGURE 9

FIGURE 10

It should be noted that there are several complex physical behaviours associated with the dynamics of gas flow, but these results highlight the impacts of variable compressibility that should not be neglected to avoid any potential misinterpretations.

4 Conclusion

This research presented two analytical solutions to transient highly compressible gas flow through porous media. Quadratic pressure gradient term and variable compressibility terms were considered in the governing equations to propose more realistic solutions that apply to a broad range of cases dealing with high-pressure gas flows. These two terms convert the linear diffusivity equation into a nonlinear Partial Differential Equation (PDE). This nonlinear PDE was solved analytically, and an analytical solution in the time domain was presented.

The results of the developed models were used to determine the slope of logarithmic pressure decay versus time, which is required for permeability quantifications. Compared to prior simplified models, an error in the range of 30% for very high-pressure scenarios (i.e., 80–10 bar) was detected. In addition, neglecting quadratic pressure gradient term and variable compressibility for lower-pressure cases (e.g., 3–1 bar) caused around 15% error, proving those nonlinear terms cannot be ignored for transient compressible gas flows under higher pressure differences.

The developed equations and presented results improve our understanding of transient highly compressible gas flow through porous media. However, additional work is required to consolidate these findings. Several designed PPD experiments have been performed in our laboratory to investigate the reliability of the presented solutions. To avoid having a lengthy paper, the results of these experiments are included in a separate publication (i.e., Part 2 of this research). In addition, Part 2 will present some core-scale numerical simulations, which are complementary to these analytical and experimental investigations.

References

• 1

  AbabouR.BahlouliM. H.SaâdiZ.ValeraI. C. (2025). Single Phase Compressible Gas Flow in Porous Media: review and Advances: R. Ababou et al. Transp. Porous Media152 (11), 89. 10.1007/s11242-025-02226-9

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AbbasiM.IzadmehrM.KarimiM.SharifiM.KazemiA. (2017). Analytical study of fluid flow modeling by diffusivity equation including the quadratic pressure gradient term. Comput. Geotechnics89, 1–8. 10.1016/j.compgeo.2017.04.001

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AbdelmalekB.KarpynZ. T.LiuS.YoonH.DewersT. (2018). Gas permeability measurements from pressure pulse decay laboratory data using pseudo-pressure and pseudo-time transformations. J. Petroleum Explor. Prod. Technol.8 (3), 839–847. 10.1007/s13202-017-0376-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AndersonD.MattarL. (2007). An improved pseudo-time for gas reservoirs with significant transient flow. J. Can. Petroleum Technol.46 (07). 10.2118/07-07-05

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BraceW. F.WalshJ.FrangosW. (1968). Permeability of granite under high pressure. J. Geophys. Research73 (6), 2225–2236. 10.1029/jb073i006p02225

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  BraeuningS.JelmertT. A.VikS. A. (1998). The effect of the quadratic gradient term on variable-rate well-tests. J. Petroleum Sci. Eng.21 (3-4), 203–222. 10.1016/s0920-4105(98)00073-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  CuiX.BustinA. M. M.BustinR. M. (2009). Measurements of gas permeability and diffusivity of tight reservoir rocks: different approaches and their applications. Geofluids9 (3), 208–223. 10.1111/j.1468-8123.2009.00244.x

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  DickerA. I.SmitsR. M. (1988). A Practical Approach for Determining Permeability from Laboratory Pressure-Pulse Decay Measurements.

  • Google Scholar
  • View reference in article

• 9

  FengR. (2017). An optimized transient technique and flow modeling for laboratory permeability measurements of unconventional gas reservoirs with tight structure. J. Nat. Gas Sci. Eng.46, 603–614. 10.1016/j.jngse.2017.08.032

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  FengR.HarpalaniS.PandeyR. (2017a). Evaluation of various pulse-decay laboratory permeability measurement techniques for highly stressed coals. Rock Mech. Rock Eng.50 (2), 297–308. 10.1007/s00603-016-1109-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  FengR.LiuJ.HarpalaniS. (2017b). Optimized pressure pulse-decay method for laboratory estimation of gas permeability of sorptive reservoirs: part 1 – background and numerical analysis. Fuel191, 555–564. 10.1016/j.fuel.2016.11.079

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  FinjordJ.AadnoyB. S. (1989). Effects of the quadratic gradient term in steady-state and semisteady-state solutions for reservoir pressure. SPE Form. Eval.4 (03), 413–417. 10.2118/15969-pa

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  HaskettS. E.NaraharaG. M.HolditchS. A. (1988). A method for simultaneous determination of permeability and porosity in low-permeability cores. SPE Form. Eval.3 (03), 651–658. 10.2118/15379-pa

  • CrossRef
  • Google Scholar

• 14

  HeB.XieL.ZhaoP.RenL.ZhangY. (2020). Highly efficient and simplified method for measuring the permeability of ultra-low permeability rocks based on the pulse-decay technique. Rock Mech. Rock Eng.53 (1), 291–303. 10.1007/s00603-019-01911-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  JadhawarP.DoryaniH. (2025). Geochemical and Geomechanical Synergies During Underground Hydrogen Storage Schemes, Subsurface Hydrogen Energy Storage. Elsevier, 239–263.

  • Google Scholar
  • View reference in article

• 16

  JamiolahmadyM.DaneshA.TehraniD.DuncanD. (2003). Positive effect of flow velocity on gas–condensate relative permeability: network modelling and comparison with experimental results. Transp. Porous Media52 (2), 159–183. 10.1023/a:1023529300395

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  JelmertT. A.VikS. A. (1996). Analytic solution to the non-linear diffusion equation for fluids of constant compressibility. J. Petroleum Sci. Eng.14 (3-4), 231–233. 10.1016/0920-4105(95)00050-x

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  JiaB.TsauJ.-S.BaratiR. (2018). Experimental and numerical investigations of permeability in heterogeneous fractured tight porous media. J. Nat. Gas Sci. Eng.58, 216–233. 10.1016/j.jngse.2018.08.011

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  JonesS. C. (1997). A technique for faster pulse-decay permeability measurements in tight rocks. SPE Form. Eval.12 (01), 19–26. 10.2118/28450-pa

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  LarkiS.SkaugeA.SorbieK.MackayE. (2025). Numerical simulation of viscous fingering in CO₂ storage: addressing the limitations of laboratory-derived relative permeabilities. Int. J. Greenh. Gas Control148, 104523. 10.1016/j.ijggc.2025.104523

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  LarkiS.SorbieK.SkaugeA.MackayE. (2026). The generation of flow-representative relative permeabilities from experimental data to model viscous fingering in CO2 storage. Geoenergy Sci. Eng.265, 214568. 10.1016/j.geoen.2026.214568

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  LeeW. J.HolditchS. A. (1982). Application of pseudotime to buildup Test analysis of low-permeability gas Wells with long-duration wellbore storage distortion. J. Petroleum Technol.34 (12), 2877–2887. 10.2118/9888-pa

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  LiJ.-T.NieR.-S.JiaY.-L.WangD.-L. (2015). Semianalytical solution of the nonlinear dual-porosity flow model with the quadratic pressure gradient term. Math. Problems Eng.2015, 1–8. 10.1155/2015/387368

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  LiangY.PriceJ. D.WarkD. A.WatsonE. B. (2001). Nonlinear pressure diffusion in a porous medium: approximate solutions with applications to permeability measurements using transient pulse decay method. J. Geophys. Res. Solid Earth106 (B1), 529–535. 10.1029/2000jb900344

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  LiuZ.Emami-MeybodiH. (2021). A unified approach to the nonlinearity of the diffusivity equation and assessment of pseudotime. SPE J.26 (01), 241–261. 10.2118/204219-pa

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  LiuW.YaoJ.ChenZ.LiuY. (2016). Effect of quadratic pressure gradient term on a one-dimensional moving boundary problem based on modified Darcy’s law. Acta Mech. Sin.32 (1), 38–53. 10.1007/s10409-015-0526-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  MitrinovicD. S.KeckicJ. D. (1984). The Cauchy Method of Residues: Theory and Applications, 9. Springer Science and Business Media.

  • Google Scholar
  • View reference in article

• 28

  NabizadehA. (2024). High-Velocity Compressible Gas Flow Modelling in Porous Media. Heriot-Watt University.

  • Google Scholar
  • View reference in article

• 29

  NabizadehA.JamiolahmadyM. (2025). A methodology for predicting choked compressible gas flow in pores using dimensional analysis. Phys. Fluids37 (9), 096124. 10.1063/5.0292017

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  NabizadehA.SharifiM.AminshahidiB. (2020). An analytical solution for spherical fluid flow modeling in porous media with non-uniform initial pressure considering the supercharging effect. Petroleum.

  • Google Scholar
  • View reference in article

• 31

  NabizadehA.AbbasiM.SiavashiJ.SharifiM.MovagharM. R. K. (2022). Fluid flow modeling through pressure-dependent porous media: an analytical solution and a computational fluid dynamics approach. Groundw. Sustain. Dev.18, 100776. 10.1016/j.gsd.2022.100776

  • CrossRef
  • Google Scholar

• 32

  NabizadehA.DoryaniH.JamiolahmadyM. (2024). High-Velocity compressible gas flow modelling through porous media considering the quadratic pressure gradient: part 2. Experimental investigations and numerical analysis. Geoenergy Sci. Eng.

  • Google Scholar
  • View reference in article

• 33

  NasieefA. (2021). “Condensate recovery from tight gas-condensate reservoirs using a new method based on a gas-based chemical,” in SPE Annual Technical Conference and Exhibition? SPE, D012S065R001.

  • Google Scholar
  • View reference in article

• 34

  Nazari MoghaddamR.JamiolahmadyM. (2016). Slip flow in porous media. Fuel173, 298–310. 10.1016/j.fuel.2016.01.057

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  NieR.-S.DingY. (2010). Research on the nonlinear spherical percolation model with quadratic pressure gradient and its percolation characteristics. Nat. Sci.02 (02), 98–105. 10.4236/ns.2010.22016

  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  OdehA. S.BabuD. K. (1988). Comparison of solutions of the nonlinear and linearized diffusion equations. SPE Reserv. Eng.3 (04), 1202–1206. 10.2118/17270-pa

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  ScibekJ.KuboT.KoikeK.Achtziger-ZupančičP. (2025). Permeability of granitic drill core tested by steady flow and transient flow gas permeameter probes: improvements to methods and applications. Eng. Geol.351, 108023. 10.1016/j.enggeo.2025.108023

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  SeteyeobotI.JamiolahmadyM.DoryanidaryuniH.MolokwuV. (2022). A Coreflood Simulation Study of a Systematic CO2 Huff-N-Puff Injection for Storage and Enhanced Gas Condensate Recovery. SPE Annual Technical Conference and Exhibition? SPE, D031S059R004.

  • Google Scholar
  • View reference in article

• 39

  ShojaeeA.WangG.GregoryS.MackayE.DopffelN. (2025). New flow simulation framework for underground hydrogen storage modelling considering microbial and geochemical reactions. Int. J. Hydrogen Energy159, 150453. 10.1016/j.ijhydene.2025.150453

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  ShojaeeA.GhanbariS.WangG.MackayE. (2026). Modelling coupled flow with key biogeochemical reactions in subsurface hydrogen systems. Int. J. Hydrogen Energy224, 154478. 10.1016/j.ijhydene.2026.154478

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  WangC.PanL.ZhaoY.ZhangY.ShenW. (2019). Analysis of the pressure-pulse propagation in rock: a new approach to simultaneously determine permeability, porosity, and adsorption capacity. Rock Mech. Rock Eng.52 (11), 4301–4317. 10.1007/s00603-019-01874-w

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  WangY.NolteS.TianZ.Amann‐HildenbrandA.KroossB.WangM. (2021). A modified pulse‐decay approach to simultaneously measure permeability and porosity of tight rocks. Energy Sci. and Eng.9 (12), 2354–2363. 10.1002/ese3.989

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  WidderD. V. (2015). Laplace Transform (PMS-6), 61. Princeton university press.

  • Google Scholar

• 44

  YangZ.SangQ.DongM.ZhangS.LiY.GongH. (2015). A modified pressure-pulse decay method for determining permeabilities of tight reservoir cores. J. Nat. Gas Sci. Eng.27, 236–246. 10.1016/j.jngse.2015.08.058

  • CrossRef
  • Google Scholar

• 45

  ZhangJ. J.LiuH.-H.BoudjatitM. (2020). Matrix permeability measurement from fractured unconventional source-rock samples: method and application. J. Contam. Hydrology233, 103663. 10.1016/j.jconhyd.2020.103663

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  ZhaoY.WangC.ZhangY.LiuQ. (2019). A method of differentiating the early-time and late-time behavior in pressure-pulse decay permeametry. Geofluids.

  • Google Scholar

• 47

  ZhaoY.WangC.BiJ. (2020). A method for simultaneously determining axial permeability and transverse permeability of tight reservoir cores by canister degassing test. Energy Sci. and Eng.8 (4), 1220–1230. 10.1002/ese3.580

  • CrossRef
  • Google Scholar

• 48

  ZhengL.BaiX.XiaoW.ZhuB.ShiJ.RenJ.et al (2025). Advanced Co2 injection in tight oil reservoirs: pore-scale mechanisms from high-temperature, high-pressure visualization and nmr. Energy and Fuels39 (42), 20315–20329. 10.1021/acs.energyfuels.5c04001

  • CrossRef
  • Google Scholar