The effect of gas compressibility during flow through a porous medium modifies the governing equation for pressure, and for low pressures, this leads to the so-called Klinkenberg effect. Darcy’s law remains unchanged, and accounting for the equation of state leads to a parabolic governing equation for pressure. The increased flow rate at low pressure modifies the nonlinearity of this governing equation. We present two broad approaches for constructing exact solutions to the steady state problem. First, geometric reduction is employed to construct exact solutions in Cartesian and polar coordinates. Next, the governing equation is interpreted so that a stream function exists for the flow, and this is used to demonstrate that solutions can only be found when the flow is irrotational. A second, rather broad, class of exact solutions is thus constructed from potential flows and their generalization for variable permeability cases. The latter leads to a non-constant coefficient problem, and we provide both an algorithm illustrating how to use an existing linear numerical solver to solve the nonlinear problem and an explicit exact solution for an annular domain.

Modeling the flow of gas through porous media has long been an area of interest for both pure and applied research. There are numerous natural and industrial examples of studies of gas flow through porous media. A recent series of papers highlights the natural behavior of porous media convection in the high Rayleigh number limit with applications toward sequestration of CO2 from deep saline aquifers.1,2 In the context of the mining industry, researchers have modeled the advection of oxygen through waste rock piles for the examination of sulfuric acid runoff3–5 as well as flow in tight shale-gas reservoirs.6–8 

Regardless of the application, a salient feature of the flow through a porous medium is that the hydraulic conductance is a function of the permeability of the medium. In many naturally occurring circumstances, the permeability is an intrinsic property of the medium, often being a function of geometry of the complex pore network9 (see the work by Chastanet, Royer, and Auriault10 for examples of the permeability for several different media). However, in the context of low pressure gas flow, experimental evidence suggests that the permeability is also a function of the nature of the flow itself. As early as 1931, Muskat and Botset11 performed experimental studies of gas flow through columns of glass beads and noted that the permeability of the medium was higher than expected. Since then, it has been noted that when pore sizes are only 1–1000 times larger than the mean free path of gas molecules, contact with pore walls is more likely than contact with other molecules.12,13 Thus, gas molecules slip along the pore boundaries, and this slippage effect adds an additional flux, thereby leading to a larger estimated permeability than liquids.

In the literature, this slippage effect is included in models by adding a term with permeability that depends on the inverse gas pressure. This is commonly known as the Klinkenberg correction.10,14 Under the Klinkenberg correction, as the local gas pressure decreases, the permeability coefficient increases. The slip-flow regime described by the Klinkenberg effect has practical importance within (among other contexts) the petroleum industry.10 For example, gas injection methods are often used as a means to extract hydrocarbons from low permeability shale gas reservoirs. Corrections to the permeability are required as these reservoirs have pore sizes on the scale of the mean free path of hydrocarbon molecules or could be under low reservoir pressure conditions.15 

While the effects of compressibility make the governing equation for pressure parabolic, the inclusion of the Klinkenberg effect modifies the equation’s nonlinearity. This implies that even in steady state, standard methods in classical linear partial differential equations cannot be applied verbatim. Thus, in order to obtain analytical solutions, either some simplifications must be made or classical methods must be modified. Analytical solutions to simplified versions of this equation have been proposed for the cases of steady-state and transient gas flow.16,17 When steady-state flow is assumed, exact solutions for one-dimensional gas flow have been obtained for constant, equivalent permeability (which is an approximation).16 In addition, exponential-like solutions have been proposed for transient pulse-decay methods to estimate the Klinkenberg coefficient.18 

While multiple methods have been presented in the literature, the approach has largely been piecemeal. The aim of this manuscript is thus twofold. First, we describe several classes of exact solutions to the Darcy–Klinkenberg equations when geometric simplification is possible. We begin by deriving the analytical pressure distribution for a variable permeability one-dimensional porous strip, followed by a two-dimensional porous annulus, incorporating Klinkenberg effects in both cases. Solutions for a porous strip and axisymmetric flow in annular domains have been obtained in the literature before16,17,19 but have assumed constant permeability; an unnecessary assumption as shown below. Second, we show that the Klinkenberg effect does not induce vorticity in the velocity field, but in some cases, a spatially variable permeability can. From there, we make the link with classical potential flow theory. This link allows us to represent the solution of the steady Darcy–Klinkenberg equations, with and without a spatially variable permeability, in closed form in terms of the solution of a standard elliptic problem. When the permeability is constant, potential flow solutions from classical fluid mechanics may be used to explicitly construct solutions, and we illustrate with several examples. When the permeability varies in space, the elliptic problem may be solved numerically, although we develop one example of an exact solution with variable permeability. We close the article with some brief conclusions and future directions.

