# Effective transmissivity for slip flow in a fracture

Tony Zaouter, Francisco J Valdés-Parada, Marc Prat, Didier Lasseux

► To cite this version:

Tony Zaouter, Francisco J Valdés-Parada, Marc Prat, Didier Lasseux. Effective transmissivity for slip flow in a fracture. Journal of Fluid Mechanics, 2023, 969, pp.A9. 10.1017/jfm.2023.564. hal-04181152v2

HAL Id: hal-04181152

https://hal.science/hal-04181152v2

Submitted on 12 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effective transmissivity for slip flow in a fracture

Tony Zaouter, Francisco J Valdés-Parada, Marc Prat, Didier Lasseux

► To cite this version:

Tony Zaouter, Francisco J Valdés-Parada, Marc Prat, Didier Lasseux. Effective transmissivity for slip flow in a fracture. Journal of Fluid Mechanics, 2023, 969, pp.A9. 10.1017/jfm.2023.564. hal-04181152v2

HAL Id: hal-04181152

https://hal.science/hal-04181152v2

Submitted on 12 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effective transmissivity for slip flow in a fracture

Tony Zaouter<fnref n="1" />, Francisco J. Valdés-Parada<fnref n="2" />, Marc Prat<fnref n="3" /> and Didier Lasseux<fnref n="4" />,†

(Received 13 March 2023; revised 7 July 2023; accepted 7 July 2023)

A simple efficient method is presented for the determination of the intrinsic transmissivity tensor, as well as the intrinsic correction tensors at successive orders in the dimensionless slip parameter, that predicts the effective transmissivity tensorial coefficient for steady, one-phase, isothermal, creeping flow of a Newtonian fluid with slip boundary condition in a rough fracture. It is demonstrated that the solution of the first *N* ancillary closure problems provides the slip correction tensors up to the 2*N* − 1 order, hence reducing the computational requirements by a factor of ~2 compared with the conventional approach. In particular, the first-order correction tensor (i.e. a Klinkenberg-like tensor) can be obtained by solving the closure problem required for the computation of the intrinsic transmissivity tensor. In addition, symmetry and definiteness (positiveness or negativeness) properties of the individual tensors are analysed. It is shown that a Padé approximant, built on the correction tensors at the first three orders, outperforms the predictions for the effective transmissivity tensor. The new approach is illustrated and validated with numerical examples on model rough fractures.

**Key words:** lubrication theory, porous media

## 1. Introduction

Flow in fractures is a situation widely encountered in practical applications such as the leakage rate of seals (e.g. Pérez-Ràfols, Larsson & Almqvist 2016; Rui *et al.* 2023), gas leak through reservoir caps or gas recovery in fractured rocks (see Zhou *et al.* (2023) and references therein), among many others. Modelling flow in rough fractures

† Email address for correspondence: didier.lasseux@u-bordeaux.fr

<footnote n="1">¹CEA, DES, ISEC, DPME, SEME, Laboratoire d'Étanchéité, Université de Montpellier, 30207 Marcoule, France</footnote>

<footnote n="2">²División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, Av. Ferrocarril San Rafael Atlixco 186, Col. Leyes de Reforma, 09310, CDMX, Mexico</footnote>

<footnote n="3">³Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, 31400 Toulouse, France</footnote>

<footnote n="4">⁴Université de Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400, Talence, France</footnote>
<figure><img src="image_2.png" /></figure>

# Effective transmissivity for slip flow in a fracture

Tony Zaouter<fnref n="1" />, Francisco J. Valdés-Parada<fnref n="2" />, Marc Prat<fnref n="3" /> and Didier Lasseux<fnref n="4" />,†

(Received 13 March 2023; revised 7 July 2023; accepted 7 July 2023)

A simple efficient method is presented for the determination of the intrinsic transmissivity tensor, as well as the intrinsic correction tensors at successive orders in the dimensionless slip parameter, that predicts the effective transmissivity tensorial coefficient for steady, one-phase, isothermal, creeping flow of a Newtonian fluid with slip boundary condition in a rough fracture. It is demonstrated that the solution of the first *N* ancillary closure problems provides the slip correction tensors up to the 2*N* − 1 order, hence reducing the computational requirements by a factor of ~2 compared with the conventional approach. In particular, the first-order correction tensor (i.e. a Klinkenberg-like tensor) can be obtained by solving the closure problem required for the computation of the intrinsic transmissivity tensor. In addition, symmetry and definiteness (positiveness or negativeness) properties of the individual tensors are analysed. It is shown that a Padé approximant, built on the correction tensors at the first three orders, outperforms the predictions for the effective transmissivity tensor. The new approach is illustrated and validated with numerical examples on model rough fractures.

**Key words:** lubrication theory, porous media

## 1. Introduction

Flow in fractures is a situation widely encountered in practical applications such as the leakage rate of seals (e.g. Pérez-Ràfols, Larsson & Almqvist 2016; Rui *et al.* 2023), gas leak through reservoir caps or gas recovery in fractured rocks (see Zhou *et al.* (2023) and references therein), among many others. Modelling flow in rough fractures

† Email address for correspondence: didier.lasseux@u-bordeaux.fr

<footnote n="1">¹CEA, DES, ISEC, DPME, SEME, Laboratoire d'Étanchéité, Université de Montpellier, 30207 Marcoule, France</footnote>

<footnote n="2">²División de Ciencias Básicas e Ingeniería, Universidad Autónoma Metropolitana-Iztapalapa, Av. Ferrocarril San Rafael Atlixco 186, Col. Leyes de Reforma, 09310, CDMX, Mexico</footnote>

<footnote n="3">³Institut de Mécanique des Fluides de Toulouse (IMFT), Université de Toulouse, CNRS, 31400 Toulouse, France</footnote>

<footnote n="4">⁴Université de Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400, Talence, France</footnote>
has been carried out using the average Reynolds equations over several decades. In early
works, surface roughness was treated as a stochastic variable in one dimension (Tzeng
& Saibel 1967); however, in the work by Patir & Cheng (1978), an average Reynolds
equation was proposed in terms of the macroscopic pressure gradient and shear flow
factors, thus making it possible to handle three-dimensional fractures. Using the volume
averaging method, Prat, Plouraboué & Letalleur (2002) later showed that these factors
were predictable from closure problem solutions in a representative periodic unit cell.

While many studies of flow through rough fractures have adopted no-slip boundary
conditions at the solid-fluid interface, this assumption is not always justified in situations
of practical interest such as non-wetting fluid flow (Lee, Yo & Lee 2013) or rarefied
gas flow (Wang, Tang & Jing 2019). In the slip regime, characterised by a Knudsen
number (i.e. the ratio of the gas molecules mean-free path to a characteristic constriction
length) smaller than ~0.1, gas flow can be modelled, at the local scale, with the Reynolds
(lubrication) equation that includes gas slippage effects at the walls (Zaouter, Lasseux &
Prat 2018). The same type of model operates as well for slip mechanisms different from
Knudsen effects.

Since fractures often exhibit a roughness that is multiscale in nature, upscaled flow
models are favoured in practice. Such models also have a Reynolds-like form involving
the macroscopic effective transmissivity tensor that lumps the geometrical and slip effects.
To separate the latter from the pure viscous contribution, a Knudsen number power-series
expansion can be performed, which provides the intrinsic transmissivity tensor as well as
a series of intrinsic correction tensors at multiple orders. These tensors can be considered
to be intrinsic as they only depend on the fracture geometry. They can be computed
from solving a set of sequentially coupled ancillary closure problems on a periodic unit
cell representative of the system under consideration. These problems can be relatively
difficult to solve and computationally costly, and the sequential coupling is prone to the
propagation of numerical errors in the solutions at the successive orders. Moreover, the
use of power-series expansions does not necessarily guarantee convergence and this poses
a serious problem in terms of accuracy.

To address this issue, a procedure is derived that considerably simplifies the
determination of the successive correction tensors as the solution of the first *N* closure
problems provides the correction tensors up to the 2*N* − 1 order. The analysis allows
determination of the symmetry and definiteness properties of the individual tensors. In
addition, it is shown that the power-series expansion can be employed to form a Padé
approximant, which, in its simplest form, is built using the first three correction tensors,
that are obtained from the solution of only the first two ancillary closure problems.

The article is organised as follows. In § 2, the local and macroscale models derived in Zaouter et al. (2018) for gas slip flow in a rough fracture are recalled. In addition, the associated closure problems resulting from power-series expansions in the Knudsen number are presented. In § 3, an integral approach is used to express the correction tensors at the different orders in the Knudsen number, showing that the computational requirements are reduced by a factor of almost two. The analysis also provides the proof of symmetry and semi-positive definiteness of the tensors. In addition, a Padé approximant is constructed requiring the minimum number of correction tensors to provide an alternative estimate of the effective transmissivity tensor. In § 4, numerical examples are used to validate the method derived here and illustrate the dependence of the transmissivity on the average Knudsen number, considering random rough fractures. The performance of the first-order (i.e. Klinkenberg) and higher-order corrections is explored. The Padé approximant is shown to outperform the predictions from the original power-series expansion. Finally, conclusions are presented in § 5.
has been carried out using the average Reynolds equations over several decades. In early
works, surface roughness was treated as a stochastic variable in one dimension (Tzeng
& Saibel 1967); however, in the work by Patir & Cheng (1978), an average Reynolds
equation was proposed in terms of the macroscopic pressure gradient and shear flow
factors, thus making it possible to handle three-dimensional fractures. Using the volume
averaging method, Prat, Plouraboué & Letalleur (2002) later showed that these factors
were predictable from closure problem solutions in a representative periodic unit cell.

While many studies of flow through rough fractures have adopted no-slip boundary
conditions at the solid-fluid interface, this assumption is not always justified in situations
of practical interest such as non-wetting fluid flow (Lee, Yo & Lee 2013) or rarefied
gas flow (Wang, Tang & Jing 2019). In the slip regime, characterised by a Knudsen
number (i.e. the ratio of the gas molecules mean-free path to a characteristic constriction
length) smaller than ~0.1, gas flow can be modelled, at the local scale, with the Reynolds
(lubrication) equation that includes gas slippage effects at the walls (Zaouter, Lasseux &
Prat 2018). The same type of model operates as well for slip mechanisms different from
Knudsen effects.

Since fractures often exhibit a roughness that is multiscale in nature, upscaled flow
models are favoured in practice. Such models also have a Reynolds-like form involving
the macroscopic effective transmissivity tensor that lumps the geometrical and slip effects.
To separate the latter from the pure viscous contribution, a Knudsen number power-series
expansion can be performed, which provides the intrinsic transmissivity tensor as well as
a series of intrinsic correction tensors at multiple orders. These tensors can be considered
to be intrinsic as they only depend on the fracture geometry. They can be computed
from solving a set of sequentially coupled ancillary closure problems on a periodic unit
cell representative of the system under consideration. These problems can be relatively
difficult to solve and computationally costly, and the sequential coupling is prone to the
propagation of numerical errors in the solutions at the successive orders. Moreover, the
use of power-series expansions does not necessarily guarantee convergence and this poses
a serious problem in terms of accuracy.

To address this issue, a procedure is derived that considerably simplifies the
determination of the successive correction tensors as the solution of the first *N* closure
problems provides the correction tensors up to the 2*N* − 1 order. The analysis allows
determination of the symmetry and definiteness properties of the individual tensors. In
addition, it is shown that the power-series expansion can be employed to form a Padé
approximant, which, in its simplest form, is built using the first three correction tensors,
that are obtained from the solution of only the first two ancillary closure problems.

The article is organised as follows. In § 2, the local and macroscale models derived
in Zaouter et al. (2018) for gas slip flow in a rough fracture are recalled. In addition,
the associated closure problems resulting from power-series expansions in the Knudsen
number are presented. In § 3, an integral approach is used to express the correction
tensors at the different orders in the Knudsen number, showing that the computational
requirements are reduced by a factor of almost two. The analysis also provides the proof of
symmetry and semi-positive definiteness of the tensors. In addition, a Padé approximant is
constructed requiring the minimum number of correction tensors to provide an alternative
estimate of the effective transmissivity tensor. In § 4, numerical examples are used to
validate the method derived here and illustrate the dependence of the transmissivity
on the average Knudsen number, considering random rough fractures. The performance
of the first-order (i.e. Klinkenberg) and higher-order corrections is explored. The Padé
approximant is shown to outperform the predictions from the original power-series
expansion. Finally, conclusions are presented in § 5.
## 2. Reynolds models at the local and macroscopic scales

The developments that follow are carried out considering the (slightly) compressible, isothermal, Newtonian and steady creeping flow of a gas within a fracture, assuming that Knudsen effects may be significant but in a rarefied regime corresponding to slip flow. Nevertheless, generalisation to slip flow in a broader context is straightforward, as it would require assuming a diffuse reflection at the solid–fluid interfaces and replacing the Knudsen number by the slip length to characteristic aperture size ratio (Bolaños & Vernescu 2017).

### 2.1. Local description

Derivation of the local Reynolds equation for flow in a rough fracture in the presence of slip was carried out in a previous work (see § 2 in Zaouter et al. 2018). The flow model at the roughness scale can be written in the following form:

$$ \nabla \cdot q = 0, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.1a) $$

$$ q = -\rho \frac{h^3}{12\mu} \left(1 + 6\xi \frac{\lambda}{h}\right) \nabla p, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.1b) $$

$$ \rho = \varphi(p), \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.1c) $$

$$ q \cdot n = 0, \quad \text{at } \mathcal{C}_{\beta\sigma}. \qquad (2.1d) $$

In the above equations, $q$ denotes the local mass flow rate per unit width of the fracture, and can be referred to as the local mass flux. In (2.1b), $h = h(x, y)$ represents the local aperture at a point $(x, y)$ in the mid-plane of the fracture (see figure 1) where the pressure is $p$, the density is $\rho$ and the mean-free path is $\lambda$, $\mu$ being the dynamic fluid viscosity, that is assumed constant. In addition, $\xi = (2 - \sigma_v)/\sigma_v$ is a parameter that accounts for the molecular reflection process at the solid–fluid interfaces, $\sigma_v$ being the tangential momentum accommodation coefficient (Maxwell 1879). Moreover, $\mathcal{A}_\beta$ represents the portion of the fracture that is occupied by the fluid phase, $\beta$, and $\mathcal{C}_{\beta\sigma}$ denotes the contours of the solid contact zones, $\sigma$ (where $h = 0$), that are eventually present between the two rough surfaces forming the fracture (see figure 1). The state equation, (2.1c) relates the density only to the fluid pressure, under the assumption of isothermal conditions. To arrive at this flow model, it is assumed that the roughness at both walls is slowly varying. This means that the slope of the two surfaces is everywhere small compared with unity, i.e. that $\varepsilon = h_\beta / \ell_\beta \ll 1$, $h_\beta$ being the characteristic aperture size and $\ell_\beta$ the characteristic length scale in the x- and y-directions of the fracture over which $h_\beta$ varies significantly. The model is $O(\varepsilon^2(1 + \xi Kn))$ accurate and its use is constrained by the Knudsen number, $Kn = \lambda/h_\beta \lesssim 0.1$.

### 2.2. Macroscopic model

