A stochastic two-scale model for pressure-driven flow between rough surfaces

Seal surface topography typically consists of global-scale geometric features as well as local-scale roughness details and homogenization-based approaches are, therefore, readily applied. These provide for resolving the global scale (large domain) with a relatively coarse mesh, while resolving the local scale (small domain) in high detail. As the total flow decreases, however, the flow pattern becomes tortuous and this requires a larger local-scale domain to obtain a converged solution. Therefore, a classical homogenization-based approach might not be feasible for simulation of very small flows. In order to study small flows, a model allowing feasibly-sized local domains, for really small flow rates, is developed. Realization was made possible by coupling the two scales with a stochastic element. Results from numerical experiments, show that the present model is in better agreement with the direct deterministic one than the conventional homogenization type of model, both quantitatively in terms of flow rate and qualitatively in reflecting the flow pattern.


Introduction
Pressure-driven flow through preamble structures such as the percolation of fluids in static seals, see e.g. [1][2][3] and the flow through fractured porous media, in e.g. [4,5], is a problem that has been addressed in many works before. The numerical solution of this problem is, however, complicated due to its multi-scale prevents the usage of a deterministic solution to such problem. Therefore, a two-scale approach is used.
In order to compute the flow rate and the flow pattern, it is first necessary to compute the deformed shape of the gap between the two surfaces (contact mechanics problem). The deformed gap is due to fluid pressure and asperity contact. In this case, we consider the fluid pressure contribution to be small and it is therefore neglected. This permits separating the problem into two smaller problems that can be solved sequentially: compute the deformed gap and compute the flow rate through that gap.
In the following, the method used in this work is presented. First, an introduction to the HMM framework, used to develop the model, is given. After that the two-scale stochastic model is introduced. This includes the theory behind the two-scale formulation of the flow model and the contact mechanics problem and the introduction of the stochastic element. Also, an overview of the solution procedure is given. Finally, the two-scale flow formulation presented in this work is compared against the well-established homogenization technique.

(a) Heterogeneous multiscale method
The HMM is a general framework used to build two-scale models [27] which is flexible in the sense that global and local-scale models do not need to be of the same nature. This will allow us to introduce the stochastic element in a more natural way than the more rigid multi-scale homogenization technique, this is, without the imposition of periodicity in roughness. Under the HMM framework, the model is constructed into steps: (i) definition of a global-scale model, which will include some parameter or variable containing the information of the local-scale model and (ii) definition of a local-scale model that permits finding those variables or parameters. The first step is crucial as a wrong definition of the global-scale method can lead to a wrong specification of the local-scale constraints. Although it can be derived from the local-scale data, it is better to use previous (experimental or theoretical) knowledge about the global-scale behaviour. The localscale model is usually easier to identify, as one can use a more fundamental representation of the reality. However, one must be careful on the selection of boundary conditions for the local-scale problem and the data processing used to transfer information to the global scale. In particular, the boundary conditions in the local-scale model must be specified so that the solution is consistent with the global-scale model. One should understand by consistent that the local scale should be constructed so that the global scale can be seen as a coarse representation of it.

(b) The two-scale stochastic model
In this section, the two-scale stochastic model developed is presented. We start by introducing the two-scale methods for both the flow model and the contact mechanics problem. After this, the stochastic element is introduced in the two-scale formulation. Finally, the solution procedure is outlined.