In the porous flow regime where the continuum approximation is valid, Darcy’s law, in various guises, has been used as the governing equation for over a century. Following the study by Wu et al.,16 Darcy’s law reads as20,21

(1)

where u is the Darcy velocity of the gas phase, P is the gas pressure, kg is the permeability (here the subscript “g” has been added to suggest we are considering gas flow), and η is the dynamic viscosity of the fluid (gas in the present case). As noted in Introduction, empirical studies of porous gas flow11,14 have found discrepancies between the flow of liquids and gases. As the pore size decreases, the gas flow enters a transitional regime where boundary slip effects become significant. In these instances, the gas permeability is corrected via a term that is inversely proportional to pressure;14 the so called Klinkenberg correction is given by

(2)

where b is referred to as the Klinkenberg factor and k is the gas-phase permeability when Klinkenberg effects are negligible.

By applying the correction to Darcy’s law and using conservation of mass, it has been shown16 that the resulting governing equation for gas flow is a non-linear partial differential equation of parabolic type for pressure P(x,t), namely,

(3)

where φ is the porosity. It is worth noting that even if the Klinkenberg correction is neglected (b = 0), this equation is nonlinear. A non-zero b modifies the nonlinearity of the equation and adds complications to the task of seeking exact solutions. Throughout this article, we will consider the steady state version of the above, a nonlinear partial differential equation of elliptic type.

The first set of exact solutions to Eq. (3) that will be presented considers geometrically simplified situations. Figure 1 displays the two configurations of interest to this study (the case of a pressurized hemisphere would follow analogously in spherical polar coordinates): the case of a drop from constant atmospheric pressure across a rectangular slab of a porous medium and the corresponding problem for an annular region of a porous medium. Variable permeability, k=k(x), can also be incorporated into the exact solution procedure.

FIG. 1.

Diagram of the setup for simplified problem geometries: (a) one-dimensional porous strip and (b) two-dimensional porous annulus.

FIG. 1.

Diagram of the setup for simplified problem geometries: (a) one-dimensional porous strip and (b) two-dimensional porous annulus.

Close modal

The central idea is to use existing expressions for the divergence and gradient in Cartesian and polar coordinates and integrate them step-wise. The procedure is demonstrated in detail for one dimensional Cartesian coordinates on the interval [0, L], and the results are presented for both Cartesian and polar coordinates. Computer algebra packages such as Maple22 can prove to be very useful in evaluating the various integrals that appear as a part of the solution procedure.

The one-dimensional form of Eq. (3) is

(4)

Integrating once, we find

The left hand side can be treated via the chain rule, while the right hand side can be integrated formally to yield the solution in the quadrature,

(5)

Note that when the boundary conditions at x = 0 are of Dirichlet type (i.e., P is specified at x = 0), c2 can be immediately determined, and the boundary condition at x = L can be used to compute c1.

At this point, the formal solution is a quadratic equation for pressure. To determine a unique solution, we apply the quadratic formula and choose the root that guarantees that pressure is non-negative. This yields the solution

(6)

While trivial from an algebra point of view, this step is intimately linked with the fact that the governing equation is nonlinear. As such, we cannot simply multiply the solution by a scaling factor and expect to retain the result as a solution, as we would for a linear governing equation.

Several sample solutions are given in Fig. 2. The parameters used to generate the plots are given in Table I. Panel (a) shows the effect of exponentially varying permeability, k(x) = exp(λx), compared to the reference, constant coefficient solution (black curve). The red (blue) curve shows the solution for an exponential increase (decrease) in permeability with x. It can be seen that increasing (decreasing) permeability yields a concave up (concave down) solution. In panel (b), we show the results for a piecewise constant permeability,