The Reynolds model given in (2.1) operates locally in the two-dimensional domain corresponding to the mid-plane of the fracture. A macroscopic description over a representative element of the fracture is nevertheless of major interest and was developed in Zaouter et al. (2018) (cf. § 3). Under the assumption of a separation of length scales,
## 2. Reynolds models at the local and macroscopic scales

The developments that follow are carried out considering the (slightly) compressible, isothermal, Newtonian and steady creeping flow of a gas within a fracture, assuming that Knudsen effects may be significant but in a rarefied regime corresponding to slip flow. Nevertheless, generalisation to slip flow in a broader context is straightforward, as it would require assuming a diffuse reflection at the solid–fluid interfaces and replacing the Knudsen number by the slip length to characteristic aperture size ratio (Bolaños & Vernescu 2017).

### 2.1. Local description

Derivation of the local Reynolds equation for flow in a rough fracture in the presence of slip was carried out in a previous work (see § 2 in Zaouter et al. 2018). The flow model at the roughness scale can be written in the following form:

$$ \nabla \cdot \mathbf{q} = 0, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.1a) $$

$$ \mathbf{q} = -\rho \frac{h^3}{12\mu} \left(1 + 6\xi \frac{\lambda}{h}\right) \nabla p, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.1b) $$

$$ \rho = \varphi(p), \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.1c) $$

$$ \mathbf{q} \cdot \mathbf{n} = 0, \quad \text{at } \mathcal{C}_{\beta\sigma}. \qquad (2.1d) $$

In the above equations, $q$ denotes the local mass flow rate per unit width of the fracture, and can be referred to as the local mass flux. In (2.1b), $h = h(x, y)$ represents the local aperture at a point $(x, y)$ in the mid-plane of the fracture (see figure 1) where the pressure is $p$, the density is $\rho$ and the mean-free path is $\lambda$, $\mu$ being the dynamic fluid viscosity, that is assumed constant. In addition, $\xi = (2 - \sigma_v)/\sigma_v$ is a parameter that accounts for the molecular reflection process at the solid–fluid interfaces, $\sigma_v$ being the tangential momentum accommodation coefficient (Maxwell 1879). Moreover, $\mathcal{A}_{\beta}$ represents the portion of the fracture that is occupied by the fluid phase, $\beta$, and $\mathcal{C}_{\beta\sigma}$ denotes the contours of the solid contact zones, $\sigma$ (where $h = 0$), that are eventually present between the two rough surfaces forming the fracture (see figure 1). The state equation, (2.1c) relates the density only to the fluid pressure, under the assumption of isothermal conditions. To arrive at this flow model, it is assumed that the roughness at both walls is slowly varying. This means that the slope of the two surfaces is everywhere small compared with unity, i.e. that $\varepsilon = h_{\beta}/\ell_{\beta} \ll 1$, $h_{\beta}$ being the characteristic aperture size and $\ell_{\beta}$ the characteristic length scale in the x- and y-directions of the fracture over which $h_{\beta}$ varies significantly. The model is $O(\varepsilon^2(1 + \xi Kn))$ accurate and its use is constrained by the Knudsen number, $Kn = \lambda/h_{\beta} \lesssim 0.1$.

### 2.2. Macroscopic model

The Reynolds model given in (2.1) operates locally in the two-dimensional domain corresponding to the mid-plane of the fracture. A macroscopic description over a representative element of the fracture is nevertheless of major interest and was developed in Zaouter et al. (2018) (cf. § 3). Under the assumption of a separation of length scales,
<figure><img src="image_3.png" /><figcaption>Figure 1. Top: sketch of the system consisting of a fracture between two rough surfaces. Note that $h(x, y)$ represents the local aperture at any point in the mid-plane of the fracture. Bottom left: top view of part of the fracture. Bottom right: detail of the two-dimensional domain in which the Reynolds model applies, including the corresponding phases and characteristic lengths, the periodic unit cell, the axes and the notations.</figcaption></figure>

i.e. $\ell_{\beta} \ll L$, $L$ being the size of the system, it is given by

$$ \nabla \cdot \langle q \rangle = 0, \qquad (2.2a) $$

$$ \langle q \rangle = -\langle \rho \rangle^{\beta} \frac{\mathbf{K}}{\mu} \cdot \nabla \langle p \rangle^{\beta}, \qquad (2.2b) $$

$$ \langle \rho \rangle^{\beta} = \varphi(\langle p \rangle^{\beta}). \qquad (2.2c) $$

In this model, the superficial and intrinsic surface averages of a variable, $\psi$, taking values in $\mathcal{A}_{\beta}$ are respectively defined as

$$ \langle \psi \rangle = \frac{1}{A} \int_{\mathcal{A}_{\beta}} \psi \, dA, \qquad (2.3a) $$

$$ \langle \psi \rangle^{\beta} = \frac{1}{A_{\beta}} \int_{\mathcal{A}_{\beta}} \psi \, dA. \qquad (2.3b) $$

In these definitions, $A$ and $A_{\beta}$ are the areas of $\mathcal{A}$ and $\mathcal{A}_{\beta}$, respectively, $\mathcal{A}$ being the averaging surface identified as the representative (periodic) unit cell of the fracture aperture field, and $\mathcal{A}_{\beta}$ the portion of $\mathcal{A}$ occupied by the $\beta$-phase, i.e. where $h \neq 0$ (cf. figure 1). Equation (2.2c) is the formal state equation at the macroscopic scale that follows from (2.1c). The above macroscopic model is obtained assuming a slightly compressible flow characterised by $\rho - \langle \rho \rangle^{\beta} \ll \langle \rho \rangle^{\beta}$. In the momentum equation (2.2b), $\mathbf{K}$ is the apparent (effective) transmissivity tensor that is obtained from the solution of the
<figure><img src="image_3.png" /><figcaption>Figure 1. Top: sketch of the system consisting of a fracture between two rough surfaces. Note that $h(x, y)$ represents the local aperture at any point in the mid-plane of the fracture. Bottom left: top view of part of the fracture. Bottom right: detail of the two-dimensional domain in which the Reynolds model applies, including the corresponding phases and characteristic lengths, the periodic unit cell, the axes and the notations.</figcaption></figure>

i.e. $l_β \ll L$, $L$ being the size of the system, it is given by

$$ \nabla \cdot \langle q \rangle = 0, \qquad (2.2a) $$

$$ \langle q \rangle = -\langle \rho \rangle^{\beta} \frac{\mathbf{K}}{\mu} \cdot \nabla \langle p \rangle^{\beta}, \qquad (2.2b) $$

$$ \langle \rho \rangle^{\beta} = \varphi(\langle p \rangle^{\beta}). \qquad (2.2c) $$

In this model, the superficial and intrinsic surface averages of a variable, $\psi$, taking values in $\mathcal{A}_β$ are respectively defined as

$$ \langle \psi \rangle = \frac{1}{A} \int_{\mathcal{A}_\beta} \psi \, dA, \qquad (2.3a) $$

$$ \langle \psi \rangle^{\beta} = \frac{1}{A_{\beta}} \int_{\mathcal{A}_{\beta}} \psi \, dA. \qquad (2.3b) $$

In these definitions, $A$ and $A_β$ are the areas of $\mathcal{A}$ and $\mathcal{A}_β$, respectively, $\mathcal{A}$ being the averaging surface identified as the representative (periodic) unit cell of the fracture aperture field, and $\mathcal{A}_β$ the portion of $\mathcal{A}$ occupied by the $β$-phase, i.e. where $h \neq 0$ (cf. figure 1). Equation (2.2c) is the formal state equation at the macroscopic scale that follows from (2.1c). The above macroscopic model is obtained assuming a slightly compressible flow characterised by $\rho - \langle \rho \rangle^{\beta} \ll \langle \rho \rangle^{\beta}$. In the momentum equation (2.2b), $\mathbf{K}$ is the apparent (effective) transmissivity tensor that is obtained from the solution of the
following closure problem in a periodic (two-dimensional) unit cell:

$$ \nabla \cdot (k(\nabla b + \mathbf{l})) = \mathbf{0}, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.4a) $$

$$ \mathbf{n} \cdot k(\nabla \mathbf{b} + \mathbf{l}) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (2.4b) $$

$$ \mathbf{b}(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}(\mathbf{r}), \quad i = 1, 2, \qquad (2.4c) $$

$$ \mathbf{K} = \langle k(\nabla \mathbf{b} + \mathbf{l}) \rangle. \qquad (2.4d) $$

Here, the local transmissivity $k$ is defined as

$$ k = k_0 + \alpha k_1, \quad k_0 = \frac{h^3}{12}, \quad k_1 = \frac{h^2}{12}, \quad \alpha = 6\xi\bar{\lambda}, \qquad (2.5) $$

and $\bar{\lambda}$ is the mean-free path corresponding to the average density, $\langle \rho \rangle^\beta$. Note that $\xi\bar{\lambda}$ can be identified as the slip length in the general case, beyond the context of gas flow in the presence of Knudsen effects. Moreover, $\mathbf{l}$ is the identity tensor, and in (2.4c), $\mathbf{l}_i$ represents the lattice vector of the two-dimensional unit cell in the $i$th-direction (i.e. $\mathbf{l}_1 = \ell_x \mathbf{e}_x$, $\mathbf{l}_2 = \ell_y \mathbf{e}_y$, ($\mathbf{e}_x, \mathbf{e}_y$) forming the system of coordinates in the plane of the fracture, see figure 1). Note that the problem defined in (2.4) may appear as ill posed since $\mathbf{b}$ is defined to within an additive constant. However, from (2.4d), it is clear that this constant plays no role in the prediction of $\mathbf{K}$ and, hence, it is not required to specify it. The field of $\mathbf{b}$, compliant with the closure procedure, is nevertheless constrained by $\langle \mathbf{b} \rangle^\beta = \mathbf{0}$.

The second-order effective transmissivity tensor, $\mathbf{K}$, is apparent (i.e. non-intrinsic) in the sense that it depends on the flow characteristic through $\bar{\lambda}$. This means that if one is interested in studying the dependence of $\mathbf{K}$ on $\bar{\lambda}$, it is necessary to solve the boundary-value problem given in (2.4) for any value of this parameter. To elucidate this dependence, $\mathbf{b}$ can be expanded in power series of $\bar{K}n = \bar{\lambda}/h_\beta$, which, in dimensional form, translates into $\mathbf{b} = \sum_{j=0}^{m} \alpha^j \mathbf{b}_j$. Accordingly, the power-series expansion, $\hat{\mathbf{K}}_m$, of $\mathbf{K}$ at the $m$th-order in $\bar{K}n$, is obtained as

$$ \hat{\mathbf{K}}_m = \sum_{j=0}^{m} \alpha^j \mathbf{K}_j. \qquad (2.6) $$

Note that, although the flux to force linearity is preserved during the averaging procedure, the linearity between $\mathbf{K}$ and $\xi\bar{\lambda}$ does not persist, in general, due to the disordered nature of the fracture aperture field that is captured in the upscaling process. For this reason, it is relevant to consider the above expansion beyond the first order. Nevertheless, the use of a power-series expansion, as expressed in (2.6), although of practical value, shall be considered with caution regarding the radius of convergence with respect to $\bar{K}n$. A thorough analysis of this convergence criterion, which may depend on the fracture aperture field, is certainly of utmost importance but goes beyond the scope of the present work. To circumvent this issue, the correction tensors series can be used to build a Padé approximant as will be detailed in § 3.5.

With the above expansions, the second-order tensors, $\mathbf{K}_j$ for $j \ge 0$, are obtained from the solution of the closure problem at the corresponding order given by

Order $j = 0$

$$ \nabla \cdot (k_0(\nabla \mathbf{b}_0 + \mathbf{l})) = \mathbf{0}, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.7a) $$

$$ \mathbf{n} \cdot k_0(\nabla \mathbf{b}_0 + \mathbf{l}) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (2.7b) $$

$$ \mathbf{b}_0(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}_0(\mathbf{r}), \quad i = 1, 2, \qquad (2.7c) $$

$$ \mathbf{K}_0 = \langle k_0(\nabla \mathbf{b}_0 + \mathbf{l}) \rangle. \qquad (2.7d) $$
following closure problem in a periodic (two-dimensional) unit cell:

$$ \nabla \cdot (k(\nabla b + \mathbf{l})) = \mathbf{0}, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.4a) $$

$$ \mathbf{n} \cdot k(\nabla \mathbf{b} + \mathbf{l}) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (2.4b) $$

$$ \mathbf{b}(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}(\mathbf{r}), \quad i = 1, 2, \qquad (2.4c) $$

$$ \mathbf{K} = \langle k(\nabla \mathbf{b} + \mathbf{l}) \rangle. \qquad (2.4d) $$

Here, the local transmissivity $k$ is defined as

$$ k = k_0 + \alpha k_1, \quad k_0 = \frac{h^3}{12}, \quad k_1 = \frac{h^2}{12}, \quad \alpha = 6\xi\bar{\lambda}, \qquad (2.5) $$

and $\bar{\lambda}$ is the mean-free path corresponding to the average density, $\langle \rho \rangle^\beta$. Note that $\xi\bar{\lambda}$ can be identified as the slip length in the general case, beyond the context of gas flow in the presence of Knudsen effects. Moreover, $\mathbf{l}$ is the identity tensor, and in (2.4c), $\mathbf{l}_i$ represents the lattice vector of the two-dimensional unit cell in the $i$th-direction (i.e. $\mathbf{l}_1 = \ell_x e_x$, $\mathbf{l}_2 = \ell_y e_y$, ($e_x$, $e_y$) forming the system of coordinates in the plane of the fracture, see figure 1). Note that the problem defined in (2.4) may appear as ill posed since $\mathbf{b}$ is defined to within an additive constant. However, from (2.4d), it is clear that this constant plays no role in the prediction of $\mathbf{K}$ and, hence, it is not required to specify it. The field of $\mathbf{b}$, compliant with the closure procedure, is nevertheless constrained by $\langle \mathbf{b} \rangle^\beta = \mathbf{0}$.

The second-order effective transmissivity tensor, $\mathbf{K}$, is apparent (i.e. non-intrinsic) in the sense that it depends on the flow characteristic through $\bar{\lambda}$. This means that if one is interested in studying the dependence of $\mathbf{K}$ on $\bar{\lambda}$, it is necessary to solve the boundary-value problem given in (2.4) for any value of this parameter. To elucidate this dependence, $\mathbf{b}$ can be expanded in power series of $\bar{K}n = \bar{\lambda}/h_\beta$, which, in dimensional form, translates into $\mathbf{b} = \sum_{j=0}^{m} \alpha^j \mathbf{b}_j$. Accordingly, the power-series expansion, $\hat{\mathbf{K}}_m$, of $\mathbf{K}$ at the $m$th-order in $\bar{K}n$, is obtained as

$$ \hat{\mathbf{K}}_m = \sum_{j=0}^{m} \alpha^j \mathbf{K}_j. \qquad (2.6) $$

