Gaseous flow through heterogeneous, partially connected networks of pipes

Simulations of flow of an ideal gas through heterogeneous simple cubic pipe networks with different pipe radius distributions and variable bond coordination numbers were performed. Networks with monomodal and bimodal radius distributions were constructed. A very wide range of Knudsen numbers was achieved. Flow simulations of purely viscous gases and incompressible liquids were also carried out for comparison. The permeability to gas in the purely viscous regime was larger than the permeability to an incompressible liquid. Based on a variety of computational tests, this result was likely not a numerical artifact. The simulated macroscopic flow behavior differed from the underlying single pipe model, depending on the radius distribution, network connectivity and magnitude of the externally applied pressure gradient, and was compatible with the Klinkenberg analysis only when the maximum Knudsen number used in each simulation was lower than 1. In this condition, the Klinkenberg coefficient was nearly proportional to the inverse of the network hydraulic radius while the effect of the radius distribution was weak and that of the network connectivity essentially negligible. The bimodal simulations displayed a typical percolation behavior, with the Klinkenberg coefficient remaining constant as long as the large pipe population was connected.

Gas flow through pipes and ducts has been a problem of great interest to physicists for nearly 150 years 1 . It was early recognized that, at high gas pressures, gas flow obeys identical laws as liquid flow, although the compressibility of gas is much higher than that of liquids. But gas flow appeared utterly different at low gas pressures. Careful experiments showed that, with all other variables being held constant, the apparent transmissivity of long pipes increased with decreasing mean gas pressure. Moreover, the magnitude of this effect was observed to increase when thinner pipes were used.
Analysis of the experimental evidence demonstrated that gas flow depends on the dimensionless Knudsen number K n = λ/R, where R is the pipe radius and λ = (μ/p) (πR g T/2) 1/2 the mean free path of the gas molecules, with p denoting the gas pressure, μ the viscosity, R g the specific gas constant and T the temperature. Two extreme flow regimes can be described depending on K n . At very high gas pressures or, equivalently, very low K n , λ is infinitesimally small and the gas motion is controlled by the collisions of gas molecules among themselves. In this limit, the gas can be described as a viscous fluid subject to the usual equations of continuum fluid dynamics, supplemented with a no-slip boundary condition along the internal pipe wall. For long pipes, the mass flow rate is derived from the compressible Poiseuille flow formula: where L is the pipe length, Z the gas compressibility factor (Z = 1 for ideal gases), p the mean gas pressure and Δp the pressure difference driving the flow. At very low gas pressures (i.e., very large K n ), λ is many times greater than the pipe radius, intermolecular collisions are essentially inexistent and the gas motion is controlled by the collisions of gas molecules with the internal pipe wall. In these conditions, continuum fluid dynamics and the concept of viscosity become inapplicable. The proper description is instead provided by Boltzmann's equation, which can be shown to reduce to a diffusion equation where the independent variable is gas pressure 2 . This flow regime is usually called free molecular or Knudsen diffusion. In long pipes, the Knudsen diffusion mass flow rate is expressed as 3  One implication of equation 2 is that, for ideal gases, Knudsen diffusion is insensitive to the mean gas pressure. The limits of the viscous flow regime and Knudsen diffusion are usually estimated to be K n < 0.001 and K n > 10, respectively. The problem is to discover what happens in the transition range, 0.001 < K n < 10. It was early realized that, for K n slightly greater than 0.001, the concept of viscosity and the Navier-Stokes equation can be retained. The net effect of the increasing number of molecule-wall collisions is merely to change the gas velocity boundary condition along the internal pipe wall (see Klinkenberg 4 and references therein). The finite tangential gas velocity developing along the internal pipe wall can be expressed based on the kinetic theory of gases as: where σ, the tangential momentum accommodation coefficient (TMAC), is the fraction of wall-colliding molecules that undergo diffusive (as opposed to specular) reflections. The TMAC is usually assumed equal to 1, although lower values have been inferred from experimental data and numerical simulations [5][6][7][8][9] . By application of a perturbation method and equation 3 to long cylindrical pipes 10 , the mass flow rate is found to be equal to: S g 4 n where K n is the Knudsen number evaluated at the mean gas pressure (for simplicity, the mean Knudsen number will be denoted K n hereafter). Equation 4 describes the so-called slip flow regime and is normally deemed to hold only up to K n ≈ 0.1, although the Klinkenberg analysis, which is based on the slip flow concept, is often used for much higher values of K n . How to model the transitional gas flow regime, i.e., for 0.1 < K n < 10, is not totally settled at present. According to one school of thought,  m K is negligible in comparison to  m V when K n approaches zero while the reverse is true for very high values of K n . Thus, the sum +   m m V K approaches the correct limits at high and low Knudsen numbers and can be assumed to provide an estimate of the mass flow rate with relatively small errors for any value of K n in the transitional flow regime 11 . This notion is theoretically justified if viscous flow and molecular diffusion can be considered to be independent and uncoupled 12 . A number of variations of this model were used to study gas flow in porous media 3,[13][14][15][16][17] .
An alternative idea for modeling transitional flow is to note that the right-hand side of equation 3 can be viewed as the first term of a Taylor expansion in λ or, equivalently, K n . The idea is then to model transitional flow by adding higher-order terms 18 . A variety of second-order slip models have been proposed in recent years 1,8,9,19,20 . Among them, the Beskok and Karniadakis's model 7 (BK) is presently receiving the most attention [21][22][23][24][25] . According to the BK model, the mass flow rate is given by: In the BK original work, the ideal gas law and total momentum accommodation were assumed, implying Z = 1 and σ = 1 7 . The function F defined by equation 6 was then empirically determined by comparison with experimental and numerical simulation data 7 . The empirical values b BK = −1, α 0 = 6/5, α 1 = 4 and α 2 = 2/5 were inferred from the Loyalka and Hamoodi's computational data 26  when K n approaches zero, while it becomes asymptotically linear with a slope equal to 6, when K n grows to infinity. Although the value of α 0 mentioned above provides the best fit to the Loyalka and Hamoodi's data 26 , it is not consistent with the Knudsen diffusion limit. Civan 22 proposed to use α 0 = 64/15π instead (he also modified the first factor in the right-hand side of equation 6 to avoid the ArcTan function). This change of α 0 only causes minor changes to F in the low Knudsen number limit (i.e., the leading terms are unchanged) and produces a steeper slope (32/3) in the large Knudsen number limit. The BK input values (including Z = 1 and σ = 1) will be used hereafter. According to Klinkenberg 4 , the ratio of the gas permeability to the intrinsic permeability of porous media, k g /k, is a linear function of 1/p of the form, k g /k = 1 + b K /p, where the Klinkenberg coefficient b K is a positive constant. Since =   F m m / T V for a single pipe is analogous to k g /k, it is interesting to assess how different F is from the Klinkenberg linear model, which, in the case of F, can be rewritten as F = A + B K n , where A should be equal to unity. I constructed discrete representations of F spanning intervals of Knudsen  higher maximum values K max as shown in Fig. 1 (I used the log-log scale for a better visibility of the low K n region). Straight lines did fit these discrete sets of data-points very well in all cases (the goodness-of-fit R 2 coefficients were all equal to 0.9999 or better), but, as can be expected from the non-linearity of equation 6, A and B varied with K max . The intercept A was always lower than 1, increasing with decreasing K max and asymptotically reaching 1 only when K max dropped below 1. The slope B decreased with K max , presumably approaching 4, the value predicted by equation 4, for vanishing values of K max (Fig. 1). These observations suggest that standard linear regression analysis applied to limited discrete datasets cannot easily detect the slight upward curvature of F even for K n as high as 100. The implication is that the Klinkenberg linear extrapolation of experimental gas permeabilities to zero inverse gas pressure 1/p may lead to significant underestimation of the intrinsic permeability k in porous media with a behavior analogous to that expressed by equation 6. The obvious solution to this problem is, of course, to measure the intrinsic permeability directly by performing gas flow tests at very high gas pressures, i.e., in conditions where equation 1 applies. But this method cannot work in practice owing to the deformability of real porous media and the ensuing pressure sensitivity of intrinsic permeability.
One problem of major interest is to extend the single pipe models discussed above to general porous media. Among many possible approaches, one of the most popular is to follow the equivalent channel model (ECM) method 27,28 . According to ECM, the complex features of the pore space can be lumped into a few parameters characterizing a single representative pore channel (e.g., volume fraction, cross-sectional size and tortuosity factor). If the equivalent channel is idealized as a simple cylindrical pipe, it is easy to demonstrate that the intrinsic permeability is given by 27,28 : where φ denotes the porosity, τ the tortuosity (i.e., ratio of the pipe length to its projection in the nominal flow direction) and the characteristic pipe radius is defined as the hydraulic radius R H (i.e., twice the volume-to-surface ratio of the pore space). Note that φ and R H can be physically measured in most porous media, while τ is not directly accessible. Variations of the ECM approach have been frequently used to model gas flow, sometimes in combination with numerical simulations of gas flow through three-dimensional reconstructions of the pore space 1,14,20,[22][23][24]29 . The main drawback of the ECM approach is that it does not account for the very large variability in pore shape and size typical of most porous media. One simple and versatile alternative approach is to simulate gas flow through heterogeneous networks of pipes by incorporating the single pipe models discussed above in the Kirchoff equations 3,15-17 . Imperfectly connected networks can also be constructed, allowing investigation of the effect of pore connectivity on gas flow 17 .

