Compressible gas flow in porous media characterized through a compressibility number and pressure curvature

Compressible gas flow in porous media characterized through a1 compressibility number Π and pressure curvature2 Alejandro F. Visentini1,*, Juan J. Hidalgo1, and Marco Dentz1 3 1Institute of Environmental Assessment and Water Research, Barcelona, Spain4 Corresponding author: *Alejandro F. Visentini (afv[email protected]; [email protected])5 Abstract6 We analyze single-phase, isothermal gas flow in 1-D homogeneous porous media through a7 pressure-squared porous-medium equation (PME). Non-dimensionalizing the PME with a steady8 advective time scale exposes a compressibility number Π that measures relative pressure change9 across the domain and parametrizes a family of globally attracting equilibria. We examine steady10 and transient gas flow behavior across six decades in Π and three inlet boundary flow data, namely,11 pressure, volumetric flux or mass flux, while keeping the pressure fixed at the outlet. At steady12 state, we derive explicit relations linking different inlet data to a target Π, thereby selecting the13 equilibrium pressure profile, the latter featuring a boundary layer of width ∝Π−1, toward the14 low-pressure outlet boundary, where the volumetric flux overshoots and the gas expands at a15 rate ∝Π2. For transients, high-resolution numerical simulations provide diffusion times based16 on the L2-norm decay of the pressure toward its steady-state. We compare these observed times17 against the slowest-decaying mode time of the eigenfunction expansion of the diffusion operator,18 evaluated using a time–dependent effective diffusivity obtained as the harmonic mean of the local19 diffusion coefficient. Close agreement identifies regimes where this scalar reduction captures the20 variable–coefficient dynamics while systematic deviations quantify when spatial variability signifi-21 cantly alters the path to the attracting equilibrium and a constant–coefficient surrogate fails. We22 finally assess transient gas compression via pressure curvature. Differentiating the PME we obtain23 an advection–diffusion–reaction (ADR) evolution for the curvature that explains the emergence,24 migration, and sharpening of a curvature front. The proposed Π–curvature framework provides a25 compact and operational lens to characterize compressible gas flows.26 Keywords Flow compression, gas flow, porous medium equation27 1 Introduction28 Gas flow through porous media occurs in countless natural and industrial settings, including volcanic29 degassing (e.g., Massol and Jaupart, 1999; Oppenheimer et al., 2013), lake sediment venting (e.g.,30 1 Scandella et al., 2016), underground storage of natural gas (e.g., Tek, 1989), hydrogen (e.g., Heine-31 mann et al., 2021) and carbon dioxide (e.g., Bachu, 2008), air-based groundwater remediation (e.g.,32 Hunt et al., 1988; Hu et al., 2010), and operation of membrane fuel cells (e.g. Lee et al., 2019). In most33 such settings, gas compressibility strongly affects transient and steady state flow behavior, leading to34 a nonlinear relationship between the gas injection or pumping rate and the spatio-temporal evolution35 of pressure and flow velocity. A notable example is large-scale gas reservoir operations, where gas is36 periodically injected and pumped at different rates and through different prescribed boundary con-37 ditions, namely, pressure, volumetric flux or mass flux, which causes qualitatively distinct transient38 paths of the reservoir’s pressure toward its steady-state. How such paths differ as a function of the39 type and intensity of the forcing through gas compression is not well characterized. Further, there is40 no single operational measure quantifying the relevance of compression effects at the reservoir-scale.41 Gas flow in porous media is generally described by combining mass and momentum (e.g., Darcy’s42 law) conservation equations, and an equation of state (e.g., Liepmann and Roshko, 2001). The resulting43 system is highly nonlinear, admitting analytical treatments only for few special cases (e.g., Bear, 1972).44 This is because of the many gas flow parameters depending on pressure, namely, gas density, viscosity,45 temperature and compressibility factor, the latter measuring departure of the gas from ideal behavior46 (e.g., Chen et al., 2006; Chen, 2007). However, for many practical cases, the pressure dependence of the47 gas density dominates over the other parameters, which can be safely assumed constant. For such case,48 the governing system reduces to a porous-medium equation (PME), written in terms of the pressure49 squared and featuring a pressure-dependent diffusivity (e.g., Al-Hussainy et al., 1966; V´azquez, 2007).50 The PME is well known across disciplines concerned with the study of gas flow, including hydroge-51 ology (e.g., Baehr and Joss, 1995), low (tight) permeable media, (e.g., Wu et al., 1998; Monteiro et al.,52 2012; Wu et al., 2014; Zhao et al., 2016), fluid mechanics (e.g., De Socio and Marino, 2006), petroleum53 engineering (e.g., Chierici, 2012) and pure mathematics (e.g., V´azquez, 2007). Notable advances come,54 on the one hand, from petroleum engineering literature, that considered the PME early on (e.g. Muskat,55 1934, 1937), in combination with transformations such as the pseudo-pressure or pseudo-time to model56 the flows of real gases (e.g., Al-Hussainy et al., 1966; Aanonsen, 1985). Here, the PME is typically57 linearized at some representative value of gas diffusivity (e.g., at mean pressure), enabling analyti-58 cal treatment of well tests interpretation, average pressure forecasting and pressure buildup analyses,59 among other numerous applications (e.g., Matthews et al., 1967; Tek, 1989; Chierici, 2012; Patzek60 et al., 2013; Ahmed, 2018). However, these approaches are typically not concerned with resolving the61 detailed spatial structure of the pressure field and its evolution, but focused on aggregated reservoir62 behavior, and, consequently, do not explicitly address the nonlinear impact of gas compressibility in63 the flow behavior. On the other hand, the mathematical analysis of non-linear diffusion equations64 continuously investigates the PME from a theoretical viewpoint, reporting classical results concerning65 finite-speed propagation (e.g., Martinson, 1976), self-similar solutions (e.g., Aronson and Graveleau,66 1993; Barenblatt, 2003), degenerate diffusion (e.g., Burgers, 2013), geometrical properties of solutions67 of the PME (e.g., Lee and V´azquez, 2003), and regularity of free boundaries (e.g., Daskalopoulos and68 2 Hamilton, 1998). However, these analyses are rarely connected to operational contexts of applied fluid69 mechanics or subsurface engineering, offering little practical guidance or physical insights on how the70 global compression state affects gas flow evolution as a function of input boundary flow controls. De-71 spite extensive work across petroleum, fluid mechanics and PME-specialized mathematical literature,72 there is, to our knowledge, no framework systematically quantifying the impact of gas compression73 effects on emergent gas flow behavior across operational boundary conditions at the reservoir scale.74 Here, to progress toward the goal of having a global framework to interpret gas flow under variable75 compression effects and boundary conditions, we introduce a compressibility number, Π, which non-76 dimensionalises the pressure-squared PME and measures the strength of gas compression across the77 domain via the relative pressure change between the initial and steady states. The parameter Π78 is conceptually related to the compressibility number recently introduced by Cuttle and MacMinn79 (2023), used to model gas–liquid displacement in a capillary tube and measuring gas compression80 relative to viscous dissipation effects. In contrast, the present work focuses on single-phase gas flow at81 the continuum scale, where compressibility affects transient storage and pressure diffusion rather than82 interfacial displacement. The parameter Π provides a single yardstick linking steady and transient83 responses across boundary conditions. We first revisit steady solutions through the lens of Π to expose84 an outlet boundary layer with width ∝Π−1where the gas increasingly expands at a rate proportional85 to Π2, and we derive explicit inlet–data relations that select a targeted Π (i.e. a specific attracting86 equilibrium). We then examine transients with high-resolution MRST simulations (Lie, 2019; Lie87 and Møyner, 2021), extracting diffusion times from the L2-norm decay of pressure toward its steady-88 state, and comparing them against the slowest-decaying mode time corresponding to the eigenfunction89 expansion of the diffusion operator evaluated using a time–dependent effective diffusivity, obtained90 as the harmonic mean of the (spatio-temporally variable) local diffusion coefficient.The comparison91 delineates regimes where a constant-coefficient surrogate is adequate and where spatial variability in the92 pressure-dependent diffusion coefficient reshapes the approach to steady-state. We finally use pressure93 curvature as a diagnostic of flow compression. Differentiating the PME twice with respect to space94 we obtain a curvature evolution equation with the form of an advection-diffusion-reaction equation95 (ADRE) that explains the emergence, downstream migration, and sharpening of a curvature front.96 We track the location, amplitude and integrated curvature associated with this front, demonstrating97 how these metrics provide compact diagnostics for identifying gas compression and expansion flow98 regimes.99 The article is organized as follows. In Section 2 we present the governing equations for gas flow,100 including the porous medium equation (PME). In Section 3 we introduce the compressibility number101 used to non-dimensionalize the PME and the pressure curvature evolution equation used as diagnostic102 of compression regimes. Section 5 analyzes steadyand transient gas flow behavior. Section 6 concludes103 the article.104 3 2 Governing equations105 The laminar flow of real gases in porous media is described by the following set of equations with106 adequate initial and boundary conditions:107 ∂ϕρ ∂t =−∇·(ρq),(1) q=−k µ∇p−ρgz,(2) ρ=pW RTZ .(3) Equation 1 is a statement of fluid mass conservation in the absence of sources or sinks, with ϕ(x, t) the108 porosity ρ(x, t) the fluid mass density (kg m−3) and q(x, t) the (volumetric) specific discharge (m s−1).109 The fluid velocity u(x, t) = q(x, t)/ϕ(x, t). The product of ρqdefines the mass flux, qm=ρq(in110 kg m−2s−1). x= (xx, zz)Trepresents the position vector in the 2-D domains used here. Equation 2,111 Darcy’s law, is a statement of momentum conservation at the continuum scale, with k(x), µ(x, t),112 p(x, t) and gdenoting, respectively, the intrinsic permeability of the medium (in m2), the dynamic113 viscosity of the fluid (in Pa s), the phase pressure (in Pa) and the constant of gravity acceleration114 (in m s−2). For gases, the permeability is in general dependent on the pressure (Klinkenberg, 1941),115 although here we consider it constant given the large confining pressures (e.g., Wu et al., 1998).116 The equation of state 3 relates the gas density with its pressure, via the (average) molecular weight117 of the pure gas (mixture) W(kg mol−1), the universal R= 8.314 J mol−1K−1and the compressibility118 factor Z(p), the latter measuring the deviation of a real gas (or mixture) from ideal behavior. Here,119 we use the Peng-Robinson relationship (Peng and Robinson, 1976) for establishing the value of Z(p)120 as a function of the gas properties. This is detailed in Appendix A.121 2.1 Gas flow equation122 Combination of Equations 1-3 and a series of manipulations and assumptions, detailed in Appendix A,123 leads to the non-linear partial differential equation on the pressure-squared (e.g., Al-Hussainy et al.,124 1966; Chierici, 2012) describing gas flow in homogeneous porous media:125 ∂p2 ∂t =D∇2p2,(4) where the non-linearity stems from the pressure dependence of the diffusion coefficient126 D=k ϕµγ ,(5) 4 via its thermodynamic compressibility γggiven in isothermal conditions as:127 γg=1 ρ ∂ρ ∂pT=T0 (6) =1 p−1 Z ∂Z ∂p . Equation 4 assumes negligible gravity effects and that the product µZ is approximately constant128 within the pressure range considered (e.g., Al-Hussainy et al., 1966). As we show in Subsection A,129 the second assumption holds for pressures below 15 MPa for hydrogen. Note that Equation 4 is an130 instantiation of the Porous Medium Equation (PME) for m= 2 (e.g., V´azquez, 2007).131 3 Compressibility number Π132 We now non-dimensionalize the PME Equation 4 as follows:133 ∂p2 ∂t =D∇2p2,(7) where ·denotes dimensionless variables or operators, defined as:134 ∇2=L2∇2, p2=p2/p2 L,x=x/L, t =t/τ∞ adv,(8) with La characteristic spatial scale, taken as the horizontal length of the reservoir, and the char-135 acteristic pressure-squared p2 Ltaken as the initial mean reservoir pressure. The advective time-scale136 at steady-state τ∞ adv =u∞/L is given in terms of the characteristic flow velocity at steady-state,137 u∞= (k/µ)(∆∞p/L), where ∆p∞=p∞ 0−pLis the pressure difference between inlet and outlet at138 steady-state. From now on, we drop ·since we work in non-dimensional space throughout.139 The non-dimensional diffusion coefficient Dis given as:140 D(x, t) = p(x, t) Π,(9) where the parameter Π = τD/τ∞ adv (= u∞L/D) is the P´eclet number for pressure, comparing τ∞ adv 141 against the characteristic time for pressure to diffuse across L,τD=L2/D.142 We show now that Π can be used as a measure of the degree of compressibility effects experienced143 by the system at steady-state. Indeed, substituting u∞= (k/µ)(∆∞p/L) and D=k/γµϕ , one144 obtains:145 Π=∆∞pγ. (10) Since the intrinsic thermodynamic compressibility of the gas γhas units of inverse pressure, Π expresses146 a relative pressure (or density) change at steady-state, hence giving an average compression degree147 at that state. Small (<1) and large (>1) values of Π indicate mild and strong compression effects,148 5 respectively. Consequently, we label Π as a compressibility number, which is to be distinguished from149 γ(Eq. 6). In Subsection 5.1.1 we give operational definition for Π as a function of applied inlet flow150 data.151 3.1 Temporal re-scaling using effective diffusivity for 1-D media152 We consider a temporal re-scaling accounting for time-varying diffusivity as follows (e.g., Crank, 1975):153 t(τ) = Zτ 0 Deq(s)ds, (11) where the time-varying equivalent diffusivity Deq(t) is calculated as154 Deq(t) = 1 ⟨R(x, t)⟩x ,(12) with ⟨R(x, t)⟩xthe (spatial) arithmetic mean of the resistivity R(x, t), the inverse of diffusivity D(x, t) =155 p(x, t)/Π. Time-series and transient results are considered in terms of the diffusion-scaled time variable156 t(τ).157 4 Evolution equation for the pressure curvature κ=∂xxp158 We now derive an evolution equation for the pressure curvature κ:=∂xxp, which is zero in the159 incompressible flow limit and increases in amplitude with the degree of flow compression. Further, the160 sign of κindicates either, ongoing expansion, for κ < 0 (i.e., ∂xqv>0), or ongoing gas compression,161 for κ > 0 (i.e., ∂xqv<0). Thus, we use κas a diagnostic (signed) scalar of the compression state of162 the gas over the domain. Using the fact that ∂t(p2)=2p ∂tpand ∂2 x(p2)=2(px)2+p pxx, we re-write163 Equation (7) as164 pt=1 Π(px)2+p pxx.(13) Differentiating twice with respect to x, we obtain165 κt=4px Πκx+p Πκxx +3 Πκ2.(14) Equation 14 has the form of an advective diffusion reaction equation (ADRE), thereby featuring three166 distinct mechanisms that govern the transient evolution of κ. First, an advection term 4px Πκx, with167 velocity uκgiven as168 uκ=−4px Π,(15) 6 a diffusion term p Πκxx, with coefficient Dκgiven as:169 Dκ=p Π,(16) which has the same form as in Equation 7, and a nonlinear reaction term,170 Rκ= (3/Π) κ2(17) (self-)amplifying κ= 0. In Section 5.2 we use Equation 14 to interpret transient evolution under171 different boundary conditions.172 Since κbehaves as a signed transported scalar (e.g., Acheson, 1990; Saffman, 1992), we consider173 the evolution of the zeroth-order moment (e.g., Taylor, 1950; Aris, 1956) of its positive and negative174 parts, respectively, κ+(x, t) and κ−(x, t), to characterize emergence and decay of compression and175 expansion regimes across the system:176 M + − i(t):=Zdκ+ −.(18) We also consider the point of maximum κ+(x, t):177 xpeak(t):= arg max κ+(x, t),(19) which we use to track the compression front migrating toward the outlet while the system converges178 to its steady-state. We do not consider moments of order >0 (e.g., a center of mass or spread)179 because curvature mass is not conserved but lost early on toward the outlet boundary, leading to180 systematic biases of higher-order moments. This is not an artifact but an intended feature of our181 setting, designed to investigate gas flow behavior that converges toward a boundary-modulated steady-182 state (Subsec. 5.1). A moment-based analysis of curvature evolution away from the influence of outlet183 boundaries remains the topic of further research.184 Finally, we emphasize that, unlike classical porous medium flows with compact support, where185 the interface between regions of zero and positive ”mathematical” pressure constitutes a moving free186 boundary (e.g., Daskalopoulos and Hamilton, 1998; V´azquez, 2007), the present problem features187 p(x, t)>0 throughout and, consequently, it remains smooth. As a result, classical spatial derivatives188 up to fourth order are well defined, permitting formal derivation of Equation (14). We also stress that189 Equation 14 is a diagnostic evolution law, that is, derived under assumed knowledge of pand px, which190 are predicted from the governing equation of the system, Equation 7.191 7 5 Analysis of gas flow in terms of the compressibility number192 Πand pressure curvature κ193 We focus on 1-D linear gas flow over a domain of length L= 1, since our goal is to understand the194 relationship between imposed inlet boundary data and emergent flow behavior via the Porous Medium195 Equation (PME) across different Π-values. Thus, at x= 0 we will apply either pressure p0, volumetric196 flux qv 0or mass flux qm 0, while at x= 1 we fix the pressure at p1= 1 throughout.197 5.1 Steady-state behavior198 We start by examining steady-state behavior. Setting ∂p2/∂t = 0 in Equation 7 the governing Laplace199 equation in 1-D reads:200 d2p2 dx2(x)=0.(20) Integrating Equation 20 once and using the fact that Π = ∇p2L/2pLγ(Eq. ??), we obtain:201 dp2 dx =−2Π,(21) Integrating again, the affine-form solution in terms of Π reads as:202 p(x) = p2Π(1 −x) + 1,(22) and the corresponding expressions for the non-dimensional volumetric and mass flux, respectively,203 qv(x) and qm(x), read:204 qv(x)=−Π p(x),(23) qm(x)=Π.(24) where we have used the equality between non-dimensional pressure and density, that is, p=ρ, we have205 set the characteristic volumetric and mass flux scales as qv c= (k/µ)pL/L and qm c= (WKp2 L)/(RTµL),206 respectively (Subsec. ??). Differentiating Equation 23 we obtain the spatial gradient of qvor, equiva-207 lently, the pressure curvature κ(x):208 κ(x)=−Π2/p3(x).(25) Figure 1 features profiles of non-dimensional pressure, flow velocity and modulus of pressure cur-209 vature for Π ∈[0.01,0.1,1,10,100].210 8 Figure 1: Profiles of (a) pressure p, (b) flow velocity qvand (c) modulus of pressure curvature κfor Π-values of (light to dark blue) (0.01,0.1,1,10,100). The pressure profiles feature a concave shape, with near-constant gradient toward the inlet and a211 steeper gradient approaching the outlet (Fig. 1a). Correspondingly, a near-uniform flow velocity qv 212 (Fig. 1b) is observed for x≪1, and an overshoot develops as x→1. The spatial heterogeneity in213 qvaccentuates with increasing Π and toward the outlet, which is reflected in a negative curvature κ214 of increasing amplitude (Fig. 1c) toward the outlet. Physically, it reflects that gas parcels undergo215 increasing volumetric expansion as they flow downstream, with qvincreasing toward the outlet to216 conserve the constant mass flux qm=ρqvover the profile at steady-state while the pressure decays217 due to the imposed outlet condition. To understand how the gas expansion increases and distributes218 over space, we study the curvature scaling for the limit Π ≫1. We partition the domain into an219 interior region [0,1−δBL,] and a thin outlet boundary layer (BL) [1 −δBL,1] of width δBL. We define220 δBL by requiring that the pressure at the inner edge of the layer remains O(1) relative to its outlet221 value (e.g., Bender and Orszag, 2013), that is:222 p(1 −δBL) = c, (26) 9 becomes Mκ− ≈Mκ+(colored dots), indicating the transition from compressionto an expansion-313 dominated steady-state. Imposing fluxes, that is, px(red and green, Fig. 4a), leads to a growth of Mκ+ 314 ∝t1via the reaction term, until t∼10−5. A delay is also observed before the curvature peak begins315 its downstream motion (red and green, Fig. 4c) because, despite the advection term being active from316 t= 0, via imposing px, at such early stage the curvature field is still being formed by the reaction317 term so there is little κ+ xor κ+ xx to advect or diffuse, respectively. Once a sufficiently sharp κ+lobe318 has developed the advective term kicks in and xmax κ+evolves ∝t1. For the case Π = 100, Mκ+grows319 at a non-monotonic rate, approximately ∝t1for t<10−6and ∝t0.5for 10−6<t<10−5, maximizing320 again at t∼10−5and decaying ∝t0.5until t∼1, when the transition to negative curvature occurs.321 Non-negligible negative curvature mass is observed well before t∼1 (dashed line in Fig. 4b), indicating322 gas expansion regions during the transient period. Higher inlet pressure leads to stronger curvature323 diffusing rapidly away from the inlet, since the diffusion coefficient increases with pressure. At the324 same time, such an incoming curvature generated upstream must be accommodated at the outlet to325 enforce the Dirichlet boundary condition. For small Π this process is weaker and delayed. Using the326 signed curvature field and associated PDE (Eq. 14) provides a complementary interpretative frame-327 work to understand how gas compression and expansion distributes over space and time, for instance328 disentangling the concurrent effects of gas compression and expansion, which are not directly apparent329 using mass conservation and pressure-based analysis only.330 5.6 Impact of spatially-variable pressure diffusivity on gas flow331 We now examine the impact on pressure p(x, t) evolution arising from spatial variability in the diffusion332 coefficient D(x, t)=p(x, t)/Π. For this we consider distributions of pressure calculated analytically,333 p2 calc(x, t), through the classical solution of the linearized PME (Eq. 7), which is exactly valid for the334 case p∼1 (e.g., Carslaw and Jaeger, 1959; Crank, 1975):335 p2 calc(x, t) = p2 ss(x) + X n Anϕn(x) exp −λnt,(39) where p2 ss(x) is the steady-state p2distribution and the second term is the eigenfunction expansion336 of the residual satisfying homogeneous boundary conditions. The coefficients Anare the projections337 of such residual onto the Laplacian eigenfunctions ϕn(x), with BC-dependent eigenvalues λngiven338 by λn= (nπ/L)2and λn= [(n+ 1/2)π/L)]2for imposed inlet pressure or flux, respectively. For339 more details please refer to Subsection A. Note that the variable tcorresponds to time scaled by the340 equivalent diffusivity Deq(t) (Eq. 11). Consequently, p2 calc(x, t) accounts for temporal but not spatial341 variability in diffusivity D(x, t). Figure 5 features p2 calc(x, t) (dashed lines) against the simulated342 p2-profiles.343 16 Figure 5: Pressure profiles as a function of distance xat three different diffusion-scaled times t(Eq. 11) equal to (dark to light) 0.1, 0.2 and 1, for (left column) Π = 0.01 and (right column) Π = 100 for imposed (first row, blue) pressure, (second row, green) volumetric flux and (third row, red) mass flux. Solid and dashed lines indicate pressure profiles calculated numerically and analytically, respectively. The numerically-calculated profiles evolve according to a spatiallyand temporally-variable diffusion coefficient D(x, t) = p(x, t)/Π, while the analytically-calculated ones evolve according to a time-variable equivalent diffusion coefficient Deq(t) given as the harmonic mean of D(x, t). As expected, the discrepancy between the simulated p2-profiles and p2 calc(x, t) is negligible for344 Π=0.01, but increases for Π = 100. Further, the error behaves differently depending on the applied345 inlet flow data. The analytical solution predominantly underand overestimates the true distribution346 when applying pressure (blue) or volumetric flux (green), respectively, while for applied mass flux347 17 (red), the bias is predominantly positive at early times and features mixed positive and negative biases348 toward later times.349 To quantify the aggregated error in the analytical solution due to unaccounted spatial heterogeneity350 in the diffusivity D(x, t), we consider the time-averaged of the residual energy, µe(t), defined as:351 µe(t) = 1 τZτ 0 e∆u(t),(40) where the energy e∆u(t) of the residual difference ∆u(x, t)=u(x, t)−ucalc(x, t) is given by (e.g.,352 Evans, 2022):353 e∆u(t) = 1 2Z1 0 [u(x, t)−ucalc(x, t)]2dx, (41) with the residuals defined as u(x, t)=p2(x, t)−p2 ss(x, t) and ucalc(x, t)=p2 calc(x, t)−p2 ss(x, t).354 The variable τindicates a generic time over which the integration is performed, selected such that355 e∆u(τ)/e∆u(0) < ϵ, where we choose ϵ= 0.01. Assuming identical steady-states Equation 41 reads:356 e∆u(t) = 1 2Z1 0 p2(x, t)−p2 calc(x, t)dx. (42) Figure 6 reports µe(t) as a function of Π and applied inlet boundary flow data.357 Figure 6: Time-average µe(t) of the residual energy e(t) between pressure distributions calculated numerically (p2(x, t)) and analytically (p2 calc(x, t)), the former accounting for spatiallyand temporally-variable diffusivity D(x, t), and the later using a time-variable equivalent diffusivity Deq(t), obtained as the harmonic mean of D(x, t). The average is performed over the interval [0, τ], with τdefined in the main text. The black line indicates scalings. The time-averaged residual energy increases with approximately ∝Π3/2, indicating more cumulated358 error due to neglected spatial variability in diffusivity D(x, t) with increasing Π.359 As an additional diagnostic of cumulated error, we consider the characteristic diffusion time-scale360 18 needed for the pressure to reach its steady-state, τobs D, estimated from the numerical simulations as361 follows:362 τobs D:= inf nt≥0 : eu(t) eu(0) ≤ϵo(43) with eu(t)=1/2R1 0u(x, t)dx, again choosing ϵ= 0.01. That is, τobs Dmarks the minimum time at363 which the energy residual, normalized by its initial (maximum) value, is at a l2-distance smaller than364 ϵ. We compare it against a calculated diffusion time-scale τcalc D, satisfying the same criteria given365 by Equation 43 but using eucalc (t)=1/2R1 0ucalc(x, t)dx. Considering the analytical expression of366 p2 calc(x, t) via (Eq. 39), it is straightforward to show that:367 t(τcalc D)≈ln (1/ϵ) λ,(44) which corresponds to the time for which the slowest-decaying mode of p2 calc(x, t), exp −λt(τ)(Eq. 39),368 is ≤ϵand λcorresponds to the lowest-order eigenvalue of the the expansion, which is λ1=π2 369 and λ0= (π/2)2for imposed pressure or flux, respectively. A more detailed derivation is given in370 Subsection A.371 Figure 7 reports τobs Dand τcalc Dnormalized by tlinear D= Π, the latter obtained by setting p= 1 in372 Equation 9.373 Figure 7: Characteristic times τDfor gas pressure to diffuse across a domain of length L= 1, normalized by the reference diffusive time τpL D= (kpL/µϕ)−1at initial domain pressure pL, as a function of Π. Two different diffusive times are reported, (dots) τobs D, calculated empirically from numerical simulations, satisfying criteria of Equation 43 and (dashed) τcalc D, calculated such that it satisfies criteria given by Equation 44, based on the slowest-decaying mode of the eigenfuction expansion of the linear diffusion operator (Subsec. B). For small targeted steady-state compressibility, that is, Π ≪1, the predicted τcalc Dsits at 0.467τpL D 374 for applied inlet pressure p0(dashed blue, Fig. 7) and at 1.866tpL Dfor imposed volumetric qv 0or mass flux375 qm 0(dashed green and red). Indeed, for Π ≪1 the diffusivity Dis D≈DpL, with DpL=kpL/µ and so376 19 from the definition of the diffusivity re-scaled time (Eq. 44) we have that t(τ)≈τDpL. Combining this377 with t(τcalc D) = ln(1/ϵ)/λ, we have τcalc D= ln(1/ϵ)/DpLλ, yielding then τcalc D/τpL D≈ln(1/ϵ)/λ, which is378 equal to 0.467 and 1.866 when λ=π2(imposed p0)orλ= (π/2)2(imposed qv 0or qm 0), respectively.379 For Π ≲102, a good agreement is observed between τcalc Dand τobs Dfor the cases of imposed (blue)380 pressure p0and (red) mass flux qm 0, indicating that the slowest–mode formula using a time-variable381 equivalent diffusivity Deq(t) can be used to predict diffusion time-scales over a domain with spatially382 variable local diffusivity D(x, t). As Π >102, stronger pressure gradients lead to highly non-uniform383 D(x, t) and the slowest–mode estimate fed by Deq(t) increasingly over predicts τobs D. For the case of384 imposed (green) volumetric flux qv 0the discrepancy between τcalc Dand τobs Dremains negligible until385 Π∼1, it increases and decays again for Π >600. Such non-monotonic dependence on Π highlights the386 nonlinear dependence of the pressure evolution with the applied injection, which appears accentuated387 for the qv 0-case. This remains to be further investigated.388 6 Conclusions389 We developed a compact framework for compressible, single-phase gas flow in homogeneous porous me-390 dia that (i) casts the dynamics in a pressure-squared porous medium equation (PME), (ii) introduces391 a compressibility number Π that provides a single scale for comparing flow regimes across different392 inlet flow data, namely, pressure p0, volumetric flux qv 0or mass flux qm 0, and (iii) uses pressure cur-393 vature as a diagnostic of gas compression. At steady state, closed-form solutions in terms of Π lead394 to simple operational maps from a targeted Π to inlet controls and show that compressibility concen-395 trates the pressure drop into an outlet boundary layer of thickness δBL ∼Π−1, with curvature scaling396 |κ|∼Π2, and integrated curvature growing linearly with Π; these scalings quantify volumetric-flux397 overshoot. For transients, we measure the time to steady by tracking the decay of the L2-norm of the398 deviation of p2from its steady state, and we compare these measured times against predicted times399 based on the slowest-decaying mode fed by a scalar effective diffusivity (the harmonic mean of the400 field at each time). The comparison is accurate over wide ranges of Π and deviates in regimes where401 spatial variability in D(x, t)∝p(x, t) localizes gradients and invalidates a constant-coefficient surro-402 gate. Differentiating the governing PME yields an advection–diffusion–reaction law for the curvature,403 which explains the emergence, sharpening, and downstream migration of a curvature front; tracking404 its position, amplitude, and mass provides concise diagnostics of the transient. High-resolution MRST405 simulations across several decades of Π corroborate the steady scalings, the operational maps, and406 the conditions under which the slowest-mode reduction suffices. The proposed Π–curvature frame-407 work offers insights of compressible gas flows and suggests extensions to heterogeneous media, higher408 dimensions, and non-ideal gas properties, as well as opportunities to use curvature-based metrics for409 control and optimization of injection–withdrawal protocols.410 20 References411 Aanonsen, S. (1985), Application of pseudotime to estimate average reservoir pressure, in SPE Annual412 Technical Conference and Exhibition?, pp. SPE–14,256, SPE.413 Acheson, D. J. (1990), Elementary fluid dynamics, Oxford University Press.414 Ahmed, T. (2018), Reservoir Engineering Handbook, Gulf Professional Publishing.415 Al-Hussainy, R., H. Ramey Jr, and P. Crawford (1966), The flow of real gases through porous media,416 Journal of Petroleum Technology,18(05), 624–636.417 Ambrosio, L., N. Gigli, and G. Savar´e (2005), Gradient flows: in metric spaces and in the space of418 probability measures, Springer.419 Aris, R. (1956), On the dispersion of a solute in a fluid flowing through a tube, Proceedings of the420 Royal Society of London. Series A. Mathematical and Physical Sciences,235 (1200), 67–77.421 Aronson, D. G., and J. Graveleau (1993), A selfsimilar solution to the focusing problem for the porous422 medium equation, European Journal of Applied Mathematics,4(1), 65–81.423 Bachu, S. (2008), Co2 storage in geological media: Role, means, status and barriers to deployment,424 Progress in energy and combustion science,34 (2), 254–273.425 Baehr, A. L., and C. J. Joss (1995), An updated model of induced airflow in the unsaturated zone,426 Water Resources Research,31 (2), 417–421.427 Barenblatt, G. I. (2003), Scaling, vol. 34, Cambridge University Press.428 Bear, J. (1972), Dynamics of Fluids in Porous Media, Dover Publications.429 Bender, C. M., and S. A. Orszag (2013), Advanced mathematical methods for scientists and engineers430 I: Asymptotic methods and perturbation theory, Springer Science & Business Media.431 Burgers, J. M. (2013), The nonlinear diffusion equation: asymptotic solutions and statistical problems,432 Springer Science & Business Media.433 Carslaw, H., and J. Jaeger (1959), Conduction of heat in solids, clarendon.434 Chen, Z. (2007), Reservoir simulation: mathematical techniques in oil recovery, SIAM.435 Chen, Z., G. Huan, and Y. Ma (2006), Computational methods for multiphase flows in porous media,436 SIAM.437 Chierici, G. L. (2012), Principles of Petroleum Reservoir Engineering: Volume 2, vol. 2, Springer438 Science & Business Media.439 Crank, J. (1975), The mathematics of diffusion. oxford, England: Clarendon.440 21 Cuttle, C., and C. W. MacMinn (2023), Dynamics of compression-driven gas-liquid displacement in a441 capillary tube, Physical Review Letters,130 (11), 114,001.442 Daskalopoulos, P., and R. Hamilton (1998), Regularity of the free boundary for the porous medium443 equation, Journal of the American mathematical society,11(4), 899–965.444 De Socio, L., and L. Marino (2006), Gas flow in a permeable medium, Journal of Fluid Mechanics,445 557, 119–133.446 Evans, L. C. (2022), Partial differential equations, vol. 19, American mathematical society.447 Heinemann, N., J. Alcalde, J. M. Miocic, S. J. Hangx, J. Kallmeyer, C. Ostertag-Henning, A. Has-448 sanpouryouzband, E. M. Thaysen, G. J. Strobel, C. Schmidt-Hattenberger, et al. (2021), Enabling449 large-scale hydrogen storage in porous media–the scientific challenges, Energy & Environmental450 Science,14(2), 853–864.451 Hu, L., X. Wu, Y. Liu, J. N. Meegoda, and S. Gao (2010), Physical modeling of air flow during air452 sparging remediation, Environmental science & technology,44(10), 3883–3888.453 Hunt, J. R., N. Sitar, and K. S. Udell (1988), Nonaqueous phase liquid transport and cleanup: 1.454 analysis of mechanisms, Water resources research,24(8), 1247–1258.455 Klinkenberg, L. (1941), The permeability of porous media to liquids and gases, Drilling and Production456 Practice, pp. 200–213.457 Lee, C., B. Zhao, R. Abouatallah, R. Wang, and A. Bazylak (2019), Compressible-gas invasion into458 liquid-saturated porous media: application to polymer-electrolyte-membrane electrolyzers, Physical459 Review Applied,11(5), 054,029.460 Lee, K.-A., and J. L. V´azquez (2003), Geometrical properties of solutions of the porous medium461 equation for large times, Indiana University Mathematics Journal, pp. 991–1016.462 Lie, K.-A. (2019), An introduction to reservoir simulation using MATLAB/GNU Octave: User guide463 for the MATLAB Reservoir Simulation Toolbox (MRST), Cambridge University Press.464 Lie, K.-A., and O. Møyner (2021), Advanced modelling with the MATLAB reservoir simulation toolbox,465 Cambridge University Press.466 Liepmann, H. W., and A. Roshko (2001), Elements of Gas Dynamics, Courier Corporation.467 Martinson, L. K. (1976), Finite velocity of propagation of thermal perturbations in media with constant468 heat conduction coefficient, Zhurnal Vychislitelnoi Matematiki i Matematicheskoi Fiziki,16, 1233–469 1241.470 Massol, H., and C. Jaupart (1999), The generation of gas overpressure in volcanic eruptions, Earth471 and Planetary Science Letters,166(1-2), 57–70.472 22 Matthews, C. S., D. G. Russell, et al. (1967), Pressure Buildup and Flow Tests in Wells, vol. 1, Henry473 L. Doherty Memorial Fund of AIME New York.474 Monteiro, P. J., C. H. Rycroft, and G. I. Barenblatt (2012), A mathematical model of fluid and gas475 flow in nanoporous media, Proceedings of the National Academy of Sciences,109(50), 20,309–20,313.476 Muskat, M. (1934), The flow of compressible fluids through porous media and some problems in heat477 conduction, Physics,5(3), 71–94.478 Muskat, M. (1937), Use of data oil the build-up of bottom-hole pressures, Transactions of the AIME,479 123(01), 44–48.480 Oppenheimer, C., T. P. Fischer, and B. Scaillet (2013), Volcanic degassing: Process and impact,481 Treatise on Geochemistry: Second Edition, pp. 111–179.482 Otto, F. (2001), The geometry of dissipative evolution equations: the porous medium equation.483 Patzek, T. W., F. Male, and M. Marder (2013), Gas production in the barnett shale obeys a simple484 scaling theory, Proceedings of the National Academy of Sciences,110(49), 19,731–19,736.485 Peng, D.-Y., and D. B. Robinson (1976), A new two-constant equation of state, Industrial & Engi-486 neering Chemistry Fundamentals,15(1), 59–64.487 Saffman, P. G. (1992), Vortex dynamics, in Theoretical Approaches to Turbulence, pp. 263–277,488 Springer.489 Scandella, B. P., L. Pillsbury, T. Weber, C. Ruppel, H. F. Hemond, and R. Juanes (2016), Ephemerality490 of discrete methane vents in lake sediments, Geophysical Research Letters,43(9), 4374–4381.491 Taylor, G. I. (1950), The instability of liquid surfaces when accelerated in a direction perpendicular492 to their planes. i, Proceedings of the Royal Society of London. Series A. Mathematical and Physical493 Sciences,201(1065), 192–196.494 Tek, M. R. (1989), Underground Storage of Natural Gas: Theory and Practice, 171, Springer Science495 & Business Media.496 V´azquez, J. L. (2007), The Porous Medium Equation: Mathematical Theory, Oxford University Press.497 Villani, C. (2021), Topics in optimal transportation, vol. 58, American Mathematical Soc.498 Wu, Y.-S., K. Pruess, and P. Persoff (1998), Gas flow in porous media with Klinkenberg effects,499 Transport in Porous Media,32, 117–137.500 Wu, Y.-S., J. Li, D.-Y. Ding, C. Wang, and Y. Di (2014), A generalized framework model for the501 simulation of gas production in unconventional gas reservoirs, Spe Journal,19(05), 845–857.502 23 Zhao, J., J. Yao, M. Zhang, L. Zhang, Y. Yang, H. Sun, S. An, and A. Li (2016), Study of gas flow503 characteristics in tight porous media with a microscale lattice boltzmann model, Scientific reports,504 6(1), 32,393.505 A Derivation of the porous medium equation (PME) and an-506 alytical solution for the linear case507 Combining Equations 1-3 and assuming g= 0, constant values for the temperature T, permeability k508 and porosity ϕ, it is straightforward to obtain (e.g., Al-Hussainy et al., 1966):509 ϕ k ∂ ∂t p Z(p)=∇·p µ(p)Z(p)∇p,(A.1) which, assuming Z(p)≈Z0and µ(p)≈µ0, and noting that ∂tp≡∂tp2/2pand further that p∇p≡510 p2/2p(for p= 0), transforms into the PME. In non-dimensional form it reads:511 ∂p2 ∂t =D∇2p2,(A.2) with diffusion coefficient512 D(x, t) = p(x, t) Π.(A.3) For negligible spatial variability of p(x, t), and hence D(x, t), Equation A.2 admits an exact ana-513 lytical solution in terms of an eigenfunction expansion of the residual ucalc(x, t)≡p2 calc(x, t)−p2 ss(x),514 between a given transient state p2 calc(x, t) and the steady-state p2 ss(x) (e.g., Carslaw and Jaeger, 1959;515 Crank, 1975):516 ucalc(x, t) = X n Anϕn(x) exp −λnt,(A.4) where the projections Anof u(x, t) onto the Laplacian eigenfunctions ϕn(x) are given as:517 An=Z1 0 u(x, t)ϕn(x)dx Z1 0 ϕn(x)2dx .(A.5) with518 ϕn(x, t) = sin pλnx(A.6) and λn= (nπ/L)2for imposed pressure at the inlet, and519 ϕn(x, t) = cos pλnx(A.7) 24 and λn= [(n+ 1/2)π/L]2, for imposed volumetric or mass flux at the inlet (n= 1,2, ...). The variable520 tis a diffusivity-scaled time to account for time-varying diffusivity (e.g., Crank, 1975):521 t=Zτ 0 D(s)ds. (A.8) We now define the calculated diffusion time-scale τcalc Din terms of the slowest-decaying modes of522 the eigenfunction expansion of Equation A.4, given for the eigenvalues λ1=π2and λ0=π2/42, for523 the inlet Dirichlet and Neumann conditions, respectively. For this, let the residual energy be defined524 by (e.g., Evans, 2022)525 ecalc(t) = 1 2Z1 0 u2 calc(x, t)dx (A.9) From Equation A.4 we have, at late times, that ecalc(t)/ecalc(0) ≈exp −λ1t. Thus, defining τcalc D 526 such that ecalc(t)/ecalc(0) = ϵ,τcalc Dsatisfies:527 t(τcalc D) = ln (1/ϵ) λ1 ,(A.10) or λ0for the Neumann condition.528 For non-negligible spatial variations of p(x, t) and D(x, t), D(t) in Equation 11 is replaced by an529 equivalent value at each time, Deff (t), given by the harmonic mean over x, that is:530 Deq(t) = 1 Rdx D(x,t) .(A.11) B Numerical simulations531 We simulate hydrogen injection into hydrogen-saturated 1-D reservoirs by solving numerically the532 coupled system of Equations 1, 2 and 3 using the Matlab Reservoir Simulation Toolbox (MRST) (Lie,533 2019; Lie and Møyner, 2021). The reservoir has a length and height of 600 m and 1 m, respectively.534 The top and bottom boundaries are assumed impervious. We vary the boundary condition at the535 right side while we keep a constant pressure of p0= 10 bar, equal to the initial pressure distribution.536 The viscosity µ, temperature Tand compressibility factor Zare assumed constant throughout, with537 values µ= 8.1×10−6Pa s−1,T= 320 K and Z= 1. The porosity is ϕ= 1 throughout. The domain538 is discretized into 2400 ×1 cells, giving a cell size of 0.025 m.539 25