where H(·) is the Heaviside step function. It can be seen that a 90% decrease in permeability (blue curve) has a far more drastic effect on the pressure distribution than a ten-fold increase in permeability (red curve).

FIG. 2.

Pressure profiles in Cartesian coordinates scaled by the pressure at the left boundary. (a) Exponential profiles with the exponential parameter λ. (b) Piecewise constant profiles with the magnitude of jump set by β.

FIG. 2.

Pressure profiles in Cartesian coordinates scaled by the pressure at the left boundary. (a) Exponential profiles with the exponential parameter λ. (b) Piecewise constant profiles with the magnitude of jump set by β.

Close modal
TABLE I.

Table of example values used throughout the article. b represents the Klinkenberg factor, and P1 represents the atmospheric pressure (P2 is not explicitly included as all results are in terms of a scaled pressure difference), a higher boundary pressure. L is the width of the porous strip, and Rin and Rout are the inner and outer radii of the porous annulus, respectively. These values are varied, and both values for each are included in the table. λ represents the rate with which the permeability changes with space, and β represents the magnitude of the instantaneous change in the permeability.

ParameterValueDimensions
b 4.75 × 104 Pa 
P1 1.013 25 × 105 Pa 
L 
Rin 0.1, 1 
Rout 1.1, 2 
λ 2, −2 m−1 
β 9, −0.9 ⋯ 
ParameterValueDimensions
b 4.75 × 104 Pa 
P1 1.013 25 × 105 Pa 
L 
Rin 0.1, 1 
Rout 1.1, 2 
λ 2, −2 m−1 
β 9, −0.9 ⋯ 

The corresponding problem for polar coordinates, in which the pressure is assumed to be steady and varying only in the radial direction, reads

(7)

A general solution can be acquired analogously to the Cartesian case and reads

(8)

Several sample solutions for the polar case are given in Fig. 3. The parameters used to generate the plots are given in Table I. The pressure drop is taken to be 10% of the reference atmospheric pressure, except for the test of nonlinearity in panel (b). Panel (a) shows the effect of exponentially varying permeability, k(r) = exp(λr), compared to the reference, constant coefficient solution (black curve). As for the Cartesian case, the red (blue) curve shows the solution for an exponential increase (decrease) in permeability with r. It can be seen that increasing (decreasing) permeability yields a concave up (concave down) solution. In panel (b), we test the effective nonlinearity of the solution. We do this by reducing the pressure decrease by 50% (green curve) and show this solution along with both the original solution (black curve) and a scaled version of the green curve (pink dashed curve). For this test, the pressure perturbation was increased to 50% of the reference pressure. While some differences are visible, they are small.

FIG. 3.

Pressure profiles for the annular region in polar coordinates scaled by the pressure at the outer boundary. (a) Pressure profiles with power law permeability, with parameter λ (Rin, Rout) = (0.1, 1.1), (b) pressure profiles for the nonlinearity test constant permeability, (c) pressure profiles with power law permeability, with parameter λ (Rin, Rout) = (1, 2) showing the effect of larger interior radius, and (d) pressure profiles with piecewise constant permeability profiles with the magnitude of jump set by β (Rin, Rout) = (0.1, 1.1).

FIG. 3.

Pressure profiles for the annular region in polar coordinates scaled by the pressure at the outer boundary. (a) Pressure profiles with power law permeability, with parameter λ (Rin, Rout) = (0.1, 1.1), (b) pressure profiles for the nonlinearity test constant permeability, (c) pressure profiles with power law permeability, with parameter λ (Rin, Rout) = (1, 2) showing the effect of larger interior radius, and (d) pressure profiles with piecewise constant permeability profiles with the magnitude of jump set by β (Rin, Rout) = (0.1, 1.1).

Close modal

Panel (c) of Fig. 3 revisits the exponential permeability profile, but for an annular region with a larger interior radius (Rin = 1). It can be seen that the decrease in curvature leads to a far smaller effect for the same change in permeability. Finally, in panel (d), we reconsider the piecewise constant case

