# Impact of the fracture contact area on macro-dispersion in single rough fractures

Anthony Beaudoin, Mohamad Farhat

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

## To cite this version:

Anthony Beaudoin, Mohamad Farhat. Impact of the fracture contact area on macro-dispersion in single rough fractures. Comptes Rendus. Mécanique, 2021, 349 (2), pp.203-224. 10.5802/crmeca.82 . hal-04077882

HAL Id: hal-04077882

https://univ-poitiers.hal.science/hal-04077882

Submitted on 18 Dec 2023

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Mécanique

Anthony Beaudoin and Mohamad Farhat

Impact of the fracture contact area on macro-dispersion in single rough
fractures

Volume 349, issue 2 (2021), p. 203-224

Published online: 27 April 2021

https://doi.org/10.5802/crmeca.82

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

©

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org

e-ISSN : 1873-7234
Short paper / *Note*

# Impact of the fracture contact area on macro-dispersion in single rough fractures

Anthony Beaudoin*,<sup>◎ a</sup> and Mohamad Farhat<sup>a</sup>

<sup>a</sup> Institut Pprime, SP2MI - Téléport 2, boulevard Marie et Pierre Curie, BP 30179, 86962 Futuroscope Chasseneuil Cedex, France

*E-mails*: anthony.beaudoin@univ-poitiers.fr (A. Beaudoin), mohamad.farhat@univ-poitiers.fr (M. Farhat)

**Abstract.** In the scientific literature, the study of the impact of the fracture contact area on macro-dispersion in single rough fractures is still an open question. In this work, we study numerically the combined effects of the fracture roughness and the fracture contact area on the non-Fickian transport in single rough fractures. In particular, we quantify the contribution of the fracture contact area on macro-dispersion. These objectives are achieved by estimating the macro-dispersion coefficient from Monte Carlo parallel numerical simulations in pure advection and advection-diffusion cases. When the fractional void $S_O$ is equal to 1 (i.e., for $\sigma_{lnb} < 0.25$), the Monte Carlo simulations show that macro-dispersion results of two contributions, dispersion caused by the heterogeneity of fracture apertures that induces a channelization of flow paths and molecular diffusion, as shown by the analytical solution proposed by Gelhar in 1993. When the fraction void $S_O$ is different from 1 (i.e., for $\sigma_{lnb} > 0.25$), a third mechanism plays in macro-dispersion. The presence of contacts or obstacles causes a disruption of flow paths. This mechanism is identical to that induced by the fracture roughness with a lower amplitude. Its amplitude is the function of the fractional void $S_O$.

**Keywords.** Single rough fractures, Fracture contact area, Macro-dispersion, Non-Fickian transport, Monte Carlo simulations, RW model.

Manuscript received 20th December 2020, revised 9th February 2021 and 17th March 2021, accepted 26th March 2021.

## 1. Introduction

The nearly ubiquitous presence of fractures on the crustal rocks is the direct consequence of tectonic processes. The discrete fracture networks dominate the hydromechanical and hydrological behavior of subsurface rocks [1, 2]. To understand the effect of discrete fracture networks on this behavior is essential for a number of environmental applications such as petroleum

<footnote>* Corresponding author.</footnote>
recovery, geothermal energy, contaminated groundwater remediation, and geological storage of carbon [3, 4]. With the current technology, the three-dimensional (3D) measurement of geometrical characteristics of fractures in the fields is impossible [2, 5]. Thus the study of the behavior of fractures should be made with single fractures [3, 4, 6, 7].

In the scientific literature, it was shown that the particle transport has a non-Fickian behavior in a single rough fracture. The particle transport cannot then be adequately described by the advection-dispersion equation [8]. This was explained by the heterogeneity of fracture apertures that induces complex flow fields [3, 6, 9–12]. Numerous numerical studies have been made to investigate particle transport in single fractures with different physical complexities, solutes or dense particles transport in a uniform fracture or in a rough fracture with spatially variable apertures [6, 9, 13–23].

In the aforementioned numerical studies, the flow field was generated from the Poiseuille law or from the classical cubic law [6, 9, 13–23]. Solving the Navier-Stokes equations and the advection-dispersion equation in a 3D single fracture is very expensive in CPU time. That is why few numerical studies have solved the Navier-Stokes equations directly in a 3D single fracture [6, 24, 25]. Thus many researchers prefer to use the Lagrangian random walk particle tracking (RWPT) algorithm with two-dimensional (2D) flow fields generated from the classical local cubic law [9, 26–28].

Another difficulty is to model the geometrical properties of a single rough fracture. Due to causes like shear solicitations or thermal stresses or chemical processes, the two surfaces of a single rough fracture do not have identical shapes [2, 29–33]. This difficulty is also accentuated by the fact that the two surfaces of a single rough fracture can move in response to boundary stresses or displacements [2, 34, 35]. This results in a modification of the fracture roughness and the fracture contact area [36–39].

In the scientific literature, the study of the impact of the fracture contact area on macro-dispersion in single rough fractures is still an open question. The researchers have observed that the increase of the fracture contact area with normal stresses or normal deformations induces a channeling of the flow which in turn leads to a non-Fickian transport [12, 40–43]. In this paper, we study numerically the combined effects of the fracture roughness and the fracture contact area on the non-Fickian transport. In particular, we quantify the contribution of the fracture contact area on macro-dispersion.

In order to reach our objective, we use a methodology to determine the macro-dispersion coefficient from Monte Carlo parallel numerical simulations in pure advection and advection-diffusion cases [44–48]. This paper is organized according to the description of this methodology: Section 2—Fracture aperture generation; Section 3—Flow velocity calculation; Section 4—Particle tracking simulation; and Section 5—Macro-dispersion estimation. In these four sections, the numerical methods used are described and validated. The numerical results are presented and discussed in Sections 6 and 7, respectively.

## 2. Fracture aperture generation

The generation of single rough fractures is based on the work presented by Mourzenko et al. in 1996 [30, 31, 49]. The model, developed by Mourzenko et al., allows to generate a single rough fracture presenting a fracture contact area.

As shown in Figure 1 [30, 31, 49], the single rough fracture is delimited by two surfaces $z^{\pm}$ oscillating randomly around two average planes $h_0^{\pm}$, parallel to the xy-plane

$$z^{\pm}(\mathbf{x}) = h_0^{\pm} + h^{\pm}(\mathbf{x}), \qquad (1)$$
<figure><img src="image_6.png" /><figcaption>**Figure 1.** Description and symbols for the fracture geometry [30].</figcaption></figure>

where $h^{\pm}$ represents the fluctuations of the surface around average planes. From the previous equation, the geometrical properties of the single rough fracture studied in this paper are given by

$$b_m = \langle \omega(\mathbf{x}) \rangle \quad \text{with} \quad \omega(\mathbf{x}) = z^{+}(\mathbf{x}) - z^{-}(\mathbf{x}), \qquad (2)$$

$$b(\mathbf{x}) = \begin{cases} \omega(\mathbf{x}), & \omega(\mathbf{x}) > 0, \\ 0, & \omega(\mathbf{x}) \le 0, \end{cases} \qquad (3)$$

where $b_m$ is the distance between the two average planes, $\omega$ is the separation between the two random surfaces, and $b$ is the fracture aperture. $\langle \cdot \rangle$ denotes the statistical expectation.

In (1), the fluctuations $h^{\pm}$ of two surfaces of the single rough fracture are characterized by a Gaussian probability distribution $\phi$ and a spatial covariance function $C$, which are defined by [50–53]

$$\phi(h^{\pm}) = \frac{1}{\sqrt{2\pi\sigma_h^2}} \exp\left[-\frac{h^{\pm 2}}{2\sigma_h^2}\right] \quad \text{with} \quad \sigma_h = \sigma_h^+ = \sigma_h^-, \qquad (4)$$

