Exact solutions for flow through porous media with the Klinkenberg effect

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 CO₂ 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 runoff^3–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 network⁹ (see the work by Chastanet, Royer, and Auriault¹⁰ 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 Botset¹¹ 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.¹⁰ 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.¹⁵ 

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).¹⁶ In addition, exponential-like solutions have been proposed for transient pulse-decay methods to estimate the Klinkenberg coefficient.¹⁸ 

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 before^16,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.,¹⁶ Darcy’s law reads as^20,21

where $u⃗$ is the Darcy velocity of the gas phase, P is the gas pressure, k_g 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 flow^11,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;¹⁴ the so called Klinkenberg correction is given by

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 shown¹⁶ that the resulting governing equation for gas flow is a non-linear partial differential equation of parabolic type for pressure $P(x⃗,t)$⁠, namely,

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.

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 Maple²² 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

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,

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

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

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).

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

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

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.

Panel (c) of Fig. 3 revisits the exponential permeability profile, but for an annular region with a larger interior radius (R_in = 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,²³ 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

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 Ψ₁ = x, the governing equations read

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

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

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

The governing equations read

Integrating both equations, we find the candidate solutions,

However, there is no possible choice of c₁(y) and c₂(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

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

We immediately note that we must have c₁(x) = c₂(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, ∇²ψ. From the governing equations, we find that

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

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

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

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,

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

where subscripts on ϕ denote partial derivatives. Since many ϕ are available in textbooks (that by Kundu and Cohen,²⁴ 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,²⁴ 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.

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

or alternatively

where

A quick calculation shows this can be rewritten as

We can thus simplify the expression for vorticity as

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.²⁴ 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

where ξ = P²/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(R_out, θ) = P₁ + a₁ 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,

some simple algebra shows that the equation for f_n reads

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

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

then the solution can be written as

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 P₁. 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.

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,²⁵ 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 Aguilera²⁶ or Tang, Tao, and He²⁷) 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 Aguilera²⁶ 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 Civan²⁸ for an example of the behavior of the hydraulic conductance as a function of the Knudsen number.

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.