Figure 1.
Three discrete representations of the function F spanning intervals of Knudsen numbers reaching a maximum value K max equal to 74 (black dots), 7.4 (blue dots) and 0.74 (red dots). The diagram shows the function F multiplied by 10 (blue) and 100 (red) for better visibility. The straight lines best fitting these discrete sets of data-points are indicated as solid curves (the curvature is a result of the log-log scale) in matching colors. The intercept A is clearly lower than 1, increased with decreasing K max and asymptotically reached 1 when K max dropped below 1.

Numerical Procedures
The aim of this work is not to verify/falsify the single pipe models mentioned above, but to investigate their extension to heterogeneous and partially connected porous media. The network simulation approach is most certainly indicated for this purpose. Since the BK model has received greater support from experimental and numerical simulation data than other contending models, I elected to use it (i.e., equations 5 and 6) in my network simulations. The numerical procedures employed here are based on the Bernabé et al. 's 30 workflow template for studying liquid flow and ionic conduction through heterogeneous networks. The numerical procedures followed in the present study are briefly sketched below (more details about the Bernabé et al. 30 procedures can also be found in previous related papers 17,31 ).
Construction of network realizations. Two kinds of network realizations were constructed, corresponding to monomodal and bimodal radii distributions. In monomodal realizations, the radii R i of individual pipes were randomly assigned according to log-uniform distributions such that the hydraulic radius of the network R H honored pre-selected values (note that the mean pipe radius is not equal to the hydraulic radius; see Bernabé et al. 30 for details). Since the Knudsen number depends on pore size, I considered five values of R H , namely 30, 10, 3, 1 and 0.3 10 −6 m (the pipe length was L = 10 R H ). The coefficient of variation of the distribution, CV, is a convenient measure of the pore size variability. Here, I prepared network realizations ranging from nearly homogeneous (CV = 0.05) to extremely heterogeneous (CV = 1.05). I also varied pore connectivity by randomly removing a fraction of the pipes and thus decreasing the coordination number z (mean number of connected pipes per nodes in the network). The coordination numbers considered here ranged from 3 to 6, corresponding to half and fully filled SC networks, respectively. The critical coordination number z c at the percolation threshold is approximately equal to 1.5 in all three-dimensional networks 32 . For all combinations of R H , CV and z, multiple 12 × 12 × 12 and 16 × 16 × 16 network realizations were constructed. For the bimodal realizations, two populations of pipe with very narrow radii distributions (CV = 0.05) were considered. The set of large pipes had a hydraulic radius ten times greater than that of the small pipes (the pipe length was equal to ten times the hydraulic radius of the large pipes population). Five combinations of large and small hydraulic radii was used, i.e., 30/3, 10/1, 3/0.3, 1/0.1 and 0.3/0.03 10 −6 m. In each realization, the individual pipes were randomly assigned one or the other distribution as specified by pre-selected number fractions of large and small pipes, w L and w S = 1 − w L . The pipe radii were then accordingly drawn. No pipes were removed, implying z = 6 in all bimodal realizations. However, the separate coordination numbers of large and small pipes, z L and z S , varied between 0 and 6 depending on w L .