(i) The two-scale flow model
Before presenting the two-scale formulation, we introduce the deterministic flow model. In this work, we consider the rectangular domain covering the full seal area. In this domain, the gap between the two rough surfaces, h (x), is resolved at a high level of detail. The subscript indicates that the gap oscillates rapidly due to the roughness, thus imposing a high resolution in the representation of the domain. In order to solve the fluid pressure distribution through the gap, the lubrication approximation is used.
where p is the fluid pressure, which will oscillate rapidly in accordance with h . The problem is completed by posing the boundary conditions. In order to obtain a pressure-driven flow, a pressure drop, p, is enforced in the x 1 -direction by setting Dirichlet boundary conditions. In the x 2 -direction, periodic boundary conditions are imposed. Once the fluid pressure distribution is known, the total leak rate can be computed by integrating the mass flow over the domain Ω s . In accordance with the Reynolds equation, this is, were η is the viscosity of the fluid. As stated previously, it often becomes too computationally expensive to resolve the roughness at the global scale. Provided that the period of the oscillations is sufficiently small, one can separate the problem into two scales: one considering the detail of the oscillations (local scale) and the other accounting for the slow variations in the whole domain (global scale). Following the HMM framework, we start by presenting the global-scale model. Then we will present the local-scale one.
The global scale is solved in a domain Ω c , which is a coarse representation of Ω s . As the Reynolds equation is a mass-conservation law, it is natural to define the global-scale model also as a mass conservation law. Following a finite volume formulation, this can be posed as where j 0 are the global-scale fluxes and (i, j) identify any point of the global scale. For each grid point, the fluxes must, therefore, be estimated from the solutions of the local-scale model. We will show that, in agreement with Darcy's law, the flux is proportional to the local pressure drop, δp, and we will refer to these constants as permeability, K. Therefore, (2.4) can be rewritten as and Here the first superscript indicates the direction of the corresponding flux and the second one the direction of the local pressure drop. Note that, in order to express the flux in m 3 s −1 and the permeability in m 3 , (2.5) is scaled by the (constant) viscosity. Together with the boundary conditions, posed equally to those in the deterministic problem, (2.5) permits computing the global-scale fluid pressure. The total leak rate is then where ∂ o Ω c is the part of the boundary of Ω c regarded as an outlet.
In the global-scale model, the permeabilities are the parameters containing the local-scale information. We, therefore, define the local scale in order to compute them. Let us start by deriving the fluxes in (2.4). As an example, consider j 0 i+1/2,j ; the others can be defined analogously. to do so, we define a local domain between the neighbour global-scale points, which has a size x 1 × x 2 . Similarly as in the deterministic formulation, we solve the fluid pressure distribution by means of the Reynolds equation (2.2). We, however, need to pose the boundary conditions in a way that the local-scale model is consistent with the global-scale one. In order to see which conditions must be fulfilled, we assume that the global-scale component of p , p 0 , varies linearly between the global-scale points, i.e.
where the capital letters refer to the discrete representation in the global scale. In order for it to be consistent, the local-scale solution should be in agreement with the above representation. Therefore, it should satisfy and · indicates average over the local domain. Other conditions over the mean value of p could be introduced, but the flux would not be affected. Several boundary conditions can be posed so that the constraint (2.9) is satisfied. Among them, the one providing better results is the near periodic condition, i.e. periodicity of p − p 0 [28]. The boundary condition selected, however, introduces an unwanted coupling between the global and the local scale. In order to remove this coupling, we start by defining p = p 1 + p 2 , with p 1 , near periodic in the x 1 -direction and periodic in the x 2 -direction, and p 2 , defined analogously. Furthermore, we define the non-dimensional variables x 1 , x 2 , Therefore, (2.2) can be written in ω as We can now find a solution by setting the two terms equal to zero independently and thus obtain two problems, both of which are independent from the global scale, i.e.
∇ · (h (y 1 ) 3 ∇p 1 (y 1 )) = 0, p 1 (y) − y 1 1 periodic, (2.11a) and ∇ · (h (y 2 ) 3 ∇p 2 (y 2 )) = 0, p 2 (y) − y 2 2 periodic. It can be seen that the solution obtained by combining this two problems satisfies the condition (2.9). Once the local problems are solved, we can compute the local flux, i.e. We can now identify the local-and the global-scale contribution to the flux. By defining the permeability, K as the local-scale contribution, we can write (ii) The two-scale contact mechanics model The fluid flows through the gaps, h , between the contacting surfaces. This gap depends on the applied load. The flow defined previously is solved on a clearance between the surfaces, h , which varies depending on the applied load. Therefore, the deformed shape of the gap must be computed before the flow problem can be assessed. In order to do so, a model based on the one presented in [15], is used. The equations governing this method can be posed as where the clearance, h , is defined as the original (no pressure) gap, h 0 , plus the deformation, u, caused due to the contact pressure, p c , and a rigid body movement, g 00 . The deformation, u, is computed by the convolution of a coefficients kernel, k, and the contact pressure, p c . This can be seen as adding the deformation caused by the pressure at all points. The equivalent Young modulus, E * , is defined as where E i and ν i are the Young modulus and Poisson ratio of each surface. Equation (2.14c) establishes a complementarity relation between h and p c , meaning that when there is contact, the clearance is zero and when there is no contact, the pressure is zero. Finally (2.14d) imposes that the average contact pressure is equal to a given total load, W. A perfect plasticity condition is defined by imposing p c ≤ H, where H is the hardness of the softer material. This problem can be solved by using the algorithm presented in [15], modified to account for the non-periodicity in the x 1 -direction (e.g. [29]).
In the same way as in the fluid computation, the computational time for the contact mechanics problem rises as the domain size grows. Therefore, a two-scale solution is also desired. In order to justify the two-scale formulation of this problem, we need to justify that (i) it is a correct approximation to use a coarser representation of the surface to capture the global-scale trends in the contact mechanics results and (ii) that it is a correct approximation to compute the contact mechanics problem in the local cells by using the nominal pressure p nom = W local as the only information from the global scale. We start by justifying the point (i). The goal is to show that the utilization of a coarse representation of h 0 gives as a result a coarse representation of h and p c . To do so, we apply a Gaussian filter, g, to (2.14a)-(2.14d). We definê (2.16) as the Fourier transform of a function f (x). The equations then read where the capital letters identify the filtered functions. In order to obtain (2.17c), the property of the filter to maintain the mean value has been used. Equation (2.14c) has not been added because it does not hold exactly for the coarsened problem. If one think on the complementarity condition and how the filter works, one can see that it holds approximately provided that the coarsening does not alter the geometry at the global scale. Regarding the plastic limit, the pressure is expected to be lower, and in a more spread area due to the filtering. Therefore, there should be less plastic deformation. This, however, should not significantly affect the result. Unlike the flow case, where each local cell corresponds to one point in the global scale, this is not the case in this coarse representation. In order to obtain the nominal pressure, p nom , of a local cell, the global-scale result should be averaged over a region of equivalent size.
We now discuss about the correctness of the second statement (ii). It has been shown that the elastic deformation caused by pressure is a long-range effect and it cannot be neglected [30]. However, the requirement here is smaller. The deformation on the local domain needs to be computed only up to a constant (g 00 ). Therefore, by assuming that the elastic deformation caused by pressures far-away from the local domain has a constant value over the considered domain, the local scale contact mechanics problem can be solved independently from the global domain.
In order to assess the validity of this assumption one should note that, according to the St. Vennant's principle, the contribution of the distant regions is smooth. Then, by assuming that the local cell is reasonably away from the boundaries, similar deformation can be assumed in all directions, leading to a flat overall deformation. Although crude, it is an adequate first approximation. The assumption of flat deformation is clearly not adequate if one refers to the deformation on points near the cell boundaries caused by pressures in points neighbouring the local cell. In such close points, however, it is a reasonable claim that the roughness is similar to that on the cell. Therefore, a periodic roughness assumption is also reasonable.

