# A Stochastic Study of Flow Anisotropy and Channelling in Open Rough Fractures

Sophie Marchand, Olivier Mersch, Michael Selzer, Fabian Nitschke, Martin

Schoenball, Jean Schmittbuhl, Britta Nestler, Thomas Kohl

► To cite this version:

Sophie Marchand, Olivier Mersch, Michael Selzer, Fabian Nitschke, Martin Schoenball, et al. A Stochastic Study of Flow Anisotropy and Channelling in Open Rough Fractures. Rock Mechanics and Rock Engineering, 2020, 53, pp.233-249. 10.1007/s00603-019-01907-4. hal-02323211

HAL Id: hal-02323211

https://hal.science/hal-02323211v1

Submitted on 29 Sep 2025

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
**Title**

A Stochastic Study of Flow Anisotropy and Channelling in Open Rough Fractures

**Permalink**

https://escholarship.org/uc/item/9th1g9t7

**Journal**

Rock Mechanics and Rock Engineering, 53(1)

**ISSN**

0723-2632

**Authors**

Marchand, Sophie
Mersch, Olivier
Selzer, Michael
et al.

**Publication Date**

2020

**DOI**

10.1007/s00603-019-01907-4

Peer reviewed
# A Stochastic Study of Flow Anisotropy and Channelling in Open Rough Fractures

Sophie Marchand¹ · Olivier Mersch² · Michael Selzer² · Fabian Nitschke¹ · Martin Schoenball³ · Jean Schmittbuhl⁴ · Britta Nestler² · Thomas Kohl¹

¹ Karlsruhe Institute of Technology Institute of Applied Geosciences, Kaiserstr. 12, 76131 Karlsruhe, Germany ² Karlsruhe Institute of Technology Institute for Applied Materials, Straße am Forum 7, 76131 Karlsruhe, Germany ³ Lawrence Berkeley National Laboratory, 1 Cyclotron Road, MS74R316C, Berkeley, CA 94720, USA ⁴ School and Observatory of Earth Sciences, 5 Rue René Descartes, 67084 Strasbourg, France

Sophie Marchand soph.march@gmail.com

## Abstract

The quantification of fluid flow in rough fractures is of high interest for reservoir engineering, especially for deep geothermal applications. Herein, rough self-affine fractures are stochastically generated with incremental shear displacement and geometrically described by two aperture definitions, the vertical aperture $a_{\text{vert}}$ and the effective aperture $a_{\text{eff}}$. In order to compare their effect on fracture flow, such as anisotropy and channelling, Local Cubic Law (LCL) model-based 2D fluid flow is simulated. The particularity of this approach is the combination of a stochastic generation of self-affine fractures with a statistical analysis (560 individual realizations) of the impact of the LCL's aperture constraint on fracture flow. The results show that aperture definition affects the quantitative interpretation of flow anisotropy and channeling as well as the aperture distribution of the fractures with shearing. Higher values of mean aperture for a given fracture are found using $a_{\text{vert}}$, whereas the aperture standard deviation is larger with $a_{\text{eff}}$. In addition, flow anisotropy is significantly sensitive to aperture definition for small shear displacements and shows a relative higher dispersion with $a_{\text{eff}}$. Thus, LCL prediction models based on $a_{\text{vert}}$ are expected to lead to higher dispersion of anisotropy results with a higher uncertainty (factor ~ 2).

Realizations based on $a_{\text{vert}}$ lead to an enhanced clustering of high flow rates for higher shearing displacements. This channeling development results in higher total flow rates for these simulations. These findings support the direct calibration of pre-existing LCL anisotropy simulations based on $a_{\text{vert}}$ towards more representative results using $a_{\text{eff}}$.

**Keywords:** Fracture · Roughness · Aperture · Stochastic · Anisotropy · Channelling
List of Symbols

$\alpha$ Slope of the linear trend lines for $I_2$ in terms of $I_1$ (unitless)

$\beta$ Intercept of the linear trend lines for $I_2$ in terms of $I_1$ (unitless)

$\Delta h$ Height difference used to define $p$ (in m)

$\Delta r$ Vector distance used to define $p$ (in m)

$\Delta x$ Discretisation element length along the x-axis (in mm)

$\Delta y$ Discretisation element length along the y-axis (in mm)

$\lambda$ Arbitrary scaling factor used to define $p$ (unitless)

$\nu$ Kinematic viscosity of the fluid (in $m^2 s^{-1}$)

$\mu$ Dynamic viscosity of the fluid (in kg m$^{-1}$ s$^{-1}$)

$\rho$ Fluid density (in kg m$^{-3}$)

$\sigma^a$ Standard deviation of the aperture within a fracture (in mm)

$\sigma_{eff}^a$ Standard deviation of the aperture within a fracture based on the effective aperture definition (in mm)

$\sigma_{vert}^a$ Standard deviation of the aperture within a fracture based on the vertical aperture definition (in mm)