$$C(\mathbf{r}) = \sigma_h^2 \exp[-(\|\mathbf{r}\|/\lambda)^{2H}] \quad \text{with} \quad 0 \le H \le 1, \qquad (5)$$

where $\sigma_h$ is the standard deviation of $h^{\pm}$, $\lambda$ is the correlation length, and $H$ is the Hurst exponent. The two surfaces of the single rough fracture are self-affine for $0 < H < 1$ and Gaussian for $H = 1$. The mean of $h^{\pm}$ is assumed to be zero. In this paper, we have considered Gaussian surfaces, that is, $H = 1$.

The fracture contact area can be quantified by the fractional void $S_O$, which is defined by [30, 31, 49]

$$S_O = \frac{A_O}{A}, \qquad (6)$$

where $A_O$ and $A$ are the projection of the total surface of the void volume of the fracture and the cross-sectional area of the fracture plane, respectively. $S_O$ is also equal to the portion of fracture apertures where the fracture aperture $b$ is strictly positive

$$S_O = \frac{1}{A} \int_A B(\mathbf{x}) \, dA \quad \text{with} \quad B(\mathbf{x}) = \begin{cases} 1, & \omega(\mathbf{x}) > 0, \\ 0, & \omega(\mathbf{x}) \le 0, \end{cases} \qquad (7)$$
<figure><img src="image_7.png" /><figcaption>**Figure 2.** Schematic diagram of the coordinate system and boundary conditions for the flow and solute transport simulations through a single rough fracture.</figcaption></figure>

where B is the phase function of void zones. In this paper, the two surfaces of the single rough fracture are uncorrelated. This implies that the statistical expectation $\langle S_O \rangle$ of the fractional void area is given by

$$ \langle S_O \rangle = \frac{1}{2} \operatorname{erfc} \left( -\frac{b_m}{2\sigma_h} \right). \qquad (8) $$

The assumption of two uncorrelated surfaces of the single rough fracture also allows to find an expression for the statistical expectation and standard deviation of the fracture aperture, $\langle b \rangle$ and $\sigma_b^2$ [30, 31, 49]

$$ \langle b \rangle = S_O b_m + R \quad \text{with} \quad R = \frac{\sigma_h}{\sqrt{\pi}} \exp \left( -\frac{b_m^2}{4\sigma_h^2} \right), \qquad (9) $$

$$ \sigma_b^2 = S_O(b_m^2 + 2\sigma_h^2) + Rb_m - \langle b \rangle^2. \qquad (10) $$

The fracture aperture field is generated within a computational domain of dimensions $L_x \times L_y$ in the x- and y-directions with a spatial step $\Delta$ identical in all the directions (see Figure 2). The Fourier method is applied here for evaluating the fracture aperture field [54, 55]. The Fourier method is the most common of the spectral methods. The representation of random fields in Fourier spaces makes these spectral methods more efficient than other methods such as sequential Gaussian methods [56]. The parallel library FFTW is used for performing the Fourier method [57, 58]. Figure 3 shows the fluctuations $h^\pm$ and fracture aperture b fields for $\sigma_h = 3$ mm, $\lambda = 4$ mm and $b_m = 12$ mm. The transverse dimension of the computational domain, $L_y = 1024$ mm, has been chosen to simulate correctly the flow. Its longitudinal dimension $L_x$ has been fixed following a convergence study. In Figure 4, the standard deviation $\sigma_h$ of fluctuations $h^\pm$ has been plotted as a function of the longitudinal dimension $L_x$ with $b_m = 4$ mm and $\lambda = 4$ mm. The exact value of $\sigma_h$ is 3, represented by the horizontal line in Figure 4. The convergence is reached when $L_x = 4096$ mm. The generation of the fracture aperture b has been verified by comparing the numerical values of $\sigma_b$ and $\langle b \rangle$ with the analytical solutions given by (9) and (10). In Figure 5, $\sigma_b$ and $\langle b \rangle$ have been plotted as functions of $b_m$ with $\sigma_h = 1$ mm and $\lambda = 4$ mm. We can observe that a computational domain defined with $L_x \times L_y = 4096$ mm $\times$ 1024 mm and $\Delta = 1$ mm allows to generate correctly the geometrical properties of the single rough fracture.
<figure><img src="image_8.png" /><figcaption>**Figure 3.** Example of fluctuations $h^±$ (top) and fracture aperture $b$ (bottom) fields in a computational domain of dimensions $L_x \times L_y = 4096$ mm × 1024 mm with $\Delta = 1$ mm for $\sigma_h = 3$ mm, $\lambda = 4$ mm, and $b_m = 12$ mm.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4.** Standard deviation of fluctuations $h^±$, $\sigma_h$ as a function of the computational domain size in the x-direction $L_x$ with $b_m = 4$ mm, $\lambda = 4$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm. (The exact value of $\sigma_h$ is 3, represented by the horizontal line.)</figcaption></figure>

**3. Flow velocity calculation**

When the Reynolds number $Re$ is lower than 1 (or the flow velocity is low) and the mean plane of the single rough fracture presents weak fluctuations, the Reynolds equation (or the local cubic law) can be used, in lieu of the full Navier-Stokes equations, to describe the flow in fractures [1, 25, 41, 59].
<figure><img src="image_10.png" /><figcaption>**Figure 5.** Standard deviation $\sigma_b$ and mean $\langle b \rangle$ of the aperture $b$ as functions of the distance between the two average planes $b_m$ with $\sigma_h = 1$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm. (The analytical solutions of $\sigma_b$ and $\langle b \rangle$ are given by (9) and (10), respectively.)</figcaption></figure>

Assuming that the fluid is incompressible, in a steady state, the governing equation can be written as

$$ \frac{\partial}{\partial x} \left( T_{xx} \frac{\partial H}{\partial x} \right) + \frac{\partial}{\partial y} \left( T_{yy} \frac{\partial H}{\partial y} \right) = 0, \quad (11) $$

where $H$ is the hydraulic head, $T_{xx}$ and $T_{yy}$ are the fracture transmissivities in the x- and y-directions, respectively. In this paper, the transmissivity is assumed to be isotropic and is defined by the Cubic law

$$ T_{xx} = T_{yy} = T = \frac{\rho g b^3}{12\mu}, \quad (12) $$

where $\rho$ is the fluid density, $\mu$ the dynamic viscosity of the fluid, and $g$ the gravity acceleration. To ensure the uniqueness of the solution, the boundary conditions are fixed hydraulic heads on the vertical sides of computational domain and periodic conditions on the horizontal sides of computational domain (see Figure 2). Thus the main flow direction is parallel to the x-direction because of these boundary conditions. The value of hydraulic heads is chosen in such a way as to keep the Reynolds number below 1.

The components of the mean flow velocity in both the x- and y-directions are calculated by [9, 17]

$$ v_x = - \frac{\rho g b^2}{12\mu} \frac{\partial H}{\partial x} \quad \text{and} \quad v_y = - \frac{\rho g b^2}{12\mu} \frac{\partial H}{\partial y}. \quad (13) $$

Equation (11) is discretized by means of a finite volume scheme with a harmonic composition rule for the transmissivity between adjacent mesh cells [60, 61]. The linear system, giving the hydraulic head $H$, is solved by using the algebraic multigrid iterative method implemented in HYPRE [62, 63]. The two components of the mean velocity are then computed on each grid face from (13).

Figure 6 shows the transmissivity $T$ and the horizontal flow velocity component $v_x$ in the computational domain for $\sigma_h = 3$ mm, $\lambda = 4$ mm, and $b_m = 12$ mm. For all the numerical simulations performed in this paper, the fluid density, the dynamic viscosity of the fluid, and the
<figure><img src="image_11.png" /><figcaption>**Figure 6.** Example of transmissivity $T$ (top) and horizontal flow velocity component $v_x$ (bottom) fields in a 2D computational domain of dimensions $L_x \times L_y = 4096$ mm $\times$ 1024 mm with $\Delta = 1$ mm for $\sigma_h = 3$ mm, $\lambda = 4$ mm, and $b_m = 12$ mm.</figcaption></figure>