Gas flow simulations.
Simulating fluid flow through pipe networks consists in solving the Kirchoff equations, i.e., to impose local conservation of mass at each node i as ∑ α  m a i = 0, where the subscript α refers to the nodes connected to i and the mass fluxes can be expressed using the appropriate one among equations 1-5 (note that, for the sake of simplicity, pure Poiseuille flow was assumed and entrance flow corrections were not included in the Kirchoff equations; the pipe length was generally chosen long enough to validate this assumption). One important difference of gas and incompressible liquid flow is that the laws of transport through individual bonds are non-linear in the case of gas flow, whereas they are linear for incompressible liquids. As a consequence of the non-linearity, the Kirchoff equations are much harder to solve for gas than liquid flow simulations. Limitations in computer power prevented me from utilizing the high-performance, parallelized, successive over-relaxation (SOR) iterative solver of Li et al. 17 . Instead, I implemented a simple, basic version of the relaxation method. The slowness of my iterative non-linear solver imposed strong restrictions on the number and size of the simulations that I could be run in a realistic time frame. In particular, I could not carry out the Bernabé et al. 30 technique of using different types of networks (i.e., simple cubic, BCC and FCC) to test the generality of the results. Only simple cubic simulations were performed here (note that Li et al. 17 obtained results from BCC and FCC simulations quite consistent with the simple cubic ones).
The so-called permeameter boundary conditions were used (constant pressures applied to the entry and exit faces, no flow allowed through the sides). To test the accuracy of the solver I compared the intrinsic permeability simulated using a gas in the purely viscous regime (equation 1) and the permeability k L of an incompressible liquid through identical network realizations (in liquid flow simulations, the linear Kirchoff equations were solved using the Krylov method). In principle, one would expect k L and k to be equal, but I found that the solver yielded values of k greater than k L . The difference between k and k L depended on the pore radius variability measure CV and, to a lesser extent, to the pore coordination number z. The ratio k/k L ranged from 1 to 1. Although the solver used here was accurate in nearly homogeneous networks, the accuracy could be deteriorating in networks with increasing heterogeneity, thus casting doubts about the validity of the k/k L > 1 inequality. Alternatively, the increase of k/k L with increasing heterogeneity may reflect a real gas flow property. This issue is difficult to resolve. Since the accuracy of the solver depends on the assumed initial pressure field (IPF), I tried several possible IPF's, including the pressure field predicted by equation 1 for a perfectly homogeneous network (R i = R H for all pipes) or that produced by the liquid flow simulation. The diverse solutions obtained did not differ by more than 1 or 2 percent. I also tried to change the solver arrest condition in purely viscous gas flow simulations and found that the intrinsic permeability increased with tightening of the arrest condition. A tighter arrest condition presumably leads to improved accuracy, albeit at a great cost in CPU time. Thus, the observed increase of k/k L with increasing heterogeneity is in fact more likely to be a valid result than not. In any event, to avoid including potentially spurious results in the subsequent analyses, I discarded the highly heterogeneous network realizations, for which k/k L was larger than 2 (i.e., the orange data-points in Fig. 2).
Since the Knudsen number varies with gas pressure, gas flow for each network realization was simulated at 6 different mean macroscopic gas pressures P (namely, 10, 0.1, 0.01, 0.003, 0.001 and 0.0003 MPa). These variations in mean gas pressure combined with those of the hydraulic radius produced more than 6 orders of magnitude changes in Knudsen numbers, K n = (μ/P R H ) (πR g T/2) 1/2 . The maximum Knudsen numbers K max achieved for each networks realization were 74, 22, 7.4, 2.2 and 0.74 corresponding to R H = 0.3, 1, 3, 10 and 30 10 −6 m, respectively. Varying P also allowed performing the Klinkenberg analysis for each network realization (i.e., examining the dependence of the gas permeability k gas on the inverse gas pressure 1/P ).
Two different rules were used to set the macroscopic pressure difference ΔP, the small-gradient rule, ΔP = P /10, and the large-gradient rule, ΔP = 2(P − 0.0001). The large-gradient rule obviously produced much greater variations in local Knudsen numbers than the small-gradient rule, especially when the mean pressure P was high (it also yielded slightly higher k/k L ratios). Only the small gradient-rule was used with the bimodal network realizations. All the simulations assumed that the saturating gas was nitrogen at room temperature. For simplicity, the compressibility factor Z was taken to be equal to 1 and full momentum accommodation was considered (σ = 1).