$\sigma_{N_{offset}}^I$ Standard deviation of the indicator $I$ associated with the shearing $N_{offset}$ (unitless)
<table>
  <tr>
    <td>σ&#770;<sup>I</sup><sub>N offset</sub></td>
    <td>Relative standard deviation for the indicator I at a given N<sub>offset</sub> (unitless)</td>
    <td>l</td>
    <td>Scaling of the maximum surface amplitude of the fractures (in mm)</td>
  </tr>
  <tr>
    <td>σ&#770;<sup>I<sub>1</sub></sup><sub>N offset</sub></td>
    <td>Relative standard deviation for the indicator I<sub>1</sub> at a given N<sub>offset</sub> (unitless)</td>
    <td>l<sub>w</sub></td>
    <td>Channel width (in mm)</td>
  </tr>
  <tr>
    <td>σ&#770;<sup>I<sub>2</sub></sup><sub>N offset</sub></td>
    <td>Relative standard deviation for the indicator I<sub>2</sub> at a given N<sub>offset</sub> (unitless)</td>
    <td>m(x, y, N<sub>offset</sub>)</td>
    <td>Height of the median between the top and bottom surface passing by the point (x, y, N<sub>offset</sub>) (in mm)</td>
  </tr>
  <tr>
    <td>AF</td>
    <td>Anisotropy factor (unitless)</td>
    <td>ṁ<sub>x</sub></td>
    <td>Finite derivative of m(x, y, N<sub>offset</sub>) along the x-axis (unitless)</td>
  </tr>
  <tr>
    <td>AF<sub>eff</sub></td>
    <td>Anisotropy factor based on the effective aperture definition (unitless)</td>
    <td>ṁ<sub>y</sub></td>
    <td>Finite derivative of m(x, y, N<sub>offset</sub>) along the y-axis (unitless)</td>
  </tr>
  <tr>
    <td>AF<sub>vert</sub></td>
    <td>Anisotropy factor based on the vertical aperture definition (unitless)</td>
    <td>N</td>
    <td>Total number of elements constituting a fracture (unitless)</td>
  </tr>
  <tr>
    <td>a</td>
    <td>Arithmetic mean aperture (in m)</td>
    <td>N<sub>offset</sub></td>
    <td>Value of the shear displacement along the x-axis (in mm or in number of element)</td>
  </tr>
  <tr>
    <td>ā</td>
    <td>Mean aperture within a fracture (in mm)</td>
    <td>n</td>
    <td>Number of elements constituting the side length of a fracture (in number of element)</td>
  </tr>
  <tr>
    <td>ā<sub>eff</sub></td>
    <td>Mean aperture within a fracture based on the effective aperture definition (in mm)</td>
    <td>n(m)</td>
    <td>Normal vector to the middle plane (unitless)</td>
  </tr>
  <tr>
    <td>ā<sub>vert</sub></td>
    <td>Mean aperture within a fracture based on the vertical aperture definition (in mm)</td>
    <td>OF</td>
    <td>Projection of the outgoing flow perpendicular to the pressure gradient (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>a<sub>e</sub></td>
    <td>Equivalent fracture aperture in the parallel-plate model (in mm)</td>
    <td>OF<sub>||</sub></td>
    <td>Projection of the outgoing flow perpendicular to the pressure gradient parallel to the shearing direction (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>a<sub>eff</sub></td>
    <td>Effective aperture (in mm)</td>
    <td>OF⊥</td>
    <td>Projection of the outgoing flow perpendicular to the pressure gradient perpendicular to the shearing direction (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>a<sub>max</sub></td>
    <td>Maximum local aperture observed among the fractures (in mm)</td>
    <td>p<sub>d</sub></td>
    <td>Probability density function of the roughness distribution (in m<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>a<sub>i</sub></td>
    <td>Aperture value of the ith element of a considered fracture (in mm)</td>
    <td>p</td>
    <td>Hydrodynamic pressure (in kg m<sup>-1</sup> s<sup>-2</sup>)</td>
  </tr>
  <tr>
    <td>a<sub>vert</sub></td>
    <td>Vertical aperture (in mm)</td>
    <td>Q</td>
    <td>Volumetric flow rate (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>a<sub>x,y</sub></td>
    <td>Local aperture at (x, y) (in m)</td>
    <td>Q1</td>
    <td>Lower quartile value (units of the variable considered)</td>
  </tr>
  <tr>
    <td>C</td>
    <td>Set of elements belonging to one of the n/2 + 1 channels (unitless)</td>
    <td>Q3</td>
    <td>Upper quartile value (units of the variable considered)</td>
  </tr>
  <tr>
    <td>C<sup>I</sup><sub>N offset</sub></td>
    <td>Centroid value for all indicators I of one offset (unitless)</td>
    <td>Q<sub>i</sub></td>
    <td>Volumetric flow rate vector at the ith element of the fracture (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>f</td>
    <td>Body force vector acting on the fluid (in kg m<sup>-2</sup> s<sup>-2</sup>)</td>
    <td>Q<sub>x</sub></td>
    <td>Projected volumetric flow rate on the x-axis (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>g</td>
    <td>Acceleration of gravity vector (in kg m<sup>-2</sup> s<sup>-2</sup>)</td>
    <td>Q<sub>y</sub></td>
    <td>Total volumetric flow rate for each element of the fracture identified by the coordinates (x, y) in the direction j where j = {X, Y} (in m<sup>3</sup> s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>grad(m)</td>
    <td>Gradient of the median m(x, y, N<sub>offset</sub>) (unitless)</td>
    <td>Re</td>
    <td>Reynolds number (unitless)</td>
  </tr>
  <tr>
    <td>H</td>
    <td>Hurst roughness exponent (unitless)</td>
    <td>Re<sup>local</sup><sub>max</sub></td>
    <td>Maximum local Reynolds number observed (unitless) element surface of a considered fracture (in m<sup>2</sup>)</td>
  </tr>
  <tr>
    <td>h<sub>bot</sub></td>
    <td>Height of a fracture element on its bottom surface (in mm)</td>
    <td>U</td>
    <td>Flow velocity (in m s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>h<sup>0</sup><sub>top</sub></td>
    <td>Height of a fracture element on its initial top surface before shearing (in mm)</td>
    <td>u</td>
    <td>Flow velocity vector (in m s<sup>-1</sup>)</td>
  </tr>
  <tr>
    <td>h<sup>N offset</sup><sub>top</sub></td>
    <td>Height of a fracture element on its top surface after shearing of N<sub>offset</sub> (in mm)</td>
    <td>x</td>
    <td>Fracture coordinate along the x-axis (in mm)</td>
  </tr>
  <tr>
    <td>I<sub>1</sub></td>
    <td>Relative proportion of channel area in a fracture (unitless)</td>
    <td>y</td>
    <td>Fracture coordinate along the y-axis (in mm.)</td>
  </tr>
  <tr>
    <td>I<sub>2</sub></td>
    <td>Maximum channel length in a given fracture normalized by n (unitless)</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>L</td>
    <td>Scaling for the total length of the fractures (in m)</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>L<sub>X</sub></td>
    <td>Fracture width along the x-axis (in m)</td>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td>L<sub>Y</sub></td>
    <td>Fracture width along the y-axis (in m)</td>
    <td></td>
    <td></td>
  </tr>
</table>
# 1 Introduction

Permeability prediction in natural media like fractured rocks is a stepping stone for the development of projects such as disposal sites (Bear et al. 1993) and underground tunnelling (Evans et al. 2013) where low permeability is aimed. However, for reservoir exploitation (Schmittbuhl et al. 2008; Zimmerman and Bodvarsson 1996) such as geothermal sites, high permeability is pursued using hydraulic stimulation to shear the fractured rocks. Through the Darcy's law for laminar flow, permeability is defined from the proportionality factor linking the flow rate and the fluid viscosity to the pressure gradient. In fractured rocks, fluid flow is affected by the local aperture which is typically linked to the local roughness distribution. Moreover, flow rate depends on mechanical deformations, i.e. shearing or normal opening. It is well known that this can lead to a dependency of flow on orientation, yielding effects of anisotropy and channelling which have been intensively studied by multiple authors (Berkowitz 2002). Auradou et al. (2001) illustrated flow anisotropy with laboratory studies using dyed fluid in a self-affine rough fracture. They observed a dependency of anisotropy on the lateral displacement of the fracture surfaces. Based on natural fracture replica, Gentier et al. (1997) also verified this phenomenon on the permeability field through laboratory experiments. They concluded that hydraulic permeability of a fracture depends clearly on the shear direction at displacements below 0.5 mm with significant changes in flow direction.

Méheust and Schmittbuhl (2001) quantified anisotropy by taking into account the direction of the pressure gradient and the geometrical heterogeneities. Auradou et al. (2006) performed an experimental and numerical study demonstrating that shear displacement induces anisotropy with enhanced permeability perpendicular to the shear direction. They extended their model to include flow channels developing perpendicularly to the shear displacement. Channelling is referring to the phenomenon of flow concentration along preferential pathways. Silliman (1989) demonstrated the development of channelling structures by laboratory experiments and underlined the associated presence of flow anisotropy. Channelling was experimentally established when small parts of a fracture plane can concentrate 90% of the fluid flow in single fractures (Rasmuson and Neretnieks 1986). By numerical studies, the importance of channelling for transport phenomena in a strongly heterogeneous medium, such as a fracture, was investigated (Tsang and Tsang 1989) and identified (Moreno and Tsang 1994).

These studies have an important impact on fractured geothermal reservoir systems which tend to be situated in tectonically active areas with an individual history of shearing events. Shearing is noticed through seismic events with a relative movement of fracture surfaces. This can take place under natural or under man-made conditions (“natural vs. induced earthquake”). In a larger context of a fracture network, these studies highlight that anisotropy and channelling are crucial to understand the
permeability patterns in reservoirs. Clearly, its quantification depends on the flow law applied and on the geometrical characterization of aperture. Considering a laminar flow and an individual fracture, the widely used equation to evaluate the influence of aperture variability on fluid flow computation is the cubic law (Brown et al. 1995, 1998; Witherspoon et al. 1980). With this approach, the fracture surface is simplified to parallel plates and inertial as well as nonlinear effects are neglected. Under purely tensile conditions when the fracture surfaces are displaced normally to the fracture plane even very rough surfaces could imply a constant aperture (Méheust and Schmittbuhl 2003). However, under normal or strike-slip faulting conditions fracture surfaces will be displaced, and aperture distribution will change locally. Given the crucial importance of mechanical interaction in a fractured geothermal system, its hydraulic impact needs to be investigated as function of the changing local aperture with displacement. In the following, we account for the aperture by defining a local transmissivity derived from the Local Cubic Law (LCL) that is a function of a stochastically generated roughness on the fracture surface. Utilization of the LCL is attractive for a broad range of engineering applications, at least as a first approach, due to its simplicity and computational efficiency. However, one of the critical issues related to LCL is the definition of its geometrical constraint which is the local aperture used in this equation (Konzuk and Kueper 2004; Oron and Berkowitz 1998; Wang Lichun et al. 2015; Zimmerman and Bodvarsson 1996).

The goal of the present work is to establish a recalibration between the LCL simulation results based on the commonly used vertical aperture and the more realistic effective aperture which considers the flow directionality. The limitation of the LCL are not put into question as our aim is to present the first results documenting an alternative to balance computational time and accuracy based on this law. To do so, we need to quantify the accuracy of the generally applied LCL models by accounting on stochastically generated fracture surface geometry. By its small-scale stochastic nature, it concludes on the bandwidth of uncertainty at large-scale numerical model approaches. To provide this quantification, a stochastic approach is presented to generate the roughness distribution of fracture surfaces. This process results in many realizations using the same fractal parameter extracted from a field observation. On these non-identical fractures, different shearing displacements are applied to simulate the individual hydraulic impact assuming viscous Poiseuille flow. The statistical dispersion of the generated datasets is analysed to highlight the uncertainty related to the definition of a LCL aperture. The stochastic investigations are concentrating on anisotropy and channelling, representing two principal hydraulic phenomena in fractured media. In the context of complex fractured systems, this study aims therewith also at improving the reliability of the hydraulically coupled processes, i.e. in geothermal models.
First, we describe the numerical framework and the procedure to generate the stochastic database in Sect. 2. In Sect. 3, the aperture distributions and the statistical dispersion of the hydraulic results are presented to analyse the associated anisotropy variations. Then, in Sect. 4, indicators are defined to quantify channelling and its dependence on shearing under the chosen stochastic approach.

## 2 Methods

### 2.1 Generation of Rough Self-Affine Surfaces

The procedure to generate the stochastic database of a sheared single fracture is conducted by constructing 3D surfaces with a software developed by Schmittbuhl and described in (Méheust and Schmittbuhl 2001). Herein, it is assumed that fractal geometries represent a good approximation of the roughness distribution of real fracture surfaces (Bouchaud 1997; Schmittbuhl et al. 1993, 1995). This roughness distribution is defined from self-affine surfaces with isotropic correlation functions. Its associated probability density function $p_d$ can be stated as follows:

$$p_d(\Delta h, \Delta \mathbf{r}) = \lambda^H p_d(\lambda^H \Delta h, \lambda \Delta \mathbf{r}), \quad (1)$$

assuming an arbitrary scaling factor $\lambda$ and the Hurst roughness exponent $H$. It indicates that $p_d$ is invariant regarding the transformations on the height difference $\Delta h$ and the distance $\Delta r$ defined in Eq. (1) (Méheust and Schmittbuhl 2001; Talon et al. 2010).

A roughness exponent of $H=0.8$ is chosen in this study being widely observed for granite (Amitrano and Schmittbuhl 2002; Schmittbuhl et al. 1993, 1995). Moreover, every surface generated (Fig. 1) is distinct as it is, respectively, based on a distribution extracted from a white noise generator (Press et al. 1992). The generated surfaces have a resolution of 2048 by 2048 discretization elements. In addition, we measure that the topography of the generated surfaces is $4 \times 10^{-6}$ mm. This value is relatively higher but comparable with the experimental measurements in Schmittbuhl et al. (2008) where the value of $2 \times 10^{-7}$ mm is found. Moreover, we verify that the Power Spectral Density (PSD) of the surfaces has a slope of $-2.6$ which is coherent with $H=0.8$ as this slope for self-affine surfaces is of $-1-2H$ (Méheust and Schmittbuhl 2001). To implement the stochastic approach a total of 560 fractures were generated. This number of fractures has been chosen to be able to quantify the accuracy of anisotropy and channelling phenomena.
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Example of surface generated with a roughness exponent of *H* = 0.8 after scaling. Note that *x*, *y*-axis are in m and *z*-axis is in mm. The colours refer to the height in the direction of the *z*-axis</figcaption></figure>

## 2.2 Sheared Fracture Generation

In a geoscientific context, fractured systems represent discontinuities in the solid rock. Under increasing stress loading, their surfaces are likely to be displaced which is a typical behaviour in deep geothermal reservoirs or in active tectonic systems. To obtain numerically a sheared synthetic fracture from a generated surface, we adopt the methodology of the study realized by Auradou et al. (2006) using sheared fractures made initially of complementary walls. After having duplicated the initial surface, its duplicate is translated horizontally by *N*<sub>offset</sub> with a periodic boundary due to the small displacements (<1% of the fracture length). Then, a vertical translation results automatically when both surfaces have at least one contact point (Fig. 2). With this last step, we place the study in the limit of a rigid approximation, and thus no plastic deformation and no rigid rotation are considered. The constructed fractures are scaled on the measurements from Schmittbuhl et al. (1993). Also, the scaling of the surfaces considers a total length *L* of 0.64m0.64m and a maximum surface amplitude of *l*=30mm. Based on a self-affine geometry, the fracture-scale studied can be included in larger reservoir models through upscaling procedures. Indeed, this approximation is enabled by the self-affine geometry which is based on the self-similarity property which is featured by scale invariance.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Illustration of the fracture construction from a generated surface. For clarity, the construction steps are shown projected in 2D. Note that the top surface refers to the duplicate of the bottom surface</figcaption></figure>

In total, we consider 88 shear displacements along the x-axis, noted $N_{offset}$, and for each of them 70 fractures are, respectively, sheared, resulting in a total number of fractures to be generated of 560. The value of $N_{offset}$ controls the shearing displacement and is given in number of elements (i.e. each length has 2048 elements). The equivalent physical lengths of the different shear displacements are presented in Table 1. The shear displacement is linearly increasing from 1 to 8 elements of a surface with the aim to study the early development of anisotropy and channelling processes. Due to the relatively small shearing displacement compared to the overall fracture size, the rigid approximation is applied as stated before. Indeed, a shearing displacement in the millimetre range [0.31;2.50] is resulting for a fracture of 0.64m. In a mechanically brittle medium, this could be associated with micro-earthquakes of magnitude $M<1$ (Wells and Coppersmith 1994) typically observed during enhanced geothermal exploitations.

<table><caption>**Table 1** Shear displacement expressed in number of elements and physical length (mm)</caption><thead><tr><th>N<sub>offset</sub> (number of elements)</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th></tr></thead><tbody><tr><td>Shear replacement</td><td>0.31</td><td>0.63</td><td>0.94</td><td>1.25</td><td>1.56</td><td>1.88</td><td>2.19</td><td>2.50</td></tr></tbody></table>

## 2.3 Local Aperture Distributions

The common practice to compute LCL is to measure the aperture normal to the fracture plane, meaning that a constant direction is considered (Brown 1989). Since this practice leads to an overestimation of transmissivity (Berkowitz 2002), a variety of aperture definitions have been proposed for fluid flow simulations. For the definition of an effective aperture herein, we follow the approach from Ge (1997) suggesting to use a “true aperture” calculated from the normal to the local orientation of the centerline between
both fracture surfaces. Other authors proposed to define the aperture as the largest sphere fitting into the fracture at a given node (Mourzenko et al. 1995) or defined flow-oriented apertures (Lichun et al. 2015).

For each of the 560 generated fractures, two definitions of aperture distributions were used: (1) $a_{\text{vert}}$ the vertical aperture and (2) $a_{\text{eff}}$ the effective (“true”) aperture (Fig. 3). This extraction is performed through the framework integrated in Pace3D (Hötzer et al. 2018). First, to obtain $a_{\text{vert}}$, we consider a constant vertical direction with respect to the horizontal mean plane at each element of a generated fracture. Also, taking the vertical distance between the lower and upper surfaces of the fracture at a fixed position, we obtain the local vertical aperture distribution. Then, a local $a_{\text{vert}}$ can be formulated as follows:

$$a_{\text{ver}}(x, y, N_{\text{offset}}) = h_{\text{top}}^{N_{\text{offset}}}(x, y) - h_{\text{bot}}(x, y) = h_{\text{top}}^{0}(x - N_{\text{offset}}, y) - h_{\text{bot}}(x, y), \quad (2)$$

with $h_{\text{top}}^{N_{\text{offset}}}$ and $h_{\text{bot}}$, respectively, the height of a fracture element on its top and bottom surfaces as labelled on Fig. 3. Note that the shearing is applied on the top surface of the fracture along the x-axis. Thus, the superscript $N_{\text{offset}}$

of $h_{\text{top}}^{N_{\text{offset}}}$ indicates the displacement of the top surface relatively to the bottom one. Moreover, the vertical displacement of the top surface is indirectly defined through the common contact point of the sheared surfaces. In addition, $h_{\text{top}}^{0}$ refers to the initial top surface that has not been sheared.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Schematization of vertical $a_{\text{vert}}$ and effective aperture $a_{\text{eff}}$ definitions projected onto the plane (x, z) in a portion of an open fracture sheared by $N_{\text{offset}}$. The grey lines represent the fracture surfaces and the red line is the centerline between the fracture surfaces</figcaption></figure>
Second, the effective aperture is defined through the gradient of the median $m$, noted **grad(m)**, between the top and bottom surfaces (Selzer 2014). Considering the shear displacement $N_{offset}$, we establish for a given fracture the local height of the median curve at the point $(x,y)$ as follows:

$$m(x, y, N_{offset}) = \frac{h_{bot}(x, y) + h_{top}^{N_{offset}}(x, y)}{2}, \qquad (3)$$

Then, $a_{eff}$ is determined as the distance between the intersection points of

$h_{top}^{N_{offset}}$ and $h_{bot}$ surfaces with the line of slope **grad(m)** passing by the point $(x,y,m)$. Depending on **grad(m)**, $a_{eff}$ is perpendicular with the flow direction making it more realistic than the $a_{vert}$ in the context of LCL application. Using finite differences, we calculate the gradient as follows:

$$\textbf{grad}(\mathbf{m}) = \begin{bmatrix} \frac{m(x+1,y,N_{offset}) - m(x-1,y,N_{offset})}{2\Delta x} \\ \frac{m(x,y+1,N_{offset}) - m(x,y-1,N_{offset})}{2\Delta y} \\ 0 \end{bmatrix} = \begin{bmatrix} \dot{m}_x \\ \dot{m}_y \\ 0 \end{bmatrix}, \qquad (4)$$

Then, we apply a rotation of $\pi/2$ on it and we obtain **n(m)**, the normal vector to the middle plane:

$$\mathbf{n}(\mathbf{m}) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & -1 \\ 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \dot{m}_x \\ \dot{m}_y \\ 0 \end{bmatrix} = \begin{bmatrix} \dot{m}_x \\ 0 \\ \dot{m}_y \end{bmatrix}, \qquad (5)$$

Thus, $a_{eff}$ is defined as follows:

$$a_{eff}(x, y, N_{offset}) = BT, \qquad (6)$$

where $B$ is the intersection point between the line of slope defined at Eq. (5) and passing by the point $(x,y,m)$ with the bottom surface of the fracture. Similarly, $T$ is the intersection between this line and the top surface.

To document the difference between the aperture definitions, we compute the aperture ratios $a_{eff}/a_{vert}$ in each fracture and present the distribution associated with one fracture where $N_{offset}=5$ (Fig. 4). Aperture ratios superior to one confirm that $a_{eff}$ is not only a projection of $a_{vert}$ onto the direction aligned with **BT**. We observed for each fracture of every shear displacement similar histograms. Also, the previous conclusions can be extended to all study.
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Histogram of the probability density in terms of aperture ratio $a_{\text{eff}}/a_{\text{vert}}$ for a given fracture at $N_{\text{offset}} = 5$. Values superior to one indicate that $a_{\text{eff}}$ is not only a projection of $a_{\text{vert}}$ onto the perpendicular direction given by **BT** shown in Fig. 3</figcaption></figure>

Furthermore, we compute the PSD for each of the simulated fractures and we present the results for a given fracture with $N_{\text{offset}}=6$ in Fig. 5. Similar shapes of the PSD for $a_{\text{vert}}$ and $a_{\text{eff}}$ are observed regardless of the shear displacement and the fracture considered. For spatial frequencies lower than $10\text{m}^{-1}$ (i.e. high spatial variations), both aperture definitions lead to the same weighting of these frequencies. From $10\text{m}^{-1}$ to $100\text{m}^{-1}$, the PSD of $a_{\text{eff}}$ presents the same trend than the one based on $a_{\text{vert}}$ but with higher values showing that $a_{\text{eff}}$ weights more this spatial range. Finally, from $100\text{m}^{-1}$ to the largest spatial frequencies, we notice that $a_{\text{eff}}$ has globally a decreasing weighting of this range compared to $a_{\text{vert}}$ whose variations are heterogeneous with smaller values than $a_{\text{eff}}$ of at least one order of magnitude. Thus, we have shown that $a_{\text{eff}}$ differs from $a_{\text{vert}}$ notably in the weighting of low spatial variations approximately below 10mm which include the aperture values. Also, $a_{\text{vert}}$ distributions present higher averaged apertures as the weight on higher spatial variations is decreasing more rapidly in the range starting from $100\text{m}^{-1}$ compared to $a_{\text{eff}}$ weighting over this range. Note that to investigate the resolution impact on the PSD curves, we computed it for a resolution smaller by factor 2 and we obtain similar shapes and relative variation according to the aperture definition (Fig. 6). Thus, the relationship between both aperture definitions is transposable to other resolutions.
<figure><img src="image_6.png" /><figcaption>Fig. 5 Power spectral density of $a_{\text{vert}}$ and $a_{\text{eff}}$ for a fixed fracture whom $N_{\text{offset}} = 6$ with a resolution of $2048 \times 2048$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6 Power spectral density of $a_{\text{vert}}$ and $a_{\text{eff}}$ for a fixed fracture whom $N_{\text{offset}} = 6$ with a resolution of $1024 \times 1024$ which is twice smaller compared to the Fig. 5</figcaption></figure>

To conclude, both extractions of local aperture distributions are projections where we obtain 2D local aperture distributions from the 3D fractures. Note that the horizontal mean plane resolution is identical to the resolution of the fracture (2048 × 2048). This numerical resolution is chosen for this study as the finest resolution reachable with an acceptable computational time (less than 10 h of simulation by fracture). Moreover, to analyse the variability of the local aperture distribution for each of the generated fractures, we define $\bar{a}$ representing the mean aperture defined as follows:

$$ \bar{a} = \frac{1}{N} \sum_{i=1}^{N} a_i, \qquad (7) $$

where $a_i$ is the aperture value of the ith element taken among the N elements that constitute the fracture. Note that this definition is identically

used for each aperture definition giving $\bar{a}_{\text{vert}}$ and $\bar{a}_{\text{eff}}$. Similarly, the distribution of the local aperture is also characterized by its dispersion through $\sigma^a$ which is the standard deviation of the aperture computed for both apertures as follows:

$$ \sigma^a = \sqrt{\frac{\sum_{i=1}^{N} (a_i - \bar{a})^2}{N - 1}}, \qquad (8) $$

## 2.4 Fluid Flow Governing Equations

Fluid flow through a single fracture is governed by the Navier-Stokes (NS) equations. In the case of incompressible and steady Newtonian flow, the NS equations can be expressed in vector form as follows (Foias et al. 2001):

$$ (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \cdot \mathbf{u} + \mathbf{f}, \qquad (9) $$
with **u** the flow velocity vector, ρ the fluid density, ρ the hydrodynamic pressure, ν the kinematic viscosity of the fluid and **f** the body force acting on the fluid, typically **f**=**g** where **g** is the acceleration of gravity. Solving the NS equations requires to solve nonlinear partial differential equations. The source of nonlinearity is in the advective term (u·∇)u which shows the inertial forces acting on the fluid. The Stokes equation can be obtained by neglecting the inertial term in the NS equation if the Reynolds number Re<<1, where Re is defined as ρQa/μ with ρ=10³kg m⁻³ the density of water, Q the volumetric flow rate, a the arithmetic mean aperture and μ=1.306×10⁻³Pa s the viscosity of water. The simulations in this work show a maximum local flow rate of approximately 10⁻⁸m³s⁻¹. Given a maximum local aperture of amax=5mm and a maximum velocity of Q=U×S~10⁻⁸×amax×L/2048, we obtain a maximum local Reynolds number

Re<sub>max</sub><sup>local</sup> ~ 10<sup>-11</sup> . This allows for reducing the NS equations to the linear Stokes equation:

$$0 = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \cdot \mathbf{u} + \mathbf{f}, \qquad (10)$$

The validity of the Stokes equation in rough fracture profiles has been studied in Brown et al. (1995); Brush and Thomson (2003); Mourzenko et al. (1995). Nevertheless, the complexity of the Stokes equation to solve fluid flow in rock fracture requires another level of simplification through geometric and kinematic assumptions allowed in the lubrication approximation context. Then, neglecting the kinematic forces, the following standard simplification is to use the Reynolds equation (Eq. (11)) to compute the pressure distribution with a “cubic law” for the flux (Brown 1987; Mourzenko et al. 1995; Tsang and Tsang 1989; Zimmerman et al. 1991):

$$\nabla \cdot \left( \frac{a_e^3(x,y)}{12\mu} \nabla p \right) = 0, \qquad (11)$$

$$Q_x = -L_y \frac{a_e^3}{12\mu} \frac{\partial p}{\partial x}, \qquad (12)$$

where a_e is the equivalent fracture aperture in the parallel-plate model, μ the dynamic viscosity of the fluid, Q_x is the projected volumetric flow rate on the x-axis and L_y is the fracture width. Due to their roughness, rock fractures present aperture variations. Thus, to simply consider this roughness, an approximation is to consider the “Local Cubic Law” with the hypothesis of an isotropic media:

$$Q_{j}^{x,y} = -L_{j} \frac{a_{e}^{3} \partial p_{x,y}}{12\mu \partial j}, \qquad (13)$$
$Q_j^{x,y}$

with $Q_j^{x,y}$ the total volumetric flow rate for each element of the fracture identified by the coordinates $(x,y)$ in the direction $j$ where $j=\{X,Y\}$ and $a_{x,y}$ the local aperture at $(x,y)$. Then, Eq. (13) is identically used to determine

$Q_X^{x,y}$

$Q_Y^{x,y}$

and $Q_j^{x,y}$. Widely applied for fluid flow and solute transport (Zimmerman and Bodvarsson 1996; Zimmerman et al. 1991), this simplified model has been investigated theoretically and numerically (Brown 1987; Brown et al. 1995; Brush and Thomson 2003; Mourzenko et al. 1995; Zimmerman and Bodvarsson 1996; Zimmerman et al. 1991). Finally, the flow rate is computed with the software Abaqus using 4-node linear quadrilateral elements to solve the Eq. (13) where the characteristic length is set to $L_Y=L_X=L=0.64\text{m}$ (see Sect. 2.1) and the viscosity of water at 10 °C is taken as $\mu=1.306\times10^{-3}\text{Pa s}$. Note that the temperature is considered constant at

$Q_j^{x,y}$

10 °C during the simulation. Note that herein $Q_j^{x,y}$ is based on the global length scale $L$. A discrete volumetric definition of it can be directly obtained

$Q_j^{x,y}$

by dividing $Q_j^{x,y}$ by the fracture dimensional resolution 2048.

Using Eq. (13), we want to quantify the evolution of flow anisotropy with shearing by considering the velocity in the outlet layer according to the pressure gradient direction as illustrated in Fig. 7. With the outlet plane as reference, the perpendicular projection of the outgoing flow (OF) is a measure for the total outgoing flow obtained and the potential fluid recovery associated. For reservoir exploitation, this flow rate is a crucial parameter; mimicking it by OF enables us to analyse the anisotropy with a field application perspective. Also, we define the OF for both pressure cases as follows:

$$ \mathrm{OF}_{\parallel} = \sum_{\substack{x=N; y=1 \\ \rightarrow Q_X \cdot \vec{x} > 0}}^{x=N; y=N} Q_X^{x,y}, \qquad (14) $$

$$ \mathrm{OF}_{\perp} = \sum_{\substack{x=1; y=N \\ \rightarrow Q_Y \cdot \vec{y} > 0}}^{x=N; y=N} Q_Y^{x,y}, \qquad (15) $$
<figure><img src="image_8.png" /><figcaption>Fig. 7 a, b Illustration of the outgoing flow (OF) definitions and boundary conditions for both pressure gradient directions {‖, ⊥}. The colored lines represent the outlet layer considered for each case of pressure gradient direction. Note that the shear displacement is made</figcaption></figure>

along the x-axis and the pressure gradient is parallel or perpendicular to this shearing direction. c, d Schematization of the boundary conditions applied in each cases of pressure gradient directions {‖, ⊥} (color figure online)

Based on the quantitative characterization of anisotropy established by Auradou et al. (2005), we quantify the dependence of the OF with the orientation of the hydraulic pressure gradient through the anisotropy factor (AF) defined as follows:

$$AF = \frac{OF_{\perp}}{OF_{\parallel}}, \qquad (16)$$

## 2.5 Boundary Conditions

A fixed hydraulic pressure is imposed at the inlet and outlet of the fracture as shown in Fig. 7. According to the pressure gradient direction considered, the inlet is set to 1bar1bar and the outlet pressure is equal to 0 bar. Then, the two other walls form a closed loop by being virtually connected to obtain an artificial infinite domain. Finally, the rough fracture surfaces are impermeable to flow. Thus, the flow is 2D in the (x,y)-plane.

## 3 Results

### 3.1 Local Aperture Distribution with Shear Displacement

Both aperture definitions {$a_{\text{vert}}, a_{\text{eff}}$} are successively applied to determine the local aperture distribution as illustrated in Fig. 8. As expected, the local $a_{\text{eff}}$ distribution reveals lower values compared to the distribution of $a_{\text{vert}}$. This observation is confirmed by the mean aperture $\bar{a}$ distribution over the
fracture replications. In order to evaluate the variations of $\bar{a}$, we display the boxplots of its distribution at each shearing offset in Fig. 9. Note, that for all boxplots shown here, the box extends from the lower (Q1) to the upper (Q3) quartile value with a line at the median. The whiskers show the range of data defined as the following interval: $[Q1-1.5×(Q3-Q1);Q3+1.5×(Q3-Q1)]$. Finally, the circles at the extremes represent the outliers whom values are outside of the previous range. Using this visualization, we can extend the previous remark by observing that $\bar{a}_{\text{vert}}$ presents globally higher values than $\bar{a}_{\text{eff}}$ at each offset. Moreover, the spreading of the distribution is larger for $\bar{a}_{\text{vert}}$ compared to $\bar{a}_{\text{eff}}$ and this spreading increases for both cases with higher shear displacements.

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Example of local aperture distribution obtained for a generated fracture whom $N_{\text{offset}} = 5$ according to the aperture definition $\{a_{\text{vert}}, a_{\text{eff}}\}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Boxplots of the mean aperture $\bar{a}$ in terms of shear displacement for each simulation configurations based on $\{a_{\text{vert}}, a_{\text{eff}}\}$. Note that the box extends from the lower (Q1) to the upper (Q3) quartile</figcaption></figure>

values with a line at the median. The whiskers show the range of data defined at Sect. 3.1. The circles at the extremes represent the outliers

Thus, the two definitions of aperture lead to a significant difference of the fracture characterization by $\bar{a}$, but also by the standard deviation $\sigma^a$ (Fig. 10). In the case of $\bar{a}$, we observe that $\sigma^a$ distributions show higher values for
higher shear displacement, regardless of the aperture definition. It is shown
that $\sigma_{\text{eff}}^a$ distribution presents higher spreading with more extreme values
than $\sigma_{\text{vert}}^a$. Furthermore, the spreading of $\sigma_{\text{vert}}^a$ increases with shearing
values, whereas $\sigma_{\text{eff}}^a$ varies heterogeneously. Finally, the distribution of local
vertical apertures reveals higher dispersion of $\bar{a}$ with shearing than $\bar{a}_{\text{eff}}$ but
concerning their standard deviation, $\sigma_{\text{vert}}^a$ is less dispersed than $\sigma_{\text{eff}}^a$.

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Boxplots of the aperture standard deviation $\sigma^a$ in terms of shear displacement for each simulation configurations based on $\{a_{\text{vert}}, a_{\text{eff}}\}$. Note that the box extends from the lower (Q1) to the</figcaption></figure>

upper (Q3) quartile values with a line at the median. The whiskers show the range of data defined at Sect. 3.1. The circles at the extremes represent the outliers

## 3.2 Fluid Flow and Total Outgoing Flow

Using $a_{\text{vert}}$ and $a_{\text{eff}}$ to spatially characterize the fractures, we study their impact on the fluid flow computed with Eq. (13). For one generated fracture, we analyse four different fluid flow models as presented in Fig. 11 which correspond to the combination of possible simulation configurations parametrized by $\{a_{\text{vert}}, a_{\text{eff}}\}$ and $\{\parallel, \perp\}$. As expected, this visualization qualitatively depicts an alignment of the high flow values with the pressure gradient direction. Moreover, the flow rate values calculated with $a_{\text{vert}}$ are higher than those calculated with $a_{\text{eff}}$. This observation is consistent with the previous finding showing that $\bar{a}_{\text{eff}}$ values are lower than $\bar{a}_{\text{vert}}$ values (Fig. 9).
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Example of local flow rate distributions in a given fracture for each simulation configurations based on {$a_{\text{vert}}$, $a_{\text{eff}}$} and {$\|\cdot\|, \perp\}$. The flow rate values are given for each simulation element in $m^3 s^{-1}$ which</figcaption></figure>

is equivalent to $10^9 \text{mm}^3 \text{s}^{-1}$. Note that the white arrows indicate the direction of the pressure gradient

In order to quantify the evolution of fluid flow anisotropy, the variations of $\text{OF}_{\|}$ and $\text{OF}_{\perp}$ (Eqs. 14, 15) are plotted against the shear displacement (Fig. 12). The spreading of the distributions increases with larger whiskers of the boxplots and further outlier values while shearing increases. $\text{OF}_{\|}$ and $\text{OF}_{\perp}$ distributions trends are relatively similar for a given aperture definition, whereas differences are significant of $\text{OF}_{\|}$ and $\text{OF}_{\perp}$ values for $a_{\text{vert}}$ compared to $a_{\text{eff}}$ (factor 2 larger).
<figure><img src="image_13.png" /><figcaption>**Fig. 12** Boxplots of the outgoing flow OF in terms of shear displacement for each simulation configurations based on {$a_{\text{vert}}, a_{\text{eff}}$} and {$\|\|, \perp$}. Note that the box extends from the lower (Q1) to the upper</figcaption></figure>

(Q3) quartile values with a line at the median. The whiskers show the range of data defined at Sect. 3.1. The circles at the extremes represent the outliers

<figure><img src="image_14.png" /><figcaption>**Fig. 13** Boxplots of AF in terms of shear displacement for each simulation configurations based on {$a_{\text{vert}}, a_{\text{eff}}$}. Note that the box extends from the lower (Q1) to the upper (Q3) quartile values with a line at</figcaption></figure>

the median. The whiskers show the range of data defined at Sect. 3.1. The circles at the extremes represent the outliers

## 3.3 Flow Anisotropy Results

To highlight and compare the variations of $OF_{\perp}$ and $OF_{\|}$ for both aperture definitions, we use their ratio AF (Eq. (16)) which facilitates the quantification of flow anisotropy. Figure 13 shows the evolution of AF over shearing for both aperture definitions. AF is computed for each fracture. It is shown that AF globally increases with increasing shear displacement. Compared to cases based on $a_{\text{vert}}$, the distributions associated with $a_{\text{eff}}$ present a larger dispersion independent of the shearing offset. This can be observed through the length of the whiskers as well as the number of outliers. Moreover, $AF_{\text{vert}}$
distributions reach higher values than AF<sub>eff</sub>. For small shearing displacements (up to 0.94mm) a relative fast increase of anisotropic behaviour with shearing can be observed. Towards higher values of displacement, this increase is saturating. The dispersion of AF values for both cases differ significantly. The spread (interquartile range) for a<sub>eff</sub> is by factor two larger than for a<sub>vert</sub>.

## 4 Discussions

The previous simulations demonstrate that anisotropy increases relatively fast with increasing shearing for small displacements. For larger displacements, this growth decreases. However, AF variations depend on the aperture definition with higher values observed for a<sub>vert</sub> and more dispersion (σ<sup>a</sup><sub>eff</sub) related to a<sub>eff</sub>. We also note that the range of values obtained for σ<sup>a</sup><sub>eff</sub> is of the same order of magnitude than the experimental values obtained by Auradou et al. (2005) on a granite replicated fracture sheared until 1mm1mm. However, replication of such experiment would be needed to confirm that the fractures simulated are geometrically coherent with rock samples.

Concerning the potential fracture scale effects on the previous results, we can consider the experimental study of Kumar Singh et al. (2016). One main result of this previous article is that on real samples increasing the fracture scale leads to a decrease of the flow rate. Also, in the present work, the measured flow rates might decrease if we consider a larger fracture scale. However, the impact on the anisotropy factor is not direct and will be investigated in an additional study. Finally, channelling appears to be fracture scale independent as described in Watanabe et al. 2015.

In order to deepen the investigation of the shearing effects on the hydraulic properties of a fracture, the study is expanded on flow channelling, a phenomenon that can be directly linked to anisotropy (Auradou et al. 2005, 2006).

Channelling is described as a spatial concentration of flow along preferential pathways. Indeed, flow tends to focus on high permeable zones and this creates preferred paths with higher velocity (Knudby and Ramírez 2005; Koltermann and Gorelick 1996). Also, for each simulation, channels are identified as largest connected area of high flow rate values (Fig. 14). In this work, we defined the identification threshold to be the third quartile (Q3) of flow rate values. The channel identification method consists in applying this threshold on the fluid flow and detecting the connected zones. Herein, connected areas are defined as areas showing flow rates above the threshold (Q3) and being connected, respectively. The maximum number of individual parallel channels developing in a fracture of a side length of *n* elements, with *n* being an even number, is *n*<sup>2</sup>+1.
<figure><img src="image_15.png" /><figcaption>**Fig. 14** Identification of channels (b) and its associated flow rate distribution (a), example based on $a_{eff}$ and in the $\perp$ case. The flow rate values are in $m^3 s^{-1}$ which is equivalent to $10^9 mm^3 s^{-1}$. Note that each colour on b labels a different channel and the background colour is set to black corresponding to the flow rates blow the Q75 quartile</figcaption></figure>

For this channel identification, two indicators are applied to follow channelling evolution inspired from Le Goc et al. (2010). The first one ($I_1$) quantifies the relative proportion of channel area in a fracture and the second one ($I_2$) quantifies the continuity of the flow path by measuring the maximum channel length in a given fracture. $I_1$ is defined as:

$$I_1 = \frac{1}{N} \sum_{i=1}^{N} \mathbb{1}_{\{Q_i > Q_{75}, i \in C\}}, \qquad (17)$$

where $Q_i$ is the flow rate vector at the $i$th element of the fracture, $C$ is the set of elements belonging to one of the $n/2+1$ channels and $N$ the total number of elements such that $N=n \times n$. If the flow rate is homogeneous over the fracture area, $I_1$ is equal to zero as no element verifies $Q_i > Q_{75}$. In the case of a flow rate superior to Q3 distributed among $n/2+1$ wide and $n$ long channel(s), we obtain $I_1 \sim (1/N) \times (n^2/2) = 1/2$. More generally, for channel(s) of width $I_w$ and length $n$, we have $I_1 = (1/N) \times n \times I_w = I_w/n$ for $I_w < n$. Thus, $I_1$ tends toward 1 for increasing channelling which is defined as a combined increase of the channel(s) size as well as their number.

$I_2$ is defined as the major axis of an ellipse that encapsulates the largest channel normalized by $n$. $I_1$ and $I_2$ are expressed in percentage where $I_2$ can exceed 100% as the orientation of the largest channel is not necessarily

parallel to $\vec{x} \to \vec{y}$ direction. Figure 15 displays the values $I_2$ over $I_1$ for each simulation and for each initial configuration of {$a_{vert}, a_{eff}$} and {$\|, \perp$}. Moreover, realizations of the same shearing displacement are shown with identical colour. Thus, the effect of shearing on channelling behaviour, represented by $I_1$ and $I_2$, can be studied for all cases.
<figure><img src="image_16.png" /><figcaption>**Fig. 15** Scatter plot of $I_2$ in terms of $I_1$ for each simulation configurations based on {$a_{\text{vert}}, a_{\text{eff}}$} and {$\|, \perp$}. The colors refer to the offset displacement and points of same color therefore represent the 70</figcaption></figure>

simulation replications of each offset. The black crosses represent the centroid of each offset cloud and the red line is the linear trend line fitted on them (color figure online)

By plotting the centroid for every shearing set of realizations, the general variation of $I_1$ and $I_2$ with shearing can be observed and quantified through the slope $\alpha$ and the intercept $\beta$ of the linear trend lines (Table 2). Larger values of $\alpha$ indicate a faster increase of $I_2$ compared to $I_1$. It is shown that $I_2$ increases faster than $I_1$ for all simulation configurations. Thus, the increasing connectivity of channelling represented by $I_2$ is more fostered by shearing than the proportion of channelled high flow rates indicated by $I_1$. It can also be observed that increasing shearing leads to significantly higher values of $I_1$ and higher or equivalent values of $I_2$ with higher associated uncertainty, defined as the area of the data distribution for a given shearing displacement. Additionally, we notice as well that the data sets are significantly spread along the $I_2$ axis reflecting a higher uncertainty for $I_2$ compared to $I_1$. Higher dispersion is also observed for small shear displacements and in particular for $N_{\text{offset}}=\{0.31,0.62\}$.

<table><caption>**Table 2** Coefficients of the linear trend lines for the distribution of centroids (Fig. 15)</caption><thead><tr><th>y = α × x + β</th><th>a<sub>vert</sub> and ||</th><th>a<sub>vert</sub> and ⊥</th><th>a<sub>eff</sub> and ||</th><th>a<sub>eff</sub> and ⊥</th></tr></thead><tbody><tr><td>α</td><td>7.22</td><td>4.39</td><td>3.52</td><td>9.16</td></tr><tr><td>β</td><td>-78.70</td><td>-4.97</td><td>-11.59</td><td>-115.32</td></tr></tbody></table>
The analysis of the channelling dependence on the aperture definition exhibits that the evolution of $\alpha$ for the configurations $\{\parallel, \perp\}$ is inverse compared to $\{a_{\text{vert}}, a_{\text{eff}}\}$ With $a_{\text{vert}}$, we observe that $\alpha$ is smaller by a factor of 2 comparing $\parallel$ to $\perp$ simulation cases. For $a_{\text{eff}}$, however, $\alpha$ more than doubles from $\parallel$ to $\perp$. Ultimately, we analysed here an extreme case in terms of the directions of pressure gradient and shearing. Herein, the relative evolution of $I_1$ and $I_2$ defined by $\alpha$ are inverse according to the aperture considered (Table 2). This reinforces the major impact of aperture definition in channelling process analysis with shearing on dispersion and relative evolution of channelling indicators.

To quantify the dispersion of the channelling indicators, we define the relative standard deviation for both values at a given $N_{\text{offset}}$ as follows:

$$ \tilde{\sigma}_{N_{\text{offset}}}^{I} = \frac{\sigma_{N_{\text{offset}}}^{I}}{C_{N_{\text{offset}}}^{I}}, \qquad (18) $$

where $\sigma_{N_{\text{offset}}}^{I}$ is the standard deviation of the indicator $I$ associated with the shearing $N_{\text{offset}}$. $C_{N_{\text{offset}}}^{I}$ is the centroid value for all indicators $I$ of one offset.

We observe that $\tilde{\sigma}_{N_{\text{offset}}}^{I_1}$ and $\tilde{\sigma}_{N_{\text{offset}}}^{I_2}$ differ by more than one order of magnitude (Fig. 16). This quantifies the previous observation that the dispersion of $I_2$ is globally larger than the one of $I_1$ (Fig. 15). Nevertheless, no homogenous tendency of the dispersion with increasing shearing can be observed regardless of the aperture definition for $I_1$. Considering the aperture definitions, only $I_1$ shows variations related to it, whereas $I_2$ varies more significantly with the pressure gradient direction.

<figure><img src="image_17.png" /><figcaption>Fig. 16 Plots of the relative standard deviations in percentage of $I_1$ and $I_2$ in terms of shear displacement for each simulation configurations based on $\{a_{\text{vert}}, a_{\text{eff}}\}$ and $\{\parallel, \perp\}$. See Eq. (18) for the literal expression of the relative standard deviation</figcaption></figure>

After considering the variations of $I_1$ and $I_2$ with shear displacement, we study these variations with the outgoing flow OF values (Fig. 17). Maximum
values of OF are obtained for high $I_1$ values corresponding mainly to the largest three shearing displacements spreading widely along the $I_2$ axis. If we analyse these graphs regarding aperture definition, it is shown that two clusters of high OF fractures can be identified (encircled grey in Fig. 17) for $a_{\text{vert}}$. This clustering may indicate that for $a_{\text{vert}}$ two types of high OF fractures can be differentiated, showing similar values of $I_1$ and significant differences of $I_2$. However, this typology is not shown in case of $a_{\text{eff}}$. Also, further experimental investigations should be carried out to determine if channelling clustering is observed. These results then can substantially contribute to identify which aperture definition represents the reality best.

<figure><img src="image_18.png" /><figcaption>**Fig. 17** Scatter plot of $I_2$ in terms of $I_1$ for each simulation configurations based on {$a_{\text{vert}}$; $a_{\text{eff}}$} and {$\|\cdot\|, \perp\|\cdot\|\}$}. Note that the colors refer to the outgoing flow value OF as defined by the Eq. (14) and (15). Grey</figcaption></figure>

circles indicate possible clustering of the fractures according to their channelling behaviour (color figure online)

## 5 Conclusions

Although fluid flow in fractured rock is investigated since many years, the impact of the complex interplay between fracture geometry and hydraulic flow under field scale scenarios remains an important research subject. At larger scales and in the geothermal context, this impact of the geometrical characterization of the fracture may lead to a re-examination of the principles used to predict the fluid circulation in the reservoir and its dynamic
with the stress field. The present paper presents an important step towards bridging small-scale to large-scale applications. The current investigation becomes, therefore, most important when the uncertainty range of large reservoir models needs to be quantified. By specifying a stochastic approach, both LCL models, based on *a*<sub>eff</sub> and on *a*<sub>vert</sub>, yield a large dispersion range of anisotropy (AF) with a doubled dispersion of AF for *a*<sub>eff</sub> compared to *a*<sub>vert</sub> models. Moreover, for both aperture definitions, we observe that higher permeability is obtained perpendicular to the shearing direction of the fracture. Being coherent with literature, this statement underlines that both LCL models based on *a*<sub>eff</sub> and *a*<sub>vert</sub> capture the anisotropy phenomenon. In addition, *a*<sub>eff</sub> is considered to be more representative to describe hydraulic phenomena than *a*<sub>vert</sub> as this first is locally perpendicular to the flow directionality. Also, former prediction LCL models with *a*<sub>vert</sub> must take into account that results for anisotropy may be overestimated and present higher dispersion range in the case where they will be based on the more representative *a*<sub>eff</sub> aperture. Finally, regarding the evolution trend of AF, we evaluated that small shear displacements up to 0.94mm0.94mm are more reliable to assure an increase of AF regardless of the aperture definition. These observations aim to support the direct calibration of pre-existing LCL anisotropy simulations based on *a*<sub>vert</sub> toward more representative results using *a*<sub>eff</sub>.

For the channelling indicators ($l_1, l_2$), we observe a higher dispersion of $l_1$ in the case of $a_{\text{eff}}$ and no specific trend for $l_2$ between the two aperture definitions. The channelling common trends to both apertures are a growing proportion of channels ($l_1$) with increasing shearing and a significant enhancement of $l_1$ with early shearing compared to later ones. Furthermore, the channel continuity ($l_2$) is largely dispersed and do not present a trend with shearing. Nevertheless, the variations of $l_2$ with outgoing flow (OF) values indicate a potential typology of the fracture channelling behaviour. Indeed, for a same range of $l_1$, two clusters are observed in the case of $a_{\text{vert}}$ along $l_2$. Similarly, to the observations made on anisotropy phenomena, the ones related to channelling can adjust existing LCL simulations using $a_{\text{vert}}$. In addition, the specific clustering of high OF observed with $a_{\text{vert}}$ is another line of research to deepen the understanding the deviations obtained from LCL models for a crucial parameter for reservoir exploitation.

In future studies, we will weigh our observations with experimental results based on similar synthetic fracture replications using the 3D-printing technology. Moreover, this experimental set up associated with additional numerical simulation will help us to increase the complexity of our model by considering the mechanical effects of shearing on the walls of the fracture and by considering the possible scale effects. This present work can also be a starting point to establish a more accurate prediction of the processes participating into the extension of the flow circulation in reservoir through shearing. In the geothermal context, it could enable us to forecast the dynamic of the heat exchange area with the stress field. Indeed, the
shearing displacements studied can be associated with earthquakes of magnitude 1–2 which are representative of a geothermal seismic activity. This paper may also be concluded by remaking that aperture is a determinant factor for fluid flow equations and thus, we estimate that the sensibility of fluid flow interpretation with aperture definition presented here can be found in more refined flow simulations. Finally, the goal to raise attention on the spatial fracture descriptions impacts on averaged and dispersion values of fluid flow phenomena is embodied by this study.

## Acknowledgements

The study is part of the Helmholtz portfolio project “Geoenergy”. The support from the program “Renewable Energies”, under the topic “Geothermal Energy Systems”, is gratefully acknowledged. We also thank the EnBW Energie Baden-Württemberg AG for supporting geothermal research at KIT.

## References

Amitrano D, Schmittbuhl J (2002) Fracture roughness and gouge distribution of a granite shear band. J Geophys Res 107:2375.
https://doi.org/10.1029/2002JB001761

Auradou H, Hulin JP, Roux S (2001) Experimental study of miscible displacement fronts in rough self-affine fractures. Phys Rev E Stat Nonlin Soft Matter Phys 63:066306. https://doi.org/10.1103/PhysRevE.63.066306

Auradou H, Drazer G, Hulin JP, Koplik J (2005) Permeability anisotropy induced by the shear displacement of rough fracture walls. Water Resour Res. https://doi.org/10.1029/2005WR003938Google Scholar

Auradou H, Drazer G, Boschan A, Hulin J-P, Koplik J (2006) Flow channeling in a single fracture induced by shear displacement. Geothermics 35:576–588.
https://doi.org/10.1016/j.geothermics.2006.11.004

Bear J, Tsang CF, De Marsily G (1993) Flow and contaminant transport in fractured rock. Academic Press, London

Berkowitz B (2002) Characterizing flow and transport in fractured geological media: a review. Adv Water Resour 25:861–884.
https://doi.org/10.1016/S0309-1708(02)00042-8

Bouchaud E (1997) Scaling properties of cracks. J Phys Condens Matter 9:4319

Brown SR (1987) Fluid flow through rock joints: the effect of surface roughness. J Geophys Res 92:1337–1347.
https://doi.org/10.1029/JB092iB02p01337

Brown S (1989) Transport of fluid and electric current in a single fracture. J Geophys Res 94:9429–9438. https://doi.org/10.1029/JB094iB07p09429
Brown S, Stockman HW, Reeves SJ (1995) Applicability of the Reynolds Equation for modeling fluid flow between rough surfaces. Geophys Res Lett 22:2537-2540. https://doi.org/10.1029/95GL02666

Brown S, Arvind C, Robert H (1998) Experimental observation of fluid flow channels in a single fracture. J Geophys Res 103:5125-5132. https://doi.org/10.1029/97JB03542

Brush DJ, Thomson NR (2003) Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations. Water Resour Res. https://doi.org/10.1029/2002WR001346

Evans DD, Rasmussen TC, Nicholson TJ (2013) Flow and transport through unsaturated fractured rock, vol 42. American Geophysical Union, New York

Foias C, Manley O, Rosa R, Temam R (2001) Navier-Stokes equations and turbulence. Cambridge University Press, Cambridge

Ge S (1997) A governing equation for fluid flow in rough fractures. Water Resour Res 33:53-61. https://doi.org/10.1029/96WR02588

Gentier S, Lamontagne E, Archambault G, Riss J (1997) Anisotropy of flow in a fracture undergoing shear and its relationship to the direction of shearing and injection pressure. Int J Rock Mech Min 34:94.e1-94.e12. https://doi.org/10.1016/S1365-1609(97)00085-3

Hötzer J, Reiter A, Hierl H, Steinmetz P, Selzer M, Nestler B (2018) The parallel multi-physics phase-field framework Pace3D. J Comput Sci Neth 26:1-12. https://doi.org/10.1016/j.jocs.2018.02.011

Knudby C, Ramírez J (2005) On the relationship between indicators of geostatistical, flow and transport connectivity. Adv Water Resour 28:405-421. https://doi.org/10.1016/j.advwatres.2004.09.001

Koltermann CE, Gorelick SM (1996) Heterogeneity in sedimentary deposits: a review of structure-imitating, process-imitating, and descriptive approaches. Water Resour Res 32:2617-2658. https://doi.org/10.1029/96WR00025

Konzuk JS, Kueper BH (2004) Evaluation of cubic law based models describing single-phase flow through a rough-walled fracture. Water Resour Res. https://doi.org/10.1029/2003WR002356

Kumar Singh K, Narain Singh D, Pathegama Gamage R (2016) Effect of sample size on the fluid flow through a single fractured granitoid. J. Rock Mech Geo Eng 8-3:329-340. https://doi.org/10.1016/j.jrmge.2015.12.004

Le Goc R, de Dreuzy J-R, Davy P (2010) Statistical characteristics of flow as indicators of channeling in heterogeneous porous and fractured media. Adv Water Resour 33:257-269. https://doi.org/10.1016/j.advwatres.2009.12.002

Lichun W, Bayani CM, Slottke DT, Ketcham RA, Sharp JM (2015) Modification of the Local Cubic Law of fracture flow for weak inertia, tortuosity, and
roughness. Water Resour Res 51:2064-2080.
https://doi.org/10.1002/2014WR015815

Méheust Y, Schmittbuhl J (2001) Geometrical heterogeneities and permeability anisotropy of rough fractures. J Geophys Res 106:2089-2102.
https://doi.org/10.1029/2000JB900306

Méheust Y, Schmittbuhl J (2003) Scale Effects Related To Flow In Rough Fractures. In: Kümpel H-J (ed) Thermo-hydro-mechanical coupling in fractured rock. Basel, Birkhäuser Basel, pp 1023-1050

Moreno L, Tsang C-F (1994) Flow channeling in strongly heterogeneous porous media: a numerical study. Water Resour Res 30:1421-1430.
https://doi.org/10.1029/93wr02978

Mourzenko V, Thovert J-F, Adler P (1995) Permeability of a single fracture; validity of the Reynolds equation. J Phys Paris II 5:465-482.
https://doi.org/10.1051/jp2:1995133

Oron AP, Berkowitz B (1998) Flow in rock fractures. Water Resour Res 34:2811-2825. https://doi.org/10.1029/98WR02285

Press WH, Teukolsky SA, Vetterling WT, Flannery BP (1992) Numerical recipes in C: the art of scientific computing. Cambridge University Press, Cambridge

Rasmuson A, Neretnieks I (1986) Radionuclide transport in fast channels in crystalline rock. Water Resour Res 22:1247-1256.
https://doi.org/10.1029/WR022i008p01247

Schmittbuhl J, Gentier S, Roux S (1993) Field measurements of the roughness of fault surfaces. Geophys Res Lett 20:639-641.
https://doi.org/10.1029/93GL00170

Schmittbuhl J, Schmitt F, Scholz C (1995) Scale invariance of crack surfaces. J Geophys Res 100:5953-5973. https://doi.org/10.1029/94JB02885

Schmittbuhl J, Steyer A, Jouniaux L, Toussaint R (2008) Fracture morphology and viscous transport. Int J Rock Mech Min 45:422-430.
https://doi.org/10.1016/j.ijrmms.2007.07.007

Selzer M (2014) Mechanische und Strömungsmechanische Topologieoptimierung mit der Phasenfeldmethode. Dissertation. Karlsruhe Institute of Technology

Silliman SE (1989) An interpretation of the difference between aperture estimates derived from hydraulic and tracer tests in a single fracture. Water Resour Res 25:2275-2283. https://doi.org/10.1029/WR025i010p02275

Talon L, Auradou H, Hansen A (2010) Permeability of self-affine aperture fields. Phys Rev E Stat Nonlinear Soft Matter Phys 82:046108.
https://doi.org/10.1103/PhysRevE.82.046108
Tsang YW, Tsang CF (1989) Flow channeling in a single fracture as a two-dimensional strongly heterogeneous permeable medium. Water Resour Res 25:2076-2080. https://doi.org/10.1029/WR025i009p02076

Watanabe N, Ishibashi T, Tsuchiya N (2015) Predicting the channeling flows through fractures at various scales. In: Proceedings world geothermal congress 2015

Wells DL, Coppersmith KJ (1994) New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. B Seismol Soc Am 84:974-1002

Witherspoon P, Wang J, Iwai K, Gale JE (1980) Validity of Cubic Law for fluid flow in a deformable rock fracture. Water Resour Res 16:1016-1024. https://doi.org/10.1029/WR016i006p01016

Zimmerman RW, Bodvarsson GS (1996) Hydraulic conductivity of rock fractures. Transp Porous Med 23:1-30. https://doi.org/10.1007/BF00145263

Zimmerman RW, Kumar S, Bodvarsson GS (1991) Lubrication theory analysis of the permeability of rough-walled fractures. Int J Rock Mech Min 28:325-331. https://doi.org/10.1016/0148-9062(91)90597-F