(iii) The stochastic approach
In the previous sections, a two-scale model has been developed. The stochastic element has not, however, been introduced yet. Indeed, following the presented procedure, four cell problems (two for each direction) must be solved for each grid node in the global scale. Moreover, the results can be sensitive to the measurement used for the computations. Because the surface topography is a random process, two different measurements of the same surface (or two equivalent ones) are expected to give different results. It is, therefore, of greater interest to obtain the results in the form of an expected value plus a confidence interval rather than simply a given value.
If one now takes a look at the permeabilities computed one notices that they can be modelled as a random variable following a log-normal distribution. This distribution is, however, rather broad. Moreover, it is expected to see some spatial correlation between the different permeabilities, caused by fluctuations in the global scale. Therefore, a reference parameter, either average interfacial separation,h, or nominal pressure, p nom , is taken from the global-scale contact mechanics computation and the permeabilities are fitted to different log-normal distributions as a function of these parameters. These distributions have a narrower permeability range and, therefore, are more meaningful.
Following this approach, the global scale is therefore not constructed by solving the local problem on each coarse grid point but by randomly generating a permeability value for each of those points. In order to ensure sufficient accuracy, numerous realizations of the global scale are computed following a Monte Carlo approach. The output is, therefore, given as a probability distribution instead of a single value.
It is worth noticing that modelling a gap between two rough surfaces by randomly assign permeabilities from a log-normal distribution (or other similar distributions) is a common practice in the porous media literature (e.g. [20,31]). With the present approach, the parameters for such distribution are no longer obtained experimentally but from modelling the problem at the local scale.
(iv) Solution procedure A flow chart of the solution procedure for the model is depicted in figure 1. In the following, a more detailed description is given. Both scales are effectively separated during the computations and therefore they are treated separately here also. For clarity, the different domains considered are represented in figure 2 We start by considering the local-scale solution procedure. In order to solve the local problems, the whole measured domain is divided into a set of local cells of size x 1 × x 2 . For the following numerical analysis, these cells are mirrored in order to obtain a periodic roughness. We note that this implies that the permeabilities K 12 and K 21 will be zero (see [32] for a justification). However, even without mirroring, these are expected to be small in most cases, also in the surface textures used in this work (figure 3; this should be reconsidered the model is applied to surfaces with different textures). The solution procedure is then as follows: (i) The deformed gap in the local domain is computed for a range of nominal pressures p nom , as described in §2b(ii). It is advised to compute the deformation for a very small nominal pressure, which shall serve as a zero value during the interpolation. (ii) The fluid pressure distribution is computed by (2.11) and the permeability is computed according to (2.12) and (2.13). This is done for the range of nominal pressures previously computed and for a specified range of average interfacial separationsh. A combination of Dirichlet boundary condition in the direction of the pressure drop and homogeneous Neumann boundary conditions in the transverse direction is used instead of the near periodic condition. The reason for this is that this combination is easier to implement and yet equivalent (see [32] for details). (iii) The distribution of permeabilities for the given set of cells, and for each reference parameter (p nom orh) are fitted to a log-normal distribution. The set of cells should be large enough in order to obtain a good estimate for the parameters of the log-normal distribution. A random selection of cells is desired to increase robustness against any possible spacial variation of permeability.
Once the permeabilities are obtained from the local-scale, the global-scale model is implemented in the following steps: (i) A measurement of the topography is used as the input. This measurement domain, Ω, is filtered by using a Gaussian filter and re-sampled to a coarser grid. filter and coarsen measurement domain solve contact mechanics problem by (14) compute deformed gap by (14) compute fluid pressure distribution by (11) and permeability K by (12)  (v) Permeabilities are assigned to each point in the global scale. To do so, the parameters of a log-normal distribution are first assigned to each point by interpolation based on the reference parameter (p nom orh). Then a permeability value is generated using that distribution. (vi) The global scale fluid pressure distribution is computed by means of (2.5) and the total leak rate is computed by means of (2.6). (vii) Points 5 and 6 are repeated until good estimates for the leak rate and its variance are achieved. If the leak rate does not follow a normal distribution, the 95% CI can be computed by performing a large-enough number of realizations so that 95% of them have a leak rate which consistently falls inside a fixed range.
It is important to note that there will be cells with zero local permeability, occurring when there is no available path connecting the inlet and the outlet in the local domain. Those cells cannot be included in the log-normal distribution. Instead, the fraction of cells with zero permeability is stored. Then a number of cells corresponding to this fraction is set to zero in the global scale. The reference parameter should be chosen between the nominal pressure and the average interfacial separation as the one that better correlates with the permeability. For example, for the turned surface with the topographies depicted in figure 3, permeabilities in the x 1 -direction correlate well with nominal pressure (except for when contact is lost) while the one in the x 2 -direction correlate better with the average interfacial separation.