Results
Monomodal networks. For each network realization, values of k and k g were calculated using equation 1 and equations 5 and 6, respectively. The simulated k g 's corresponding to different mean gas pressures P can conveniently be expressed in normalized form, k g /k, and displayed in Klinkenberg-type plots (Fig. 3). Linear functions, k g /k = A* + B*/P , appeared to fit the simulated data very well in all cases, although a subtle underlying non-linearity was indicated by the dependence of A* and B* on the maximum Knudsen number reached for each k g /k curve. For comparison with F, the normalized gas permeability can be recast as a function of the Knudsen number where A and B depend on K max . The values of A and B in nearly homogeneous networks (CV = 0.05) were very close to those estimated from the function F and, furthermore, were not significantly affected by changes in connectivity (Fig. 4). This good agreement of k g /k and F in homogeneous networks was not preserved, however, for higher levels of heterogeneity. In network realizations with CV ≥ 0.55, A and B for k g /k increasingly diverged from their F counterparts when CV was increased and/or z decreased. But the most striking result was that the intersect A had completely different forms when the small-and large-gradient rules were used. With the small-gradient rule, the results were analogous to those shown in Fig. 1, namely, A was lower than 1 and decreased with increasing K max . In contrast, values of A larger than 1 and increasing with K max were obtained when the large-gradient rule was applied (Fig. 4). An increase of A with K max indicates a non-linear Klinkenberg curve with a downward curvature and implies an overestimation of k if the Klinkenberg extrapolation is used. The implication is thus that experimental Klinkenberg plots may yield incorrect underestimates or overestimates of the intrinsic permeability if the maximum Knudsen number investigated is larger than about 1 (a limit suggested by Fig. 4), depending on the heterogeneity level and the magnitude of the pressure gradient.
Since the Klinkenberg analysis can be properly performed when K max is lower than 1, I calculated the Klinkenberg coefficient b K for the various network realizations excluding the k g /k data obtained for K n > 1 (Fig. 5). In simulations using the small-gradient rule, b K was approximately proportional to the inverse hydraulic radius, b K ≈ C*/R H , with C* a relatively weak linear function of CV (specifically, C* = 0.034 + 0.2CV). Importantly, there was no discernable effect of connectivity (the data-points fell on nearly horizontal lines in constant R H planes; Fig. 5). The large-gradient rule produced similar results with slightly larger and more variable values of b K (Fig. 5).
Bimodal networks. The intrinsic permeability k generally decreased with increasing w S and a sharp drop occurred when the set of large pipes became disconnected (i.e., for w L ≈ 0.25 or w S ≈ 0.75). The variations of k with w S , including the singularity observed above the large pipes percolation threshold (for w L ≥ 0.25 or w S ≤ 0.75; Fig. 6a), were similar for all hydraulic radii combinations and were well reproduced by a simple modification of  1.05 and z = 3). The small-gradient rule was used in (a) and the large-gradient rule in (b). The values of R H and K max corresponding to each set of simulations are indicated in the insets in colors matching the data-points. The color-matching solid lines represent the best fitting linear functions, k g /k = A* + B*/P . These lines show that the intercept A* was lower than 1 when the small gradient was used and greater than 1 otherwise. the percolation-based binary mixing model of Bernabé et al. 33 . The Bernabé et al. 33 model uses different mixing laws above and below the percolation threshold of the high-permeability phase, namely, the upper and lower Hashin-Shtrikman bounds for w L > 0.25 and w L < 0.25, respectively. A better fit was obtained here by replacing the lower Hashin-Shtrikman bound with the geometric average (see details in Appendix A).
The Klinkenberg coefficient b K also displayed a typical percolation behavior. For each hydraulic radii combination, b K remained nearly constant until w S approached the percolation threshold, where a sudden jump occurred (Fig. 6b). Strikingly, the singularity was located below the large pipes percolation threshold (for w L ≤ 0.25 or   Fig. 6b) and not above as normally expected. Thus, the connectivity of the large pipes significantly affected k but had a negligible influence on b K (this is similar to the lack of sensitivity of b K to z − z c in Fig. 5). Below the percolation threshold, the (disconnected) large pipes still strongly affected b K while their effect on k was weak.

