A Stochastic Two-Scale Model for Rarefied Gas Flow in Highly Heterogeneous Porous Media

This paper presents the development of a model enabling the analysis of rarefied gas flow through highly heterogeneous porous media. To capture the characteristics associated with the global- and the local-scale topology of the permeable phase in a typical porous medium, the heterogeneous multi-scale method, which is a flexible framework for constructing two-scale models, was employed. The rapid spatial variations associated with the local-scale topology are accounted for stochastically, by treating the permeability of different local-scale domains as a random variable. The results obtained with the present model show that an increase in the spatial variability in the heterogeneous topology of the porous medium significantly reduces the relevance of rarefaction effects. This clearly shows the necessity of considering a realistic description of the pore topology and questions the applicability of the results obtained for topologies exhibiting regular pore patterns. Although the present model is developed to study low Knudsen number flows, i.e. the slip-flow regime, the same development procedure could be readily adapted for other regimes as well.


Introduction
When the density of a gas is very small, or when it is forced through a very small channel, the gas can no longer be described using a continuum representation. In this case, one says that the flow is rarefied. Knudsen (1909) observed this rarefaction effect already in 1909, when he studied the flow of a gas driven through a tube at low pressure. He found that the flow rate differed from the one predicted by the Poiseuille formula, which is based on theory governing flow of continuous media. Another situation where a gas does not behave as a continuum is when it flows through channels where the mean free path of the gas molecules is comparable to the dimensions of the channel. The classification of rarefaction, in terms of the Knudsen number, K n , which is defined as the ratio of mean free path to the characteristic dimensions of the channel, presented in Barber and Emerson (2006) tells us (i) that the flow is rarefied when K n ≳ 10 −3 , (ii) when 10 −3 ≲ K n ≲ 10 −1 , the flow is in the slip-flow regime and continuum theory incorporating the slip flow boundary condition is applicable, (iii) for 10 −1 ≲ K n ≲ 10 1 , the flow is in the transition regime in which continuum theory is not applicable, and (iv) for K n ≳ 10 1 , the flow is said to be in the free molecular regime, which admits the collisionless Boltzmann equation to govern the flow behaviour.
In a porous medium, the fluid percolates through narrow channels and even more narrow constrictions, at which the pressure, consequentially, drops markedly. This means that gas rarefaction likely occurs when gas flows through a porous medium, and this is an event occurring in many applications such as gas production, storage and separation, catalytic reactions in low-porosity ceramics and filtration in medical applications (Johansson et al. 2018). In general, quantified prediction of rarefied flow through a given porous medium requires a realistic and thus highly detailed representation of the flow domain. Recently, new imaging techniques have made it possible to capture the topology of such domains (Maire and Withers 2014). Unfortunately, the models used to study rarefied gas flow become computationally expensive to solve even for well-described and simplistic types of topologies. This is also even more pronounced if the transition to rarefied flow is resolved with high accuracy by using kinetic theories based on the Boltzmann equation (Kalarakis et al. 2012;Pavan and Oxarango 2007) and its derived methods, such as the Bhatnagar-Gross-Krook kinetic equation (Germanou et al. 2018;Wu et al. 2017). Examples of some topologies that have been investigated can be found in Lasseux et al. (2014) where regularly packed balls were considered and in Colin (2005), where well-defined patterns of micro-channels are considered. The relevance of constrictions was highlighted in Zheng et al. (2016), studying the regular multi-scale Sierpinski carpet type of topology. Related to this, several studies have investigated the effect of rarefaction in topologies defined by simple regular patterns (Yang and Weigand 2018) and topologies with self-affine fractal nature (Zheng et al. 2013;Wu and Chen 2016). It has also been shown that the effect of rarefaction, i.e. how much does the flow differ from what would be expected from a continuum theory, cannot be explained by a simple parameter such as a representative pore size. Instead, it also depends on porosity (Meghdadi Isfahani 2017), heterogeneity (Zheng et al. 2017;Germanou et al. 2018) and anisotropy Germanou et al. 2018), which doubtlessly complicates the analysis of rarefied gas in porous media. This underlines the importance of considering the structure of the pore topology, which significantly affects the behaviour of the flow (Yang and Weigand 2018), preventing a universal description to be developed (Kalarakis et al. 2012).
There are a few recent publications considering rarefied gas flow through porous media with heterogeneous topology. For instance, Germanou et al. (2018), Wu et al. (2017) and Lasseux et al. (2016) considered topologies consisting of a randomly packed collection of two-dimensional balls, Kawagoe et al. (2016) considered a random packing of spheres, Naraghi and Javadpour (2015) studied randomly generated topologies, and Wang et al. (2016) considered topologies reconstructed from measurements. The objectives and main contributions of the present work have, therefore, been specified as (i) to investigate in detail the importance of considering the random variations of the topology of the porous medium, both for qualitative and for quantitative predictions, when studying rarefied gas flow through it, and (ii) to propose a modelling framework that permits adapting existing models to account for these variations. To this end, a stochastic two-scale model is proposed here.
Conceptually, a good starting point for developing a two-scale model can be found in the successful use of the volume averaging method in the study of rarefied gas flow in porous media; see e.g. Lasseux et al. (2014Lasseux et al. ( , 2016 and Jin and Chen (2018). To find examples in which a heterogeneous porous medium has been explicitly considered, however, one must resort to the more prolific studies of the random nature of the non-rarefied (pressure driven) flows in porous media; see, e.g. Ricciardi et al. (2005), Nordlund et al. (2006), Selvadurai and Selvadurai (2010), Guilleminot et al. (2012) and Selvadurai et al. (2014). The advantage of using these techniques is that they are well tested, and thus we know that they are correct. As most of the available two-scale models, these are formulated as a combination of a global problem considering the macro-scale features and a local problem considering the micro-scale features. In this way, the global problem can be discretised using a coarse grid, which only needs to resolve the long wavelength variations induced by the macro-scale topology. This does, however, require permeability values obtained by solving problems defined at small, highly resolved domains, including the micro-scale details.
When working with highly heterogeneous porous media, however, these "small" domains must actually be rather large if they are expected to be representative of the whole domain of study. Unfortunately, when discretised, this may lead to a system of equations too large to be manageable in practice. This can be circumvented, e.g. by applying the approach presented in Nordlund et al. (2006), Guilleminot et al. (2012) and Pérez-Ràfols et al. (2016), which also is adopted here. In this approach, not one, but a number of small local-scale domains are sampled even though these, individually, are not considered to represent the statistics of the original topology. The variation within the sampled set of domains is included by considering their permeabilities as a random variable. This is the way the influence of the micro-scale topology which, spatially, varies significantly in permeability is treated in this type of approach. This means that the effect of the highly heterogeneous porous medium can be efficiently homogenised and included in the global problem, which solution can be easily computed. This does, however, imply that a quite large number of local-scale domains must be considered before convergence of the distribution can be observed, and this makes it a rather demanding problem numerically.
Since the objective of this work is to investigate the importance of considering the random variations of the topology of the porous medium, and not to improve the fundamental description of rarefied gas flow, only the most simplistic underlying physics will be considered. For the model developed here, we will, therefore, consider the slip-flow regime 10 −3 ≲ K n ≲ 10 −1 . Indeed, with a model constructed in a similar way as the one in Lasseux et al. (2016), we can obtain coefficients, e.g. permeability, intrinsic to the topology, which can be interpreted in a stochastic manner. We stress, however, that a similar twoscale stochastic framework could be constructed from other local-scale models that, for instance, describe the slip flow differently, concern other types of rarefied flow or even include more complex physics; see e.g. Moyne and Christian (2003), Charrier and Dubroca (2004), Wu et al. (2017) and Germanou et al. (2018).
The remainder of this paper is organised as follows: firstly, in Sect. 2, a reference problem to aid the presentation is introduced. Then, the underlying physics and the mathematical model will be presented in Sect. 3. Thereafter, in Sects. 4 and 5, the proposed two-scale model and its stochastic extension will be introduced. Lastly, the interplay between rarefaction and variability of the topology of the porous medium will be discussed in Sect. 6, and a closure of the paper will be given in Sect. 7.

Benchmark Topology
For the numerical simulations performed by means of the stochastic two-scale model developed in this work, we have chosen a collection of randomly packed balls, to represent a benchmark for topology of heterogeneous porous media. To this end, we have chosen to use the topology depicted in Fig. 1.
This particular collection of balls was measured with X-ray microtomography, using a Zeiss Xradia Versa 510 system. The sample consists of 200 steel spheres with diameter 2.0 mm randomly packed inside a PMMA tube with inner diameter 14 mm. The full height of the bed is 10.2 mm, resulting in a packing ratio of 0.53. The steel spheres (balls) are designed for use in bearings and has a well-determined uniformity and size distribution specified by the manufacturer. The scanning was carried out with an X-ray tube voltage and output of 160 kV and of 10 W, respectively. The field of view (FOV) was 15.2 mm, and the spatial resolution in terms of voxel size 14.9 μ m. The total number of projection images was 1601, acquired with and exposure time of 14 s, resulting in a total scan time of approximately 7 h. As stated above, the diameter of the balls is 2.0 mm. We have, however, applied a scaling of the topology to match the requirements of the current investigation. More precisely, we have taken samples from the measurement depicted in Fig. 1, to construct the flow domains for the local-scale problem presented in Sect. 4.2. For the numerical analysis of the effects of rarefaction and variability of the topology, the samples were scaled so that the side l of the cubical local-scale domain is equal to the radius of the balls. A motivation for this choice will be given in Sect. 4.
A collection of randomly packed balls was chosen as a benchmark topology because it is simple, but still exhibits the random spatial variation expected in a heterogeneous porous medium. As such, it can be used to illustrate the applicability of the model and it also makes it possible to study the influence of the spatial variability of the topology on the effect of rarefaction. It should, however, be emphasised that the modelling framework developed herein can be applied to study the flow of rarefied gases through other more complex types of heterogeneous porous media, provided that the due assumptions are met. Moreover, when considering more complex geometries, the random variations exhibited by the porous medium can only be expected to be more influential than in the simpler geometry considered here.

Modelling Rarefied Gas Flow
Let us consider the steady-state regime of (slightly) rarefied gas flow with Knudsen numbers 10 −3 ≲ K n ≲ 10 −1 . To this end, the Stokes momentum equations for compressible flow and the continuity equation will here be combined with a first-order slip-flow boundary condition. Indeed, where we have assumed that inertia forces are small compared to the viscous ones, that the flow is in steady-state, isothermal conditions and that the fluid can be modelled as an ideal gas with constant viscosity ( ) and obeying the well-known constitutive relationship (2) between density and pressure where M is the molar mass in kg/mol, T is the temperature in K and R = 8.314J∕(molK) is the universal gas constant. In (1), p, v , and are the pressure, velocity, density and dynamic viscosity of the fluid, which occupies the volume V f , and = 2 − v ∕ v , where v is the accommodation coefficient. In the first-order slip-flow boundary condition, at the interface sf between the solid and the fluid (1c), is the mean free path of the molecules given by where k B is the Boltzmann constant, T is the temperature and is the collision diameter. Moreover, the unit normal vector at sf is denoted as n and indicates the identity matrix.
Together with adequate boundary conditions, (1) can be solved for the pressure (p) and the velocity (v), which approximately describe the state of the rarefied gas flow. Due to the small size of the pores, however, a very fine mesh would be required and solving the problem becomes extremely computationally intensive, even when considering fairly small domains. For this kind of situations, a two-scale formulation can provide for the solution. This will be the topic of the following section.

The Two-Scale Model
In this section, a two-scale model is built upon the rarefied gas flow model, introduced in the previous section. This two-scale model will, as we will see, form the basis for the twoscale stochastic model, presented in Sect. 5. The present model will be based on the heterogeneous multi-scale method (HMM), which is a fairly general framework that can be used to construct two-scale models for various different physics; see e.g. E and Engquist (2003) for details. The reason for selecting this framework before more rigorous and well-established methods such as the volume averaging method employed in Lasseux et al. (2016) or the two-scale homogenisation technique, see e.g. Lukkassen et al. (2002), is that it is as readily applied and more flexible at the same time. Indeed, the proposed modelling framework can be adapted to other underlying base models, by following the same procedure as we present here. Also, it avoids certain assumptions that would be in contradiction with the stochastic part of the model introduced in Sect. 5. For instance, the two-scale homogenisation technique demands periodicity and the volume averaging method requires the localscale domain to be large compared with a typical size of the pores. We note, moreover, that this type of framework has previously been applied to a similar problem; see Pérez-Ràfols et al. (2016). We will, however, show that the resulting model in many ways is equivalent to the model derived by Lasseux et al. (2016), which is based on the more rigorous volume averaging method.
A schematic illustration of how the present two-scale model is constructed is depicted in Fig. 2. As usual, when constructing a two-scale model, global-and local-scale domains are defined. The global-scale domain covers the whole topology of interest, and, as indicated by the red grid, a coarse mesh is applied. Effectively, this means that only the large features of the topology are considered here, whereas the fine details of the topology are considered at the significantly smaller local-scale domain, with mesh resolution high enough to capture them.
Following the HMM framework, two steps are taken here to couple the global and the local scale. That is, (i) the global-scale problem is defined with a set of coupling variables-in this case permeabilities which summarise the effect of the fine details of the topology in an averaged sense and (ii) these permeabilities are computed based on the solution of the local-scale problem, defined on a small highly resolved local-scale domain. The global-scale model governs the fundamental physical nature, and in the case of fluid flow, it is a standard practice to employ a mass conservation criterion. When constructing the local-scale problem for the present work, the gas rarefaction physics is described by combining the Stokes momentum equation with the slip-flow boundary condition and the continuity equation, i.e. (1). One must, however, pay special attention when describing the boundary conditions for the local-scale domain so that the interfaces between the local-and the global-scale domains are described in a consistent manner.
In the following subsections, the two aforementioned steps will be developed to model rarefied gas flow through heterogeneous porous media. To this end, it will be assumed that L >> l , where L is the characteristic size of the global-scale and l characterises the size of the local-scale domain. This results in that the global-scale variables vary slowly and that they, therefore, can be assumed constant inside each localscale domain. In the majority of the available two-scale models for flow through porous media, it is also required that l >> d p , where d p represents the characteristic diameter of the pores. This is to ensure that the results obtained from two different local-scale domains would be similar. One of the advantages with the present stochastic two-scale approach, which is based on the HMM, is that rather large variations in the local-scale response (due to the stochastic nature of the topology) can be considered and having l ≳ d p is, therefore, sufficient. The numerical investigation includes a study of the rarefaction effect associated with the size of the channels that the voids between the balls constitute. This is done by analysing the effect of l, and the results are presented in Sect. 6. When zooming in, however, one sees that there is a rapid variation of the topology of the solid phase (dark areas) and fluid phase (light areas). In this zoomed image, a coarse global-scale mesh is defined, as depicted by the coarse grid in red. An example of a local-scale domain, V f , centred at ijk , is highlighted in green. A three-dimensional representation of this domain is depicted on the right. The mass fluxes Q I , in all three directions I = 1 , I = 2 and I = 3 , from and towards the neighbouring domains, are depicted with yellow arrows. Note that the directions and magnitudes of the arrows in the figure exemplify a given flow situation

The Global-Scale Problem
As stated previously, in the case of fluid flow, the associated global-scale problem is normally expressed as an equation describing conservation of mass. This is also the case here, and the discretised form is formulated by means of the finite volume method. To this end, a coarse grid is marked in red in Fig. 2, in which the elements are the local-scale cells where the mass-flow balance is imposed. More precisely, for the steady-state case considered here, mass-flow balance for a local-scale cell V f , centred at ijk ∶= X 1i , X 2j , X 3k , will be equated as where Q I indicate the net in-and out mass fluxes of the local-scale domain V f centred at ijk , in each direction I = 1, 2, 3 . In order to compute these mass flow components, a law of the form of Darcy's law is used. There are, however, various ways that can be applied to discretise such law, and for the present analysis, we have adopted the one in (5), in which only the diagonal terms of the permeability matrix are included, i.e.
where P ijk is the global-scale pressure defined as P ijk ∶= P ijk , K I ijk is the permeability in the X I -direction of the local-scale domain V f , centred at ijk , coloured in green in Fig. 2, ijk ∶= (P ijk ) the fluid's density at the pressure P ijk and l is the size of the domain. We note that we have here assumed that the domain is a cube with all three dimensions equal. Otherwise, the term l should be replaced by A o ∕l p , where A o is the area of the outlet and l p is the dimension in the direction of the pressure gradient.
The particular choice of law in (5) does not contain the off-diagonal terms ( K IJ , I ≠ J ) of the permeability tensor, meaning that a pressure drop in a given direction does not contribute to flow in the directions perpendicular to it. The reason for only considering the permeabilities in the main directions (by including only the diagonal terms K I ∶= K II ) is related to the construction of the local-scale domains, which is the same as in the model presented in Patir and Cheng (1978). A more detailed justification and an explanation of the implications of this choice is given in the end of Sect. 4.2.2. Note, moreover, that, despite we use a law of the form given in (5), we do not assume that Darcy's law per se applies, as we will allow the permeability to be vary with pressure due to rarefaction effects.
By approximating the values of the pressure at the mid-grid nodes i± 1 2 ij with the following expression and analogously for the other directions, the fluxes in each direction in (5) become An equation for the global-scale pressure can now be obtained from (4) and (7) In the following section, we will describe and define the local-scale problem and how to obtain the fluxes Q I ijk and, thereafter the permeabilities K I ijk , so that we can solve (8) for the global-scale pressure P ijk . We note that, since we will explicitly allow permeability to vary from cell to cell, the effects of porosity and velocity variations between cells will be implicitly considered, and Brinkman correction terms will not be necessary.

The Local-Scale Problem
In this section, the local-scale problem, that needs to be solved in order to obtain the permeabilities K I ijk from (7), which are required to solve (8) for the global-scale pressure P ijk , will be defined. But, before we present the details for this, we will elaborate on the requirements for obtaining K I ijk .
At the global scale, each component K I ijk of K I represents the permeability in the X I -direction of the local-scale cell V f centred at ijk . The spatial variation of the topology implies that obtaining K I ijk would require solving as many local-scale problems as there are local-scale cells V f . In Sect. 5, we will present the stochastic model, which presents an effective way of avoiding this problem. Moreover, it should be noted that (7) is not a simple linear equation, and an approximation that the density's variation is small enough to consider it constant within each local-scale domain V f , i.e. ijk will, therefore, be applied. This leads to that K I ijk (approximately) can be determined as where I P ijk = P i−1jk − P i+1jk , with Q I ijk computed from the solution obtained by solving the local-scale problem (12). This implies that K I ijk would depend on the pressure P ijk as well as ijk and . For a fluid which is iso-viscous and almost incompressible, a combination of a non-dimensionalisation and a Taylor expansion of the permeability can be adopted to alleviate the need to compute K I ijk for all the values of I P ijk and ijk for the global-scale grid employed.

3
Let us now consider a full three-dimensional local-scale domain, such as the one illustrated to the right in Fig 2, which would be used to compute K I ijk . Note that the boundaries i , o and l , as specified in Fig 2 (right), are for the flow in the x 1 -direction. Since the permeabilities at other locations and in other directions are computed analogously by choosing the appropriate domain, the presentation will in the following be exemplified for K 1ijk . At the local scale, the flow in the x 1 -direction is given by where v 1 = v 1 ⋅ 1 ( 1 being the unit vector in x 1 -direction) is the velocity of the fluid in the x 1 -direction and v 1 is the solution to (12). Note that, in (10), the mass flow has been defined through integration of the velocity over the whole volume instead of along the output surface, hence the 1/l term. Due to the topology of the domain and the chosen boundary conditions, the two formulations are equivalent. The latter formulation is, however, more robust when employed for numerical computation, and it also compares better with other two-scale models. By (10) and by employing the approximation (6), the permeability as defined in (9) can be computed as In order to obtain v 1 , a model governing the flow through the local-scale domain is required. For this purpose, we have chosen to adopt the incompressible form of the model posed in (1), which together with the appropriate boundary is defined as We do note that assuming that the density does not vary is not technically correct, since it is coupled to the pressure, which varies linearly through the domain. Despite that, provided the local-scale domains are much smaller than the global-scale geometry, the pressure drop inside a single domain can be expected to be small compared to its mean value. As a result, the density variation would also be small as compared to its mean value within the domain. Therefore, the error induced by assuming that the density is constant is expected to be small. The great simplification introduced by this extra assumption then justifies taking it. Also, ̄ is defined as the mean free path at the average pressure, p , i.e. (10) The pressure boundary conditions at the inlet i and the outlet o ensure that the v 1 component of the solution v is determined correctly, but note that the approximation (6) is used to determine the values for P i± 1 2 jk . At the sides of the cell that are neither the inlet nor the outlet, l , a no-flow symmetry condition has been specified. A discussion on the validity of this boundary condition is given in Sect. 4.2.2.

Dimensionless Formulation of the Local-Scale Problem
The local problem as posed in (12), and its variants in x 2 -and x 3 -directions, can be used to compute the permeabilities K n directly. However, since the global pressure (at each side of V f ) enters the problem via the boundary conditions, the local-scale flow components obtained from its solution would depend on the global pressure too. Since neither of them are known a priori, the local problems would then have to be solved for a whole range of mean pressures ( p ) and pressure pairs (P i− 1 2 jk , P i+ 1 2 jk ) , requiring a tremendous computational effort. To relax this problem, the independent-and dependent variables in (12) are made dimensionless by the scaling (14), which also include a translation of the pressure solution. Indeed, by introducing the scaling into (12), we obtain where ∇ * indicates that the derivatives are taken with respect to x * . We note that the averaged dimensionless mean free path can be expressed as according to (13) and (14). Thus, although P i− 1 2 jk and P i+ 1 2 jk , which are approximated by means of (6), do not appear explicitly in (15), they affect the solution through ̄ * . Note that ̄ * also contributes to the nonlinearity of the local-scale problem as it includes the localscale variable p * , which represents the dimensionless mean value of the pressure within V f . Having solved (15), the permeability can be computed from the dimensionless form, as K 1 = l 2 K * 1 . The reason for the variation in K 1 with X ijk resides in the definition of ̄ * , and the way in which this dependency will be treated is presented in Sect. 4.2.3.

Comparison with the Volume Averaging Method
As we have stated previously, the method used here is less rigorous than other methods used for constructing two-scale methods such as homogenisation or the volume averaging method. The choice of HMM was motivated by the fact that some assumptions of those methods are incompatible with the construction of the present model. For example, demands such as periodicity of the local-scale domain (needed for homogenisation) or that the typical pore size is much smaller than the local-scale domain (needed for the volume averaging method) are hard to reconcile with the large variation in permeability between neighbouring local-scale domains assumed in this work. As discussed in this section, however, the resulting local-scale problem is very similar to those obtained with the more rigorous methods. To show that, we compare the local-scale problem (15) and the problem (3.7) in Lasseux et al. (2016), which they obtained by using the volume averaging method.
Note first that our problem (15) can be seen as a projection of theirs in x * 1 -direction, so that only K * 1 , and not the whole permeability tensor, is considered. This goes for the problems defined to obtain K * 2 and K * 3 too. From this point of view, the equations are equal if one interprets their D as v * and their d as p * − x * 1 . (a more explicit comparison for a similar case can be found in Pérez-Ràfols et al. (2016)). Note also that (3.7e) in Lasseux et al. (2016) is equivalent to (17) and that the periodicity condition (3.7d) is replaced by the combined symmetry and Dirichlet boundary conditions imposed here. The main differences between our formulation and the one in Lasseux et al. (2016) are, therefore, that we omit the off-diagonal terms of the permeability tensor and that the faces other than the inlet and outlet have no-flow condition instead of being periodic boundaries as in Lasseux et al. (2016). We shall now see, however, that these differences are not as notable as they seem.
In the present analysis, the local-scale domains, acquired from the X-ray microtomography scan, are not periodic, as clearly seen in Fig. 3a. If this domain is not manipulated, this would impede imposing periodic boundaries in the faces other than the inlet and the outlet. Indeed, the fluid domain would not match at the periodic boundary. Due to the small size of the local-scale domain, as compared to the size of the balls, this is hard to solve unless the local-scale domain is mirrored. Such a mirrored domain was studied in Almqvist et al. (2011), where the authors compared local-scale problems including no-flow conditions with the periodic local-scale problem obtained by homogenisation. They could show that, if the local-scale domain was mirrored, the two problems were equivalent. Moreover, the off-diagonal terms in their "flow-factor matrix" vanish, provided that the local-scale domain has been obtained by mirroring. This directly translates to the present three-dimensional problem, and it means that the solutions to (15) and (3.7) in Lasseux et al. (2016) are equivalent, provided that the local-scale domain in the periodic problem is constructed by mirroring. Moreover, we note that, because of the rapid spatial change in permeability, flow perpendicular to the global-scale pressure drop will arise in all global-scale simulations, thus reducing the possible error coming from the zero off-diagonal terms. Also, as shown in Almqvist  Almqvist et al. (2007Almqvist et al. ( , 2008 by means of another approach related to homogenisation, the off-diagonal terms do not significantly influence the global-scale solution, if the local-scale topologies do not exhibit distinct directionality. Actually, topologies with a directionality aligned with or perpendicular to the flow direction are intrinsically symmetric; thus, the off-diagonal terms are identical to zero in these cases.

Expansion of the Local-Scale Problem in *
In the previous section, we have shown that the permeability's dependence on global-scale conditions resides in ̄ * defined in (16). More precisely, these global-scale conditions are the global-scale pressure P and the temperature T. Since we are considering ideal gas flow under iso-thermal conditions, the temperature is just an input, the viscosity is constant and density is given by the pressure, following the constitutive relationship (2). This effectively means that the global-scale pressure is the only variable that couples the scales. Next, we will present a methodology for how to fully decouple the scales. We start by noting that, in (15), ̄ * appears in combination with , i.e. as the product ̄ * . In order to estimate the magnitude of this product, one should recall that, in this work, only low Knudsen number flows are considered, i.e. 10 −3 ≲ K n ≲ 10 −1 . The Knudsen number can be expressed as where d p , in this case, characterises the size of the voids between the balls. Now since ̄ * =̄∕l , where l is chosen such that l ≳ d p , it is clear that ̄ * ≲ K n . Since is of order unity, we have that ̄ * ≲ 10 −1 . Therefore, ̄ * is small enough to allow for expansions of the dependent variables v * 1 and p * in terms of it and is also what they did in Lasseux et al. (2016). Indeed, where R v m and R p m are the residuals of the expansions up to order m. We really want to stress here that the expansion considered here intends to capture only the effect of ̄ * on the topology. Obviously, as the original equations are of first order in K n , we do not claim that the model provides for an approximation of higher order for the behaviour of the gas.
Introducing (19) into (15) and requiring that the system must be valid for each order of ̄ * independently, we find for J = 0 and, for J ≥ 1 These problems are similar to (15), but as they do not include ̄ * , they are independent of the global-scale pressure, making the global and the local scale fully decoupled. Thus, having solved the virtual problems (20) and (21), the auxiliary dimensionless permeabilities K * I J , for each ijk , can then be computed as Thereafter, an approximation of the permeability K 1 can be obtained by as K I = l 2 K * I . Note that K I (0) = l 2 K * I 0 represents the permeability without slip and the terms K I (J) = l 2 ̄ * J K * I J , for J = 1, 2, … , m , corrects for slip flow. We can, therefore, define rarefaction correction factors S I (J) = K I (J) ∕K I (0) and write the permeability K I as The advantage of expressing the effect of rarefaction in terms S I (J) instead of K I (J) is that the correlation between K I (J) and the former is less strong than with the latter. It is also on this basis that we neglect the correlation completely, while constructing the resulting stochastic model, which then becomes much simpler.

A Realisation of a Local-Scale Solution
Before presenting the stochastic two-scale model, we present an example of a local-scale result, with topology and solutions shown in Fig. 3. The topology of the permeable phase of the local-scale domain depicted in grey in Fig. 3a has been randomly sampled from the original measurement shown in Fig. 1 and its size is comparable to the radius of the balls. In Fig. 3b-e, the pressure p * 0 and velocity v * 0 of the virtual problem (20) and p * J and v * J of (21), for J = 1 , 2 and 3, are depicted. These were obtained using COMSOL Multiphysics. A maximum element size of 0.033 has been used on the boundaries of the dimensionless local-scale domains, while the elements within the bulk were restrained to a maximum size of 0.13. The maximum growth rate in the transition between bulk and boundary was chosen as 1.18. This has led, in this particular case, to around 4 × 10 5 elements, although some of the domains were resolved with the double amount of elements. The reason for the smaller elements on the boundary surfaces and for the (relatively) small growth rate is to capture a more accurate representation of the velocity gradient close to the boundaries of the local-scale domain. Note that this is necessary in order to accurately capture the boundary condition (21c), which is the driving force of the flow. It has been shown to be quite challenging to accurately resolve this condition for several reasons, that is, (i) due to the fine mesh needed, (ii) the effect of measurement errors and irregularities on the surface (note that the curvature is particularly affected by those) and (iii) the cumulative nature of the error as J increases. One can, in fact, identify some oscillations in the resulting pressures when inspecting closely Fig. 3, particularly in Fig. 3d.
Focusing now on the flow pattern, the case J = 0 depicted in Fig. 3b is the most intuitive as it has pressure-drop type of boundary conditions and as it reflects flow without rarefaction. It thus makes sense that the flow occurs from left to right, with the streamlines getting closer to each other when moving from the larger inlet to the smaller outlet. The others are harder to interpret as the driving force behind the flow is given by the boundary condition (21c). In Fig. 3d, for instance, one can clearly see that the main portion of the flow occurs along the boundary. One should notice that the absolute value of p * J increases as J increases since, in general, | which then translates into an increased pressure. The same occurs with the magnitude of the velocity. In the non-scaled version of the problem, this increase can be shown to be of order 1/l. The scaling, not being able to capture the average velocity, cannot either fully erase this increase. The relevance of this is that, in order for any convergence to be expected in (23), we must have ̄ * < 1 . This, however, should not be a problem in general, since we need 10 −3 ≲ K n ≲ 10 −1 for the localscale problem to be valid. Recall that ̄ * differs from the Knudsen number by a factor d p ∕l , where d p is the characteristic size of the voids between the balls. In our case, d p ∕l ≃ 10 −1 , and the difference is a factor of ten, but in other arrangements or porous media, this is a concern that one should be aware of.

The Stochastic Two-Scale Model
The model described in Sect. 4 allows for computing the flow in a two-scale manner. Following it literally, however, would imply computing the permeability for each of the localscale domains V f , defined by the elements of the global-scale grid. This is, obviously, not practical from a computational point of view, but most importantly, it is not necessary. This can be explained as follows. Despite that the permeability can vary greatly between each local-scale domain, the values will follow a certain probability distribution, governed by the structure of the porous medium. The approach adopted for the stochastic two-scale model here is, therefore, to compute the permeability for as many local-scale domains as it takes to obtain a representative probability density function (PDF), which describes the variability with respect to the topology. Subsequently, the PDF is used to randomly generate a permeability distribution occupying each element of the global-scale domain, which in a statistical sense represents the original one. Note that this is actually more informative than a deterministic solution, since it describes how the flow through a whole set of statistically equivalent realisations varies and not just how the flow through one particular realisation varies. These ideas are elaborated upon within this section, which first describes how the local-scale permeability variability can be characterised and how the associated PDF's can be obtained. Thereafter, the recipe for constructing the global-scale model is presented.

Characterisation of Local-Scale Variability
As stated previously, the permeability is regarded as a random variable given by a PDF, which describes its variability with respect to the local-scale topology. We remark that the size (l) of scaled domain, obtained by scaling the original measurement depicted in Fig. 1, is of the same order as the diameter of the balls ( d p ). This implies that there is a rather high possibility that the resulting permeabilities exhibit a strong correlation with each other. There are, however, cases where the correlation is expected to be small. For example, if the centre of one of the balls is located close to the centre of one of the boundaries, then the permeability in the direction normal to the boundary surface would be minimal, while the permeabilities in the other direction may be rather large in comparison.
In the present work, nine local-scale domains (such as the one exemplified in Fig. 3) were extracted from the original domain in Fig. 1 and then scaled to one of the four different sizes l = 18 , 32, 56 and 100μ m. Thereafter, 27 permeabilities (in all three directions for each of the nine local-scale domains) were computed, as described in (24). Because of not considering the correlation, a single PDF for the permeability distributions in all three  Fig. 4. It is clear that the spread is quite large (notice that the only one with linear x-axis is S I (2) ), which, of course, is to be expected due to the small size of the domain. In order to utilise these results to generate permeabilities and distribute them over the global-scale mesh, we have chosen to fit normal (Gaussian) distributions whenever possible. In cases where normality could not be assumed, i.e. for K * I (0) , S I (1) and S I (3) , log-normal distributions were used instead. This choice is motivated by the fact that these variables can only take positive values, and the criteria for not using normal distributions were based on the Lilliefors' 1967 statistical test. As indicated in Fig. 4, the uncertainty in the coefficients of the fitted distributions is now much smaller and thus can be used in the global-scale simulation. The log-normal distribution and other options as well have been used to model the variability of the permeability of porous media previously; see e.g. (Selvadurai et al. 2014;Selvadurai and Selvadurai 2010) for the former and Ricciardi et al. (2005) which employed the beta distribution. One advantage with the log-normal distribution is that the unintentional generation of negative values of permeability is avoided. Albeit being out of the scope of this paper, an interesting, and more complex, approach would be to follow the one in Guilleminot et al. (2012), where the maximum entropy principle is used to find the distribution that introduces the smallest bias. At last, it should be noted that for a flow close to the percolation limit, many of the local-scale domains would exhibit zero permeability. Such a distortion does, in turn, mean that it will no longer be possible to describe it using a log-normal distribution.

Construction of the Global-Scale Model
At this point, the apparent permeability K I (0) , as well as the correction functions S I (J) , has been described by log-normal and normal distributions. We will now proceed with the construction of the global-scale model, for which a mesh, with nodes at each ijk , of the whole domain of study is first defined. Then, values of K I (0) and S I (J) are randomly sampled from the corresponding distribution and, thereafter, assigned to each element V f (centre points at ijk ) of the mesh. With the values of K I (0) and S I (J) generated, the permeabilities K 1 , K 2 and K 3 can be obtained and (8) can be solved for the global-scale pressure distribution P ijk . Once the pressure distribution has been obtained, the total flow can simply be computed by summing up the flows of the cells located at the outlet of the global-scale domain. Computed in this way, the resulting total flow is a direct result of one realisation of a randomly generated permeability map. Realisations of randomly generated permeability maps and subsequent computations of the total flow must and are, therefore, repeated until a converged mean value and standard deviation are obtained.

The Interplay Between Rarefaction and Variability
In order to study the interplay between rarefaction and the variability of the local-scale permeability, with respect to the local-scale topology, numerous global-scale simulations were performed. The global-scale domain was discretised using a mesh consisting of 100 cells in the x 1 -direction and 10 cells in both the x 2 -and the x 3 -directions. To drive the flow, a pressure drop ( p i = 10 5 Pa and p o = 10 4 Pa) was applied in the x 1 -direction and periodic boundaries are applied in the x 2 -and x 3 -directions. This configuration, with the cell aspect ratio aligned with the direction of the pressure gradient, ensures that the pressure drop over each cell would lead to density variations small enough for the local-scale flow to be considered as incompressible and modelled by (12). The size of the local-scale domain, l, was used to vary the topology and four different sizes were considered, i.e. l = 18 , 32, 56 and 100 μ m. In this way, the size of the balls and the voids between are also scaled. These choices of l were taken to vary the Knudsen number (18), while keeping it within the range of applicability of the presented model. The fluid is assumed to be nitrogen and in this work modelled as an ideal gas = pM∕(RT) , with M = 28 g/mol , R = 8.314 J/molK) , viscosity ( = 10μPas ) and the temperature was specified as T = 273 K.
Before presenting the outcome of the study, let us define what we hereafter will refer to as the rarefaction effect R and the variability effect V , which will be used in order to interpret the results. To this end, we have chosen the following definitions: and where Q r v denotes the flow including both variability of the permeability and gas rarefaction, Q 0 v with variability but without rarefaction. Similarly, Q r 0 denotes the total flow without considering the variability but including rarefaction and Q 0 0 when both variability and rarefaction are excluded. Indeed, when rarefaction is not considered, the Knudsen number is zero (so is ̄ * ) and the total flow Q 0 v is determined by only K * -computed by (24), and not the correction functions S I (J) . When including the variability of the permeability, it is included via the PDF according to Fig. 4, only. When variability is not considered, the flow ( Q r 0 ) is determined by the permeabilities K I , which are computed by (24), assuming that the permeability K I (0) and the correction functions S I (J) are constant and equal to the mean values; log K I (0) = −6.137 , log S I (1) = 2.746 , S I (2) = 22.651 and log S I (3) = 11.395 of the corresponding probability density function.
The effects of gas rarefaction and the effect of variability of the local-scale permeability with respect to the local-scale topology, with varying domain size (l), are visualised in Fig. 5.
As expected, V 0 is constant, since, in the absence of rarefaction, scaling the geometry only changes the flow by a multiplicative factor. Therefore, Q 0 v and Q 0 0 are changed by the same factor and their ratio remains unchanged. Of course, the variability still does have a significant effect and V 0 is not unity. When there is rarefaction, however, the variability effect V r increases as l increases and that the influence of gas rarefaction decreases at the same time. This is also confirmed by the results depicted in Fig. 5 (right), which also shows that including the variability reduces the rarefaction effect. These results are consistent with the ones presented by Germanou et al. (2018) when studying a two-dimensional random pack of spheres with a deterministic model. Let us now elaborate on the possible explanation for them. To this end, we start by revisiting the work by Pérez-Ràfols et al. (2016), where it was shown that the total flow will be significantly lower in the case when there is no gas rarefaction. Indeed, due to the variability in the local-scale permeability, the fluid does not follow a direct path between the inlet and the outlet. Instead, it meanders through the local-scale cells, looking for the path offering the least resistance. This does, however, increase the effective tortuosity, making the course of flow longer. In turn, this lowers the pressure drops and reduces the total flow. Moreover, the probability of passing through a low-permeability region on this meandering path is increased, and this reduces the total flow. In practice, this acts as a flow restriction, increasing the pressure upstream. When the flowing gas is rarefied, the presence of these restrictions will cause a cascade effect, as the increase in pressure upstream leads to a decrease in ̄ * which reduces the effective permeability there (24). In turn, this will cause a further decrease on flow and a further increase in pressure (and thus permeability) further upstream. The relevance of this effect is shown in Fig. 6, in terms of the global-scale pressure and the averaged Fig. 5 The relationship between flow and the size of the local-scale domain including the variability of the permeability with respect to the local-scale topology, both with and without gas rarefaction (left) and including gas rarefaction both with and without including the variability of the permeability with respect to the local-scale topology (right) Fig. 6 Non-dimensional global pressure ( P∕p o ) over the local cells as a function of normalised length, from inlet (0) to outlet (1) (left) and the corresponding increase in ̄ * , which is related to the Knudsen number K n (right). Results obtained without including the effect of rarefaction are depicted by the continuous black lines, and those that do are depicted in colour, including error bars for the variability. The results are presented for the four local-scale cell sizes l = 18 , 32, 56 and 100 μm and the legend is the same for both the figures dimensionless mean free path ̄ * (which is directly proportional to the Knudsen number K n ). Clearly, the global-scale pressure becomes larger (left) at locations far from the outlet and ̄ * becomes smaller (right). It should be noticed, however, that this reduction likely vanishes for a flow situation close to the percolation threshold. In such case, it is the critical constriction that controls the flow and what happens far away from there is rather irrelevant. Thus, although still present, the cascade effect would only very slightly influence the total flow.
One might note that through the discussion of the results, there has been no need to invoke the precise nature of the problem. Indeed, we have only used the fact that permeability increases due to rarefaction but not the precise way in which this occurs. This suggests that the reduction in the rarefaction effect induces in the total flow-due to the variability of the local-scale permeability-can be expected to be relevant in a more general case, not necessarily involving slip flow nor connected to the particular model used here. Of course, the strength of the reduction due to the rarefaction effect can vary notably in other cases, and this requires more detailed investigations. A critical implication of the result presented here is that, in general, a universal relation between the rarefaction effect and the intrinsic permeability, as presented in Bravo (2007) and Kalarakis et al. (2012), cannot be established. Indeed, rarefaction would affect the flow through two different porous media with the same intrinsic permeability differently, depending on the characteristics of their localscale topology. Of course, if the porous medium studied is close enough in terms of their spatial variation, as in Bravo (2007) and Kalarakis et al. (2012), such general relations can be found even for fairly different topologies.
We note that the heterogeneous multi-scale method is feasible and that it may be used when constructing two-scale models. This implies that the same procedure followed here can be applied to models with second-order corrections  or more complex models describing the transition regime, e.g. Moyne and Christian (2003) and Charrier and Dubroca (2004). Although the computational effort to run sufficiently many localscale cases would still be high, it would be markedly lower compared to if a domain large enough to obtain converged results would be considered. A difficulty may arise in case when the local-scale permeability depends on both the pressure drop and the mean pressure in a more complex way. In such a case, covering a sufficient range of conditions may be very computationally expensive. The expenses could, however, be significantly reduced, by realising only the conditions that are strictly necessary. An efficient design of experiments, which takes into consideration the variation of the permeability with varying conditions and using knowledge-based interpolation schemes, would be required for this. An example of such approach can be found in the work by de Boer et al. (2016), albeit applied to a problem originating in another field of research.

Concluding Remarks
The modelling framework developed and presented here was applied to study fluid flow through porous media confined within the slip-flow regime, although it could be adapted to study other cases as well. This framework is based on a two-scale approach in which the highly heterogeneous nature of the porous medium is modelled by treating the permeability of local-scale domains as a random variable. By utilising the present model, it was shown that the rarefaction effect is significantly reduced because of the variability of permeability with respect to the local-scale topology. This underpins the importance of considering the characteristics of the local-scale topology and that results obtained by studying the flow through topologies exhibiting regular patterns cannot be used to draw useful conclusions pertaining to the effects of rarefied flows through highly heterogeneous porous media.
Funding Open access funding provided by Lulea University of Technology.

Compliance with Ethical Standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.