(c) An average permeability approach
Most of the most successful two-scale models for fluid flow are based on the homogenization technique. Therefore, we benchmark our flow model to it. We first note that the proposed boundary conditions make the problem equivalent to that of standard homogenization. Indeed, by making the change of variables χ 1 = p 1 + y 1 1 and χ 2 = p 2 + y 2 2 , the local problems read ∇ · (h 3 ∇χ 1 ) = which is equivalent to the homogenized results by reinterpreting the flow factors a 11 and a 21 as the local permeabilities K 11 and K 21 (e.g. [10]). We note that this resemblance was already studied in [32] and that a more rigorous derivation of the consistency between the HMM formulation and the standard homogenization is given in [33] for the more general time-dependent case.
The main difference lies then in the assumptions regarding the local-scale roughness. In homogenization it is assumed to be periodic, therefore, it is natural to assume that the flow factors (or permeability) are constant or vary smoothly in the global scale or, at least, in subdomains of the global scale much larger than the global-scale size. One then obtains the values by averaging several cell realizations. Hereafter, we will refer to this procedure as the averaged permeability approach. By contrast, the periodic boundaries in the given approach represent only a consistent way to couple the two scales. This allows accounting naturally for spacial variation of permeability. We have used this feature to introduce a stochastic representation of the surface in §2b(iii).
In the average permeability approach, the pressure-driven flow problem considered will lead to a constant pressure gradient in the x 1 -direction and a zero pressure gradient in the x 2 -direction. The total leak rate, Q H , is therefore whereK H is the average permeability in the x 1 -direction.