Discussion
Upscaling of the BK single pipe model. The relationship of k g /k and K n obtained in the simulations, in effect, represents the result of upscaling the BK model from the scale of a single pore to that of a porous body containing many pores. One manifest result is that the simulated macroscopic behavior was not always similar to that predicted by the BK model. The simulated data revealed very weakly non-linear Klinkenberg curves, asymptotically becoming linear at high Knudsen numbers, as expected from the BK model. But both upwardand downward-bending of the Klinkenberg curves were observed whereas equations 5 and 6 only show upward curvature. Indeed, rare but clear instances of upward or downward bending of experimental Klinkenberg curves have indeed been observed [34][35][36] . However, the subtle non-linearity predicted by the simulations could often be obscured by experimental uncertainties.
The main factor producing weakly non-linear Klinkenberg curves in the simulations was the heterogeneity of the field of local Knudsen numbers, which itself resulted from the geometrical heterogeneity of the networks combined with the gas pressure difference imposed externally. The origin of the strong non-linearity occasionally observed in rocks is less clear. Rarefaction effects such as those predicted by the BK model must have played a role but were probably not the only cause since Knudsen numbers larger than 1 were rarely reached in these experiments (see Fig. 9 of Sinha et al. 34 or Figures 14 and 15 of Wang et al. 36 ). Other causal factors include deformability of the porous material, which may have a strong effect at high gas pressures (for example, in Fig. 7 of Wang et al. 36 the increase of k g as the inverse gas pressure approaches 0 was likely due to a decrease of the effective pressure; see also Letham and Bustin 37 ), the non-ideal constitutive behavior of the gas 3 and gas adsorption in porous media with a very small pore size 1 .