gravity acceleration are fixed to $\rho = 1000$ kg/m$^3$, $\mu = 10^{-3}$ Pa·s, and $g = 9.81$ m/s$^2$, respectively, for a temperature of 20 °C. We have estimated the Reynolds number $Re = \rho b_m \langle |v| \rangle / \mu = 0.6$ by taking the distance between the two average planes $b_m = 12$ mm as the characteristic length and the statistical expectation of the mean flow velocity $\langle |v| \rangle = 0.05$ mm/s as the characteristic velocity. In Figure 7, the dimensionless parameter $T_R/T_e$ has been plotted as a function of the mean separation $b_m/\sigma_h$. This dimensionless parameter corresponds to the ratio between the Reynolds transmissivity $T_R = \mu \langle v_x b \rangle / \langle \nabla H \rangle$ and the effective transmissivity $T_e = \langle b \rangle^3 / 12$, as defined in the work by Mourzenko et al. in 1995 [64]. $\langle \nabla H \rangle$ corresponds to the statistical expectation of the hydraulic head gradient. In the following of this work, the fixed parameters are $L_x = 4096$ mm, $L_y = 1024$ mm, $\Delta = 1$ mm and $\lambda = 4$ mm. The numerical results are in good agreement with the numerical results obtained by Mourzenko et al. in 1995 [64].

## 4. Particle tracking simulation

In this paper, the RWPT method is used to simulate the solute transport in the single rough fracture. The RWPT method has proven its efficiency in numerous studies [9, 15, 27, 28, 65].

Considering only the mechanisms of advection and diffusion in the single rough fracture and neglecting the matrix diffusion in the rock formations [12, 66], the solute transport equation is given in the x- and y-directions by [67]

$$x^{n+1} = x^n + v_x(x^n, y^n)\delta t + \mathcal{N}(0, 1)\sqrt{2D_m\delta t}, \quad (14)$$

$$y^{n+1} = y^n + v_y(x^n, y^n)\delta t + \mathcal{N}(0, 1)\sqrt{2D_m\delta t}, \quad (15)$$

where $(x^n, y^n)$ are the particle location in the Cartesian space at time $n\delta t$, $n$ refers to the current time step, $n+1$ is the succeeding time step, $\mathcal{N}(0, 1)$ is a random selection from the standard normal distribution with zero mean and unit variance, and $D_m$ is the molecular diffusion coefficient.

Within grid cells, the mean flow velocity $v$ is derived from linear interpolations in both directions, as it is the sole interpolation scheme that verifies the continuity equation [68–72].
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Dimensionless transmissivity $T_R/T_e$ as a function of the mean separation $b_m/\sigma_h$ with $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm.</figcaption></figure>

$N_p$ particles are injected in an injection window of size $0.8 \ L_y$, orthogonal to the flow direction and located at least five correlation lengths $\lambda$ downstream from the side of the system (see Figure 2). The ratio between the mean advection rate and the molecular diffusion is called the Péclet number $Pe$ given by

$$Pe = \frac{b_m \langle |v| \rangle}{D_m} \qquad (16)$$

In the numerical simulations performed for this paper, we have studied two values of the diffusion coefficient, $D_m = 0$ and $0.01 \ mm^2/s$. It allows to study the purely advective regime ($Pe = \infty$) and the advective-diffusive regime ($Pe < 100$).

As in the scientific literature [44–48], the Courant-Friedrichs-Levy condition is used to maintain the stability and accuracy of the Euler scheme by adapting the time step $\delta t$

$$\delta t = \frac{1}{N_{\delta t}} \min \left( \frac{\Delta}{|v(x(t))|}, \frac{\Delta^2}{2D_m} \right), \qquad (17)$$

where $N_{\delta t}$ is a positive integer representing the number of time steps performed by the particle in a cell and $|v(x(t))|$ is the norm of the mean flow velocity at the location $x(t)$.

A parametric study has been performed to fix the value of $N_{\delta t}$ in order to prevent the particles from entering the fracture contact areas. This allows to verify the assumption that the matrix diffusion is neglected in the rock formation. In this parameter study, the single rough fracture is replaced by a single smooth fracture conserving the same values of the fractional void $S_O$ and distance between the two average planes $b_m$. In the following work, this single smooth fracture is called a single equivalent fracture. Figure 8 shows the transmissivity $T$ and horizontal flow velocity component $v_x$ fields of a single equivalent fracture obtained from a single rough fracture with $\sigma_h = 3$ mm and $b_m = 12$ mm. In this single equivalent fracture, the motion of a cloud of 10,000 particles has been simulated by considering the advection and molecular diffusion effects. The statistical expectation $\langle |v| \rangle$ of the mean flow velocity and the diffusion coefficient $D_m$ are fixed to 0.05 mm/s and 0.01 mm²/s, respectively. The injection process of particles is identical between the single rough and equivalent fractures. Figure 9 shows the cloud of particles around a fracture contact area selected in the single equivalent fracture. The value of $N_{\delta t}$ used has been
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Transmissivity *T* and horizontal flow velocity component *v*<sub>x</sub> fields in a single equivalent fracture obtained from a single rough fracture with *b*<sub>m</sub> = 12 mm, *σ*<sub>h</sub> = 3 mm, *λ* = 4 mm, *L*<sub>x</sub> = 4096 mm, *L*<sub>y</sub> = 1024 mm, and *Δ* = 1 mm.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 9.** Cloud of particles around a fracture contact area in the single equivalent fracture generated with *b*<sub>m</sub> = 12 mm, *σ*<sub>h</sub> = 3 mm, *λ* = 4 mm, *L*<sub>x</sub> = 4096 mm, *L*<sub>y</sub> = 1024 mm, and *Δ* = 1 mm.</figcaption></figure>

chosen from Figure 10 in which the number *N*<sub>pp</sub> of particles trapped in the same fracture contact area *N*<sub>pp</sub> is plotted as a function of *N*<sub>δt</sub>. There are no particles entering the fracture contact area from *N*<sub>δt</sub> = 20.
<figure><img src="image_15.png" /><figcaption>Figure 10. Number $N_{pp}$ of particles trapped in the fracture contact area of the single equivalent fracture as a function of the number of time steps $N_{\delta t}$ performed by the particles in a cell.</figcaption></figure>

## 5. Macro-dispersion estimation

The estimation of longitudinal and transverse dispersivities is performed in two steps [44–48]. First, the longitudinal and transverse dispersivities are evaluated on one Monte Carlo simulation by means of the following expressions

$$ \alpha_L^i(t) = \frac{1}{2\langle|v|\rangle} \frac{d}{dt} (\langle x^2(t) \rangle_i - \langle x(t) \rangle_i^2), \quad (18) $$

$$ \alpha_T^i(t) = \frac{1}{2\langle|v|\rangle} \frac{d}{dt} (\langle y^2(t) \rangle_i - \langle y(t) \rangle_i^2). \quad (19) $$

The $k$th moments of the solute cloud of the $i$th simulation, $\langle x^k(t) \rangle_i$ and $\langle y^k(t) \rangle_i$, are evaluated from

$$ \langle x^k(t) \rangle_i = \frac{1}{N_p} \sum_{j=1}^{j=N_p} (x_j^i(t))^k \quad \text{and} \quad \langle y^k(t) \rangle_i = \frac{1}{N_p} \sum_{j=1}^{j=N_p} (y_j^i(t))^k. \quad (20) $$

$x_j^i$ and $y_j^i$ are the coordinates of the $j$th particle within the $i$th simulation. All time-dependent dispersion results will be presented against $t/t^*$ defined as the time $t$ normalized by the characteristic time $t^*$ necessary to cross a correlation length, $t^* = \lambda/\langle|v|\rangle$. Second, the average over a number MC of Monte Carlo simulations is performed