Results
The performance of the model is shown in two steps. First, the two-scale model is validated against a deterministic solution. Then, an example of its usage, the model is applied to a test case. In both cases, two different topographies are used, as shown in figure 3. Topography 1 is highly anisotropic and corresponds to a turned surface. Topography 2 has been sand-blasted to obtain a more isotropic surface (although still with a clear anisotropic component). These topographies have also been chosen in order to obtain two clearly different total flow rates. The leakage studied is the one occurring between these two surfaces and   in the x 1 -direction, in order to impose a pressure drop p. The topography used is mirrored in the x 2 -direction in order to avoid a jump in topography when applying the periodic boundaries.
The material, assumed to be linear elastic-perfectly plastic, is described by the parameters E 1 = E 2 = 206 GPa, ν 1 = ν 2 = 0.3 and H = 2.75 GPa. The leak rate is computed by (2.3), (2.6) or (2.20), depending on the methodology. In the presented results, however, it is scaled by the factor η/ p. The Gaussian filter used for coarsening have a cut-off length 0.01 times the total length of the measurement and the filtered surface is re-sampled to one-sixteenth of the previous number of points. The original lateral resolution for the contact mechanics (both for the deterministic and the local-scale cells) is of 0.896 µm, which corresponds to a grid size of 600 × 1440 nodes for the topographies presented in figure 3. The tolerance for the load is set to 1 × 10 −3 % and the one for the contact plane is set to 10 −10 m (see points 4 and 8 of the algorithm presented in [15] for reference).