Note that, although the flux to force linearity is preserved during the averaging procedure, the linearity between $\mathbf{K}$ and $\xi\bar{\lambda}$ does not persist, in general, due to the disordered nature of the fracture aperture field that is captured in the upscaling process. For this reason, it is relevant to consider the above expansion beyond the first order. Nevertheless, the use of a power-series expansion, as expressed in (2.6), although of practical value, shall be considered with caution regarding the radius of convergence with respect to $\bar{K}n$. A thorough analysis of this convergence criterion, which may depend on the fracture aperture field, is certainly of utmost importance but goes beyond the scope of the present work. To circumvent this issue, the correction tensors series can be used to build a Padé approximant as will be detailed in § 3.5.

With the above expansions, the second-order tensors, $\mathbf{K}_j$ for $j \ge 0$, are obtained from the solution of the closure problem at the corresponding order given by

Order $j = 0$

$$ \nabla \cdot (k_0(\nabla \mathbf{b}_0 + \mathbf{l})) = \mathbf{0}, \quad \text{in } \mathcal{A}_{\beta}, \qquad (2.7a) $$

$$ \mathbf{n} \cdot k_0(\nabla \mathbf{b}_0 + \mathbf{l}) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (2.7b) $$

$$ \mathbf{b}_0(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}_0(\mathbf{r}), \quad i = 1, 2, \qquad (2.7c) $$

$$ \mathbf{K}_0 = \langle k_0(\nabla \mathbf{b}_0 + \mathbf{l}) \rangle. \qquad (2.7d) $$
Order $j \ge 1$

$$ \nabla \cdot (k_0 \nabla \mathbf{b}_j + k_1 (\nabla \mathbf{b}_{j-1} + \delta_{j1}^K \mathbf{I})) = \mathbf{0}, \quad \text{in } \mathcal{A}_\beta, \qquad (2.8a) $$

$$ \mathbf{n} \cdot (k_0 \nabla \mathbf{b}_j + k_1 (\nabla \mathbf{b}_{j-1} + \delta_{j1}^K \mathbf{I})) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (2.8b) $$

$$ \mathbf{b}_j(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}_j(\mathbf{r}), \quad i = 1, 2, \qquad (2.8c) $$

$$ \mathbf{K}_j = \langle k_0 \nabla \mathbf{b}_j + k_1 (\nabla \mathbf{b}_{j-1} + \delta_{j1}^K \mathbf{I}) \rangle. \qquad (2.8d) $$

In (2.8), $\delta_{j1}^K$ represents the Kronecker delta. As in the problem for $\mathbf{b}$, no constraint is assigned for $\mathbf{b}_j$, which is hence defined to within an arbitrary additive constant that has no influence on $\mathbf{K}_j$.

It must be noted that the zeroth-order problem yields the intrinsic transmissivity tensor that is the only relevant effective coefficient at the macroscale in the absence of slip. The successive orders $j \ge 1$ provide the effective correction tensors due to slip effects. Conversely to $\mathbf{K}$, the tensors $\mathbf{K}_0, \mathbf{K}_1, \dots$, are all intrinsic to the fracture under consideration as they only depend on its aperture field, making $\hat{\mathbf{K}}_m$ fully predictive. Moreover, they allow distinguishing of the contributions from slip effects from purely viscous effects. In particular, at the first order in $\overline{\mathbf{K}}n$, the correction leads to the Klinkenberg relationship when the gas can be considered as ideal (see end of § 3, p. 429 in Zaouter et al. 2018), similar to the classical correction to Darcy’s law in the context of slip flow in porous media (Lasseux et al. 2014; Lasseux, Valdés-Parada & Porter 2016).

From the above expansion in $\overline{\mathbf{K}}n$, it appears that the determination of the $j$th-order correction tensor requires the solution of the series of $j+1$ problems (from order 0 to $j$). This can be extremely costly in terms of computational time and resources to achieve enough accuracy, in particular because the closure problem at order $j$ involves the gradient of the solution at the $j-1$ order, hence contributing to a propagation of errors detrimental to the final solution. This motivates exploring of simplifications in the procedure to obtain the correction tensors, which are presented in the following section. The approach is based upon an integral formula, that results from the divergence theorem, which is used recursively to express $\mathbf{K}_{2N-1}$ in terms of the first $N$ closure problems solution. Moreover, this approach allows determination of the symmetry and definiteness properties of the correction tensors in a very simple way. In addition, the use of a simple Padé approximant is explored, that only requires the solution of the first two closure problems.

### 3. An efficient procedure to obtain the correction tensors $\mathbf{K}_j, j \ge 1$

#### 3.1. Integral formula

In order to progress towards a simplified method of determining $\mathbf{K}_j$, it is of interest to derive an integral formula that will be extensively used in the following. Consider an arbitrary second-order tensor field, $\mathbf{A}$, and a vector field $\mathbf{a}$ defined in a domain $\Omega_\beta$ of boundary $\partial\Omega_\beta$ and having sufficient regularity. On the basis of the divergence theorem, they satisfy the following identity:

$$ \int_{\Omega_\beta} (\nabla \cdot \mathbf{A}) \mathbf{a} \, d\Omega + \int_{\Omega_\beta} \mathbf{A}^T \cdot \nabla \mathbf{a} \, d\Omega = \int_{\partial\Omega_\beta} \mathbf{n} \cdot \mathbf{A} \mathbf{a} \, dS. \qquad (3.1) $$

Here, $\mathbf{n}$ is the unit normal vector at $\partial\Omega_\beta$ pointing out of $\Omega_\beta$. When this formula is considered in the periodic unit cell, $\mathcal{A}$, representative of the fracture ($\Omega_\beta \equiv \mathcal{A}_\beta$) with
Order $j \ge 1$

$$ \nabla \cdot (k_0 \nabla \mathbf{b}_j + k_1 (\nabla \mathbf{b}_{j-1} + \delta_{j1}^K \mathbf{I})) = \mathbf{0}, \quad \text{in } \mathcal{A}_\beta, \qquad (2.8a) $$

$$ \mathbf{n} \cdot (k_0 \nabla \mathbf{b}_j + k_1 (\nabla \mathbf{b}_{j-1} + \delta_{j1}^K \mathbf{I})) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (2.8b) $$

$$ \mathbf{b}_j(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}_j(\mathbf{r}), \quad i = 1, 2, \qquad (2.8c) $$

$$ \mathbf{K}_j = \langle k_0 \nabla \mathbf{b}_j + k_1 (\nabla \mathbf{b}_{j-1} + \delta_{j1}^K \mathbf{I}) \rangle. \qquad (2.8d) $$

In (2.8), $\delta_{j1}^K$ represents the Kronecker delta. As in the problem for $\mathbf{b}$, no constraint is assigned for $\mathbf{b}_j$, which is hence defined to within an arbitrary additive constant that has no influence on $\mathbf{K}_j$.

It must be noted that the zeroth-order problem yields the intrinsic transmissivity tensor that is the only relevant effective coefficient at the macroscale in the absence of slip. The successive orders $j \ge 1$ provide the effective correction tensors due to slip effects. Conversely to $\mathbf{K}$, the tensors $\mathbf{K}_0, \mathbf{K}_1, \dots$, are all intrinsic to the fracture under consideration as they only depend on its aperture field, making $\hat{\mathbf{K}}_m$ fully predictive. Moreover, they allow distinguishing of the contributions from slip effects from purely viscous effects. In particular, at the first order in $\overline{\mathbf{K}}n$, the correction leads to the Klinkenberg relationship when the gas can be considered as ideal (see end of § 3, p. 429 in Zaouter et al. 2018), similar to the classical correction to Darcy’s law in the context of slip flow in porous media (Lasseux et al. 2014; Lasseux, Valdés-Parada & Porter 2016).

From the above expansion in $\overline{\mathbf{K}}n$, it appears that the determination of the $j$th-order correction tensor requires the solution of the series of $j+1$ problems (from order 0 to $j$). This can be extremely costly in terms of computational time and resources to achieve enough accuracy, in particular because the closure problem at order $j$ involves the gradient of the solution at the $j-1$ order, hence contributing to a propagation of errors detrimental to the final solution. This motivates exploring of simplifications in the procedure to obtain the correction tensors, which are presented in the following section. The approach is based upon an integral formula, that results from the divergence theorem, which is used recursively to express $\mathbf{K}_{2N-1}$ in terms of the first $N$ closure problems solution. Moreover, this approach allows determination of the symmetry and definiteness properties of the correction tensors in a very simple way. In addition, the use of a simple Padé approximant is explored, that only requires the solution of the first two closure problems.

### 3. An efficient procedure to obtain the correction tensors $\mathbf{K}_j, j \ge 1$

#### 3.1. Integral formula

In order to progress towards a simplified method of determining $\mathbf{K}_j$, it is of interest to derive an integral formula that will be extensively used in the following. Consider an arbitrary second-order tensor field, $\mathbf{A}$, and a vector field $\mathbf{a}$ defined in a domain $\Omega_\beta$ of boundary $\partial\Omega_\beta$ and having sufficient regularity. On the basis of the divergence theorem, they satisfy the following identity:

$$ \int_{\Omega_{\beta}} (\nabla \cdot \mathbf{A}) \mathbf{a} \, d\Omega + \int_{\Omega_{\beta}} \mathbf{A}^T \cdot \nabla \mathbf{a} \, d\Omega = \int_{\partial\Omega_{\beta}} \mathbf{n} \cdot \mathbf{A} \mathbf{a} \, dS. \qquad (3.1) $$

Here, $\mathbf{n}$ is the unit normal vector at $\partial\Omega_\beta$ pointing out of $\Omega_\beta$. When this formula is considered in the periodic unit cell, $\mathcal{A}$, representative of the fracture ($\Omega_\beta \equiv \mathcal{A}_\beta$) with
periodic fields **A** and **a**, further considering **A** to be solenoidal and satisfying **n**·**A** = **0** at $\mathcal{C}_{\beta\sigma}$, the above identity simplifies to give

$$ \langle \mathbf{A}^T \cdot \nabla \mathbf{a} \rangle = \mathbf{0}, \qquad (3.2) $$

where $\langle \psi \rangle$ is the surface superficial average of $\psi$ defined in (2.3a). With this result at hand, the transmissivity tensors can be reformulated. Important properties for **K** and **K**<sub>0</sub> are first analysed before focusing on **K**<sub>j</sub>, j ≥ 1.

## 3.2. Properties of K and K<sub>0</sub>

In Lasseux & Valdés-Parada (2017) (see section II C), **K** and **K**<sub>0</sub> were shown to be symmetric tensors. The analysis is now completed with a proof of their positiveness. To do so, the relationship given in (3.2) is considered taking **A** = k(∇**b** + **I**) (see (2.4)) and **a** = **b**. This leads to $\langle k\nabla b \rangle + \langle k\nabla b^T \cdot \nabla b \rangle = \mathbf{0}$. Adding the transpose of this result to the expression of **K** given in (2.4d) yields $\mathbf{K} = \langle k\nabla b + k\mathbf{I} + k\nabla b^T + k\nabla b^T \cdot \nabla b \rangle$, that is

$$ \mathbf{K} = \langle k(\nabla b^T + \mathbf{I}) \cdot (\nabla b + \mathbf{I}) \rangle. \qquad (3.3) $$

Similarly, taking **A** = k<sub>0</sub>(∇**b**<sub>0</sub> + **I**), (see (2.7)) and **a** = **b**<sub>0</sub> in relation (3.2) gives

$$ \mathbf{K}_0 = \langle k_0(\nabla b_0^T + \mathbf{I}) \cdot (\nabla b_0 + \mathbf{I}) \rangle. \qquad (3.4) $$

The last two expressions confirm that both **K** and **K**<sub>0</sub> are symmetric tensors and, moreover, that they are semi-definite positive. Indeed, for any arbitrary constant non-zero vector, $\omega$, the quantity $\omega \cdot \mathbf{K} \cdot \omega = \langle k((\nabla b + \mathbf{I}) \cdot \omega)^2 \rangle$ is positive (and similarly for $\omega \cdot \mathbf{K}_0 \cdot \omega$). Therefore, **K** and **K**<sub>0</sub> admit an inverse, which indicates that (2.2b) can be used in a reciprocal form, regardless flow is in the slip regime or not. These results, in terms of symmetry, positiveness (and consequently on the existence of an inverse) can straightforwardly be generalised to many diffusive processes in porous media, for instance, for the effective diffusivity tensor for passive mass diffusion (i.e. with no adsorption or chemical reaction) in homogeneous (Whitaker 1999, chap. 1) and heterogeneous porous media, the effective conductivity tensor for thermal conduction in a porous medium (Whitaker 1999, chap. 2) and the effective permeability for one-phase creeping flow in heterogeneous porous media (Whitaker 1999, chap. 5). Finally, the fact that **K** is symmetric and positive is consistent with the Clausius–Duhem inequality.

## 3.3. Reformulation of K<sub>1</sub>

The interest is now focused on the first-order correction tensor, **K**<sub>1</sub>. Applying the integral formula given in (3.2) with **A** = k<sub>0</sub>∇**b**<sub>1</sub> + k<sub>1</sub>(∇**b**<sub>0</sub> + **I**), (see (2.8) with j = 1) and **a** = **b**<sub>0</sub>, and taking the transpose of the result, yields $\langle k_0\nabla b_0^T \cdot \nabla b_1 \rangle = -\langle k_1\nabla b_0^T \cdot (\nabla b_0 + \mathbf{I}) \rangle$. In addition, using again (3.2) with **A** = k<sub>0</sub>(∇**b**<sub>0</sub> + **I**), (see (2.7)) and **a** = **b**<sub>1</sub>, gives $\langle k_0\nabla b_0^T \cdot \nabla b_1 \rangle = -\langle k_0\nabla b_1 \rangle$. When the combination of these two expressions is substituted back into the definition of **K**<sub>1</sub> given in (2.8d) (taking j = 1), it follows that

$$ \mathbf{K}_1 = \langle k_1(\nabla b_0^T + \mathbf{I}) \cdot (\nabla b_0 + \mathbf{I}) \rangle. \qquad (3.5) $$

This last result has many important consequences. First, it shows that the first-order slip correction tensor can be obtained from the solution of the zeroth-order closure problem that hence provides both **K**<sub>0</sub> and **K**<sub>1</sub>. The above result shows that, if one is restricted to the classical first-order correction at the macroscopic scale, the zeroth-order closure problem
periodic fields **A** and **a**, further considering **A** to be solenoidal and satisfying **n**·**A** = **0** at $\mathcal{C}_{\beta\sigma}$, the above identity simplifies to give

$$ \langle \mathbf{A}^T \cdot \nabla \mathbf{a} \rangle = 0, \qquad (3.2) $$

where $\langle \psi \rangle$ is the surface superficial average of $\psi$ defined in (2.3a). With this result at hand, the transmissivity tensors can be reformulated. Important properties for **K** and **K**<sub>0</sub> are first analysed before focusing on **K**<sub>j</sub>, j ≥ 1.

## 3.2. Properties of K and K<sub>0</sub>

In Lasseux & Valdés-Parada (2017) (see section II C), **K** and **K**<sub>0</sub> were shown to be symmetric tensors. The analysis is now completed with a proof of their positiveness. To do so, the relationship given in (3.2) is considered taking **A** = k(∇**b** + **I**) (see (2.4)) and **a** = **b**. This leads to $\langle k\nabla b \rangle + \langle k\nabla b^T \cdot \nabla b \rangle = 0$. Adding the transpose of this result to the expression of **K** given in (2.4d) yields $\mathbf{K} = \langle k\nabla b + k\mathbf{I} + k\nabla b^T + k\nabla b^T \cdot \nabla b \rangle$, that is