$$ \alpha_L(t) = \langle \alpha_L^i(t) \rangle_{i=1,\text{MC}} \quad \text{and} \quad \alpha_T(t) = \langle \alpha_T^i(t) \rangle_{i=1,\text{MC}}. \quad (21) $$

In all the numerical simulations, the number $N_p$ of particles injected is fixed to 10,000. A single convergence study on the number MC of Monte Carlo simulations has been realized in order to obtain the asymptotic regime of macro-dispersion. Figures 11 and 12 show the convergence to the asymptotic regime of the longitudinal dispersivity $\alpha_L$ with the number MC of Monte Carlo simulations performed in the advective–diffusive case. The transverse dispersivity $\alpha_T$ has not been plotted because all the numerical simulations give a zero value. Taking $\sigma_h = 3$ mm, the single rough fractures are generated with $b_m = 4$ mm (see Figure 11) and 12 mm (see Figure 12). The diffusion coefficient $D_m$ is fixed to 0.01 mm²/s. The statistical expectation of the mean flow velocity $\langle|v|\rangle$ is equal to 0.244 mm/s for $b_m = 4$ mm giving $Pe = 97.6$ (see Figure 11) and 0.066 mm/s
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$ for different values of the number MC of Monte Carlo simulations in single rough fractures generated with $b_m = 4$ mm, $\sigma_h = 3$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm in the advective-diffusive case $Pe = 97.6$. (The horizontal line represents the asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity.)</figcaption></figure>

for $b_m = 12$ mm giving $Pe = 79.2$ (see Figure 12). In the two cases, we can observe that the conver-
gence of the longitudinal dispersivity $\alpha_L$ on the asymptotic regime is very fast with the number
MC of Monte Carlo simulations performed. The asymptotic regime seems to be established for
$MC > 20$. In the following work, we have selected $MC = 50$ in order to obtain correctly the asymp-
totic regime. This gives an asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity estimated
from 5000,000 particle trajectories in the single rough fractures. The asymptotic value $\alpha_{\text{Lasymp}}^R$
of the longitudinal dispersivity, indicated in Figures 11 and 12 with the horizontal line, is deter-
mined by averaging the time evolution curve of the longitudinal dispersivity $\alpha_L(t/t^*)$ over the
time range $[0.5 t_{fb}, t_{fb}]$ over which $\alpha_L(t/t^*)$ is observed to have reached its asymptotic regime,
where $t_{fb}$ is the first breakthrough time (time for which the first particle arrives at a distance of
$0.05 L_x$ from the output border) [44–48].

**6. Numerical results**

In this work, two values of $b_m$ are tested, 4 mm and 12 mm. The standard deviation $\sigma_h$ of $h^\pm$ is varied between 0.1 and 1 for $b_m = 4$ mm and between 0.5 and 3 for $b_m = 12$ mm. The hydraulic heads, in the inlet and outlet of the computational domain, are fixed such that the Reynolds assumption is satisfied. The statistical expectation of the mean flow velocity $\langle|v|\rangle$ is equal to 0.244 mm/s with $b_m = 4$ mm and 0.066 mm/s with $b_m = 12$ mm. It corresponds to a Reynolds number $Re = 0.976$ with $b_m = 4$ mm and 0.792 with $b_m = 12$ mm. The coefficient of molecular diffusion is fixed to $D_m = 0$ ($Pe = \infty$) and 0.01 mm$^2$/s ($Pe = 97.6$ for $b_m = 4$ mm and 79.2 for $b_m = 12$ mm). The correlation length is fixed to $\lambda = 4$ mm. The characteristics of all the numerical simulations performed in this work are summarized in the Table 1.

From all the numerical simulations performed in this work, Figures 13, 14, 15, and 16 show the
longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$ for different values of
$\sigma_h$. Figures 13 and 14 correspond to $D_m = 0$ and 0.01 mm$^2$/s with $b_m = 4$ mm. Figures 15 and 16
correspond to $D_m = 0$ and 0.01 with $b_m = 12$ mm. We can observe that the asymptotic regime is
<figure><img src="image_17.png" /><figcaption>**Figure 12.** Longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$ for different values of the number MC of Monte Carlo simulations in single rough fractures generated with $b_m = 12$ mm, $\sigma_h = 3$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm in the advective–diffusive case $Pe = 79.2$. (The horizontal line represents the asymptotic value $\alpha_{Lasymp}^R$ of the longitudinal dispersivity.)</figcaption></figure>

<table><caption>**Table 1.** Characteristics of numerical simulations performed in this work</caption><thead><tr><th>Parameters</th><th colspan="2">Values</th></tr></thead><tbody><tr><td>$b_m$ (mm)</td><td>4</td><td>12</td></tr><tr><td>$\sigma_h$ (mm)</td><td>0.1, 0.25, 0.5, 0.75, 0.875, 1</td><td>0.5, 1, 1.5, 2, 2.5, 3</td></tr><tr><td>$\langle|\mathbf{v}|\rangle$ (mm/s)</td><td>0.244</td><td>0.066</td></tr><tr><td>$Re$</td><td>0.976</td><td>0.792</td></tr><tr><td>$D_m$ (mm<sup>2</sup>/s)</td><td>0 &rarr; ($Pe = \infty$) 0.01 &rarr; ($Pe = 97.6$)</td><td>0 &rarr; ($Pe = \infty$) 0.01 &rarr; ($Pe = 79.2$)</td></tr><tr><td>$\lambda$ (mm)</td><td colspan="2">4</td></tr><tr><td>$L_x \times L_y$ (mm × mm)</td><td colspan="2">4096 × 1024</td></tr><tr><td>$\Delta$ (mm)</td><td colspan="2">1</td></tr><tr><td>$N_p$</td><td colspan="2">10,000</td></tr><tr><td>MC</td><td colspan="2">50</td></tr></tbody></table>

reached for all the cases studied. In Figures 13, 14, 15, and 16, the horizontal line corresponds still to the asymptotic value $\alpha_{Lasymp}^R$ of the longitudinal dispersivity.

This asymptotic value $\alpha_{Lasymp}^R$ is presented in Figures 17 and 18 as a function of $\sigma_{lnb}$ for $b_m = 4$ mm and 12 mm, respectively. In the two figures, the analytical solution proposed by Gelhar in 1993 has been plotted [73]. Based on a stochastic analysis of solute transport in single rough fractures, the expression of this analytical solution is given by

$$ \alpha_{\text{anal}} = \frac{[3 + G(\sigma_{\text{lnb}}^2)]\sigma_{\text{lnb}}^2 \lambda + D_m}{\langle|\mathbf{v}|\rangle} \quad \text{with } \sigma_{\text{lnb}} = \sqrt{\ln\left(\frac{\sigma_b^2}{b_m^2} + 1\right)}. \qquad (22) $$

In the previous equation, the function G is defined by

$$ G(x) = 1 + 0.205x^2 + 0.16x^4 + 0.045x^6 + 0.0115x^8 \quad \text{with } 0 < x < 5. \qquad (23) $$
<figure><img src="image_18.png" /><figcaption>**Figure 13.** Longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$ for different values of $\sigma_h$ with $D_m = 0 \text{ mm}^2/\text{s}$ (or $Pe = \infty$), $b_m = 4$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm. (The horizontal lines present the asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity.)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$ for different values of $\sigma_h$ with $D_m = 0.01 \text{ mm}^2/\text{s}$ (or $Pe = 97.6$), $b_m = 4$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, and $\Delta = 1$ mm. (The horizontal lines present the asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity.)</figcaption></figure>