(a) Two-scale model validation
In order to do the validation the measured domain Ω is separated into a set of local cells, ω. The size of the local cells in the x 1 -direction, x 1 is taken so that each cell corresponds to one wavelength of the main frequency in Topography 1 (figure 3). In order to facilitate comparison, the same sizes have been also used for Topography 2. This means that the measurement domain is divided into six cells, each of which with a length of 0.21 mm. In order to verify convergence with the domain size, three widths in the x 2 -direction, x 2 , are taken. These are 0.18, 0.09 and 0.045 mm respectively, giving 6, 12 and 24 local-scale cells in that direction. A typical two-scale formulation requires the local domain to be at least one order of magnitude smaller that the global domain [34], which is not fulfilled in the present validation. This is because of the restriction on the measurement domain size. We note, however, that the error obtained when using a larger domain is expected to be smaller. Therefore, the validation can be directly extended to the more common situation where the two-scale separation is clearer.
The three techniques presented in § §2b and c are compared. In all three cases, the contact mechanics model described in §2 is used to compute the deformed gap. The deterministic approach for contact mechanics is used to compute the deterministic flow, while the two-scale contact mechanics approach is used for the other two techniques. The results for a range of total load W and for the three different widths x 2 are depicted in figure 4.
The present two-scale model gives a solution that is in good agreement with the deterministic one for both tested topographies, specially for the two wider local cells. The average permeability approach predicts generally a higher leak rate. Although it could be acceptable for Topography 2 when using the higher width, this is not true for Topography 1. Moreover, for a given cell size, the present two-scale model gives significantly better accuracy. The error between the present two-scale formulation and the deterministic solution, defined as is depicted in figure 5. It can be observed that reasonable values are obtained for the two wider widths (maximum error for Topography 1 using cells of width 0.18 and 0.09 mm is 7 and 11%, respectively, and 5.5 and 19% for Topography 2 and the same widths). In order to identify the main source of error, the error between the deterministic solution and the solution with the present two-scale fluid model using the deterministic computation for the gap as the flow domain is also depicted. For Topography 1, it can be seen that it is much smaller, which identifies the two-scales contact mechanics model as the main source of error for this case. For Topography 2, however, this is not the case. The reason given is that due to the cell selection made for Topography 1, the boundary conditions in the flow model are particularly well suited, as the pressure becomes nearly constant at the edges of the domain due to the larger gap. This is, however, not the case for Topography 2. Thus, it becomes clear that the error and its main source is topography dependent. Despite that, given a careful choice of domain size, a sufficiently good approximation can always be made. The average permeability approach based on the homogenization technique has been proved to produce good results for similar pressure-driven conditions [16]. The discrepancy observed in this work is, therefore, attributed to the small size of the local cells, which makes them not representative of the topography. This is further enhanced when the leak rate is reduced. In order to see why the local-scale cells are not representative when using the average permeability approach but are representative in the present two-scale formulation, one can analyse the flow pattern. The flow pattern obtained using Topography 1, in terms of deterministic flux, is depicted in figure 6a for a high total load W, i.e. for a small leak rate. It can be seen that the fluid follows a particular pattern through the gap. Whenever a path is available, it advances in the x 1 -direction, which is the direction of the pressure gradient. When this is not possible, or when the path is too small, the flow advances in the x 2 -direction until a large path is found and flow in the x 1direction is possible again. Furthermore, the advance in the x 2 -direction is as large as the half of the width of Ω, which is the maximum value due to the surface has been mirrored. In fact, one should expect even larger advances in this direction when computing over a larger domain. This implies that a representative cell for the average permeability approach would need to be at least as large as Ω for it to be possible to capture the flow pattern correctly. In the present two-scale formulation, however, the heterogeneity in the permeability distribution can enforce such flows in the x 2 -direction even when using small local cells. This can be seen in figure 6b, 1.00 where it can be observed that the flow pattern is satisfactorily captured by the present two-scales approach. Another fundamental difference of the two approaches lies in how much large values of permeability affect the total leak rate. In the average permeability approach the cells with high permeability will increase the average permeability and therefore significantly affect the total leak rate. In the present two-scale formulation, the leak rate is controlled by the narrowest constriction that must be crossed. This is, in fact, a more correct representation of the problem [26]. Therefore, the influence of the high permeability values is not as large as in the average permeability approach. This is also the explanation for the observed differences in the leak rate prediction shown in figure 4. In figure 7, a similar comparison, using Topography 2, is depicted. The present two-scale approach captures, again the correct flow pattern. In this case, however, the flow in the x 2 -direction is much less important (although significant enough to affect the results) and the constrictions are more evenly spread. This, in turn, explains why the average permeability approach does not deviate as much as in the previous case.
It is important to emphasize that the average permeability approach only predicts a too high leak rate when it is small, i.e. at high total load. Otherwise, the variation in the permeabilities are not expected to be large and the flow perpendicular to the direction of the pressure drop is expected to be less important. In those cases, the average permeability approach can be considered a good approximation. This can be seen in figure 4, where the prediction from the average  permeability approach is close to the deterministic one at the lowest total load and the deviation observed when using Topography 2, which exhibits a higher leakage, is smaller.