Difference of liquid and intrinsic gas permeability.
Another important result is that the intrinsic permeability k estimated from purely viscous gas flow simulations was apparently larger than the liquid permeability k L . As discussed in section 2, there is a reasonable possibility that the enhanced efficiency of compressible flow with respect to incompressible flow was not a numerical artifact but resulted from genuine differences between the pressure fields produced in identical networks by flow of a compressible gas and an incompressible liquid. To illustrate this point, the normalized gas and liquid pressure fields, P i and Π i , obtained for the same network realization (R H = 30 10 −6 m, CV = 0.55 and z − z c = 4.5) were plotted against each other in Fig. 7a. Two versions of P i corresponding to the small-and large-gradient rules are shown in the same diagram. The effect of gas compressibility is evident for the large-gradient dataset (quasi-quadratic P i versus Π i trend) while it is almost invisible in the small gradient case (quasi-linear trend). Most importantly for this discussion, the relationship between P i and Π i is blurred by strong fluctuations with respect to the overall trends. These fluctuations are irregular but not random as they resulted deterministically from the specific R i field of this network realization. Under close examination, identical very strong fluctuations can easily be identified in the large-and small-gradient data-point clusters, implying that the geometrical heterogeneity of the network affected the two pressure fields in a nearly indistinguishable way. This view is supported by a cross-plot of the large-and small-gradient gas pressure fields, precisely outlining a fluctuation-free, compressibility-controlled relationship among the two pressure fields (Fig. 7b). Similar observations were made for network realizations with larger CV and/or smaller z -z c , although the level of numerical noise was higher than in Fig. 7.
Ratios of Klinkenberg-corrected gas permeability to liquid permeability, k/k L , significantly greater than 1 have been experimentally observed in crystalline rocks (Christian David and Jerome Wasserman, personal communication), in tight sandstones 38,39 and other materials such as clay-bearing fault gouge 40 or mortar 41 . The origin of the experimental k > k L discrepancy remains unclear. The nature of the saturating liquid or gas appears to be very important. Permeability to water tends to be lower than permeability to non-polar liquids 41 , suggesting that water flow tests can be strongly affected by water-solid interactions (note that distilled water was observed to reduce permeability more than brine 38,40 ). Gas flow experiments showed that permeability to nearly ideal gases such as helium is substantially higher than the permeability to non-ideal gases such as methane and that the difference increases with non-ideality 3,34 . Effect of the pressure gradient. Besides the effects discussed above, the magnitude of the pressure gradient also affected the gas pressure fields simulated in conditions corresponding to moderate to high Knudsen numbers. Cross-plots of the large-and small-gradient pressure fields simulated in the same realizations are shown in Fig. 8 for P = 0.01 and 0.0001 MPa and various input parameters (CV = 0.55, z -z c = 4.5 and R H = 30, 10, 3, 1 and 0.3 10 −6 m). The Knudsen numbers achieved in these simulations increased from 0.022 to 22 with decreasing R H and P . As in Fig. 7b, the cross-plotted curves show very small fluctuations about the well-delineated overall trends. The important observation is that these curves reveal a regular transition with decreasing K n from quadratic to nearly linear trends (Fig. 8). Thus rarefaction and compressibility effects tend to balance each other in gas flow at large Knudsen numbers. should be much smaller than the confining pressure to reduce variations in effective pressure and the associated pore space deformations, (c) the lowest mean gas pressure should correspond to a Knudsen number no greater than 1, and finally, (d) substantially non-ideal gases should be avoided (e.g., methane, carbon dioxide).