This analytical solution is based on two assumptions. The logarithm of the fracture aperture is a statistically stationary, isotropic, 2D Gaussian random field. The flow velocity through a single rough fracture can be modeled by the Reynolds equation. These two assumptions are verified by the present numerical model, making comparison possible. In the right-hand side of (22), the first member represents the effect of the fracture roughness on the macro-dispersion. In Figures 17 and 18, we can observe that the numerical results have a good agreement with the analytical
<figure><img src="image_20.png" /><figcaption>**Figure 15.** Longitudinal dispersivity $a_L$ as a function of the dimensionless time $t/t^*$ for different values of $\sigma_h$ with $D_m = 0 \text{ mm}^2/\text{s}$ (or $Pe = \infty$), $b_m = 12 \text{ mm}$, $\lambda = 4 \text{ mm}$, $L_x = 4096 \text{ mm}$, $L_y = 1024 \text{ mm}$, and $\Delta = 1 \text{ mm}$. (The horizontal lines present the asymptotic value $a_{\text{Lasymp}}^R$ of the longitudinal dispersivity.)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 16.** Longitudinal dispersivity $a_L$ as a function of the dimensionless time $t/t^*$ for different values of $\sigma_h$ with $D_m = 0.01 \text{ mm}^2/\text{s}$ (or $Pe = 79.2$), $b_m = 12 \text{ mm}$, $\lambda = 4 \text{ mm}$, $L_x = 4096 \text{ mm}$, $L_y = 1024 \text{ mm}$, and $\Delta = 1 \text{ mm}$. (The horizontal lines present the asymptotic value $a_{\text{Lasymp}}^R$ of the longitudinal dispersivity.)</figcaption></figure>

solution for $\sigma_{\text{Inb}} < 0.25$. In the contrary case, a difference between the numerical results and the analytical solution increases with $\sigma_{\text{Inb}}$. This shows that the macro-dispersion does indeed result from the sum of effects of the fracture roughness and molecular diffusion up to $\sigma_{\text{Inb}} < 0.25$, where the fraction void $S_O$ remains equal to 1. Tables 2 and 3 summarize the asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity in the single rough fractures and the absolute difference $|\Delta\alpha|$ between the numerical and analytical asymptotic values of the longitudinal dispersivity in the single rough fractures.
<figure><img src="image_22.png" /><figcaption>**Figure 17.** Asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity in the single rough fractures and analytical solution $\alpha_{\text{anal}}$ proposed by Gelhar in 1993, as functions of $\sigma_{\text{lnb}}$ for $b_m = 4$ mm and two values of number Péclet $Pe = \infty$ and 97.6.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 18.** Asymptotic value $\alpha_{\text{Lasymp}}^R$ of the longitudinal dispersivity in the single rough fractures and analytical solution $\alpha_{\text{anal}}$ proposed by Gelhar in 1993, as functions of $\sigma_{\text{lnb}}$ for $b_m = 12$ mm and two values of number Péclet $Pe = \infty$ and 79.2.</figcaption></figure>

**7. Discussion**

Figures 19 and 20 show the fractional void $S_O$ as a function of $\sigma_{\text{lnb}}$ for $b_m = 4$ and 12 mm, respectively. The numerical results are obtained by counting the number of mesh cells with $b = 0$ mm in the computational domain to estimate (6). Tables 2 and 3 summarize the value of the fractional void $S_O$ in the single rough fractures. The analytical results are given by (8). The numerical and analytical results are relatively close. The value of $S_O$ becomes different to 1 for $\sigma_{\text{lnb}} > 0.25$. From this limit $\sigma_{\text{lnb}} = 0.25$, the roughness of two fracture surfaces is strong enough to induce a non-zero surface contact area. In Figure 3, this is what happens for the single rough fracture with $b_m = 12$ mm and $\sigma_h = 3$ mm.
<table><caption>**Table 2.** Fractional void $S_O$, asymptotic value $\alpha_{\text{Lasymp}}^i$ of the longitudinal dispersivity in the single rough ($i = R$) and equivalent ($i = E$) fractures, and absolute difference $|\Delta\alpha|$ between the numerical and analytical asymptotic values of the longitudinal dispersivity in the single rough fractures as functions of $\sigma_{\text{Inb}}$ for $b_m = 4$ mm and two values of $Pe$</caption><thead><tr><th rowspan="2">σ<sub>Inb</sub> (mm)</th><th rowspan="2">S<sub>O</sub></th><th colspan="3">Pe = ∞</th><th colspan="3">Pe = 97.6</th></tr><tr><th>α<sup>R</sup><sub>Lasymp</sub> (mm)</th><th>|Δα| (mm)</th><th>α<sup>E</sup><sub>Lasymp</sub> (mm)</th><th>α<sup>R</sup><sub>Lasymp</sub> (mm)</th><th>|Δα| (mm)</th></tr></thead><tbody><tr><td>0.032</td><td>1</td><td>0.009</td><td>0.007</td><td>0</td><td>0.062</td><td>0.005</td></tr><tr><td>0.085</td><td>1</td><td>0.078</td><td>0.036</td><td>0</td><td>0.16</td><td>0.002</td></tr><tr><td>0.172</td><td>1</td><td>0.486</td><td>0.013</td><td>0</td><td>0.556</td><td>0.036</td></tr><tr><td>0.257</td><td>0.9999</td><td>1.114</td><td>0.052</td><td>0.028</td><td>1.157</td><td>0.049</td></tr><tr><td>0.299</td><td>0.9997</td><td>1.581</td><td>0.14</td><td>0.12</td><td>1.628</td><td>0.141</td></tr><tr><td>0.342</td><td>0.9985</td><td>2.488</td><td>0.609</td><td>0.64</td><td>2.5</td><td>0.571</td></tr></tbody></table>

<table><caption>**Table 3.** Fractional void $S_O$, asymptotic value $\alpha_{\text{Lasymp}}^i$ of the longitudinal dispersivity in the single rough ($i = R$) and equivalent ($i = E$) fractures, and absolute difference $|\Delta\alpha|$ between the numerical and analytical asymptotic values of the longitudinal dispersivity in the single rough fractures as functions of $\sigma_{\text{Inb}}$ for $b_m = 12$ mm and two values of $Pe$</caption><thead><tr><th rowspan="2">σ<sub>Inb</sub> (mm)</th><th rowspan="2">S<sub>O</sub></th><th colspan="3">Pe = ∞</th><th colspan="3">Pe = 79.2</th></tr><tr><th>α<sup>R</sup><sub>Lasymp</sub> (mm)</th><th>|Δα| (mm)</th><th>α<sup>E</sup><sub>Lasymp</sub> (mm)</th><th>α<sup>R</sup><sub>Lasymp</sub> (mm)</th><th>|Δα| (mm)</th></tr></thead><tbody><tr><td>0.059</td><td>1</td><td>0.046</td><td>0.01</td><td>0</td><td>0.235</td><td>0.028</td></tr><tr><td>0.117</td><td>1</td><td>0.2</td><td>0.021</td><td>0</td><td>0.439</td><td>0.064</td></tr><tr><td>0.175</td><td>1</td><td>0.469</td><td>0.024</td><td>0</td><td>0.703</td><td>0.054</td></tr><tr><td>0.233</td><td>0.9999</td><td>0.917</td><td>0.024</td><td>0.022</td><td>1.079</td><td>0.053</td></tr><tr><td>0.289</td><td>0.9998</td><td>1.467</td><td>0.125</td><td>0.12</td><td>1.707</td><td>0.205</td></tr><tr><td>0.345</td><td>0.9985</td><td>2.65</td><td>0.734</td><td>0.8</td><td>2.848</td><td>0.764</td></tr></tbody></table>

The surface contact area results then in contacts or obstacles that the flow has to bypass. The flow through a single rough fracture occurs in a flow path system that is separated by the surface contact area, as shown in Figure 6. This mechanism can be compared to the one observed in the evolution of the surface contact area during significant overload or large displacement events in the fractures. Experimental and numerical works have shown that the increase of the fracture roughness and the fracture contact area leads to tortuous flow paths in the single rough fracture because the flow has to bypass the contacts or obstacles. This leads to a non-Fickian transport with an increase of macro-dispersion [12, 37, 39–41, 74].