(b) Two-scale stochastic model
Once having validated the two-scale part of the model, we show in this section the performance of the model by applying it to a case example. The measurement domain used for this is the same as in the previous section. Let's start by considering the local-scale results. Having computed the permeabilities for the set of local cells, these values must be fitted to log-normal distributions. In order to decide which reference parameter to use, we study the correlation between the permeability and the two reference parameters (nominal pressure and average interfacial separation). These relations are shown in figure 8. Focusing first on the correlation between average interfacial separation and permeability, one can observe that the permeabilities in the x 2 -direction can be approximated by K ≈h 3 . The same trend is followed by permeabilities in the x 1 -direction when there is no contact and the separation is large. This approximation will be used for simplicity when obtaining the permeability distributions in the global scale. When the nominal pressure increases, the permeability in the x 1 -direction is reduced significantly (note here that a large number of cells have zero permeability at high nominal pressures), while the average interfacial separation remains nearly constant. Therefore, we use the nominal pressure instead of average interfacial separation as the reference parameter. As seen in figure   -15 x 2 -direction (mm) log 10 K 11 log 10 K 11 log 10 K 22 log 10 K 22 correlation between permeability and nominal pressure, but one can clearly observe a trend in both mean value (decreasing) and spread (increasing) of permeability. The percentage of cells with zero permeability also increases with increasing nominal pressure. In order to describe the permeability in the x 1 -direction log-normal distributions are fitted for each value of the nominal pressure. The log-normal distributions relating permeability and nominal pressure and the approximation K ≈h 3 are used to build the permeability distributions and the global scale leak rate is computed. An example of a permeability distribution is shown in figure 9 and the resulting flow pattern for a unitary pressure drop in the x 1 -direction is depicted in figure 10. It can be seen in figure 10a how the flow advances in the x 2 -direction until an easy path is found in the x 1direction. For the case of Topography 2 (figure 10b) the leakage is much higher, resulting in a less important flow in the x 2 -direction. One could define this flow as advancing in the x 1 -direction until a blockage is found, then doing a small shift in the x 2 -direction to avoid it and continue in the x 1 -direction.
One can think on these flow patterns, particularly the one corresponding to Topography 1, similarly to the description given by Persson & Yang [2], where the flow domain is described as a network of paths with critical constrictions randomly located. In the present model, however, no assumption on either the distribution or the size of the paths and the constrictions needs to be imposed.
The leak rate is computed for several permeability distribution realizations. When a sufficient number of realizations have been computed, the leak rate can be described statistically. The probability distribution for the leak rate for three global domain sizes can be seen in figure 11a for Topography 1 and in figure 11c   the leak rate or the global-scale domain become smaller, the variance of the leak rate distribution increases and it becomes skewed. If a normal distribution can be assumed for the leak rate, it can be characterized by using only the mean value and the standard deviation of that distribution.
In the more general case where the distribution is not known, a 95% CI for the leak rate can be computed, as shown in figure 11b and d. Without knowing the distribution, this interval is obtained by performing a large number of computations until 95% of the computed leak rates fall consistently inside a fixed range. Some general trends can be observed with increase of globalscale domain. As it is increased in the x 2 -direction, a reduction of variability of the leak rate is observed. Both the reduction of variability and the tendency to become normal when increasing the size in the x 2 -direction are expected results as the effect of extreme values will be smaller and all realizations will be more similar due to the presence of more local cells. It is also notable that the leak rate per unit width is slightly increased in the case of Topography 1. This is because of the larger domain allow areas with high permeability to appear more often in the realizations. This effect is not observed in Topography 2 because its larger leakage allows for a better averaging even at the smallest domain. An increase of the global domain in the x 1 -direction has the opposite effect and mean leak rate is reduced. The reason for this is that the fluid must cover a longer distance and a narrower constriction is more likely to be encountered. One should expected the leak rate per unit width to converge to a certain value as the size of the domain in the x 1 and x 2 -directions increases. It is noticeable, however, that convergence with size in the x 1 -direction is far from being reached even at realistic seal dimensions.

Concluding remarks
In pressure-driven flows, as the leak rate decreases, the fluid follows a pattern which becomes wide in the direction transverse to the pressure gradient. This obliges utilizing a large local domain in the common two-scale formulations. However, the required size of the local domain can become too large and approach the global domain, which makes it not possible to perform a scale separation of the flow. In order to avoid this problem, a two-scale model based on the HMM framework has been developed. This model does not assume periodic repetition of the topography and, therefore, can capture this wide flow pattern via local variation of permeability. This allows using much smaller local-scale domains. It has been shown that in this way a two-scales model which captures correctly the flow pattern can be developed. Moreover, the presented two-scale formulation permits the construction of the local scale by considering the local permeabilities a random variable. This feature allows inclusion of the inherent uncertainty coming from the surface topography in the model and estimating the uncertainty on the leak rate due to this cause.
In the case study presented, it has been shown that smaller global-scale domains, as well as smaller leak rates, lead to more uncertainty on the predicted leak rate. It has also been shown that, even for relatively large global domains, the leak rate per unit width is dependent on the domain size. This has been explained as an effect of the random construction of the permeability distributions.
Data accessibility. The datasets supporting this article have been uploaded as part of the electronic supplementary material.