Implications for the
The relationship between k and b K is also an issue of great interest in rock physics. Experimental data suggest a power law relationship, b K ∝ k −n . Values of the exponent n between 0.33 and 0.39 were obtained from large rock datasets 38,42 . However, these values can be considered uncertain owing to the very large scatter clouding the data. Wang et al. 36 measured k and b K in two gneiss samples subjected to increasing confining pressures. They found very well defined power laws with n = 0.20 and 0.53, demonstrating that the relationship between k and b K is highly variable for individual rocks. Civan 22 derived a model predicting a slightly different power law, b K ∝ (k/φ) −n , where φ denotes porosity and n = 1/2 (Sampath and Keighin 39 applied this model to their data and found n = 0.53). The values of k and b K obtained here in monomodal simulations are represented in a three-dimensional diagram against each other and the hydraulic radius R H of the network realizations (Fig. 9). The diagram shows that, as already mentioned in sections 3 and 4, k and b K are very sensitive to the hydraulic radius and generally obey the power laws, k ∝ R H 2 and b K ∝ 1/R H (the dotted lines shown in Fig. 9 for comparison were calculated such that their projections on the k − R H and b K − R H planes obeyed these power laws). The intrinsic permeability is also significantly affected by pore size heterogeneity CV and pore connectivity z − z c while the Klinkenberg coefficient is weakly affected by CV and almost totally insensitive to z − z c (Li et al. 17 obtained very similar results using the Javadpour pipe model). Combining these power laws, one easily finds b K ∝ k −1/2 . However exponents lower than ½ as observed in natural rocks 38,42 can be obtained if subsets of the simulated data are appropriately selected. Reduction of n is easily achieved by introducing some sort of correlation of k to z − z c and anti-correlation to CV, a reasonable assumption since highly permeable rocks, indeed, possess better pore connectivity and tend to show less pore size relative variability than those with very low permeability.  Bimodal networks. The bimodal networks investigated here are not realistic representations of rocks containing micro-porosity (a more truthful image is given by dual porosity networks 43 ). The results obtained are nevertheless quite instructive. Both the intrinsic permeability and Klinkenberg coefficient displayed typical percolation singularities when w L crossed the percolation threshold of the large pipe population (Fig. 6), although the behaviors of k and b K were very different. The intrinsic permeability simulated data are very well reproduced using two different binary mixing laws, the upper Hashin-Shtrikman bound above (i.e., to the left of the grey line in Fig. 6a) and geometric averaging below the percolation threshold (see Appendix A). In this model, the percolation singularity results from the critical power law describing the number fraction of large pipes that belong to the connected cluster (here, I used the classic value, 0.41, of the critical exponent). In the case of the Klinkenberg coefficient, however, the singularity is located below the percolation threshold (i.e., to the right of the grey line in Fig. 6b). Remarkably, b K appears to be constant, independent of w L , above the percolation threshold. This peculiar behavior may be related to the insensitivity of b K to z -z c observed in monomodal simulations. This result suggests that, in dual porosity rocks, b K may be insensitive to the presence of micro-porosity as long as the macro-pore population is connected.