$$ \mathbf{K} = \langle k(\nabla b^T + \mathbf{I}) \cdot (\nabla b + \mathbf{I}) \rangle. \qquad (3.3) $$

Similarly, taking **A** = k<sub>0</sub>(∇**b**<sub>0</sub> + **I**) (see (2.7)) and **a** = **b**<sub>0</sub> in relation (3.2) gives

$$ \mathbf{K}_0 = \langle k_0(\nabla b_0^T + \mathbf{I}) \cdot (\nabla b_0 + \mathbf{I}) \rangle. \qquad (3.4) $$

The last two expressions confirm that both **K** and **K**<sub>0</sub> are symmetric tensors and, moreover, that they are semi-definite positive. Indeed, for any arbitrary constant non-zero vector, $\omega$, the quantity $\omega \cdot \mathbf{K} \cdot \omega = \langle k((\nabla b + \mathbf{I}) \cdot \omega)^2 \rangle$ is positive (and similarly for $\omega \cdot \mathbf{K}_0 \cdot \omega$). Therefore, **K** and **K**<sub>0</sub> admit an inverse, which indicates that (2.2b) can be used in a reciprocal form, regardless flow is in the slip regime or not. These results, in terms of symmetry, positiveness (and consequently on the existence of an inverse) can straightforwardly be generalised to many diffusive processes in porous media, for instance, for the effective diffusivity tensor for passive mass diffusion (i.e. with no adsorption or chemical reaction) in homogeneous (Whitaker 1999, chap. 1) and heterogeneous porous media, the effective conductivity tensor for thermal conduction in a porous medium (Whitaker 1999, chap. 2) and the effective permeability for one-phase creeping flow in heterogeneous porous media (Whitaker 1999, chap. 5). Finally, the fact that **K** is symmetric and positive is consistent with the Clausius–Duhem inequality.

## 3.3. Reformulation of K<sub>1</sub>

The interest is now focused on the first-order correction tensor, **K**<sub>1</sub>. Applying the integral formula given in (3.2) with **A** = k<sub>0</sub>∇**b**<sub>1</sub> + k<sub>1</sub>(∇**b**<sub>0</sub> + **I**), (see (2.8) with j = 1) and **a** = **b**<sub>0</sub>, and taking the transpose of the result, yields $\langle k_0\nabla b_0^T \cdot \nabla b_1 \rangle = -\langle k_1\nabla b_0^T \cdot (\nabla b_0 + \mathbf{I}) \rangle$. In addition, using again (3.2) with **A** = k<sub>0</sub>(∇**b**<sub>0</sub> + **I**), (see (2.7)) and **a** = **b**<sub>1</sub>, gives $\langle k_0\nabla b_0^T \cdot \nabla b_1 \rangle = -\langle k_0\nabla b_1 \rangle$. When the combination of these two expressions is substituted back into the definition of **K**<sub>1</sub> given in (2.8d) (taking j = 1), it follows that

$$ \mathbf{K}_1 = \langle k_1(\nabla b_0^T + \mathbf{I}) \cdot (\nabla b_0 + \mathbf{I}) \rangle. \qquad (3.5) $$

This last result has many important consequences. First, it shows that the first-order slip correction tensor can be obtained from the solution of the zeroth-order closure problem that hence provides both **K**<sub>0</sub> and **K**<sub>1</sub>. The above result shows that, if one is restricted to the classical first-order correction at the macroscopic scale, the zeroth-order closure problem
(which physically corresponds to flow under no-slip conditions) contains all the physical information that is necessary to predict both the intrinsic transmissivity and first-order correction tensors. Second, it readily proves that $K_1$ is a symmetric, semi-definite positive tensor. Therefore, if a first-order correction is considered reasonable, which is the case when the Knudsen number is small enough (see some examples in § 4 of Zaouter et al. 2018), only the solution of (2.7) is required and $\hat{K}_1$ (see (2.6)) is symmetric, semi-definite positive. Consequently, it admits an inverse. These conclusions are similar to those obtained for slip flow in porous media analysed in Lasseux, Zaouter & Valdés-Parada (2023).

At this point, it is of interest to investigate the slip correction tensors of higher orders.

### 3.4. Alternative expression of $K_j$, $j \ge 2$

A simpler expression of $K_j$ compared with (2.8d), $j \ge 2$, can be obtained by first considering (3.2), taking $\mathbf{A} \equiv k_0(\nabla b_0 + \mathbf{I})$ (see (2.7)) and $\mathbf{a} \equiv b_j$, which gives $\langle k_0 \nabla b_j \rangle = -\langle k_0 \nabla b_0^T \cdot \nabla b_j \rangle$. Similarly, taking $\mathbf{A} \equiv k_0 \nabla b_1 + k_1(\nabla b_0 + \mathbf{I})$ (see (2.8) with $j=1$) and $\mathbf{a} \equiv b_{j-1}$ yields $\langle k_1 \nabla b_{j-1} \rangle = -\langle (k_0 \nabla b_1^T + k_1 \nabla b_0^T) \cdot \nabla b_{j-1} \rangle$. When these two results are substituted back into the expression of $K_j$ given in (2.8d) ($j > 1$), the following alternative expression of this tensor is obtained:

$$ K_j = -\langle k_0 \nabla b_1^T \cdot \nabla b_{j-1} \rangle. \quad (3.6) $$

To arrive at this result, the fact that $\langle \nabla b_0^T \cdot (k_0 \nabla b_j + k_1 \nabla b_{j-1}) \rangle = \mathbf{0}$ was employed. This follows from an additional use of (3.2) with $\mathbf{A} \equiv k_0 \nabla b_j + k_1 \nabla b_{j-1}$ (see (2.8)) and $\mathbf{a} \equiv b_0$ and the application of a transpose.

If $j=2$, no further step is required as (3.6) gives

$$ K_2 = -\langle k_0 \nabla b_1^T \cdot \nabla b_1 \rangle. \quad (3.7) $$

If $j \ge 2$, the next step, taking (3.6) as the initial stage, is to further employ the integral formula (3.2) repeatedly taking $\mathbf{A} \equiv k_0 \nabla b_{j-n} + k_1 \nabla b_{j-n-1}$ and $\mathbf{a} \equiv b_n$ with $n$ increasing from 1 to $(j-1)/2$ if $j$ is odd or from 1 to $j/2$ if $j$ is even. After substituting the ensuing relationships into the successive expressions of $K_j$ at each value of $n$, this finally yields (note that this procedure is indeed also valid for $j=2$)

$$ K_j = \begin{cases} \langle k_1 \nabla b_{(j-1)/2}^T \cdot \nabla b_{(j-1)/2} \rangle, & j \text{ odd,} \\ -\langle k_0 \nabla b_{j/2}^T \cdot \nabla b_{j/2} \rangle, & j \text{ even,} \end{cases} \quad j \ge 2. \quad (3.8) $$

This recurrent procedure is illustrated in Appendix A for $j=3$ and $j=4$.

The result given in (3.8) has important consequences. First, it shows that the $j$th-order correction tensor is obtained from the solution of the closure problems from order 0 to $(j-1)/2$ if $j$ is odd or $j/2$ if $j$ is even, instead of the $j+1$ closure problems from order 0 to $j$, as could be inferred from the expanded closure scheme (see (2.8d)). In other words, the solution of the first $N$ closure problems provides all the correction tensors, $K_j$, for $j$ up to $2N-1$, i.e. $\hat{K}_{2N-1}$. This represents a considerable simplification since the computational requirements are divided by a factor of roughly two for the same targeted order of the power-series expansion. Second, it shows that all the correction tensors are symmetric and, moreover, that the odd-order correction tensors are semi-definite positive whereas the even-order ones are semi-definite negative. This further proves that the correction tensors form an alternate series. Furthermore, the fact that $K_2$ is negative indicates that the diagonal terms of $K$ are concave functions of $\alpha$ near $\alpha=0$. All these conclusions are similar to those reached in Lasseux et al. (2023) for slip flow in porous media.
(which physically corresponds to flow under no-slip conditions) contains all the physical information that is necessary to predict both the intrinsic transmissivity and first-order correction tensors. Second, it readily proves that $K_1$ is a symmetric, semi-definite positive tensor. Therefore, if a first-order correction is considered reasonable, which is the case when the Knudsen number is small enough (see some examples in § 4 of Zaouter et al. 2018), only the solution of (2.7) is required and $\hat{K}_1$ (see (2.6)) is symmetric, semi-definite positive. Consequently, it admits an inverse. These conclusions are similar to those obtained for slip flow in porous media analysed in Lasseux, Zaouter & Valdés-Parada (2023).

At this point, it is of interest to investigate the slip correction tensors of higher orders.

### 3.4. Alternative expression of $K_j$, $j \ge 2$

A simpler expression of $K_j$ compared with (2.8d), $j \ge 2$, can be obtained by first considering (3.2), taking $\mathbf{A} = k_0(\nabla\mathbf{b}_0 + \mathbf{I})$ (see (2.7)) and $\mathbf{a} = \mathbf{b}_j$, which gives $\langle k_0\nabla\mathbf{b}_j \rangle = -\langle k_0\nabla\mathbf{b}_0^T \cdot \nabla\mathbf{b}_j \rangle$. Similarly, taking $\mathbf{A} = k_0\nabla\mathbf{b}_1 + k_1(\nabla\mathbf{b}_0 + \mathbf{I})$ (see (2.8) with $j = 1$) and $\mathbf{a} = \mathbf{b}_{j-1}$ yields $\langle k_1\nabla\mathbf{b}_{j-1} \rangle = -\langle (k_0\nabla\mathbf{b}_1^T + k_1\nabla\mathbf{b}_0^T) \cdot \nabla\mathbf{b}_{j-1} \rangle$. When these two results are substituted back into the expression of $K_j$ given in (2.8d) ($j > 1$), the following alternative expression of this tensor is obtained:

$$ K_j = -\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_{j-1} \rangle. \qquad (3.6) $$

To arrive at this result, the fact that $\langle \nabla\mathbf{b}_0^T \cdot (k_0\nabla\mathbf{b}_j + k_1\nabla\mathbf{b}_{j-1}) \rangle = \mathbf{0}$ was employed. This follows from an additional use of (3.2) with $\mathbf{A} = k_0\nabla\mathbf{b}_j + k_1\nabla\mathbf{b}_{j-1}$ (see (2.8)) and $\mathbf{a} = \mathbf{b}_0$ and the application of a transpose.

If $j = 2$, no further step is required as (3.6) gives

$$ K_2 = -\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_1 \rangle. \qquad (3.7) $$

If $j \ge 2$, the next step, taking (3.6) as the initial stage, is to further employ the integral formula (3.2) repeatedly taking $\mathbf{A} = k_0\nabla\mathbf{b}_{j-n} + k_1\nabla\mathbf{b}_{j-n-1}$ and $\mathbf{a} = \mathbf{b}_n$ with $n$ increasing from 1 to $(j-1)/2$ if $j$ is odd or from 1 to $j/2$ if $j$ is even. After substituting the ensuing relationships into the successive expressions of $K_j$ at each value of $n$, this finally yields (note that this procedure is indeed also valid for $j = 2$)

$$ K_j = \begin{cases} \langle k_1 \nabla \mathbf{b}_{(j-1)/2}^T \cdot \nabla \mathbf{b}_{(j-1)/2} \rangle, & j \text{ odd,} \\ -\langle k_0 \nabla \mathbf{b}_{j/2}^T \cdot \nabla \mathbf{b}_{j/2} \rangle, & j \text{ even,} \end{cases} \quad j \ge 2. \qquad (3.8) $$

This recurrent procedure is illustrated in Appendix A for $j = 3$ and $j = 4$.

The result given in (3.8) has important consequences. First, it shows that the $j$th-order correction tensor is obtained from the solution of the closure problems from order 0 to $(j-1)/2$ if $j$ is odd or $j/2$ if $j$ is even, instead of the $j+1$ closure problems from order 0 to $j$, as could be inferred from the expanded closure scheme (see (2.8d)). In other words, the solution of the first $N$ closure problems provides all the correction tensors, $K_j$, for $j$ up to $2N-1$, i.e. $\hat{K}_{2N-1}$. This represents a considerable simplification since the computational requirements are divided by a factor of roughly two for the same targeted order of the power-series expansion. Second, it shows that all the correction tensors are symmetric and, moreover, that the odd-order correction tensors are semi-definite positive whereas the even-order ones are semi-definite negative. This further proves that the correction tensors form an alternate series. Furthermore, the fact that $K_2$ is negative indicates that the diagonal terms of $K$ are concave functions of $\alpha$ near $\alpha = 0$. All these conclusions are similar to those reached in Lasseux et al. (2023) for slip flow in porous media.
3.5. Padé approximant

Since the asymptotic behaviour of **K** in the limit $\xi \overline{Kn} \to +\infty$ is known and the power-series expansion in $\xi \overline{Kn}$ is alternate, another use of this expansion is proposed to construct a Padé approximant, $\tilde{\mathbf{K}}_{(n,d)}$, which is given by (Baker & Graves-Morris 1996)

$$ \tilde{\mathbf{K}}_{(n,d)} = \left( \sum_{i=0}^{n} \alpha^i \mathbf{G}_i \right) \cdot \left( \sum_{j=0}^{d} \alpha^j \mathbf{F}_j \right)^{-1}, \qquad (3.9) $$

where $\mathbf{G}_i$ and $\mathbf{F}_j$ are two series of second-order tensors that will be identified later on. Note that (3.9) corresponds to the ‘right’ approximant in the sense that the inverse term is placed at the right of the expression. Equivalently, a ‘left’ approximant can be constructed and it can be shown that the two are in fact equal (Zinn-Justin 1971).

To determine the values of $n$ and $d$, an asymptotic analysis of the closure problem (2.4) can be carried out. The limit $\alpha \to +\infty$, shows that $\mathbf{K} \sim \alpha \mathbf{K}'_1$ (i.e. $\mathbf{K}$ scales linearly with $\alpha$ at sufficiently large values of this parameter). The intrinsic tensor $\mathbf{K}'_1$ providing the asymptotic slope is given by the solution of the following intrinsic closure problem for $\mathbf{b}'$