where H(·) is the Heaviside step function. Once again, it can be seen that a 90% decrease in permeability (blue curve) has a far more drastic effect on the pressure distribution than a ten-fold increase in permeability (red curve). The results can be compared to the work presented in the study by Penney and Stastna,23 who demonstrated the effects on the fluid pressure distribution of the bands of low and high permeability (blockage and ducting, respectively). Their work assumed that the saturating fluid was incompressible or, in other words, either fluid or air was incompressible at near atmospheric pressures. Hence, they did not include Klinkenberg effects. They found that modifications of the surface flux into the domain were dominated by the magnitude of the permeability change instead of the width and location of the change. Blockages (a band of lower permeability) caused a significant decrease in the inward fluid flux, while ducting (a band of higher permeability) allowed a significant annular flux within the ducting region, but with very minimal differences in the non-ducting outer region, when compared to the flux of a constant permeability case.

The basic Eq. (3) in steady state and with a constant permeability may be rewritten as

where

In two dimensions, this expression guarantees the existence of a streamfunction so that in Cartesian coordinates, we have

(9)
(10)

The nonlinearity of the above-mentioned equations along with the requirement that pressure is non-negative makes the problem of determining a valid pressure distribution non-trivial.

First, consider a simple case, a uniform stream, to illustrate the solution procedure for P. For the stream-function Ψ1 = x, the governing equations read

(11)
(12)

Integrating Eq. (11) with respect to x and Eq. (12) with respect to y results in

(13)
(14)

In order to be able to find a unique solution, we must choose c1(y) and c2(x) so that the two equations are the same. Here, the choice c1(y) = −y and c2(x) = 0 clearly yields a single equation. Solving for P(x, y) gives

(15)

While it is possible to carry out the procedure for a broad class of streamfunctions, the procedure does not always work. A cautionary example of this is given by the class of closed streamlines, including circular streamlines. Consider the family of streamfunctions

(16)

The governing equations read

(17)
(18)

Integrating both equations, we find the candidate solutions,

(19)
(20)

However, there is no possible choice of c1(y) and c2(x) that reduces these two equations to a single equation, and hence, a unique solution equation for p is impossible to find.

While moderately useful for constructing solutions, the integration procedure is more useful as a means to restrict classes of streamfunctions. For example, consider streamfunctions of the form

(21)

for differentiable functions f(·) and g(·) whose anti-derivatives we will denoted by capital letters [F(·) and G(·)]. Substitution and integration give

(22)
(23)

We immediately note that we must have c1(x) = c2(y) = 0, and by equating the remaining portions of the right hand sides, we find that either

or

Taking the derivative, we find the standard second order, constant coefficient ordinary differential equation

implying either f is sinusoidal and g is exponential or vice versa.

A restriction with a broader set of consequences can be obtained by forming the equation for the vorticity, ∇2ψ. From the governing equations, we find that

(24)

Thus, neither the nonlinearity due compressible effects nor the Klinkenberg correction can lead to the possibility of rotational flow.

We have demonstrated in Sec. III that any flow governed by the Darcy–Klinkenberg Eq. (3) is irrotational. This is in addition to the condition that our vector V is divergence-free by construction. It is important to note that V is not the physical velocity, which is compressible. The two are related via the relation

(25)

which we will employ below.

From the divergence-free and irrotational restrictions on V, we are thus guaranteed the existence of a potential ϕ. Recall that, for now, we have assumed that permeability is constant.

The governing equations read

(26)
(27)

Integration is now relatively trivial, with the choice c1(y) = c2(x) = 0 guaranteeing a solution. Choosing the positive root gives the expression

(28)

Since the potential is only defined up to an additive constant, we can choose the arbitrary constant to guarantee a physical (i.e., non-negative) pressure distribution.

The solution may be expanded in b so that the case of no Klinkenberg effects is the leading order term,

(29)

This is a useful expression since it shows that at first order, Klinkenberg effects merely shift the value of pressure by a constant, with geometric effects evident only at second order.

Once the pressure is known, an expression for the velocity field is given by

(30)

where subscripts on ϕ denote partial derivatives. Since many ϕ are available in textbooks (that by Kundu and Cohen,24 for example), e.g., as the real part of any analytic function F(z) of a complex variable, we are able to explicitly construct broad classes of solutions for pressure and the velocity field.