Conclusions
(a) The flow of an ideal gas through heterogeneous and imperfectly connected simple cubic networks of pipes was simulated numerically. The simulations included both rarefaction and compressibility effects. The permeability to gas in the purely viscous regime (K n = 0) was found to be greater than the permeability to an incompressible liquid. This result suggests that a compressible fluid may flow through a heterogeneous porous media more efficiently than an incompressible one. Non-ideal constitutive laws may have an important additional effect. (b) The functional dependence of gas flow on the macroscopic Knudsen number differed from that predicted by the BK single pipe model and was contingent on the pipe radius distribution, the pore connectivity and the magnitude of the externally applied pressure gradient. The implication in terms of the Klinkenberg analysis was that extrapolation of a Klinkenberg curve cannot be trusted to provide an accurate estimate of the intrinsic permeability unless the maximum Knudsen number investigated is lower than 1. (c) The Klinkenberg analysis was applied to appropriate subsets of the simulated data (K max < 1). The Klinkenberg coefficient b K is almost completely insensitive to pore connectivity, only moderately affected by the width of the pipe radius distribution, and thus nearly proportional to the inverse hydraulic radius of the network. (d) In the bimodal simulations, the intrinsic permeability k displayed a typical percolation behavior, with a singularity occurring immediately above the percolation threshold of the large pipe population. The Klinkenberg coefficient b K showed an unusual behavior. It remained constant above the percolation threshold and developed a typical percolation singularity below it.