$$ \nabla \cdot (k_1(\nabla \mathbf{b}' + \mathbf{I})) = \mathbf{0}, \quad \text{in } \mathcal{A}_{\beta}, \qquad (3.10a) $$

$$ \mathbf{n} \cdot k_1(\nabla \mathbf{b}' + \mathbf{I}) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (3.10b) $$

$$ \mathbf{b}'(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}'(\mathbf{r}), \quad i = 1, 2, \qquad (3.10c) $$

$$ \mathbf{K}'_1 = \langle k_1(\nabla \mathbf{b}' + \mathbf{I}) \rangle. \qquad (3.10d) $$

Consequently, a Padé approximant also exhibiting a linear asymptotic behaviour in the limit $\alpha \to +\infty$ shall be preferentially chosen, hence suggesting $n = d + 1$ in (3.9). For the sake of brevity in presentation, the developments are carried out for the case $d = 1$, which provides the simplest Padé approximant $\tilde{\mathbf{K}}_{(2,1)}$. So as to determine the tensorial coefficients $\mathbf{G}_i$, $i = 0, \dots, 2$ and $\mathbf{F}_j$, $j = 0, 1$, (3.9) is equated to $\hat{\mathbf{K}}_3$, which gives

$$ (\mathbf{G}_0 + \alpha \mathbf{G}_1 + \alpha^2 \mathbf{G}_2) - (\mathbf{K}_0 + \alpha \mathbf{K}_1 + \alpha^2 \mathbf{K}_2 + \alpha^3 \mathbf{K}_3) \cdot (\mathbf{F}_0 + \alpha \mathbf{F}_1) = \mathbf{O}(\alpha^4). \qquad (3.11) $$

Identification of the different terms at the successive orders in $\alpha$ leads to

$$ \mathbf{G}_0 = \mathbf{K}_0 \cdot \mathbf{F}_0, \qquad (3.12a) $$

$$ \mathbf{G}_1 = \mathbf{K}_1 \cdot \mathbf{F}_0 + \mathbf{K}_0 \cdot \mathbf{F}_1, \qquad (3.12b) $$

$$ \mathbf{G}_2 = \mathbf{K}_2 \cdot \mathbf{F}_0 + \mathbf{K}_1 \cdot \mathbf{F}_1, \qquad (3.12c) $$

$$ \mathbf{0} = \mathbf{K}_3 \cdot \mathbf{F}_0 + \mathbf{K}_2 \cdot \mathbf{F}_1. \qquad (3.12d) $$

Without loss of generality, $\mathbf{F}_0$ can be chosen as $\mathbf{F}_0 = \mathbf{I}$, and this allows one to solve the preceding system of equations to obtain

$$ \mathbf{G}_0 = \mathbf{K}_0, \qquad (3.13a) $$

$$ \mathbf{G}_1 = \mathbf{K}_1 - \mathbf{K}_0 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3, \qquad (3.13b) $$

$$ \mathbf{G}_2 = \mathbf{K}_2 - \mathbf{K}_1 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3, \qquad (3.13c) $$

$$ \mathbf{F}_1 = -\mathbf{K}_2^{-1} \cdot \mathbf{K}_3. \qquad (3.13d) $$
### 3.5. Padé approximant

Since the asymptotic behaviour of $K$ in the limit $\xi \bar{K}n \to +\infty$ is known and the power-series expansion in $\xi \bar{K}n$ is alternate, another use of this expansion is proposed to construct a Padé approximant, $\tilde{K}_{(n,d)}$, which is given by (Baker & Graves-Morris 1996)

$$ \tilde{K}_{(n,d)} = \left( \sum_{i=0}^{n} \alpha^i \mathbf{G}_i \right) \cdot \left( \sum_{j=0}^{d} \alpha^j \mathbf{F}_j \right)^{-1}, \quad (3.9) $$

where $\mathbf{G}_i$ and $\mathbf{F}_j$ are two series of second-order tensors that will be identified later on. Note that (3.9) corresponds to the 'right' approximant in the sense that the inverse term is placed at the right of the expression. Equivalently, a 'left' approximant can be constructed and it can be shown that the two are in fact equal (Zinn-Justin 1971).

To determine the values of $n$ and $d$, an asymptotic analysis of the closure problem (2.4) can be carried out. The limit $\alpha \to +\infty$, shows that $K \sim \alpha K'_1$ (i.e. $K$ scales linearly with $\alpha$ at sufficiently large values of this parameter). The intrinsic tensor $K'_1$ providing the asymptotic slope is given by the solution of the following intrinsic closure problem for $b'$

$$ \nabla \cdot (k_1(\nabla b' + \mathbf{I})) = \mathbf{0}, \quad \text{in } \mathcal{A}_{\beta}, \qquad (3.10a) $$

$$ \mathbf{n} \cdot k_1(\nabla b' + \mathbf{I}) = \mathbf{0}, \quad \text{at } \mathcal{C}_{\beta\sigma}, \qquad (3.10b) $$

$$ \mathbf{b}'(\mathbf{r} + \mathbf{l}_i) = \mathbf{b}'(\mathbf{r}), \quad i = 1, 2, \qquad (3.10c) $$

$$ K'_1 = \langle k_1(\nabla \mathbf{b}' + \mathbf{I}) \rangle. \qquad (3.10d) $$

Consequently, a Padé approximant also exhibiting a linear asymptotic behaviour in the limit $\alpha \to +\infty$ shall be preferentially chosen, hence suggesting $n = d + 1$ in (3.9). For the sake of brevity in presentation, the developments are carried out for the case $d = 1$, which provides the simplest Padé approximant $\tilde{K}_{(2,1)}$. So as to determine the tensorial coefficients $\mathbf{G}_i$, $i = 0, \dots, 2$ and $\mathbf{F}_j$, $j = 0, 1$, (3.9) is equated to $\hat{\mathbf{K}}_3$, which gives

$$ (\mathbf{G}_0 + \alpha \mathbf{G}_1 + \alpha^2 \mathbf{G}_2) - (\mathbf{K}_0 + \alpha \mathbf{K}_1 + \alpha^2 \mathbf{K}_2 + \alpha^3 \mathbf{K}_3) \cdot (\mathbf{F}_0 + \alpha \mathbf{F}_1) = \mathbf{O}(\alpha^4). \quad (3.11) $$

Identification of the different terms at the successive orders in $\alpha$ leads to

$$ \mathbf{G}_0 = \mathbf{K}_0 \cdot \mathbf{F}_0, \qquad (3.12a) $$

$$ \mathbf{G}_1 = \mathbf{K}_1 \cdot \mathbf{F}_0 + \mathbf{K}_0 \cdot \mathbf{F}_1, \qquad (3.12b) $$

$$ \mathbf{G}_2 = \mathbf{K}_2 \cdot \mathbf{F}_0 + \mathbf{K}_1 \cdot \mathbf{F}_1, \qquad (3.12c) $$

$$ \mathbf{0} = \mathbf{K}_3 \cdot \mathbf{F}_0 + \mathbf{K}_2 \cdot \mathbf{F}_1. \qquad (3.12d) $$

Without loss of generality, $\mathbf{F}_0$ can be chosen as $\mathbf{F}_0 = \mathbf{I}$, and this allows one to solve the preceding system of equations to obtain

$$ \mathbf{G}_0 = \mathbf{K}_0, \qquad (3.13a) $$

$$ \mathbf{G}_1 = \mathbf{K}_1 - \mathbf{K}_0 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3, \qquad (3.13b) $$

$$ \mathbf{G}_2 = \mathbf{K}_2 - \mathbf{K}_1 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3, \qquad (3.13c) $$

$$ \mathbf{F}_1 = -\mathbf{K}_2^{-1} \cdot \mathbf{K}_3. \qquad (3.13d) $$
<figure><img src="image_4.png" /><figcaption>Figure 2. Unit cell aperture field of an anisotropic Gaussian fracture. White areas represent contact spots of zero aperture.</figcaption></figure>

This yields the final result

$$ \tilde{\mathbf{K}}_{(2,1)} = (\mathbf{K}_0 + \alpha(\mathbf{K}_1 - \mathbf{K}_0 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3) + \alpha^2(\mathbf{K}_2 - \mathbf{K}_1 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3)) \cdot (I - \alpha \mathbf{K}_2^{-1} \cdot \mathbf{K}_3)^{-1}. \quad (3.14) $$

It can be shown that the Padé approximant tensor $\tilde{\mathbf{K}}_{(2,1)}$ given by (3.14) is symmetric. This follows from the symmetry of all the tensors $\mathbf{K}_j, j \ge 0$ and the fact that the transpose of $\tilde{\mathbf{K}}_{(2,1)}$ would be equal to the ‘left’ Padé approximant, which itself is equal to the ‘right’ approximant $\tilde{\mathbf{K}}_{(2,1)}$ (Zinn-Justin 1971). In addition, since the expression of $\tilde{\mathbf{K}}_{(2,1)}$ solely involves the tensors $\mathbf{K}_j$ up to $j = 3$, this approximant can be obtained from the solution of only the first two closure problems, as shown in the previous section. The performance of this approach is assessed by means of numerical simulations in what follows.

## 4. Illustrative example

In this section, a numerical example illustrating the dependence of the transmissivity upon the average Knudsen number is, reported, considering synthetic random rough fractures. The asymptotic expansion of the transmissivity given by (2.6), truncated up to the third order, as well as the Padé approximant given by (3.14), are investigated as well.

The first fracture under consideration was obtained by generating an anisotropic Gaussian shortly correlated surface employing the procedure described in Bergström (2012), with a discretisation of 1025 × 1025 points. A size $\ell_x = \ell_y = \ell$ for the unit cell, correlation lengths $\sigma_{cx} = 0.04\ell$ and $\sigma_{cy} = 0.02\ell$ (in the x- and y-directions, respectively) and a root mean square roughness $R_q = 10^{-3}\ell$ were used. The surface was then brought into contact with a flat plane by vertical overlap, and the aperture was defined as the height difference between the two surfaces. When this difference was negative, the aperture was set to zero, which produces contact zones. The relative position of the two surfaces was chosen to obtain a contact area equal to 7% of the total area. The resulting periodic unit cell aperture field is represented in figure 2.

To begin with, it is of interest to investigate the local flow fields for the particular unit cell configuration reported in figure 2 and make comparisons between the fields of the
<figure><img src="image_4.png" /><figcaption>Figure 2. Unit cell aperture field of an anisotropic Gaussian fracture. White areas represent contact spots of zero aperture.</figcaption></figure>

This yields the final result

$$ \tilde{\mathbf{K}}_{(2,1)} = (\mathbf{K}_0 + \alpha(\mathbf{K}_1 - \mathbf{K}_0 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3) + \alpha^2(\mathbf{K}_2 - \mathbf{K}_1 \cdot \mathbf{K}_2^{-1} \cdot \mathbf{K}_3)) \cdot (I - \alpha \mathbf{K}_2^{-1} \cdot \mathbf{K}_3)^{-1}. \quad (3.14) $$

It can be shown that the Padé approximant tensor $\tilde{\mathbf{K}}_{(2,1)}$ given by (3.14) is symmetric. This follows from the symmetry of all the tensors $\mathbf{K}_j, j \ge 0$ and the fact that the transpose of $\tilde{\mathbf{K}}_{(2,1)}$ would be equal to the ‘left’ Padé approximant, which itself is equal to the ‘right’ approximant $\tilde{\mathbf{K}}_{(2,1)}$ (Zinn-Justin 1971). In addition, since the expression of $\tilde{\mathbf{K}}_{(2,1)}$ solely involves the tensors $\mathbf{K}_j$ up to $j = 3$, this approximant can be obtained from the solution of only the first two closure problems, as shown in the previous section. The performance of this approach is assessed by means of numerical simulations in what follows.

## 4. Illustrative example

In this section, a numerical example illustrating the dependence of the transmissivity upon the average Knudsen number is, reported, considering synthetic random rough fractures. The asymptotic expansion of the transmissivity given by (2.6), truncated up to the third order, as well as the Padé approximant given by (3.14), are investigated as well.

The first fracture under consideration was obtained by generating an anisotropic Gaussian shortly correlated surface employing the procedure described in Bergström (2012), with a discretisation of 1025 × 1025 points. A size $\ell_x = \ell_y = \ell$ for the unit cell, correlation lengths $\sigma_{cx} = 0.04\ell$ and $\sigma_{cy} = 0.02\ell$ (in the x- and y-directions, respectively) and a root mean square roughness $R_q = 10^{-3}\ell$ were used. The surface was then brought into contact with a flat plane by vertical overlap, and the aperture was defined as the height difference between the two surfaces. When this difference was negative, the aperture was set to zero, which produces contact zones. The relative position of the two surfaces was chosen to obtain a contact area equal to 7% of the total area. The resulting periodic unit cell aperture field is represented in figure 2.

To begin with, it is of interest to investigate the local flow fields for the particular unit cell configuration reported in figure 2 and make comparisons between the fields of the
<figure><img src="image_5.png" /><figcaption>Figure 3. Magnitudes of the microscale flux fields (a) $q$ and (b) $\hat{q}_1$ normalised by $\|\langle q \rangle\|$. The corresponding flow streamlines are superimposed in green colour. The results were obtained on the periodic unit cell of figure 2, taking $\xi \bar{K} \bar{n}_x = 0.15$ and $\xi \bar{K} \bar{n}_y = 0.19$ and a pressure gradient along $-e_x$.</figcaption></figure>

microscale flux, $q$, and its first-order expansion, $\hat{q}_1$. To this end, it is necessary to derive an expression for $q$ and $\hat{q}_1$. This is accomplished by using Gray's decomposition $p = \tilde{p} + \langle p \rangle^\beta$ (Gray 1975) and the representation of $\tilde{p}$ in terms of the macroscopic pressure gradient applied on the unit cell according to the closure procedure reported in § 3.2 in Zaouter et al. (2018) in (2.1b). For $q$, the closure for $\tilde{p}$ reads $\tilde{p} = b \cdot \nabla \langle p \rangle^\beta$ and this leads to

$$q = -\frac{\rho}{\mu}k(l + \nabla b) \cdot \nabla \langle p \rangle^{\beta}. \quad (4.1a)$$

For $\hat{q}_1$, $b$ is expanded at the first order to yield $\tilde{p} = (b_0 + \alpha b_1) \cdot \nabla \langle p \rangle^\beta$. This leads to

$$\hat{q}_1 = -\frac{\rho}{\mu}[k_0(l + \nabla b_0) + \alpha(k_1(l + \nabla b_0) + k_0\nabla b_1)] \cdot \nabla \langle p \rangle^{\beta}. \quad (4.1b)$$

In the above equations, $b$, $b_0$ and $b_1$ result from the solutions of (2.4), (2.7) and (2.8) (with $j = 1$), respectively. They are carried out with the same finite volume numerical method as that employed in Zaouter et al. (2018), that is second-order accurate in space. For the sake of brevity, numerical simulations are restricted to the closure problems solution up to the first order so that the comparison between only $q$ and $\hat{q}_1$ is reported.

In figure 3, the magnitudes and streamlines of the microscale flux fields, $q$ and $\hat{q}_1$ (normalised by $\|\langle q \rangle\|$), imposing incompressible flow, are reported for $\xi \bar{K} \bar{n}_x = 0.15$ and $\xi \bar{K} \bar{n}_y = 0.19$. For the sake of simplicity, the macroscopic pressure gradient is taken to be along $-e_x$. As can be observed in this figure, the flux magnitudes do not seem to exhibit significant differences. However, this can be misleading since the flux magnitudes are strongly influenced by the $x$-component values, which may not markedly differ between themselves, and this may damp the differences between the $y$-components. A much more significant mismatch is noticeable in the streamline patterns. These remarks suggest that significant differences may be observed between $K$ and $\hat{K}_1$ as will be shown later. For the moment, it is pertinent to gain a more quantitative idea about the differences between the microscale flux predictions. In figure 4 the fields of the relative differences, defined as

$$\epsilon_i = \frac{\hat{q}_{1i} - q_i}{\langle q_i \rangle}, \quad i = x, y, \quad (4.2)$$

are reported for both components of the flux, taking the same conditions as those considered in figure 3. The results show that the largest differences are found for the
<figure><img src="image_5.png" /><figcaption>Figure 3. Magnitudes of the microscale flux fields (a) $q$ and (b) $\hat{q}_1$ normalised by $\|\langle q \rangle\|$. The corresponding flow streamlines are superimposed in green colour. The results were obtained on the periodic unit cell of figure 2, taking $\xi \bar{K} n_x = 0.15$ and $\xi \bar{K} n_y = 0.19$ and a pressure gradient along $-e_x$.</figcaption></figure>

microscale flux, $q$, and its first-order expansion, $\hat{q}_1$. To this end, it is necessary to derive an expression for $q$ and $\hat{q}_1$. This is accomplished by using Gray's decomposition $p = \tilde{p} + \langle p \rangle^\beta$ (Gray 1975) and the representation of $\tilde{p}$ in terms of the macroscopic pressure gradient applied on the unit cell according to the closure procedure reported in § 3.2 in Zaouter et al. (2018) in (2.1b). For $q$, the closure for $\tilde{p}$ reads $\tilde{p} = b \cdot \nabla \langle p \rangle^\beta$ and this leads to

$$q = -\frac{\rho}{\mu}k(l + \nabla b) \cdot \nabla \langle p \rangle^{\beta}. \quad (4.1a)$$

For $\hat{q}_1$, $b$ is expanded at the first order to yield $\tilde{p} = (b_0 + \alpha b_1) \cdot \nabla \langle p \rangle^\beta$. This leads to

$$\hat{q}_1 = -\frac{\rho}{\mu}[k_0(l + \nabla b_0) + \alpha(k_1(l + \nabla b_0) + k_0\nabla b_1)] \cdot \nabla \langle p \rangle^{\beta}. \quad (4.1b)$$

In the above equations, $b$, $b_0$ and $b_1$ result from the solutions of (2.4), (2.7) and (2.8) (with $j=1$), respectively. They are carried out with the same finite volume numerical method as that employed in Zaouter et al. (2018), that is second-order accurate in space. For the sake of brevity, numerical simulations are restricted to the closure problems solution up to the first order so that the comparison between only $q$ and $\hat{q}_1$ is reported.

In figure 3, the magnitudes and streamlines of the microscale flux fields, $q$ and $\hat{q}_1$ (normalised by $\|\langle q \rangle\|$), imposing incompressible flow, are reported for $\xi \bar{K} n_x = 0.15$ and $\xi \bar{K} n_y = 0.19$. For the sake of simplicity, the macroscopic pressure gradient is taken to be along $-e_x$. As can be observed in this figure, the flux magnitudes do not seem to exhibit significant differences. However, this can be misleading since the flux magnitudes are strongly influenced by the $x$-component values, which may not markedly differ between themselves, and this may damp the differences between the $y$-components. A much more significant mismatch is noticeable in the streamline patterns. These remarks suggest that significant differences may be observed between $K$ and $\hat{K}_1$ as will be shown later. For the moment, it is pertinent to gain a more quantitative idea about the differences between the microscale flux predictions. In figure 4 the fields of the relative differences, defined as

$$\epsilon_i = \frac{\hat{q}_{1i} - q_i}{\langle q_i \rangle}, \quad i = x, y, \quad (4.2)$$

are reported for both components of the flux, taking the same conditions as those considered in figure 3. The results show that the largest differences are found for the
<figure><img src="image_6.png" /><figcaption>Figure 4. Differences between the (a) x and (b) y components of the microscale flux $q$ and $\hat{q}_1$ reported in figure 3, respectively normalised by the corresponding averages $\langle q_i \rangle$ ($i = x, y$).</figcaption></figure>

<table><caption>Table 1. Dimensionless components of the intrinsic transmissivity tensor, $K_0$, and slip correction tensor, $K_1$, computed from the conventional method ((2.7d) and (2.8d) with $j = 1$, respectively) and the improved method ((3.4) and (3.5), respectively). Here, $\Delta$ is the relative difference between the values obtained from the conventional and improved methods, taking the former as the reference.</caption><thead><tr><th rowspan="2">Component</th><th colspan="3"><i>K</i><sub>0</sub>/<i>ℓ</i><sup>3</sup></th><th colspan="3"><i>K</i><sub>1</sub>/<i>ℓ</i><sup>2</sup></th></tr><tr><th>Conventional</th><th>Improved</th><th>Δ (%)</th><th>Conventional</th><th>Improved</th><th>Δ (%)</th></tr></thead><tbody><tr><td>xx</td><td>2.533 × 10<sup>-10</sup></td><td>2.532 × 10<sup>-10</sup></td><td>0.049</td><td>1.741 × 10<sup>-7</sup></td><td>1.740 × 10<sup>-7</sup></td><td>0.102</td></tr><tr><td>yy</td><td>1.149 × 10<sup>-10</sup></td><td>1.149 × 10<sup>-10</sup></td><td>0.056</td><td>9.661 × 10<sup>-8</sup></td><td>9.651 × 10<sup>-8</sup></td><td>0.106</td></tr><tr><td>xy</td><td>1.169 × 10<sup>-12</sup></td><td>1.160 × 10<sup>-12</sup></td><td>0.821</td><td>2.356 × 10<sup>-9</sup></td><td>2.335 × 10<sup>-9</sup></td><td>0.854</td></tr><tr><td>yx</td><td>1.169 × 10<sup>-12</sup></td><td>1.160 × 10<sup>-12</sup></td><td>0.821</td><td>2.356 × 10<sup>-9</sup></td><td>2.335 × 10<sup>-9</sup></td><td>0.854</td></tr></tbody></table>

predictions of the y-component of the flux, especially in throats formed between close solid contact patches. This suggests that the first-order predictions of $K$ may be less accurate for the off-diagonal terms than for the diagonal ones, as will be further explored.

With the solutions for $b_0$ and $b_1$ at hand, $K_0$ and $K_1$ can be computed either from (2.7d) and (2.8d) ($j = 1$), respectively, (this is referred to as the conventional method hereafter) or from (3.4) and (3.5), respectively, (a procedure called the improved method). The components of $K_0$ and $K_1$, computed with the two different methods, are reported in table 1. These results confirm the symmetry of these two tensors and show an excellent agreement for all the components, with a maximum relative difference, $\Delta$, around 0.8% with respect to the conventional method. For $K_0$, this difference only originates from signatures of the numerical inaccuracies on $\nabla b_0$ between formulas (2.7d) and (3.4). Note that, in practice, (2.7d) should be preferred to (3.4) to compute $K_0$ as it is expected to be, however, more accurate. For $K_1$, a supplementary effect comes from the fact that the solution for $b_1$ is tainted by the inaccuracies on $\nabla b_0$ when the conventional method is employed. In any case, the differences remain small and decrease as the computational mesh size is decreased (i.e. for sufficiently converged solutions). Such small error appears more than acceptable compared with the computational burden of having to solve multiple closure problems up to the desired order. Nevertheless, it is difficult to assess which method is more accurate to compute $K_j$ for $j \ge 1$. These results validate the approach presented in this work.
<figure><img src="image_6.png" /><figcaption>Figure 4. Differences between the (a) x and (b) y components of the microscale flux $q$ and $\hat{q}_1$ reported in figure 3, respectively normalised by the corresponding averages $\langle q_i \rangle$ ($i = x, y$).</figcaption></figure>

<table><caption>Table 1. Dimensionless components of the intrinsic transmissivity tensor, $K_0$, and slip correction tensor, $K_1$, computed from the conventional method ((2.7d) and (2.8d) with $j = 1$, respectively) and the improved method ((3.4) and (3.5), respectively). Here, $\Delta$ is the relative difference between the values obtained from the conventional and improved methods, taking the former as the reference.</caption><thead><tr><th rowspan="2">Component</th><th colspan="3"><i>K</i><sub>0</sub>/<i>ℓ</i><sup>3</sup></th><th colspan="3"><i>K</i><sub>1</sub>/<i>ℓ</i><sup>2</sup></th></tr><tr><th>Conventional</th><th>Improved</th><th>Δ (%)</th><th>Conventional</th><th>Improved</th><th>Δ (%)</th></tr></thead><tbody><tr><td>xx</td><td>2.533 × 10<sup>-10</sup></td><td>2.532 × 10<sup>-10</sup></td><td>0.049</td><td>1.741 × 10<sup>-7</sup></td><td>1.740 × 10<sup>-7</sup></td><td>0.102</td></tr><tr><td>yy</td><td>1.149 × 10<sup>-10</sup></td><td>1.149 × 10<sup>-10</sup></td><td>0.056</td><td>9.661 × 10<sup>-8</sup></td><td>9.651 × 10<sup>-8</sup></td><td>0.106</td></tr><tr><td>xy</td><td>1.169 × 10<sup>-12</sup></td><td>1.160 × 10<sup>-12</sup></td><td>0.821</td><td>2.356 × 10<sup>-9</sup></td><td>2.335 × 10<sup>-9</sup></td><td>0.854</td></tr><tr><td>yx</td><td>1.169 × 10<sup>-12</sup></td><td>1.160 × 10<sup>-12</sup></td><td>0.821</td><td>2.356 × 10<sup>-9</sup></td><td>2.335 × 10<sup>-9</sup></td><td>0.854</td></tr></tbody></table>

predictions of the y-component of the flux, especially in throats formed between close solid contact patches. This suggests that the first-order predictions of $K$ may be less accurate for the off-diagonal terms than for the diagonal ones, as will be further explored.

With the solutions for $b_0$ and $b_1$ at hand, $K_0$ and $K_1$ can be computed either from (2.7d) and (2.8d) ($j = 1$), respectively, (this is referred to as the conventional method hereafter) or from (3.4) and (3.5), respectively, (a procedure called the improved method). The components of $K_0$ and $K_1$, computed with the two different methods, are reported in table 1. These results confirm the symmetry of these two tensors and show an excellent agreement for all the components, with a maximum relative difference, $\Delta$, around 0.8% with respect to the conventional method. For $K_0$, this difference only originates from signatures of the numerical inaccuracies on $\nabla b_0$ between formulas (2.7d) and (3.4). Note that, in practice, (2.7d) should be preferred to (3.4) to compute $K_0$ as it is expected to be, however, more accurate. For $K_1$, a supplementary effect comes from the fact that the solution for $b_1$ is tainted by the inaccuracies on $\nabla b_0$ when the conventional method is employed. In any case, the differences remain small and decrease as the computational mesh size is decreased (i.e. for sufficiently converged solutions). Such small error appears more than acceptable compared with the computational burden of having to solve multiple closure problems up to the desired order. Nevertheless, it is difficult to assess which method is more accurate to compute $K_j$ for $j \ge 1$. These results validate the approach presented in this work.
Effective transmissivity for slip flow in a fracture

<figure><img src="image_7.png" /><figcaption>Figure 5. (a) The xx-components, (b) yy-components, (c) xy-components of the effective transmissivity tensor, **K**, of the anisotropic Gaussian fracture, sketched in figure 2, normalised by the corresponding components of the intrinsic transmissivity tensor, **K**₀ (open symbols), as a function of the Knudsen number; components of the power-series expansions, $\tilde{K}_{mij}$ (for *m* up to 3, see (2.6)), and the Padé approximant, $\tilde{K}_{(2,1)ij}$, *i*, *j* = *x*, *y* (lines). (d–f) Corresponding relative errors.</figcaption></figure>

The focus is now laid upon the series expansions of the transmissivity tensor. From the zeroth- and first-order closure problem solutions, the correction tensors $K_j$ in (2.6) can be computed up to $j = 3$ with the improved method. The analysis is carried out in terms of the Knudsen number, which is considered to be dependent on the flow direction following the definition proposed in Zaouter *et al*. (2018), i.e. $\overline{Kn}_i = \bar{\lambda}/h_{\beta i}$, $i = x, y$, where $h_{\beta i}$ is estimated as $h_{\beta i} = \sqrt[3]{12K_{0ii}}$. In this way, $h_{\beta i}$ represents the aperture of a fracture made of two plane parallel plates that would exhibit the same intrinsic transmissivity $K_{0ii}$ in the *i*-th-direction. In agreement with the slip flow regime, results are reported for $\xi \overline{Kn}_i \le 0.15$, with the idea that $\xi \approx 1$. With these estimates of $h_{\beta i}$, and taking the correlation length, $\sigma_{ci}$, as a measure of $\ell_{\beta i}$, this leads to $\varepsilon \sim 0.036$ in the *x*-direction and $\varepsilon \sim 0.055$ in the *y*-direction, confirming that the Reynolds approximation, which is the basis of the model employed here, is perfectly justified.

In figure 5, the components of the transmissivity tensors, normalised by the corresponding intrinsic transmissivity component, are represented as a function of $\xi \overline{Kn}_i$, $i = x, y$. For $K_{xy}/K_{0xy}$, the Knudsen number is taken as $\overline{Kn}_x$. Since the tensors are all symmetric, values of only one of the off-diagonal components are reported. From this figure, it can be clearly seen that the components of **K** (open symbols in figure 5) are increasing nonlinear functions (in magnitude) of the Knudsen number. Nonlinearity was already discussed in Zaouter *et al*. (2018) (see end of p. 431). Moreover, the concave dependence of the diagonal terms upon $\xi \overline{Kn}_i$ is confirmed. A similar behaviour is also observed for the off-diagonal terms, the formal proof of which, in the general case, remains yet unclear. Successive expansions of the transmissivity at increasing orders (lines in
Effective transmissivity for slip flow in a fracture

<figure><img src="image_7.png" /><figcaption>Figure 5. (a) The xx-components, (b) yy-components, (c) xy-components of the effective transmissivity tensor, **K**, of the anisotropic Gaussian fracture, sketched in figure 2, normalised by the corresponding components of the intrinsic transmissivity tensor, **K**₀ (open symbols), as a function of the Knudsen number; components of the power-series expansions, $\tilde{K}_{mij}$ (for *m* up to 3, see (2.6)), and the Padé approximant, $\tilde{K}_{(2,1)ij}$, *i*, *j* = *x*, *y* (lines). (d–f) Corresponding relative errors.</figcaption></figure>

The focus is now laid upon the series expansions of the transmissivity tensor. From the zeroth- and first-order closure problem solutions, the correction tensors $K_j$ in (2.6) can be computed up to $j = 3$ with the improved method. The analysis is carried out in terms of the Knudsen number, which is considered to be dependent on the flow direction following the definition proposed in Zaouter *et al*. (2018), i.e. $\overline{Kn}_i = \bar{\lambda}/h_{\beta i}$, $i = x, y$, where $h_{\beta i}$ is estimated as $h_{\beta i} = \sqrt[3]{12K_{0ii}}$. In this way, $h_{\beta i}$ represents the aperture of a fracture made of two plane parallel plates that would exhibit the same intrinsic transmissivity $K_{0ii}$ in the *i*-th-direction. In agreement with the slip flow regime, results are reported for $\xi \overline{Kn}_i \le 0.15$, with the idea that $\xi \approx 1$. With these estimates of $h_{\beta i}$, and taking the correlation length, $\sigma_{ci}$, as a measure of $\ell_{\beta i}$, this leads to $\varepsilon \sim 0.036$ in the *x*-direction and $\varepsilon \sim 0.055$ in the *y*-direction, confirming that the Reynolds approximation, which is the basis of the model employed here, is perfectly justified.

In figure 5, the components of the transmissivity tensors, normalised by the corresponding intrinsic transmissivity component, are represented as a function of $\xi \overline{Kn}_i$, $i = x, y$. For $K_{xy}/K_{0xy}$, the Knudsen number is taken as $\overline{Kn}_x$. Since the tensors are all symmetric, values of only one of the off-diagonal components are reported. From this figure, it can be clearly seen that the components of **K** (open symbols in figure 5) are increasing nonlinear functions (in magnitude) of the Knudsen number. Nonlinearity was already discussed in Zaouter *et al*. (2018) (see end of p. 431). Moreover, the concave dependence of the diagonal terms upon $\xi \overline{Kn}_i$ is confirmed. A similar behaviour is also observed for the off-diagonal terms, the formal proof of which, in the general case, remains yet unclear. Successive expansions of the transmissivity at increasing orders (lines in
<table><thead><tr><th>Component</th><th>K/ℓ³</th><th>K̂₁/ℓ³</th><th>K̃<sub>(2,1)</sub>/ℓ³</th><th>ΔK̂₁ (%)</th><th>ΔK̃<sub>(2,1)</sub> (%)</th></tr></thead><tbody><tr><td>xx</td><td>1.690 × 10<sup>-9</sup></td><td>1.821 × 10<sup>-9</sup></td><td>1.726 × 10<sup>-9</sup></td><td>7.711</td><td>2.104</td></tr><tr><td>yy</td><td>9.047 × 10<sup>-10</sup></td><td>9.845 × 10<sup>-10</sup></td><td>9.303 × 10<sup>-10</sup></td><td>8.823</td><td>2.831</td></tr><tr><td>xy</td><td>1.634 × 10<sup>-11</sup></td><td>2.237 × 10<sup>-11</sup></td><td>1.762 × 10<sup>-11</sup></td><td>36.883</td><td>7.803</td></tr></tbody><tfoot><tr><td colspan="6">Table 2. Comparison of the dimensionless effective transmissivity tensor with K̂₁/ℓ³ and K̃<sub>(2,1)</sub>/ℓ³ for the fracture depicted in figure 2 taking ξK̄nₓ = 1.03 and ξK̄nᵧ = 1.35. Here, ΔK̂₁ = |Kᵢⱼ - K̂₁ᵢⱼ|/Kᵢⱼ, ΔK̃<sub>(2,1)</sub> = |Kᵢⱼ - K̃<sub>(2,1)</sub>|<sub>ij</sub>/Kᵢⱼ, i, j = x, y.</td></tr></tfoot></table>

figure 5) provide greater accuracy at small Knudsen numbers. This is further supported by figure 5(d-f) showing the relative error, $\epsilon_{ij}$, between the ij-component of the effective transmissivity tensor, $K_{ij}$, $i, j = x, y$, and its prediction either by the series expansion, $\hat{K}_{m_{ij}}$, $m = 1, 2, 3$, or the Padé approximant, $\tilde{K}_{(2,1)ij}$, taking $K_{ij}$ as the reference.

The results in figure 5 show that a first-order approximation is only reasonable for sufficiently small values of $\xi\overline{Kn}_i$ (i.e. smaller than 0.05 in the case under study). Note that a third-order approximation contributes to reduce the error in the predictions up until a cross-over value of $\xi\overline{Kn}_i$ that is approximately 0.08. For larger values, a third-order power-series expansion actually decreases the accuracy of the predictions. This behaviour is certainly related to the radius of convergence of the power-series expansion. Furthermore, the Padé approximant clearly outperforms the predictions from the power-series expansion. Indeed, the maximum relative error resulting from this approach is approximately 1% for the off-diagonal term, whereas for the diagonal terms this value decreases by about one order of magnitude. These results are in agreement with the observations made on the flux fields reported in figures 3 and 4. These conclusions are indeed dependent on the unit cell geometry. In Appendix B, the results corresponding to an isotropic Gaussian geometry are reported. In this case, the cross-over value of $\xi\overline{Kn}_i$ is approximately twice the one deduced from figure 5. Nevertheless, the Padé approximant remains in excellent agreement. This confirms the relevance of this approach.

It is worth pointing out that the formalism of rarefied flow used in this work limits the Knudsen number values that can be used. Still, in the context of flow with slip effects induced by other physical mechanisms (for example, an effective slip boundary condition resulting from replacing a rough surface by an effective one), it is admissible to consider much larger values of the dimensionless slip length. Consequently, significantly larger errors in the effective transmissivity predictions from a first-order expansion can be expected. For the purpose of illustration, results for the fracture depicted in figure 2 corresponding to $\xi\overline{Kn}_x = 1.03$ and $\xi\overline{Kn}_y = 1.35$ are reported in table 2, which confirm the expected behaviour and further justify the relevance of the Padé approximant. Note that the relative errors given in this table are expected to remain more or less constant as $\xi\overline{Kn}$ keeps increasing. This is due to the fact that $K$, $\hat{K}_1$ and $\tilde{K}_{(2,1)}$ all behave linearly in the limit $\xi\overline{Kn} \to +\infty$.

## 5. Conclusions

In this work, an efficient method to determine the intrinsic and slip correction transmissivity (second-order) tensors, involved in the macroscopic model for steady, creeping, single-phase, isothermal, Newtonian slip flow in a rough fracture is reported.
<table><thead><tr><th>Component</th><th>K/ℓ³</th><th>K̂₁/ℓ³</th><th>K̃<sub>(2,1)</sub>/ℓ³</th><th>ΔK̂₁ (%)</th><th>ΔK̃<sub>(2,1)</sub> (%)</th></tr></thead><tbody><tr><td>xx</td><td>1.690 × 10<sup>-9</sup></td><td>1.821 × 10<sup>-9</sup></td><td>1.726 × 10<sup>-9</sup></td><td>7.711</td><td>2.104</td></tr><tr><td>yy</td><td>9.047 × 10<sup>-10</sup></td><td>9.845 × 10<sup>-10</sup></td><td>9.303 × 10<sup>-10</sup></td><td>8.823</td><td>2.831</td></tr><tr><td>xy</td><td>1.634 × 10<sup>-11</sup></td><td>2.237 × 10<sup>-11</sup></td><td>1.762 × 10<sup>-11</sup></td><td>36.883</td><td>7.803</td></tr></tbody></table>
Table 2. Comparison of the dimensionless effective transmissivity tensor with K̂₁/ℓ³ and K̃<sub>(2,1)</sub>/ℓ³ for the fracture depicted in figure 2 taking ξK̄nₓ = 1.03 and ξK̄nᵧ = 1.35. Here, ΔK̂₁ = |Kᵢⱼ - K̂₁ᵢⱼ|/Kᵢⱼ, ΔK̃<sub>(2,1)</sub> = |Kᵢⱼ - K̃<sub>(2,1)</sub>|/Kᵢⱼ, i, j = x, y.

figure 5) provide greater accuracy at small Knudsen numbers. This is further supported by figure 5(d-f) showing the relative error, $\epsilon_{ij}$, between the ij-component of the effective transmissivity tensor, $K_{ij}$, $i, j = x, y$, and its prediction either by the series expansion, $\hat{K}_{m_{ij}}$, $m = 1, 2, 3$, or the Padé approximant, $\tilde{K}_{(2,1)ij}$, taking $K_{ij}$ as the reference.

The results in figure 5 show that a first-order approximation is only reasonable for sufficiently small values of $\xi\overline{Kn}_i$ (i.e. smaller than 0.05 in the case under study). Note that a third-order approximation contributes to reduce the error in the predictions up until a cross-over value of $\xi\overline{Kn}_i$ that is approximately 0.08. For larger values, a third-order power-series expansion actually decreases the accuracy of the predictions. This behaviour is certainly related to the radius of convergence of the power-series expansion. Furthermore, the Padé approximant clearly outperforms the predictions from the power-series expansion. Indeed, the maximum relative error resulting from this approach is approximately 1% for the off-diagonal term, whereas for the diagonal terms this value decreases by about one order of magnitude. These results are in agreement with the observations made on the flux fields reported in figures 3 and 4. These conclusions are indeed dependent on the unit cell geometry. In Appendix B, the results corresponding to an isotropic Gaussian geometry are reported. In this case, the cross-over value of $\xi\overline{Kn}_i$ is approximately twice the one deduced from figure 5. Nevertheless, the Padé approximant remains in excellent agreement. This confirms the relevance of this approach.

It is worth pointing out that the formalism of rarefied flow used in this work limits the Knudsen number values that can be used. Still, in the context of flow with slip effects induced by other physical mechanisms (for example, an effective slip boundary condition resulting from replacing a rough surface by an effective one), it is admissible to consider much larger values of the dimensionless slip length. Consequently, significantly larger errors in the effective transmissivity predictions from a first-order expansion can be expected. For the purpose of illustration, results for the fracture depicted in figure 2 corresponding to $\xi\overline{Kn}_x = 1.03$ and $\xi\overline{Kn}_y = 1.35$ are reported in table 2, which confirm the expected behaviour and further justify the relevance of the Padé approximant. Note that the relative errors given in this table are expected to remain more or less constant as $\xi\overline{Kn}$ keeps increasing. This is due to the fact that $K, \hat{K}_1$ and $\tilde{K}_{(2,1)}$ all behave linearly in the limit $\xi\overline{Kn} \to +\infty$.

## 5. Conclusions

In this work, an efficient method to determine the intrinsic and slip correction transmissivity (second-order) tensors, involved in the macroscopic model for steady, creeping, single-phase, isothermal, Newtonian slip flow in a rough fracture is reported.
Effective transmissivity for slip flow in a fracture

<figure><img src="image_8.png" /><figcaption>Figure 6. Unit cell aperture field of an isotropic Gaussian fracture.</figcaption></figure>

These coefficients result from a power-series expansion of the closure problem that yields the effective transmissivity tensor in terms of the dimensionless slip parameter (i.e. the Knudsen number for rarefied gas flow or, more generally, the slip length to characteristic aperture size ratio). They are obtained from the solution of closure problems at the successive orders in the slip parameter within a periodic unit cell representative of the system. The reported method shows that the solution of the first *N* closure problems provides the intrinsic and correction transmissivity tensors up to the order $2N - 1$. In particular, the classical first-order (Klinkenberg-like) correction tensor is obtained from the solution of the same (zeroth-order) closure problem that is required to compute the intrinsic transmissivity tensor (that is the single effective coefficient in the absence of slip effects). This result allows one to conclude that the no-slip flow pore-scale fields (i.e. the solution of the zeroth-order closure problem) also contain all the physical information characterising the first-order slip correction. The formulation allows one to also conclude that the effective and intrinsic transmissivity tensors, together with the odd-order correction tensors, are semi-definite positive, whereas the even-order correction tensors are semi-definite negative, and further, that all the tensors are symmetric. In addition, the power-series expansion was used to construct the simplest Padé approximant for the effective transmissivity tensor, which only requires the solution of the zeroth- and first-order closure problems. The method is validated through numerical examples of slip flow within random rough fractures showing its relevance. In particular, the numerical results show that the Padé approximant outperforms the predictions from the power-series expansion.

**Funding.** T.Z. acknowledges the support from the joint CEA-TECHNETICS Group France maestral sealing laboratory. D.L. thanks the RRI BEST-Usine du futur (factory of the future) of the University of Bordeaux and the GdR HydroGEMM of CNRS for their support. F.J.V.P. thanks ENSAM for the financial support for a visiting professor position.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

D Tony Zaouter https://orcid.org/0000-0003-1330-9130;

D Francisco J. Valdés-Parada https://orcid.org/0000-0003-4208-1075;

D Didier Lasseux https://orcid.org/0000-0002-6080-8226.
Effective transmissivity for slip flow in a fracture

<figure><img src="image_8.png" /><figcaption>Figure 6. Unit cell aperture field of an isotropic Gaussian fracture.</figcaption></figure>

These coefficients result from a power-series expansion of the closure problem that yields the effective transmissivity tensor in terms of the dimensionless slip parameter (i.e. the Knudsen number for rarefied gas flow or, more generally, the slip length to characteristic aperture size ratio). They are obtained from the solution of closure problems at the successive orders in the slip parameter within a periodic unit cell representative of the system. The reported method shows that the solution of the first *N* closure problems provides the intrinsic and correction transmissivity tensors up to the order $2N - 1$. In particular, the classical first-order (Klinkenberg-like) correction tensor is obtained from the solution of the same (zeroth-order) closure problem that is required to compute the intrinsic transmissivity tensor (that is the single effective coefficient in the absence of slip effects). This result allows one to conclude that the no-slip flow pore-scale fields (i.e. the solution of the zeroth-order closure problem) also contain all the physical information characterising the first-order slip correction. The formulation allows one to also conclude that the effective and intrinsic transmissivity tensors, together with the odd-order correction tensors, are semi-definite positive, whereas the even-order correction tensors are semi-definite negative, and further, that all the tensors are symmetric. In addition, the power-series expansion was used to construct the simplest Padé approximant for the effective transmissivity tensor, which only requires the solution of the zeroth- and first-order closure problems. The method is validated through numerical examples of slip flow within random rough fractures showing its relevance. In particular, the numerical results show that the Padé approximant outperforms the predictions from the power-series expansion.

**Funding.** T.Z. acknowledges the support from the joint CEA-TECHNETICS Group France maestral sealing laboratory. D.L. thanks the RRI BEST-Usine du futur (factory of the future) of the University of Bordeaux and the GdR HydroGEMM of CNRS for their support. F.J.V.P. thanks ENSAM for the financial support for a visiting professor position.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.**

D Tony Zaouter https://orcid.org/0000-0003-1330-9130;

D Francisco J. Valdés-Parada https://orcid.org/0000-0003-4208-1075;

D Didier Lasseux https://orcid.org/0000-0002-6080-8226.
<figure><img src="image_9.png" /><figcaption>Figure 7. (a) The xx-components, (b) yy-components, (c) xy-components of the effective transmissivity tensor, **K**, of the isotropic Gaussian fracture, sketched in figure 2, normalised by the corresponding components of the intrinsic transmissivity tensor, **K**₀ (open symbols), as a function of the Knudsen number; components of the power-series expansions, ̂*K*<sub>*mij*</sub> (for *m* up to 3, see (2.6)), and the Padé approximant ̃*K*<sub>(2,1)ij</sub>, *i*, *j* = *x*, *y* (lines). (d–f) Corresponding relative errors.</figcaption></figure>

**Appendix A. Procedure to express K_j, j ≥ 2**

The purpose of this appendix is to illustrate the recurrent procedure leading to the expressions of $K_j$ given in (3.8). This is carried out for $j = 3$ and $j = 4$.

Procedure for $\mathbf{K}_3$ ($j=3$)

In this case, (3.6) gives $\mathbf{K}_3 = -\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle$ and the recurrence procedure is carried out once (for $n=1$) with the identity given in (3.2) taking $\mathbf{A} = k_0 \nabla \mathbf{b}_{j-n} + k_1 \nabla \mathbf{b}_{j-n-1}$ and $\mathbf{a} = \mathbf{b}_n$, i.e. $\mathbf{A} = k_0 \nabla \mathbf{b}_2 + k_1 \nabla \mathbf{b}_1$ and $\mathbf{a} = \mathbf{b}_1$. This yields $\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle = -\langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_1 \rangle$. Upon substitution of this result into the above expression of $\mathbf{K}_3$, it follows that $\mathbf{K}_3 = \langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_1 \rangle$.

Procedure for $\mathbf{K}_4$ ($j=4$)

The starting point is (3.6), which gives $\mathbf{K}_4 = -\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_3 \rangle$, whereas the recurrence is employed with $n=1, 2$. Consequently, the integral formula (3.2) is used twice, first with $\mathbf{A} = k_0 \nabla \mathbf{b}_3 + k_1 \nabla \mathbf{b}_2$ and $\mathbf{a} = \mathbf{b}_1$ ($n=1$), and second, with $\mathbf{A} = k_0 \nabla \mathbf{b}_2 + k_1 \nabla \mathbf{b}_1$ and $\mathbf{a} = \mathbf{b}_2$ ($n=2$). This leads to $\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_3 \rangle = -\langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle$ for $n=1$ and $\langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle = -\langle k_0 \nabla \mathbf{b}_2^T \cdot \nabla \mathbf{b}_2 \rangle$ for $n=2$. Combining these two results leads to $\mathbf{K}_4 = -\langle k_0 \nabla \mathbf{b}_2^T \cdot \nabla \mathbf{b}_2 \rangle$.

**Appendix B. Results for an isotropic Gaussian surface**

In this section, results are presented for the isotropic Gaussian fracture illustrated in
**figure 6**. This geometry was constructed employing the same procedure described in the
<figure><img src="image_9.png" /><figcaption>Figure 7. (a) The xx-components, (b) yy-components, (c) xy-components of the effective transmissivity tensor, **K**, of the isotropic Gaussian fracture, sketched in figure 2, normalised by the corresponding components of the intrinsic transmissivity tensor, **K**₀ (open symbols), as a function of the Knudsen number; components of the power-series expansions, ̂*K*<sub>*mij*</sub> (for *m* up to 3, see (2.6)), and the Padé approximant ̃*K*<sub>(2,1)ij</sub>, *i*, *j* = *x*, *y* (lines). (d–f) Corresponding relative errors.</figcaption></figure>

**Appendix A. Procedure to express K_j, j ≥ 2**

The purpose of this appendix is to illustrate the recurrent procedure leading to the expressions of $K_j$ given in (3.8). This is carried out for $j = 3$ and $j = 4$.

Procedure for $\mathbf{K}_3$ ($j=3$)

In this case, (3.6) gives $\mathbf{K}_3 = -\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle$ and the recurrence procedure is carried out once (for $n=1$) with the identity given in (3.2) taking $\mathbf{A} = k_0 \nabla \mathbf{b}_{j-n} + k_1 \nabla \mathbf{b}_{j-n-1}$ and $\mathbf{a} = \mathbf{b}_n$, i.e. $\mathbf{A} = k_0 \nabla \mathbf{b}_2 + k_1 \nabla \mathbf{b}_1$ and $\mathbf{a} = \mathbf{b}_1$. This yields $\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle = -\langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_1 \rangle$. Upon substitution of this result into the above expression of $\mathbf{K}_3$, it follows that $\mathbf{K}_3 = \langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_1 \rangle$.

Procedure for $\mathbf{K}_4$ ($j=4$)

The starting point is (3.6), which gives $\mathbf{K}_4 = -\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_3 \rangle$, whereas the recurrence
is employed with $n=1, 2$. Consequently, the integral formula (3.2) is used twice, first
with $\mathbf{A} = k_0 \nabla \mathbf{b}_3 + k_1 \nabla \mathbf{b}_2$ and $\mathbf{a} = \mathbf{b}_1$ ($n=1$), and second, with $\mathbf{A} = k_0 \nabla \mathbf{b}_2 + k_1 \nabla \mathbf{b}_1$
and $\mathbf{a} = \mathbf{b}_2$ ($n=2$). This leads to $\langle k_0 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_3 \rangle = -\langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle$ for $n=1$ and
$\langle k_1 \nabla \mathbf{b}_1^T \cdot \nabla \mathbf{b}_2 \rangle = -\langle k_0 \nabla \mathbf{b}_2^T \cdot \nabla \mathbf{b}_2 \rangle$ for $n=2$. Combining these two results leads to $\mathbf{K}_4 =$
$-\langle k_0 \nabla \mathbf{b}_2^T \cdot \nabla \mathbf{b}_2 \rangle$.

**Appendix B. Results for an isotropic Gaussian surface**

In this section, results are presented for the isotropic Gaussian fracture illustrated in
**figure 6**. This geometry was constructed employing the same procedure described in the
main text for the Gaussian anisotropic fracture, albeit the correlation lengths are now $\sigma_{cx} = \sigma_{cy} = 0.05\ell$ maintaining the same root mean square roughness value. Two rough surfaces were generated in this way and the relative position of the two surfaces was chosen in order for the contact area to occupy now 30% of the total area.

The predictions of the transmissivity tensor components using both the power-series expansion and the Padé approximant as functions of $\xi Kn_i$, $i = x, y$ are reported in figure 7 along with the corresponding relative errors. Contrary to the case presented in the main text, for the present geometry and the ranges of $\xi Kn_i$, the predictions from the power-series expansion improve as more terms are included in the series. Moreover, the values of the relative errors are, overall, smaller than in the anisotropic case. Furthermore, the cross-over value of $\xi Kn_i$ beyond which the third-order correction term in the power-series expansion no longer provides more accuracy is almost twice that in the anisotropic case for the xx and yy components. This is a clear indication that the performance of the power-series expansion depends on the aperture field structure. In addition, the predictions from the Padé approximant are in excellent agreement (the relative error remains less than approximately $10^{-3}$) with the results for K, as in the anisotropic case reported in the main text. The differences between this case and the latter are obviously attributed to the geometries considered in the unit cells.

REFERENCES

BAKER, G.A. & GRAVES-MORRIS, P. 1996 *Padé Approximants*, 2nd edn, Encyclopedia of Mathematics and its Applications, Series Number 59. Cambridge University Press.

BERGSTRÖM, D. 2012 Rough surface generation & analysis. Available at: http://www.mysimlabs.com/surface_generation.html.

BOLAÑOS, S.J. & VERNESCU, B. 2017 Derivation of the Navier slip and slip length for viscous flows over a rough boundary. *Phys. Fluids* **29**, 057103.

GRAY, W.G. 1975 A derivation of the equations for multi-phase transport. *Chem. Engng Sci.* **30** (2), 229-233.

LASSEUX, D. & VALDÉS-PARADA, F.J. 2017 Symmetry properties of macroscopic transport coefficients in porous media. *Phys. Fluids* **29** (4), 043303.

LASSEUX, D., VALDÉS-PARADA, F.J., OCHOA-TAPIA, J.A. & GOYEAU, B. 2014 A macroscopic model for slightly compressible gas slip-flow in homogeneous porous media. *Phys. Fluids* **26** (5), 053102.

LASSEUX, D., VALDÉS-PARADA, F.J. & PORTER, M.L. 2016 An improved macroscale model for gas slip flow in porous media. *J. Fluid Mech.* **805**, 118–146.

LASSEUX, D., ZAOUTER, T. & VALDÉS-PARADA, F.J. 2023 Determination of Klinkenberg and higher-order correction tensors for slip flow in porous media. *Phys. Rev. Fluids* **8**, 053401.

LEE, H.B., YO, I.W. & LEE, K.K. 2013 The modified Reynolds equation for non-wetting fluid flow through a rough-walled rock fracture. *Adv. Water Resour.* **53**, 242–249.

MAXWELL, J.C. 1879 On stresses in rarified gases arising from inequalities of temperature. *Phil. Trans. R. Soc. Lond.* **170** (0), 231–256.

PATIR, N. & CHENG, H.S. 1978 An average flow model for determining effects of three dimensional roughness on partial hydrodynamic lubrication. *Trans. ASME* **100**, 12–17.

PÉREZ-RÀFOLS, F., LARSSON, R. & ALMQVIST, A. 2016 Modelling of leakage on metal-to-metal seals. *Tribol. Intl* **94**, 421–427.

PRAT, M., PLOURABOUÉ, F. & LETALLEUR, N. 2002 Averaged Reynolds equation for flows between rough surfaces in sliding motion. *Transp. Porous Med.* **48** (3), 291–313.

RUI, M.X., WEI, W., HEHUA, Z., YAJI, J., RUI, J. & XINBIN, T. 2023 Gas sealing behavior of gasketed segmental joints in shield tunnels: an experimental and computational study. *Tunnel. Undergr. Space Technol.* **132**, 104906.

TZENG, S.T. & SAIBEL, E. 1967 Surface roughness effect on slider bearing lubrication. *ASLE Trans.* **10**, 334–348.

WANG, J., TANG, D. & JING, Y. 2019 Analytical solution of gas flow in rough-walled microfracture at in situ conditions. *Water Resour. Res.* **55**, 6001–6017.

WHITAKER, S. 1999 *The Method of Volume Averaging*. Springer.
Effective transmissivity for slip flow in a fracture

main text for the Gaussian anisotropic fracture, albeit the correlation lengths are now $\sigma_{cx} = \sigma_{cy} = 0.05\ell$ maintaining the same root mean square roughness value. Two rough surfaces were generated in this way and the relative position of the two surfaces was chosen in order for the contact area to occupy now 30% of the total area.

The predictions of the transmissivity tensor components using both the power-series expansion and the Padé approximant as functions of $\xi Kn_i$, $i = x, y$ are reported in figure 7 along with the corresponding relative errors. Contrary to the case presented in the main text, for the present geometry and the ranges of $\xi Kn_i$, the predictions from the power-series expansion improve as more terms are included in the series. Moreover, the values of the relative errors are, overall, smaller than in the anisotropic case. Furthermore, the cross-over value of $\xi Kn_i$ beyond which the third-order correction term in the power-series expansion no longer provides more accuracy is almost twice that in the anisotropic case for the xx and yy components. This is a clear indication that the performance of the power-series expansion depends on the aperture field structure. In addition, the predictions from the Padé approximant are in excellent agreement (the relative error remains less than approximately $10^{-3}$) with the results for K, as in the anisotropic case reported in the main text. The differences between this case and the latter are obviously attributed to the geometries considered in the unit cells.

REFERENCES

BAKER, G.A. & GRAVES-MORRIS, P. 1996 *Padé Approximants*, 2nd edn, Encyclopedia of Mathematics and its Applications, Series Number 59. Cambridge University Press.

BERGSTRÖM, D. 2012 Rough surface generation & analysis. Available at: http://www.mysimlabs.com/surface_generation.html.

BOLAÑOS, S.J. & VERNESCU, B. 2017 Derivation of the Navier slip and slip length for viscous flows over a rough boundary. *Phys. Fluids* **29**, 057103.

GRAY, W.G. 1975 A derivation of the equations for multi-phase transport. *Chem. Engng Sci.* **30** (2), 229-233.

LASSEUX, D. & VALDÉS-PARADA, F.J. 2017 Symmetry properties of macroscopic transport coefficients in porous media. *Phys. Fluids* **29** (4), 043303.

LASSEUX, D., VALDÉS-PARADA, F.J., OCHOA-TAPIA, J.A. & GOYEAU, B. 2014 A macroscopic model for slightly compressible gas slip-flow in homogeneous porous media. *Phys. Fluids* **26** (5), 053102.

LASSEUX, D., VALDÉS-PARADA, F.J. & PORTER, M.L. 2016 An improved macroscale model for gas slip flow in porous media. *J. Fluid Mech.* **805**, 118–146.

LASSEUX, D., ZAOUTER, T. & VALDÉS-PARADA, F.J. 2023 Determination of Klinkenberg and higher-order correction tensors for slip flow in porous media. *Phys. Rev. Fluids* **8**, 053401.

LEE, H.B., YO, I.W. & LEE, K.K. 2013 The modified Reynolds equation for non-wetting fluid flow through a rough-walled rock fracture. *Adv. Water Resour.* **53**, 242–249.

MAXWELL, J.C. 1879 On stresses in rarified gases arising from inequalities of temperature. *Phil. Trans. R. Soc. Lond.* **170** (0), 231–256.

PATIR, N. & CHENG, H.S. 1978 An average flow model for determining effects of three dimensional roughness on partial hydrodynamic lubrication. *Trans. ASME* **100**, 12–17.

PÉREZ-RÀFOLS, F., LARSSON, R. & ALMQVIST, A. 2016 Modelling of leakage on metal-to-metal seals. *Tribol. Intl* **94**, 421–427.

PRAT, M., PLOURABOUÉ, F. & LETALLEUR, N. 2002 Averaged Reynolds equation for flows between rough surfaces in sliding motion. *Transp. Porous Med.* **48** (3), 291–313.

RUI, M.X., WEI, W., HEHUA, Z., YAJI, J., RUI, J. & XINBIN, T. 2023 Gas sealing behavior of gasketed segmental joints in shield tunnels: an experimental and computational study. *Tunnel. Undergr. Space Technol.* **132**, 104906.

TZENG, S.T. & SAIBEL, E. 1967 Surface roughness effect on slider bearing lubrication. *ASLE Trans.* **10**, 334–348.

WANG, J., TANG, D. & JING, Y. 2019 Analytical solution of gas flow in rough-walled microfracture at in situ conditions. *Water Resour. Res.* **55**, 6001–6017.

WHITAKER, S. 1999 *The Method of Volume Averaging*. Springer.
ZAOUTER, T., LASSEUX, D. & PRAT, M. 2018 Gas slip flow in a fracture: local Reynolds equation and upscaled macroscopic model. *J. Fluid Mech.* **837**, 413–442.

ZHOU, C.-B., CHEN, Y.-F., HU, R. & YANG, Z. 2023 Groundwater flow through fractured rocks and seepage control in geotechnical engineering: theories and practices. *J. Rock Mech. Geotech. Engng* **15**, 1–36.

ZINN-JUSTIN, J. 1971 Strong interactions dynamics with Padé approximants. *Phys. Rep.* **1**, 55–102.
ZAOUTER, T., LASSEUX, D. & PRAT, M. 2018 Gas slip flow in a fracture: local Reynolds equation and upscaled macroscopic model. *J. Fluid Mech.* **837**, 413–442.

ZHOU, C.-B., CHEN, Y.-F., HU, R. & YANG, Z. 2023 Groundwater flow through fractured rocks and seepage control in geotechnical engineering: theories and practices. *J. Rock Mech. Geotech. Engng* **15**, 1–36.

ZINN-JUSTIN, J. 1971 Strong interactions dynamics with Padé approximants. *Phys. Rep.* **1**, 55–102.