Here, we consider two examples from classical potential flow theory, namely, the corner flow,24 with the complex potential

where z is a complex number, M is a scaling factor, and n is an integer. When n = 2, this potential yields the flow in the first quadrant bounded by walls at the axes. When n = 3, this potential yields the flow in the 60° wedge in the first quadrant. The real part of ϕ(c) is the potential used to construct the pressure in the Klinkenberg theory. The imaginary part gives the streamlines.

The sample solutions are shown in Fig. 4 with the classical solutions on the left and the new Klinkenberg theory solutions on the right. It can be seen that the flow field remains the same but the pressure distribution is significantly modified.

FIG. 4.

Exact solutions constructed from potential flow solutions for flow in a corner. (a) The classical solution for a first quadrant flow ϕ = Mz2 + c with the potential shaded and the streamlines overlaid in white, (b) the pressure (shaded) constructed from the potential in (a) and the velocity field superimposed as a quiver plot, (c) the classical solution for a 60° flow ϕ = Mz3 + c with the potential shaded and the streamlines overlaid in white, and (d) the pressure (shaded) constructed from the potential in (b) and the velocity field superimposed as a quiver plot. The dividing streamline from (c) has been added for visual clarity. Note that c is the additive constant required to ensure real-valued pressure.

FIG. 4.

Exact solutions constructed from potential flow solutions for flow in a corner. (a) The classical solution for a first quadrant flow ϕ = Mz2 + c with the potential shaded and the streamlines overlaid in white, (b) the pressure (shaded) constructed from the potential in (a) and the velocity field superimposed as a quiver plot, (c) the classical solution for a 60° flow ϕ = Mz3 + c with the potential shaded and the streamlines overlaid in white, and (d) the pressure (shaded) constructed from the potential in (b) and the velocity field superimposed as a quiver plot. The dividing streamline from (c) has been added for visual clarity. Note that c is the additive constant required to ensure real-valued pressure.

Close modal

While both the explicit solutions of Sec. II, and those derived from potential solutions in Sec. IV, have their uses, an interesting point to explore is the following question: “To what extent is the assumption of a constant permeability required?” If we group the ratio of permeability and viscosity as the general function a(x, y), the problem can be written as

(31)

or alternatively

where

A quick calculation shows this can be rewritten as

We can thus simplify the expression for vorticity as

(32)

because the curl of a gradient is identically zero. In the absence of Klinkenberg effects, this condition is analogous to the baroclinic vorticity generation term in density stratified flow.24 When Klinkenberg effects are present, the condition is geometrically non-trivial. Thus, while the Klinkenberg effect in itself is not sufficient to induce vorticity, in some circumstances, a variable permeability is.

Regardless of whether vorticity is present or not, Eq. (31) can be rewritten as

(33)

where ξ = P2/2 + bP. While Eq. (33) does not have a massive number of exact solutions that Laplace’s equation does, it does have the advantage of being in standard form. This means almost every standard finite element package, for example, will have well validated solvers for it. The algorithm for solution would proceed as follows:

As an aside, polar coordinates are ideal for demonstrating how to choose the boundary conditions on ξ that are consistent with a physically valid set of boundary conditions on P. Say P at the outer boundary is written as a truncated cosine series

Then the quadratic equation relating ξ and P gives

The right hand side may be simplified using cosine addition formulas. Thus, in polar coordinates, the procedure for solving P in terms of ξ is explicit, effectively a recipe. As a simple example, consider P(Rout, θ) = P1 + a1 cos(θ). After some simple algebra, we find that

Thus, in order to impose a boundary condition with a constant value perturbed by a single azimuthal, mode-1 sinusoid on P, we must impose a boundary condition on ξ that has non-zero values for azimuthal mode-0, mode-1, and mode-2.

As noted above, exact solutions for general a(x, y) are rare, but in some examples, they can be constructed. Consider an annular region, as shown in Fig. 1. In polar coordinates, and assuming a single mode solution in the azimuthal,

(34)

some simple algebra shows that the equation for fn reads

(35)

Choosing a(r) = 1/r allows for a significant simplification to

(36)

This is an Euler-type differential equation, and if we define

(37)

then the solution can be written as

