Abstract
Accurate evaluation of transmissibility is central to reliable numerical simulation of gas flow in porous media, because of its pressure-dependent properties such as density and viscosity. In commercial reservoir simulators, inter-block transmissibility is calculated using upstream averaging and in the well-containing grid-block, it is evaluated at the grid-block pressure. This procedure does not capture the physics of flow because gas properties change gradually one grid-block to another and significantly within the well-containing grid-block. This study presents a novel integral-based formulation in which pressure-dependent properties are integrated over a pressure range which provides exact values to be used in transmissibility terms, without invoking assumptions or mathematical averaging. The proposed method is fully physics-based, robust, and straightforward to implement within finite-volume reservoir simulators. Owing to its general derivation, the method is applicable to natural gas reservoirs, CO2 injection and storage for CCUS applications which involves multiphase flow, and hydrogen subsurface storage. This methodology helps obtain pressure and production values more precisely by avoiding unjustified mathematical averaging. Comparison of the results obtained from the proposed formulation with the ones obtained from commercial softwares shows significant improvement when grid-blocks are relatively large and absolute permeability is low.
1 Introduction
Accurate representation of gas flow in porous media is essential for the numerical simulation of natural gas reservoirs, CO2 sequestration, and underground hydrogen storage (UHS). In finite-volume formulations, transmissibility governs mass transfer between adjacent grid blocks and between the reservoir and wells, therefore plays a central role in the accuracy of the numerical solution. Because gas density and viscosity vary strongly and nonlinearly with pressure, transmissibility becomes a pressure-dependent quantity that must be treated carefully to preserve the physics of gas flow ().
In most reservoir simulators, such as ECLIPSE or CMG, inter-block transmissibility is evaluated using upstream averaging, and in well-containing grid-blocks it is evaluated at the grid-block pressure. It should be noted that pressure does not change discretely from one block to another, but varies continuously along the flow path, and the pressure drop within a well-block is steep. As a result, evaluating gas mobility at a single pressure misrepresents the true physics, especially in low-permeability formations or when grid blocks are large.
The need to rigorously handle pressure-dependent gas properties was first recognized in analytical studies. Al-Hussainy et al. showed that gas-property variation must be treated through integration over pressure rather than by evaluating at a single point, establishing the foundation of real-gas pseudo-pressure theory (), which was explained in detail through numerous applications by (). Subsequent simulation practices provided ways to calculate flow within well-containing grid-blocks for slightly compressible fluids, not for compressible ones (; ).
As reservoir simulation matured, several researchers (; ) extended the treatment of multiphase flow and anisotropy, yet they continued to compute gas mobility at single pressures rather than along pressure gradients. While novel discretization schemes, such as Multi-Point Flux Approximation (MPFA), addressed geometric and anisotropic complexities, they still relied on pointwise mobility sampling for transmissibility evaluations ().
In unconventional reservoirs, the issue becomes more pronounced. Artus and Fructus showed that nonlinear pressure fields near fractures can significantly distort production forecasts if transmissibility is not corrected, yet their correction did not incorporate the explicit pressure dependence of gas viscosity and density within the transmissibility term itself, leaving a key physical gap (). Recent application areas further amplify the need for physically accurate transmissibility. In CO2 geological storage, especially near injection wells, CO2 can approach the critical point, where properties vary sharply with pressure and temperature. Nghiem et al. showed that conventional simulators often struggle in these regimes unless this strong pressure dependence is represented with high fidelity—reinforcing the importance of evaluating mobility along a continuous pressure path rather than at a single discrete value ().
Similarly, in the area of underground hydrogen storage () identified reservoir and operational controls where pressure gradients govern withdrawal performance () benchmarked simulators and showed that hydrogen mixing and recovery efficiency are sensitive to how pressure-dependent mobility is treated; and broad reviews by (; ) emphasize cyclic pressure behaviour, geomechanics, and storage efficiency, all of which demand physically consistent mobility treatment.
() provided direct laboratory measurements showing that the mobility of common reservoir gases (CO2, CH4, N2, and air) varies strongly with pressure, confirming that mobility cannot be reliably represented using a single pressure point and must instead be evaluated across the actual pressure range. These findings reinforce that evaluating mobility at a single pressure is physically inconsistent with real-gas behaviour ().
Multidimensional upwinding analyses show grid-orientation bias and modern near-well modelling demonstrates that pointwise property evaluation can fail when pressure changes rapidly around the well, prompting local corrections (). In particular (), showed that Peaceman-type well-block treatments become inaccurate when strong pressure variations occur near the wellbore, requiring local corrections to capture near-well physics more faithfully. Also (), demonstrated that classical single-point upstream weighting introduces significant grid-orientation effects and systematic bias in multiphase transport, even on structured grids, highlighting that pointwise mobility evaluation can be physically inadequate even before considering strong gas-property nonlinearity. Prior transmissibility corrections in tight/shale contexts provide useful context (), but they did not integrate pressure-dependent gas properties.
This study presents an integral-based formulation for the computation of gas-phase transmissibility that directly accounts for its pressure dependence without relying on assumptions or mathematical averaging, for both grid-block boundaries and within well-containing grid blocks. Because the approach is independent of reservoir geometry or gas type, it is applicable to a broad class of porous media flow problems involving compressible gases in multi-phase flow applications.
2 Mathematical development
In this work, primary production from a gas reservoir is studied in which production is due to gas expansion in the reservoir as the pressure depletes with time. Using a mass balance approach, it is known that . Considering a finite element (with a cross-sectional area of and length of ) in a 1D system, (Equation 1) holds,
Where, is the gas density, is the gas volumetric flowrate and is the porosity of the reservoir rock. Using Darcy’s law and definition of compressibility, the final gas flow equation is obtained shown in Equation 2,
Discretization of the equation above in -direction and multiplying both sides by gives,
The derivation of Equation 3 is presented in Supplementary file 1. An extra equation is required to describe gas flow within the well-containing grid-block.
Where, is well-index and equal to and . The main concern is to obtain an accurate value for in the equations above, which is the primary goal on this study.
Considering two adjacent grid-blocks, gas flows from higher-pressure grid-block to lower pressure grid-block . In numerical simulation, due to discretization, pressure changes suddenly from one grid-block to another. But, in nature this pressure change is gradual. For this purpose, imaginary points are taken between the centers of two grid-blocks, as shown in Figure 1.
FIGURE 1
Now the LHS of Equation 3, i.e., can be re-written in the following form,
Darcy’s law works only when fluid viscosity and compressibility are constant between two points, i.e., it works in a differential form, where viscosity and compressibility are almost constant over a very short distance. Therefore, the RHS of Equation 5 can be expressed in an integral form,
The RHS of Equation 6 is the area under the curve if vs. pressure was plotted, therefore the quantity averaged over a pressure range is obtained in the following way,
2.1 Inter-block transmissibility calculation
Equation 7 is used to obtain at the boundary between each two grid-blocks for all timesteps during the numerical simulation calculation, as shown by the large red dot in Figure 2. The details are presented in Supplementary file 1. It is noteworthy to compare this value with that is obtained at the upstream pressure (Figure 2). The larger the difference, the larger the error that would be introduced into the numerical simulation. Note that the schematic shown in Figure 2 is not to scale and adjacent grid-blocks have much closer pressures in the process of numerical simulation. Also, it is important to know that the area under the red dashed line is equal to the area of the blue dashed area.
FIGURE 2
In Equation 7, pressure values at timestep are used, but in the IMPES formulation used for this study, pressure values at timestep are used. Pressure difference between two grid-blocks varies only very slightly from timestep to . To ensure that, smaller timesteps are chosen during early times of primary production hence, during the numerical simulation run.
2.2 Transmissibility calculation in the well equation
This part is the most important part of this study. Flow of gas from the well-containing grid-block to the wellbore is described by Equation 4. To obtain the transmissibility term, i.e., , the same method explained above is used, as graphically shown in Figure 3. A more detailed step-by-step procedure is presented in Supplementary file 1. There is a noticeable difference between the value the proposed integral method provides and the value that is calculated at the well-containing grid-block pressure. This difference introduces a large error into the numerical simulation. This error is more significant 1) during early times of production where the pressure in the well-containing grid-block is much larger than the wellbore pressure (which is constant during production), 2) in low permeability reservoirs, and 3) when the well is in a large cubic grid-block.
FIGURE 3
Iterative methods update the pressure in the well-containing grid-block. Computing the transmissibility at the updated pressure does not provide a physics-based value–it gives almost the same value. Taking a smaller timestep does not solve the problem either because gas still flows from a high-pressure region in the well-containing grid-block towards the wellbore with a fixed pressure. The proposed integral method accounts for the gradual pressure drop in this grid-block.
2.3 RHS of the mass balance equation (i.e., equation 3)
On the RHS of Equation 3, and are pressure dependent. Should they be calculated at or ? Pressure in the well-containing grid-block does not change much in one timestep. Therefore, these parameters can be calculated at the current timestep. Strictly speaking, pressure drop in this grid-block is significantly sharper during the early times of production. Therefore, it is best to use relatively smaller timesteps during early times of production to avoid very large pressure drop in this grid-block, then these parameters can safely be calculated at the current timestep.
3 Results and discussion
A simulation example is presented here using a sample data taken from (). In Table 1, values of gas formation volume factor and viscosity at different pressures are given, from which is calculated. Best fit equations for these parameters are then obtained, as shown in Equations 8–10,
TABLE 1
| 5000 | 0.000511 | 0.002867 | 0.0381 | 9155 |
| 4500 | 0.000570 | 0.003200 | 0.0364 | 8594 |
| 4000 | 0.000644 | 0.003618 | 0.0346 | 7982 |
| 3500 | 0.000741 | 0.004159 | 0.0329 | 7311 |
| 3000 | 0.000870 | 0.004884 | 0.0312 | 6572 |
| 2500 | 0.001052 | 0.005907 | 0.0294 | 5755 |
| 2000 | 0.001328 | 0.007455 | 0.0277 | 4846 |
| 1500 | 0.001792 | 0.010064 | 0.0259 | 3830 |
Gas formation volume factor and viscosity at different pressures.
Integrating the equation for prior to simulation provides equation that is used to calculate the inter-block and well transmissibility terms during the numerical simulation, shown in Equation 11.
System properties and flow conditions for simulation runs are given in Table 2. Results from the following three different tools are presented and compared: 1) CMG, 2) Eclipse, and 3) The proposed integral equation. CMG and Eclipse employ “upstream averaging scheme” to obtain mobility and viscosity terms at the boundaries between grid-block and compute these parameters at the grid pressure for the well-containing grid-block. We have created an in-house code to numerically simulate gas flow in porous media in which the integral method is incorporated. A detailed procedure in presented in Supplementary file 1.
TABLE 2
| Reservoir properties and flow conditions | Value | Unit |
|---|---|---|
| Length in x-direction, | 5000 | |
| Length in y-direction, | 50 | |
| Length in z-direction, | 50 | |
| Porosity, ϕ | 0.15 | |
| Permeability, | Variable | |
| Production wellbore radius, | 0.2 | |
| Initial pressure in the layer, | 5000 | |
| Production well pressure, | 1500 | |
| Number of grid-blocks, | 1000 | |
| Rock compressibility, | 0 | |
| Timestep size, | 0.1 |
Reservoir properties and flow conditions for the simulation.
Note that this is a dry gas reservoir, i.e., the reservoir is entirely filled with dry gas prior to production. And there will be no condensation around the wellbore during production. Rock compressibility is assumed to be zero to see the effect of gas expansion only.
For the set of data presented above, two identical reservoirs are considered which differed only in the value of absolute permeability. The following important points have been observed.
In both cases, 1000 grid-blocks are used, therefore grid-blocks have very close pressures during the numerical simulation and inter-block transmissibility terms obtained from the proposed integral method are almost the same as those calculated at upstream pressures.
The transmissibility term in the well equation obtained from the proposed integral method is noticeably different from the one computed at the well pressure as long as there is a large difference between the grid pressure and the constant well pressure, as shown graphically in Figure 3. Hence, production profiles given by CMG or Eclipse will be significantly different from the one given by our in-house code. At later times of production, all three tools give very close results.
When the reservoir permeability is higher, reservoir depletes faster, so does the pressure in the well-containing grid-block. Therefore, the three tools give different results only for a short period of time during the early times of production. In a reservoir with a lower permeability, well-containing grid-block pressure declines at a slower rate. During this time there is a significant difference between well pressure and grid pressure, hence a noticeable difference between the results given by CMG or Eclipse and our in-house code.
In large cubic grid-blocks, well-containing pressure takes a longer time to drop, therefore during this time the proposed integral formulation gives more precise results. In very small redial well-containing grid-blocks this integral formulation is not necessary because pressure in the grid-block drops quickly.
Case 1
Gas production profiles for the reservoir with this choice of permeability, obtained from three tools, differ only during very early periods of production, as shown in Figure 4. The difference is due to different treatments of transmissibility terms in the numerical simulation. CMG and Eclipse give identical results which confirms they both employ upstream averaging (identical blue and green curves in Figure 4). The difference disappears with time as the pressure in the well-containing grid-block approaches the well pressure. The cumulative error in the calculation of production becomes negligible compared with the total gas production in 1 year (Figure 5). This necessitates using the integral method when the early-time behavior of the reservoir is studied.
FIGURE 4
FIGURE 5
Case 2
A significantly lower permeability is considered here. At a lower permeability, pressure in the well-containing grid-block drops more slowly, therefore results from CMG/Eclipse and our in-house code will differ more significantly and for a longer time. Results are presented for the first 100 days in Figure 6, which are significantly different. In Figure 7, it is shown that even after 27 years of production these three tools give considerably different results, for the reason presented above. This very low permeability case confirms that a robust and physics-based methodology must be employed to calculate in the equations. It is again worthwhile noting that CMG and Eclipse give identical results.
FIGURE 6
FIGURE 7
4 Conclusion
In nature, pressure changes gradually from one grid-block to another, and in the well-containing grid-block, it changes gradually from the edges of the grid-block to the wellbore. To account for this behavior, a simple integral method is proposed to calculate in transmissibility terms which provides an averaged value over a pressure range either between two grid-blocks or within well-containing grid-block.
It is not necessary to use this method to obtain inter-block transmissibility if there are a large number of grid-blocks in the system due to very close pressures in adjacent grid-blocks. But it is imperative to use this method in the well equation especially when well-containing grid-block is large. When this grid-block is small or a radial discretization is used (with a very small radial grid-block at the centre), can safely be calculated at the well pressure.
This integral method is more necessary to be employed when the permeability is low. Taking smaller timestep or iterative procedures do not help with getting correct results because these are numerics, not physics. The proposed method is fully physics-based, robust, and straightforward to implement within finite-volume reservoir simulators. Owing to its general derivation, the method is applicable to a wide range of gas-dominated reservoir systems, including natural gas reservoirs, CO2 injection and storage for CCUS applications which involve multiphase flow, and hydrogen subsurface storage. The approach provides a unified and reliable framework for gas-flow simulation.
Statements
Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.
Author contributions
SG: Validation, Conceptualization, Writing – review and editing, Writing – original draft, Methodology, Investigation, Visualization. SN: Methodology, Writing – original draft, Data curation. KS: Formal Analysis, Writing – review and editing, Supervision, Validation, Investigation.
Funding
The author(s) declared that financial support was not received for this work and/or its publication.
Conflict of interest
Author SN was employed by the company TotalEnergies.
The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2026.1841365/full#supplementary-material.
References
1
AavatsmarkI. (2002). An introduction to multipoint flux approximations for quadrilateral grids. Comput. Geosci.6, 405–432. 10.1023/a:1021291114475
2
AbunumahO.OgunludeP.GobinaE. (2021). Experimental evaluation of the mobility profile of enhanced oil recovery gases. Adv. Chem. Eng. Sci.11, 154–164. 10.4236/aces.2021.112010
3
AL-HussainyR.RameyJ. R. H.CrawfordP. (1966). The flow of real gases through porous media. Journal Petroleum Technol.18, 624–636. 10.2118/1243-a-pa
4
AliM.IsahA.YekeenN.HassanpouryouzbandA.SarmadivalehM.OkoroaforE. R.et al (2025). Recent progress in underground hydrogen storage. Energy and Environ. Sci.18, 5740–5810. 10.1039/d4ee04564e
5
ArtusV.FructusD. (2012). Transmissibility corrections and grid control for shale gas numerical production forecasts. Oil and Gas Sci. Technology–Revue D’IFP Energies Nouvelles67, 805–821. 10.2516/ogst/2012054
6
AzizK. (1979). Petroleum Reservoir Simulation, 476. London, United Kingdom: Applied Science Publishers.
7
BabuD.OdehA. S. (1989). Productivity of a horizontal well. SPE Reservoir Engineering4, 417–421. 10.2118/18298-PA
8
ChenZ.ZhangY. (2009). Well flow models for various numerical methods. Int. Journal Numerical Analysis Modeling6, 375–388.
9
CuiH. (2018). “Improving peaceman well flow modelling with 1D radial simulation,” in ECMOR XVI-16th European Conference on the Mathematics of Oil Recovery. European Association of Geoscientists and Engineers, 1–9.
10
DakeL. P. (1983). Fundamentals of Reservoir Engineering. Elsevier.
11
Do RosarioR. C. D.DE SouzaG.Amaral SoutoH. (2020). A comparative study of some well-reservoir coupling models in the numerical simulation of oil reservoirs. Int. J. Adv. Eng. Res. Sci.7, 126–148. 10.22161/ijaers.79.16
12
HogewegS.StrobelG.HagemannB. (2022). Benchmark study for the simulation of underground hydrogen storage operations. Comput. Geosci.26, 1367–1378. 10.1007/s10596-022-10163-5
13
KeilegavlenE.KozdonJ. E.MallisonB. T. (2012). Multidimensional upstream weighting for multiphase transport on general grids. Comput. Geosci.16, 1021–1042. 10.1007/s10596-012-9301-7
14
MccainJ. R. W. D. (1972). Properties of Petroleum Fluids.
15
MiocicJ.HeinemannN.EdlmannK.ScafidiJ.MolaeiF.AlcaldeJ. (2023). Underground Hydrogen Storage: A Review.
16
Mohd RazakS.CornelioJ.ChoY.LiuH.-H.VaidyaR.JafarpourB. (2022). Transfer learning with recurrent neural networks for long-term production forecasting in unconventional reservoirs. Spe J.27, 2425–2442. 10.2118/209594-pa
17
NghiemL.GorucuS. E.ShrivastavaV.ZhangG.KohseB. (2025). “Carbon storage simulation near the saturation pressure and temperature, and critical point of CO2,” in SPE Reservoir Simulation Conference. SPE, D031S011R001.
18
OkoroaforE. R.SaltzerS. D.KovscekA. R. (2022). Toward underground hydrogen storage in porous media: reservoir engineering insights. Int. J. Hydrogen Energy47, 33781–33802. 10.1016/j.ijhydene.2022.07.239
19
PeacemanD. W. (1978). Interpretation of well-block pressures in numerical reservoir simulation. Soc. Petroleum Eng. J.18, 183–194. (includes associated paper 6988). 10.2118/6893-pa
20
PletcherJ. L. (2002). Improvements to reservoir material-balance methods. SPE Reserv. Eval. and Eng.5, 49–59. 10.2118/75354-pa
21
VAN PoollenH.BreitenbachE.ThurnauD. (1968). Treatment of individual wells and grids in reservoir modeling. Soc. Petroleum Eng. J.8, 341–346. 10.2118/2022-pa
Summary
Keywords
gas reservoir simulation, transmissibility, pressure-dependent properties, upstream averaging, integral methodology
Citation
Goodarzian S, Naderi S and Sorbie K (2026) An integral method for accurate calculation of the transmissibility in the numerical simulation of gas flow in porous media. Front. Earth Sci. 14:1841365. doi: 10.3389/feart.2026.1841365
Received
28 March 2026
Revised
24 May 2026
Accepted
03 June 2026
Published
27 July 2026
Volume
14 - 2026
Edited by
Mehrdad Vasheghani Farahani, The University of Manchester, United Kingdom
Reviewed by
Kwame Sarkodie, Kwame Nkrumah University of Science and Technology, Ghana
Sina Omrani, The University of Manchester, United Kingdom
Updates
Copyright
© 2026 Goodarzian, Naderi and Sorbie.
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Saleh Goodarzian, sg370@hw.ac.uk
Disclaimer
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.