To estimate the impact of the fracture contact area on macro-dispersion, we have estimated the longitudinal dispersivity $\alpha_L$ in single equivalent fractures generated from single rough fractures by conserving the same values of the fractional void $S_O$ and distance between the two average planes $b_m$. Figure 8 shows an example of a single equivalent fracture generated from a single rough fracture with $b_m = 12$ mm and $\sigma_h = 3$ mm. As for the single rough fractures, each value of the longitudinal dispersivity $\alpha_L$ is estimated from 500,000 particle trajectories (10,000 particles injected for each Monte Carlo simulation with MC = 50) in the same injection conditions. The mechanism of advection is only considered in the single equivalent fracture.

From these numerical simulations, Figures 21 and 22 show the longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$ for $b_m = 4$ mm with $\sigma_h = 0.75, 0.875$, and 1, and $b_m = 12$ mm with $\sigma_h = 2, 2.5$ and 3. In all the cases studied, the asymptotic regime of the longitudinal dispersivity $\alpha_L$ is established. The averaging method is still used for the estimation of
<figure><img src="image_24.png" /><figcaption>**Figure 19.** Fractional void $S_O$ as a function of $\sigma_{\text{Inb}}$ for a single rough fracture with $L_x = 4096$ mm, $L_y = 1024$ mm, $\Delta = 1$ mm, $\lambda = 4$ mm and $b_m = 4$ mm. (Triangles: analytical value given by (8) and diamonds: numerical value given by the present numerical model.)</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 20.** Fractional void $S_O$ as a function of $\sigma_{\text{Inb}}$ for a single rough fracture with $L_x = 4096$ mm, $L_y = 1024$ mm, $\Delta = 1$ mm, $\lambda = 4$ mm and $b_m = 12$ mm. (Triangles: analytical value given by (8) and diamonds: numerical value given by the present numerical model.)</figcaption></figure>

$\alpha_{\text{Lasymp}}^E$ in the single equivalent fractures. Tables 2 and 3 summarize the asymptotic value $\alpha_{\text{Lasymp}}^E$ of the longitudinal dispersivity in the single equivalent fractures.

In Figures 23 and 24, we have reported the asymptotic value $\alpha_{\text{Lasymp}}^E$ of the longitudinal dispersivity in the single equivalent fractures and the absolute difference $|\Delta\alpha|$ between the numerical and analytical asymptotic values of the longitudinal dispersivity in the single rough fractures as functions of $\sigma_{\text{Inb}}$ in the purely advective case ($Pe = \infty$) and the advective-diffusive case ($Pe < 100$) for $b_m = 4$ and 12 mm, respectively. We can observe that the two quantities, $\alpha_{\text{Lasymp}}^E$ and $|\Delta\alpha|$ are relatively close. This shows that the absolute difference $|\Delta\alpha|$ between the numerical and analytical asymptotic values of the longitudinal dispersivity in the single rough fractures is caused by the presence of contacts disrupting the flow velocity fields.
<figure><img src="image_26.png" /><figcaption>**Figure 21.** Longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$, in a single equivalent fracture generated from a single rough fracture with $b_m = 4$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, $\Delta = 1$ mm, and different values of $\sigma_h$ for the purely advective case, $D_m = 0 \text{ mm}^2/\text{s}$ or $Pe = \infty$. (The horizontal lines present the asymptotic value $\alpha_{\text{Lasymp}}^E$ of the longitudinal dispersivity.)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 22.** Longitudinal dispersivity $\alpha_L$ as a function of the dimensionless time $t/t^*$, in a single equivalent fracture generated from a single rough fracture with $b_m = 12$ mm, $\lambda = 4$ mm, $L_x = 4096$ mm, $L_y = 1024$ mm, $\Delta = 1$ mm, and different values of $\sigma_h$ for the purely advective case, $D_m = 0 \text{ mm}^2/\text{s}$ or $Pe = \infty$. (The horizontal lines present the asymptotic value $\alpha_{\text{Lasymp}}^E$ of the longitudinal dispersivity.)</figcaption></figure>

**8. Conclusions**

The impact of the fracture contact area on macro-dispersion in single rough fractures has been highlighted by means of Monte Carlo simulations in which four steps are performed:
<figure><img src="image_28.png" /><figcaption>**Figure 23.** Asymptotic value $\alpha_{\text{Lasymp}}^E$ of the longitudinal dispersivity in the single equiv- alent fractures and absolute difference $|\Delta\alpha|$ between the numerical and analytical asymp- totic values of the longitudinal dispersivity in the single rough fractures, as functions of $\sigma_{\text{Inb}}$ for $b_m = 4$ mm and two values of $Pe$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 24.** Asymptotic value $\alpha_{\text{Lasymp}}^E$ of the longitudinal dispersivity in the single equiv- alent fractures and absolute difference $|\Delta\alpha|$ between the numerical and analytical asymp- totic values of the longitudinal dispersivity in the single rough fractures, as functions of $\sigma_{\text{Inb}}$ for $b_m = 12$ mm and two values of $Pe$.</figcaption></figure>

step 1—fracture aperture generation; step 2—flow velocity calculation; step 3—particle tracking simulation; and step 4—macro-dispersion estimation.

When the fractional void $S_O$ is equal to 1 (i.e., for $\sigma_{\text{Inb}} < 0.25$), the Monte Carlo simulations show that macro-dispersion results of two contributions, dispersion caused by the heterogeneity of fracture apertures that induces a channelization of flow paths and molecular diffusion, as shown by the analytical solution proposed by Gelhar in 1993 (22) [73] (see Figures 17 and 18). When the fraction void $S_O$ is different from 1 (i.e., for $\sigma_{\text{Inb}} > 0.25$), a third mechanism plays in macro-dispersion. The presence of contacts or obstacles causes a disruption of flow paths. This
mechanism is identical to that induced by the fracture roughness with a lower amplitude (see Figures 23 and 24). Its amplitude is function of the fractional void $S_O$.

These conclusions are specific to the model conditions studied. The single rough fractures are made with two Gaussian surfaces. The generation of single rough fractures is based on the model developed by Mourzenko et al. [30, 31, 49] (see (1)–(10)). This model allows also to generate self-affine fractures. It has not been discussed how macro-dispersion varies if the generation of single rough fractures is changed. Macro-dispersion is estimated from moments that need to evaluate particle trajectories (see (18)–(21)). The particle trajectories are evaluated by considering only advection and molecular diffusion (see (14) and (15)). The Taylor-Aris dispersion is also present in the fractures. The impact of Taylor-Aris dispersion on macro-dispersion has not been studied in this work.

## Acknowledgments

This work was granted access to the high performance computing resources of IDRIS under the allocation 2021-A0092A10954 made by GENCI. The action 2—Dynamique des contaminants in the axis 3 Chimie verte et préservation des ressources (eau, sol, biodiversité, carbone renouvelable) of the ECONAT program—gestion durable des ECOsystèmes et des ressources NATurelles (CPER 2015–2020) is acknowledged for its financial funding.

## References

[1] R. W. Zimmerman, I. Main, “Hydromechanical behavior of fractured rocks”, in Mechanics of Fluid-saturated Rocks (Y. Gueguen, M. Bouteca, eds.), Elsevier, London, 2004, p. 363-421.

[2] Q. Lei, J. P. Latham, C. F. Tsang, “The use of discrete fracture networks for modelling coupled geomechanical and hydrological behaviour of fractured rocks”, Comput. Geotech. 85 (2017), p. 151-176.

[3] Z. Dou, B. Sleep, H. Zhan, Z. Zhou, J. Wang, “Multiscale roughness influence on conservative solute transport in self-affine fracture”, Int. J. Heat Mass Trans. 133 (2019), p. 606-618.