(38)

where subscripts on γ denote the positive or negative root.

We have used this solution procedure and the algorithm mentioned above to construct an exact solution for P using the parameters referenced in Table I. The solution is shown in Fig. 5. Panel (a) shows the shaded pressure field, scaled by the constant reference pressure P1. Panel (b) shows the solution of the linear problem, ξ, that is used to construct P in panel (a). Panels (c) and (d) show azimuthal and radial profiles of pressure, respectively. The imposed change in pressure is 20% of the reference value. The green dashed curve in panel (c) shows the single cosine fit of the solution. The nonlinear effect is appreciable, although not dominant.

FIG. 5.

(a) Shaded contours of P for the exact, variable k(r) solution scaled by the reference value of pressure P1. (b) Shaded contours of ξ (scaled by b2), the solution of the linear problem that is used to construct the pressure in panel (a). (c) Azimuthal profile of the pressure (scaled by P1) on the interior (blue) and exterior (black) boundaries; a single cosine fit to the solution on the exterior boundary is shown as a dashed green line. (d) Radial profile of pressure (scaled by P1) at a fixed value of θ = 0.

FIG. 5.

(a) Shaded contours of P for the exact, variable k(r) solution scaled by the reference value of pressure P1. (b) Shaded contours of ξ (scaled by b2), the solution of the linear problem that is used to construct the pressure in panel (a). (c) Azimuthal profile of the pressure (scaled by P1) on the interior (blue) and exterior (black) boundaries; a single cosine fit to the solution on the exterior boundary is shown as a dashed green line. (d) Radial profile of pressure (scaled by P1) at a fixed value of θ = 0.

Close modal

It is worth noting also that the standard form Eq. (33) allows theorems from the theory of elliptic partial differential equations, such as the maximum principle,25 to be applied. These results suggest that pressure maxima must occur on the boundary of a closed domain under only mild technical conditions, which are satisfied in any laboratory or field setting.

We have presented several classes of solutions for the problem of gas flow through a porous medium, including the Klinkenberg correction. The geometrically simple cases could be used in laboratory settings for material characterization or as a springboard to study the inverse problem of whether pressure measurements can be used to infer spatial distributions of permeability. The solutions constructed from potential flow provide a surprisingly broad class of solutions to the nonlinear Darcy–Klinkenberg problem and unify the many solutions presented in the literature on an ad hoc basis.

It is well known that solutions of Laplace’s equation with purely Neumann (flux) boundary conditions are only specified up to an additive constant. For gas saturated porous media, the governing steady state equation is nonlinear, and hence, an arbitrary constant cannot be added. The boundary value problem outlined in Algorithm 1 demonstrates how a constant added to the linear problem for ξ is expressed in the solution for P. However, the potential used to construct the solution for P in (28) covers an unbounded region (i.e., the “corner” of the first quadrant that is of interest), and hence, it is only defined up to an additive constant. Practical use of the potential thus has the additional problem of specifying the value of this constant a priori. We can estimate a bound on the constant that ensures that pressures (and therefore fluid velocities) are real-valued. This is done simply by showing the values of the additive constant that ensure the argument under the radical in (28) is greater than zero. However, in the absence of sufficient boundary conditions, the nature of Klinkenberg correction precludes us from bounding the constant in such a way that we can find an upper bound on the pressure field itself. This is because, given a particular porous medium, as the pressure approaches zero, the Klinkenberg model, (2), states that the permeability must continue to grow. Several studies (for example, that by Ziarani and Aguilera26 or Tang, Tao, and He27) have demonstrated that higher order corrections are required as the Knudsen number (i.e., the ratio of the mean free path of molecules to the pore width) becomes large. For example, see the study by Ziarani and Aguilera26 for details on higher order terms in the Klinkenberg correction. There is potential for future work in developing simple mathematical models in the high Knudsen number limit. However, as the Knudsen number becomes very large (greater than about ten), the continuum assumption breaks down, and molecules act ballistically. See the work by Civan28 for an example of the behavior of the hydraulic conductance as a function of the Knudsen number.

ALGORITHM 1.

Algorithm for constructing a solution of the nonlinear boundary value problem due to Klinkenberg effects, from a linear elliptic problem in standard form.

1. Use P to find boundary conditions on ξ 
2. Solve (33) 
3. Construct P from ξ via the quadratic formula 
4. Construct u from P
1. Use P to find boundary conditions on ξ 
2. Solve (33) 
3. Construct P from ξ via the quadratic formula 
4. Construct u from P

When the permeability is variable, the systematic algorithm (Algorithm 1) we have provided can be used to build a solver for the nonlinear problem from any well-validated solver for linear elliptic problems. The exact solution we have presented for radially dependent permeability of the power law type is to the best extent of our knowledge unique in the literature, encompassing all nonlinearities—a nontrivial azimuthal structure and variable permeability.

In order to establish exact solutions based on potential flow, we were able to demonstrate that compressibility, including the Klinkenberg correction, is not sufficient to induce vorticity in the flow. However, variable permeability is able to induce vorticity, as shown by Eq. (32).

Future work could follow several avenues in addition to the high Knudsen number limit discussed above. A general numerical solver that implements the above-mentioned methodology could be written for community use, the geometrically simplified cases could be used to explore the inverse problem of determining permeability distribution from pressure measurements, and the physical assumptions of the model could be loosened to consider variable temperature/buoyancy effects.

This research project was funded in part by the Ontario Ministry of Research, Innovation and Science through an Ontario Research Fund-Research Excellence grant (Grant No. RE09-061). The authors would also like to thank Jason Olsthoorn whose suggestions and criticisms have led to substantial improvements to this paper.

The authors have no conflicts to disclose.

Marek Stastna: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). Andrew P. Grace: Formal analysis (equal); Supervision (supporting); Writing – original draft (supporting); Writing – review & editing (equal). Travis Robinson: Formal analysis (supporting); Writing – original draft (supporting).

All codes are available through the first author’s GitHub page.