[4] Z. Yang, D. Li, S. Xue, R. Hu, Y. F. Chen, “Effect of aperture field anisotropy on two-phase in rough fractures”, Adv. Water Res. 132 (2019), article no. 103390.

[5] W. A. Illman, “Hydraulic tomography offers improved imaging of heterogeneity in fractured rocks”, Groundwater 52 (2014), p. 659-684.

[6] L. Wang, M. Bayani Cardenas, “An efficient quasi-3D particle tracking-based approach for transport through fractures with application to dynamic dispersion calculation”, J. Contam. Hydrol. 179 (2015), p. 47-54.

[7] L. Wang, M. Bayani Cardenas, “Transition from non-Fickian to Fickian longitudinal transport through 3-D rough fractures: Scale-(in) sensitivity and roughness dependence”, J. Contam. Hydrol. 198 (2017), p. 1-10.

[8] S. P. Neuman, D. M. Tartakovsky, “Perspective on theories of non-Fickian transport in heterogeneous media”, Adv. Water Res. 32-5 (2009), p. 670-680.

[9] Q. Zheng, S. Dickson, Y. Guo, “Influence of aperture field heterogeneity and anisotropy on dispersion regimes and dispersivity in single fractures”, J. Geophys. Res. 114 (2009), no. B3, p. 12.

[10] L. Wang, M. B. Cardenas, “Non-Fickian transport through two-dimensional rough fractures: assessment and prediction”, Water Resour. Res. 50 (2014), no. 2, p. 871-884.

[11] Y. Ding, X. Meng, D. Yang, “Numerical simulation of polydisperse dense particles transport in a random-orientated fracture with spatially variable apertures”, Colloids Surf. A: Physicochem. Eng. Asp. 610 (2020), article no. 125729.

[12] Y. Hu, W. Xu, L. Zhan, Z. Ye, Y. Chen, “Non-Fickian solute transport in rough-walled fractures: the effect of contact area”, Water 12 (2020), no. 7, p. 2049-2074.

[13] S. C. James, C. V. Chrysikopoulos, “Transport of polydisperse colloid suspensions in a single fracture”, Water Resour. Res. 35 (1999), no. 3, p. 707-718.

[14] S. C. James, C. V. Chrysikopoulos, “Transport of polydisperse colloids in a saturated fracture with spatially variable aperture”, Water Resour. Res. 36 (2000), no. 6, p. 1457-1465.

[15] S. C. James, C. V. Chrysikopoulos, “An efficient particle tracking equation with specified spatial step for the solution of the diffusion equation”, Chem. Eng. Sci. 56 (2001), no. 23, p. 6535-6543.

[16] P. W. Reimus, S. C. James, “Determining the random time step in a constant spatial step particle tracking algorithm”, Chem. Eng. Sci. 57 (2002), no. 21, p. 4429-4434.
[17] C. V. Chrysikopoulos, S. C. James, "Transport of neutrally buoyant and dense variably sized colloids in a two-dimensional fracture with anisotropic aperture", *Transp. Porous Med.* **51** (2003), p. 191-210.

[18] S. C. James, C. V. Chrysikopoulos, "Effective velocity and effective dispersion coefficient for finite-sized particles flowing in a uniform fracture", *J. Colloid Interface Sci.* **263** (2003), no. 1, p. 288-295.

[19] S. C. James, C. V. Chrysikopoulos, "Monodisperse and polydisperse colloid transport in water-saturated fractures with various orientations: gravity effects", *Adv. Water Resour.* **34** (2011), no. 10, p. 1249-1255.

[20] X. Meng, D. Yang, "Determination of dynamic dispersion coefficient for particles flowing in a parallel-plate fracture", *Colloids Surf. A: Physicochem. Eng. Asp.* **509** (2016), p. 259-278.

[21] S. C. James, L. Wang, C. V. Chrysikopoulos, "Modeling colloid transport in fractures with spatially variable aperture and surface attachment", *J. Hydrol.* **566** (2018), p. 735-742.

[22] L. Zheng, L. Wang, C. James, "When can the local advection-dispersion equation simulate non-Fickian transport through fractures?", *Stoch. Environ. Res. Risk Assess.* **33** (2019), no. 3, p. 931-938.

[23] Y. Ding, X. Meng, D. Yang, "Determination of dynamic dispersion coefficient for solid particles flowing in a fracture with consideration of gravity effect", *J. Energy Resour. Technol.* **145** (2020), no. 2, article no. 053101.

[24] D. J. Brush, N. R. Thomson, "Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations", *Water Resour. Res.* **39** (2003), no. 4, p. 1085-1100.

[25] A. H. Al-Yaarubi, C. C. Pain, C. A. Grattoni, R. W. Zimmerman, "Navier-Stokes simulations of fluid flow through a rock fracture", in *Dynamics of Fluids and Transport in Fractured Rock*, Geophysical Monograph Series, vol. 162, AGU, Washington, DC, 2020, p. 55-64.

[26] P. W. Reimus, "The use of synthetic colloids in tracer transport experiments in saturated rock fractures", 1995, Thesis, University of New Mexico, Albuquerque, New Mexico.

[27] R. L. Detwiler, H. Rajaram, R. J. Glass, "Solute transport in variable-aperture fractures: an investigation of the relative importance of Taylor dispersion and macrodispersion", *Water Resour. Res.* **36** (2000), no. 7, p. 1611-1625.

[28] S. C. James, T. Bilezikjian, C. Chrysikopoulos, "Contaminant transport in a fracture with spatially variable aperture in the presence of monodisperse and polydisperse colloids", *Stoch. Env. Res. Risk A* **19** (2005), no. 4, p. 266-279.

[29] D. D. Pollard, A. Aydin, "Progress in understanding jointing over the past century", *Geol. Soc. Am. Bull.* **100** (1988), p. 1181-1204.

[30] V. V. Mourzenko, J. F. Thovert, P. M. Adler, "Geometry of simulated fractures", *Phys. Rev. E* **53** (1996), no. 6, p. 5606-5626.

[31] P. M. Adler, J. F. Thovert, *Fractures and Fracture Networks*, vol. 15, Springer Science & Business Media, 1999.

[32] R. Schultz, "Growth of geologic fractures into large-strain populations: review of nomenclature, subcritical crack growth, and some implications for rock engineering", *Int. J. Rock Mech. Min. Sci.* **37** (2000), p. 403-411.

[33] Q. Lei, X. Wang, "Tectonic interpretation of the connectivity of a multiscale fracture system in limestone", *Geophys. Res. Lett.* **43** (2016), p. 1551-1558.

[34] J. Rutqvist, O. Stephansson, "The role of hydromechanical coupling in fractured rock engineering", *Hydrogeol. J.* **11** (2003), p. 7-40.

[35] C. F. Tsang, R. Bernier, C. Davies, "Geohydromechanical processes in the excavation damaged zone in crystalline rock, rock salt, and indurated and plastic clays—in the context of radioactive waste disposal", *Int. J. Rock Mech. Min. Sci.* **42** (2005), p. 109-125.

[36] S. P. Bertels, D. A. DiCarlo, M. J. Blunt, "Measurement of aperture distribution, capillary pressure, relative permeability, and in situ saturation in a rock fracture using computed tomography scanning", *Water Resour. Res.* **37** (2001), p. 649-662.

[37] J. Zhou, S. Hu, S. Fang, Y. Chen, C. Zhou, "Nonlinear flow behavior at low Reynolds numbers through rough-walled fractures subjected to normal compressive loading", *Int. J. Rock Mech. Min. Sci.* **80** (2015), p. 202-218.

[38] D. Huo, S. M. Benson, "Experimental investigation of stress-dependency of relative permeability in rock fractures", *Transp. Porous Med.* **113** (2016), p. 567-590.