1.
S.
Pirozzoli
,
M.
De Paoli
,
F.
Zonta
, and
A.
Soldati
, “
Towards the ultimate regime in Rayleigh–Darcy convection
,”
J. Fluid Mech.
911
,
R4
(
2021
).
2.
M.
De Paoli
,
S.
Pirozzoli
,
F.
Zonta
, and
A.
Soldati
, “
Strong Rayleigh–Darcy convection regime in three-dimensional porous media
,”
J. Fluid Mech.
943
,
A51
(
2022
).
3.
X.
Chi
,
R. T.
Amos
,
M.
Stastna
,
D. W.
Blowes
,
D. C.
Sego
, and
L.
Smith
, “
The Diavik Waste Rock Project: Implications of wind-induced gas transport
,”
Appl. Geochem.
36
,
246
255
(
2013
).
4.
R.
Lefebvre
,
D.
Hockley
,
J.
Smolensky
, and
P.
Gélinas
, “
Multiphase transfer processes in waste rock piles producing acid mine drainage: 1: Conceptual model and system characterization
,”
J. Contam. Hydrol.
52
,
137
164
(
2001
).
5.
R. T.
Amos
,
D. W.
Blowes
,
B. L.
Bailey
,
D. C.
Sego
,
L.
Smith
, and
A. I. M.
Ritchie
, “
Waste-rock hydrogeology and geochemistry
,”
Appl. Geochem.
57
,
140
156
(
2015
).
6.
Y.-S.
Wu
,
J.
Li
,
D.-Y.
Ding
,
C.
Wang
, and
Y.
Di
, “
A generalized framework model for the simulation of gas production in unconventional gas reservoirs
,”
SPE J.
19
,
845
857
(
2014
).
7.
N.
Alharthy
,
M.
Al Kobaisi
,
M. A.
Torcuk
,
H.
Kazemi
, and
R.
Graves
, “
Physics and modeling of gas flow in shale reservoirs
,” in
Abu Dhabi International Petroleum Conference and Exhibition
(
OnePetro
,
2012
).
8.
Y. J.
Zhuang
,
H. Z.
Yu
, and
Q. Y.
Zhu
, “
Experimental and numerical investigations on the flow around and through the fractal soft rocks with water vapor absorption
,”
Int. J. Heat Mass Transfer
96
,
413
429
(
2016
).
9.
M.
Honarpour
and
S. M.
Mahmood
, “
Relative-permeability measurements: An overview
,”
J. Pet. Technol.
40
,
963
966
(
1988
).
10.
J.
Chastanet
,
P.
Royer
, and
J.-L.
Auriault
, “
Does Klinkenberg’s law survive upscaling?
,”
Transp. Porous Media
56
,
171
198
(
2004
).
11.
M.
Muskat
and
H. G.
Botset
, “
Flow of gas through porous materials
,”
Physics
1
,
27
47
(
1931
).
12.
H.
Adzumi
, “
Studies on the flow of gaseous mixtures through capillaries. II. The molecular flow of gaseous mixtures
,”
Bull. Chem. Soc. Jpn.
12
,
285
291
(
1937
).
13.
F.
Javadpour
,
D.
Fisher
, and
M.
Unsworth
, “
Nanoscale gas flow in shale gas sediments
,”
J. Can. Pet. Technol.
46
,
55
61
(
2007
).
14.
L.
Klinkenberg
, “
The permeability of porous media to liquids and gases
,”
Am. Petrol. Inst., Drill. Proj. Pract.
2
,
200
213
(
1941
).
15.
F. P.
Wang
,
R. M.
Reed
,
A.
John
, and
G.
Katherine
, “
Pore networks and fluid flow in gas shales
,” in
SPE Annual Technical Conference and Exhibition
(
OnePetro
,
2009
).
16.
Y. S.
Wu
,
K.
Pruess
, and
P.
Persoff
, “
Gas flow in porous media with Klinkenberg effects
,”
Transp. Porous Media
32
,
117
137
(
1998
).
17.
A. L.
Baehr
and
C. J.
Joss
, “
An updated model of induced airflow in the unsaturated zone
,”
Water Resour. Res.
31
,
417
421
, (
1995
).
18.
Y.
Liang
,
J. D.
Price
,
D. A.
Wark
, and
E. B.
Watson
, “
Nonlinear pressure diffusion in a porous medium: Approximate solutions with applications to permeability measurements using transient pulse decay method
,”
J. Geophys. Res.: Solid Earth
106
,
529
535
, (
2001
).
19.
A. L.
Baehr
and
M. F.
Hult
, “
Evaluation of unsaturated zone air permeability through pneumatic tests
,”
Water Resour. Res.
27
,
2605
2617
, (
1991
).
20.
H.
Darcy
,
Les Fontaines Publiques de la Ville de Dijon: Exposition et Application des Principes à Suivre et des Formules à Employer Dans les Questions de Distribution D’eau: Ouvrage Terminé Par un Appendice Relatif Aux Fournitures D’eau de Plusieurs Villes, au Filtrage des Eaux et à la Fabrication des Tuyaux de Fonte de Plomb, de Tôle et de Bitume
(
V. Dalmont
,
1856
), Vol. 2.
21.
S.
Whitaker
, “
Flow in porous media I: A theoretical derivation of Darcy’s law
,”
Transp. Porous Media
1
,
3
25
(
1986
).
22.
Maplesoft, a division of Waterloo Maple, Inc., Maple.
23.
J.
Penney
and
M.
Stastna
, “
Numerical simulations of flow through a variable permeability circular cylinder
,”
Phys. Fluids
33
,
117113
(
2021
).
24.
P. K.
Kundu
and
I.
Cohen
,
Fluid Mechanics
, 4th ed. (
Academic Press
,
Amsterdam
,
2008
).
25.
M. H.
Protter
and
H. F.
Weinberger
,
Maximum Principles in Differential Equations
(
Springer Science & Business Media
,
2012
).
26.
A. S.
Ziarani
and
R.
Aguilera
, “
Knudsen’s permeability correction for tight porous media
,”
Transp. Porous Media
91
,
239
260
(
2012
).
27.
G. H.
Tang
,
W. Q.
Tao
, and
Y. L.
He
, “
Gas slippage effect on microscale porous flow using the lattice Boltzmann method
,”
Phys. Rev. E
72
,
056301
(
2005
).
28.
F.
Civan
, “
Effective correlation of apparent gas permeability in tight porous media
,”
Transp. Porous Media
82
,
375
384
(
2010
).