[39] Y. Chen, H. Lian, W. Liang, J. Yang, V. P. Nguyen, S. P. A. Bordas, "The influence of fracture geometry variation on non-Darcy flow in fractures under confining stresses", *Int. J. Rock Mech. Min. Sci.* **113** (2019), p. 59-71.

[40] W. Jeong, J. Song, "Numerical investigations for flow and transport in a rough fracture with a hydromechanical effect", *Energy Sources B* **27** (2005), p. 997-1011.

[41] T. Koyama, B. Li, Y. Jiang, L. Jing, "Numerical simulations for the effects of normal loading on particle transport in rock fractures during shear", *Int. J. Rock Mech. Min. Sci.* **45** (2008), p. 1403-1419.

[42] P. K. Kang, S. Brown, R. Juanes, "Emergence of anomalous transport in stressed rough fractures", *Earth Planet Sci. Lett.* **454** (2016), p. 46-54.

[43] L. Zou, L. Jing, V. Cvetkovic, "Modeling of solute transport in a 3D rough-walled fracture–matrix system", *Transp. Porous Med.* **116** (2017), p. 1005-1029.

[44] J. R. de Dreuzy, A. Beaudoin, J. Erhel, "Asymptotic dispersion in 2D heterogeneous porous media determined by parallel numerical simulations", *Water Resour. Res.* **43** (2007), no. 10, p. 13.

[45] J. R. de Dreuzy, A. Beaudoin, J. Erhel, "Reply to comment by A. Fiori et al. on Asymptotic dispersion in 2D
heterogeneous porous media determined by parallel numerical simulations", *Water Resour. Res.* **44** (2008), no. 6, article no. W06604.

[46] A. Beaudoin, J. R. de Dreuzy, J. Erhel, "Numerical Monte Carlo analysis of the influence of pore-scale dispersion on macrodispersion in 2-D heterogeneous porous media", *Water Resour. Res.* **46** (2010), no. 12, p. 9.

[47] A. Beaudoin, J. R. de Dreuzy, "Numerical assessment of 3-D macrodispersion in heterogeneous porous media", *Water Resour. Res.* **49** (2013), no. 5, p. 2489-2496.

[48] A. Beaudoin, A. Dartois, S. Huberson, "Analysis of the influence of averaged positive second invariant Qav of deformation tensor u on the maximum dilution index Emax in steady Darcy flows through isotropic heterogeneous porous media", *Adv. Water Res.* **128** (2019), p. 39-47.

[49] V. V. Mourzenko, J. F. Thovert, P. M. Adler, "Conductivity and transmissivity of a single fracture", *Tranp. Porous Med.* **112** (2016), no. 3, p. 235-256.

[50] S. R. Brown, R. L. Kranz, B. P. Bonner, "Correlation between the surfaces of natural rock joints", *Geophys. Res. Lett.* **13** (1986), p. 1430-1433.

[51] S. K. Sinha, E. B. Sirota, S. Garoff, H. B. Stanley, "X-ray and neutron scattering from rough surfaces", *Phys. Rev. B* **38** (1988), p. 2297-2311.

[52] B. B. Mandelbrot, *The Fractal Geometry of Nature*, Freeman, San Francisco, 1983.

[53] E. Bouchaud, J. P. Bouchaud, "Fracture surfaces: apparent roughness, relevant length scales, and fracture toughness", *Phys. Rev. B* **17** (1994), p. 17752.

[54] A. L. Gutjahr, "Fast fourier transforms for random field generation", 1989, New Mexico Institute of Mining and Technology Project Report for Los Alamos Grant to New Mexica Tech (Doctoral Dissertation).

[55] E. Pardo-Iguzquiza, M. Chica-Olmo, "The fourier integral method: an efficient spectral method for simulation of random fields", *Math. Geol.* **25** (1993), no. 2, p. 177-217.

[56] F. Hebe, V. Prykhodko, S. Schluter, S. Attinger, "Generating random fields with a truncated power-law variogram: a comparison of several numerical methods", *Environ. Modell. Softw.* **55** (2014), p. 32-48.

[57] M. Frigo, S. G. Johnson, "The design and implementation of FFTW3", *Proc. IEEE* **93** (2005), no. 2, p. 216-231.

[58] M. Roberts, J. C. Bowman, "Multithreaded implicitly deadliased convolutions", *J. Comp. Phys.* **356** (2018), p. 98-114.

[59] I. W. Yeo, M. H. De Freitas, R. W. Zimmerman, "Effect of shear displacement on the aperture and permeability of rock", *Int. J. Rock Mech. Min. Sci.* **35** (1998), p. 1051-1070.

[60] R. Eymard, T. Gallouet, R. Herbin, "Finite volume methods", in *Handbook of Numerical Analysis*, vol. 7, Elsevier, 2000, p. 713-1018.

[61] R. Eymard, T. Gallouet, R. Herbin, "A new finite volume scheme for anisotropic diffusion problems on general grids: convergence analysis", *C. R. Math. 344* (2007), no. 6, p. 403-406.

[62] J. Erhel, J. R. de Dreuzy, A. Beaudoin, E. Bresciani, D. Tromeur-Dervout, "A parallel scientific software for heterogeneous hydrogeology", in *Paralle Computational Fluid Dynamics*, Springer, Berlin, Heidelberg, 2009.

[63] R. D. Falgout, J. E. Jones, U. M. Yang, "Pursuing scalability for hypre's conceptual interfaces", *ACM Trans. Math. Softw.* (*TOMS*) **31** (2005), no. 3, p. 326-350.

[64] V. V. Mourzenko, J. F. Thovert, P. M. Adler, "Permeability of a single fracture: validity of the Reynolds equation", *J. Phys. II* **5** (1995), no. 3, p. 465-482.

[65] L. Moreno, Y. W. Tsang, C. F. Tsang, F. Hale, I. Neretnieks, "Flow and tracer transport in a single fracture: A stochastic model and its relation to some field observations", *Water Resour. Res.* **24** (1988), no. 12, p. 2033-2048.

[66] M. W. Becker, A. M. Shapiro, "Tracer transport in fractured crystalline rock: evidence of nondiffusive breakthrough tailing", *Water Resour. Res.* **36** (2000), p. 1677-1686.

[67] A. F. B. Tompson, L. W. Gelhar, "Numerical simulation of solute transport in three-dimensional randomly heterogeneous porous media", *Water Resour. Res.* **26** (1990), no. 10, p. 2541-2562.

[68] M. Prevost, M. G. Edwards, M. J. Blunt, "Streamline tracing on curvilinear structured and unstructured grids", *SPE J.* **7** (2002), no. 2, p. 139-148.

[69] D. W. Pollock, "Semianalytical computation of path lines for finite-difference models", *Groundwater* **26** (1988), no. 6, p. 743-750.

[70] Y. Kuznetsov, S. Repin, "New mixed finite element method on polygonal and polyhedral meshes", *Russ. J. Numer. Anal. Math. Modell.* **18** (2003), p. 261-278.

[71] E. Jimenez, K. Sabir, A. Datta-Gupta, M. J. King, "Spatial error and convergence in streamline simulation", in *SPE Reservoir Simulation Symposium*, Society of Petroleum Engineers, 2005.

[72] H. M. Cheng, J. Droniou, "An HMM-ELLAM scheme on generic polygonal meshes for miscible incompressible flows in porous media", *J. Pet. Sci. Eng.* **172** (2019), p. 707-723.

[73] L. W. Gelhar, *Stochastic Subsurface Hydrology*, Prentice-Hall, Old Tappan, NJ, 1993, 272-279 pages.

[74] D. Vogler, R. R. Settgast, C. Annavarapu, C. Madonna, P. Bayer, F. Amann, "Experiments and simulations of fully hydro-mechanically coupled response of rough fractures exposed to high-pressure fluid injection", *J. Geophys. Res.* **123** (2018), p. 1186-1200.